C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_init_fixed.F,v 1.19 2014/11/25 01:08:20 jmc Exp $
C $Name:  $

#include "OBCS_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: OBCS_INIT_FIXED

C     !INTERFACE:
      SUBROUTINE OBCS_INIT_FIXED( myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | SUBROUTINE OBCS_INIT_FIXED
C     | o Initialise OBCs fixed arrays
C     *==========================================================*

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"

C     !INPUT/OUTPUT PARAMETERS:
C     myThid   :: my Thread Id. number
      INTEGER myThid
CEOP

#ifdef ALLOW_OBCS
C     !LOCAL VARIABLES:
C     msgBuf   :: Informational/error message buffer
C     OB_ApplX :: number of grid points (in X dir) overwritten by obcs_apply
C     OB_ApplY :: number of grid points (in Y dir) overwritten by obcs_apply
C     bi,bj    :: tile indices
C     i, j     :: Loop counters
      CHARACTER*(MAX_LEN_MBUF) msgBuf, errMsg
      INTEGER OB_ApplX
      INTEGER OB_ApplY
      INTEGER bi, bj
      INTEGER  i, j
      INTEGER im, jm
      INTEGER iB, jB
      LOGICAL flag
      INTEGER ioUnit

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_ENTER('OBCS_INIT_FIXED',myThid)
#endif

C==   Set Interior mask at Cell Center:

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
           OBCS_insideMask(i,j,bi,bj) = 1.
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      IF ( insideOBmaskFile.EQ.' ' ) THEN
C--   If no maskFile is provided, set Interior mask from OB list of indices

       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
C-    Eastern boundary
          IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN
           flag = .TRUE.
           DO i=OB_Ie(j,bi,bj),sNx
             flag = flag .AND.
     &              kSurfC(i,j,bi,bj).LE.Nr .AND. i.NE.OB_Iw(j,bi,bj)
             IF ( flag ) OBCS_insideMask(i,j,bi,bj) = 0.
           ENDDO
          ENDIF
C-    Western boundary
          IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN
           flag = .TRUE.
           DO i=OB_Iw(j,bi,bj),1,-1
             flag = flag .AND.
     &              kSurfC(i,j,bi,bj).LE.Nr .AND. i.NE.OB_Ie(j,bi,bj)
             IF ( flag ) OBCS_insideMask(i,j,bi,bj) = 0.
           ENDDO
          ENDIF
         ENDDO
         DO i=1,sNx
C-    Northern boundary
          IF ( OB_Jn(i,bi,bj).NE.OB_indexNone ) THEN
           flag = .TRUE.
           DO j=OB_Jn(i,bi,bj),sNy
             flag = flag .AND.
     &              kSurfC(i,j,bi,bj).LE.Nr .AND. j.NE.OB_Js(i,bi,bj)
             IF ( flag ) OBCS_insideMask(i,j,bi,bj) = 0.
           ENDDO
          ENDIF
C-    Southern boundary
          IF ( OB_Js(i,bi,bj).NE.OB_indexNone ) THEN
           flag = .TRUE.
           DO j=OB_Js(i,bi,bj),1,-1
             flag = flag .AND.
     &              kSurfC(i,j,bi,bj).LE.Nr .AND. j.NE.OB_Jn(i,bi,bj)
             IF ( flag ) OBCS_insideMask(i,j,bi,bj) = 0.
           ENDDO
          ENDIF
         ENDDO

C--   end bi,bj loops
        ENDDO
       ENDDO

      ELSE
C--   Read in Interior mask from file :

       CALL READ_FLD_XY_RS( insideOBmaskFile, ' ', OBCS_insideMask,
     &                      0, myThid )

       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           IF ( OBCS_insideMask(i,j,bi,bj).NE.0. )
     &          OBCS_insideMask(i,j,bi,bj) = 1.
          ENDDO
         ENDDO
        ENDDO
       ENDDO

C--   end computing/reading Interior mask
      ENDIF

C--   Fill in the overlap:
      _EXCH_XY_RS( OBCS_insideMask, myThid )

C==   Set interior mask at U & V location (grid-cell Wester & Southern edges)
C     leave OB edges inside (maskIn=1) (e.g., Eastern OB: maskInW(OB_Ie)=1 )
C     so that velocity normal-component at OB is still in Interior region.
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=2-OLy,sNy+OLy
          DO i=2-OLx,sNx+OLx
            maskInW(i,j,bi,bj) = maskInW(i,j,bi,bj)
     &                     *MAX( OBCS_insideMask(i-1,j,bi,bj),
     &                           OBCS_insideMask(i,j,bi,bj) )
            maskInS(i,j,bi,bj) = maskInS(i,j,bi,bj)
     &                     *MAX( OBCS_insideMask(i,j-1,bi,bj),
     &                           OBCS_insideMask(i,j,bi,bj) )
          ENDDO
         ENDDO
       ENDDO
      ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   After exchange, set tiled index arrays OB_Jn/Js/Ie/Iw in overlap region
C     issue a warning if not consistent (similar to OBCS_CHECK but for overlap)

c     IF ( .TRUE. ) THEN
      IF ( OBCS_indexStatus .LT. 2 ) THEN
       ioUnit = standardMessageUnit
       WRITE(msgBuf,'(2A)')
     &      'OBCS_INIT_FIXED: Setting OB indices in Overlap'
       CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
       WRITE(errMsg,'(2A)') 'S/R OBCS_INIT_FIXED: ',
     &              'Inside Mask and OB locations disagree :'
       flag = .TRUE.
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)

C-    Eastern boundary
         DO j=1-OLy,sNy+OLy
          DO i=1,sNx+1
           IF ( OBCS_insideMask(i,j,bi,bj).LT.
     &          OBCS_insideMask(i-1,j,bi,bj)
     &          .AND. ( j.LT.1 .OR. j.GT.sNy )
     &          .AND. kSurfW(i,j,bi,bj).LE.Nr ) THEN
            IF ( OB_Ie(j,bi,bj).EQ.OB_indexNone ) THEN
              OB_Ie(j,bi,bj) = i
              WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
     &        ' Sets OBE(j,bi,bj=',j,',',bi,',',bj,')=', OB_Ie(j,bi,bj)
              CALL PRINT_MESSAGE( msgBuf,ioUnit,SQUEEZE_RIGHT,myThid )
            ELSEIF ( OB_Ie(j,bi,bj).NE.i ) THEN
              IF ( flag ) CALL PRINT_ERROR( errMsg, myThid )
              flag = .FALSE.
              WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
     &        ' OBE(j,bi,bj=',j,',',bi,',',bj,')=', OB_Ie(j,bi,bj),
     &        ' but from insideMask expects I=', i
              CALL PRINT_ERROR( msgBuf, myThid )
            ENDIF
           ENDIF
          ENDDO
         ENDDO
C-    Western boundary
         DO j=1-OLy,sNy+OLy
          DO i=0,sNx
           IF ( OBCS_insideMask(i,j,bi,bj).LT.
     &          OBCS_insideMask(i+1,j,bi,bj)
     &          .AND. ( j.LT.1 .OR. j.GT.sNy )
     &          .AND. kSurfW(i+1,j,bi,bj).LE.Nr ) THEN
            IF ( OB_Iw(j,bi,bj).EQ.OB_indexNone ) THEN
              OB_Iw(j,bi,bj) = i
              WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
     &        ' Sets OBW(j,bi,bj=',j,',',bi,',',bj,')=', OB_Iw(j,bi,bj)
              CALL PRINT_MESSAGE( msgBuf,ioUnit,SQUEEZE_RIGHT,myThid )
            ELSEIF ( OB_Iw(j,bi,bj).NE.i ) THEN
              IF ( flag ) CALL PRINT_ERROR( errMsg, myThid )
              flag = .FALSE.
              WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
     &        ' OBW(j,bi,bj=',j,',',bi,',',bj,')=', OB_Iw(j,bi,bj),
     &        ' but from insideMask expects I=', i
              CALL PRINT_ERROR( msgBuf, myThid )
            ENDIF
           ENDIF
          ENDDO
         ENDDO
C-    Northern boundary
         DO j=1,sNy+1
          DO i=1-OLx,sNx+OLx
           IF ( OBCS_insideMask(i,j,bi,bj).LT.
     &          OBCS_insideMask(i,j-1,bi,bj)
     &          .AND. ( i.LT.1 .OR. i.GT.sNx )
     &          .AND. kSurfS(i,j,bi,bj).LE.Nr ) THEN
            IF ( OB_Jn(i,bi,bj).EQ.OB_indexNone ) THEN
              OB_Jn(i,bi,bj) = j
              WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
     &        ' Sets OBN(i,bi,bj=',i,',',bi,',',bj,')=', OB_Jn(i,bi,bj)
              CALL PRINT_MESSAGE( msgBuf,ioUnit,SQUEEZE_RIGHT,myThid )
            ELSEIF ( OB_Jn(i,bi,bj).NE.j ) THEN
              IF ( flag ) CALL PRINT_ERROR( errMsg, myThid )
              flag = .FALSE.
              WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
     &        ' OBN(i,bi,bj=',i,',',bi,',',bj,')=', OB_Jn(i,bi,bj),
     &        ' but from insideMask expects J=', j
              CALL PRINT_ERROR( msgBuf, myThid )
            ENDIF
           ENDIF
          ENDDO
         ENDDO
C-    Southern boundary
         DO j=0,sNy
          DO i=1-OLx,sNx+OLx
           IF ( OBCS_insideMask(i,j,bi,bj).LT.
     &          OBCS_insideMask(i,j+1,bi,bj)
     &          .AND. ( i.LT.1 .OR. i.GT.sNx )
     &          .AND. kSurfS(i,j+1,bi,bj).LE.Nr ) THEN
            IF ( OB_Js(i,bi,bj).EQ.OB_indexNone ) THEN
              OB_Js(i,bi,bj) = j
              WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
     &        ' Sets OBS(i,bi,bj=',i,',',bi,',',bj,')=', OB_Js(i,bi,bj)
              CALL PRINT_MESSAGE( msgBuf,ioUnit,SQUEEZE_RIGHT,myThid )
            ELSEIF ( OB_Js(i,bi,bj).NE.j ) THEN
              IF ( flag ) CALL PRINT_ERROR( errMsg, myThid )
              flag = .FALSE.
              WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
     &        ' OBS(i,bi,bj=',i,',',bi,',',bj,')=', OB_Js(i,bi,bj),
     &        ' but from insideMask expects J=', j
              CALL PRINT_ERROR( msgBuf, myThid )
            ENDIF
           ENDIF
          ENDDO
         ENDDO

        ENDDO
       ENDDO
       WRITE(msgBuf,'(2A)')
     &      'OBCS_INIT_FIXED: Setting OB indices in Overlap <= done'
       CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C==   After EXCH: reset Interior mask to zero beyond OB: this is necessary
C     when EXCH are not disabled (e.g. with EXCH1) between tile Edges
C     that are closed by OB.
C     Do it over OLx,OLy grid points beyond OB, in agreement with OBCS code
C     (apply_tracer) which overwrites tracer over the same width.
      OB_ApplX = OLx
      OB_ApplY = OLy

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)

C--   Set Interior mask to zero beyond OB
C-    Eastern boundary
        DO j=1-OLy,sNy+OLy
         jm = MAX( j-1, 1-OLy )
         iB = OB_Ie(j,bi,bj)
         IF ( iB.NE.OB_indexNone ) THEN
           DO i=iB,iB+OB_ApplX-1
             OBCS_insideMask(i,j,bi,bj) = 0.
           ENDDO
           DO i=iB+1,iB+OB_ApplX-1
             maskInW(i,j,bi,bj) = 0.
           ENDDO
           IF ( OB_Ie(jm,bi,bj).NE.OB_indexNone ) THEN
            iB = MAX( iB, OB_Ie(jm,bi,bj) )
            DO i=iB,iB+OB_ApplX-1
             maskInS(i,j,bi,bj) = 0.
            ENDDO
           ENDIF
         ENDIF
        ENDDO
C-    Western boundary
        DO j=1-OLy,sNy+OLy
         jm = MAX( j-1, 1-OLy )
         iB = OB_Iw(j,bi,bj)
         IF ( iB.NE.OB_indexNone ) THEN
           DO i=1-OB_ApplX+iB,iB
             OBCS_insideMask(i,j,bi,bj) = 0.
           ENDDO
           DO i=2-OB_ApplX+iB,iB
             maskInW(i,j,bi,bj) = 0.
           ENDDO
           IF ( OB_Iw(jm,bi,bj).NE.OB_indexNone ) THEN
            iB = MIN( iB, OB_Iw(jm,bi,bj) )
            DO i=1-OB_ApplX+iB,iB
             maskInS(i,j,bi,bj) = 0.
            ENDDO
           ENDIF
         ENDIF
        ENDDO
C-    Northern boundary
        DO i=1-OLx,sNx+OLx
         im = MAX( i-1, 1-OLx )
         jB = OB_Jn(i,bi,bj)
         IF ( jB.NE.OB_indexNone ) THEN
           DO j=jB,jB+OB_ApplY-1
             OBCS_insideMask(i,j,bi,bj) = 0.
           ENDDO
           DO j=jB+1,jB+OB_ApplY-1
             maskInS(i,j,bi,bj) = 0.
           ENDDO
           IF ( OB_Jn(im,bi,bj).NE.OB_indexNone ) THEN
            jB = MAX( jB, OB_Jn(im,bi,bj) )
            DO j=jB,jB+OB_ApplY-1
             maskInW(i,j,bi,bj) = 0.
            ENDDO
           ENDIF
         ENDIF
        ENDDO
C-    Southern boundary
        DO i=1-OLx,sNx+OLx
         im = MAX( i-1, 1-OLx )
         jB = OB_Js(i,bi,bj)
         IF ( jB.NE.OB_indexNone ) THEN
           DO j=1-OB_ApplY+jB,jB
             OBCS_insideMask(i,j,bi,bj) = 0.
           ENDDO
           DO j=2-OB_ApplY+jB,jB
             maskInS(i,j,bi,bj) = 0.
           ENDDO
           IF ( OB_Js(im,bi,bj).NE.OB_indexNone ) THEN
            jB = MIN( jB, OB_Js(im,bi,bj) )
            DO j=1-OB_ApplY+jB,jB
             maskInW(i,j,bi,bj) = 0.
            ENDDO
           ENDIF
         ENDIF
        ENDDO

C--   Apply mask to maskInC :
        DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
            maskInC(i,j,bi,bj) = maskInC(i,j,bi,bj)
     &                          *OBCS_insideMask(i,j,bi,bj)
          ENDDO
        ENDDO

C--   end bi,bj loops
       ENDDO
      ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   Set OB active tiles:
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
         tileHasOBE(bi,bj) = .FALSE.
         tileHasOBW(bi,bj) = .FALSE.
         tileHasOBN(bi,bj) = .FALSE.
         tileHasOBS(bi,bj) = .FALSE.
         DO j=1-OLy,sNy+OLy
          tileHasOBE(bi,bj) = tileHasOBE(bi,bj) .OR.
     &                      ( OB_Ie(j,bi,bj).NE.OB_indexNone )
          tileHasOBW(bi,bj) = tileHasOBW(bi,bj) .OR.
     &                      ( OB_Iw(j,bi,bj).NE.OB_indexNone )
         ENDDO
         DO i=1-OLx,sNx+OLx
          tileHasOBN(bi,bj) = tileHasOBN(bi,bj) .OR.
     &                      ( OB_Jn(i,bi,bj).NE.OB_indexNone )
          tileHasOBS(bi,bj) = tileHasOBS(bi,bj) .OR.
     &                      ( OB_Js(i,bi,bj).NE.OB_indexNone )
         ENDDO
       ENDDO
      ENDDO

C--   Set domain connected-piece Id for OB grid points:
      CALL OBCS_SET_CONNECT( myThid )

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_LEAVE('OBCS_INIT_FIXED',myThid)
#endif

#endif /* ALLOW_OBCS */
      RETURN
      END