model/src/ini_masks_etc.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_masks_etc.F,v 1.57 2017/04/11 18:07:15 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_MASKS_ETC
C     !INTERFACE:
      SUBROUTINE INI_MASKS_ETC( myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_MASKS_ETC
C     | o Initialise masks and topography factors
C     *==========================================================*
C     | These arrays are used throughout the code and describe
C     | the topography of the domain through masks (0s and 1s)
C     | and fractional height factors (0<hFac<1). The latter
C     | distinguish between the lopped-cell and full-step
C     | topographic representations.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef NONLIN_FRSURF
# include "SURFACE.h"
#endif /* NONLIN_FRSURF */

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

C     !LOCAL VARIABLES:
C     bi, bj    :: tile indices
C     i, j, k   :: Loop counters
C     tmpFld    :: Temporary array used to compute & write Total Depth
C     tmpVar*   :: Temporary array used to integrate column thickness
      INTEGER bi, bj
      INTEGER i, j, k
      _RS tmpFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL tmpVar1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL tmpVar2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL hFacMnSz, hFacCtmp
      _RL hFac1tmp, hFac2tmp
CEOP

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

#ifdef ALLOW_SHELFICE
      IF ( useShelfIce ) THEN
C--   Modify  ocean upper boundary position according to ice-shelf topography
       CALL SHELFICE_INIT_DEPTHS(
     U               R_low, Ro_surf,
     I               myThid )
      ENDIF
#endif /* ALLOW_SHELFICE */

      IF ( selectSigmaCoord.EQ.0 ) THEN
C---  r-coordinate with partial-cell or full cell mask

C--   Initialise rLow & reference rSurf at Western & Southern edges (U & V pts)
C     Note: not final value since these estimates ignore hFacMin constrain
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        i = 1-OLx
        DO j=1-OLy,sNy+OLy
           rLowW (i,j,bi,bj) = rF(1)
           rSurfW(i,j,bi,bj) = rF(1)
        ENDDO
        j = 1-OLy
        DO i=1-OLx,sNx+OLx
           rLowS (i,j,bi,bj) = rF(1)
           rSurfS(i,j,bi,bj) = rF(1)
        ENDDO
        DO j=1-OLy,sNy+OLy
         DO i=2-OLx,sNx+OLx
           rLowW(i,j,bi,bj)  =
     &           MAX(   R_low(i-1,j,bi,bj),   R_low(i,j,bi,bj) )
           rSurfW(i,j,bi,bj) =
     &           MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
           rSurfW(i,j,bi,bj) =
     &           MAX( rSurfW(i,j,bi,bj), rLowW(i,j,bi,bj) )
         ENDDO
        ENDDO
        DO j=2-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
           rLowS(i,j,bi,bj)  =
     &           MAX(   R_low(i,j-1,bi,bj),   R_low(i,j,bi,bj) )
           rSurfS(i,j,bi,bj) =
     &           MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
           rSurfS(i,j,bi,bj) =
     &           MAX( rSurfS(i,j,bi,bj), rLowS(i,j,bi,bj) )
         ENDDO
        ENDDO
       ENDDO
      ENDDO

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

C--   Calculate lopping factor hFacC : over-estimate the part inside of the domain
C     taking into account the lower_R Boundary (Bathymetry / Top of Atmos)
        DO k=1, Nr
         hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
C      o Non-dimensional distance between grid bound. and domain lower_R bound.
           hFacCtmp = (rF(k)-R_low(i,j,bi,bj))*recip_drF(k)
C      o Select between, closed, open or partial (0,1,0-1)
           hFacCtmp = MIN( MAX( hFacCtmp, zeroRL ) , oneRL )
C      o Impose minimum fraction and/or size (dimensional)
           IF ( hFacCtmp.LT.hFacMnSz ) THEN
            IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
             hFacC(i,j,k,bi,bj) = 0.
            ELSE
             hFacC(i,j,k,bi,bj) = hFacMnSz
            ENDIF
           ELSE
             hFacC(i,j,k,bi,bj) = hFacCtmp
           ENDIF
          ENDDO
         ENDDO
        ENDDO

C-    Re-calculate lower-R Boundary position, taking into account hFacC
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
           tmpVar1(i,j) = 0.
         ENDDO
        ENDDO
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           tmpVar1(i,j) = tmpVar1(i,j) + drF(k)*hFacC(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
           R_low(i,j,bi,bj) = rF(1) - tmpVar1(i,j)
         ENDDO
        ENDDO

C--   Calculate lopping factor hFacC : Remove part outside of the domain
C     taking into account the Reference (=at rest) Surface Position Ro_surf
        DO k=1, Nr
         hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
C      o Non-dimensional distance between grid boundary and model surface
           hFacCtmp = (rF(k)-Ro_surf(i,j,bi,bj))*recip_drF(k)
C      o Reduce the previous fraction : substract the outside part.
           hFacCtmp = hFacC(i,j,k,bi,bj) - MAX( hFacCtmp, zeroRL )
C      o set to zero if empty Column :
           hFacCtmp = MAX( hFacCtmp, zeroRL )
C      o Impose minimum fraction and/or size (dimensional)
           IF ( hFacCtmp.LT.hFacMnSz ) THEN
            IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
             hFacC(i,j,k,bi,bj) = 0.
            ELSE
             hFacC(i,j,k,bi,bj) = hFacMnSz
            ENDIF
           ELSE
             hFacC(i,j,k,bi,bj) = hFacCtmp
           ENDIF
          ENDDO
         ENDDO
        ENDDO

C-    Re-calculate Reference surface position, taking into account hFacC
C     initialize Total column fluid thickness and surface k index
C       Note: if no fluid (continent) ==> kSurf = Nr+1
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          tmpVar2(i,j) = 0.
          tmpFld(i,j,bi,bj) = 0.
          kSurfC(i,j,bi,bj) = Nr+1
          kLowC (i,j,bi,bj) = 0
         ENDDO
        ENDDO
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           tmpVar2(i,j) = tmpVar2(i,j) + drF(k)*hFacC(i,j,k,bi,bj)
           tmpFld(i,j,bi,bj) = tmpFld(i,j,bi,bj) + 1.
           IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) kLowC(i,j,bi,bj) = k
          ENDDO
         ENDDO
        ENDDO
        DO k=Nr,1,-1
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) kSurfC(i,j,bi,bj) = k
          ENDDO
         ENDDO
        ENDDO
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          Ro_surf(i,j,bi,bj) = R_low(i,j,bi,bj) + tmpVar2(i,j)
          maskInC(i,j,bi,bj) = 0.
          IF ( kSurfC(i,j,bi,bj).LE.Nr ) maskInC(i,j,bi,bj) = 1.
         ENDDO
        ENDDO

C-    end bi,bj loops.
       ENDDO
      ENDDO

      IF ( plotLevel.GE.debLevB ) THEN
c       CALL PLOT_FIELD_XYRS( tmpFld,
c    &           'Model Depths K Index' , -1, myThid )
        CALL PLOT_FIELD_XYRS(R_low,
     &           'Model R_low (ini_masks_etc)', -1, myThid )
        CALL PLOT_FIELD_XYRS(Ro_surf,
     &           'Model Ro_surf (ini_masks_etc)', -1, myThid )
      ENDIF

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

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

C--   Calculate quantities derived from XY depth map
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
C         Total fluid column thickness (r_unit) :
          tmpVar1(i,j) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
          tmpFld(i,j,bi,bj) = tmpVar1(i,j)
C         Inverse of fluid column thickness (1/r_unit)
          IF ( tmpVar1(i,j) .LE. zeroRL ) THEN
           recip_Rcol(i,j,bi,bj) = 0.
          ELSE
           recip_Rcol(i,j,bi,bj) = 1. _d 0 / tmpVar1(i,j)
          ENDIF
         ENDDO
        ENDDO

C- Method-1 (useMin4hFacEdges = T):
C    compute hFacW,hFacS as minimum of adjacent hFacC factor
C- Method-2 (useMin4hFacEdges = F):
C    compute hFacW,hFacS from rSurfW,S and rLowW,S by applying
C    same rules as for hFacC
C Note: Currently, no difference between methods except when useShelfIce=T and
C       if, in adjacent columns, ice-draft and bathy are within the same level k

        IF ( useMin4hFacEdges ) THEN
C--   hFacW and hFacS (at U and V points):
C-    Method-1: use simply minimum of adjacent hFacC factor

         DO k=1, Nr
          DO j=1-OLy,sNy+OLy
           hFacW(1-OLx,j,k,bi,bj) = 0.
           DO i=2-OLx,sNx+OLx
            hFacW(i,j,k,bi,bj) =
     &        MIN( hFacC(i,j,k,bi,bj), hFacC(i-1,j,k,bi,bj) )
           ENDDO
          ENDDO
          DO i=1-OLx,sNx+OLx
            hFacS(i,1-OLy,k,bi,bj) = 0.
          ENDDO
          DO j=2-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            hFacS(i,j,k,bi,bj) =
     &        MIN( hFacC(i,j,k,bi,bj), hFacC(i,j-1,k,bi,bj) )
           ENDDO
          ENDDO
         ENDDO

        ELSE
C--   hFacW and hFacS (at U and V points):
C-    Method-2: compute new hFacW,S from rSurfW,S and rLowW,S
C               by applying same rules as for hFacC

         DO k=1, Nr
          hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
C      o Non-dimensional distance between grid bound. and domain lower_R bound.
            hFac1tmp = ( rF(k) - rLowW(i,j,bi,bj) )*recip_drF(k)
            hFacCtmp = MIN( hFac1tmp, oneRL )
c           hFacCtmp = MAX( hFacCtmp, zeroRL )
C      o Impose minimum fraction and/or size (dimensional)
            IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
              hFac1tmp = 0.
            ELSE
              hFac1tmp = MAX( hFacCtmp, hFacMnSz )
            ENDIF
C      o Reduce the previous fraction : substract the outside fraction
C        (i.e., beyond reference (=at rest) surface position rSurfW)
            hFac2tmp = ( rF(k) -rSurfW(i,j,bi,bj) )*recip_drF(k)
            hFacCtmp = hFac1tmp - MAX( hFac2tmp, zeroRL )
C      o Impose minimum fraction and/or size (dimensional)
            IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
              hFacW(i,j,k,bi,bj) = 0.
            ELSE
              hFacW(i,j,k,bi,bj) = MAX( hFacCtmp, hFacMnSz )
            ENDIF
           ENDDO
          ENDDO
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
C      o Non-dimensional distance between grid bound. and domain lower_R bound.
            hFac1tmp = ( rF(k) - rLowS(i,j,bi,bj) )*recip_drF(k)
            hFacCtmp = MIN( hFac1tmp, oneRL )
c           hFacCtmp = MAX( hFacCtmp, zeroRL )
C      o Impose minimum fraction and/or size (dimensional)
            IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
              hFac1tmp = 0.
            ELSE
              hFac1tmp = MAX( hFacCtmp, hFacMnSz )
            ENDIF
C      o Reduce the previous fraction : substract the outside fraction
C        (i.e., beyond reference (=at rest) surface position rSurfS)
            hFac2tmp = ( rF(k) -rSurfS(i,j,bi,bj) )*recip_drF(k)
            hFacCtmp = hFac1tmp - MAX( hFac2tmp, zeroRL )
C      o Impose minimum fraction and/or size (dimensional)
            IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
              hFacS(i,j,k,bi,bj) = 0.
            ELSE
              hFacS(i,j,k,bi,bj) = MAX( hFacCtmp, hFacMnSz )
            ENDIF
           ENDDO
          ENDDO
         ENDDO
        ENDIF

C--   Update rLow & reference rSurf at Western & Southern edges (U & V pts):
C     account for adjusted R_low & Ro_surf due to hFacMin constrain on hFacC.
C     Might need further adjustment (e.g., if useShelfIce=T) to match
C     integrated level thickness ( =Sum_k(drF*hFac) )
        DO j=1-OLy,sNy+OLy
         DO i=2-OLx,sNx+OLx
           rLowW(i,j,bi,bj)  =
     &           MAX(   R_low(i-1,j,bi,bj),   R_low(i,j,bi,bj) )
           rSurfW(i,j,bi,bj) =
     &           MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
           rSurfW(i,j,bi,bj) =
     &           MAX( rSurfW(i,j,bi,bj), rLowW(i,j,bi,bj) )
         ENDDO
        ENDDO
        DO j=2-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
           rLowS(i,j,bi,bj)  =
     &           MAX(   R_low(i,j-1,bi,bj),   R_low(i,j,bi,bj) )
           rSurfS(i,j,bi,bj) =
     &           MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
           rSurfS(i,j,bi,bj) =
     &           MAX( rSurfS(i,j,bi,bj), rLowS(i,j,bi,bj) )
         ENDDO
        ENDDO

c       IF ( useShelfIce ) THEN
C--   Adjust rLow & reference rSurf at Western & Southern edges (U & V pts)
C     to get consistent column thickness from Sum_k(hFac*drF) and rSurf-rLow

C-    Total column thickness at Western & Southern edges
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
            tmpVar1(i,j) = 0. _d 0
            tmpVar2(i,j) = 0. _d 0
          ENDDO
         ENDDO
         DO k=1,Nr
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            tmpVar1(i,j) = tmpVar1(i,j) + drF(k)*hFacW(i,j,k,bi,bj)
            tmpVar2(i,j) = tmpVar2(i,j) + drF(k)*hFacS(i,j,k,bi,bj)
           ENDDO
          ENDDO
         ENDDO

         IF ( useMin4hFacEdges ) THEN
C-    Adjust only rSurf at W and S edges (correct for the difference)
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             rSurfW(i,j,bi,bj) = rLowW(i,j,bi,bj) + tmpVar1(i,j)
             rSurfS(i,j,bi,bj) = rLowS(i,j,bi,bj) + tmpVar2(i,j)
            ENDDO
           ENDDO
         ELSE
C-    Adjust both rLow and rSurf at W & S edges (split correction by half)
C     adjust rSurfW and rLowW:
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             tmpVar1(i,j) = rLowW(i,j,bi,bj) + tmpVar1(i,j)
             tmpVar1(i,j) = ( tmpVar1(i,j) -rSurfW(i,j,bi,bj) )*halfRL
            ENDDO
           ENDDO
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             rSurfW(i,j,bi,bj) = rSurfW(i,j,bi,bj) + tmpVar1(i,j)
             rLowW (i,j,bi,bj) = rLowW (i,j,bi,bj) - tmpVar1(i,j)
            ENDDO
           ENDDO
C     Adjust rSurfS and rLowS:
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             tmpVar2(i,j) = rLowS(i,j,bi,bj) + tmpVar2(i,j)
             tmpVar2(i,j) = ( tmpVar2(i,j) -rSurfS(i,j,bi,bj) )*halfRL
            ENDDO
           ENDDO
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             rSurfS(i,j,bi,bj) = rSurfS(i,j,bi,bj) + tmpVar2(i,j)
             rLowS (i,j,bi,bj) = rLowS (i,j,bi,bj) - tmpVar2(i,j)
            ENDDO
           ENDDO
         ENDIF

C-    end if useShelfIce
c       ENDIF

C-    end bi,bj loops.
       ENDDO
      ENDDO

      CALL EXCH_UV_XYZ_RS( hFacW,  hFacS, .FALSE., myThid )
      CALL EXCH_UV_XY_RS( rSurfW, rSurfS, .FALSE., myThid )
      CALL EXCH_UV_XY_RS(  rLowW,  rLowS, .FALSE., myThid )

C--   Additional closing of Western and Southern grid-cell edges: for example,
C     a) might add some "thin walls" in specific location
C--   b) close non-periodic N & S boundaries of lat-lon grid at the N/S poles.
      CALL ADD_WALLS2MASKS( myThid )

C--   Calculate surface k index for interface W & S (U & V points)
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          kSurfW(i,j,bi,bj) = Nr+1
          kSurfS(i,j,bi,bj) = Nr+1
          DO k=Nr,1,-1
           IF (hFacW(i,j,k,bi,bj).NE.zeroRS) kSurfW(i,j,bi,bj) = k
           IF (hFacS(i,j,k,bi,bj).NE.zeroRS) kSurfS(i,j,bi,bj) = k
          ENDDO
          maskInW(i,j,bi,bj)= 0.
          IF ( kSurfW(i,j,bi,bj).LE.Nr ) maskInW(i,j,bi,bj)= 1.
          maskInS(i,j,bi,bj)= 0.
          IF ( kSurfS(i,j,bi,bj).LE.Nr ) maskInS(i,j,bi,bj)= 1.
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      ELSE
#ifndef DISABLE_SIGMA_CODE
C---  Sigma and Hybrid-Sigma set-up:
        CALL INI_SIGMA_HFAC( myThid )
#endif /* DISABLE_SIGMA_CODE */
      ENDIF

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

C--   Write to disk: Total Column Thickness & hFac(C,W,S):
C     This I/O is now done in write_grid.F
c     CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpFld,0,myThid)
c     CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
c     CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
c     CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)

      IF ( plotLevel.GE.debLevB ) THEN
        CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 0, myThid )
        CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 0, myThid )
        CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 0, myThid )
      ENDIF

C--   Masks and reciprocals of hFac[CWS]
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) THEN
            recip_hFacC(i,j,k,bi,bj) = 1. _d 0 / hFacC(i,j,k,bi,bj)
            maskC(i,j,k,bi,bj) = 1.
           ELSE
            recip_hFacC(i,j,k,bi,bj) = 0.
            maskC(i,j,k,bi,bj) = 0.
           ENDIF
           IF ( hFacW(i,j,k,bi,bj).NE.zeroRS ) THEN
            recip_hFacW(i,j,k,bi,bj) = 1. _d 0 / hFacW(i,j,k,bi,bj)
            maskW(i,j,k,bi,bj) = 1.
           ELSE
            recip_hFacW(i,j,k,bi,bj) = 0.
            maskW(i,j,k,bi,bj) = 0.
           ENDIF
           IF ( hFacS(i,j,k,bi,bj).NE.zeroRS ) THEN
            recip_hFacS(i,j,k,bi,bj) = 1. _d 0 / hFacS(i,j,k,bi,bj)
            maskS(i,j,k,bi,bj) = 1.
           ELSE
            recip_hFacS(i,j,k,bi,bj) = 0.
            maskS(i,j,k,bi,bj) = 0.
           ENDIF
          ENDDO
         ENDDO
        ENDDO
#ifdef NONLIN_FRSURF
C--   Save initial geometrical hFac factor into h0Fac (fixed in time):
C     Note: In case 1 pkg modifies hFac (from packages_init_fixed, called
C     later in sequence of calls) this pkg would need also to update h0Fac.
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           h0FacC(i,j,k,bi,bj) = _hFacC(i,j,k,bi,bj)
           h0FacW(i,j,k,bi,bj) = _hFacW(i,j,k,bi,bj)
           h0FacS(i,j,k,bi,bj) = _hFacS(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO
#endif /* NONLIN_FRSURF */
C-    end bi,bj loops.
       ENDDO
      ENDDO

c #ifdef ALLOW_NONHYDROSTATIC
C--   Calculate "recip_hFacU" = reciprocal hfac distance/volume for W cells
C NOTE:  not used ; computed locally in CALC_GW
c #endif

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/model/src/ini_masks_etc.F,v 1.57 2017/04/11 18:07:15 jmc Exp $
2	C $Name: $
3
4	#include "PACKAGES_CONFIG.h"
5	#include "CPP_OPTIONS.h"
6
7	CBOP
8	C !ROUTINE: INI_MASKS_ETC
9	C !INTERFACE:
10	SUBROUTINE INI_MASKS_ETC( myThid )
11	C !DESCRIPTION: \bv
12	C ==========================================================
13	C \| SUBROUTINE INI_MASKS_ETC
14	C \| o Initialise masks and topography factors
15	C ==========================================================
16	C \| These arrays are used throughout the code and describe
17	C \| the topography of the domain through masks (0s and 1s)
18	C \| and fractional height factors (0<hFac<1). The latter
19	C \| distinguish between the lopped-cell and full-step
20	C \| topographic representations.
21	C ==========================================================
22	C \ev
23
24	C !USES:
25	IMPLICIT NONE
26	C === Global variables ===
27	#include "SIZE.h"
28	#include "EEPARAMS.h"
29	#include "PARAMS.h"
30	#include "GRID.h"
31	#ifdef NONLIN_FRSURF
32	# include "SURFACE.h"
33	#endif /* NONLIN_FRSURF */
34
35	C !INPUT/OUTPUT PARAMETERS:
36	C myThid :: my Thread Id number
37	INTEGER myThid
38
39	C !LOCAL VARIABLES:
40	C bi, bj :: tile indices
41	C i, j, k :: Loop counters
42	C tmpFld :: Temporary array used to compute & write Total Depth
43	C tmpVar* :: Temporary array used to integrate column thickness
44	INTEGER bi, bj
45	INTEGER i, j, k
46	_RS tmpFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
47	_RL tmpVar1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48	_RL tmpVar2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49	_RL hFacMnSz, hFacCtmp
50	_RL hFac1tmp, hFac2tmp
51	CEOP
52
53	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
54
55	#ifdef ALLOW_SHELFICE
56	IF ( useShelfIce ) THEN
57	C-- Modify ocean upper boundary position according to ice-shelf topography
58	CALL SHELFICE_INIT_DEPTHS(
59	U R_low, Ro_surf,
60	I myThid )
61	ENDIF
62	#endif /* ALLOW_SHELFICE */
63
64	IF ( selectSigmaCoord.EQ.0 ) THEN
65	C--- r-coordinate with partial-cell or full cell mask
66
67	C-- Initialise rLow & reference rSurf at Western & Southern edges (U & V pts)
68	C Note: not final value since these estimates ignore hFacMin constrain
69	DO bj=myByLo(myThid), myByHi(myThid)
70	DO bi=myBxLo(myThid), myBxHi(myThid)
71	i = 1-OLx
72	DO j=1-OLy,sNy+OLy
73	rLowW (i,j,bi,bj) = rF(1)
74	rSurfW(i,j,bi,bj) = rF(1)
75	ENDDO
76	j = 1-OLy
77	DO i=1-OLx,sNx+OLx
78	rLowS (i,j,bi,bj) = rF(1)
79	rSurfS(i,j,bi,bj) = rF(1)
80	ENDDO
81	DO j=1-OLy,sNy+OLy
82	DO i=2-OLx,sNx+OLx
83	rLowW(i,j,bi,bj) =
84	& MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
85	rSurfW(i,j,bi,bj) =
86	& MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
87	rSurfW(i,j,bi,bj) =
88	& MAX( rSurfW(i,j,bi,bj), rLowW(i,j,bi,bj) )
89	ENDDO
90	ENDDO
91	DO j=2-OLy,sNy+OLy
92	DO i=1-OLx,sNx+OLx
93	rLowS(i,j,bi,bj) =
94	& MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
95	rSurfS(i,j,bi,bj) =
96	& MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
97	rSurfS(i,j,bi,bj) =
98	& MAX( rSurfS(i,j,bi,bj), rLowS(i,j,bi,bj) )
99	ENDDO
100	ENDDO
101	ENDDO
102	ENDDO
103
104	DO bj=myByLo(myThid), myByHi(myThid)
105	DO bi=myBxLo(myThid), myBxHi(myThid)
106
107	C-- Calculate lopping factor hFacC : over-estimate the part inside of the domain
108	C taking into account the lower_R Boundary (Bathymetry / Top of Atmos)
109	DO k=1, Nr
110	hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
111	DO j=1-OLy,sNy+OLy
112	DO i=1-OLx,sNx+OLx
113	C o Non-dimensional distance between grid bound. and domain lower_R bound.
114	hFacCtmp = (rF(k)-R_low(i,j,bi,bj))*recip_drF(k)
115	C o Select between, closed, open or partial (0,1,0-1)
116	hFacCtmp = MIN( MAX( hFacCtmp, zeroRL ) , oneRL )
117	C o Impose minimum fraction and/or size (dimensional)
118	IF ( hFacCtmp.LT.hFacMnSz ) THEN
119	IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
120	hFacC(i,j,k,bi,bj) = 0.
121	ELSE
122	hFacC(i,j,k,bi,bj) = hFacMnSz
123	ENDIF
124	ELSE
125	hFacC(i,j,k,bi,bj) = hFacCtmp
126	ENDIF
127	ENDDO
128	ENDDO
129	ENDDO
130
131	C- Re-calculate lower-R Boundary position, taking into account hFacC
132	DO j=1-OLy,sNy+OLy
133	DO i=1-OLx,sNx+OLx
134	tmpVar1(i,j) = 0.
135	ENDDO
136	ENDDO
137	DO k=1,Nr
138	DO j=1-OLy,sNy+OLy
139	DO i=1-OLx,sNx+OLx
140	tmpVar1(i,j) = tmpVar1(i,j) + drF(k)*hFacC(i,j,k,bi,bj)
141	ENDDO
142	ENDDO
143	ENDDO
144	DO j=1-OLy,sNy+OLy
145	DO i=1-OLx,sNx+OLx
146	R_low(i,j,bi,bj) = rF(1) - tmpVar1(i,j)
147	ENDDO
148	ENDDO
149
150	C-- Calculate lopping factor hFacC : Remove part outside of the domain
151	C taking into account the Reference (=at rest) Surface Position Ro_surf
152	DO k=1, Nr
153	hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
154	DO j=1-OLy,sNy+OLy
155	DO i=1-OLx,sNx+OLx
156	C o Non-dimensional distance between grid boundary and model surface
157	hFacCtmp = (rF(k)-Ro_surf(i,j,bi,bj))*recip_drF(k)
158	C o Reduce the previous fraction : substract the outside part.
159	hFacCtmp = hFacC(i,j,k,bi,bj) - MAX( hFacCtmp, zeroRL )
160	C o set to zero if empty Column :
161	hFacCtmp = MAX( hFacCtmp, zeroRL )
162	C o Impose minimum fraction and/or size (dimensional)
163	IF ( hFacCtmp.LT.hFacMnSz ) THEN
164	IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
165	hFacC(i,j,k,bi,bj) = 0.
166	ELSE
167	hFacC(i,j,k,bi,bj) = hFacMnSz
168	ENDIF
169	ELSE
170	hFacC(i,j,k,bi,bj) = hFacCtmp
171	ENDIF
172	ENDDO
173	ENDDO
174	ENDDO
175
176	C- Re-calculate Reference surface position, taking into account hFacC
177	C initialize Total column fluid thickness and surface k index
178	C Note: if no fluid (continent) ==> kSurf = Nr+1
179	DO j=1-OLy,sNy+OLy
180	DO i=1-OLx,sNx+OLx
181	tmpVar2(i,j) = 0.
182	tmpFld(i,j,bi,bj) = 0.
183	kSurfC(i,j,bi,bj) = Nr+1
184	kLowC (i,j,bi,bj) = 0
185	ENDDO
186	ENDDO
187	DO k=1,Nr
188	DO j=1-OLy,sNy+OLy
189	DO i=1-OLx,sNx+OLx
190	tmpVar2(i,j) = tmpVar2(i,j) + drF(k)*hFacC(i,j,k,bi,bj)
191	tmpFld(i,j,bi,bj) = tmpFld(i,j,bi,bj) + 1.
192	IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) kLowC(i,j,bi,bj) = k
193	ENDDO
194	ENDDO
195	ENDDO
196	DO k=Nr,1,-1
197	DO j=1-OLy,sNy+OLy
198	DO i=1-OLx,sNx+OLx
199	IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) kSurfC(i,j,bi,bj) = k
200	ENDDO
201	ENDDO
202	ENDDO
203	DO j=1-OLy,sNy+OLy
204	DO i=1-OLx,sNx+OLx
205	Ro_surf(i,j,bi,bj) = R_low(i,j,bi,bj) + tmpVar2(i,j)
206	maskInC(i,j,bi,bj) = 0.
207	IF ( kSurfC(i,j,bi,bj).LE.Nr ) maskInC(i,j,bi,bj) = 1.
208	ENDDO
209	ENDDO
210
211	C- end bi,bj loops.
212	ENDDO
213	ENDDO
214
215	IF ( plotLevel.GE.debLevB ) THEN
216	c CALL PLOT_FIELD_XYRS( tmpFld,
217	c & 'Model Depths K Index' , -1, myThid )
218	CALL PLOT_FIELD_XYRS(R_low,
219	& 'Model R_low (ini_masks_etc)', -1, myThid )
220	CALL PLOT_FIELD_XYRS(Ro_surf,
221	& 'Model Ro_surf (ini_masks_etc)', -1, myThid )
222	ENDIF
223
224	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
225
226	DO bj = myByLo(myThid), myByHi(myThid)
227	DO bi = myBxLo(myThid), myBxHi(myThid)
228
229	C-- Calculate quantities derived from XY depth map
230	DO j=1-OLy,sNy+OLy
231	DO i=1-OLx,sNx+OLx
232	C Total fluid column thickness (r_unit) :
233	tmpVar1(i,j) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
234	tmpFld(i,j,bi,bj) = tmpVar1(i,j)
235	C Inverse of fluid column thickness (1/r_unit)
236	IF ( tmpVar1(i,j) .LE. zeroRL ) THEN
237	recip_Rcol(i,j,bi,bj) = 0.
238	ELSE
239	recip_Rcol(i,j,bi,bj) = 1. _d 0 / tmpVar1(i,j)
240	ENDIF
241	ENDDO
242	ENDDO
243
244	C- Method-1 (useMin4hFacEdges = T):
245	C compute hFacW,hFacS as minimum of adjacent hFacC factor
246	C- Method-2 (useMin4hFacEdges = F):
247	C compute hFacW,hFacS from rSurfW,S and rLowW,S by applying
248	C same rules as for hFacC
249	C Note: Currently, no difference between methods except when useShelfIce=T and
250	C if, in adjacent columns, ice-draft and bathy are within the same level k
251
252	IF ( useMin4hFacEdges ) THEN
253	C-- hFacW and hFacS (at U and V points):
254	C- Method-1: use simply minimum of adjacent hFacC factor
255
256	DO k=1, Nr
257	DO j=1-OLy,sNy+OLy
258	hFacW(1-OLx,j,k,bi,bj) = 0.
259	DO i=2-OLx,sNx+OLx
260	hFacW(i,j,k,bi,bj) =
261	& MIN( hFacC(i,j,k,bi,bj), hFacC(i-1,j,k,bi,bj) )
262	ENDDO
263	ENDDO
264	DO i=1-OLx,sNx+OLx
265	hFacS(i,1-OLy,k,bi,bj) = 0.
266	ENDDO
267	DO j=2-OLy,sNy+OLy
268	DO i=1-OLx,sNx+OLx
269	hFacS(i,j,k,bi,bj) =
270	& MIN( hFacC(i,j,k,bi,bj), hFacC(i,j-1,k,bi,bj) )
271	ENDDO
272	ENDDO
273	ENDDO
274
275	ELSE
276	C-- hFacW and hFacS (at U and V points):
277	C- Method-2: compute new hFacW,S from rSurfW,S and rLowW,S
278	C by applying same rules as for hFacC
279
280	DO k=1, Nr
281	hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
282	DO j=1-OLy,sNy+OLy
283	DO i=1-OLx,sNx+OLx
284	C o Non-dimensional distance between grid bound. and domain lower_R bound.
285	hFac1tmp = ( rF(k) - rLowW(i,j,bi,bj) )*recip_drF(k)
286	hFacCtmp = MIN( hFac1tmp, oneRL )
287	c hFacCtmp = MAX( hFacCtmp, zeroRL )
288	C o Impose minimum fraction and/or size (dimensional)
289	IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
290	hFac1tmp = 0.
291	ELSE
292	hFac1tmp = MAX( hFacCtmp, hFacMnSz )
293	ENDIF
294	C o Reduce the previous fraction : substract the outside fraction
295	C (i.e., beyond reference (=at rest) surface position rSurfW)
296	hFac2tmp = ( rF(k) -rSurfW(i,j,bi,bj) )*recip_drF(k)
297	hFacCtmp = hFac1tmp - MAX( hFac2tmp, zeroRL )
298	C o Impose minimum fraction and/or size (dimensional)
299	IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
300	hFacW(i,j,k,bi,bj) = 0.
301	ELSE
302	hFacW(i,j,k,bi,bj) = MAX( hFacCtmp, hFacMnSz )
303	ENDIF
304	ENDDO
305	ENDDO
306	DO j=1-OLy,sNy+OLy
307	DO i=1-OLx,sNx+OLx
308	C o Non-dimensional distance between grid bound. and domain lower_R bound.
309	hFac1tmp = ( rF(k) - rLowS(i,j,bi,bj) )*recip_drF(k)
310	hFacCtmp = MIN( hFac1tmp, oneRL )
311	c hFacCtmp = MAX( hFacCtmp, zeroRL )
312	C o Impose minimum fraction and/or size (dimensional)
313	IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
314	hFac1tmp = 0.
315	ELSE
316	hFac1tmp = MAX( hFacCtmp, hFacMnSz )
317	ENDIF
318	C o Reduce the previous fraction : substract the outside fraction
319	C (i.e., beyond reference (=at rest) surface position rSurfS)
320	hFac2tmp = ( rF(k) -rSurfS(i,j,bi,bj) )*recip_drF(k)
321	hFacCtmp = hFac1tmp - MAX( hFac2tmp, zeroRL )
322	C o Impose minimum fraction and/or size (dimensional)
323	IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
324	hFacS(i,j,k,bi,bj) = 0.
325	ELSE
326	hFacS(i,j,k,bi,bj) = MAX( hFacCtmp, hFacMnSz )
327	ENDIF
328	ENDDO
329	ENDDO
330	ENDDO
331	ENDIF
332
333	C-- Update rLow & reference rSurf at Western & Southern edges (U & V pts):
334	C account for adjusted R_low & Ro_surf due to hFacMin constrain on hFacC.
335	C Might need further adjustment (e.g., if useShelfIce=T) to match
336	C integrated level thickness ( =Sum_k(drF*hFac) )
337	DO j=1-OLy,sNy+OLy
338	DO i=2-OLx,sNx+OLx
339	rLowW(i,j,bi,bj) =
340	& MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
341	rSurfW(i,j,bi,bj) =
342	& MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
343	rSurfW(i,j,bi,bj) =
344	& MAX( rSurfW(i,j,bi,bj), rLowW(i,j,bi,bj) )
345	ENDDO
346	ENDDO
347	DO j=2-OLy,sNy+OLy
348	DO i=1-OLx,sNx+OLx
349	rLowS(i,j,bi,bj) =
350	& MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
351	rSurfS(i,j,bi,bj) =
352	& MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
353	rSurfS(i,j,bi,bj) =
354	& MAX( rSurfS(i,j,bi,bj), rLowS(i,j,bi,bj) )
355	ENDDO
356	ENDDO
357
358	c IF ( useShelfIce ) THEN
359	C-- Adjust rLow & reference rSurf at Western & Southern edges (U & V pts)
360	C to get consistent column thickness from Sum_k(hFac*drF) and rSurf-rLow
361
362	C- Total column thickness at Western & Southern edges
363	DO j=1-OLy,sNy+OLy
364	DO i=1-OLx,sNx+OLx
365	tmpVar1(i,j) = 0. _d 0
366	tmpVar2(i,j) = 0. _d 0
367	ENDDO
368	ENDDO
369	DO k=1,Nr
370	DO j=1-OLy,sNy+OLy
371	DO i=1-OLx,sNx+OLx
372	tmpVar1(i,j) = tmpVar1(i,j) + drF(k)*hFacW(i,j,k,bi,bj)
373	tmpVar2(i,j) = tmpVar2(i,j) + drF(k)*hFacS(i,j,k,bi,bj)
374	ENDDO
375	ENDDO
376	ENDDO
377
378	IF ( useMin4hFacEdges ) THEN
379	C- Adjust only rSurf at W and S edges (correct for the difference)
380	DO j=1-OLy,sNy+OLy
381	DO i=1-OLx,sNx+OLx
382	rSurfW(i,j,bi,bj) = rLowW(i,j,bi,bj) + tmpVar1(i,j)
383	rSurfS(i,j,bi,bj) = rLowS(i,j,bi,bj) + tmpVar2(i,j)
384	ENDDO
385	ENDDO
386	ELSE
387	C- Adjust both rLow and rSurf at W & S edges (split correction by half)
388	C adjust rSurfW and rLowW:
389	DO j=1-OLy,sNy+OLy
390	DO i=1-OLx,sNx+OLx
391	tmpVar1(i,j) = rLowW(i,j,bi,bj) + tmpVar1(i,j)
392	tmpVar1(i,j) = ( tmpVar1(i,j) -rSurfW(i,j,bi,bj) )*halfRL
393	ENDDO
394	ENDDO
395	DO j=1-OLy,sNy+OLy
396	DO i=1-OLx,sNx+OLx
397	rSurfW(i,j,bi,bj) = rSurfW(i,j,bi,bj) + tmpVar1(i,j)
398	rLowW (i,j,bi,bj) = rLowW (i,j,bi,bj) - tmpVar1(i,j)
399	ENDDO
400	ENDDO
401	C Adjust rSurfS and rLowS:
402	DO j=1-OLy,sNy+OLy
403	DO i=1-OLx,sNx+OLx
404	tmpVar2(i,j) = rLowS(i,j,bi,bj) + tmpVar2(i,j)
405	tmpVar2(i,j) = ( tmpVar2(i,j) -rSurfS(i,j,bi,bj) )*halfRL
406	ENDDO
407	ENDDO
408	DO j=1-OLy,sNy+OLy
409	DO i=1-OLx,sNx+OLx
410	rSurfS(i,j,bi,bj) = rSurfS(i,j,bi,bj) + tmpVar2(i,j)
411	rLowS (i,j,bi,bj) = rLowS (i,j,bi,bj) - tmpVar2(i,j)
412	ENDDO
413	ENDDO
414	ENDIF
415
416	C- end if useShelfIce
417	c ENDIF
418
419	C- end bi,bj loops.
420	ENDDO
421	ENDDO
422
423	CALL EXCH_UV_XYZ_RS( hFacW, hFacS, .FALSE., myThid )
424	CALL EXCH_UV_XY_RS( rSurfW, rSurfS, .FALSE., myThid )
425	CALL EXCH_UV_XY_RS( rLowW, rLowS, .FALSE., myThid )
426
427	C-- Additional closing of Western and Southern grid-cell edges: for example,
428	C a) might add some "thin walls" in specific location
429	C-- b) close non-periodic N & S boundaries of lat-lon grid at the N/S poles.
430	CALL ADD_WALLS2MASKS( myThid )
431
432	C-- Calculate surface k index for interface W & S (U & V points)
433	DO bj=myByLo(myThid), myByHi(myThid)
434	DO bi=myBxLo(myThid), myBxHi(myThid)
435	DO j=1-OLy,sNy+OLy
436	DO i=1-OLx,sNx+OLx
437	kSurfW(i,j,bi,bj) = Nr+1
438	kSurfS(i,j,bi,bj) = Nr+1
439	DO k=Nr,1,-1
440	IF (hFacW(i,j,k,bi,bj).NE.zeroRS) kSurfW(i,j,bi,bj) = k
441	IF (hFacS(i,j,k,bi,bj).NE.zeroRS) kSurfS(i,j,bi,bj) = k
442	ENDDO
443	maskInW(i,j,bi,bj)= 0.
444	IF ( kSurfW(i,j,bi,bj).LE.Nr ) maskInW(i,j,bi,bj)= 1.
445	maskInS(i,j,bi,bj)= 0.
446	IF ( kSurfS(i,j,bi,bj).LE.Nr ) maskInS(i,j,bi,bj)= 1.
447	ENDDO
448	ENDDO
449	ENDDO
450	ENDDO
451
452	ELSE
453	#ifndef DISABLE_SIGMA_CODE
454	C--- Sigma and Hybrid-Sigma set-up:
455	CALL INI_SIGMA_HFAC( myThid )
456	#endif /* DISABLE_SIGMA_CODE */
457	ENDIF
458
459	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
460
461	C-- Write to disk: Total Column Thickness & hFac(C,W,S):
462	C This I/O is now done in write_grid.F
463	c CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpFld,0,myThid)
464	c CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
465	c CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
466	c CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
467
468	IF ( plotLevel.GE.debLevB ) THEN
469	CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 0, myThid )
470	CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 0, myThid )
471	CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 0, myThid )
472	ENDIF
473
474	C-- Masks and reciprocals of hFac[CWS]
475	DO bj = myByLo(myThid), myByHi(myThid)
476	DO bi = myBxLo(myThid), myBxHi(myThid)
477	DO k=1,Nr
478	DO j=1-OLy,sNy+OLy
479	DO i=1-OLx,sNx+OLx
480	IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) THEN
481	recip_hFacC(i,j,k,bi,bj) = 1. _d 0 / hFacC(i,j,k,bi,bj)
482	maskC(i,j,k,bi,bj) = 1.
483	ELSE
484	recip_hFacC(i,j,k,bi,bj) = 0.
485	maskC(i,j,k,bi,bj) = 0.
486	ENDIF
487	IF ( hFacW(i,j,k,bi,bj).NE.zeroRS ) THEN
488	recip_hFacW(i,j,k,bi,bj) = 1. _d 0 / hFacW(i,j,k,bi,bj)
489	maskW(i,j,k,bi,bj) = 1.
490	ELSE
491	recip_hFacW(i,j,k,bi,bj) = 0.
492	maskW(i,j,k,bi,bj) = 0.
493	ENDIF
494	IF ( hFacS(i,j,k,bi,bj).NE.zeroRS ) THEN
495	recip_hFacS(i,j,k,bi,bj) = 1. _d 0 / hFacS(i,j,k,bi,bj)
496	maskS(i,j,k,bi,bj) = 1.
497	ELSE
498	recip_hFacS(i,j,k,bi,bj) = 0.
499	maskS(i,j,k,bi,bj) = 0.
500	ENDIF
501	ENDDO
502	ENDDO
503	ENDDO
504	#ifdef NONLIN_FRSURF
505	C-- Save initial geometrical hFac factor into h0Fac (fixed in time):
506	C Note: In case 1 pkg modifies hFac (from packages_init_fixed, called
507	C later in sequence of calls) this pkg would need also to update h0Fac.
508	DO k=1,Nr
509	DO j=1-OLy,sNy+OLy
510	DO i=1-OLx,sNx+OLx
511	h0FacC(i,j,k,bi,bj) = _hFacC(i,j,k,bi,bj)
512	h0FacW(i,j,k,bi,bj) = _hFacW(i,j,k,bi,bj)
513	h0FacS(i,j,k,bi,bj) = _hFacS(i,j,k,bi,bj)
514	ENDDO
515	ENDDO
516	ENDDO
517	#endif /* NONLIN_FRSURF */
518	C- end bi,bj loops.
519	ENDDO
520	ENDDO
521
522	c #ifdef ALLOW_NONHYDROSTATIC
523	C-- Calculate "recip_hFacU" = reciprocal hfac distance/volume for W cells
524	C NOTE: not used ; computed locally in CALC_GW
525	c #endif
526
527	RETURN
528	END