/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F
ViewVC logotype

Diff of /MITgcm/pkg/generic_advdiff/gad_calc_rhs.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.44 by adcroft, Sat Jan 20 21:20:11 2007 UTC revision 1.58 by jmc, Thu Dec 1 14:16:30 2011 UTC
# Line 12  C !INTERFACE: ========================== Line 12  C !INTERFACE: ==========================
12       I           xA, yA, maskUp, uFld, vFld, wFld,       I           xA, yA, maskUp, uFld, vFld, wFld,
13       I           uTrans, vTrans, rTrans, rTransKp1,       I           uTrans, vTrans, rTrans, rTransKp1,
14       I           diffKh, diffK4, KappaR, TracerN, TracAB,       I           diffKh, diffK4, KappaR, TracerN, TracAB,
15       I           tracerIdentity, advectionScheme, vertAdvecScheme,       I           deltaTLev, tracerIdentity,
16         I           advectionScheme, vertAdvecScheme,
17       I           calcAdvection, implicitAdvection, applyAB_onTracer,       I           calcAdvection, implicitAdvection, applyAB_onTracer,
18         I           trUseGMRedi, trUseKPP,
19       U           fVerT, gTracer,       U           fVerT, gTracer,
20       I           myTime, myIter, myThid )       I           myTime, myIter, myThid )
21    
# Line 77  C vertAdvecScheme  :: advection scheme t Line 79  C vertAdvecScheme  :: advection scheme t
79  C calcAdvection    :: =False if Advec computed with multiDim scheme  C calcAdvection    :: =False if Advec computed with multiDim scheme
80  C implicitAdvection:: =True if vertical Advec computed implicitly  C implicitAdvection:: =True if vertical Advec computed implicitly
81  C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)  C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)
82    C trUseGMRedi      :: true if this tracer uses GM-Redi
83    C trUseKPP         :: true if this tracer uses KPP
84  C myTime           :: current time  C myTime           :: current time
85  C myIter           :: iteration number  C myIter           :: iteration number
86  C myThid           :: thread number  C myThid           :: thread number
# Line 96  C myThid           :: thread number Line 100  C myThid           :: thread number
100        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
102        _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
103          _RL deltaTLev(Nr)
104        INTEGER tracerIdentity        INTEGER tracerIdentity
105        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
106        LOGICAL calcAdvection        LOGICAL calcAdvection
107        LOGICAL implicitAdvection, applyAB_onTracer        LOGICAL implicitAdvection, applyAB_onTracer
108          LOGICAL trUseGMRedi, trUseKPP
109        _RL     myTime        _RL     myTime
110        INTEGER myIter, myThid        INTEGER myIter, myThid
111    
# Line 124  C locABT           :: local copy of (AB- Line 130  C locABT           :: local copy of (AB-
130        EXTERNAL    GAD_DIAG_SUFX        EXTERNAL    GAD_DIAG_SUFX
131  #endif  #endif
132        INTEGER i,j        INTEGER i,j
133          _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134          _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 132  C locABT           :: local copy of (AB- Line 140  C locABT           :: local copy of (AB-
140        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141        _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
142        _RL advFac, rAdvFac        _RL advFac, rAdvFac
143    #ifdef GAD_SMOLARKIEWICZ_HACK
144          _RL outFlux, trac, fac, gTrFac
145    #endif
146  CEOP  CEOP
147    
148  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 150  C--   Set diagnostic suffix for the curr Line 161  C--   Set diagnostic suffix for the curr
161        advFac  = 0. _d 0        advFac  = 0. _d 0
162        IF (calcAdvection) advFac = 1. _d 0        IF (calcAdvection) advFac = 1. _d 0
163        rAdvFac = rkSign*advFac        rAdvFac = rkSign*advFac
164        IF (implicitAdvection) rAdvFac = 0. _d 0        IF (implicitAdvection) rAdvFac = rkSign
165    
166        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
167         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
# Line 207  C--   Initialize net flux in X direction Line 218  C--   Initialize net flux in X direction
218  C-    Advective flux in X  C-    Advective flux in X
219        IF (calcAdvection) THEN        IF (calcAdvection) THEN
220          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
221            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)             CALL GAD_C2_ADV_X( bi,bj,k, uTrans, locABT, af, myThid )
222          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
223       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
224            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
225       I            dTtracerLev(k), uTrans, uFld, locABT,       I              deltaTLev(k), uTrans, uFld, locABT,
226       O            af, myThid )       O              af, myThid )
227          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSE
228            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),           DO j=1-OLy,sNy+OLy
229       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,            DO i=1-OLx,sNx+OLx
230       O            af, myThid )  #ifdef ALLOW_OBCS
231          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN             maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
232            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)  #else /* ALLOW_OBCS */
233          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN             maskLocW(i,j) = _maskW(i,j,k,bi,bj)
234            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)  #endif /* ALLOW_OBCS */
235          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            ENDDO
236            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),           ENDDO
237       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
238       O            af, myThid )             CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
239          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       I              uTrans, uFld, maskLocW, locABT,
240           IF ( inAdMode ) THEN       O              af, myThid )
241             ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
242               CALL GAD_U3_ADV_X( bi,bj,k, uTrans, maskLocW, locABT,
243         O              af, myThid )
244             ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
245               CALL GAD_C4_ADV_X( bi,bj,k, uTrans, maskLocW, locABT,
246         O              af, myThid )
247             ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
248               CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
249         I              uTrans, uFld, maskLocW, locABT,
250         O              af, myThid )
251             ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
252              IF ( inAdMode ) THEN
253  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
254  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
255  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
256            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),             CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
257       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I              uTrans, uFld, maskLocW, locABT,
258       O           af, myThid )       O              af, myThid )
259              ELSE
260               CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
261         I              uTrans, uFld, maskLocW, locABT,
262         O              af, myThid )
263              ENDIF
264    #ifndef ALLOW_AUTODIFF_TAMC
265             ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
266               CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
267         I              uTrans, uFld, maskLocW, locABT,
268         O              af, myThid )
269    #endif
270           ELSE           ELSE
271            CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),            STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
      I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,  
      O           af, myThid )  
272           ENDIF           ENDIF
         ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN  
           CALL GAD_OS7MP_ADV_X( bi,bj,k, dTtracerLev(k),  
      I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,  
      O            af, myThid )  
         ELSE  
          STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'  
273          ENDIF          ENDIF
274          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
275           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
# Line 276  C-    Add bi-harmonic diffusive flux in Line 302  C-    Add bi-harmonic diffusive flux in
302    
303  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
304  C-    GM/Redi flux in X  C-    GM/Redi flux in X
305        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
306  C *note* should update GMREDI_XTRANSPORT to set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
307          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
308            CALL GMREDI_XTRANSPORT(            CALL GMREDI_XTRANSPORT(
# Line 304  C     anelastic: advect.fluxes are scale Line 330  C     anelastic: advect.fluxes are scale
330  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
331  C       excluding advective terms:  C       excluding advective terms:
332        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
333       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
334            diagName = 'DFxE'//diagSufx            diagName = 'DFxE'//diagSufx
335            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
336        ENDIF        ENDIF
# Line 320  C--   Initialize net flux in Y direction Line 346  C--   Initialize net flux in Y direction
346  C-    Advective flux in Y  C-    Advective flux in Y
347        IF (calcAdvection) THEN        IF (calcAdvection) THEN
348          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
349            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)             CALL GAD_C2_ADV_Y( bi,bj,k, vTrans, locABT, af, myThid )
350          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
351       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
352            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
353       I            dTtracerLev(k), vTrans, vFld, locABT,       I              deltaTLev(k), vTrans, vFld, locABT,
354       O            af, myThid )       O              af, myThid )
355          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSE
356            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),           DO j=1-OLy,sNy+OLy
357       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,            DO i=1-OLx,sNx+OLx
358       O            af, myThid )  #ifdef ALLOW_OBCS
359          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN             maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
360            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)  #else /* ALLOW_OBCS */
361          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN             maskLocS(i,j) = _maskS(i,j,k,bi,bj)
362            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)  #endif /* ALLOW_OBCS */
363          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            ENDDO
364            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),           ENDDO
365       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
366       O            af, myThid )             CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
367          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       I              vTrans, vFld, maskLocS, locABT,
368           IF ( inAdMode ) THEN       O              af, myThid )
369             ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
370               CALL GAD_U3_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
371         O              af, myThid )
372             ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
373               CALL GAD_C4_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
374         O              af, myThid )
375             ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
376               CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
377         I              vTrans, vFld, maskLocS, locABT,
378         O              af, myThid )
379             ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
380              IF ( inAdMode ) THEN
381  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
382  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
383  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
384            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),             CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
385       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I              vTrans, vFld, maskLocS, locABT,
386       O           af, myThid )       O              af, myThid )
387              ELSE
388               CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
389         I              vTrans, vFld, maskLocS, locABT,
390         O              af, myThid )
391              ENDIF
392    #ifndef ALLOW_AUTODIFF_TAMC
393             ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
394               CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
395         I              vTrans, vFld, maskLocS, locABT,
396         O              af, myThid )
397    #endif
398           ELSE           ELSE
399            CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),             STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
      I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,  
      O           af, myThid )  
400           ENDIF           ENDIF
         ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN  
           CALL GAD_OS7MP_ADV_Y( bi,bj,k, dTtracerLev(k),  
      I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,  
      O            af, myThid )  
         ELSE  
           STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'  
401          ENDIF          ENDIF
402          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
403           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
# Line 389  C-    Add bi-harmonic flux in Y Line 430  C-    Add bi-harmonic flux in Y
430    
431  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
432  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
433        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
434  C *note* should update GMREDI_YTRANSPORT to set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
435          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
436            CALL GMREDI_YTRANSPORT(            CALL GMREDI_YTRANSPORT(
# Line 417  C     anelastic: advect.fluxes are scale Line 458  C     anelastic: advect.fluxes are scale
458  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
459  C       excluding advective terms:  C       excluding advective terms:
460        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
461       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
462            diagName = 'DFyE'//diagSufx            diagName = 'DFyE'//diagSufx
463            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
464        ENDIF        ENDIF
# Line 439  C-    Compute vertical advective flux in Line 480  C-    Compute vertical advective flux in
480          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
481       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
482            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
483       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
484       O         af, myThid )       O         af, myThid )
485          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
486            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
487       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
488       O         af, myThid )       O         af, myThid )
489          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
490            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
# Line 451  C-    Compute vertical advective flux in Line 492  C-    Compute vertical advective flux in
492            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
493          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
494            CALL GAD_DST3_ADV_R( bi,bj,k,            CALL GAD_DST3_ADV_R( bi,bj,k,
495       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
496       O         af, myThid )       O         af, myThid )
497          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
498  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
# Line 459  cph IF inAdExact=.FALSE., we want to use Line 500  cph IF inAdExact=.FALSE., we want to use
500  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
501            IF ( inAdMode ) THEN            IF ( inAdMode ) THEN
502             CALL GAD_DST3_ADV_R( bi,bj,k,             CALL GAD_DST3_ADV_R( bi,bj,k,
503       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
504       O         af, myThid )       O         af, myThid )
505            ELSE            ELSE
506             CALL GAD_DST3FL_ADV_R( bi,bj,k,             CALL GAD_DST3FL_ADV_R( bi,bj,k,
507       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
508       O         af, myThid )       O         af, myThid )
509            ENDIF            ENDIF
510    #ifndef ALLOW_AUTODIFF_TAMC
511          ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
512             CALL GAD_OS7MP_ADV_R( bi,bj,k,             CALL GAD_OS7MP_ADV_R( bi,bj,k,
513       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
514       O         af, myThid )       O         af, myThid )
515    #endif
516          ELSE          ELSE
517            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
518          ENDIF          ENDIF
# Line 509  C           boundary condition. Line 552  C           boundary condition.
552    
553  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
554  C-    GM/Redi flux in R  C-    GM/Redi flux in R
555        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
556  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
557          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
558            CALL GMREDI_RTRANSPORT(            CALL GMREDI_RTRANSPORT(
# Line 537  C *note* should update GMREDI_RTRANSPORT Line 580  C *note* should update GMREDI_RTRANSPORT
580  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
581  C       Explicit terms only & excluding advective terms:  C       Explicit terms only & excluding advective terms:
582        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
583       &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN       &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN
584            diagName = 'DFrE'//diagSufx            diagName = 'DFrE'//diagSufx
585            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
586        ENDIF        ENDIF
# Line 545  C       Explicit terms only & excluding Line 588  C       Explicit terms only & excluding
588    
589  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
590  C-    Set non local KPP transport term (ghat):  C-    Set non local KPP transport term (ghat):
591        IF ( useKPP .AND. k.GE.2 ) THEN        IF ( trUseKPP .AND. k.GE.2 ) THEN
592         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
593          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
594           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
# Line 553  C-    Set non local KPP transport term ( Line 596  C-    Set non local KPP transport term (
596         ENDDO         ENDDO
597         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
598          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
599       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
600       O     df )       O           df,
601         I           myTime, myIter, myThid )
602         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
603          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
604       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
605       O     df )       O           df,
606         I           myTime, myIter, myThid )
607  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
608         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
609          CALL KPP_TRANSPORT_PTR(          CALL KPP_TRANSPORT_PTR(
610       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
611       I     tracerIdentity-GAD_TR1+1,       I           tracerIdentity-GAD_TR1+1,
612       O     df )       O           df,
613         I           myTime, myIter, myThid )
614  #endif  #endif
615         ELSE         ELSE
616          PRINT*,'invalid tracer indentity: ', tracerIdentity          WRITE(errorMessageUnit,*)
617          STOP 'GAD_CALC_RHS: Ooops'       &    'tracer identity =', tracerIdentity, ' is not valid => STOP'
618            STOP 'ABNORMAL END: S/R GAD_CALC_RHS: invalid tracer identity'
619         ENDIF         ENDIF
620         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
621          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
# Line 576  C-    Set non local KPP transport term ( Line 623  C-    Set non local KPP transport term (
623       &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)       &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
624          ENDDO          ENDDO
625         ENDDO         ENDDO
626    #ifdef ALLOW_DIAGNOSTICS
627    C-    Diagnostics of Non-Local Tracer (vertical) flux
628           IF ( useDiagnostics ) THEN
629             diagName = 'KPPg'//diagSufx
630             CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
631    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
632    C        does it only if k=1 (never the case here)
633             IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
634           ENDIF
635    #endif
636          ENDIF
637    #endif /* ALLOW_KPP */
638    
639    #ifdef GAD_SMOLARKIEWICZ_HACK
640    coj   Hack to make redi (and everything else in this s/r) positive
641    coj   (see Smolarkiewicz MWR 1989 and Bott MWR 1989).
642    coj   Only works if 'down' is k+1 and k loop in thermodynamics is k=Nr,1,-1
643    coj
644    coj   Apply to all tracers except temperature
645          IF (tracerIdentity.NE.GAD_TEMPERATURE .AND.
646         &    tracerIdentity.NE.GAD_SALINITY) THEN
647           DO j=1-Oly,sNy+Oly-1
648            DO i=1-Olx,sNx+Olx-1
649    coj   Add outgoing fluxes
650             outFlux=deltaTLev(k)*
651         &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
652         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
653         &    *( MAX(0. _d 0,fZon(i+1,j)) + MAX(0. _d 0,-fZon(i,j))
654         &      +MAX(0. _d 0,fMer(i,j+1)) + MAX(0. _d 0,-fMer(i,j))
655         &      +MAX(0. _d 0,fVerT(i,j,kDown)*rkSign)
656         &      +MAX(0. _d 0,-fVerT(i,j,kUp)*rkSign)
657         &     )
658             IF ( applyAB_onTracer ) THEN
659               trac=TracerN(i,j,k,bi,bj)
660             ELSE
661               trac=TracAB(i,j,k,bi,bj)
662             ENDIF
663    coj   If they would reduce tracer by a fraction of more than
664    coj   SmolarkiewiczMaxFrac, scale them down
665             IF (outFlux.GT.0. _d 0 .AND.
666         &       outFlux.GT.SmolarkiewiczMaxFrac*trac) THEN
667    coj   If tracer is already negative, scale flux to zero
668               fac = MAX(0. _d 0,SmolarkiewiczMaxFrac*trac/outFlux)
669    
670               IF (fZon(i+1,j).GT.0. _d 0) fZon(i+1,j)=fac*fZon(i+1,j)
671               IF (-fZon(i,j) .GT.0. _d 0) fZon(i,j)  =fac*fZon(i,j)
672               IF (fMer(i,j+1).GT.0. _d 0) fMer(i,j+1)=fac*fMer(i,j+1)
673               IF (-fMer(i,j) .GT.0. _d 0) fMer(i,j)  =fac*fMer(i,j)
674               IF (-fVerT(i,j,kUp)*rkSign .GT.0. _d 0)
675         &       fVerT(i,j,kUp)=fac*fVerT(i,j,kUp)
676    
677               IF (k.LT.Nr .AND. fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN
678    coj   Down flux is special: it has already been applied in lower layer,
679    coj   so we have to readjust this.
680    coj   Note: for k+1, gTracer is now the updated tracer, not the tendency!
681    coj   thus it has an extra factor deltaTLev(k+1)
682                 gTrFac=deltaTLev(k+1)
683    coj   Other factors that have been applied to gTracer since the last call:
684    #ifdef NONLIN_FRSURF
685                 IF (nonlinFreeSurf.GT.0) THEN
686                  IF (select_rStar.GT.0) THEN
687    #ifndef DISABLE_RSTAR_CODE
688                    gTrFac = gTrFac/rStarExpC(i,j,bi,bj)
689    #endif /* DISABLE_RSTAR_CODE */
690                  ENDIF
691                 ENDIF
692    #endif /* NONLIN_FRSURF */
693    coj   Now: undo down flux, ...
694                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
695         &        +gTrFac
696         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
697         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
698         &         *recip_rhoFacC(k+1)
699         &         *( -fVerT(i,j,kDown)*rkSign )
700    coj   ... scale ...
701                 fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)
702    coj   ... and reapply
703                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
704         &        +gTrFac
705         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
706         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
707         &         *recip_rhoFacC(k+1)
708         &         *( fVerT(i,j,kDown)*rkSign )
709               ENDIF
710    
711             ENDIF
712            ENDDO
713           ENDDO
714        ENDIF        ENDIF
715  #endif  #endif
716    
# Line 589  C     Anelastic: scale vertical fluxes b Line 724  C     Anelastic: scale vertical fluxes b
724       &   *( (fZon(i+1,j)-fZon(i,j))       &   *( (fZon(i+1,j)-fZon(i,j))
725       &     +(fMer(i,j+1)-fMer(i,j))       &     +(fMer(i,j+1)-fMer(i,j))
726       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
727       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
728       &                   +(vTrans(i,j+1)-vTrans(i,j))       &                   +(vTrans(i,j+1)-vTrans(i,j))*advFac
729       &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac       &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
730       &                  )*advFac       &                  )
731       &    )       &    )
732         ENDDO         ENDDO
733        ENDDO        ENDDO
734    
735  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
736        IF ( debugLevel .GE. debLevB        IF ( debugLevel .GE. debLevC
737       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
738       &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0       &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
739       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
Legend:
Removed from v.1.44  
changed lines
  Added in v.1.58
  ViewVC Help
Powered by ViewVC 1.1.22