Back to home page

MITgcm
 
 

    

File indexing completed on 2020-02-28 06:11:24 UTC

view on githubraw file Latest commit 58240850 on 2019-12-26 20:51:18 UTC
7155c3b3df Jean*0001 #include "SHELFICE_OPTIONS.h"
                0002 CBOP
                0003 C     !ROUTINE: SHELFICE_REMESH_UV_MASK
                0004 C     !INTERFACE:
                0005       SUBROUTINE SHELFICE_REMESH_UV_MASK(
                0006      O                    k1SurfW, k1SurfS, mrgFacW, mrgFacS,
                0007      I                    myTime, myIter, myThid )
                0008 
                0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
                0011 C     | SUBROUTINE SHELFICE_REMESH_UV_MASK
                0012 C     | o Update masks and geom factors at U and V points
                0013 C     | o save surface level index and thickness factors for
                0014 C     |   update of horizontal velocity state variables
                0015 C     *==========================================================*
                0016 C     \ev
                0017 
                0018 C     !USES:
                0019       IMPLICIT NONE
                0020 C     === Global variables ===
                0021 #include "SIZE.h"
                0022 #include "EEPARAMS.h"
                0023 #include "PARAMS.h"
                0024 #include "GRID.h"
                0025 #ifdef NONLIN_FRSURF
                0026 # include "SURFACE.h"
                0027 #endif /* NONLIN_FRSURF */
                0028 
                0029 C     !INPUT/OUTPUT PARAMETERS:
                0030 C     k1SurfW   :: surface level index (at cell W.Edge) before remeshing
                0031 C     k1SurfS   :: surface level index (at cell S.Edge) before remeshing
                0032 C     mrgFacW   :: store hFacW to compute merging weight for U velocity
                0033 C     mrgFacS   :: store hFacS to compute merging weight for V velocity
                0034 C     myTime    :: Current time in simulation
                0035 C     myIter    :: Current iteration number
                0036 C     myThid    :: my Thread Id number
                0037       INTEGER k1SurfW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0038       INTEGER k1SurfS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0039       _RL mrgFacW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,nSx,nSy)
                0040       _RL mrgFacS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,nSx,nSy)
                0041       _RL myTime
                0042       INTEGER myIter
                0043       INTEGER myThid
                0044 
                0045 #ifdef ALLOW_SHELFICE_REMESHING
                0046 #ifdef NONLIN_FRSURF
                0047 C     !LOCAL VARIABLES:
                0048 C     bi, bj    :: tile indices
                0049 C     i, j, k   :: Loop counters
                0050 C     tmpVar    :: Temporary array used to integrate column thickness
                0051       INTEGER bi, bj
                0052       INTEGER i, j, k
                0053       INTEGER ks, k1, k2
                0054       _RL tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0055       _RL hFacMnSz, hFacCtmp
                0056       _RL hFac1tmp, hFac2tmp
                0057       _RL Rmin_tmp, hfacInfMom
                0058 CEOP
                0059 
                0060 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0061 
                0062       DO bj = myByLo(myThid), myByHi(myThid)
                0063        DO bi = myBxLo(myThid), myBxHi(myThid)
                0064 
                0065         DO j=1-OLy,sNy+OLy
                0066          DO i=1-OLx,sNx+OLx
                0067            k1SurfW(i,j,bi,bj) = kSurfW(i,j,bi,bj)
                0068            k1SurfS(i,j,bi,bj) = kSurfS(i,j,bi,bj)
                0069            mrgFacW(i,j,1,bi,bj) = 0.
                0070            mrgFacW(i,j,2,bi,bj) = 0.
                0071            mrgFacS(i,j,1,bi,bj) = 0.
                0072            mrgFacS(i,j,2,bi,bj) = 0.
                0073          ENDDO
                0074         ENDDO
                0075 
                0076         DO j=1-OLy,sNy+OLy
                0077          DO i=2-OLx,sNx+OLx
                0078            rSurfW(i,j,bi,bj) =
                0079      &           MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
                0080            rSurfW(i,j,bi,bj) =
                0081      &           MAX( rSurfW(i,j,bi,bj), rLowW(i,j,bi,bj) )
                0082          ENDDO
                0083         ENDDO
                0084         DO j=2-OLy,sNy+OLy
                0085          DO i=1-OLx,sNx+OLx
                0086            rSurfS(i,j,bi,bj) =
                0087      &           MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
                0088            rSurfS(i,j,bi,bj) =
                0089      &           MAX( rSurfS(i,j,bi,bj), rLowS(i,j,bi,bj) )
                0090          ENDDO
                0091         ENDDO
                0092 
                0093        ENDDO
                0094       ENDDO
                0095 
                0096       DO bj = myByLo(myThid), myByHi(myThid)
                0097        DO bi = myBxLo(myThid), myBxHi(myThid)
                0098 
                0099 C- Method-1 (useMin4hFacEdges = T):
                0100 C    compute hFacW,hFacS as minimum of adjacent hFacC factor
                0101 C- Method-2 (useMin4hFacEdges = F):
                0102 C    compute hFacW,hFacS from rSurfW,S and rLowW,S by applying
                0103 C    same rules as for hFacC
                0104 C Note: Currently, no difference between methods except when useShelfIce=T and
                0105 C       if, in adjacent columns, ice-draft and bathy are within the same level k
                0106 
                0107         IF ( useMin4hFacEdges ) THEN
                0108 C--   hFacW and hFacS (at U and V points):
                0109 C-    Method-1: use simply minimum of adjacent hFacC factor
                0110 
                0111          DO k=1, Nr
                0112           DO j=1-OLy,sNy+OLy
                0113            h0FacW(1-OLx,j,k,bi,bj) = 0.
                0114            DO i=2-OLx,sNx+OLx
                0115             h0FacW(i,j,k,bi,bj) =
                0116      &        MIN( h0FacC(i,j,k,bi,bj), h0FacC(i-1,j,k,bi,bj) )
                0117            ENDDO
                0118           ENDDO
                0119           DO i=1-OLx,sNx+OLx
                0120             h0FacS(i,1-OLy,k,bi,bj) = 0.
                0121           ENDDO
                0122           DO j=2-OLy,sNy+OLy
                0123            DO i=1-OLx,sNx+OLx
                0124             h0FacS(i,j,k,bi,bj) =
                0125      &        MIN( h0FacC(i,j,k,bi,bj), h0FacC(i,j-1,k,bi,bj) )
                0126            ENDDO
                0127           ENDDO
                0128          ENDDO
                0129 
                0130         ELSE
                0131 C--   hFacW and hFacS (at U and V points):
                0132 C-    Method-2: compute new hFacW,S from rSurfW,S and rLowW,S
                0133 C               by applying same rules as for hFacC
                0134 
                0135          DO k=1, Nr
                0136           hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
                0137           DO j=1-OLy,sNy+OLy
                0138            DO i=1-OLx,sNx+OLx
                0139 C      o Non-dimensional distance between grid bound. and domain lower_R bound.
                0140             hFac1tmp = ( rF(k) - rLowW(i,j,bi,bj) )*recip_drF(k)
                0141             hFacCtmp = MIN( hFac1tmp, oneRL )
                0142 c           hFacCtmp = MAX( hFacCtmp, zeroRL )
                0143 C      o Impose minimum fraction and/or size (dimensional)
                0144             IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
                0145               hFac1tmp = 0.
                0146             ELSE
                0147               hFac1tmp = MAX( hFacCtmp, hFacMnSz )
                0148             ENDIF
                0149 C      o Reduce the previous fraction : substract the outside fraction
                0150 C        (i.e., beyond reference (=at rest) surface position rSurfW)
                0151             hFac2tmp = ( rF(k) -rSurfW(i,j,bi,bj) )*recip_drF(k)
                0152             hFacCtmp = hFac1tmp - MAX( hFac2tmp, zeroRL )
                0153 C      o Impose minimum fraction and/or size (dimensional)
                0154             IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
                0155               h0FacW(i,j,k,bi,bj) = 0.
                0156             ELSE
                0157               h0FacW(i,j,k,bi,bj) = MAX( hFacCtmp, hFacMnSz )
                0158             ENDIF
                0159            ENDDO
                0160           ENDDO
                0161           DO j=1-OLy,sNy+OLy
                0162            DO i=1-OLx,sNx+OLx
                0163 C      o Non-dimensional distance between grid bound. and domain lower_R bound.
                0164             hFac1tmp = ( rF(k) - rLowS(i,j,bi,bj) )*recip_drF(k)
                0165             hFacCtmp = MIN( hFac1tmp, oneRL )
                0166 c           hFacCtmp = MAX( hFacCtmp, zeroRL )
                0167 C      o Impose minimum fraction and/or size (dimensional)
                0168             IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
                0169               hFac1tmp = 0.
                0170             ELSE
                0171               hFac1tmp = MAX( hFacCtmp, hFacMnSz )
                0172             ENDIF
                0173 C      o Reduce the previous fraction : substract the outside fraction
                0174 C        (i.e., beyond reference (=at rest) surface position rSurfS)
                0175             hFac2tmp = ( rF(k) -rSurfS(i,j,bi,bj) )*recip_drF(k)
                0176             hFacCtmp = hFac1tmp - MAX( hFac2tmp, zeroRL )
                0177 C      o Impose minimum fraction and/or size (dimensional)
                0178             IF ( hFacCtmp.LT.hFacMnSz*halfRL ) THEN
                0179               h0FacS(i,j,k,bi,bj) = 0.
                0180             ELSE
                0181               h0FacS(i,j,k,bi,bj) = MAX( hFacCtmp, hFacMnSz )
                0182             ENDIF
                0183            ENDDO
                0184           ENDDO
                0185          ENDDO
                0186         ENDIF
                0187 
                0188 C--   Update rLow & reference rSurf at Western & Southern edges (U & V pts):
                0189 C     account for adjusted R_low & Ro_surf due to hFacMin constrain on hFacC.
                0190 C     Might need further adjustment (e.g., if useShelfIce=T) to match
                0191 C     integrated level thickness ( =Sum_k(drF*hFac) )
                0192 
                0193 C--   Adjust rLow & reference rSurf at Western & Southern edges (U & V pts)
                0194 C     to get consistent column thickness from Sum_k(hFac*drF) and rSurf-rLow
                0195 
                0196 C-    Total column thickness at Western edges
                0197         DO j=1-OLy,sNy+OLy
                0198          DO i=1-OLx,sNx+OLx
                0199            tmpVar(i,j) = 0. _d 0
                0200          ENDDO
                0201         ENDDO
                0202         DO k=1,Nr
                0203          DO j=1-OLy,sNy+OLy
                0204           DO i=1-OLx,sNx+OLx
                0205             tmpVar(i,j) = tmpVar(i,j) + drF(k)*h0FacW(i,j,k,bi,bj)
                0206           ENDDO
                0207          ENDDO
                0208         ENDDO
                0209 C-    Adjust only rSurf at W (correct for the difference)
                0210         DO j=1-OLy,sNy+OLy
                0211          DO i=1-OLx,sNx+OLx
                0212            rSurfW(i,j,bi,bj) = rLowW(i,j,bi,bj) + tmpVar(i,j)
                0213          ENDDO
                0214         ENDDO
                0215 
                0216 C-    Total column thickness at Southern edges
                0217         DO j=1-OLy,sNy+OLy
                0218          DO i=1-OLx,sNx+OLx
                0219            tmpVar(i,j) = 0. _d 0
                0220          ENDDO
                0221         ENDDO
                0222         DO k=1,Nr
                0223          DO j=1-OLy,sNy+OLy
                0224           DO i=1-OLx,sNx+OLx
                0225             tmpVar(i,j) = tmpVar(i,j) + drF(k)*h0FacS(i,j,k,bi,bj)
                0226           ENDDO
                0227          ENDDO
                0228         ENDDO
                0229 C-    Adjust only rSurf at edges (correct for the difference)
                0230         DO j=1-OLy,sNy+OLy
                0231          DO i=1-OLx,sNx+OLx
                0232            rSurfS(i,j,bi,bj) = rLowS(i,j,bi,bj) + tmpVar(i,j)
                0233          ENDDO
                0234         ENDDO
                0235 
                0236 C-    end bi,bj loops.
                0237        ENDDO
                0238       ENDDO
                0239 
                0240       CALL EXCH_UV_XYZ_RS( h0FacW,  h0FacS, .FALSE., myThid )
                0241       CALL EXCH_UV_XY_RS( rSurfW, rSurfS, .FALSE., myThid )
                0242 
                0243       DO bj=myByLo(myThid), myByHi(myThid)
                0244        DO bi=myBxLo(myThid), myBxHi(myThid)
                0245 
                0246 C--   Account for additional closing of Western and Southern grid-cell edges
                0247 C     ( which jave been reccorded as kSurfW/S = Nr+2 )
                0248 C     and calculate surface k index for interface W & S (U & V points)
                0249         DO j=1-OLy,sNy+OLy
                0250          DO i=1-OLx,sNx+OLx
                0251           kSurfW(i,j,bi,bj) = MAX( k1SurfW(i,j,bi,bj), Nr+1 )
                0252           kSurfS(i,j,bi,bj) = MAX( k1SurfS(i,j,bi,bj), Nr+1 )
                0253          ENDDO
                0254         ENDDO
                0255         DO k=1,Nr
                0256          DO j=1-OLy,sNy+OLy
                0257           DO i=1-OLx,sNx+OLx
                0258            IF ( kSurfW(i,j,bi,bj).EQ.Nr+2 ) THEN
                0259              h0FacW(i,j,k,bi,bj) = zeroRS
                0260            ELSEIF ( kSurfW(i,j,bi,bj).EQ.Nr+1 ) THEN
                0261              IF (h0FacW(i,j,k,bi,bj).NE.zeroRS) kSurfW(i,j,bi,bj) = k
                0262            ENDIF
                0263            IF ( kSurfS(i,j,bi,bj).EQ.Nr+2 ) THEN
                0264              h0FacS(i,j,k,bi,bj) = zeroRS
                0265            ELSEIF ( kSurfS(i,j,bi,bj).EQ.Nr+1 ) THEN
                0266              IF (h0FacS(i,j,k,bi,bj).NE.zeroRS) kSurfS(i,j,bi,bj) = k
                0267            ENDIF
                0268           ENDDO
                0269          ENDDO
                0270         ENDDO
                0271 
                0272 C-    end bi,bj loops.
                0273        ENDDO
                0274       ENDDO
                0275 
                0276 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0277 
                0278       DO bj = myByLo(myThid), myByHi(myThid)
                0279        DO bi = myBxLo(myThid), myBxHi(myThid)
                0280 
                0281 C--  Compute the mimimum value of r_surf (used for computing hFac_surfC)
                0282         hFacInfMOM = hFacInf
                0283         DO j=1,sNy
                0284          DO i=1,sNx
                0285            ks = kSurfC(i,j,bi,bj)
                0286            IF (ks.LE.Nr) THEN
                0287              Rmin_tmp = rF(ks+1)
                0288              IF ( ks.EQ.kSurfW(i,j,bi,bj))
                0289      &          Rmin_tmp = MAX(Rmin_tmp, R_low(i-1,j,bi,bj))
                0290              IF ( ks.EQ.kSurfW(i+1,j,bi,bj))
                0291      &          Rmin_tmp = MAX(Rmin_tmp, R_low(i+1,j,bi,bj))
                0292              IF ( ks.EQ.kSurfS(i,j,bi,bj))
                0293      &          Rmin_tmp = MAX(Rmin_tmp, R_low(i,j-1,bi,bj))
                0294              IF ( ks.EQ.kSurfS(i,j+1,bi,bj))
                0295      &          Rmin_tmp = MAX(Rmin_tmp, R_low(i,j+1,bi,bj))
                0296 
                0297              Rmin_surf(i,j,bi,bj) =
                0298      &        MAX( MAX(rF(ks+1),R_low(i,j,bi,bj)) + hFacInf*drF(ks),
                0299      &                                Rmin_tmp + hFacInfMOM*drF(ks)
                0300      &           )
                0301            ENDIF
                0302          ENDDO
                0303         ENDDO
                0304 
                0305 C--  Prepare merging weights for both components of horizontal velocity;
                0306 C    they will get updated once new thickness factor are available
                0307         DO j=1-OLy,sNy+OLy
                0308          DO i=1-OLx,sNx+OLx
                0309           IF ( kSurfW(i,j,bi,bj).NE.k1SurfW(i,j,bi,bj) ) THEN
                0310             k1 = k1SurfW(i,j,bi,bj)
                0311             k2 =  kSurfW(i,j,bi,bj)
                0312             IF ( k2.GT.k1 .AND. k2.LE.Nr ) THEN
                0313 C-      merging former (=k1) into new (=k2) surface grid cell
                0314               mrgFacW(i,j,1,bi,bj) = hFacW(i,j,k1,bi,bj)
                0315               mrgFacW(i,j,2,bi,bj) = hFacW(i,j,k2,bi,bj)
                0316             ENDIF
                0317           ENDIF
                0318          ENDDO
                0319         ENDDO
                0320         DO j=1-OLy,sNy+OLy
                0321          DO i=1-OLx,sNx+OLx
                0322           IF ( kSurfS(i,j,bi,bj).NE.k1SurfS(i,j,bi,bj) ) THEN
                0323             k1 = k1SurfS(i,j,bi,bj)
                0324             k2 =  kSurfS(i,j,bi,bj)
                0325             IF ( k2.GT.k1 .AND. k2.LE.Nr ) THEN
                0326 C-      merging former (=k1) into new (=k2) surface grid cell
                0327               mrgFacS(i,j,1,bi,bj) = hFacS(i,j,k1,bi,bj)
                0328               mrgFacS(i,j,2,bi,bj) = hFacS(i,j,k2,bi,bj)
                0329             ENDIF
                0330           ENDIF
                0331          ENDDO
                0332         ENDDO
                0333 
                0334 C--   Masks and reciprocals of hFac[CWS]
                0335         DO k=1,Nr
                0336          DO j=1-OLy,sNy+OLy
                0337           DO i=1-OLx,sNx+OLx
                0338            IF ( h0FacW(i,j,k,bi,bj).NE.zeroRS ) THEN
                0339             recip_hFacW(i,j,k,bi,bj) = 1. _d 0 / h0FacW(i,j,k,bi,bj)
                0340             maskW(i,j,k,bi,bj) = 1.
                0341            ELSE
                0342             recip_hFacW(i,j,k,bi,bj) = 0.
                0343             maskW(i,j,k,bi,bj) = 0.
                0344            ENDIF
                0345            IF ( h0FacS(i,j,k,bi,bj).NE.zeroRS ) THEN
                0346             recip_hFacS(i,j,k,bi,bj) = 1. _d 0 / h0FacS(i,j,k,bi,bj)
                0347             maskS(i,j,k,bi,bj) = 1.
                0348            ELSE
                0349             recip_hFacS(i,j,k,bi,bj) = 0.
                0350             maskS(i,j,k,bi,bj) = 0.
                0351            ENDIF
                0352           ENDDO
                0353          ENDDO
                0354         ENDDO
                0355 
                0356 C--   Copy initial geometrical h0FacW & h0FacW into hFac factor.
                0357         DO k=1,Nr
                0358          DO j=1-OLy,sNy+OLy
                0359           DO i=1-OLx,sNx+OLx
                0360            hFacW(i,j,k,bi,bj) = h0FacW(i,j,k,bi,bj)
                0361            hFacS(i,j,k,bi,bj) = h0FacS(i,j,k,bi,bj)
                0362           ENDDO
                0363          ENDDO
                0364         ENDDO
                0365 
                0366 C-    end bi,bj loops.
                0367        ENDDO
                0368       ENDDO
                0369 
                0370       CALL EXCH_XY_RL( Rmin_surf, myThid )
                0371 
                0372 #endif /* NONLIN_FRSURF */
                0373 #endif /* ALLOW_SHELFICE_REMESHING */
                0374 
                0375       RETURN
                0376       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!