Back to home page

MITgcm

 
 

    


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

view on githubraw file Latest commit e0cf625c on 2020-02-18 17:29:07 UTC
e0cf625cb7 Jean*0001 #include "SHELFICE_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C     !ROUTINE: SHELFICE_REMESH_STATE
                0005 C     !INTERFACE:
                0006       SUBROUTINE SHELFICE_REMESH_STATE(
                0007      I                    k1SurfC, k1SurfW, k1SurfS, mrgFacC,
                0008      U                    mrgFacW, mrgFacS,
                0009      I                    myTime, myIter, myThid )
                0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
                0012 C     | SUBROUTINE SHELFICE_REMESH_STATE
                0013 C     | o Update state variables near the surface
                0014 C     |   where remeshing occurs
                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 #include "DYNVARS.h"
                0026 
                0027 C     !INPUT/OUTPUT PARAMETERS:
                0028 C     k1SurfC   :: surface level index (at cell center) before remeshing
                0029 C     k1SurfW   :: surface level index (at cell W.Edge) before remeshing
                0030 C     k1SurfS   :: surface level index (at cell S.Edge) before remeshing
                0031 C     mrgFacC   :: merging weight for tracers (cell center)
                0032 C     mrgFacW   :: merging weight for U component velocity
                0033 C     mrgFacS   :: merging weight for V component velocity
                0034 C     myTime    :: Current time in simulation
                0035 C     myIter    :: Current iteration number
                0036 C     myThid    :: my Thread Id number
                0037       INTEGER k1SurfC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0038       INTEGER k1SurfW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0039       INTEGER k1SurfS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0040       _RL mrgFacC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,nSx,nSy)
                0041       _RL mrgFacW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,nSx,nSy)
                0042       _RL mrgFacS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,nSx,nSy)
                0043       _RL myTime
                0044       INTEGER myIter
                0045       INTEGER myThid
                0046 CEOP
                0047 
                0048 #ifdef ALLOW_SHELFICE_REMESHING
                0049 C     !LOCAL VARIABLES:
                0050 C     bi,bj   :: tile indices
                0051 C     i,j,k   :: Loop counters
                0052 C     k1, k2  :: previous and new surface level
                0053 C     ioUnit  :: temp for writing msg unit
                0054 C     msgBuf  :: Informational/error message buffer
                0055       INTEGER bi, bj
                0056       INTEGER i, j
                0057       INTEGER k1, k2
                0058       _RL r_newDz
                0059 
                0060 C--  Update model state: 1) tracers
                0061       DO bj = myByLo(myThid), myByHi(myThid)
                0062        DO bi = myBxLo(myThid), myBxHi(myThid)
                0063         DO j = 1-OLy,sNy+OLy
                0064          DO i = 1-OLx,sNx+OLx
                0065           k1 = k1SurfC(i,j,bi,bj)
                0066           k2 = kSurfC(i,j,bi,bj)
                0067           IF ( k2.LT.k1 .AND. k1.LE.Nr ) THEN
                0068 C-    Splitting top cell (=k1):
                0069 C     just copy tracer value from former top (=k1) into new top cell (=k2)
                0070             theta(i,j,k2,bi,bj) = theta(i,j,k1,bi,bj)
                0071             salt(i,j,k2,bi,bj)  = salt(i,j,k1,bi,bj)
                0072 #ifdef ALLOW_ADAMSBASHFORTH_3
                0073             gtNm(i,j,k2,bi,bj,1) = gtNm(i,j,k1,bi,bj,1)
                0074             gtNm(i,j,k2,bi,bj,2) = gtNm(i,j,k1,bi,bj,2)
                0075             gsNm(i,j,k2,bi,bj,1) = gsNm(i,j,k1,bi,bj,1)
                0076             gsNm(i,j,k2,bi,bj,2) = gsNm(i,j,k1,bi,bj,2)
                0077 #else
                0078             gtNm1(i,j,k2,bi,bj) = gtNm1(i,j,k1,bi,bj)
                0079             gsNm1(i,j,k2,bi,bj) = gsNm1(i,j,k1,bi,bj)
                0080 #endif
                0081           ENDIF
                0082           IF ( k2.GT.k1 .AND. k2.LE.Nr ) THEN
                0083 C-    Merging top cell (=k1) into the one just below (=k2):
                0084 C     average tracer into new top cell (=k2) using weight factor "mrgFacC"
                0085             theta(i,j,k2,bi,bj) =
                0086      &            theta(i,j,k1,bi,bj)*mrgFacC(i,j,1,bi,bj)
                0087      &          + theta(i,j,k2,bi,bj)*mrgFacC(i,j,2,bi,bj)
                0088             theta(i,j,k1,bi,bj) = 0.0
                0089             salt (i,j,k2,bi,bj) =
                0090      &            salt(i,j,k1,bi,bj)*mrgFacC(i,j,1,bi,bj)
                0091      &          + salt(i,j,k2,bi,bj)*mrgFacC(i,j,2,bi,bj)
                0092             salt(i,j,k1,bi,bj) = 0.0
                0093 #ifdef ALLOW_ADAMSBASHFORTH_3
                0094             gtNm(i,j,k2,bi,bj,1) =
                0095      &            gtNm(i,j,k1,bi,bj,1)*mrgFacC(i,j,1,bi,bj)
                0096      &          + gtNm(i,j,k2,bi,bj,1)*mrgFacC(i,j,2,bi,bj)
                0097             gtNm(i,j,k2,bi,bj,2) =
                0098      &            gtNm(i,j,k1,bi,bj,2)*mrgFacC(i,j,1,bi,bj)
                0099      &          + gtNm(i,j,k2,bi,bj,2)*mrgFacC(i,j,2,bi,bj)
                0100             gtNm(i,j,k1,bi,bj,1)= 0.0
                0101             gtNm(i,j,k1,bi,bj,2)= 0.0
                0102 C-
                0103             gsNm(i,j,k2,bi,bj,1) =
                0104      &            gsNm(i,j,k1,bi,bj,1)*mrgFacC(i,j,1,bi,bj)
                0105      &          + gsNm(i,j,k2,bi,bj,1)*mrgFacC(i,j,2,bi,bj)
                0106             gsNm(i,j,k2,bi,bj,2) =
                0107      &            gsNm(i,j,k1,bi,bj,2)*mrgFacC(i,j,1,bi,bj)
                0108      &          + gsNm(i,j,k2,bi,bj,2)*mrgFacC(i,j,2,bi,bj)
                0109             gsNm(i,j,k1,bi,bj,1)= 0.0
                0110             gsNm(i,j,k1,bi,bj,2)= 0.0
                0111 #else
                0112             gtNm1(i,j,k2,bi,bj) =
                0113      &            gtNm1(i,j,k1,bi,bj)*mrgFacC(i,j,1,bi,bj)
                0114      &          + gtNm1(i,j,k2,bi,bj)*mrgFacC(i,j,2,bi,bj)
                0115             gtNm1(i,j,k1,bi,bj) = 0.0
                0116             gsNm1(i,j,k2,bi,bj) =
                0117      &            gsNm1(i,j,k1,bi,bj)*mrgFacC(i,j,1,bi,bj)
                0118      &          + gsNm1(i,j,k2,bi,bj)*mrgFacC(i,j,2,bi,bj)
                0119             gsNm1(i,j,k1,bi,bj) = 0.0
                0120 #endif
                0121           ENDIF
                0122          ENDDO
                0123         ENDDO
                0124        ENDDO
                0125       ENDDO
                0126 
                0127       DO bj = myByLo(myThid), myByHi(myThid)
                0128        DO bi = myBxLo(myThid), myBxHi(myThid)
                0129 
                0130 C--  Update merging weights for both components of horizontal velocity
                0131         DO j = 1-OLy,sNy+OLy
                0132          DO i = 1-OLx,sNx+OLx
                0133           IF ( kSurfW(i,j,bi,bj).NE.k1SurfW(i,j,bi,bj) ) THEN
                0134             k1 = k1SurfW(i,j,bi,bj)
                0135             k2 =  kSurfW(i,j,bi,bj)
                0136             IF ( k2.GT.k1 .AND. k2.LE.Nr ) THEN
                0137 C-      merging former (=k1) into new (=k2) surface grid cell
                0138 c             r_newDz = recip_drF(k2)/hFacW(i,j,k2,bi,bj)
                0139               r_newDz = recip_drF(k2)*recip_hFacW(i,j,k2,bi,bj)
                0140               mrgFacW(i,j,1,bi,bj) = mrgFacW(i,j,1,bi,bj)
                0141      &                              *drF(k1)*r_newDz
                0142               mrgFacW(i,j,2,bi,bj) = mrgFacW(i,j,2,bi,bj)
                0143      &                              *drF(k2)*r_newDz
                0144             ENDIF
                0145           ENDIF
                0146          ENDDO
                0147         ENDDO
                0148         DO j = 1-OLy,sNy+OLy
                0149          DO i = 1-OLx,sNx+OLx
                0150           IF ( kSurfS(i,j,bi,bj).NE.k1SurfS(i,j,bi,bj) ) THEN
                0151             k1 = k1SurfS(i,j,bi,bj)
                0152             k2 =  kSurfS(i,j,bi,bj)
                0153             IF ( k2.GT.k1 .AND. k2.LE.Nr ) THEN
                0154 C-      merging former (=k1) into new (=k2) surface grid cell
                0155 c             r_newDz = recip_drF(k2)/hFacS(i,j,k2,bi,bj)
                0156               r_newDz = recip_drF(k2)*recip_hFacS(i,j,k2,bi,bj)
                0157               mrgFacS(i,j,1,bi,bj) = mrgFacS(i,j,1,bi,bj)
                0158      &                              *drF(k1)*r_newDz
                0159               mrgFacS(i,j,2,bi,bj) = mrgFacS(i,j,2,bi,bj)
                0160      &                              *drF(k2)*r_newDz
                0161             ENDIF
                0162           ENDIF
                0163          ENDDO
                0164         ENDDO
                0165 
                0166 C--  Update model state: 2) horizontal velocity components U & V
                0167         DO j = 1-OLy,sNy+OLy
                0168          DO i = 1-OLx,sNx+OLx
                0169           IF ( k1SurfW(i,j,bi,bj).NE.kSurfW(i,j,bi,bj) ) THEN
                0170            k1 = k1SurfW(i,j,bi,bj)
                0171            k2 =  kSurfW(i,j,bi,bj)
                0172            IF ( k2.LT.k1 .AND. k1.LE.Nr ) THEN
                0173 C-      spliting case: just copy former surface value (=k1) into new one (=k2)
                0174             uVel(i,j,k2,bi,bj) = uVel(i,j,k1,bi,bj)
                0175 #ifdef ALLOW_ADAMSBASHFORTH_3
                0176             guNm(i,j,k2,bi,bj,1) = guNm(i,j,k1,bi,bj,1)
                0177             guNm(i,j,k2,bi,bj,2) = guNm(i,j,k1,bi,bj,2)
                0178 #else
                0179             guNm1(i,j,k2,bi,bj) = guNm1(i,j,k1,bi,bj)
                0180 #endif
                0181            ENDIF
                0182            IF ( k2.GT.k1 .AND. k2.LE.Nr ) THEN
                0183 C-      merging former (=k1) into new (=k2) surface grid cell
                0184             uVel(i,j,k2,bi,bj) =
                0185      &            uVel(i,j,k1,bi,bj)*mrgFacW(i,j,1,bi,bj)
                0186      &          + uVel(i,j,k2,bi,bj)*mrgFacW(i,j,2,bi,bj)
                0187             uVel(i,j,k1,bi,bj) = 0.0
                0188 #ifdef ALLOW_ADAMSBASHFORTH_3
                0189             guNm(i,j,k2,bi,bj,1) =
                0190      &            guNm(i,j,k1,bi,bj,1)*mrgFacW(i,j,1,bi,bj)
                0191      &          + guNm(i,j,k2,bi,bj,1)*mrgFacW(i,j,2,bi,bj)
                0192             guNm(i,j,k2,bi,bj,2) =
                0193      &            guNm(i,j,k1,bi,bj,2)*mrgFacW(i,j,1,bi,bj)
                0194      &          + guNm(i,j,k2,bi,bj,2)*mrgFacW(i,j,2,bi,bj)
                0195             guNm(i,j,k1,bi,bj,1)= 0.0
                0196             guNm(i,j,k1,bi,bj,2)= 0.0
                0197 #else
                0198             guNm1(i,j,k2,bi,bj) =
                0199      &            guNm1(i,j,k1,bi,bj)*mrgFacW(i,j,1,bi,bj)
                0200      &          + guNm1(i,j,k2,bi,bj)*mrgFacW(i,j,2,bi,bj)
                0201             guNm1(i,j,k1,bi,bj) = 0.0
                0202 #endif
                0203            ENDIF
                0204           ENDIF
                0205          ENDDO
                0206         ENDDO
                0207         DO j = 1-OLy,sNy+OLy
                0208          DO i = 1-OLx,sNx+OLx
                0209           IF ( k1SurfS(i,j,bi,bj).NE.kSurfS(i,j,bi,bj) ) THEN
                0210            k1 = k1SurfS(i,j,bi,bj)
                0211            k2 =  kSurfS(i,j,bi,bj)
                0212            IF ( k2.LT.k1 .AND. k1.LE.Nr ) THEN
                0213 C-      spliting case: just copy former surface value (=k1) into new one (=k2)
                0214             vVel(i,j,k2,bi,bj) = vVel(i,j,k1,bi,bj)
                0215 #ifdef ALLOW_ADAMSBASHFORTH_3
                0216             gvNm(i,j,k2,bi,bj,1) = gvNm(i,j,k1,bi,bj,1)
                0217             gvNm(i,j,k2,bi,bj,2) = gvNm(i,j,k1,bi,bj,2)
                0218 #else
                0219             gvNm1(i,j,k2,bi,bj) = gvNm1(i,j,k1,bi,bj)
                0220 #endif
                0221            ENDIF
                0222            IF ( k2.GT.k1 .AND. k2.LE.Nr ) THEN
                0223 C-      merging former (=k1) into new (=k2) surface grid cell
                0224             vVel(i,j,k2,bi,bj) =
                0225      &            vVel(i,j,k1,bi,bj)*mrgFacS(i,j,1,bi,bj)
                0226      &          + vVel(i,j,k2,bi,bj)*mrgFacS(i,j,2,bi,bj)
                0227             vVel(i,j,k1,bi,bj) = 0.0
                0228 #ifdef ALLOW_ADAMSBASHFORTH_3
                0229             gvNm(i,j,k2,bi,bj,1) =
                0230      &            gvNm(i,j,k1,bi,bj,1)*mrgFacS(i,j,1,bi,bj)
                0231      &          + gvNm(i,j,k2,bi,bj,1)*mrgFacS(i,j,2,bi,bj)
                0232             gvNm(i,j,k2,bi,bj,2) =
                0233      &            gvNm(i,j,k1,bi,bj,2)*mrgFacS(i,j,1,bi,bj)
                0234      &          + gvNm(i,j,k2,bi,bj,2)*mrgFacS(i,j,2,bi,bj)
                0235             gvNm(i,j,k1,bi,bj,1)= 0.0
                0236             gvNm(i,j,k1,bi,bj,2)= 0.0
                0237 #else
                0238             gvNm1(i,j,k2,bi,bj) =
                0239      &            gvNm1(i,j,k1,bi,bj)*mrgFacS(i,j,1,bi,bj)
                0240      &          + gvNm1(i,j,k2,bi,bj)*mrgFacS(i,j,2,bi,bj)
                0241             gvNm1(i,j,k1,bi,bj) = 0.0
                0242 #endif
                0243            ENDIF
                0244           ENDIF
                0245          ENDDO
                0246         ENDDO
                0247 
                0248        ENDDO
                0249       ENDDO
                0250 
                0251 #endif /* ALLOW_SHELFICE_REMESHING */
                0252       RETURN
                0253       END