model/src/impldiff.F

73b66b887d Jean*0001 #include "PACKAGES_CONFIG.h"
1dbaea09ee Chri*0002 #include "CPP_OPTIONS.h"
779cd6d73d Alis*0003
9366854e02 Chri*0004 CBOP
                0005 C     !ROUTINE: IMPLDIFF
                0006 C     !INTERFACE:
779cd6d73d Alis*0007       SUBROUTINE IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
73b66b887d Jean*0008      I                     tracerId, KappaRX, recip_hFac,
d64c4d306c Jean*0009      U                     gTracer,
779cd6d73d Alis*0010      I                     myThid )
9366854e02 Chri*0011 C     !DESCRIPTION: \bv
                0012 C     *==========================================================*
d64c4d306c Jean*0013 C     | S/R IMPLDIFF
                0014 C     | o Solve implicit diffusion equation for vertical
                0015 C     |   diffusivity.
9366854e02 Chri*0016 C     *==========================================================*
d64c4d306c Jean*0017 C     | o Recoded from 2d intermediate fields to 3d to reduce
b7b61e618a Mart*0018 C     |   TAF storage
d64c4d306c Jean*0019 C     | o Fixed missing masks for fields a(), c()
9366854e02 Chri*0020 C     *==========================================================*
                0021 C     \ev
                0022
                0023 C     !USES:
a0b25fcf44 Chri*0024       IMPLICIT NONE
                0025 C     == Global data ==
779cd6d73d Alis*0026 #include "SIZE.h"
                0027 #include "DYNVARS.h"
81bc00c2f0 Chri*0028 #include "EEPARAMS.h"
779cd6d73d Alis*0029 #include "PARAMS.h"
                0030 #include "GRID.h"
4e66ab0b67 Oliv*0031 #ifdef ALLOW_GENERIC_ADVDIFF
02d90fb24c Jean*0032 # include "GAD.h"
4e66ab0b67 Oliv*0033 #endif
                0034 #ifdef ALLOW_LONGSTEP
02d90fb24c Jean*0035 # include "LONGSTEP_PARAMS.h"
4e66ab0b67 Oliv*0036 #endif
                0037 #ifdef ALLOW_PTRACERS
02d90fb24c Jean*0038 # include "PTRACERS_SIZE.h"
                0039 # include "PTRACERS_PARAMS.h"
8adf9f02ba Patr*0040 #endif
9b89fcf692 antn*0041 #ifdef ALLOW_LAYERS
                0042 # include "LAYERS_P2SHARE.h"
                0043 #endif
02d90fb24c Jean*0044 c#ifdef ALLOW_AUTODIFF_TAMC
                0045 c#endif
8adf9f02ba Patr*0046
9366854e02 Chri*0047 C     !INPUT/OUTPUT PARAMETERS:
d64c4d306c Jean*0048 C     tracerId   :: tracer Identificator (if > 0) ; = -1 or -2 when
                0049 C                   solving vertical viscosity implicitly for U or V
                0050 C     KappaRk    :: vertical diffusion coefficient
                0051 C     recip_hFac :: Inverse of cell open-depth factor
                0052 C     gTracer    :: future tracer field
779cd6d73d Alis*0053       INTEGER bi,bj,iMin,iMax,jMin,jMax
73b66b887d Jean*0054       INTEGER tracerId
d64c4d306c Jean*0055       _RL KappaRX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0056       _RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
23a7f3050f Jean*0057       _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
779cd6d73d Alis*0058       INTEGER myThid
a0b25fcf44 Chri*0059
9366854e02 Chri*0060 C     !LOCAL VARIABLES:
779cd6d73d Alis*0061       INTEGER i,j,k
698b6992ee Jean*0062       _RL deltaTX(Nr), locUpdate
d64c4d306c Jean*0063       _RL locTr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0064       _RL a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0065       _RL b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0066       _RL c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0067       _RL bet(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0068       _RL gam(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
73b66b887d Jean*0069 #ifdef ALLOW_DIAGNOSTICS
                0070       CHARACTER*8 diagName
                0071       CHARACTER*4 diagSufx
                0072 #ifdef ALLOW_GENERIC_ADVDIFF
                0073       CHARACTER*4 GAD_DIAG_SUFX
                0074       EXTERNAL    GAD_DIAG_SUFX
                0075 #endif
                0076       LOGICAL     DIAGNOSTICS_IS_ON
                0077       EXTERNAL    DIAGNOSTICS_IS_ON
698b6992ee Jean*0078       _RL recip_dT
d64c4d306c Jean*0079       _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73b66b887d Jean*0080 #endif /* ALLOW_DIAGNOSTICS */
9366854e02 Chri*0081 CEOP
779cd6d73d Alis*0082
1189d93920 Patr*0083 cph(
                0084 cph Not good for TAF: may create irreducible control flow graph
                0085 cph      IF (Nr.LE.1) RETURN
                0086 cph)
1e1678964e Jean*0087
4e66ab0b67 Oliv*0088 #ifdef ALLOW_PTRACERS
                0089       IF ( tracerId.GE.GAD_TR1) THEN
                0090         DO k=1,Nr
                0091          deltaTX(k) = PTRACERS_dTLev(k)
                0092         ENDDO
                0093       ELSEIF ( tracerId.GE.1 ) THEN
                0094 #else
73b66b887d Jean*0095       IF ( tracerId.GE.1 ) THEN
4e66ab0b67 Oliv*0096 #endif
73b66b887d Jean*0097         DO k=1,Nr
                0098          deltaTX(k) = dTtracerLev(k)
                0099         ENDDO
                0100       ELSE
                0101         DO k=1,Nr
02d90fb24c Jean*0102          deltaTX(k) = deltaTMom
73b66b887d Jean*0103         ENDDO
                0104       ENDIF
                0105
8016bde8f3 Patr*0106 C--   Initialise
b8452ee69a Patr*0107       DO k=1,Nr
6f96b2277c Mart*0108        DO j=1-OLy,sNy+OLy
                0109         DO i=1-OLx,sNx+OLx
d64c4d306c Jean*0110          locTr(i,j,k) = 0. _d 0
b8452ee69a Patr*0111         ENDDO
                0112        ENDDO
                0113       ENDDO
8b6a578407 Patr*0114
8016bde8f3 Patr*0115 C--   Old aLower
46da898ec4 Jean*0116       DO j=1-OLy,sNy+OLy
                0117        DO i=1-OLx,sNx+OLx
                0118          a(i,j,1) = 0. _d 0
                0119        ENDDO
                0120       ENDDO
8016bde8f3 Patr*0121       DO k=2,Nr
6f96b2277c Mart*0122 #ifdef TARGET_NEC_SX
                0123        DO j=1-OLy,sNy+OLy
                0124         DO i=1-OLx,sNx+OLx
                0125 #else
bcd7bce512 Patr*0126        DO j=jMin,jMax
                0127         DO i=iMin,iMax
6f96b2277c Mart*0128 #endif
d64c4d306c Jean*0129           a(i,j,k) = -deltaTX(k)*recip_hFac(i,j,k)*recip_drF(k)
4606c28752 Jean*0130      &               *recip_deepFac2C(k)*recip_rhoFacC(k)
8016bde8f3 Patr*0131      &               *KappaRX(i,j, k )*recip_drC( k )
4606c28752 Jean*0132      &               *deepFac2F(k)*rhoFacF(k)
d64c4d306c Jean*0133           IF (recip_hFac(i,j,k-1).EQ.0.) a(i,j,k)=0.
8016bde8f3 Patr*0134         ENDDO
                0135        ENDDO
                0136       ENDDO
                0137
                0138 C--   Old aUpper
                0139       DO k=1,Nr-1
6f96b2277c Mart*0140 #ifdef TARGET_NEC_SX
                0141        DO j=1-OLy,sNy+OLy
                0142         DO i=1-OLx,sNx+OLx
                0143 #else
bcd7bce512 Patr*0144        DO j=jMin,jMax
                0145         DO i=iMin,iMax
6f96b2277c Mart*0146 #endif
d64c4d306c Jean*0147           c(i,j,k) = -deltaTX(k)*recip_hFac(i,j,k)*recip_drF(k)
4606c28752 Jean*0148      &               *recip_deepFac2C(k)*recip_rhoFacC(k)
8016bde8f3 Patr*0149      &               *KappaRX(i,j,k+1)*recip_drC(k+1)
4606c28752 Jean*0150      &               *deepFac2F(k+1)*rhoFacF(k+1)
d64c4d306c Jean*0151           IF (recip_hFac(i,j,k+1).EQ.0.) c(i,j,k)=0.
779cd6d73d Alis*0152         ENDDO
                0153        ENDDO
8016bde8f3 Patr*0154       ENDDO
46da898ec4 Jean*0155       DO j=1-OLy,sNy+OLy
                0156        DO i=1-OLx,sNx+OLx
                0157          c(i,j,Nr) = 0. _d 0
                0158        ENDDO
                0159       ENDDO
8b6a578407 Patr*0160
8016bde8f3 Patr*0161 C--   Old aCenter
                0162       DO k=1,Nr
46da898ec4 Jean*0163 #ifdef TARGET_NEC_SX
6f96b2277c Mart*0164        DO j=1-OLy,sNy+OLy
                0165         DO i=1-OLx,sNx+OLx
46da898ec4 Jean*0166 #else
                0167        DO j=jMin,jMax
                0168         DO i=iMin,iMax
                0169 #endif
93a010d2da Jean*0170           b(i,j,k) = 1. _d 0 - ( a(i,j,k) + c(i,j,k) )
                0171 C-    to recover older (prior to 2016-10-05) results:
                0172 c         b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
8016bde8f3 Patr*0173         ENDDO
                0174        ENDDO
                0175       ENDDO
                0176
                0177 C--   Old and new gam, bet are the same
                0178       DO k=1,Nr
6f96b2277c Mart*0179        DO j=1-OLy,sNy+OLy
                0180         DO i=1-OLx,sNx+OLx
74e25897f7 Jean*0181           bet(i,j,k) = 1. _d 0
8016bde8f3 Patr*0182           gam(i,j,k) = 0. _d 0
                0183         ENDDO
                0184        ENDDO
                0185       ENDDO
8adf9f02ba Patr*0186
8016bde8f3 Patr*0187 C--   Only need do anything if Nr>1
46da898ec4 Jean*0188       IF (Nr.GT.1) THEN
8adf9f02ba Patr*0189
8016bde8f3 Patr*0190        k = 1
                0191 C--    Beginning of forward sweep (top level)
46da898ec4 Jean*0192 #ifdef TARGET_NEC_SX
6f96b2277c Mart*0193        DO j=1-OLy,sNy+OLy
                0194         DO i=1-OLx,sNx+OLx
46da898ec4 Jean*0195 #else
                0196        DO j=jMin,jMax
                0197         DO i=iMin,iMax
                0198 #endif
8016bde8f3 Patr*0199          IF (b(i,j,1).NE.0.) bet(i,j,1) = 1. _d 0 / b(i,j,1)
8adf9f02ba Patr*0200         ENDDO
                0201        ENDDO
                0202
46da898ec4 Jean*0203       ENDIF
                0204
a0b25fcf44 Chri*0205 C--   Middle of forward sweep
46da898ec4 Jean*0206       IF (Nr.GE.2) THEN
8b6a578407 Patr*0207
8016bde8f3 Patr*0208 CADJ loop = sequential
                0209        DO k=2,Nr
8adf9f02ba Patr*0210
46da898ec4 Jean*0211 #ifdef TARGET_NEC_SX
6f96b2277c Mart*0212         DO j=1-OLy,sNy+OLy
                0213          DO i=1-OLx,sNx+OLx
46da898ec4 Jean*0214 #else
                0215         DO j=jMin,jMax
                0216          DO i=iMin,iMax
                0217 #endif
8016bde8f3 Patr*0218           gam(i,j,k) = c(i,j,k-1)*bet(i,j,k-1)
d64c4d306c Jean*0219           IF ( ( b(i,j,k) - a(i,j,k)*gam(i,j,k) ) .NE. 0.)
8016bde8f3 Patr*0220      &        bet(i,j,k) = 1. _d 0 / ( b(i,j,k) - a(i,j,k)*gam(i,j,k) )
779cd6d73d Alis*0221          ENDDO
                0222         ENDDO
8adf9f02ba Patr*0223
779cd6d73d Alis*0224        ENDDO
8b6a578407 Patr*0225
46da898ec4 Jean*0226       ENDIF
8b6a578407 Patr*0227
6f96b2277c Mart*0228 #ifdef TARGET_NEC_SX
                0229       DO j=1-OLy,sNy+OLy
                0230        DO i=1-OLx,sNx+OLx
                0231 #else
8016bde8f3 Patr*0232       DO j=jMin,jMax
                0233        DO i=iMin,iMax
6f96b2277c Mart*0234 #endif
23a7f3050f Jean*0235         locTr(i,j,1) = gTracer(i,j,1)*bet(i,j,1)
779cd6d73d Alis*0236        ENDDO
8016bde8f3 Patr*0237       ENDDO
                0238       DO k=2,Nr
6f96b2277c Mart*0239 #ifdef TARGET_NEC_SX
                0240        DO j=1-OLy,sNy+OLy
                0241         DO i=1-OLx,sNx+OLx
                0242 #else
779cd6d73d Alis*0243        DO j=jMin,jMax
                0244         DO i=iMin,iMax
6f96b2277c Mart*0245 #endif
d64c4d306c Jean*0246          locTr(i,j,k) = bet(i,j,k)*
23a7f3050f Jean*0247      &        (gTracer(i,j,k) - a(i,j,k)*locTr(i,j,k-1))
8adf9f02ba Patr*0248         ENDDO
                0249        ENDDO
8016bde8f3 Patr*0250       ENDDO
8adf9f02ba Patr*0251
a0b25fcf44 Chri*0252 C--    Backward sweep
8016bde8f3 Patr*0253 CADJ loop = sequential
46da898ec4 Jean*0254        DO k=Nr-1,1,-1
6f96b2277c Mart*0255 #ifdef TARGET_NEC_SX
698b6992ee Jean*0256         DO j=1-OLy,sNy+OLy
                0257          DO i=1-OLx,sNx+OLx
6f96b2277c Mart*0258 #else
46da898ec4 Jean*0259         DO j=jMin,jMax
                0260          DO i=iMin,iMax
6f96b2277c Mart*0261 #endif
d64c4d306c Jean*0262           locTr(i,j,k) = locTr(i,j,k) - gam(i,j,k+1)*locTr(i,j,k+1)
46da898ec4 Jean*0263          ENDDO
                0264         ENDDO
                0265        ENDDO
                0266
                0267        DO k=1,Nr
                0268 #ifdef TARGET_NEC_SX
                0269         DO j=1-OLy,sNy+OLy
                0270          DO i=1-OLx,sNx+OLx
                0271 #else
                0272         DO j=jMin,jMax
                0273          DO i=iMin,iMax
                0274 #endif
698b6992ee Jean*0275           locUpdate =  locTr(i,j,k) - gTracer(i,j,k)
46da898ec4 Jean*0276           gTracer(i,j,k) = locTr(i,j,k)
698b6992ee Jean*0277           locTr(i,j,k) = locUpdate
46da898ec4 Jean*0278          ENDDO
8016bde8f3 Patr*0279         ENDDO
                0280        ENDDO
779cd6d73d Alis*0281
73b66b887d Jean*0282 #ifdef ALLOW_DIAGNOSTICS
698b6992ee Jean*0283 C--   Diagnostics of momentum dissipation/viscous tendency, implicit part:
                0284       IF ( useDiagnostics .AND.
                0285      &     ( tracerId.EQ. -1 .OR. tracerId.EQ. -2 ) ) THEN
                0286         IF ( tracerId.EQ. -1 ) diagName = 'Um_ImplD'
                0287         IF ( tracerId.EQ. -2 ) diagName = 'Vm_ImplD'
                0288         IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN
                0289           recip_dT = 0. _d 0
                0290           IF ( deltaTMom.GT.zeroRL ) recip_dT = oneRL / deltaTMom
                0291           DO k=1,Nr
                0292            DO j=1-OLy,sNy+OLy
                0293             DO i=1-OLx,sNx+OLx
                0294               locTr(i,j,k) = locTr(i,j,k)*recip_dT
                0295             ENDDO
                0296            ENDDO
                0297           ENDDO
                0298           CALL DIAGNOSTICS_FILL( locTr,diagName, 0,Nr,2,bi,bj, myThid )
                0299         ENDIF
                0300       ENDIF
                0301
                0302 C--   Diagnostics of vertical diffusion flux (or vertical viscous flux):
2e7d726b27 Jean*0303       IF ( useDiagnostics .AND.tracerId.NE.0 ) THEN
                0304         IF ( tracerId.GE. 1 ) THEN
73b66b887d Jean*0305 C--   Set diagnostic suffix for the current tracer
                0306 #ifdef ALLOW_GENERIC_ADVDIFF
2e7d726b27 Jean*0307           diagSufx = GAD_DIAG_SUFX( tracerId, myThid )
73b66b887d Jean*0308 #else
2e7d726b27 Jean*0309           diagSufx = 'aaaa'
73b66b887d Jean*0310 #endif
2e7d726b27 Jean*0311           diagName = 'DFrI'//diagSufx
                0312         ELSEIF ( tracerId.EQ. -1 ) THEN
                0313           diagName = 'VISrI_Um'
                0314         ELSEIF ( tracerId.EQ. -2 ) THEN
                0315           diagName = 'VISrI_Vm'
                0316         ELSE
                0317           STOP 'IMPLIDIFF: should never reach this point !'
                0318         ENDIF
9b89fcf692 antn*0319         IF ( DIAGNOSTICS_IS_ON(diagName,myThid)
                0320 #ifdef ALLOW_LAYERS
                0321      &       .OR. layers_useThermo
                0322 #endif
                0323      &       ) THEN
73b66b887d Jean*0324          DO k= 1,Nr
                0325           IF ( k.EQ.1 ) THEN
                0326 C-  Note: Needs to call DIAGNOSTICS_FILL at level k=1 even if array == 0
                0327 C         otherwise counter is not incremented !!
                0328             DO j=1-OLy,sNy+OLy
                0329              DO i=1-OLx,sNx+OLx
                0330                df(i,j) = 0. _d 0
                0331              ENDDO
                0332             ENDDO
2e7d726b27 Jean*0333           ELSEIF ( tracerId.GE.1 ) THEN
6f96b2277c Mart*0334 #ifdef TARGET_NEC_SX
                0335             DO j=1-OLy,sNy+OLy
                0336              DO i=1-OLx,sNx+OLx
                0337 #else
73b66b887d Jean*0338             DO j=1,sNy
                0339              DO i=1,sNx
6f96b2277c Mart*0340 #endif
73b66b887d Jean*0341                df(i,j) =
4606c28752 Jean*0342      &             -rA(i,j,bi,bj)*deepFac2F(k)*rhoFacF(k)
d64c4d306c Jean*0343      &            * KappaRX(i,j,k)*recip_drC(k)*rkSign
46da898ec4 Jean*0344      &            * (gTracer(i,j,k) - gTracer(i,j,k-1))
d64c4d306c Jean*0345      &            * maskC(i,j,k,bi,bj)
                0346      &            * maskC(i,j,k-1,bi,bj)
2e7d726b27 Jean*0347              ENDDO
                0348             ENDDO
                0349           ELSEIF ( tracerId.EQ.-1 ) THEN
6f96b2277c Mart*0350 #ifdef TARGET_NEC_SX
                0351             DO j=1-OLy,sNy+OLy
                0352              DO i=1-OLx,sNx+OLx
                0353 #else
2e7d726b27 Jean*0354             DO j=1,sNy
                0355              DO i=1,sNx+1
6f96b2277c Mart*0356 #endif
2e7d726b27 Jean*0357                df(i,j) =
4606c28752 Jean*0358      &             -rAw(i,j,bi,bj)*deepFac2F(k)*rhoFacF(k)
d64c4d306c Jean*0359      &            * KappaRX(i,j,k)*recip_drC(k)*rkSign
46da898ec4 Jean*0360      &            * (gTracer(i,j,k) - gTracer(i,j,k-1))
2e7d726b27 Jean*0361      &            * _maskW(i,j,k,bi,bj)
                0362      &            * _maskW(i,j,k-1,bi,bj)
                0363              ENDDO
                0364             ENDDO
                0365           ELSEIF ( tracerId.EQ.-2 ) THEN
6f96b2277c Mart*0366 #ifdef TARGET_NEC_SX
                0367             DO j=1-OLy,sNy+OLy
                0368              DO i=1-OLx,sNx+OLx
                0369 #else
2e7d726b27 Jean*0370             DO j=1,sNy+1
                0371              DO i=1,sNx
6f96b2277c Mart*0372 #endif
2e7d726b27 Jean*0373                df(i,j) =
4606c28752 Jean*0374      &             -rAs(i,j,bi,bj)*deepFac2F(k)*rhoFacF(k)
d64c4d306c Jean*0375      &            * KappaRX(i,j,k)*recip_drC(k)*rkSign
46da898ec4 Jean*0376      &            * (gTracer(i,j,k) - gTracer(i,j,k-1))
2e7d726b27 Jean*0377      &            * _maskS(i,j,k,bi,bj)
                0378      &            * _maskS(i,j,k-1,bi,bj)
73b66b887d Jean*0379              ENDDO
                0380             ENDDO
                0381           ENDIF
2ceedfef66 Jean*0382           CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
cf336ab6c5 Ryan*0383 #ifdef ALLOW_LAYERS
9b89fcf692 antn*0384           IF ( layers_useThermo ) THEN
50d8304171 Ryan*0385            CALL LAYERS_FILL( df, tracerId, 'DFR',
cf336ab6c5 Ryan*0386      &                           k, 1, 2,bi,bj, myThid )
ee16a2cae4 Ryan*0387           ENDIF
cf336ab6c5 Ryan*0388 #endif /* ALLOW_LAYERS */
73b66b887d Jean*0389          ENDDO
                0390         ENDIF
                0391       ENDIF
                0392 #endif /* ALLOW_DIAGNOSTICS */
                0393
779cd6d73d Alis*0394       RETURN
                0395       END