/[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.63 by jmc, Fri Dec 27 15:52:17 2013 UTC
# Line 11  C !INTERFACE: ========================== Line 11  C !INTERFACE: ==========================
11       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
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, diffKr4, TracerN, TracAB,
15       I           tracerIdentity, advectionScheme, vertAdvecScheme,       I           deltaTLev, trIdentity,
16         I           advectionScheme, vertAdvecScheme,
17       I           calcAdvection, implicitAdvection, applyAB_onTracer,       I           calcAdvection, implicitAdvection, applyAB_onTracer,
18         I           trUseDiffKr4, trUseGMRedi, trUseKPP,
19         O           fZon, fMer,
20       U           fVerT, gTracer,       U           fVerT, gTracer,
21       I           myTime, myIter, myThid )       I           myTime, myIter, myThid )
22    
# Line 45  C !USES: =============================== Line 48  C !USES: ===============================
48  #include "SURFACE.h"  #include "SURFACE.h"
49  #include "GAD.h"  #include "GAD.h"
50    
51    #ifdef ALLOW_AUTODIFF
52    # include "AUTODIFF_PARAMS.h"
53    #endif /* ALLOW_AUTODIFF */
54  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
55  #include "tamc.h"  # include "tamc.h"
56  #include "tamc_keys.h"  # include "tamc_keys.h"
57  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
58    
59  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
# Line 65  C uTrans,vTrans    :: 2-D arrays of volu Line 71  C uTrans,vTrans    :: 2-D arrays of volu
71  C rTrans           :: 2-D arrays of volume transports at W points  C rTrans           :: 2-D arrays of volume transports at W points
72  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
73  C diffKh           :: horizontal diffusion coefficient  C diffKh           :: horizontal diffusion coefficient
74  C diffK4           :: bi-harmonic diffusion coefficient  C diffK4           :: horizontal bi-harmonic diffusion coefficient
75  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
76    C diffKr4          :: 1-D array for vertical bi-harmonic diffusion coefficient
77  C TracerN          :: tracer field @ time-step n (Note: only used  C TracerN          :: tracer field @ time-step n (Note: only used
78  C                     if applying AB on tracer field rather than on tendency gTr)  C                     if applying AB on tracer field rather than on tendency gTr)
79  C TracAB           :: current tracer field (@ time-step n if applying AB on gTr  C TracAB           :: current tracer field (@ time-step n if applying AB on gTr
80  C                     or extrapolated fwd in time to n+1/2 if applying AB on Tr)  C                     or extrapolated fwd in time to n+1/2 if applying AB on Tr)
81  C tracerIdentity   :: tracer identifier (required for KPP,GM)  C trIdentity       :: tracer identifier (required for KPP,GM)
82  C advectionScheme  :: advection scheme to use (Horizontal plane)  C advectionScheme  :: advection scheme to use (Horizontal plane)
83  C vertAdvecScheme  :: advection scheme to use (Vertical direction)  C vertAdvecScheme  :: advection scheme to use (Vertical direction)
84  C calcAdvection    :: =False if Advec computed with multiDim scheme  C calcAdvection    :: =False if Advec computed with multiDim scheme
85  C implicitAdvection:: =True if vertical Advec computed implicitly  C implicitAdvection:: =True if vertical Advec computed implicitly
86  C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)  C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)
87    C trUseDiffKr4     :: true if this tracer uses vertical bi-harmonic diffusion
88    C trUseGMRedi      :: true if this tracer uses GM-Redi
89    C trUseKPP         :: true if this tracer uses KPP
90  C myTime           :: current time  C myTime           :: current time
91  C myIter           :: iteration number  C myIter           :: iteration number
92  C myThid           :: thread number  C myThid           :: thread number
# Line 94  C myThid           :: thread number Line 104  C myThid           :: thread number
104        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL diffKh, diffK4        _RL diffKh, diffK4
106        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107          _RL diffKr4(Nr)
108        _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)
109        _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)
110        INTEGER tracerIdentity        _RL deltaTLev(Nr)
111          INTEGER trIdentity
112        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
113        LOGICAL calcAdvection        LOGICAL calcAdvection
114        LOGICAL implicitAdvection, applyAB_onTracer        LOGICAL implicitAdvection, applyAB_onTracer
115          LOGICAL trUseDiffKr4, trUseGMRedi, trUseKPP
116        _RL     myTime        _RL     myTime
117        INTEGER myIter, myThid        INTEGER myIter, myThid
118    
119  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
120  C gTracer          :: tendency array  C gTracer          :: tendency array
121    C fZon             :: zonal flux
122    C fMer             :: meridional flux
123  C fVerT            :: 2 1/2D arrays for vertical advective flux  C fVerT            :: 2 1/2D arrays for vertical advective flux
124        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
125          _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126          _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
128    
129    C !FUNCTIONS:       ====================================================
130    #ifdef ALLOW_DIAGNOSTICS
131          CHARACTER*4 GAD_DIAG_SUFX
132          EXTERNAL    GAD_DIAG_SUFX
133    #endif /* ALLOW_DIAGNOSTICS */
134    
135  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
136  C i,j              :: loop indices  C i,j              :: loop indices
137  C df4              :: used for storing del^2 T for bi-harmonic term  C df4              :: used for storing del^2 T for bi-harmonic term
 C fZon             :: zonal flux  
 C fMer             :: meridional flux  
138  C af               :: advective flux  C af               :: advective flux
139  C df               :: diffusive flux  C df               :: diffusive flux
140  C localT           :: local copy of tracer field  C localT           :: local copy of tracer field
141  C locABT           :: local copy of (AB-extrapolated) tracer field  C locABT           :: local copy of (AB-extrapolated) tracer field
 #ifdef ALLOW_DIAGNOSTICS  
       CHARACTER*8 diagName  
       CHARACTER*4 GAD_DIAG_SUFX, diagSufx  
       EXTERNAL    GAD_DIAG_SUFX  
 #endif  
142        INTEGER i,j        INTEGER i,j
143          _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144          _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
146        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149        _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150        _RL advFac, rAdvFac        _RL advFac, rAdvFac
151    #ifdef GAD_SMOLARKIEWICZ_HACK
152          _RL outFlux, trac, fac, gTrFac
153    #endif
154    #ifdef ALLOW_DIAGNOSTICS
155          CHARACTER*8 diagName
156          CHARACTER*4 diagSufx
157    #endif
158  CEOP  CEOP
159    
160  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 143  C--   the kDown is still required Line 166  C--   the kDown is still required
166  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
167  C--   Set diagnostic suffix for the current tracer  C--   Set diagnostic suffix for the current tracer
168        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
169          diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )          diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
170        ENDIF        ENDIF
171  #endif  #endif
172    
173        advFac  = 0. _d 0        advFac  = 0. _d 0
174        IF (calcAdvection) advFac = 1. _d 0        IF (calcAdvection) advFac = 1. _d 0
175        rAdvFac = rkSign*advFac        rAdvFac = rkSign*advFac
176        IF (implicitAdvection) rAdvFac = 0. _d 0        IF (implicitAdvection) rAdvFac = rkSign
177    
178        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
179         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
# Line 179  C--   Make local copy of tracer array Line 202  C--   Make local copy of tracer array
202          ENDDO          ENDDO
203        ENDIF        ENDIF
204    
 C--   Unless we have already calculated the advection terms we initialize  
 C     the tendency to zero.  
 C     <== now done earlier at the beginning of thermodynamics.  
 c     IF (calcAdvection) THEN  
 c      DO j=1-Oly,sNy+Oly  
 c       DO i=1-Olx,sNx+Olx  
 c        gTracer(i,j,k,bi,bj)=0. _d 0  
 c       ENDDO  
 c      ENDDO  
 c     ENDIF  
   
205  C--   Pre-calculate del^2 T if bi-harmonic coefficient is non-zero  C--   Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
206        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
207         CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)         CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)
# Line 198  C--   Pre-calculate del^2 T if bi-harmon Line 210  C--   Pre-calculate del^2 T if bi-harmon
210        ENDIF        ENDIF
211    
212  C--   Initialize net flux in X direction  C--   Initialize net flux in X direction
213        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
214         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
215          fZon(i,j) = 0. _d 0          fZon(i,j) = 0. _d 0
216         ENDDO         ENDDO
217        ENDDO        ENDDO
# Line 207  C--   Initialize net flux in X direction Line 219  C--   Initialize net flux in X direction
219  C-    Advective flux in X  C-    Advective flux in X
220        IF (calcAdvection) THEN        IF (calcAdvection) THEN
221          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
222            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)             CALL GAD_C2_ADV_X( bi,bj,k, uTrans, locABT, af, myThid )
223          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
224       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
225            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
226       I            dTtracerLev(k), uTrans, uFld, locABT,       I              deltaTLev(k), uTrans, uFld, locABT,
227       O            af, myThid )       O              af, myThid )
228          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSE
229            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),           DO j=1-OLy,sNy+OLy
230       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,            DO i=1-OLx,sNx+OLx
231       O            af, myThid )  #ifdef ALLOW_OBCS
232          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN             maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
233            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)  #else /* ALLOW_OBCS */
234          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN             maskLocW(i,j) = _maskW(i,j,k,bi,bj)
235            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)  #endif /* ALLOW_OBCS */
236          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            ENDDO
237            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),           ENDDO
238       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
239       O            af, myThid )             CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
240          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       I              uTrans, uFld, maskLocW, locABT,
241           IF ( inAdMode ) THEN       O              af, myThid )
242             ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
243               CALL GAD_U3_ADV_X( bi,bj,k, uTrans, maskLocW, locABT,
244         O              af, myThid )
245             ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
246               CALL GAD_C4_ADV_X( bi,bj,k, uTrans, maskLocW, locABT,
247         O              af, myThid )
248    #ifdef ALLOW_AUTODIFF
249             ELSEIF( advectionScheme.EQ.ENUM_DST3 .OR.
250         &          (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT .AND. inAdMode)
251         &         ) THEN
252  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
253  cph IF inAdExact=.FALSE., we want to use DST3  cph If inAdExact=.FALSE., we want to use DST3
254  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
255            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),  #else /* ALLOW_AUTODIFF */
256       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,           ELSEIF( advectionScheme.EQ.ENUM_DST3 ) THEN
257       O           af, myThid )  #endif /* ALLOW_AUTODIFF */
258               CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
259         I              uTrans, uFld, maskLocW, locABT,
260         O              af, myThid )
261             ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
262               CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
263         I              uTrans, uFld, maskLocW, locABT,
264         O              af, myThid )
265    #ifndef ALLOW_AUTODIFF_TAMC
266             ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
267               CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
268         I              uTrans, uFld, maskLocW, locABT,
269         O              af, myThid )
270    #endif
271           ELSE           ELSE
272            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 )  
273           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)'  
274          ENDIF          ENDIF
275          DO j=1-Oly,sNy+Oly  #ifdef ALLOW_OBCS
276           DO i=1-Olx,sNx+Olx          IF ( useOBCS ) THEN
277    C-      replace advective flux with 1st order upwind scheme estimate
278              CALL OBCS_U1_ADV_TRACER( .TRUE., trIdentity, bi, bj, k,
279         I                             maskW(1-OLx,1-OLy,k,bi,bj),
280         I                             uTrans, locABT,
281         U                             af, myThid )
282            ENDIF
283    #endif /* ALLOW_OBCS */
284            DO j=1-OLy,sNy+OLy
285             DO i=1-OLx,sNx+OLx
286            fZon(i,j) = fZon(i,j) + af(i,j)            fZon(i,j) = fZon(i,j) + af(i,j)
287           ENDDO           ENDDO
288          ENDDO          ENDDO
289  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
290          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
291            diagName = 'ADVx'//diagSufx            diagName = 'ADVx'//diagSufx
292            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
293          ENDIF          ENDIF
294  #endif  #endif
295        ENDIF        ENDIF
# Line 262  C-    Diffusive flux in X Line 298  C-    Diffusive flux in X
298        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
299         CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)         CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
300        ELSE        ELSE
301         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
302          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
303           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
304          ENDDO          ENDDO
305         ENDDO         ENDDO
# Line 276  C-    Add bi-harmonic diffusive flux in Line 312  C-    Add bi-harmonic diffusive flux in
312    
313  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
314  C-    GM/Redi flux in X  C-    GM/Redi flux in X
315        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
316  C *note* should update GMREDI_XTRANSPORT to set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
317          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
318            CALL GMREDI_XTRANSPORT(            CALL GMREDI_XTRANSPORT(
319       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
320       I         xA,TracerN,tracerIdentity,       I         xA,TracerN,trIdentity,
321       U         df,       U         df,
322       I         myThid)       I         myThid)
323          ELSE          ELSE
324            CALL GMREDI_XTRANSPORT(            CALL GMREDI_XTRANSPORT(
325       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
326       I         xA,TracAB, tracerIdentity,       I         xA,TracAB, trIdentity,
327       U         df,       U         df,
328       I         myThid)       I         myThid)
329          ENDIF          ENDIF
330        ENDIF        ENDIF
331  #endif  #endif
332  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
333        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
334         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
335          fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)          fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
336         ENDDO         ENDDO
337        ENDDO        ENDDO
# Line 304  C     anelastic: advect.fluxes are scale Line 340  C     anelastic: advect.fluxes are scale
340  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
341  C       excluding advective terms:  C       excluding advective terms:
342        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
343       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
344            diagName = 'DFxE'//diagSufx            diagName = 'DFxE'//diagSufx
345            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
346        ENDIF        ENDIF
347  #endif  #endif
348    
349  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
350        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
351         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
352          fMer(i,j) = 0. _d 0          fMer(i,j) = 0. _d 0
353         ENDDO         ENDDO
354        ENDDO        ENDDO
# Line 320  C--   Initialize net flux in Y direction Line 356  C--   Initialize net flux in Y direction
356  C-    Advective flux in Y  C-    Advective flux in Y
357        IF (calcAdvection) THEN        IF (calcAdvection) THEN
358          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
359            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)             CALL GAD_C2_ADV_Y( bi,bj,k, vTrans, locABT, af, myThid )
360          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
361       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
362            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
363       I            dTtracerLev(k), vTrans, vFld, locABT,       I              deltaTLev(k), vTrans, vFld, locABT,
364       O            af, myThid )       O              af, myThid )
365          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSE
366            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),           DO j=1-OLy,sNy+OLy
367       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,            DO i=1-OLx,sNx+OLx
368       O            af, myThid )  #ifdef ALLOW_OBCS
369          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN             maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
370            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)  #else /* ALLOW_OBCS */
371          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN             maskLocS(i,j) = _maskS(i,j,k,bi,bj)
372            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)  #endif /* ALLOW_OBCS */
373          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            ENDDO
374            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),           ENDDO
375       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
376       O            af, myThid )             CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
377          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       I              vTrans, vFld, maskLocS, locABT,
378           IF ( inAdMode ) THEN       O              af, myThid )
379             ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
380               CALL GAD_U3_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
381         O              af, myThid )
382             ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
383               CALL GAD_C4_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
384         O              af, myThid )
385    #ifdef ALLOW_AUTODIFF
386             ELSEIF( advectionScheme.EQ.ENUM_DST3 .OR.
387         &          (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT .AND. inAdMode)
388         &         ) THEN
389  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
390  cph IF inAdExact=.FALSE., we want to use DST3  cph If inAdExact=.FALSE., we want to use DST3
391  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
392            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),  #else /* ALLOW_AUTODIFF */
393       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,           ELSEIF( advectionScheme.EQ.ENUM_DST3 ) THEN
394       O           af, myThid )  #endif /* ALLOW_AUTODIFF */
395               CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
396         I              vTrans, vFld, maskLocS, locABT,
397         O              af, myThid )
398             ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
399               CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
400         I              vTrans, vFld, maskLocS, locABT,
401         O              af, myThid )
402    #ifndef ALLOW_AUTODIFF_TAMC
403             ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
404               CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
405         I              vTrans, vFld, maskLocS, locABT,
406         O              af, myThid )
407    #endif
408           ELSE           ELSE
409            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 )  
410           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)'  
411          ENDIF          ENDIF
412          DO j=1-Oly,sNy+Oly  #ifdef ALLOW_OBCS
413           DO i=1-Olx,sNx+Olx          IF ( useOBCS ) THEN
414    C-      replace advective flux with 1st order upwind scheme estimate
415              CALL OBCS_U1_ADV_TRACER( .FALSE., trIdentity, bi, bj, k,
416         I                             maskS(1-OLx,1-OLy,k,bi,bj),
417         I                             vTrans, locABT,
418         U                             af, myThid )
419            ENDIF
420    #endif /* ALLOW_OBCS */
421            DO j=1-OLy,sNy+OLy
422             DO i=1-OLx,sNx+OLx
423            fMer(i,j) = fMer(i,j) + af(i,j)            fMer(i,j) = fMer(i,j) + af(i,j)
424           ENDDO           ENDDO
425          ENDDO          ENDDO
426  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
427          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
428            diagName = 'ADVy'//diagSufx            diagName = 'ADVy'//diagSufx
429            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
430          ENDIF          ENDIF
431  #endif  #endif
432        ENDIF        ENDIF
# Line 375  C-    Diffusive flux in Y Line 435  C-    Diffusive flux in Y
435        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
436         CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)         CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
437        ELSE        ELSE
438         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
439          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
440           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
441          ENDDO          ENDDO
442         ENDDO         ENDDO
# Line 389  C-    Add bi-harmonic flux in Y Line 449  C-    Add bi-harmonic flux in Y
449    
450  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
451  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
452        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
453  C *note* should update GMREDI_YTRANSPORT to set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
454          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
455            CALL GMREDI_YTRANSPORT(            CALL GMREDI_YTRANSPORT(
456       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
457       I         yA,TracerN,tracerIdentity,       I         yA,TracerN,trIdentity,
458       U         df,       U         df,
459       I         myThid)       I         myThid)
460          ELSE          ELSE
461            CALL GMREDI_YTRANSPORT(            CALL GMREDI_YTRANSPORT(
462       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
463       I         yA,TracAB, tracerIdentity,       I         yA,TracAB, trIdentity,
464       U         df,       U         df,
465       I         myThid)       I         myThid)
466          ENDIF          ENDIF
467        ENDIF        ENDIF
468  #endif  #endif
469  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
470        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
471         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
472          fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)          fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
473         ENDDO         ENDDO
474        ENDDO        ENDDO
# Line 417  C     anelastic: advect.fluxes are scale Line 477  C     anelastic: advect.fluxes are scale
477  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
478  C       excluding advective terms:  C       excluding advective terms:
479        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
480       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
481            diagName = 'DFyE'//diagSufx            diagName = 'DFyE'//diagSufx
482            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
483        ENDIF        ENDIF
484  #endif  #endif
485    
# Line 428  C-    Advective flux in R Line 488  C-    Advective flux in R
488  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
489  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
490        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.
491       &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.k.LT.Nr)       &     (.NOT.useAIM .OR. trIdentity.NE.GAD_SALINITY .OR. k.LT.Nr)
492       &   ) THEN       &   ) THEN
493  #else  #else
494        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN
495  #endif  #endif
496  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
497          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
498            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)             CALL GAD_C2_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
499          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
500       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
501            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,             CALL GAD_DST2U1_ADV_R( bi,bj,k,vertAdvecScheme,deltaTLev(k),
502       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
503       O         af, myThid )       O              af, myThid )
504          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
505            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
506       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
507       O         af, myThid )       O              af, myThid )
508          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
509            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)             CALL GAD_U3_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
510          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
511            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)             CALL GAD_C4_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
512          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN  #ifdef ALLOW_AUTODIFF
513            CALL GAD_DST3_ADV_R( bi,bj,k,          ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 .OR.
514       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       &         (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT .AND. inAdMode)
515       O         af, myThid )       &        ) THEN
         ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN  
516  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
517  cph IF inAdExact=.FALSE., we want to use DST3  cph If inAdExact=.FALSE., we want to use DST3
518  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
519            IF ( inAdMode ) THEN  #else /* ALLOW_AUTODIFF */
520             CALL GAD_DST3_ADV_R( bi,bj,k,          ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
521       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),  #endif /* ALLOW_AUTODIFF */
522       O         af, myThid )             CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
523            ELSE       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
524             CALL GAD_DST3FL_ADV_R( bi,bj,k,       O              af, myThid )
525       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
526       O         af, myThid )             CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
527            ENDIF       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
528         O              af, myThid )
529    #ifndef ALLOW_AUTODIFF_TAMC
530          ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
531             CALL GAD_OS7MP_ADV_R( bi,bj,k,             CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
532       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
533       O         af, myThid )       O              af, myThid )
534    #endif
535          ELSE          ELSE
536            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
537          ENDIF          ENDIF
538  C-     add the advective flux to fVerT  C-     add the advective flux to fVerT
539          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
540           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
541            fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)            fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)*maskInC(i,j,bi,bj)
542           ENDDO           ENDDO
543          ENDDO          ENDDO
544  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
545          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
546            diagName = 'ADVr'//diagSufx            diagName = 'ADVr'//diagSufx
547            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
548  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
549  C        does it only if k=1 (never the case here)  C        does it only if k=1 (never the case here)
550            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
# Line 494  C-    Diffusive flux in R Line 556  C-    Diffusive flux in R
556  C     Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper  C     Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
557  C           boundary condition.  C           boundary condition.
558        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
559         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
560          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
561           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
562          ENDDO          ENDDO
563         ENDDO         ENDDO
# Line 507  C           boundary condition. Line 569  C           boundary condition.
569         ENDIF         ENDIF
570        ENDIF        ENDIF
571    
572          IF ( trUseDiffKr4 ) THEN
573           IF ( applyAB_onTracer ) THEN
574             CALL GAD_BIHARM_R( bi,bj,k, diffKr4, TracerN, df, myThid )
575           ELSE
576             CALL GAD_BIHARM_R( bi,bj,k, diffKr4, TracAB,  df, myThid )
577           ENDIF
578          ENDIF
579    
580  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
581  C-    GM/Redi flux in R  C-    GM/Redi flux in R
582        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
583  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
584          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
585            CALL GMREDI_RTRANSPORT(            CALL GMREDI_RTRANSPORT(
586       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
587       I         TracerN,tracerIdentity,       I         TracerN,trIdentity,
588       U         df,       U         df,
589       I         myThid)       I         myThid)
590          ELSE          ELSE
591            CALL GMREDI_RTRANSPORT(            CALL GMREDI_RTRANSPORT(
592       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
593       I         TracAB, tracerIdentity,       I         TracAB, trIdentity,
594       U         df,       U         df,
595       I         myThid)       I         myThid)
596          ENDIF          ENDIF
597        ENDIF        ENDIF
598  #endif  #endif
599    
600        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
601         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
602          fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)          fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
603         ENDDO         ENDDO
604        ENDDO        ENDDO
# Line 537  C *note* should update GMREDI_RTRANSPORT Line 607  C *note* should update GMREDI_RTRANSPORT
607  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
608  C       Explicit terms only & excluding advective terms:  C       Explicit terms only & excluding advective terms:
609        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
610       &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN       &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN
611            diagName = 'DFrE'//diagSufx            diagName = 'DFrE'//diagSufx
612            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
613        ENDIF        ENDIF
614  #endif  #endif
615    
616  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
617  C-    Set non local KPP transport term (ghat):  C-    Set non local KPP transport term (ghat):
618        IF ( useKPP .AND. k.GE.2 ) THEN        IF ( trUseKPP .AND. k.GE.2 ) THEN
619         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
620          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
621           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
622          ENDDO          ENDDO
623         ENDDO         ENDDO
624         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (trIdentity.EQ.GAD_TEMPERATURE) THEN
625          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
626       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
627       O     df )       O           df,
628         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN       I           myTime, myIter, myThid )
629           ELSEIF (trIdentity.EQ.GAD_SALINITY) THEN
630          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
631       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
632       O     df )       O           df,
633         I           myTime, myIter, myThid )
634  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
635         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN         ELSEIF (trIdentity .GE. GAD_TR1) THEN
636          CALL KPP_TRANSPORT_PTR(          CALL KPP_TRANSPORT_PTR(
637       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
638       I     tracerIdentity-GAD_TR1+1,       I           trIdentity-GAD_TR1+1,
639       O     df )       O           df,
640         I           myTime, myIter, myThid )
641  #endif  #endif
642         ELSE         ELSE
643          PRINT*,'invalid tracer indentity: ', tracerIdentity          WRITE(errorMessageUnit,*)
644          STOP 'GAD_CALC_RHS: Ooops'       &    'tracer identity =', trIdentity, ' is not valid => STOP'
645            STOP 'ABNORMAL END: S/R GAD_CALC_RHS: invalid tracer identity'
646         ENDIF         ENDIF
647         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
648          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
649           fVerT(i,j,kUp) = fVerT(i,j,kUp)           fVerT(i,j,kUp) = fVerT(i,j,kUp)
650       &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)       &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
651          ENDDO          ENDDO
652         ENDDO         ENDDO
653    #ifdef ALLOW_DIAGNOSTICS
654    C-    Diagnostics of Non-Local Tracer (vertical) flux
655           IF ( useDiagnostics ) THEN
656             diagName = 'KPPg'//diagSufx
657             CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
658    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
659    C        does it only if k=1 (never the case here)
660             IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
661           ENDIF
662    #endif
663          ENDIF
664    #endif /* ALLOW_KPP */
665    
666    #ifdef GAD_SMOLARKIEWICZ_HACK
667    coj   Hack to make redi (and everything else in this s/r) positive
668    coj   (see Smolarkiewicz MWR 1989 and Bott MWR 1989).
669    coj   Only works if 'down' is k+1 and k loop in thermodynamics is k=Nr,1,-1
670    coj
671    coj   Apply to all tracers except temperature
672          IF ( trIdentity.NE.GAD_TEMPERATURE .AND.
673         &     trIdentity.NE.GAD_SALINITY ) THEN
674           DO j=1-OLy,sNy+OLy-1
675            DO i=1-OLx,sNx+OLx-1
676    coj   Add outgoing fluxes
677             outFlux=deltaTLev(k)*
678         &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
679         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
680         &    *( MAX(0. _d 0,fZon(i+1,j)) + MAX(0. _d 0,-fZon(i,j))
681         &      +MAX(0. _d 0,fMer(i,j+1)) + MAX(0. _d 0,-fMer(i,j))
682         &      +MAX(0. _d 0,fVerT(i,j,kDown)*rkSign)
683         &      +MAX(0. _d 0,-fVerT(i,j,kUp)*rkSign)
684         &     )
685             IF ( applyAB_onTracer ) THEN
686               trac=TracerN(i,j,k,bi,bj)
687             ELSE
688               trac=TracAB(i,j,k,bi,bj)
689             ENDIF
690    coj   If they would reduce tracer by a fraction of more than
691    coj   SmolarkiewiczMaxFrac, scale them down
692             IF (outFlux.GT.0. _d 0 .AND.
693         &       outFlux.GT.SmolarkiewiczMaxFrac*trac) THEN
694    coj   If tracer is already negative, scale flux to zero
695               fac = MAX(0. _d 0,SmolarkiewiczMaxFrac*trac/outFlux)
696    
697               IF (fZon(i+1,j).GT.0. _d 0) fZon(i+1,j)=fac*fZon(i+1,j)
698               IF (-fZon(i,j) .GT.0. _d 0) fZon(i,j)  =fac*fZon(i,j)
699               IF (fMer(i,j+1).GT.0. _d 0) fMer(i,j+1)=fac*fMer(i,j+1)
700               IF (-fMer(i,j) .GT.0. _d 0) fMer(i,j)  =fac*fMer(i,j)
701               IF (-fVerT(i,j,kUp)*rkSign .GT.0. _d 0)
702         &       fVerT(i,j,kUp)=fac*fVerT(i,j,kUp)
703    
704               IF (k.LT.Nr .AND. fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN
705    coj   Down flux is special: it has already been applied in lower layer,
706    coj   so we have to readjust this.
707    coj   Note: for k+1, gTracer is now the updated tracer, not the tendency!
708    coj   thus it has an extra factor deltaTLev(k+1)
709                 gTrFac=deltaTLev(k+1)
710    coj   Other factors that have been applied to gTracer since the last call:
711    #ifdef NONLIN_FRSURF
712                 IF (nonlinFreeSurf.GT.0) THEN
713                  IF (select_rStar.GT.0) THEN
714    #ifndef DISABLE_RSTAR_CODE
715                    gTrFac = gTrFac/rStarExpC(i,j,bi,bj)
716    #endif /* DISABLE_RSTAR_CODE */
717                  ENDIF
718                 ENDIF
719    #endif /* NONLIN_FRSURF */
720    coj   Now: undo down flux, ...
721                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
722         &        +gTrFac
723         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
724         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
725         &         *recip_rhoFacC(k+1)
726         &         *( -fVerT(i,j,kDown)*rkSign )
727    coj   ... scale ...
728                 fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)
729    coj   ... and reapply
730                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
731         &        +gTrFac
732         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
733         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
734         &         *recip_rhoFacC(k+1)
735         &         *( fVerT(i,j,kDown)*rkSign )
736               ENDIF
737    
738             ENDIF
739            ENDDO
740           ENDDO
741        ENDIF        ENDIF
742  #endif  #endif
743    
744  C--   Divergence of fluxes  C--   Divergence of fluxes
745  C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged  C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
746        DO j=1-Oly,sNy+Oly-1  C     for Stevens OBC: keep only vertical diffusive contribution on boundaries
747         DO i=1-Olx,sNx+Olx-1        DO j=1-OLy,sNy+OLy-1
748           DO i=1-OLx,sNx+OLx-1
749          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
750       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
751       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
752       &   *( (fZon(i+1,j)-fZon(i,j))       &   *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)
753       &     +(fMer(i,j+1)-fMer(i,j))       &     +(fMer(i,j+1)-fMer(i,j))*maskInC(i,j,bi,bj)
754       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
755       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
756       &                   +(vTrans(i,j+1)-vTrans(i,j))       &                   +(vTrans(i,j+1)-vTrans(i,j))*advFac
757       &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac       &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
758       &                  )*advFac       &                  )*maskInC(i,j,bi,bj)
759       &    )       &    )
760         ENDDO         ENDDO
761        ENDDO        ENDDO
762    
763  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
764        IF ( debugLevel .GE. debLevB        IF ( debugLevel .GE. debLevC
765       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE       &   .AND. trIdentity.EQ.GAD_TEMPERATURE
766       &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0       &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
767       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
768       &   .AND. useCubedSphereExchange ) THEN       &   .AND. useCubedSphereExchange ) THEN
Legend:
Removed from v.1.44  
changed lines
  Added in v.1.63
  ViewVC Help
Powered by ViewVC 1.1.22