Back to home page

MITgcm
 
 

    

File indexing completed on 2020-02-14 06:10:28 UTC

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