Back to home page

MITgcm
 
 

    

File indexing completed on 2021-06-17 05:11:07 UTC

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