/[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.2 by jmc, Thu Jul 12 00:26:30 2001 UTC revision 1.58 by jmc, Thu Dec 1 14:16:30 2011 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6        SUBROUTINE GAD_CALC_RHS(  CBOP
7    C !ROUTINE: GAD_CALC_RHS
8    
9    C !INTERFACE: ==========================================================
10          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,maskUp,       I           xA, yA, maskUp, uFld, vFld, wFld,
13       I           diffKh, diffK4, KappaRT, Tracer,       I           uTrans, vTrans, rTrans, rTransKp1,
14       I           tracerIdentity,       I           diffKh, diffK4, KappaR, TracerN, TracAB,
15         I           deltaTLev, tracerIdentity,
16         I           advectionScheme, vertAdvecScheme,
17         I           calcAdvection, implicitAdvection, applyAB_onTracer,
18         I           trUseGMRedi, trUseKPP,
19       U           fVerT, gTracer,       U           fVerT, gTracer,
20       I           myThid )       I           myTime, myIter, myThid )
 C     /==========================================================\  
 C     | SUBROUTINE GAD_CALC_RHS                                  |  
 C     |==========================================================|  
 C     \==========================================================/  
       IMPLICIT NONE  
21    
22  C     == GLobal variables ==  C !DESCRIPTION:
23    C Calculates the tendency of a tracer due to advection and diffusion.
24    C It calculates the fluxes in each direction indepentently and then
25    C sets the tendency to the divergence of these fluxes. The advective
26    C fluxes are only calculated here when using the linear advection schemes
27    C otherwise only the diffusive and parameterized fluxes are calculated.
28    C
29    C Contributions to the flux are calculated and added:
30    C \begin{equation*}
31    C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
32    C \end{equation*}
33    C
34    C The tendency is the divergence of the fluxes:
35    C \begin{equation*}
36    C G_\theta = G_\theta + \nabla \cdot {\bf F}
37    C \end{equation*}
38    C
39    C The tendency is assumed to contain data on entry.
40    
41    C !USES: ===============================================================
42          IMPLICIT NONE
43  #include "SIZE.h"  #include "SIZE.h"
44  #include "EEPARAMS.h"  #include "EEPARAMS.h"
45  #include "PARAMS.h"  #include "PARAMS.h"
46  #include "GRID.h"  #include "GRID.h"
47  #include "DYNVARS.h"  #include "SURFACE.h"
48  #include "GAD.h"  #include "GAD.h"
49    
50  C     == Routine arguments ==  #ifdef ALLOW_AUTODIFF_TAMC
51        INTEGER k,kUp,kDown,kM1  #include "tamc.h"
52    #include "tamc_keys.h"
53    #endif /* ALLOW_AUTODIFF_TAMC */
54    
55    C !INPUT PARAMETERS: ===================================================
56    C bi,bj            :: tile indices
57    C iMin,iMax        :: loop range for called routines
58    C jMin,jMax        :: loop range for called routines
59    C k                :: vertical index
60    C kM1              :: =k-1 for k>1, =1 for k=1
61    C kUp              :: index into 2 1/2D array, toggles between 1|2
62    C kDown            :: index into 2 1/2D array, toggles between 2|1
63    C xA,yA            :: areas of X and Y face of tracer cells
64    C maskUp           :: 2-D array for mask at W points
65    C uFld,vFld,wFld   :: Local copy of velocity field (3 components)
66    C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points
67    C rTrans           :: 2-D arrays of volume transports at W points
68    C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
69    C diffKh           :: horizontal diffusion coefficient
70    C diffK4           :: bi-harmonic diffusion coefficient
71    C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
72    C TracerN          :: tracer field @ time-step n (Note: only used
73    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
75    C                     or extrapolated fwd in time to n+1/2 if applying AB on Tr)
76    C tracerIdentity   :: tracer identifier (required for KPP,GM)
77    C advectionScheme  :: advection scheme to use (Horizontal plane)
78    C vertAdvecScheme  :: advection scheme to use (Vertical direction)
79    C calcAdvection    :: =False if Advec computed with multiDim scheme
80    C implicitAdvection:: =True if vertical Advec computed implicitly
81    C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)
82    C trUseGMRedi      :: true if this tracer uses GM-Redi
83    C trUseKPP         :: true if this tracer uses KPP
84    C myTime           :: current time
85    C myIter           :: iteration number
86    C myThid           :: thread number
87        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
88          INTEGER k,kUp,kDown,kM1
89        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91          _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92          _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93          _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94          _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _RL diffKh, diffK4        _RL diffKh, diffK4
100        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
102          _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
103          _RL deltaTLev(Nr)
104        INTEGER tracerIdentity        INTEGER tracerIdentity
105        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        INTEGER advectionScheme, vertAdvecScheme
106          LOGICAL calcAdvection
107          LOGICAL implicitAdvection, applyAB_onTracer
108          LOGICAL trUseGMRedi, trUseKPP
109          _RL     myTime
110          INTEGER myIter, myThid
111    
112    C !OUTPUT PARAMETERS: ==================================================
113    C gTracer          :: tendency array
114    C fVerT            :: 2 1/2D arrays for vertical advective flux
115        _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)
116        INTEGER myThid        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
117    
118  C     == Local variables ==  C !LOCAL VARIABLES: ====================================================
119  C     I, J, K - Loop counters  C i,j              :: loop indices
120    C df4              :: used for storing del^2 T for bi-harmonic term
121    C fZon             :: zonal flux
122    C fMer             :: meridional flux
123    C af               :: advective flux
124    C df               :: diffusive flux
125    C localT           :: local copy of tracer field
126    C locABT           :: local copy of (AB-extrapolated) tracer field
127    #ifdef ALLOW_DIAGNOSTICS
128          CHARACTER*8 diagName
129          CHARACTER*4 GAD_DIAG_SUFX, diagSufx
130          EXTERNAL    GAD_DIAG_SUFX
131    #endif
132        INTEGER i,j        INTEGER i,j
133        LOGICAL TOP_LAYER        _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134        _RL afFacT, dfFacT        _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141          _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
142          _RL advFac, rAdvFac
143    #ifdef GAD_SMOLARKIEWICZ_HACK
144          _RL outFlux, trac, fac, gTrFac
145    #endif
146    CEOP
147    
148  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
149  C--   only the kUp part of fverT is set in this subroutine  C--   only the kUp part of fverT is set in this subroutine
150  C--   the kDown is still required  C--   the kDown is still required
151        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
152  #endif  #endif
       DO j=1-OLy,sNy+OLy  
        DO i=1-OLx,sNx+OLx  
         fZon(i,j)      = 0.0  
         fMer(i,j)      = 0.0  
         fVerT(i,j,kUp) = 0.0  
        ENDDO  
       ENDDO  
153    
154        afFacT = 1. _d 0  #ifdef ALLOW_DIAGNOSTICS
155        dfFacT = 1. _d 0  C--   Set diagnostic suffix for the current tracer
156        TOP_LAYER = K .EQ. 1        IF ( useDiagnostics ) THEN
157            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
158          ENDIF
159    #endif
160    
161          advFac  = 0. _d 0
162          IF (calcAdvection) advFac = 1. _d 0
163          rAdvFac = rkSign*advFac
164          IF (implicitAdvection) rAdvFac = rkSign
165    
 C--   Make local copy of tracer array  
166        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
167         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
168          localT(i,j)=tracer(i,j,k,bi,bj)          fZon(i,j)      = 0. _d 0
169            fMer(i,j)      = 0. _d 0
170            fVerT(i,j,kUp) = 0. _d 0
171            df(i,j)        = 0. _d 0
172            df4(i,j)       = 0. _d 0
173         ENDDO         ENDDO
174        ENDDO        ENDDO
175    
176    C--   Make local copy of tracer array
177          IF ( applyAB_onTracer ) THEN
178            DO j=1-OLy,sNy+OLy
179             DO i=1-OLx,sNx+OLx
180              localT(i,j)=TracerN(i,j,k,bi,bj)
181              locABT(i,j)= TracAB(i,j,k,bi,bj)
182             ENDDO
183            ENDDO
184          ELSE
185            DO j=1-OLy,sNy+OLy
186             DO i=1-OLx,sNx+OLx
187              localT(i,j)= TracAB(i,j,k,bi,bj)
188              locABT(i,j)= TracAB(i,j,k,bi,bj)
189             ENDDO
190            ENDDO
191          ENDIF
192    
193    C--   Unless we have already calculated the advection terms we initialize
194    C     the tendency to zero.
195    C     <== now done earlier at the beginning of thermodynamics.
196    c     IF (calcAdvection) THEN
197    c      DO j=1-Oly,sNy+Oly
198    c       DO i=1-Olx,sNx+Olx
199    c        gTracer(i,j,k,bi,bj)=0. _d 0
200    c       ENDDO
201    c      ENDDO
202    c     ENDIF
203    
204  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
205        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
# Line 88  C--   Pre-calculate del^2 T if bi-harmon Line 211  C--   Pre-calculate del^2 T if bi-harmon
211  C--   Initialize net flux in X direction  C--   Initialize net flux in X direction
212        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
213         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
214          fZon(i,j) = 0.          fZon(i,j) = 0. _d 0
215         ENDDO         ENDDO
216        ENDDO        ENDDO
217    
218  C-    Advective flux in X  C-    Advective flux in X
219        IF (gad_advection_scheme.EQ.ENUM_CENTERED_2ND) THEN        IF (calcAdvection) THEN
220         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
221        ELSEIF (gad_advection_scheme.EQ.ENUM_FLUX_LIMIT) THEN             CALL GAD_C2_ADV_X( bi,bj,k, uTrans, locABT, af, myThid )
222         CALL GAD_FLUXLIMIT_ADV_X(          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
223       &      bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
224        ELSEIF (gad_advection_scheme.EQ.ENUM_UPWIND_3RD ) THEN             CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
225         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)       I              deltaTLev(k), uTrans, uFld, locABT,
226        ELSEIF (gad_advection_scheme.EQ.ENUM_CENTERED_4TH) THEN       O              af, myThid )
227         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)          ELSE
228        ELSE           DO j=1-OLy,sNy+OLy
229         STOP 'GAD_CALC_RHS: Bad gad_advection_scheme (X)'            DO i=1-OLx,sNx+OLx
230    #ifdef ALLOW_OBCS
231               maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
232    #else /* ALLOW_OBCS */
233               maskLocW(i,j) = _maskW(i,j,k,bi,bj)
234    #endif /* ALLOW_OBCS */
235              ENDDO
236             ENDDO
237             IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
238               CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
239         I              uTrans, uFld, maskLocW, locABT,
240         O              af, myThid )
241             ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
242               CALL GAD_U3_ADV_X( bi,bj,k, uTrans, maskLocW, locABT,
243         O              af, myThid )
244             ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
245               CALL GAD_C4_ADV_X( bi,bj,k, uTrans, maskLocW, locABT,
246         O              af, myThid )
247             ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
248               CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
249         I              uTrans, uFld, maskLocW, locABT,
250         O              af, myThid )
251             ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
252              IF ( inAdMode ) THEN
253    cph This block is to trick the adjoint:
254    cph IF inAdExact=.FALSE., we want to use DST3
255    cph with limiters in forward, but without limiters in reverse.
256               CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
257         I              uTrans, uFld, maskLocW, locABT,
258         O              af, myThid )
259              ELSE
260               CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
261         I              uTrans, uFld, maskLocW, locABT,
262         O              af, myThid )
263              ENDIF
264    #ifndef ALLOW_AUTODIFF_TAMC
265             ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
266               CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
267         I              uTrans, uFld, maskLocW, locABT,
268         O              af, myThid )
269    #endif
270             ELSE
271              STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
272             ENDIF
273            ENDIF
274            DO j=1-Oly,sNy+Oly
275             DO i=1-Olx,sNx+Olx
276              fZon(i,j) = fZon(i,j) + af(i,j)
277             ENDDO
278            ENDDO
279    #ifdef ALLOW_DIAGNOSTICS
280            IF ( useDiagnostics ) THEN
281              diagName = 'ADVx'//diagSufx
282              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
283            ENDIF
284    #endif
285        ENDIF        ENDIF
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         fZon(i,j) = fZon(i,j) + af(i,j)  
        ENDDO  
       ENDDO  
286    
287  C-    Diffusive flux in X  C-    Diffusive flux in X
288        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
289         CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)         CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
290        ELSE        ELSE
291         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
292          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
293           df(i,j) = 0.           df(i,j) = 0. _d 0
294          ENDDO          ENDDO
295         ENDDO         ENDDO
296        ENDIF        ENDIF
297    
298    C-    Add bi-harmonic diffusive flux in X
299          IF (diffK4 .NE. 0.) THEN
300           CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
301          ENDIF
302    
303  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
304  C-    GM/Redi flux in X  C-    GM/Redi flux in X
305        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
306  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
307          CALL GMREDI_XTRANSPORT(          IF ( applyAB_onTracer ) THEN
308       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_XTRANSPORT(
309       I     xA,Tracer,       I         iMin,iMax,jMin,jMax,bi,bj,k,
310       U     df,       I         xA,TracerN,tracerIdentity,
311       I     myThid)       U         df,
312         I         myThid)
313            ELSE
314              CALL GMREDI_XTRANSPORT(
315         I         iMin,iMax,jMin,jMax,bi,bj,k,
316         I         xA,TracAB, tracerIdentity,
317         U         df,
318         I         myThid)
319            ENDIF
320        ENDIF        ENDIF
321  #endif  #endif
322        DO j=jMin,jMax  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
323         DO i=iMin,iMax        DO j=1-Oly,sNy+Oly
324          fZon(i,j) = fZon(i,j) + df(i,j)         DO i=1-Olx,sNx+Olx
325            fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
326         ENDDO         ENDDO
327        ENDDO        ENDDO
328    
329  C-    Bi-harmonic duffusive flux in X  #ifdef ALLOW_DIAGNOSTICS
330        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
331         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)  C       excluding advective terms:
332         DO j=jMin,jMax        IF ( useDiagnostics .AND.
333          DO i=iMin,iMax       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
334           fZon(i,j) = fZon(i,j) + df(i,j)            diagName = 'DFxE'//diagSufx
335          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
336        ENDIF        ENDIF
337    #endif
338    
339  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
340        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
341         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
342          fMer(i,j) = 0.          fMer(i,j) = 0. _d 0
343         ENDDO         ENDDO
344        ENDDO        ENDDO
345    
346  C-    Advective flux in Y  C-    Advective flux in Y
347        IF (gad_advection_scheme.EQ.ENUM_CENTERED_2ND) THEN        IF (calcAdvection) THEN
348         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
349        ELSEIF (gad_advection_scheme.EQ.ENUM_FLUX_LIMIT) THEN             CALL GAD_C2_ADV_Y( bi,bj,k, vTrans, locABT, af, myThid )
350         CALL GAD_FLUXLIMIT_ADV_Y(          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
351       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
352        ELSEIF (gad_advection_scheme.EQ.ENUM_UPWIND_3RD ) THEN             CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
353         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)       I              deltaTLev(k), vTrans, vFld, locABT,
354        ELSEIF (gad_advection_scheme.EQ.ENUM_CENTERED_4TH) THEN       O              af, myThid )
355         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)          ELSE
356        ELSE           DO j=1-OLy,sNy+OLy
357         STOP 'GAD_CALC_RHS: Bad gad_advection_scheme (Y)'            DO i=1-OLx,sNx+OLx
358    #ifdef ALLOW_OBCS
359               maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
360    #else /* ALLOW_OBCS */
361               maskLocS(i,j) = _maskS(i,j,k,bi,bj)
362    #endif /* ALLOW_OBCS */
363              ENDDO
364             ENDDO
365             IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
366               CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
367         I              vTrans, vFld, maskLocS, locABT,
368         O              af, myThid )
369             ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
370               CALL GAD_U3_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
371         O              af, myThid )
372             ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
373               CALL GAD_C4_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
374         O              af, myThid )
375             ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
376               CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
377         I              vTrans, vFld, maskLocS, locABT,
378         O              af, myThid )
379             ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
380              IF ( inAdMode ) THEN
381    cph This block is to trick the adjoint:
382    cph IF inAdExact=.FALSE., we want to use DST3
383    cph with limiters in forward, but without limiters in reverse.
384               CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
385         I              vTrans, vFld, maskLocS, locABT,
386         O              af, myThid )
387              ELSE
388               CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
389         I              vTrans, vFld, maskLocS, locABT,
390         O              af, myThid )
391              ENDIF
392    #ifndef ALLOW_AUTODIFF_TAMC
393             ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
394               CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
395         I              vTrans, vFld, maskLocS, locABT,
396         O              af, myThid )
397    #endif
398             ELSE
399               STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
400             ENDIF
401            ENDIF
402            DO j=1-Oly,sNy+Oly
403             DO i=1-Olx,sNx+Olx
404              fMer(i,j) = fMer(i,j) + af(i,j)
405             ENDDO
406            ENDDO
407    #ifdef ALLOW_DIAGNOSTICS
408            IF ( useDiagnostics ) THEN
409              diagName = 'ADVy'//diagSufx
410              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
411            ENDIF
412    #endif
413        ENDIF        ENDIF
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fMer(i,j) = fMer(i,j) + af(i,j)  
        ENDDO  
       ENDDO  
414    
415  C-    Diffusive flux in Y  C-    Diffusive flux in Y
416        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
# Line 181  C-    Diffusive flux in Y Line 418  C-    Diffusive flux in Y
418        ELSE        ELSE
419         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
420          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
421           df(i,j) = 0.           df(i,j) = 0. _d 0
422          ENDDO          ENDDO
423         ENDDO         ENDDO
424        ENDIF        ENDIF
425    
426    C-    Add bi-harmonic flux in Y
427          IF (diffK4 .NE. 0.) THEN
428           CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
429          ENDIF
430    
431  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
432  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
433        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
434         CALL GMREDI_YTRANSPORT(  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
435  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*          IF ( applyAB_onTracer ) THEN
436       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_YTRANSPORT(
437       I     yA,Tracer,       I         iMin,iMax,jMin,jMax,bi,bj,k,
438       U     df,       I         yA,TracerN,tracerIdentity,
439       I     myThid)       U         df,
440         I         myThid)
441            ELSE
442              CALL GMREDI_YTRANSPORT(
443         I         iMin,iMax,jMin,jMax,bi,bj,k,
444         I         yA,TracAB, tracerIdentity,
445         U         df,
446         I         myThid)
447            ENDIF
448        ENDIF        ENDIF
449  #endif  #endif
450    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
451        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
452         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
453          fMer(i,j) = fMer(i,j) + df(i,j)          fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
454         ENDDO         ENDDO
455        ENDDO        ENDDO
456    
457  C-    Bi-harmonic flux in Y  #ifdef ALLOW_DIAGNOSTICS
458        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
459         CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)  C       excluding advective terms:
460         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
461          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
462           fMer(i,j) = fMer(i,j) + df(i,j)            diagName = 'DFyE'//diagSufx
463          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
464        ENDIF        ENDIF
465    #endif
466    
467  C--   Initialize net flux in R  C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         fVerT(i,j,kUp) = 0.  
        ENDDO  
       ENDDO  
   
468  C-    Advective flux in R  C-    Advective flux in R
469  C     Note: wVel needs to be masked  #ifdef ALLOW_AIM
470        IF (K.GE.2) THEN  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
471          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.
472         &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.k.LT.Nr)
473         &   ) THEN
474    #else
475          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN
476    #endif
477  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
478         IF (gad_advection_scheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
479          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
480         ELSEIF (gad_advection_scheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
481          CALL GAD_FLUXLIMIT_ADV_R(       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
482       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
483         ELSEIF (gad_advection_scheme.EQ.ENUM_UPWIND_3RD ) THEN       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
484          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)       O         af, myThid )
485         ELSEIF (gad_advection_scheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
486          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
487  c       CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
488         ELSE       O         af, myThid )
489          STOP 'GAD_CALC_RHS: Bad gad_advection_scheme (R)'          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
490         ENDIF            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
491  C-    Surface "correction" term at k>1 :          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
492         DO j=1-Oly,sNy+Oly            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
493          DO i=1-Olx,sNx+Olx          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
494           af(i,j) = af(i,j)            CALL GAD_DST3_ADV_R( bi,bj,k,
495       &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
496       &             rTrans(i,j)*Tracer(i,j,k,bi,bj)       O         af, myThid )
497          ENDDO          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
498         ENDDO  cph This block is to trick the adjoint:
499        ELSE  cph IF inAdExact=.FALSE., we want to use DST3
500  C-    Surface "correction" term at k=1 :  cph with limiters in forward, but without limiters in reverse.
501         DO j=1-Oly,sNy+Oly            IF ( inAdMode ) THEN
502          DO i=1-Olx,sNx+Olx             CALL GAD_DST3_ADV_R( bi,bj,k,
503           af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)       I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
504         O         af, myThid )
505              ELSE
506               CALL GAD_DST3FL_ADV_R( bi,bj,k,
507         I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
508         O         af, myThid )
509              ENDIF
510    #ifndef ALLOW_AUTODIFF_TAMC
511            ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
512               CALL GAD_OS7MP_ADV_R( bi,bj,k,
513         I         deltaTLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
514         O         af, myThid )
515    #endif
516            ELSE
517              STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
518            ENDIF
519    C-     add the advective flux to fVerT
520            DO j=1-Oly,sNy+Oly
521             DO i=1-Olx,sNx+Olx
522              fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
523             ENDDO
524          ENDDO          ENDDO
525         ENDDO  #ifdef ALLOW_DIAGNOSTICS
526            IF ( useDiagnostics ) THEN
527              diagName = 'ADVr'//diagSufx
528              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
529    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
530    C        does it only if k=1 (never the case here)
531              IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
532            ENDIF
533    #endif
534        ENDIF        ENDIF
 C-    add the advective flux to fVerT  
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         fVerT(i,j,kUp) = fVerT(i,j,kUp) + afFacT*af(i,j)  
        ENDDO  
       ENDDO  
535    
536  C-    Diffusive flux in R  C-    Diffusive flux in R
537  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
538  C           boundary condition.  C           boundary condition.
539        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
540         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
541          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
542           df(i,j) = 0.           df(i,j) = 0. _d 0
543          ENDDO          ENDDO
544         ENDDO         ENDDO
545        ELSE        ELSE
546         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         IF ( applyAB_onTracer ) THEN
547             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracerN,df,myThid)
548           ELSE
549             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracAB, df,myThid)
550           ENDIF
551        ENDIF        ENDIF
 c     DO j=jMin,jMax  
 c      DO i=iMin,iMax  
 c       fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)  
 c      ENDDO  
 c     ENDDO  
552    
553  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
554  C-    GM/Redi flux in R  C-    GM/Redi flux in R
555        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
556  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
557         CALL GMREDI_RTRANSPORT(          IF ( applyAB_onTracer ) THEN
558       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_RTRANSPORT(
559       I     maskUp,Tracer,       I         iMin,iMax,jMin,jMax,bi,bj,k,
560       U     df,       I         TracerN,tracerIdentity,
561       I     myThid)       U         df,
562  c      DO j=jMin,jMax       I         myThid)
563  c       DO i=iMin,iMax          ELSE
564  c        fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)            CALL GMREDI_RTRANSPORT(
565  c       ENDDO       I         iMin,iMax,jMin,jMax,bi,bj,k,
566  c      ENDDO       I         TracAB, tracerIdentity,
567         U         df,
568         I         myThid)
569            ENDIF
570        ENDIF        ENDIF
571  #endif  #endif
572    
573        DO j=jMin,jMax        DO j=1-Oly,sNy+Oly
574         DO i=iMin,iMax         DO i=1-Olx,sNx+Olx
575          fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)          fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
576         ENDDO         ENDDO
577        ENDDO        ENDDO
578    
579    #ifdef ALLOW_DIAGNOSTICS
580    C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
581    C       Explicit terms only & excluding advective terms:
582          IF ( useDiagnostics .AND.
583         &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN
584              diagName = 'DFrE'//diagSufx
585              CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
586          ENDIF
587    #endif
588    
589  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
590  C-    Add non local KPP transport term (ghat) to diffusive T flux.  C-    Set non local KPP transport term (ghat):
591        IF (useKPP) THEN        IF ( trUseKPP .AND. k.GE.2 ) THEN
592         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
593          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
594           df(i,j) = 0.           df(i,j) = 0. _d 0
595          ENDDO          ENDDO
596         ENDDO         ENDDO
597         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
 C *note* should update KPP_TRANSPORT_T to set df  *aja*  
598          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
599       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
600       I     KappaRT,       O           df,
601       U     df )       I           myTime, myIter, myThid )
602         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
603          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
604       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
605       I     KappaRT,       O           df,
606       U     df )       I           myTime, myIter, myThid )
607    #ifdef ALLOW_PTRACERS
608           ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
609            CALL KPP_TRANSPORT_PTR(
610         I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
611         I           tracerIdentity-GAD_TR1+1,
612         O           df,
613         I           myTime, myIter, myThid )
614    #endif
615         ELSE         ELSE
616          STOP 'GAD_CALC_RHS: Ooops'          WRITE(errorMessageUnit,*)
617         &    'tracer identity =', tracerIdentity, ' is not valid => STOP'
618            STOP 'ABNORMAL END: S/R GAD_CALC_RHS: invalid tracer identity'
619         ENDIF         ENDIF
620         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
621          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
622           fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp)
623         &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
624            ENDDO
625           ENDDO
626    #ifdef ALLOW_DIAGNOSTICS
627    C-    Diagnostics of Non-Local Tracer (vertical) flux
628           IF ( useDiagnostics ) THEN
629             diagName = 'KPPg'//diagSufx
630             CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
631    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
632    C        does it only if k=1 (never the case here)
633             IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
634           ENDIF
635    #endif
636          ENDIF
637    #endif /* ALLOW_KPP */
638    
639    #ifdef GAD_SMOLARKIEWICZ_HACK
640    coj   Hack to make redi (and everything else in this s/r) positive
641    coj   (see Smolarkiewicz MWR 1989 and Bott MWR 1989).
642    coj   Only works if 'down' is k+1 and k loop in thermodynamics is k=Nr,1,-1
643    coj
644    coj   Apply to all tracers except temperature
645          IF (tracerIdentity.NE.GAD_TEMPERATURE .AND.
646         &    tracerIdentity.NE.GAD_SALINITY) THEN
647           DO j=1-Oly,sNy+Oly-1
648            DO i=1-Olx,sNx+Olx-1
649    coj   Add outgoing fluxes
650             outFlux=deltaTLev(k)*
651         &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
652         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
653         &    *( MAX(0. _d 0,fZon(i+1,j)) + MAX(0. _d 0,-fZon(i,j))
654         &      +MAX(0. _d 0,fMer(i,j+1)) + MAX(0. _d 0,-fMer(i,j))
655         &      +MAX(0. _d 0,fVerT(i,j,kDown)*rkSign)
656         &      +MAX(0. _d 0,-fVerT(i,j,kUp)*rkSign)
657         &     )
658             IF ( applyAB_onTracer ) THEN
659               trac=TracerN(i,j,k,bi,bj)
660             ELSE
661               trac=TracAB(i,j,k,bi,bj)
662             ENDIF
663    coj   If they would reduce tracer by a fraction of more than
664    coj   SmolarkiewiczMaxFrac, scale them down
665             IF (outFlux.GT.0. _d 0 .AND.
666         &       outFlux.GT.SmolarkiewiczMaxFrac*trac) THEN
667    coj   If tracer is already negative, scale flux to zero
668               fac = MAX(0. _d 0,SmolarkiewiczMaxFrac*trac/outFlux)
669    
670               IF (fZon(i+1,j).GT.0. _d 0) fZon(i+1,j)=fac*fZon(i+1,j)
671               IF (-fZon(i,j) .GT.0. _d 0) fZon(i,j)  =fac*fZon(i,j)
672               IF (fMer(i,j+1).GT.0. _d 0) fMer(i,j+1)=fac*fMer(i,j+1)
673               IF (-fMer(i,j) .GT.0. _d 0) fMer(i,j)  =fac*fMer(i,j)
674               IF (-fVerT(i,j,kUp)*rkSign .GT.0. _d 0)
675         &       fVerT(i,j,kUp)=fac*fVerT(i,j,kUp)
676    
677               IF (k.LT.Nr .AND. fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN
678    coj   Down flux is special: it has already been applied in lower layer,
679    coj   so we have to readjust this.
680    coj   Note: for k+1, gTracer is now the updated tracer, not the tendency!
681    coj   thus it has an extra factor deltaTLev(k+1)
682                 gTrFac=deltaTLev(k+1)
683    coj   Other factors that have been applied to gTracer since the last call:
684    #ifdef NONLIN_FRSURF
685                 IF (nonlinFreeSurf.GT.0) THEN
686                  IF (select_rStar.GT.0) THEN
687    #ifndef DISABLE_RSTAR_CODE
688                    gTrFac = gTrFac/rStarExpC(i,j,bi,bj)
689    #endif /* DISABLE_RSTAR_CODE */
690                  ENDIF
691                 ENDIF
692    #endif /* NONLIN_FRSURF */
693    coj   Now: undo down flux, ...
694                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
695         &        +gTrFac
696         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
697         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
698         &         *recip_rhoFacC(k+1)
699         &         *( -fVerT(i,j,kDown)*rkSign )
700    coj   ... scale ...
701                 fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)
702    coj   ... and reapply
703                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
704         &        +gTrFac
705         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
706         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
707         &         *recip_rhoFacC(k+1)
708         &         *( fVerT(i,j,kDown)*rkSign )
709               ENDIF
710    
711             ENDIF
712          ENDDO          ENDDO
713         ENDDO         ENDDO
714        ENDIF        ENDIF
715  #endif  #endif
716    
717  C--   Divergence of fluxes  C--   Divergence of fluxes
718        DO j=jMin,jMax  C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
719         DO i=iMin,iMax        DO j=1-Oly,sNy+Oly-1
720          gTracer(i,j,k,bi,bj)=         DO i=1-Olx,sNx+Olx-1
721            gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
722       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
723       &    *recip_rA(i,j,bi,bj)       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
724       &    *(       &   *( (fZon(i+1,j)-fZon(i,j))
725       &    +( fZon(i+1,j)-fZon(i,j) )       &     +(fMer(i,j+1)-fMer(i,j))
726       &    +( fMer(i,j+1)-fMer(i,j) )       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
727       &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
728         &                   +(vTrans(i,j+1)-vTrans(i,j))*advFac
729         &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
730         &                  )
731       &    )       &    )
732         ENDDO         ENDDO
733        ENDDO        ENDDO
734    
735    #ifdef ALLOW_DEBUG
736          IF ( debugLevel .GE. debLevC
737         &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
738         &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
739         &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
740         &   .AND. useCubedSphereExchange ) THEN
741            CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
742         &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
743          ENDIF
744    #endif /* ALLOW_DEBUG */
745    
746        RETURN        RETURN
747        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.58
  ViewVC Help
Powered by ViewVC 1.1.22