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