Back to home page

MITgcm

 
 

    


File indexing completed on 2023-09-03 05:09:51 UTC

view on githubraw file Latest commit 74487008 on 2023-09-03 01:50:18 UTC
6d54cf9ca1 Ed H*0001 #include "PACKAGES_CONFIG.h"
5646dbad6f Patr*0002 #include "CPP_OPTIONS.h"
7448700841 Mart*0003 #ifdef ALLOW_AUTODIFF
                0004 # include "AUTODIFF_OPTIONS.h"
                0005 #endif
5646dbad6f Patr*0006 
9366854e02 Chri*0007 CBOP
                0008 C     !ROUTINE: CALC_VISCOSITY
                0009 C     !INTERFACE:
64f8465da1 Jean*0010       SUBROUTINE CALC_VISCOSITY(
                0011      I           bi,bj, iMin,iMax,jMin,jMax,
efd2c352d2 Jean*0012      O           kappaRU, kappaRV,
64f8465da1 Jean*0013      I           myThid )
5646dbad6f Patr*0014 
9366854e02 Chri*0015 C     !DESCRIPTION: \bv
                0016 C     *==========================================================*
64f8465da1 Jean*0017 C     | SUBROUTINE CALC_VISCOSITY
                0018 C     | o Calculate net vertical viscosity
9366854e02 Chri*0019 C     *==========================================================*
                0020 C     \ev
5646dbad6f Patr*0021 
9366854e02 Chri*0022 C     !USES:
                0023       IMPLICIT NONE
5646dbad6f Patr*0024 C     == GLobal variables ==
                0025 #include "SIZE.h"
                0026 #include "EEPARAMS.h"
                0027 #include "PARAMS.h"
                0028 #include "DYNVARS.h"
                0029 #include "GRID.h"
7448700841 Mart*0030 #ifndef EXCLUDE_PCELL_MIX_CODE
                0031 # ifdef ALLOW_AUTODIFF
                0032 #  include "tamc.h"
                0033 # endif
                0034 #endif
5646dbad6f Patr*0035 
9366854e02 Chri*0036 C     !INPUT/OUTPUT PARAMETERS:
5646dbad6f Patr*0037 C     == Routine arguments ==
64f8465da1 Jean*0038 C     iMin,iMax,jMin,jMax :: Range of points for which calculation
                0039 C     bi,bj   :: current tile indices
efd2c352d2 Jean*0040 C     kappaRU :: Total vertical viscosity for zonal flow.
                0041 C     kappaRV :: Total vertical viscosity for meridional flow.
64f8465da1 Jean*0042 C     myThid  :: my Thread Id number
                0043       INTEGER iMin,iMax,jMin,jMax
                0044       INTEGER bi,bj
efd2c352d2 Jean*0045       _RL kappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0046       _RL kappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
5646dbad6f Patr*0047       INTEGER myThid
                0048 
9366854e02 Chri*0049 C     !LOCAL VARIABLES:
5646dbad6f Patr*0050 C     == Local variables ==
64f8465da1 Jean*0051 C     i, j, k :: Loop counters
                0052       INTEGER i,j,k
efd2c352d2 Jean*0053       INTEGER ki
2d435b47ac Jean*0054 #ifndef EXCLUDE_PCELL_MIX_CODE
                0055       INTEGER km, mixSurf, mixBott
                0056       _RL pC_kFac
                0057       _RL tmpFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
7448700841 Mart*0058 # ifdef ALLOW_AUTODIFF_TAMC
                0059 C     tkey :: tape key (tile dependent)
                0060 C     kkey :: tape key (level and tile dependent)
                0061       INTEGER tkey, kkey
                0062 # endif
2d435b47ac Jean*0063 #endif
9366854e02 Chri*0064 CEOP
5646dbad6f Patr*0065 
efd2c352d2 Jean*0066       DO k = 1,Nr+1
                0067        ki = MIN(k,Nr)
5646dbad6f Patr*0068 
d8d1486ca1 Jean*0069        DO j = 1-OLy, sNy+OLy
                0070         DO i = 1-OLx, sNx+OLx
efd2c352d2 Jean*0071          kappaRU(i,j,k) = viscArNr(ki)
                0072          kappaRV(i,j,k) = viscArNr(ki)
64f8465da1 Jean*0073         ENDDO
5646dbad6f Patr*0074        ENDDO
                0075 
                0076 #ifdef ALLOW_KPP
efd2c352d2 Jean*0077        IF ( useKPP .AND. k.LE.Nr ) THEN
                0078          CALL KPP_CALC_VISC(
64f8465da1 Jean*0079      I        bi,bj, iMin,iMax,jMin,jMax, k,
efd2c352d2 Jean*0080      O        kappaRU, kappaRV,
5646dbad6f Patr*0081      I        myThid)
efd2c352d2 Jean*0082        ENDIF
5646dbad6f Patr*0083 #endif
                0084 
e864122ae8 Mart*0085 #ifdef ALLOW_PP81
efd2c352d2 Jean*0086        IF ( usePP81 .AND. k.LE.Nr ) THEN
                0087          CALL PP81_CALC_VISC(
64f8465da1 Jean*0088      I        bi,bj, iMin,iMax,jMin,jMax, k,
efd2c352d2 Jean*0089      O        kappaRU, kappaRV,
e864122ae8 Mart*0090      I        myThid)
efd2c352d2 Jean*0091        ENDIF
e864122ae8 Mart*0092 #endif
                0093 
d8d1486ca1 Jean*0094 #ifdef ALLOW_KL10
efd2c352d2 Jean*0095        IF ( useKL10 .AND. k.LE.Nr ) THEN
                0096          CALL KL10_CALC_VISC(
d8d1486ca1 Jean*0097      I        bi,bj, iMin,iMax,jMin,jMax, k,
efd2c352d2 Jean*0098      O        kappaRU, kappaRV,
d8d1486ca1 Jean*0099      I        myThid)
efd2c352d2 Jean*0100        ENDIF
d8d1486ca1 Jean*0101 #endif
                0102 
e864122ae8 Mart*0103 #ifdef ALLOW_MY82
efd2c352d2 Jean*0104        IF ( useMY82 .AND. k.LE.Nr ) THEN
                0105          CALL MY82_CALC_VISC(
64f8465da1 Jean*0106      I        bi,bj, iMin,iMax,jMin,jMax, k,
efd2c352d2 Jean*0107      O        kappaRU, kappaRV,
e864122ae8 Mart*0108      I        myThid)
efd2c352d2 Jean*0109        ENDIF
e864122ae8 Mart*0110 #endif
                0111 
69a7b27187 Mart*0112 #ifdef ALLOW_GGL90
efd2c352d2 Jean*0113        IF ( useGGL90 .AND. k.LE.Nr ) THEN
                0114          CALL GGL90_CALC_VISC(
64f8465da1 Jean*0115      I        bi,bj, iMin,iMax,jMin,jMax, k,
efd2c352d2 Jean*0116      O        kappaRU, kappaRV,
69a7b27187 Mart*0117      I        myThid)
efd2c352d2 Jean*0118        ENDIF
69a7b27187 Mart*0119 #endif
                0120 
efd2c352d2 Jean*0121        IF ( k.EQ.Nr+1 .AND.
                0122      &     ( usePP81 .OR. useKL10 .OR. useMY82 .OR. useGGL90 )
                0123      &    ) THEN
                0124         DO j = 1-OLy, sNy+OLy
                0125          DO i = 1-OLx, sNx+OLx
                0126           kappaRU(i,j,k) = kappaRU(i,j,ki)
                0127           kappaRV(i,j,k) = kappaRV(i,j,ki)
                0128          ENDDO
                0129         ENDDO
                0130        ENDIF
                0131 
64f8465da1 Jean*0132 C--   end of k loop
                0133       ENDDO
5646dbad6f Patr*0134 
2d435b47ac Jean*0135 #ifndef EXCLUDE_PCELL_MIX_CODE
7448700841 Mart*0136 # ifdef ALLOW_AUTODIFF_TAMC
                0137       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0138 CADJ STORE kappaRU = comlev1_bibj, key = tkey, kind = isbyte
                0139 CADJ STORE kappaRV = comlev1_bibj, key = tkey, kind = isbyte
                0140 # endif
2d435b47ac Jean*0141       IF ( interViscAr_pCell ) THEN
                0142 C--   This is a hack: alter vertical viscosity (instead of changing many S/R)
                0143 C     in order to account for missing hFac in viscous term
                0144        DO k = 2,Nr
                0145          km = k - 1
                0146 C-    account for true distance (including hFac) in vertical gradient
                0147          DO j = 2-OLy, sNy+OLy
                0148           DO i = 2-OLx, sNx+OLx
                0149            IF ( k.GT.kSurfW(i,j,bi,bj) .AND.
                0150      &          k.LE.MIN( kLowC(i,j,bi,bj), kLowC(i-1,j,bi,bj) )
                0151      &        ) THEN
                0152              kappaRU(i,j,k) = kappaRU(i,j,k)
                0153      &                *twoRL/(hFacW(i,j,km,bi,bj)+hFacW(i,j,k,bi,bj))
                0154            ENDIF
                0155           ENDDO
                0156          ENDDO
                0157          DO j = 2-OLy, sNy+OLy
                0158           DO i = 2-OLx, sNx+OLx
                0159            IF ( k.GT.kSurfS(i,j,bi,bj) .AND.
                0160      &          k.LE.MIN( kLowC(i,j,bi,bj), kLowC(i,j-1,bi,bj) )
                0161      &        ) THEN
                0162              kappaRV(i,j,k) = kappaRV(i,j,k)
                0163      &                *twoRL/(hFacS(i,j,km,bi,bj)+hFacS(i,j,k,bi,bj))
                0164            ENDIF
                0165           ENDDO
                0166          ENDDO
                0167        ENDDO
                0168       ENDIF
                0169 
                0170       IF ( pCellMix_select.GT.0 ) THEN
                0171 C--   This is a hack: alter vertical viscosity (instead of changing many S/R)
                0172 C     in order to to increase mixing for too thin surface/bottom partial cell
                0173        mixSurf = pCellMix_select/10
                0174        mixBott = MOD(pCellMix_select,10)
                0175        DO k = 2,Nr
7448700841 Mart*0176 # ifdef ALLOW_AUTODIFF_TAMC
                0177          kkey = k + (tkey-1)*Nr
                0178 # endif
2d435b47ac Jean*0179          km = k - 1
                0180          pC_kFac = 1.
                0181          IF ( pCellMix_delR.LT.drF(k) )
                0182      &     pC_kFac = pCellMix_delR*recip_drF(k)
                0183 
                0184 C-    Increase KappaRU above bottom level:
                0185          IF ( mixBott.GE.1 ) THEN
                0186           DO j = 2-OLy, sNy+OLy
                0187            DO i = 2-OLx, sNx+OLx
                0188              tmpFac(i,j) = 0. _d 0
                0189              IF ( k.EQ.MIN( kLowC(i,j,bi,bj), kLowC(i-1,j,bi,bj) )
                0190      &        .AND. k.GT.kSurfW(i,j,bi,bj) ) THEN
                0191                tmpFac(i,j) = pC_kFac*_recip_hFacW(i,j,k,bi,bj)
                0192              ENDIF
                0193            ENDDO
                0194           ENDDO
f2a88c9ff8 jm-c 0195           IF ( mixBott.EQ.2 ) THEN
                0196            DO j = 2-OLy, sNy+OLy
                0197             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0198              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0199             ENDDO
2d435b47ac Jean*0200            ENDDO
f2a88c9ff8 jm-c 0201           ELSEIF ( mixBott.EQ.3 ) THEN
                0202            DO j = 2-OLy, sNy+OLy
                0203             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0204              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0205             ENDDO
2d435b47ac Jean*0206            ENDDO
f2a88c9ff8 jm-c 0207           ELSEIF ( mixBott.EQ.4 ) THEN
                0208            DO j = 2-OLy, sNy+OLy
                0209             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0210              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
                0211      &                    *tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0212             ENDDO
2d435b47ac Jean*0213            ENDDO
f2a88c9ff8 jm-c 0214           ENDIF
2d435b47ac Jean*0215 C-    increase mixing above bottom (by ~(1/hFac)^mixBott) if too thin p-cell
                0216           DO j = 2-OLy, sNy+OLy
                0217            DO i = 2-OLx, sNx+OLx
                0218              tmpFac(i,j) = MIN( tmpFac(i,j), pCellMix_maxFac )
7448700841 Mart*0219 # ifdef ALLOW_AUTODIFF_TAMC
                0220            ENDDO
                0221           ENDDO
                0222 CADJ STORE tmpFac         = comlev1_bibj_k, key = kkey, kind = isbyte
                0223 CADJ STORE kappaRU(:,:,k) = comlev1_bibj_k, key = kkey, kind = isbyte
                0224           DO j = 2-OLy, sNy+OLy
                0225            DO i = 2-OLx, sNx+OLx
                0226 # endif
2d435b47ac Jean*0227              kappaRU(i,j,k) = MAX( kappaRU(i,j,k),
                0228      &                             pCellMix_viscAr(k)*tmpFac(i,j) )
                0229            ENDDO
                0230           ENDDO
                0231          ENDIF
                0232 
                0233 C-    Increase KappaRV above bottom level:
                0234          IF ( mixBott.GE.1 ) THEN
                0235           DO j = 2-OLy, sNy+OLy
                0236            DO i = 2-OLx, sNx+OLx
                0237              tmpFac(i,j) = 0. _d 0
                0238              IF ( k.EQ.MIN( kLowC(i,j,bi,bj), kLowC(i,j-1,bi,bj) )
                0239      &        .AND. k.GT.kSurfS(i,j,bi,bj) ) THEN
                0240                tmpFac(i,j) = pC_kFac*_recip_hFacS(i,j,k,bi,bj)
                0241              ENDIF
                0242            ENDDO
                0243           ENDDO
f2a88c9ff8 jm-c 0244           IF ( mixBott.EQ.2 ) THEN
                0245            DO j = 2-OLy, sNy+OLy
                0246             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0247              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0248             ENDDO
2d435b47ac Jean*0249            ENDDO
f2a88c9ff8 jm-c 0250           ELSEIF ( mixBott.EQ.3 ) THEN
                0251            DO j = 2-OLy, sNy+OLy
                0252             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0253              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0254             ENDDO
2d435b47ac Jean*0255            ENDDO
f2a88c9ff8 jm-c 0256           ELSEIF ( mixBott.EQ.4 ) THEN
                0257            DO j = 2-OLy, sNy+OLy
                0258             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0259              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
                0260      &                    *tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0261             ENDDO
2d435b47ac Jean*0262            ENDDO
f2a88c9ff8 jm-c 0263           ENDIF
2d435b47ac Jean*0264 C-    increase mixing above bottom (by ~(1/hFac)^mixBott) if too thin p-cell
                0265           DO j = 2-OLy, sNy+OLy
                0266            DO i = 2-OLx, sNx+OLx
                0267              tmpFac(i,j) = MIN( tmpFac(i,j), pCellMix_maxFac )
7448700841 Mart*0268 # ifdef ALLOW_AUTODIFF_TAMC
                0269            ENDDO
                0270           ENDDO
                0271 CADJ STORE tmpFac         = comlev1_bibj_k, key = kkey, kind = isbyte
                0272 CADJ STORE kappaRV(:,:,k) = comlev1_bibj_k, key = kkey, kind = isbyte
                0273           DO j = 2-OLy, sNy+OLy
                0274            DO i = 2-OLx, sNx+OLx
                0275 # endif
2d435b47ac Jean*0276              kappaRV(i,j,k) = MAX( kappaRV(i,j,k),
                0277      &                             pCellMix_viscAr(k)*tmpFac(i,j) )
                0278            ENDDO
                0279           ENDDO
                0280          ENDIF
                0281 
                0282          pC_kFac = 1.
                0283          IF ( pCellMix_delR.LT.drF(km) )
                0284      &     pC_kFac = pCellMix_delR*recip_drF(km)
                0285 
                0286 C-    Increase KappaRU below surface level:
                0287          IF ( mixSurf.GE.1 ) THEN
                0288           DO j = 2-OLy, sNy+OLy
                0289            DO i = 2-OLx, sNx+OLx
                0290              tmpFac(i,j) = 0. _d 0
                0291              IF ( km.EQ.kSurfW(i,j,bi,bj) .AND.
                0292      &            km.LT.MIN( kLowC(i,j,bi,bj), kLowC(i-1,j,bi,bj) )
                0293      &          ) THEN
                0294                tmpFac(i,j) = pC_kFac*_recip_hFacW(i,j,km,bi,bj)
                0295              ENDIF
                0296            ENDDO
                0297           ENDDO
f2a88c9ff8 jm-c 0298           IF ( mixSurf.EQ.2 ) THEN
                0299            DO j = 2-OLy, sNy+OLy
                0300             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0301              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0302             ENDDO
2d435b47ac Jean*0303            ENDDO
f2a88c9ff8 jm-c 0304           ELSEIF ( mixSurf.EQ.3 ) THEN
                0305            DO j = 2-OLy, sNy+OLy
                0306             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0307              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0308             ENDDO
2d435b47ac Jean*0309            ENDDO
f2a88c9ff8 jm-c 0310           ELSEIF ( mixSurf.EQ.4 ) THEN
                0311            DO j = 2-OLy, sNy+OLy
                0312             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0313              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
                0314      &                    *tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0315             ENDDO
2d435b47ac Jean*0316            ENDDO
f2a88c9ff8 jm-c 0317           ENDIF
2d435b47ac Jean*0318 C-    increase mixing below surface (by ~(1/hFac)^mixSurf) if too thin p-cell
7448700841 Mart*0319 # ifdef ALLOW_AUTODIFF_TAMC
                0320 CADJ STORE tmpFac = comlev1_bibj_k, key = kkey, kind = isbyte
                0321 # endif
2d435b47ac Jean*0322           DO j = 2-OLy, sNy+OLy
                0323            DO i = 2-OLx, sNx+OLx
                0324              tmpFac(i,j) = MIN( tmpFac(i,j), pCellMix_maxFac )
7448700841 Mart*0325 # ifdef ALLOW_AUTODIFF_TAMC
                0326            ENDDO
                0327           ENDDO
                0328 CADJ STORE tmpFac         = comlev1_bibj_k, key = kkey, kind = isbyte
                0329 CADJ STORE kappaRU(:,:,k) = comlev1_bibj_k, key = kkey, kind = isbyte
                0330           DO j = 2-OLy, sNy+OLy
                0331            DO i = 2-OLx, sNx+OLx
                0332 # endif
2d435b47ac Jean*0333              kappaRU(i,j,k) = MAX( kappaRU(i,j,k),
                0334      &                             pCellMix_viscAr(k)*tmpFac(i,j) )
                0335            ENDDO
                0336           ENDDO
                0337          ENDIF
                0338 
                0339 C-    Increase KappaRV below surface level:
                0340          IF ( mixSurf.GE.1 ) THEN
                0341           DO j = 2-OLy, sNy+OLy
                0342            DO i = 2-OLx, sNx+OLx
                0343              tmpFac(i,j) = 0. _d 0
                0344              IF ( km.EQ.kSurfS(i,j,bi,bj) .AND.
                0345      &            km.LT.MIN( kLowC(i,j,bi,bj), kLowC(i,j-1,bi,bj) )
                0346      &          ) THEN
                0347                tmpFac(i,j) = pC_kFac*_recip_hFacS(i,j,km,bi,bj)
                0348              ENDIF
                0349            ENDDO
                0350           ENDDO
f2a88c9ff8 jm-c 0351           IF ( mixSurf.EQ.2 ) THEN
                0352            DO j = 2-OLy, sNy+OLy
                0353             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0354              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0355             ENDDO
2d435b47ac Jean*0356            ENDDO
f2a88c9ff8 jm-c 0357           ELSEIF ( mixSurf.EQ.3 ) THEN
                0358            DO j = 2-OLy, sNy+OLy
                0359             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0360              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0361             ENDDO
2d435b47ac Jean*0362            ENDDO
f2a88c9ff8 jm-c 0363           ELSEIF ( mixSurf.EQ.4 ) THEN
                0364            DO j = 2-OLy, sNy+OLy
                0365             DO i = 2-OLx, sNx+OLx
2d435b47ac Jean*0366              tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
                0367      &                    *tmpFac(i,j)*tmpFac(i,j)
f2a88c9ff8 jm-c 0368             ENDDO
2d435b47ac Jean*0369            ENDDO
f2a88c9ff8 jm-c 0370           ENDIF
2d435b47ac Jean*0371 C-    increase mixing below surface (by ~(1/hFac)^mixSurf) if too thin p-cell
7448700841 Mart*0372 # ifdef ALLOW_AUTODIFF_TAMC
                0373 CADJ STORE tmpFac = comlev1_bibj_k, key = kkey, kind = isbyte
                0374 # endif
2d435b47ac Jean*0375           DO j = 2-OLy, sNy+OLy
                0376            DO i = 2-OLx, sNx+OLx
                0377              tmpFac(i,j) = MIN( tmpFac(i,j), pCellMix_maxFac )
7448700841 Mart*0378 # ifdef ALLOW_AUTODIFF_TAMC
                0379            ENDDO
                0380           ENDDO
                0381 CADJ STORE tmpFac         = comlev1_bibj_k, key = kkey, kind = isbyte
                0382 CADJ STORE kappaRV(:,:,k) = comlev1_bibj_k, key = kkey, kind = isbyte
                0383           DO j = 2-OLy, sNy+OLy
                0384            DO i = 2-OLx, sNx+OLx
                0385 # endif
2d435b47ac Jean*0386              kappaRV(i,j,k) = MAX( kappaRV(i,j,k),
                0387      &                             pCellMix_viscAr(k)*tmpFac(i,j) )
                0388            ENDDO
                0389           ENDDO
                0390          ENDIF
                0391 
                0392 C--   end of k loop
                0393        ENDDO
                0394       ENDIF
                0395 #endif /* ndef EXCLUDE_PCELL_MIX_CODE */
                0396 
5646dbad6f Patr*0397       RETURN
                0398       END