/[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.17 by jmc, Sun Jan 26 21:08:36 2003 UTC revision 1.72 by jahn, Thu Sep 11 17:56:32 2014 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,maskUp,       I           xA, yA, maskUp, uFld, vFld, wFld,
13       I           diffKh, diffK4, KappaRT, Tracer,       I           uTrans, vTrans, rTrans, rTransKp1,
14       I           tracerIdentity, advectionScheme, calcAdvection,       I           diffKh, diffK4, KappaR, diffKr4, TracerN, TracAB,
15         I           deltaTLev, trIdentity,
16         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           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 27  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 40  C !USES: =============================== Line 45  C !USES: ===============================
45  #include "EEPARAMS.h"  #include "EEPARAMS.h"
46  #include "PARAMS.h"  #include "PARAMS.h"
47  #include "GRID.h"  #include "GRID.h"
 #include "DYNVARS.h"  
48  #include "SURFACE.h"  #include "SURFACE.h"
49  #include "GAD.h"  #include "GAD.h"
50    #ifdef ALLOW_AUTODIFF
51  #ifdef ALLOW_AUTODIFF_TAMC  # include "AUTODIFF_PARAMS.h"
52  #include "tamc.h"  #endif /* ALLOW_AUTODIFF */
 #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
56  C  iMin,iMax,jMin,jMax  :: loop range for called routines  C iMin, iMax       :: for called routines, to get valid output "gTracer"
57  C  kup                  :: index into 2 1/2D array, toggles between 1 and 2  C jMin, jMax       ::                      over this range of indices
58  C  kdown                :: index into 2 1/2D array, toggles between 2 and 1  C k                :: vertical index
59  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr  C kM1              :: =k-1 for k>1, =1 for k=1
60  C  xA,yA                :: areas of X and Y face of tracer cells  C kUp              :: index into 2 1/2D array, toggles between 1|2
61  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points  C kDown            :: index into 2 1/2D array, toggles between 2|1
62  C  maskUp               :: 2-D array for mask at W points  C xA, yA           :: areas of X and Y face of tracer cells
63  C  diffKh               :: horizontal diffusion coefficient  C maskUp           :: 2-D array for mask at W points
64  C  diffK4               :: bi-harmonic diffusion coefficient  C uFld, vFld, wFld :: Local copy of velocity field (3 components)
65  C  KappaRT              :: 3-D array for vertical diffusion coefficient  C uTrans, vTrans   :: 2-D arrays of volume transports at U,V points
66  C  Tracer               :: tracer field  C rTrans           :: 2-D arrays of volume transports at W points
67  C  tracerIdentity       :: identifier for the tracer (required only for KPP)  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
68  C  advectionScheme      :: advection scheme to use  C diffKh           :: horizontal diffusion coefficient
69  C  calcAdvection        :: =False if Advec terms computed with multiDim scheme  C diffK4           :: horizontal bi-harmonic diffusion coefficient
70  C  myThid               :: thread number  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
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 trIdentity       :: 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 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
86    C myIter           :: iteration number
87    C myThid           :: thread number
88        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
89        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
90        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92          _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93          _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94          _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95          _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _RS maskUp(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 KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL diffKr4(Nr)
103        INTEGER tracerIdentity        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104        INTEGER advectionScheme        _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
105          _RL deltaTLev(Nr)
106          INTEGER trIdentity
107          INTEGER advectionScheme, vertAdvecScheme
108        LOGICAL calcAdvection        LOGICAL calcAdvection
109        INTEGER myThid        LOGICAL implicitAdvection, applyAB_onTracer
110          LOGICAL trUseDiffKr4, trUseGMRedi, trUseKPP
111          _RL     myTime
112          INTEGER myIter, myThid
113    
114  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
115  C  gTracer              :: tendancy array  C gTracer          :: tendency array
116  C  fVerT                :: 2 1/2D arrays for vertical advective flux  C fZon             :: zonal flux
117        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  C fMer             :: meridional flux
118    C fVerT            :: 2 1/2D arrays for vertical advective flux
119          _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
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
133  C  fZon                 :: zonal flux  C af               :: advective flux
134  C  fmer                 :: meridional flux  C df               :: diffusive flux
135  C  af                   :: advective flux  C localT           :: local copy of tracer field
136  C  df                   :: diffusive flux  C locABT           :: local copy of (AB-extrapolated) tracer field
 C  localT               :: local copy of tracer field  
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)
145          _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    
161    #ifdef ALLOW_DIAGNOSTICS
162    C--   Set diagnostic suffix for the current tracer
163          IF ( useDiagnostics ) THEN
164            diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
165          ENDIF
166    #endif
167    
168          advFac  = 0. _d 0
169          IF (calcAdvection) advFac = 1. _d 0
170          rAdvFac = rkSign*advFac
171          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
175          fZon(i,j)      = 0. _d 0          fZon(i,j)      = 0. _d 0
# Line 118  C--   the kDown is still required Line 177  C--   the kDown is still required
177          fVerT(i,j,kUp) = 0. _d 0          fVerT(i,j,kUp) = 0. _d 0
178          df(i,j)        = 0. _d 0          df(i,j)        = 0. _d 0
179          df4(i,j)       = 0. _d 0          df4(i,j)       = 0. _d 0
         localT(i,j)    = 0. _d 0  
180         ENDDO         ENDDO
181        ENDDO        ENDDO
182    
183  C--   Make local copy of tracer array  C--   Make local copy of tracer array
184        DO j=1-OLy,sNy+OLy        IF ( applyAB_onTracer ) THEN
185         DO i=1-OLx,sNx+OLx          DO j=1-OLy,sNy+OLy
186          localT(i,j)=tracer(i,j,k,bi,bj)           DO i=1-OLx,sNx+OLx
187         ENDDO            localT(i,j)=TracerN(i,j,k)
188        ENDDO            locABT(i,j)= TracAB(i,j,k)
189             ENDDO
190  C--   Unless we have already calculated the advection terms we initialize          ENDDO
191  C     the tendency to zero.        ELSE
192        IF (calcAdvection) THEN          DO j=1-OLy,sNy+OLy
193         DO j=1-Oly,sNy+Oly           DO i=1-OLx,sNx+OLx
194          DO i=1-Olx,sNx+Olx            localT(i,j)=TracerN(i,j,k)
195           gTracer(i,j,k,bi,bj)=0. _d 0            locABT(i,j)=TracerN(i,j,k)
196             ENDDO
197          ENDDO          ENDDO
        ENDDO  
198        ENDIF        ENDIF
199    
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
# Line 147  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
213    
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,localT,af,myThid)             CALL GAD_C2_ADV_X( bi,bj,k, uTrans, locABT, af, myThid )
218        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
219         CALL GAD_FLUXLIMIT_ADV_X(       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
220       &      bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)             CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
221        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I              deltaTLev(k), uTrans, uFld, locABT,
222         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)       O              af, myThid )
223        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSE
224         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)           DO j=1-OLy,sNy+OLy
225        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            DO i=1-OLx,sNx+OLx
226         CALL GAD_DST3_ADV_X(  #ifdef ALLOW_OBCS
227       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)             maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
228        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN  #else /* ALLOW_OBCS */
229         CALL GAD_DST3FL_ADV_X(             maskLocW(i,j) = _maskW(i,j,k,bi,bj)
230       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)  #endif /* ALLOW_OBCS */
231        ELSE            ENDDO
232         STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'           ENDDO
233        ENDIF           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
234        DO j=1-Oly,sNy+Oly             CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
235         DO i=1-Olx,sNx+Olx       I              uTrans, uFld, maskLocW, locABT,
236          fZon(i,j) = fZon(i,j) + af(i,j)       O              af, myThid )
237         ENDDO           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
238        ENDDO             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:
248    cph If inAdExact=.FALSE., we want to use DST3
249    cph with limiters in forward, but without limiters in reverse.
250    #else /* ALLOW_AUTODIFF */
251             ELSEIF( advectionScheme.EQ.ENUM_DST3 ) THEN
252    #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
267              STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
268             ENDIF
269            ENDIF
270    #ifdef ALLOW_OBCS
271            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)
282             ENDDO
283            ENDDO
284    #ifdef ALLOW_DIAGNOSTICS
285            IF ( useDiagnostics ) THEN
286              diagName = 'ADVx'//diagSufx
287              CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
288            ENDIF
289    #endif
290        ENDIF        ENDIF
291    
292  C-    Diffusive flux in X  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
301        ENDIF        ENDIF
302    
303    C-    Add bi-harmonic diffusive flux in X
304          IF (diffK4 .NE. 0.) THEN
305           CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
306          ENDIF
307    
308  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
309  C-    GM/Redi flux in X  C-    GM/Redi flux in X
310        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
 C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  
311          CALL GMREDI_XTRANSPORT(          CALL GMREDI_XTRANSPORT(
312       I     iMin,iMax,jMin,jMax,bi,bj,K,       I         trIdentity, bi, bj, k, iMin, iMax+1, jMin, jMax,
313       I     xA,Tracer,tracerIdentity,       I         xA, TracerN,
314       U     df,       U         df,
315       I     myThid)       I         myThid )
316        ENDIF        ENDIF
317  #endif  #endif
318        DO j=1-Oly,sNy+Oly  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
319         DO i=1-Olx,sNx+Olx        DO j=1-OLy,sNy+OLy
320          fZon(i,j) = fZon(i,j) + df(i,j)         DO i=1-OLx,sNx+OLx
321            fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
322         ENDDO         ENDDO
323        ENDDO        ENDDO
324    
325  C-    Bi-harmonic duffusive flux in X  #ifdef ALLOW_DIAGNOSTICS
326        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
327         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)  C       excluding advective terms:
328         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
329          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
330           fZon(i,j) = fZon(i,j) + df(i,j)            diagName = 'DFxE'//diagSufx
331          ENDDO            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
332         ENDDO  #ifdef ALLOW_LAYERS
333              IF ( useLayers ) THEN
334               CALL LAYERS_FILL_DFX( df, trIdentity, k, 1, 2,bi,bj, myThid )
335              ENDIF
336    #endif /* ALLOW_LAYERS */
337        ENDIF        ENDIF
338    #endif
339    
340  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
341        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
342         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
343          fMer(i,j) = 0. _d 0          fMer(i,j) = 0. _d 0
344         ENDDO         ENDDO
345        ENDDO        ENDDO
346    
347  C-    Advective flux in Y  C-    Advective flux in Y
348        IF (calcAdvection) THEN        IF (calcAdvection) THEN
349        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
350         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)             CALL GAD_C2_ADV_Y( bi,bj,k, vTrans, locABT, af, myThid )
351        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
352         CALL GAD_FLUXLIMIT_ADV_Y(       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
353       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)             CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
354        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I              deltaTLev(k), vTrans, vFld, locABT,
355         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)       O              af, myThid )
356        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSE
357         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)           DO j=1-OLy,sNy+OLy
358        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            DO i=1-OLx,sNx+OLx
359         CALL GAD_DST3_ADV_Y(  #ifdef ALLOW_OBCS
360       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)             maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
361        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN  #else /* ALLOW_OBCS */
362         CALL GAD_DST3FL_ADV_Y(             maskLocS(i,j) = _maskS(i,j,k,bi,bj)
363       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)  #endif /* ALLOW_OBCS */
364        ELSE            ENDDO
365         STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'           ENDDO
366        ENDIF           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
367        DO j=1-Oly,sNy+Oly             CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
368         DO i=1-Olx,sNx+Olx       I              vTrans, vFld, maskLocS, locABT,
369          fMer(i,j) = fMer(i,j) + af(i,j)       O              af, myThid )
370         ENDDO           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
371        ENDDO             CALL GAD_U3_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
372         O              af, myThid )
373             ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
374               CALL GAD_C4_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
375         O              af, myThid )
376    #ifdef ALLOW_AUTODIFF
377             ELSEIF( advectionScheme.EQ.ENUM_DST3 .OR.
378         &          (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT .AND. inAdMode)
379         &         ) THEN
380    cph This block is to trick the adjoint:
381    cph If inAdExact=.FALSE., we want to use DST3
382    cph with limiters in forward, but without limiters in reverse.
383    #else /* ALLOW_AUTODIFF */
384             ELSEIF( advectionScheme.EQ.ENUM_DST3 ) THEN
385    #endif /* ALLOW_AUTODIFF */
386               CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
387         I              vTrans, vFld, maskLocS, locABT,
388         O              af, myThid )
389             ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
390               CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
391         I              vTrans, vFld, maskLocS, locABT,
392         O              af, myThid )
393    #ifndef ALLOW_AUTODIFF
394             ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
395               CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
396         I              vTrans, vFld, maskLocS, locABT,
397         O              af, myThid )
398    #endif
399             ELSE
400               STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
401             ENDIF
402            ENDIF
403    #ifdef ALLOW_OBCS
404            IF ( useOBCS ) THEN
405    C-      replace advective flux with 1st order upwind scheme estimate
406              CALL OBCS_U1_ADV_TRACER( .FALSE., trIdentity, bi, bj, k,
407         I                             maskS(1-OLx,1-OLy,k,bi,bj),
408         I                             vTrans, locABT,
409         U                             af, myThid )
410            ENDIF
411    #endif /* ALLOW_OBCS */
412            DO j=1-OLy,sNy+OLy
413             DO i=1-OLx,sNx+OLx
414              fMer(i,j) = fMer(i,j) + af(i,j)
415             ENDDO
416            ENDDO
417    #ifdef ALLOW_DIAGNOSTICS
418            IF ( useDiagnostics ) THEN
419              diagName = 'ADVy'//diagSufx
420              CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
421            ENDIF
422    #endif
423        ENDIF        ENDIF
424    
425  C-    Diffusive flux in Y  C-    Diffusive flux in Y
426        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
427         CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)         CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
428        ELSE        ELSE
429         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
430          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
431           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
432          ENDDO          ENDDO
433         ENDDO         ENDDO
434        ENDIF        ENDIF
435    
436    C-    Add bi-harmonic flux in Y
437          IF (diffK4 .NE. 0.) THEN
438           CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
439          ENDIF
440    
441  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
442  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
443        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
444  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*          CALL GMREDI_YTRANSPORT(
445         CALL GMREDI_YTRANSPORT(       I         trIdentity, bi, bj, k, iMin, iMax, jMin, jMax+1,
446       I     iMin,iMax,jMin,jMax,bi,bj,K,       I         yA, TracerN,
447       I     yA,Tracer,tracerIdentity,       U         df,
448       U     df,       I         myThid )
      I     myThid)  
449        ENDIF        ENDIF
450  #endif  #endif
451        DO j=1-Oly,sNy+Oly  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
452         DO i=1-Olx,sNx+Olx        DO j=1-OLy,sNy+OLy
453          fMer(i,j) = fMer(i,j) + df(i,j)         DO i=1-OLx,sNx+OLx
454            fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
455         ENDDO         ENDDO
456        ENDDO        ENDDO
457    
458  C-    Bi-harmonic flux in Y  #ifdef ALLOW_DIAGNOSTICS
459        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
460         CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)  C       excluding advective terms:
461         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
462          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
463           fMer(i,j) = fMer(i,j) + df(i,j)            diagName = 'DFyE'//diagSufx
464          ENDDO            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
465         ENDDO  #ifdef ALLOW_LAYERS
466        ENDIF            IF ( useLayers ) THEN
467               CALL LAYERS_FILL_DFY( df, trIdentity,k, 1, 2,bi,bj, myThid )
468  #ifdef NONLIN_FRSURF            ENDIF
469  C--   Compute vertical flux fVerT(kDown) at interface k+1 (between k & k+1):  #endif /* ALLOW_LAYERS */
       IF ( calcAdvection .AND. K.EQ.Nr .AND.  
      &     useRealFreshWaterFlux .AND.  
      &     buoyancyRelation .EQ. 'OCEANICP' ) THEN    
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          fVerT(i,j,kDown) = convertEmP2rUnit*PmEpR(i,j,bi,bj)  
      &     *rA(i,j,bi,bj)*maskC(i,j,k,bi,bj)*Tracer(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
470        ENDIF        ENDIF
471  #endif /* NONLIN_FRSURF */  #endif
472    
473  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):
474  C-    Advective flux in R  C-    Advective flux in R
475        IF (calcAdvection) THEN  #ifdef ALLOW_AIM
476  C     Note: wVel needs to be masked  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
477        IF (K.GE.2) THEN        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.
478  C-    Compute vertical advective flux in the interior:       &     (.NOT.useAIM .OR. trIdentity.NE.GAD_SALINITY .OR. k.LT.Nr)
479         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN       &   ) THEN
480          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)  #else
481         ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN
482          CALL GAD_FLUXLIMIT_ADV_R(  #endif
483       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)         IF ( applyAB_onTracer ) THEN
484         ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN  C-    Compute vertical advective flux in the interior using TracAB:
485          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
486         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN             CALL GAD_C2_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
487          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
488         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
489          CALL GAD_DST3_ADV_R(             CALL GAD_DST2U1_ADV_R( bi,bj,k,vertAdvecScheme,deltaTLev(k),
490       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)       I              rTrans, wFld, TracAB,
491         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       O              af, myThid )
492          CALL GAD_DST3FL_ADV_R(          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
493       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
494         I              rTrans, wFld, TracAB,
495         O              af, myThid )
496            ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
497               CALL GAD_U3_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
498            ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
499               CALL GAD_C4_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
500    #ifdef ALLOW_AUTODIFF
501            ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 .OR.
502         &         (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT .AND. inAdMode)
503         &        ) THEN
504    cph This block is to trick the adjoint:
505    cph If inAdExact=.FALSE., we want to use DST3
506    cph with limiters in forward, but without limiters in reverse.
507    #else /* ALLOW_AUTODIFF */
508            ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
509    #endif /* ALLOW_AUTODIFF */
510               CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
511         I              rTrans, wFld, TracAB,
512         O              af, myThid )
513            ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
514               CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
515         I              rTrans, wFld, TracAB,
516         O              af, myThid )
517    #ifndef ALLOW_AUTODIFF
518            ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
519               CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
520         I              rTrans, wFld, TracAB,
521         O              af, myThid )
522    #endif
523            ELSE
524              STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
525            ENDIF
526         ELSE         ELSE
527          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'  C-    Compute vertical advective flux in the interior using TracerN:
528            IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
529               CALL GAD_C2_ADV_R( bi,bj,k, rTrans, TracerN, af, myThid )
530            ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
531         &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
532               CALL GAD_DST2U1_ADV_R( bi,bj,k,vertAdvecScheme,deltaTLev(k),
533         I              rTrans, wFld, TracerN,
534         O              af, myThid )
535            ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
536               CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
537         I              rTrans, wFld, TracerN,
538         O              af, myThid )
539            ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
540               CALL GAD_U3_ADV_R( bi,bj,k, rTrans, TracerN, af, myThid )
541            ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
542               CALL GAD_C4_ADV_R( bi,bj,k, rTrans, TracerN, af, myThid )
543    #ifdef ALLOW_AUTODIFF
544            ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 .OR.
545         &         (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT .AND. inAdMode)
546         &        ) THEN
547    cph This block is to trick the adjoint:
548    cph If inAdExact=.FALSE., we want to use DST3
549    cph with limiters in forward, but without limiters in reverse.
550    #else /* ALLOW_AUTODIFF */
551            ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
552    #endif /* ALLOW_AUTODIFF */
553               CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
554         I              rTrans, wFld, TracerN,
555         O              af, myThid )
556            ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
557               CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
558         I              rTrans, wFld, TracerN,
559         O              af, myThid )
560    #ifndef ALLOW_AUTODIFF
561            ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
562               CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
563         I              rTrans, wFld, TracerN,
564         O              af, myThid )
565    #endif
566            ELSE
567              STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
568            ENDIF
569         ENDIF         ENDIF
570  C-    Surface "correction" term at k>1 :  C-     add the advective flux to fVerT
571         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
572          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
573           af(i,j) = af(i,j)            fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)*maskInC(i,j,bi,bj)
      &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  
      &             rTrans(i,j)*Tracer(i,j,k,bi,bj)  
574          ENDDO          ENDDO
        ENDDO  
       ELSE  
 C-    Surface "correction" term at k=1 :  
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
       ENDIF  
 C-    add the advective flux to fVerT  
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)  
575         ENDDO         ENDDO
576        ENDDO  #ifdef ALLOW_DIAGNOSTICS
577           IF ( useDiagnostics ) THEN
578              diagName = 'ADVr'//diagSufx
579              CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
580    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
581    C        does it only if k=1 (never the case here)
582              IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
583           ENDIF
584    #endif
585        ENDIF        ENDIF
586    
587  C-    Diffusive flux in R  C-    Diffusive flux in R
588  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
589  C           boundary condition.  C           boundary condition.
590        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
591         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
592          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
593           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
594          ENDDO          ENDDO
595         ENDDO         ENDDO
596        ELSE        ELSE
597         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)          CALL GAD_DIFF_R(bi,bj,k,KappaR,TracerN,df,myThid)
598          ENDIF
599    
600          IF ( trUseDiffKr4 ) THEN
601            CALL GAD_BIHARM_R( bi,bj,k, diffKr4, TracerN, df, myThid )
602        ENDIF        ENDIF
603    
604  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
605  C-    GM/Redi flux in R  C-    GM/Redi flux in R
606        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
607  C *note* should update GMREDI_RTRANSPORT to set df  *aja*          CALL GMREDI_RTRANSPORT(
608         CALL GMREDI_RTRANSPORT(       I         trIdentity, bi, bj, k, iMin, iMax, jMin, jMax,
609       I     iMin,iMax,jMin,jMax,bi,bj,K,       I         TracerN,
610       I     Tracer,tracerIdentity,       U         df,
611       U     df,       I         myThid )
      I     myThid)  
612        ENDIF        ENDIF
613  #endif  #endif
614    
615        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
616         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
617          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)
618         ENDDO         ENDDO
619        ENDDO        ENDDO
620    
621    #ifdef ALLOW_DIAGNOSTICS
622    C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
623    C       Explicit terms only & excluding advective terms:
624          IF ( useDiagnostics .AND.
625         &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN
626              diagName = 'DFrE'//diagSufx
627              CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
628    #ifdef ALLOW_LAYERS
629              IF ( useLayers ) THEN
630               CALL LAYERS_FILL_DFR( df, trIdentity, k, 1, 2,bi,bj, myThid )
631              ENDIF
632    #endif /* ALLOW_LAYERS */
633          ENDIF
634    #endif
635    
636  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
637  C-    Add non local KPP transport term (ghat) to diffusive T flux.  C-    Set non local KPP transport term (ghat):
638        IF (useKPP) THEN        IF ( trUseKPP .AND. k.GE.2 ) THEN
639         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
640          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
641           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
642          ENDDO          ENDDO
643         ENDDO         ENDDO
644         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (trIdentity.EQ.GAD_TEMPERATURE) THEN
 C *note* should update KPP_TRANSPORT_T to set df  *aja*  
645          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
646       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
647       I     KappaRT,       O           df,
648       U     df )       I           myTime, myIter, myThid )
649         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (trIdentity.EQ.GAD_SALINITY) THEN
650          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
651       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
652       I     KappaRT,       O           df,
653       U     df )       I           myTime, myIter, myThid )
654    #ifdef ALLOW_PTRACERS
655           ELSEIF (trIdentity .GE. GAD_TR1) THEN
656            CALL KPP_TRANSPORT_PTR(
657         I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
658         I           trIdentity-GAD_TR1+1,
659         O           df,
660         I           myTime, myIter, myThid )
661    #endif
662         ELSE         ELSE
663          STOP 'GAD_CALC_RHS: Ooops'          WRITE(errorMessageUnit,*)
664         &    'tracer identity =', trIdentity, ' is not valid => STOP'
665            STOP 'ABNORMAL END: S/R GAD_CALC_RHS: invalid tracer identity'
666           ENDIF
667           DO j=1-OLy,sNy+OLy
668            DO i=1-OLx,sNx+OLx
669             fVerT(i,j,kUp) = fVerT(i,j,kUp)
670         &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
671            ENDDO
672           ENDDO
673    #ifdef ALLOW_DIAGNOSTICS
674    C-    Diagnostics of Non-Local Tracer (vertical) flux
675           IF ( useDiagnostics ) THEN
676             diagName = 'KPPg'//diagSufx
677             CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
678    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
679    C        does it only if k=1 (never the case here)
680             IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
681    #ifdef ALLOW_LAYERS
682             IF ( useLayers ) THEN
683              CALL LAYERS_FILL_DFR( df, trIdentity, k, 1, 2,bi,bj, myThid )
684             ENDIF
685    #endif /* ALLOW_LAYERS */
686         ENDIF         ENDIF
687         DO j=1-Oly,sNy+Oly  #endif
688          DO i=1-Olx,sNx+Olx        ENDIF
689           fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)  #endif /* ALLOW_KPP */
690    
691    #ifdef GAD_SMOLARKIEWICZ_HACK
692    coj   Hack to make redi (and everything else in this s/r) positive
693    coj   (see Smolarkiewicz MWR 1989 and Bott MWR 1989).
694    coj   Only works if 'down' is k+1 and k loop in thermodynamics is k=Nr,1,-1
695    coj
696    coj   Apply to all tracers except temperature
697          IF ( trIdentity.NE.GAD_TEMPERATURE .AND.
698         &     trIdentity.NE.GAD_SALINITY ) THEN
699           DO j=1-OLy,sNy+OLy-1
700            DO i=1-OLx,sNx+OLx-1
701    coj   Add outgoing fluxes
702             outFlux=deltaTLev(k)*
703         &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
704         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
705         &    *( MAX(0. _d 0,fZon(i+1,j)) + MAX(0. _d 0,-fZon(i,j))
706         &      +MAX(0. _d 0,fMer(i,j+1)) + MAX(0. _d 0,-fMer(i,j))
707         &      +MAX(0. _d 0,fVerT(i,j,kDown)*rkSign)
708         &      +MAX(0. _d 0,-fVerT(i,j,kUp)*rkSign)
709         &     )
710             trac = localT(i,j)
711    coj   If they would reduce tracer by a fraction of more than
712    coj   SmolarkiewiczMaxFrac, scale them down
713             IF (outFlux.GT.0. _d 0 .AND.
714         &       outFlux.GT.SmolarkiewiczMaxFrac*trac) THEN
715    coj   If tracer is already negative, scale flux to zero
716               fac = MAX(0. _d 0,SmolarkiewiczMaxFrac*trac/outFlux)
717    
718               IF (fZon(i+1,j).GT.0. _d 0) fZon(i+1,j)=fac*fZon(i+1,j)
719               IF (-fZon(i,j) .GT.0. _d 0) fZon(i,j)  =fac*fZon(i,j)
720               IF (fMer(i,j+1).GT.0. _d 0) fMer(i,j+1)=fac*fMer(i,j+1)
721               IF (-fMer(i,j) .GT.0. _d 0) fMer(i,j)  =fac*fMer(i,j)
722               IF (-fVerT(i,j,kUp)*rkSign .GT.0. _d 0)
723         &       fVerT(i,j,kUp)=fac*fVerT(i,j,kUp)
724    
725               IF (k.LT.Nr .AND. fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN
726    coj   Down flux is special: it has already been applied in lower layer,
727    coj   so we have to readjust this.
728    coj   undo down flux, ...
729                 gTracer(i,j,k+1) = gTracer(i,j,k+1)
730         &        +_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
731         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
732         &         *recip_rhoFacC(k+1)
733         &         *( -fVerT(i,j,kDown)*rkSign )
734    coj   ... scale ...
735                 fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)
736    coj   ... and reapply
737                 gTracer(i,j,k+1) = gTracer(i,j,k+1)
738         &        +_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
739         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
740         &         *recip_rhoFacC(k+1)
741         &         *( fVerT(i,j,kDown)*rkSign )
742               ENDIF
743    
744             ENDIF
745          ENDDO          ENDDO
746         ENDDO         ENDDO
747        ENDIF        ENDIF
748  #endif  #endif
749    
750  C--   Divergence of fluxes  C--   Divergence of fluxes
751        DO j=1-Oly,sNy+Oly-1  C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
752         DO i=1-Olx,sNx+Olx-1  C     for Stevens OBC: keep only vertical diffusive contribution on boundaries
753          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)        DO j=1-OLy,sNy+OLy-1
754           DO i=1-OLx,sNx+OLx-1
755            gTracer(i,j,k) = gTracer(i,j,k)
756       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
757       &    *recip_rA(i,j,bi,bj)       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
758       &    *(       &   *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)
759       &    +( fZon(i+1,j)-fZon(i,j) )       &     +(fMer(i,j+1)-fMer(i,j))*maskInC(i,j,bi,bj)
760       &    +( fMer(i,j+1)-fMer(i,j) )       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
761       &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
762         &                   +(vTrans(i,j+1)-vTrans(i,j))*advFac
763         &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
764         &                  )*maskInC(i,j,bi,bj)
765       &    )       &    )
766         ENDDO         ENDDO
767        ENDDO        ENDDO
768    
769  #ifdef NONLIN_FRSURF  #ifdef ALLOW_DEBUG
770  C-- account for 3.D divergence of the flow in rStar coordinate:        IF ( debugLevel .GE. debLevC
771        IF (calcAdvection .AND. select_rStar.GT.0) THEN       &   .AND. trIdentity.EQ.GAD_TEMPERATURE
772         DO j=1-Oly,sNy+Oly-1       &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
773          DO i=1-Olx,sNx+Olx-1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
774           gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)       &   .AND. useCubedSphereExchange ) THEN
775       &     - (rStarExpC(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf          CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
776       &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)       &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
         ENDDO  
        ENDDO  
       ENDIF  
       IF (calcAdvection .AND. select_rStar.LT.0) THEN  
        DO j=1-Oly,sNy+Oly-1  
         DO i=1-Olx,sNx+Olx-1  
          gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)  
      &     - rStarDhCDt(i,j,bi,bj)  
      &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
777        ENDIF        ENDIF
778  #endif /* NONLIN_FRSURF */  #endif /* ALLOW_DEBUG */
         
779    
780        RETURN        RETURN
781        END        END
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.72
  ViewVC Help
Powered by ViewVC 1.1.22