/[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.28 by jmc, Mon Sep 27 14:45:09 2004 UTC revision 1.63 by jmc, Fri Dec 27 15:52:17 2013 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: GAD_CALC_RHS  C !ROUTINE: GAD_CALC_RHS
8    
9  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
10        SUBROUTINE GAD_CALC_RHS(        SUBROUTINE GAD_CALC_RHS(
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,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I           xA, yA, maskUp, uFld, vFld, wFld,
13       I           uVel, vVel, wVel,       I           uTrans, vTrans, rTrans, rTransKp1,
14       I           diffKh, diffK4, KappaRT, Tracer,       I           diffKh, diffK4, KappaR, diffKr4, TracerN, TracAB,
15       I           tracerIdentity, advectionScheme, vertAdvecScheme,       I           deltaTLev, trIdentity,
16       I           calcAdvection, implicitAdvection,       I           advectionScheme, vertAdvecScheme,
17         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    
23  C !DESCRIPTION:  C !DESCRIPTION:
24  C Calculates the tendancy of a tracer due to advection and diffusion.  C Calculates the tendency of a tracer due to advection and diffusion.
25  C It calculates the fluxes in each direction indepentently and then  C It calculates the fluxes in each direction indepentently and then
26  C sets the tendancy to the divergence of these fluxes. The advective  C sets the tendency to the divergence of these fluxes. The advective
27  C fluxes are only calculated here when using the linear advection schemes  C fluxes are only calculated here when using the linear advection schemes
28  C otherwise only the diffusive and parameterized fluxes are calculated.  C otherwise only the diffusive and parameterized fluxes are calculated.
29  C  C
# Line 29  C \begin{equation*} Line 32  C \begin{equation*}
32  C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}  C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
33  C \end{equation*}  C \end{equation*}
34  C  C
35  C The tendancy is the divergence of the fluxes:  C The tendency is the divergence of the fluxes:
36  C \begin{equation*}  C \begin{equation*}
37  C G_\theta = G_\theta + \nabla \cdot {\bf F}  C G_\theta = G_\theta + \nabla \cdot {\bf F}
38  C \end{equation*}  C \end{equation*}
39  C  C
40  C The tendancy is assumed to contain data on entry.  C The tendency is assumed to contain data on entry.
41    
42  C !USES: ===============================================================  C !USES: ===============================================================
43        IMPLICIT NONE        IMPLICIT NONE
# 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: ===================================================
60  C bi,bj            :: tile indices  C bi,bj            :: tile indices
61  C iMin,iMax        :: loop range for called routines  C iMin,iMax        :: loop range for called routines
62  C jMin,jMax        :: loop range for called routines  C jMin,jMax        :: loop range for called routines
63  C kup              :: index into 2 1/2D array, toggles between 1|2  C k                :: vertical index
64  C kdown            :: index into 2 1/2D array, toggles between 2|1  C kM1              :: =k-1 for k>1, =1 for k=1
65  C kp1              :: =k+1 for k<Nr, =Nr for k=Nr  C kUp              :: index into 2 1/2D array, toggles between 1|2
66    C kDown            :: index into 2 1/2D array, toggles between 2|1
67  C xA,yA            :: areas of X and Y face of tracer cells  C xA,yA            :: areas of X and Y face of tracer cells
68    C maskUp           :: 2-D array for mask at W points
69    C uFld,vFld,wFld   :: Local copy of velocity field (3 components)
70  C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points  C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points
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
 C maskUp           :: 2-D array for mask at W points  
 C uVel,vVel,wVel   :: 3 components of the velcity field (3-D array)  
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 KappaRT          :: 3-D array for vertical diffusion coefficient  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
76  C Tracer           :: tracer field  C diffKr4          :: 1-D array for vertical bi-harmonic diffusion coefficient
77  C tracerIdentity   :: tracer identifier (required for KPP,GM)  C TracerN          :: tracer field @ time-step n (Note: only used
78    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
80    C                     or extrapolated fwd in time to n+1/2 if applying AB on Tr)
81    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)
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 79  C myThid           :: thread number Line 94  C myThid           :: thread number
94        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
95        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97          _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98          _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99          _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100          _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL uVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
       _RL vVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
       _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
105        _RL diffKh, diffK4        _RL diffKh, diffK4
106        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL diffKr4(Nr)
108        INTEGER tracerIdentity        _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)
110          _RL deltaTLev(Nr)
111          INTEGER trIdentity
112        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
113        LOGICAL calcAdvection        LOGICAL calcAdvection
114        LOGICAL implicitAdvection        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          :: tendancy 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
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)
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 127  C--   the kDown is still required Line 163  C--   the kDown is still required
163        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
164  #endif  #endif
165    
166    #ifdef ALLOW_DIAGNOSTICS
167    C--   Set diagnostic suffix for the current tracer
168          IF ( useDiagnostics ) THEN
169            diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
170          ENDIF
171    #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 = rkFac*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 143  C--   the kDown is still required Line 186  C--   the kDown is still required
186        ENDDO        ENDDO
187    
188  C--   Make local copy of tracer array  C--   Make local copy of tracer array
189        DO j=1-OLy,sNy+OLy        IF ( applyAB_onTracer ) THEN
190         DO i=1-OLx,sNx+OLx          DO j=1-OLy,sNy+OLy
191          localT(i,j)=tracer(i,j,k,bi,bj)           DO i=1-OLx,sNx+OLx
192         ENDDO            localT(i,j)=TracerN(i,j,k,bi,bj)
193        ENDDO            locABT(i,j)= TracAB(i,j,k,bi,bj)
194             ENDDO
195  C--   Unless we have already calculated the advection terms we initialize          ENDDO
196  C     the tendency to zero.        ELSE
197  C     <== now done earlier at the beginning of thermodynamics.          DO j=1-OLy,sNy+OLy
198  c     IF (calcAdvection) THEN           DO i=1-OLx,sNx+OLx
199  c      DO j=1-Oly,sNy+Oly            localT(i,j)= TracAB(i,j,k,bi,bj)
200  c       DO i=1-Olx,sNx+Olx            locABT(i,j)= TracAB(i,j,k,bi,bj)
201  c        gTracer(i,j,k,bi,bj)=0. _d 0           ENDDO
202  c       ENDDO          ENDDO
203  c      ENDDO        ENDIF
 c     ENDIF  
204    
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
# Line 168  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
218    
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,localT,af,myThid)             CALL GAD_C2_ADV_X( bi,bj,k, uTrans, locABT, af, myThid )
223        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
224         CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, deltaTtracer,       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
225       I          uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,             CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
226       O          af, myThid )       I              deltaTLev(k), uTrans, uFld, locABT,
227        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       O              af, myThid )
228         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)          ELSE
229        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN           DO j=1-OLy,sNy+OLy
230         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            DO i=1-OLx,sNx+OLx
231        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN  #ifdef ALLOW_OBCS
232         CALL GAD_DST3_ADV_X( bi,bj,k, deltaTtracer,             maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
233       I          uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,  #else /* ALLOW_OBCS */
234       O          af, myThid )             maskLocW(i,j) = _maskW(i,j,k,bi,bj)
235        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN  #endif /* ALLOW_OBCS */
236         CALL GAD_DST3FL_ADV_X( bi,bj,k, deltaTtracer,            ENDDO
237       I          uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,           ENDDO
238       O          af, myThid )           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
239        ELSE             CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
240         STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'       I              uTrans, uFld, maskLocW, locABT,
241        ENDIF       O              af, myThid )
242        DO j=1-Oly,sNy+Oly           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
243         DO i=1-Olx,sNx+Olx             CALL GAD_U3_ADV_X( bi,bj,k, uTrans, maskLocW, locABT,
244          fZon(i,j) = fZon(i,j) + af(i,j)       O              af, myThid )
245         ENDDO           ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
246        ENDDO             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:
253    cph If inAdExact=.FALSE., we want to use DST3
254    cph with limiters in forward, but without limiters in reverse.
255    #else /* ALLOW_AUTODIFF */
256             ELSEIF( advectionScheme.EQ.ENUM_DST3 ) THEN
257    #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
272              STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
273             ENDIF
274            ENDIF
275    #ifdef ALLOW_OBCS
276            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)
287             ENDDO
288            ENDDO
289    #ifdef ALLOW_DIAGNOSTICS
290            IF ( useDiagnostics ) THEN
291              diagName = 'ADVx'//diagSufx
292              CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
293            ENDIF
294    #endif
295        ENDIF        ENDIF
296    
297  C-    Diffusive flux in X  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
306        ENDIF        ENDIF
307    
308    C-    Add bi-harmonic diffusive flux in X
309          IF (diffK4 .NE. 0.) THEN
310           CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
311          ENDIF
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 use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
317          CALL GMREDI_XTRANSPORT(          IF ( applyAB_onTracer ) THEN
318       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_XTRANSPORT(
319       I     xA,Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
320       U     df,       I         xA,TracerN,trIdentity,
321       I     myThid)       U         df,
322         I         myThid)
323            ELSE
324              CALL GMREDI_XTRANSPORT(
325         I         iMin,iMax,jMin,jMax,bi,bj,k,
326         I         xA,TracAB, trIdentity,
327         U         df,
328         I         myThid)
329            ENDIF
330        ENDIF        ENDIF
331  #endif  #endif
332        DO j=1-Oly,sNy+Oly  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
333         DO i=1-Olx,sNx+Olx        DO j=1-OLy,sNy+OLy
334          fZon(i,j) = fZon(i,j) + df(i,j)         DO i=1-OLx,sNx+OLx
335            fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
336         ENDDO         ENDDO
337        ENDDO        ENDDO
338    
339  C-    Bi-harmonic duffusive flux in X  #ifdef ALLOW_DIAGNOSTICS
340        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
341         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)  C       excluding advective terms:
342         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
343          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
344           fZon(i,j) = fZon(i,j) + df(i,j)            diagName = 'DFxE'//diagSufx
345          ENDDO            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
        ENDDO  
346        ENDIF        ENDIF
347    #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
355    
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,localT,af,myThid)             CALL GAD_C2_ADV_Y( bi,bj,k, vTrans, locABT, af, myThid )
360        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
361         CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, deltaTtracer,       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
362       I          vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,             CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
363       O          af, myThid )       I              deltaTLev(k), vTrans, vFld, locABT,
364        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       O              af, myThid )
365         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)          ELSE
366        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN           DO j=1-OLy,sNy+OLy
367         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            DO i=1-OLx,sNx+OLx
368        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN  #ifdef ALLOW_OBCS
369         CALL GAD_DST3_ADV_Y( bi,bj,k, deltaTtracer,             maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
370       I          vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,  #else /* ALLOW_OBCS */
371       O          af, myThid )             maskLocS(i,j) = _maskS(i,j,k,bi,bj)
372        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN  #endif /* ALLOW_OBCS */
373         CALL GAD_DST3FL_ADV_Y( bi,bj,k, deltaTtracer,            ENDDO
374       I          vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,           ENDDO
375       O          af, myThid )           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
376        ELSE             CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
377         STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'       I              vTrans, vFld, maskLocS, locABT,
378        ENDIF       O              af, myThid )
379        DO j=1-Oly,sNy+Oly           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
380         DO i=1-Olx,sNx+Olx             CALL GAD_U3_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
381          fMer(i,j) = fMer(i,j) + af(i,j)       O              af, myThid )
382         ENDDO           ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
383        ENDDO             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:
390    cph If inAdExact=.FALSE., we want to use DST3
391    cph with limiters in forward, but without limiters in reverse.
392    #else /* ALLOW_AUTODIFF */
393             ELSEIF( advectionScheme.EQ.ENUM_DST3 ) THEN
394    #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
409               STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
410             ENDIF
411            ENDIF
412    #ifdef ALLOW_OBCS
413            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)
424             ENDDO
425            ENDDO
426    #ifdef ALLOW_DIAGNOSTICS
427            IF ( useDiagnostics ) THEN
428              diagName = 'ADVy'//diagSufx
429              CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
430            ENDIF
431    #endif
432        ENDIF        ENDIF
433    
434  C-    Diffusive flux in Y  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
443        ENDIF        ENDIF
444    
445    C-    Add bi-harmonic flux in Y
446          IF (diffK4 .NE. 0.) THEN
447           CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
448          ENDIF
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 use localT and set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
454         CALL GMREDI_YTRANSPORT(          IF ( applyAB_onTracer ) THEN
455       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_YTRANSPORT(
456       I     yA,Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
457       U     df,       I         yA,TracerN,trIdentity,
458       I     myThid)       U         df,
459         I         myThid)
460            ELSE
461              CALL GMREDI_YTRANSPORT(
462         I         iMin,iMax,jMin,jMax,bi,bj,k,
463         I         yA,TracAB, trIdentity,
464         U         df,
465         I         myThid)
466            ENDIF
467        ENDIF        ENDIF
468  #endif  #endif
469        DO j=1-Oly,sNy+Oly  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
470         DO i=1-Olx,sNx+Olx        DO j=1-OLy,sNy+OLy
471          fMer(i,j) = fMer(i,j) + df(i,j)         DO i=1-OLx,sNx+OLx
472            fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
473         ENDDO         ENDDO
474        ENDDO        ENDDO
475    
476  C-    Bi-harmonic flux in Y  #ifdef ALLOW_DIAGNOSTICS
477        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
478         CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)  C       excluding advective terms:
479         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
480          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
481           fMer(i,j) = fMer(i,j) + df(i,j)            diagName = 'DFyE'//diagSufx
482          ENDDO            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
        ENDDO  
483        ENDIF        ENDIF
484    #endif
485    
486  C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):  C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
487  C-    Advective flux in R  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,tracer,af,myThid)             CALL GAD_C2_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
499         ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
500          CALL GAD_FLUXLIMIT_ADV_R(       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
501       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)             CALL GAD_DST2U1_ADV_R( bi,bj,k,vertAdvecScheme,deltaTLev(k),
502         ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
503          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)       O              af, myThid )
504         ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
505          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
506         ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
507          CALL GAD_DST3_ADV_R(       O              af, myThid )
508       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
509         ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN             CALL GAD_U3_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
510          CALL GAD_DST3FL_ADV_R(          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
511       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)             CALL GAD_C4_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
512         ELSE  #ifdef ALLOW_AUTODIFF
513          STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'          ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 .OR.
514         ENDIF       &         (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT .AND. inAdMode)
515         &        ) THEN
516    cph This block is to trick the adjoint:
517    cph If inAdExact=.FALSE., we want to use DST3
518    cph with limiters in forward, but without limiters in reverse.
519    #else /* ALLOW_AUTODIFF */
520            ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
521    #endif /* ALLOW_AUTODIFF */
522               CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
523         I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
524         O              af, myThid )
525            ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
526               CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
527         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
531               CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
532         I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
533         O              af, myThid )
534    #endif
535            ELSE
536              STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
537            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
543          ENDDO          ENDDO
544         ENDDO  #ifdef ALLOW_DIAGNOSTICS
545            IF ( useDiagnostics ) THEN
546              diagName = 'ADVr'//diagSufx
547              CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
548    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
549    C        does it only if k=1 (never the case here)
550              IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
551            ENDIF
552    #endif
553        ENDIF        ENDIF
554    
555  C-    Diffusive flux in R  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
564        ELSE        ELSE
565         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         IF ( applyAB_onTracer ) THEN
566             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracerN,df,myThid)
567           ELSE
568             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracAB, df,myThid)
569           ENDIF
570          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        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         CALL GMREDI_RTRANSPORT(          IF ( applyAB_onTracer ) THEN
585       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_RTRANSPORT(
586       I     Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
587       U     df,       I         TracerN,trIdentity,
588       I     myThid)       U         df,
589         I         myThid)
590            ELSE
591              CALL GMREDI_RTRANSPORT(
592         I         iMin,iMax,jMin,jMax,bi,bj,k,
593         I         TracAB, trIdentity,
594         U         df,
595         I         myThid)
596            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
605    
606    #ifdef ALLOW_DIAGNOSTICS
607    C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
608    C       Explicit terms only & excluding advective terms:
609          IF ( useDiagnostics .AND.
610         &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN
611              diagName = 'DFrE'//diagSufx
612              CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
613          ENDIF
614    #endif
615    
616  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
617  C-    Add non local KPP transport term (ghat) to diffusive T flux.  C-    Set non local KPP transport term (ghat):
618        IF (useKPP) 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
 C *note* should update KPP_TRANSPORT_T to set df  *aja*  
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       I     KappaRT,       O           df,
628       U     df )       I           myTime, myIter, myThid )
629         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         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       I     KappaRT,       O           df,
633       U     df )       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,KappaRT,       I           trIdentity-GAD_TR1+1,
639       U     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) + df(i,j)*maskUp(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp)
650         &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
651            ENDDO
652           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          ENDDO
740         ENDDO         ENDDO
741        ENDIF        ENDIF
742  #endif  #endif
743    
744  C--   Divergence of fluxes  C--   Divergence of fluxes
745        DO j=1-Oly,sNy+Oly-1  C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
746         DO i=1-Olx,sNx+Olx-1  C     for Stevens OBC: keep only vertical diffusive contribution on boundaries
747          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_rA(i,j,bi,bj)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
751       &   *( (fZon(i+1,j)-fZon(i,j))       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
752       &     +(fMer(i,j+1)-fMer(i,j))       &   *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)
753       &     +(fVerT(i,j,kUp)-fVerT(i,j,kDown))*rkFac       &     +(fMer(i,j+1)-fMer(i,j))*maskInC(i,j,bi,bj)
754       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
755       &                   +(vTrans(i,j+1)-vTrans(i,j))       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
756       &                   +(rTrans(i,j)-rTransKp1(i,j))*rAdvFac       &                   +(vTrans(i,j+1)-vTrans(i,j))*advFac
757       &                  )*advFac       &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
758         &                  )*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
# Line 449  C--   Divergence of fluxes Line 770  C--   Divergence of fluxes
770       &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )       &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
771        ENDIF        ENDIF
772  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
773    
774        RETURN        RETURN
775        END        END
Legend:
Removed from v.1.28  
changed lines
  Added in v.1.63
  ViewVC Help
Powered by ViewVC 1.1.22