pkg/ggl90/ggl90_mixinglength.F

b038e3cc4f Mart*0001 #include "GGL90_OPTIONS.h"
                0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
                0005
                0006 CBOP
                0007 C !ROUTINE: GGL90_MIXINGLENGTH
                0008
                0009 C !INTERFACE: ======================================================
                0010       SUBROUTINE GGL90_MIXINGLENGTH(
                0011      U     GGL90mixingLength,
                0012 #ifdef ALLOW_GGL90_LANGMUIR
                0013      O     LCmixingLength,
                0014 #endif
                0015      O     rMixingLength,
                0016      I     iMin ,iMax ,jMin ,jMax,
                0017      I     bi, bj, myTime, myIter, myThid )
                0018
                0019 C !DESCRIPTION: \bv
                0020 C     *==========================================================*
                0021 C     | SUBROUTINE GGL90_MIXINGLENGTH                            |
                0022 C     | o Compute GGL90mixingLength (and LCmixingLength)         |
                0023 C     *==========================================================*
                0024 C     | Equation numbers refer to                                |
                0025 C     |  Gaspar et al. (1990), JGR 95 (C9), pp 16,179            |
                0026 C     | Some parts of the implementation follow Blanke and       |
                0027 C     |  Delecuse (1993), JPO, and OPA code, in particular the   |
                0028 C     |  computation of the                                      |
                0029 C     |  mixing length = max(min(lk,depth),lkmin)                |
                0030 C     | Note: Only call this S/R if Nr > 1 (no use if Nr=1)      |
                0031 C     *==========================================================*
                0032
                0033 C \ev
                0034
                0035 C !USES: ============================================================
                0036       IMPLICIT NONE
                0037 #include "SIZE.h"
                0038 #include "EEPARAMS.h"
                0039 #include "PARAMS.h"
                0040 #include "GRID.h"
                0041 #include "GGL90.h"
                0042 #ifdef ALLOW_SHELFICE
                0043 # include "SHELFICE.h"
                0044 #endif
                0045 #ifdef ALLOW_AUTODIFF
                0046 # include "AUTODIFF_PARAMS.h"
                0047 #endif
                0048 #ifdef ALLOW_AUTODIFF_TAMC
                0049 # include "tamc.h"
                0050 #endif
                0051
                0052 C !INPUT PARAMETERS: ===================================================
                0053 C Routine arguments
9af873c532 Hajo*0054 C     GGL90mixingLength :: mixing length (m) following Banke+Delecuse
                0055 C     rMixingLength     :: inverse of mixing length
b038e3cc4f Mart*0056 C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
                0057 C     bi, bj :: Current tile indices
                0058 C     myTime :: Current time in simulation
                0059 C     myIter :: Current time-step number
                0060 C     myThid :: My Thread Id number
                0061       _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0062 #ifdef ALLOW_GGL90_LANGMUIR
                0063       _RL     LCmixingLength   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0064 #endif
                0065       _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
9af873c532 Hajo*0066       INTEGER iMin ,iMax ,jMin ,jMax
                0067       INTEGER bi, bj
b038e3cc4f Mart*0068       _RL     myTime
                0069       INTEGER myIter
                0070       INTEGER myThid
                0071
                0072 #ifdef ALLOW_GGL90
                0073 C !LOCAL VARIABLES: ====================================================
                0074 C     i, j, k          :: array computation indices
                0075 C     kSrf             :: vertical index of surface level
                0076 C     kTop             :: index of top interface (just below surf. level)
                0077 C
                0078 C     In general, all 3D variables are defined at W-points (i.e.,
                0079 C     between k and k-1), all 2D variables are also defined at W-points
                0080 C     or at the very surface level (like uStarSquare)
                0081       INTEGER i, j, k
                0082       INTEGER kSrf, kTop
                0083       _RL MaxLength
                0084       _RL mxLength_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0085       _RL MLtmp
                0086 C     This mixed layer model is not invariant under coordinate
                0087 C     transformation to pressure coordinates, so we need these
                0088 C     factors to scale the vertical (pressure) coordinates
                0089       _RL coordFac, recip_coordFac
                0090       INTEGER locMxlMaxFlag
9af873c532 Hajo*0091       CHARACTER*(MAX_LEN_MBUF) msgBuf
b038e3cc4f Mart*0092 #ifdef ALLOW_AUTODIFF_TAMC
                0093 C     tkey :: tape key (depends on tiles)
                0094 C     kkey :: tape key (depends on levels and tiles)
                0095       INTEGER tkey, kkey
                0096 #endif
                0097 CEOP
                0098
                0099       IF ( usingPCoords ) THEN
                0100        kSrf = Nr
                0101        kTop = Nr
                0102       ELSE
                0103        kSrf =  1
                0104        kTop =  2
                0105       ENDIF
                0106
                0107       coordFac = 1. _d 0
                0108       IF ( usingPCoords) coordFac = gravity * rhoConst
                0109       recip_coordFac = 1./coordFac
                0110
                0111       locMxlMaxFlag = mxlMaxFlag
                0112 #ifdef ALLOW_AUTODIFF
                0113       IF ( inAdMode ) locMxlMaxFlag = adMxlMaxFlag
                0114 #endif
                0115 #ifdef ALLOW_AUTODIFF_TAMC
                0116       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0117 #endif /* ALLOW_AUTODIFF_TAMC */
                0118
9af873c532 Hajo*0119 C--   Initialize local fields
b038e3cc4f Mart*0120       DO k=1,Nr
                0121        DO j=1-OLy,sNy+OLy
                0122         DO i=1-OLx,sNx+OLx
                0123          rMixingLength(i,j,k) = 0. _d 0
                0124          mxLength_Dn  (i,j,k) = 0. _d 0
                0125         ENDDO
                0126        ENDDO
                0127       ENDDO
                0128       DO j=1-OLy,sNy+OLy
                0129        DO i=1-OLx,sNx+OLx
                0130 c       rMixingLength(i,j,1)  = 0. _d 0
                0131         mxLength_Dn(i,j,1) = GGL90mixingLengthMin
                0132        ENDDO
                0133       ENDDO
                0134
                0135 #ifdef ALLOW_GGL90_LANGMUIR
                0136       IF (useLANGMUIR) THEN
                0137        DO k=1,Nr
                0138         DO j=1-OLy,sNy+OLy
                0139          DO i=1-OLx,sNx+OLx
                0140           LCmixingLength(i,j,k) = GGL90mixingLengthMin
                0141          ENDDO
                0142         ENDDO
                0143        ENDDO
                0144       ENDIF
                0145 #endif
                0146
9af873c532 Hajo*0147 C-    Ensure mixing between first and second level
b038e3cc4f Mart*0148       IF (mxlSurfFlag) THEN
                0149        DO j=jMin,jMax
                0150         DO i=iMin,iMax
                0151 #ifdef ALLOW_SHELFICE
                0152          IF ( useShelfIce ) THEN
                0153           kSrf = MAX(1,kTopC(i,j,bi,bj))
                0154           kTop = MIN(kSrf+1,Nr)
                0155          ENDIF
                0156 #endif
                0157          GGL90mixingLength(i,j,kTop)=drF(kSrf)*recip_coordFac
                0158         ENDDO
                0159        ENDDO
                0160       ENDIF
                0161
                0162 C--   Impose upper and lower bound for mixing length
                0163 #ifdef ALLOW_AUTODIFF
                0164 CADJ STORE GGL90mixingLength = comlev1_bibj, key=tkey, kind=isbyte
                0165 #endif
                0166       IF ( locMxlMaxFlag .EQ. 0 ) THEN
                0167
                0168        DO k=2,Nr
                0169         DO j=jMin,jMax
                0170          DO i=iMin,iMax
                0171 C     Use thickness of water column (inverse of recip_Rcol) as the
                0172 C     maximum length.
                0173           MaxLength =
                0174      &         ( Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) )*recip_coordFac
                0175           GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
                0176      &                                   MaxLength)
                0177          ENDDO
                0178         ENDDO
                0179        ENDDO
                0180
                0181       ELSEIF ( locMxlMaxFlag .EQ. 1 ) THEN
                0182
                0183        DO k=2,Nr
                0184         DO j=jMin,jMax
                0185          DO i=iMin,iMax
                0186           MaxLength=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
                0187      &         * recip_coordFac
                0188 c         MaxLength=MAX(MaxLength,20. _d 0)
                0189           GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
                0190      &                                   MaxLength)
                0191          ENDDO
                0192         ENDDO
                0193        ENDDO
                0194
                0195       ELSEIF ( locMxlMaxFlag .EQ. 2 .OR. locMxlMaxFlag .EQ. 3 ) THEN
                0196
                0197        IF ( usingPcoords ) THEN
                0198 #ifdef ALLOW_AUTODIFF_TAMC
                0199 CADJ STORE GGL90mixingLength = comlev1_bibj, key = tkey, kind=isbyte
                0200 #endif /* ALLOW_AUTODIFF_TAMC */
                0201 C     Downward sweep, extra treatment of k=Nr for p-coordinates
                0202 C     because level Nr+1 is not available
                0203         DO j=jMin,jMax
                0204          DO i=iMin,iMax
                0205           mxLength_Dn(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
                0206      &         GGL90mixingLengthMin+drF(Nr)*recip_coordFac)
                0207          ENDDO
                0208         ENDDO
                0209         DO k=Nr-1,2,-1
                0210 #ifdef ALLOW_AUTODIFF_TAMC
                0211          kkey = k + (tkey-1)*Nr
                0212 CADJ STORE mxLength_Dn(:,:,k+1)
                0213 CADJ &     = comlev1_bibj_k, key = kkey, kind=isbyte
                0214 #endif /* ALLOW_AUTODIFF_TAMC */
                0215          DO j=jMin,jMax
                0216           DO i=iMin,iMax
                0217            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
                0218      &          mxLength_Dn(i,j,k+1)+drF(k)*recip_coordFac)
                0219           ENDDO
                0220          ENDDO
                0221         ENDDO
                0222 C     Upward sweep
                0223         DO k=2,Nr
                0224 #ifdef ALLOW_AUTODIFF_TAMC
                0225          kkey = k + (tkey-1)*Nr
                0226 C     It is important that the two k-levels of these fields are stored
                0227 C     in one statement because otherwise taf will only store one, which
                0228 C     is wrong (i.e. was wrong in previous versions).
                0229 CADJ STORE GGL90mixingLength(:,:,k-1), GGL90mixingLength(:,:,k)
                0230 CADJ &     = comlev1_bibj_k, key = kkey, kind=isbyte
                0231 #endif /* ALLOW_AUTODIFF_TAMC */
                0232          DO j=jMin,jMax
                0233           DO i=iMin,iMax
                0234            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
                0235      &          GGL90mixingLength(i,j,k-1)+drF(k-1)*recip_coordFac)
                0236           ENDDO
                0237          ENDDO
                0238         ENDDO
9af873c532 Hajo*0239
b038e3cc4f Mart*0240        ELSE
9af873c532 Hajo*0241 C-    Z-coordinate case
b038e3cc4f Mart*0242 #ifdef ALLOW_AUTODIFF_TAMC
                0243 CADJ STORE GGL90mixingLength = comlev1_bibj, key = tkey, kind=isbyte
                0244 #endif /* ALLOW_AUTODIFF_TAMC */
                0245 C     Downward sweep
                0246         DO k=2,Nr
                0247 #ifdef ALLOW_AUTODIFF_TAMC
                0248          kkey = k + (tkey-1)*Nr
                0249 CADJ STORE mxLength_Dn(:,:,k-1)
                0250 CADJ &     = comlev1_bibj_k, key = kkey, kind=isbyte
                0251 #endif /* ALLOW_AUTODIFF_TAMC */
                0252          DO j=jMin,jMax
                0253           DO i=iMin,iMax
                0254            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
                0255      &          mxLength_Dn(i,j,k-1)+drF(k-1)*recip_coordFac)
                0256           ENDDO
                0257          ENDDO
                0258         ENDDO
9af873c532 Hajo*0259
b038e3cc4f Mart*0260 C     Upward sweep, extra treatment of k=Nr for z-coordinates
                0261 C     because level Nr+1 is not available
                0262         DO j=jMin,jMax
                0263          DO i=iMin,iMax
                0264           GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
                0265      &         GGL90mixingLengthMin+drF(Nr)*recip_coordFac)
                0266          ENDDO
                0267         ENDDO
                0268         DO k=Nr-1,2,-1
                0269 #ifdef ALLOW_AUTODIFF_TAMC
                0270          kkey = k + (tkey-1)*Nr
                0271 C     It is important that the two k-levels of these fields are stored
                0272 C     in one statement because otherwise taf will only store one, which
                0273 C     is wrong (i.e. was wrong in previous versions).
                0274 CADJ STORE GGL90mixingLength(:,:,k+1), GGL90mixingLength(:,:,k)
                0275 CADJ &     = comlev1_bibj_k, key = kkey, kind=isbyte
                0276 #endif /* ALLOW_AUTODIFF_TAMC */
                0277          DO j=jMin,jMax
                0278           DO i=iMin,iMax
                0279            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
                0280      &          GGL90mixingLength(i,j,k+1)+drF(k)*recip_coordFac)
                0281           ENDDO
                0282          ENDDO
                0283         ENDDO
9af873c532 Hajo*0284 C-    end if P/Z-coordinate
b038e3cc4f Mart*0285        ENDIF
9af873c532 Hajo*0286
b038e3cc4f Mart*0287 #ifdef ALLOW_AUTODIFF_TAMC
                0288 CADJ STORE mxLength_Dn       = comlev1_bibj, key = tkey, kind=isbyte
                0289 CADJ STORE GGL90mixingLength = comlev1_bibj, key = tkey, kind=isbyte
                0290 #endif
                0291 C     Impose minimum from downward sweep
                0292        DO k=2,Nr
                0293         DO j=jMin,jMax
                0294          DO i=iMin,iMax
                0295           GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
                0296      &                                  mxLength_Dn(i,j,k))
                0297          ENDDO
                0298         ENDDO
                0299        ENDDO
                0300
                0301       ELSE
                0302        WRITE(msgBuf,'(A,I5,A)')
                0303      &   'GGL90_MIXINGLENGTH: mxlMaxFlag=',
                0304      &   locMxlMaxFlag,' not implemented'
                0305        CALL PRINT_ERROR( msgBuf, myThid )
                0306        STOP 'ABNORMAL END: S/R GGL90_MIXINGLENGTH'
                0307       ENDIF
                0308
9af873c532 Hajo*0309 #ifdef ALLOW_GGL90_LANGMUIR
                0310 C----------------------------------
                0311 C--   Langmuir circulation effect :
                0312 C----------------------------------
                0313       IF (useLANGMUIR) THEN
                0314 # ifdef ALLOW_AUTODIFF_TAMC
                0315 CADJ STORE GGL90mixingLength = comlev1_bibj, key = tkey, kind=isbyte
                0316 # endif
                0317        IF ( locMxlMaxFlag .EQ. 1 ) THEN
                0318         DO k=2,Nr
                0319          DO j=jMin,jMax
                0320           DO i=iMin,iMax
                0321            IF ( usingPcoords ) THEN
                0322             MaxLength=(rF(k)-R_low(i,j,bi,bj))*recip_coordFac
                0323            ELSE
                0324             MaxLength=(Ro_surf(i,j,bi,bj)-rF(k)) * recip_coordFac
                0325            ENDIF
                0326            IF (GGL90mixingLength(i,j,k) .EQ. MaxLength) THEN
                0327             LCmixingLength(i,j,k) = LC_Gamma * GGL90mixingLength(i,j,k)
                0328            ELSE
                0329             LCmixingLength(i,j,k) = GGL90mixingLength(i,j,k)
                0330            ENDIF
                0331           ENDDO
                0332          ENDDO
                0333         ENDDO
                0334
                0335        ELSEIF ( locMxlMaxFlag .EQ. 2 .OR. locMxlMaxFlag .EQ. 3 ) THEN
                0336
                0337 # ifdef ALLOW_AUTODIFF_TAMC
                0338 CADJ STORE mxLength_Dn       = comlev1_bibj, key = tkey, kind=isbyte
                0339 # endif
                0340         DO k=2,Nr
                0341          DO j=jMin,jMax
                0342           DO i=iMin,iMax
                0343            IF (GGL90mixingLength(i,j,k) .EQ. mxLength_Dn(i,j,k)) THEN
                0344             LCmixingLength(i,j,k) = LC_Gamma * GGL90mixingLength(i,j,k)
                0345            ELSE
                0346             LCmixingLength(i,j,k) = GGL90mixingLength(i,j,k)
                0347            ENDIF
                0348           ENDDO
                0349          ENDDO
                0350         ENDDO
                0351
                0352        ELSE
                0353         WRITE(msgBuf,'(2A,I5,A)')
                0354      &       'GGL90_MIXINGLENGTH: ',
                0355      &       'Langmuir Circ. Parameterization with mxlMaxFlag=',
                0356      &   locMxlMaxFlag,' not implemented'
                0357         CALL PRINT_ERROR( msgBuf, myThid )
                0358         STOP 'ABNORMAL END: S/R GGL90_MIXINGLENGTH'
                0359        ENDIF
                0360
                0361 C--   Impose minimum LC-mixing length to avoid division by zero
                0362        IF ( locMxlMaxFlag .EQ. 1 .OR. locMxlMaxFlag .EQ. 2 ) THEN
                0363 # ifdef ALLOW_AUTODIFF_TAMC
                0364 CADJ STORE LCmixingLength    = comlev1_bibj, key = tkey, kind=isbyte
                0365 # endif
                0366         DO k=2,Nr
                0367          DO j=jMin,jMax
                0368           DO i=iMin,iMax
                0369            MLtmp = MAX(LCmixingLength(i,j,k),GGL90mixingLengthMin)
                0370            LCmixingLength(i,j,k) = MLtmp
                0371           ENDDO
                0372          ENDDO
                0373         ENDDO
                0374        ENDIF
                0375       ENDIF
                0376 #endif /* ALLOW_GGL90_LANGMUIR */
                0377
b038e3cc4f Mart*0378 #ifdef ALLOW_AUTODIFF_TAMC
                0379 CADJ STORE GGL90mixingLength = comlev1_bibj, key = tkey, kind=isbyte
                0380 #endif
                0381 C--   Impose minimum mixing length to avoid division by zero
                0382 C     and compute inverse
                0383       IF ( locMxlMaxFlag.EQ.3 ) THEN
                0384 #ifdef ALLOW_AUTODIFF_TAMC
                0385 CADJ STORE mxLength_Dn       = comlev1_bibj, key = tkey, kind=isbyte
                0386 #endif
                0387        DO k=2,Nr
                0388         DO j=jMin,jMax
                0389          DO i=iMin,iMax
                0390 #ifdef GGL90_REGULARIZE_MIXINGLENGTH
                0391           MLtmp = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k)
                0392      &         + GGL90mixingLengthMin**2 )
                0393 #else
                0394           MLtmp = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
                0395           MLtmp = MAX( MLtmp, GGL90mixingLengthMin )
                0396 #endif
                0397           rMixingLength(i,j,k) = 1. _d 0 / MLtmp
                0398          ENDDO
                0399         ENDDO
                0400        ENDDO
                0401       ELSE
                0402        DO k=2,Nr
                0403         DO j=jMin,jMax
                0404          DO i=iMin,iMax
                0405 #ifdef GGL90_REGULARIZE_MIXINGLENGTH
                0406           MLtmp = SQRT( GGL90mixingLength(i,j,k)**2
                0407      &         + GGL90mixingLengthMin**2 )
                0408 #else
                0409           MLtmp = MAX(GGL90mixingLength(i,j,k),GGL90mixingLengthMin)
                0410 #endif
                0411           GGL90mixingLength(i,j,k) = MLtmp
                0412           rMixingLength(i,j,k) = 1. _d 0 / MLtmp
                0413          ENDDO
                0414         ENDDO
                0415        ENDDO
                0416       ENDIF
                0417
                0418 #endif /* ALLOW_GGL90 */
                0419
                0420       RETURN
                0421       END