Back to home page

MITgcm
 
 

    

File indexing completed on 2024-05-25 05:10:18 UTC

view on githubraw file Latest commit 00f81e67 on 2024-05-24 21:00:12 UTC
6d54cf9ca1 Ed H*0001 #include "PACKAGES_CONFIG.h"
1dbaea09ee Chri*0002 #include "CPP_OPTIONS.h"
a74b31ee34 Alis*0003 
9366854e02 Chri*0004 CBOP
                0005 C     !ROUTINE: INI_MASKS_ETC
                0006 C     !INTERFACE:
a74b31ee34 Alis*0007       SUBROUTINE INI_MASKS_ETC( myThid )
9366854e02 Chri*0008 C     !DESCRIPTION: \bv
                0009 C     *==========================================================*
ab5a98a4ed Jean*0010 C     | SUBROUTINE INI_MASKS_ETC
                0011 C     | o Initialise masks and topography factors
9366854e02 Chri*0012 C     *==========================================================*
ab5a98a4ed Jean*0013 C     | These arrays are used throughout the code and describe
                0014 C     | the topography of the domain through masks (0s and 1s)
                0015 C     | and fractional height factors (0<hFac<1). The latter
                0016 C     | distinguish between the lopped-cell and full-step
                0017 C     | topographic representations.
9366854e02 Chri*0018 C     *==========================================================*
                0019 C     \ev
a74b31ee34 Alis*0020 
9366854e02 Chri*0021 C     !USES:
                0022       IMPLICIT NONE
a74b31ee34 Alis*0023 C     === Global variables ===
                0024 #include "SIZE.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027 #include "GRID.h"
f88bbe67c4 Jean*0028 #ifdef NONLIN_FRSURF
                0029 # include "SURFACE.h"
                0030 #endif /* NONLIN_FRSURF */
a74b31ee34 Alis*0031 
9366854e02 Chri*0032 C     !INPUT/OUTPUT PARAMETERS:
21ab794b7f Jean*0033 C     myThid    ::  my Thread Id number
a74b31ee34 Alis*0034       INTEGER myThid
                0035 
9366854e02 Chri*0036 C     !LOCAL VARIABLES:
21ab794b7f Jean*0037 C     bi, bj    :: tile indices
                0038 C     i, j, k   :: Loop counters
                0039 C     tmpFld    :: Temporary array used to compute & write Total Depth
dff94812d5 Jean*0040 C     tmpVar    :: Temporary array used to integrate column thickness
                0041 C     rEmpty    :: empty column r-position
a74b31ee34 Alis*0042       INTEGER bi, bj
86566ea478 Jean*0043       INTEGER i, j, k
7b03dab2b2 jm-c 0044 c     _RS tmpFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
dff94812d5 Jean*0045       _RL tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a2f3aac956 Jean*0046       _RL hFacMnSz, hFac_loc
21ab794b7f Jean*0047       _RL hFac1tmp, hFac2tmp
dff94812d5 Jean*0048       _RS rEmpty
9366854e02 Chri*0049 CEOP
fb481a83c2 Alis*0050 
86566ea478 Jean*0051 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0052 
ea3c970699 Jean*0053 #ifdef ALLOW_SHELFICE
                0054       IF ( useShelfIce ) THEN
                0055 C--   Modify  ocean upper boundary position according to ice-shelf topography
                0056        CALL SHELFICE_INIT_DEPTHS(
                0057      U               R_low, Ro_surf,
                0058      I               myThid )
                0059       ENDIF
                0060 #endif /* ALLOW_SHELFICE */
00f81e6785 Ou W*0061 #ifdef ALLOW_STEEP_ICECAVITY
                0062       IF ( useSTIC ) CALL STIC_INIT_DEPTHS( myThid )
                0063 #endif /* ALLOW_STEEP_ICECAVITY */
ea3c970699 Jean*0064 
f15994caab Jean*0065       IF ( selectSigmaCoord.EQ.0 ) THEN
                0066 C---  r-coordinate with partial-cell or full cell mask
                0067 
21ab794b7f Jean*0068 C--   Initialise rLow & reference rSurf at Western & Southern edges (U & V pts)
                0069 C     Note: not final value since these estimates ignore hFacMin constrain
dff94812d5 Jean*0070       rEmpty = rF(1)
21ab794b7f Jean*0071       DO bj=myByLo(myThid), myByHi(myThid)
                0072        DO bi=myBxLo(myThid), myBxHi(myThid)
                0073         i = 1-OLx
                0074         DO j=1-OLy,sNy+OLy
dff94812d5 Jean*0075            rLowW (i,j,bi,bj) = rEmpty
                0076            rSurfW(i,j,bi,bj) = rEmpty
21ab794b7f Jean*0077         ENDDO
                0078         j = 1-OLy
                0079         DO i=1-OLx,sNx+OLx
dff94812d5 Jean*0080            rLowS (i,j,bi,bj) = rEmpty
                0081            rSurfS(i,j,bi,bj) = rEmpty
21ab794b7f Jean*0082         ENDDO
                0083         DO j=1-OLy,sNy+OLy
                0084          DO i=2-OLx,sNx+OLx
                0085            rLowW(i,j,bi,bj)  =
                0086      &           MAX(   R_low(i-1,j,bi,bj),   R_low(i,j,bi,bj) )
                0087            rSurfW(i,j,bi,bj) =
                0088      &           MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
                0089          ENDDO
                0090         ENDDO
                0091         DO j=2-OLy,sNy+OLy
                0092          DO i=1-OLx,sNx+OLx
                0093            rLowS(i,j,bi,bj)  =
                0094      &           MAX(   R_low(i,j-1,bi,bj),   R_low(i,j,bi,bj) )
                0095            rSurfS(i,j,bi,bj) =
                0096      &           MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
                0097          ENDDO
                0098         ENDDO
                0099        ENDDO
                0100       ENDDO
                0101 
89a48da7c6 Alis*0102       DO bj=myByLo(myThid), myByHi(myThid)
                0103        DO bi=myBxLo(myThid), myBxHi(myThid)
21ab794b7f Jean*0104 
                0105 C--   Calculate lopping factor hFacC : over-estimate the part inside of the domain
                0106 C     taking into account the lower_R Boundary (Bathymetry / Top of Atmos)
86566ea478 Jean*0107         DO k=1, Nr
a5212b6221 Jean*0108          hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
b45b3dc2f0 Jean*0109          DO j=1-OLy,sNy+OLy
                0110           DO i=1-OLx,sNx+OLx
aea29c8517 Alis*0111 C      o Non-dimensional distance between grid bound. and domain lower_R bound.
a2f3aac956 Jean*0112            hFac_loc = (rF(k)-R_low(i,j,bi,bj))*recip_drF(k)
aea29c8517 Alis*0113 C      o Select between, closed, open or partial (0,1,0-1)
a2f3aac956 Jean*0114            hFac_loc = MIN( MAX( hFac_loc, zeroRL ) , oneRL )
aea29c8517 Alis*0115 C      o Impose minimum fraction and/or size (dimensional)
a2f3aac956 Jean*0116            IF ( hFac_loc.LT.hFacMnSz*halfRL .OR.
36941d6f03 Jean*0117      &          R_low(i,j,bi,bj).GE.Ro_surf(i,j,bi,bj) ) THEN
a2f3aac956 Jean*0118              hFacC(i,j,k,bi,bj) = zeroRS
aea29c8517 Alis*0119            ELSE
a2f3aac956 Jean*0120              hFacC(i,j,k,bi,bj) = MAX( hFac_loc, hFacMnSz )
89a48da7c6 Alis*0121            ENDIF
                0122           ENDDO
                0123          ENDDO
                0124         ENDDO
aea29c8517 Alis*0125 
86566ea478 Jean*0126 C-    Re-calculate lower-R Boundary position, taking into account hFacC
b45b3dc2f0 Jean*0127         DO j=1-OLy,sNy+OLy
                0128          DO i=1-OLx,sNx+OLx
a2f3aac956 Jean*0129            tmpVar(i,j) = 0. _d 0
bbebcf4c4e Mart*0130          ENDDO
                0131         ENDDO
21ab794b7f Jean*0132         DO k=1,Nr
b45b3dc2f0 Jean*0133          DO j=1-OLy,sNy+OLy
                0134           DO i=1-OLx,sNx+OLx
dff94812d5 Jean*0135            tmpVar(i,j) = tmpVar(i,j) + drF(k)*hFacC(i,j,k,bi,bj)
aea29c8517 Alis*0136           ENDDO
                0137          ENDDO
                0138         ENDDO
21ab794b7f Jean*0139         DO j=1-OLy,sNy+OLy
                0140          DO i=1-OLx,sNx+OLx
dff94812d5 Jean*0141            R_low(i,j,bi,bj) = rF(1) - tmpVar(i,j)
21ab794b7f Jean*0142          ENDDO
                0143         ENDDO
c712b7e80b Chri*0144 
86566ea478 Jean*0145 C--   Calculate lopping factor hFacC : Remove part outside of the domain
                0146 C     taking into account the Reference (=at rest) Surface Position Ro_surf
                0147         DO k=1, Nr
a5212b6221 Jean*0148          hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
b45b3dc2f0 Jean*0149          DO j=1-OLy,sNy+OLy
                0150           DO i=1-OLx,sNx+OLx
aea29c8517 Alis*0151 C      o Non-dimensional distance between grid boundary and model surface
a2f3aac956 Jean*0152            hFac_loc = (rF(k)-Ro_surf(i,j,bi,bj))*recip_drF(k)
aea29c8517 Alis*0153 C      o Reduce the previous fraction : substract the outside part.
a2f3aac956 Jean*0154            hFac_loc = hFacC(i,j,k,bi,bj) - MAX( hFac_loc, zeroRL )
aea29c8517 Alis*0155 C      o set to zero if empty Column :
a2f3aac956 Jean*0156            hFac_loc = MAX( hFac_loc, zeroRL )
aea29c8517 Alis*0157 C      o Impose minimum fraction and/or size (dimensional)
a2f3aac956 Jean*0158            IF ( hFac_loc.LT.hFacMnSz*halfRL ) THEN
                0159              hFacC(i,j,k,bi,bj) = zeroRS
aea29c8517 Alis*0160            ELSE
a2f3aac956 Jean*0161              hFacC(i,j,k,bi,bj) = MAX( hFac_loc, hFacMnSz )
aea29c8517 Alis*0162            ENDIF
                0163           ENDDO
                0164          ENDDO
                0165         ENDDO
a6cbc7a360 Mart*0166 
86566ea478 Jean*0167 C-    Re-calculate Reference surface position, taking into account hFacC
                0168 C     initialize Total column fluid thickness and surface k index
                0169 C       Note: if no fluid (continent) ==> kSurf = Nr+1
b45b3dc2f0 Jean*0170         DO j=1-OLy,sNy+OLy
                0171          DO i=1-OLx,sNx+OLx
a2f3aac956 Jean*0172           tmpVar(i,j) = 0. _d 0
86566ea478 Jean*0173           kSurfC(i,j,bi,bj) = Nr+1
21ab794b7f Jean*0174           kLowC (i,j,bi,bj) = 0
                0175          ENDDO
                0176         ENDDO
                0177         DO k=1,Nr
                0178          DO j=1-OLy,sNy+OLy
                0179           DO i=1-OLx,sNx+OLx
dff94812d5 Jean*0180            tmpVar(i,j) = tmpVar(i,j) + drF(k)*hFacC(i,j,k,bi,bj)
21ab794b7f Jean*0181            IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) kLowC(i,j,bi,bj) = k
ab42872a05 Alis*0182           ENDDO
21ab794b7f Jean*0183          ENDDO
                0184         ENDDO
                0185         DO k=Nr,1,-1
                0186          DO j=1-OLy,sNy+OLy
                0187           DO i=1-OLx,sNx+OLx
                0188            IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) kSurfC(i,j,bi,bj) = k
60c223928f Mart*0189           ENDDO
21ab794b7f Jean*0190          ENDDO
                0191         ENDDO
                0192         DO j=1-OLy,sNy+OLy
                0193          DO i=1-OLx,sNx+OLx
dff94812d5 Jean*0194           Ro_surf(i,j,bi,bj) = R_low(i,j,bi,bj) + tmpVar(i,j)
a5212b6221 Jean*0195           maskInC(i,j,bi,bj) = 0.
                0196           IF ( kSurfC(i,j,bi,bj).LE.Nr ) maskInC(i,j,bi,bj) = 1.
7b03dab2b2 jm-c 0197 c         k = MAX( 0, kLowC (i,j,bi,bj) - kSurfC(i,j,bi,bj) + 1 )
                0198 c         tmpFld(i,j,bi,bj) = k
ab42872a05 Alis*0199          ENDDO
                0200         ENDDO
21ab794b7f Jean*0201 
86566ea478 Jean*0202 C-    end bi,bj loops.
ab42872a05 Alis*0203        ENDDO
                0204       ENDDO
aea29c8517 Alis*0205 
337b46d524 Jean*0206       IF ( plotLevel.GE.debLevB ) THEN
ea3c970699 Jean*0207 c       CALL PLOT_FIELD_XYRS( tmpFld,
322e3ca99c Jean*0208 c    &           'Model Depths K Index' , -1, myThid )
137309e213 Jean*0209         CALL PLOT_FIELD_XYRS(R_low,
322e3ca99c Jean*0210      &           'Model R_low (ini_masks_etc)', -1, myThid )
137309e213 Jean*0211         CALL PLOT_FIELD_XYRS(Ro_surf,
322e3ca99c Jean*0212      &           'Model Ro_surf (ini_masks_etc)', -1, myThid )
137309e213 Jean*0213       ENDIF
ab42872a05 Alis*0214 
21ab794b7f Jean*0215 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0216 
ab42872a05 Alis*0217       DO bj = myByLo(myThid), myByHi(myThid)
                0218        DO bi = myBxLo(myThid), myBxHi(myThid)
21ab794b7f Jean*0219 
                0220 C--   Calculate quantities derived from XY depth map
b45b3dc2f0 Jean*0221         DO j=1-OLy,sNy+OLy
                0222          DO i=1-OLx,sNx+OLx
aea29c8517 Alis*0223 C         Total fluid column thickness (r_unit) :
dff94812d5 Jean*0224           tmpVar(i,j) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
                0225 c         tmpFld(i,j,bi,bj) = tmpVar(i,j)
aea29c8517 Alis*0226 C         Inverse of fluid column thickness (1/r_unit)
dff94812d5 Jean*0227           IF ( tmpVar(i,j) .LE. zeroRL ) THEN
a2f3aac956 Jean*0228            recip_Rcol(i,j,bi,bj) = zeroRS
ab42872a05 Alis*0229           ELSE
dff94812d5 Jean*0230            recip_Rcol(i,j,bi,bj) = 1. _d 0 / tmpVar(i,j)
ab42872a05 Alis*0231           ENDIF
                0232          ENDDO
                0233         ENDDO
                0234 
db15b951f3 Jean*0235 C- Method-1 (useMin4hFacEdges = T):
                0236 C    compute hFacW,hFacS as minimum of adjacent hFacC factor
                0237 C- Method-2 (useMin4hFacEdges = F):
                0238 C    compute hFacW,hFacS from rSurfW,S and rLowW,S by applying
                0239 C    same rules as for hFacC
                0240 C Note: Currently, no difference between methods except when useShelfIce=T and
                0241 C       if, in adjacent columns, ice-draft and bathy are within the same level k
                0242 
21ab794b7f Jean*0243         IF ( useMin4hFacEdges ) THEN
                0244 C--   hFacW and hFacS (at U and V points):
                0245 C-    Method-1: use simply minimum of adjacent hFacC factor
a5212b6221 Jean*0246 
21ab794b7f Jean*0247          DO k=1, Nr
                0248           DO j=1-OLy,sNy+OLy
a2f3aac956 Jean*0249            hFacW(1-OLx,j,k,bi,bj) = zeroRS
21ab794b7f Jean*0250            DO i=2-OLx,sNx+OLx
                0251             hFacW(i,j,k,bi,bj) =
                0252      &        MIN( hFacC(i,j,k,bi,bj), hFacC(i-1,j,k,bi,bj) )
                0253            ENDDO
                0254           ENDDO
                0255           DO i=1-OLx,sNx+OLx
a2f3aac956 Jean*0256             hFacS(i,1-OLy,k,bi,bj) = zeroRS
21ab794b7f Jean*0257           ENDDO
                0258           DO j=2-OLy,sNy+OLy
                0259            DO i=1-OLx,sNx+OLx
                0260             hFacS(i,j,k,bi,bj) =
                0261      &        MIN( hFacC(i,j,k,bi,bj), hFacC(i,j-1,k,bi,bj) )
                0262            ENDDO
                0263           ENDDO
                0264          ENDDO
                0265 
                0266         ELSE
                0267 C--   hFacW and hFacS (at U and V points):
                0268 C-    Method-2: compute new hFacW,S from rSurfW,S and rLowW,S
                0269 C               by applying same rules as for hFacC
                0270 
                0271          DO k=1, Nr
                0272           hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
                0273           DO j=1-OLy,sNy+OLy
                0274            DO i=1-OLx,sNx+OLx
                0275 C      o Non-dimensional distance between grid bound. and domain lower_R bound.
                0276             hFac1tmp = ( rF(k) - rLowW(i,j,bi,bj) )*recip_drF(k)
a2f3aac956 Jean*0277             hFac_loc = MIN( hFac1tmp, oneRL )
                0278 c           hFac_loc = MAX( hFac_loc, zeroRL )
21ab794b7f Jean*0279 C      o Impose minimum fraction and/or size (dimensional)
a2f3aac956 Jean*0280             IF ( hFac_loc.LT.hFacMnSz*halfRL .OR.
36941d6f03 Jean*0281      &           rLowW(i,j,bi,bj).GE.rSurfW(i,j,bi,bj) ) THEN
a2f3aac956 Jean*0282               hFac1tmp = 0. _d 0
21ab794b7f Jean*0283             ELSE
a2f3aac956 Jean*0284               hFac1tmp = MAX( hFac_loc, hFacMnSz )
21ab794b7f Jean*0285             ENDIF
                0286 C      o Reduce the previous fraction : substract the outside fraction
                0287 C        (i.e., beyond reference (=at rest) surface position rSurfW)
                0288             hFac2tmp = ( rF(k) -rSurfW(i,j,bi,bj) )*recip_drF(k)
a2f3aac956 Jean*0289             hFac_loc = hFac1tmp - MAX( hFac2tmp, zeroRL )
21ab794b7f Jean*0290 C      o Impose minimum fraction and/or size (dimensional)
a2f3aac956 Jean*0291             IF ( hFac_loc.LT.hFacMnSz*halfRL ) THEN
                0292               hFacW(i,j,k,bi,bj) = zeroRS
21ab794b7f Jean*0293             ELSE
a2f3aac956 Jean*0294               hFacW(i,j,k,bi,bj) = MAX( hFac_loc, hFacMnSz )
21ab794b7f Jean*0295             ENDIF
                0296            ENDDO
                0297           ENDDO
                0298           DO j=1-OLy,sNy+OLy
                0299            DO i=1-OLx,sNx+OLx
                0300 C      o Non-dimensional distance between grid bound. and domain lower_R bound.
                0301             hFac1tmp = ( rF(k) - rLowS(i,j,bi,bj) )*recip_drF(k)
a2f3aac956 Jean*0302             hFac_loc = MIN( hFac1tmp, oneRL )
                0303 c           hFac_loc = MAX( hFac_loc, zeroRL )
21ab794b7f Jean*0304 C      o Impose minimum fraction and/or size (dimensional)
a2f3aac956 Jean*0305             IF ( hFac_loc.LT.hFacMnSz*halfRL .OR.
36941d6f03 Jean*0306      &           rLowS(i,j,bi,bj).GE.rSurfS(i,j,bi,bj) ) THEN
a2f3aac956 Jean*0307               hFac1tmp = 0. _d 0
21ab794b7f Jean*0308             ELSE
a2f3aac956 Jean*0309               hFac1tmp = MAX( hFac_loc, hFacMnSz )
21ab794b7f Jean*0310             ENDIF
                0311 C      o Reduce the previous fraction : substract the outside fraction
                0312 C        (i.e., beyond reference (=at rest) surface position rSurfS)
                0313             hFac2tmp = ( rF(k) -rSurfS(i,j,bi,bj) )*recip_drF(k)
a2f3aac956 Jean*0314             hFac_loc = hFac1tmp - MAX( hFac2tmp, zeroRL )
21ab794b7f Jean*0315 C      o Impose minimum fraction and/or size (dimensional)
a2f3aac956 Jean*0316             IF ( hFac_loc.LT.hFacMnSz*halfRL ) THEN
                0317               hFacS(i,j,k,bi,bj) = zeroRS
21ab794b7f Jean*0318             ELSE
a2f3aac956 Jean*0319               hFacS(i,j,k,bi,bj) = MAX( hFac_loc, hFacMnSz )
21ab794b7f Jean*0320             ENDIF
                0321            ENDDO
                0322           ENDDO
                0323          ENDDO
                0324         ENDIF
                0325 
                0326 C--   Update rLow & reference rSurf at Western & Southern edges (U & V pts):
                0327 C     account for adjusted R_low & Ro_surf due to hFacMin constrain on hFacC.
                0328 C     Might need further adjustment (e.g., if useShelfIce=T) to match
                0329 C     integrated level thickness ( =Sum_k(drF*hFac) )
b45b3dc2f0 Jean*0330         DO j=1-OLy,sNy+OLy
                0331          DO i=2-OLx,sNx+OLx
a239fc1811 Jean*0332            rLowW(i,j,bi,bj)  =
                0333      &           MAX(   R_low(i-1,j,bi,bj),   R_low(i,j,bi,bj) )
f88bbe67c4 Jean*0334            rSurfW(i,j,bi,bj) =
                0335      &           MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
                0336            rSurfW(i,j,bi,bj) =
                0337      &           MAX( rSurfW(i,j,bi,bj), rLowW(i,j,bi,bj) )
a239fc1811 Jean*0338          ENDDO
                0339         ENDDO
b45b3dc2f0 Jean*0340         DO j=2-OLy,sNy+OLy
                0341          DO i=1-OLx,sNx+OLx
a239fc1811 Jean*0342            rLowS(i,j,bi,bj)  =
                0343      &           MAX(   R_low(i,j-1,bi,bj),   R_low(i,j,bi,bj) )
f88bbe67c4 Jean*0344            rSurfS(i,j,bi,bj) =
                0345      &           MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
                0346            rSurfS(i,j,bi,bj) =
                0347      &           MAX( rSurfS(i,j,bi,bj), rLowS(i,j,bi,bj) )
a239fc1811 Jean*0348          ENDDO
                0349         ENDDO
a5212b6221 Jean*0350 
db15b951f3 Jean*0351 c       IF ( useShelfIce ) THEN
dff94812d5 Jean*0352 C--   Adjust reference rSurf at Western & Southern edges (U & V pts)
21ab794b7f Jean*0353 C     to get consistent column thickness from Sum_k(hFac*drF) and rSurf-rLow
                0354 
dff94812d5 Jean*0355 C-    Total column thickness at Western edge
a5212b6221 Jean*0356          DO j=1-OLy,sNy+OLy
                0357           DO i=1-OLx,sNx+OLx
dff94812d5 Jean*0358             tmpVar(i,j) = 0. _d 0
a5212b6221 Jean*0359           ENDDO
                0360          ENDDO
21ab794b7f Jean*0361          DO k=1,Nr
a5212b6221 Jean*0362           DO j=1-OLy,sNy+OLy
                0363            DO i=1-OLx,sNx+OLx
dff94812d5 Jean*0364             tmpVar(i,j) = tmpVar(i,j) + drF(k)*hFacW(i,j,k,bi,bj)
a5212b6221 Jean*0365            ENDDO
                0366           ENDDO
21ab794b7f Jean*0367          ENDDO
dff94812d5 Jean*0368 C-    Adjust rSurf at W edge (correct for the difference)
                0369          DO j=1-OLy,sNy+OLy
                0370           DO i=1-OLx,sNx+OLx
                0371              rSurfW(i,j,bi,bj) = rLowW(i,j,bi,bj) + tmpVar(i,j)
                0372           ENDDO
                0373          ENDDO
21ab794b7f Jean*0374 
dff94812d5 Jean*0375 C-    Total column thickness at Southern edges
                0376          DO j=1-OLy,sNy+OLy
                0377           DO i=1-OLx,sNx+OLx
                0378             tmpVar(i,j) = 0. _d 0
                0379           ENDDO
                0380          ENDDO
                0381          DO k=1,Nr
                0382           DO j=1-OLy,sNy+OLy
                0383            DO i=1-OLx,sNx+OLx
                0384             tmpVar(i,j) = tmpVar(i,j) + drF(k)*hFacS(i,j,k,bi,bj)
21ab794b7f Jean*0385            ENDDO
dff94812d5 Jean*0386           ENDDO
                0387          ENDDO
                0388 C-    Adjust rSurf at S edge (correct for the difference)
                0389          DO j=1-OLy,sNy+OLy
                0390           DO i=1-OLx,sNx+OLx
                0391              rSurfS(i,j,bi,bj) = rLowS(i,j,bi,bj) + tmpVar(i,j)
                0392           ENDDO
                0393          ENDDO
21ab794b7f Jean*0394 
                0395 C-    end if useShelfIce
db15b951f3 Jean*0396 c       ENDIF
a5212b6221 Jean*0397 
86566ea478 Jean*0398 C-    end bi,bj loops.
a239fc1811 Jean*0399        ENDDO
                0400       ENDDO
21ab794b7f Jean*0401 
a5212b6221 Jean*0402       CALL EXCH_UV_XYZ_RS( hFacW,  hFacS, .FALSE., myThid )
a239fc1811 Jean*0403       CALL EXCH_UV_XY_RS( rSurfW, rSurfS, .FALSE., myThid )
ea3c970699 Jean*0404       CALL EXCH_UV_XY_RS(  rLowW,  rLowS, .FALSE., myThid )
a239fc1811 Jean*0405 
86566ea478 Jean*0406 C--   Calculate surface k index for interface W & S (U & V points)
                0407       DO bj=myByLo(myThid), myByHi(myThid)
                0408        DO bi=myBxLo(myThid), myBxHi(myThid)
b45b3dc2f0 Jean*0409         DO j=1-OLy,sNy+OLy
                0410          DO i=1-OLx,sNx+OLx
86566ea478 Jean*0411           kSurfW(i,j,bi,bj) = Nr+1
                0412           kSurfS(i,j,bi,bj) = Nr+1
                0413           DO k=Nr,1,-1
a5212b6221 Jean*0414            IF (hFacW(i,j,k,bi,bj).NE.zeroRS) kSurfW(i,j,bi,bj) = k
                0415            IF (hFacS(i,j,k,bi,bj).NE.zeroRS) kSurfS(i,j,bi,bj) = k
86566ea478 Jean*0416           ENDDO
a2f3aac956 Jean*0417           maskInW(i,j,bi,bj)= zeroRS
                0418           IF ( kSurfW(i,j,bi,bj).LE.Nr ) maskInW(i,j,bi,bj)= oneRS
                0419           maskInS(i,j,bi,bj)= zeroRS
                0420           IF ( kSurfS(i,j,bi,bj).LE.Nr ) maskInS(i,j,bi,bj)= oneRS
86566ea478 Jean*0421          ENDDO
                0422         ENDDO
                0423        ENDDO
                0424       ENDDO
                0425 
dff94812d5 Jean*0426 C--   Additional closing of Western and Southern grid-cell edges: for example,
                0427 C     a) might add some "thin walls" in specific location
                0428 C     b) close non-periodic N & S boundaries of lat-lon grid at the N/S poles.
                0429 C     new: location now reccorded as kSurfW/S = Nr+2
                0430       CALL ADD_WALLS2MASKS( rEmpty, myThid )
                0431 
f15994caab Jean*0432       ELSE
                0433 #ifndef DISABLE_SIGMA_CODE
                0434 C---  Sigma and Hybrid-Sigma set-up:
                0435         CALL INI_SIGMA_HFAC( myThid )
                0436 #endif /* DISABLE_SIGMA_CODE */
                0437       ENDIF
                0438 
86566ea478 Jean*0439 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0440 
                0441 C--   Write to disk: Total Column Thickness & hFac(C,W,S):
f4a7634227 Alis*0442 C     This I/O is now done in write_grid.F
ea3c970699 Jean*0443 c     CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpFld,0,myThid)
f4a7634227 Alis*0444 c     CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
                0445 c     CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
                0446 c     CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
27d8b2bb9c Jean*0447 
337b46d524 Jean*0448       IF ( plotLevel.GE.debLevB ) THEN
322e3ca99c Jean*0449         CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 0, myThid )
                0450         CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 0, myThid )
                0451         CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 0, myThid )
137309e213 Jean*0452       ENDIF
a74b31ee34 Alis*0453 
86566ea478 Jean*0454 C--   Masks and reciprocals of hFac[CWS]
a74b31ee34 Alis*0455       DO bj = myByLo(myThid), myByHi(myThid)
                0456        DO bi = myBxLo(myThid), myBxHi(myThid)
86566ea478 Jean*0457         DO k=1,Nr
b45b3dc2f0 Jean*0458          DO j=1-OLy,sNy+OLy
                0459           DO i=1-OLx,sNx+OLx
a5212b6221 Jean*0460            IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) THEN
86566ea478 Jean*0461             recip_hFacC(i,j,k,bi,bj) = 1. _d 0 / hFacC(i,j,k,bi,bj)
a2f3aac956 Jean*0462             maskC(i,j,k,bi,bj) = oneRS
a74b31ee34 Alis*0463            ELSE
a2f3aac956 Jean*0464             recip_hFacC(i,j,k,bi,bj) = zeroRS
                0465             maskC(i,j,k,bi,bj) = zeroRS
a74b31ee34 Alis*0466            ENDIF
a5212b6221 Jean*0467            IF ( hFacW(i,j,k,bi,bj).NE.zeroRS ) THEN
86566ea478 Jean*0468             recip_hFacW(i,j,k,bi,bj) = 1. _d 0 / hFacW(i,j,k,bi,bj)
a2f3aac956 Jean*0469             maskW(i,j,k,bi,bj) = oneRS
a74b31ee34 Alis*0470            ELSE
a2f3aac956 Jean*0471             recip_hFacW(i,j,k,bi,bj) = zeroRS
                0472             maskW(i,j,k,bi,bj) = zeroRS
a74b31ee34 Alis*0473            ENDIF
a5212b6221 Jean*0474            IF ( hFacS(i,j,k,bi,bj).NE.zeroRS ) THEN
86566ea478 Jean*0475             recip_hFacS(i,j,k,bi,bj) = 1. _d 0 / hFacS(i,j,k,bi,bj)
a2f3aac956 Jean*0476             maskS(i,j,k,bi,bj) = oneRS
a74b31ee34 Alis*0477            ELSE
a2f3aac956 Jean*0478             recip_hFacS(i,j,k,bi,bj) = zeroRS
                0479             maskS(i,j,k,bi,bj) = zeroRS
a74b31ee34 Alis*0480            ENDIF
                0481           ENDDO
                0482          ENDDO
                0483         ENDDO
adc74f71dd Jean*0484 #ifdef NONLIN_FRSURF
                0485 C--   Save initial geometrical hFac factor into h0Fac (fixed in time):
                0486 C     Note: In case 1 pkg modifies hFac (from packages_init_fixed, called
                0487 C     later in sequence of calls) this pkg would need also to update h0Fac.
                0488         DO k=1,Nr
                0489          DO j=1-OLy,sNy+OLy
                0490           DO i=1-OLx,sNx+OLx
                0491            h0FacC(i,j,k,bi,bj) = _hFacC(i,j,k,bi,bj)
                0492            h0FacW(i,j,k,bi,bj) = _hFacW(i,j,k,bi,bj)
                0493            h0FacS(i,j,k,bi,bj) = _hFacS(i,j,k,bi,bj)
                0494           ENDDO
                0495          ENDDO
                0496         ENDDO
                0497 #endif /* NONLIN_FRSURF */
86566ea478 Jean*0498 C-    end bi,bj loops.
a74b31ee34 Alis*0499        ENDDO
                0500       ENDDO
                0501 
ab5a98a4ed Jean*0502 c #ifdef ALLOW_NONHYDROSTATIC
                0503 C--   Calculate "recip_hFacU" = reciprocal hfac distance/volume for W cells
d34bd75e92 Jean*0504 C NOTE:  not used ; computed locally in CALC_GW
ab5a98a4ed Jean*0505 c #endif
                0506 
a74b31ee34 Alis*0507       RETURN
                0508       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!