/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_advection.F
ViewVC logotype

Diff of /MITgcm/pkg/generic_advdiff/gad_advection.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.4 by adcroft, Wed Sep 19 20:45:09 2001 UTC revision 1.33 by jmc, Thu Dec 16 22:28:43 2004 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
 CBOI  
 C !TITLE: pkg/generic\_advdiff  
 C !AUTHORS: adcroft@mit.edu  
 C !INTRODUCTION:  
 C \section{Generica Advection Diffusion Package}  
 C  
 C Package "generic\_advdiff" provides a common set of routines for calculating  
 C advective/diffusive fluxes for tracers (cell centered quantities on a C-grid).  
 C  
 C Many different advection schemes are available: the standard centered  
 C second order, centered fourth order and upwind biased third order schemes  
 C are known as linear methods and require some stable time-stepping method  
 C such as Adams-Bashforth. Alternatives such as flux-limited schemes are  
 C stable in the forward sense and are best combined with the multi-dimensional  
 C method provided in gad\_advection.  
 C  
 C There are two high-level routines:  
 C  \begin{itemize}  
 C  \item{GAD\_CALC\_RHS} calculates all fluxes at time level "n" and is used  
 C  for the standard linear schemes. This must be used in conjuction with  
 C  Adams-Bashforth time-stepping. Diffusive and parameterized fluxes are  
 C  always calculated here.  
 C  \item{GAD\_ADVECTION} calculates just the advective fluxes using the  
 C  non-linear schemes and can not be used in conjuction with Adams-Bashforth  
 C  time-stepping.  
 C  \end{itemize}  
 CEOI  
   
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    #undef MULTIDIM_OLD_VERSION
6    
7    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
8  CBOP  CBOP
9  C !ROUTINE: GAD_ADVECTION  C !ROUTINE: GAD_ADVECTION
10    
11  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
12        SUBROUTINE GAD_ADVECTION(bi,bj,advectionScheme,tracerIdentity,        SUBROUTINE GAD_ADVECTION(
13       U                         Tracer,Gtracer,       I     implicitAdvection, advectionScheme, vertAdvecScheme,
14       I                         myTime,myIter,myThid)       I     tracerIdentity,
15         I     uVel, vVel, wVel, tracer,
16         O     gTracer,
17         I     bi,bj, myTime,myIter,myThid)
18    
19  C !DESCRIPTION:  C !DESCRIPTION:
20  C Calculates the tendancy of a tracer due to advection.  C Calculates the tendancy of a tracer due to advection.
# Line 48  C Line 25  C
25  C The algorithm is as follows:  C The algorithm is as follows:
26  C \begin{itemize}  C \begin{itemize}
27  C \item{$\theta^{(n+1/3)} = \theta^{(n)}  C \item{$\theta^{(n+1/3)} = \theta^{(n)}
28  C      - \Delta t \partial_x (u\theta) + \theta \partial_x u$}  C      - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$}
29  C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)}  C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)}
30  C      - \Delta t \partial_y (v\theta) + \theta \partial_y v$}  C      - \Delta t \partial_y (v\theta^{(n+1/3)}) + \theta^{(n)} \partial_y v$}
31  C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)}  C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)}
32  C      - \Delta t \partial_r (w\theta) + \theta \partial_r w$}  C      - \Delta t \partial_r (w\theta^{(n+2/3)}) + \theta^{(n)} \partial_r w$}
33  C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}  C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}
34  C \end{itemize}  C \end{itemize}
35  C  C
# Line 63  C !USES: =============================== Line 40  C !USES: ===============================
40  #include "SIZE.h"  #include "SIZE.h"
41  #include "EEPARAMS.h"  #include "EEPARAMS.h"
42  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
43  #include "GRID.h"  #include "GRID.h"
44  #include "GAD.h"  #include "GAD.h"
45    #ifdef ALLOW_AUTODIFF_TAMC
46    # include "tamc.h"
47    # include "tamc_keys.h"
48    # ifdef ALLOW_PTRACERS
49    #  include "PTRACERS_SIZE.h"
50    # endif
51    #endif
52    #ifdef ALLOW_EXCH2
53    #include "W2_EXCH2_TOPOLOGY.h"
54    #include "W2_EXCH2_PARAMS.h"
55    #endif /* ALLOW_EXCH2 */
56    
57  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
58  C  bi,bj                :: tile indices  C  implicitAdvection :: implicit vertical advection (later on)
59  C  advectionScheme      :: advection scheme to use  C  advectionScheme   :: advection scheme to use (Horizontal plane)
60  C  tracerIdentity       :: identifier for the tracer (required only for OBCS)  C  vertAdvecScheme   :: advection scheme to use (vertical direction)
61  C  Tracer               :: tracer field  C  tracerIdentity    :: tracer identifier (required only for OBCS)
62  C  myTime               :: current time  C  uVel              :: velocity, zonal component
63  C  myIter               :: iteration number  C  vVel              :: velocity, meridional component
64  C  myThid               :: thread number  C  wVel              :: velocity, vertical component
65        INTEGER bi,bj  C  tracer            :: tracer field
66        INTEGER advectionScheme  C  bi,bj             :: tile indices
67    C  myTime            :: current time
68    C  myIter            :: iteration number
69    C  myThid            :: thread number
70          LOGICAL implicitAdvection
71          INTEGER advectionScheme, vertAdvecScheme
72        INTEGER tracerIdentity        INTEGER tracerIdentity
73        _RL Gtracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
74          _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
75          _RL wVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
76          _RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
77          INTEGER bi,bj
78        _RL myTime        _RL myTime
79        INTEGER myIter        INTEGER myIter
80        INTEGER myThid        INTEGER myThid
81    
82  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
83  C  gTracer              :: tendancy array  C  gTracer           :: tendancy array
84        _RL Tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
85    
86  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
87  C  maskUp               :: 2-D array for mask at W points  C  maskUp        :: 2-D array for mask at W points
88  C  iMin,iMax,jMin,jMax  :: loop range for called routines  C  maskLocW      :: 2-D array for mask at West points
89  C  i,j,k                :: loop indices  C  maskLocS      :: 2-D array for mask at South points
90  C  kup                  :: index into 2 1/2D array, toggles between 1 and 2  C  iMin,iMax,    :: loop range for called routines
91  C  kdown                :: index into 2 1/2D array, toggles between 2 and 1  C  jMin,jMax     :: loop range for called routines
92  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr  C [iMin,iMax]Upd :: loop range to update tracer field
93  C  xA,yA                :: areas of X and Y face of tracer cells  C [jMin,jMax]Upd :: loop range to update tracer field
94  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points  C  i,j,k         :: loop indices
95  C  af                   :: 2-D array for horizontal advective flux  C  kup           :: index into 2 1/2D array, toggles between 1 and 2
96  C  fVerT                :: 2 1/2D arrays for vertical advective flux  C  kdown         :: index into 2 1/2D array, toggles between 2 and 1
97  C  localTij             :: 2-D array used as temporary local copy of tracer fld  C  kp1           :: =k+1 for k<Nr, =Nr for k=Nr
98  C  localTijk            :: 3-D array used as temporary local copy of tracer fld  C  xA,yA         :: areas of X and Y face of tracer cells
99  C  kp1Msk               :: flag (0,1) to act as over-riding mask for W levels  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
100  C  calc_fluxes_X        :: logical to indicate to calculate fluxes in X dir  C  rTrans        :: 2-D arrays of volume transports at W points
101  C  calc_fluxes_Y        :: logical to indicate to calculate fluxes in Y dir  C  rTransKp1     :: vertical volume transport at interface k+1
102  C  nipass               :: number of passes to make in multi-dimensional method  C  af            :: 2-D array for horizontal advective flux
103  C  ipass                :: number of the current pass being made  C  afx           :: 2-D array for horizontal advective flux, x direction
104    C  afy           :: 2-D array for horizontal advective flux, y direction
105    C  fVerT         :: 2 1/2D arrays for vertical advective flux
106    C  localTij      :: 2-D array, temporary local copy of tracer fld
107    C  localTijk     :: 3-D array, temporary local copy of tracer fld
108    C  kp1Msk        :: flag (0,1) for over-riding mask for W levels
109    C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
110    C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
111    C  interiorOnly  :: only update the interior of myTile, but not the edges
112    C  overlapOnly   :: only update the edges of myTile, but not the interior
113    C  nipass        :: number of passes in multi-dimensional method
114    C  ipass         :: number of the current pass being made
115    C  myTile        :: variables used to determine which cube face
116    C  nCFace        :: owns a tile for cube grid runs using
117    C                :: multi-dim advection.
118    C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
119        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120          _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121          _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
123        INTEGER i,j,k,kup,kDown,kp1        INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
124          INTEGER i,j,k,kup,kDown
125        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132          _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133          _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
135        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136        _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
137        _RL kp1Msk        _RL kp1Msk
138        LOGICAL calc_fluxes_X,calc_fluxes_Y        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
139          LOGICAL interiorOnly, overlapOnly
140        INTEGER nipass,ipass        INTEGER nipass,ipass
141          INTEGER myTile, nCFace
142          LOGICAL N_edge, S_edge, E_edge, W_edge
143    #ifdef ALLOW_DIAGNOSTICS
144          INTEGER     kk
145          CHARACTER*8 diagName
146          CHARACTER*4 GAD_DIAG_SUFX, diagSufx
147          EXTERNAL    GAD_DIAG_SUFX
148    #endif
149  CEOP  CEOP
150    
151    #ifdef ALLOW_AUTODIFF_TAMC
152              act0 = tracerIdentity - 1
153              max0 = maxpass
154              act1 = bi - myBxLo(myThid)
155              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
156              act2 = bj - myByLo(myThid)
157              max2 = myByHi(myThid) - myByLo(myThid) + 1
158              act3 = myThid - 1
159              max3 = nTx*nTy
160              act4 = ikey_dynamics - 1
161              igadkey = (act0 + 1)
162         &                      + act1*max0
163         &                      + act2*max0*max1
164         &                      + act3*max0*max1*max2
165         &                      + act4*max0*max1*max2*max3
166              if (tracerIdentity.GT.maxpass) then
167                 print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
168                 STOP 'maxpass seems smaller than tracerIdentity'
169              endif
170    #endif /* ALLOW_AUTODIFF_TAMC */
171    
172    #ifdef ALLOW_DIAGNOSTICS
173    C--   Set diagnostic suffix for the current tracer
174          IF ( useDiagnostics ) THEN
175            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
176          ENDIF
177    #endif
178    
179  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
180  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
181  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
# Line 136  C     uninitialised but inert locations. Line 190  C     uninitialised but inert locations.
190          rTrans(i,j)  = 0. _d 0          rTrans(i,j)  = 0. _d 0
191          fVerT(i,j,1) = 0. _d 0          fVerT(i,j,1) = 0. _d 0
192          fVerT(i,j,2) = 0. _d 0          fVerT(i,j,2) = 0. _d 0
193            rTransKp1(i,j)= 0. _d 0
194         ENDDO         ENDDO
195        ENDDO        ENDDO
196    
197    C--   Set tile-specific parameters for horizontal fluxes
198          IF (useCubedSphereExchange) THEN
199           nipass=3
200    #ifdef ALLOW_AUTODIFF_TAMC
201           IF ( nipass.GT.maxcube ) STOP 'maxcube needs to be = 3'
202    #endif
203    #ifdef ALLOW_EXCH2
204           myTile = W2_myTileList(bi)
205           nCFace = exch2_myFace(myTile)
206           N_edge = exch2_isNedge(myTile).EQ.1
207           S_edge = exch2_isSedge(myTile).EQ.1
208           E_edge = exch2_isEedge(myTile).EQ.1
209           W_edge = exch2_isWedge(myTile).EQ.1
210    #else
211           nCFace = bi
212           N_edge = .TRUE.
213           S_edge = .TRUE.
214           E_edge = .TRUE.
215           W_edge = .TRUE.
216    #endif
217          ELSE
218           nipass=2
219           N_edge = .FALSE.
220           S_edge = .FALSE.
221           E_edge = .FALSE.
222           W_edge = .FALSE.
223          ENDIF
224    
225        iMin = 1-OLx        iMin = 1-OLx
226        iMax = sNx+OLx        iMax = sNx+OLx
227        jMin = 1-OLy        jMin = 1-OLy
# Line 146  C     uninitialised but inert locations. Line 229  C     uninitialised but inert locations.
229    
230  C--   Start of k loop for horizontal fluxes  C--   Start of k loop for horizontal fluxes
231        DO k=1,Nr        DO k=1,Nr
232    #ifdef ALLOW_AUTODIFF_TAMC
233             kkey = (igadkey-1)*Nr + k
234    CADJ STORE tracer(:,:,k,bi,bj) =
235    CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte
236    #endif /* ALLOW_AUTODIFF_TAMC */
237    
238  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
239        CALL CALC_COMMON_FACTORS (        CALL CALC_COMMON_FACTORS (
# Line 153  C--   Get temporary terms used by tenden Line 241  C--   Get temporary terms used by tenden
241       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
242       I         myThid)       I         myThid)
243    
244  C--   Make local copy of tracer array  #ifdef ALLOW_GMREDI
245    C--   Residual transp = Bolus transp + Eulerian transp
246          IF (useGMRedi)
247         &   CALL GMREDI_CALC_UVFLOW(
248         &            uTrans, vTrans, bi, bj, k, myThid)
249    #endif /* ALLOW_GMREDI */
250    
251    C--   Make local copy of tracer array and mask West & South
252        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
253         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
254          localTij(i,j)=tracer(i,j,k,bi,bj)           localTij(i,j)=tracer(i,j,k,bi,bj)
255             maskLocW(i,j)=maskW(i,j,k,bi,bj)
256             maskLocS(i,j)=maskS(i,j,k,bi,bj)
257         ENDDO         ENDDO
258        ENDDO        ENDDO
259    
260    #ifndef ALLOW_AUTODIFF_TAMC
261        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
262         nipass=3          withSigns = .FALSE.
263        ELSE          CALL FILL_CS_CORNER_UV_RS(
264         nipass=1       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )
265        ENDIF        ENDIF
266         nipass=1  #endif
267    
268  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
269    C--   For cube need one pass for each of red, green and blue axes.
270        DO ipass=1,nipass        DO ipass=1,nipass
271    #ifdef ALLOW_AUTODIFF_TAMC
272             passkey = ipass + (k-1)      *maxcube
273         &                   + (igadkey-1)*maxcube*Nr
274             IF (nipass .GT. maxpass) THEN
275              STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
276             ENDIF
277    #endif /* ALLOW_AUTODIFF_TAMC */
278    
279        IF (nipass.EQ.3) THEN        interiorOnly = .FALSE.
280         calc_fluxes_X=.FALSE.        overlapOnly  = .FALSE.
281         calc_fluxes_Y=.FALSE.        IF (useCubedSphereExchange) THEN
282         IF (ipass.EQ.1 .AND. (bi.EQ.1 .OR. bi.EQ.2) ) THEN  #ifdef MULTIDIM_OLD_VERSION
283          calc_fluxes_X=.TRUE.  C-    CubedSphere : pass 3 times, with full update of local tracer field
284         ELSEIF (ipass.EQ.1 .AND. (bi.EQ.4 .OR. bi.EQ.5) ) THEN         IF (ipass.EQ.1) THEN
285          calc_fluxes_Y=.TRUE.          calc_fluxes_X = nCFace.EQ.1 .OR. nCFace.EQ.2
286         ELSEIF (ipass.EQ.2 .AND. (bi.EQ.1 .OR. bi.EQ.6) ) THEN          calc_fluxes_Y = nCFace.EQ.4 .OR. nCFace.EQ.5
287          calc_fluxes_Y=.TRUE.         ELSEIF (ipass.EQ.2) THEN
288         ELSEIF (ipass.EQ.2 .AND. (bi.EQ.3 .OR. bi.EQ.4) ) THEN          calc_fluxes_X = nCFace.EQ.3 .OR. nCFace.EQ.4
289          calc_fluxes_X=.TRUE.          calc_fluxes_Y = nCFace.EQ.6 .OR. nCFace.EQ.1
290         ELSEIF (ipass.EQ.3 .AND. (bi.EQ.2 .OR. bi.EQ.3) ) THEN  #else /* MULTIDIM_OLD_VERSION */
291          calc_fluxes_Y=.TRUE.  C-    CubedSphere : pass 3 times, with partial update of local tracer field
292         ELSEIF (ipass.EQ.3 .AND. (bi.EQ.5 .OR. bi.EQ.6) ) THEN         IF (ipass.EQ.1) THEN
293          calc_fluxes_X=.TRUE.          overlapOnly  = MOD(nCFace,3).EQ.0
294            interiorOnly = MOD(nCFace,3).NE.0
295            calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
296            calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
297           ELSEIF (ipass.EQ.2) THEN
298            overlapOnly  = MOD(nCFace,3).EQ.2
299            calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
300            calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
301    #endif /* MULTIDIM_OLD_VERSION */
302           ELSE
303            calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
304            calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
305         ENDIF         ENDIF
306        ELSE        ELSE
307         calc_fluxes_X=.TRUE.  C-    not CubedSphere
308         calc_fluxes_Y=.TRUE.          calc_fluxes_X = MOD(ipass,2).EQ.1
309            calc_fluxes_Y = .NOT.calc_fluxes_X
310        ENDIF        ENDIF
311    
312    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
313  C--   X direction  C--   X direction
314        IF (calc_fluxes_X) THEN        IF (calc_fluxes_X) THEN
315    
316  C--   Internal exchange for calculations in X  C-     Do not compute fluxes if
317        IF (useCubedSphereExchange) THEN  C       a) needed in overlap only
318         DO j=1,Oly  C   and b) the overlap of myTile are not cube-face Edges
319          DO i=1,Olx         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
320           localTij( 1-i , 1-j )=localTij( 1-j ,    i    )  
321           localTij( 1-i ,sNy+j)=localTij( 1-j , sNy+1-i )  #ifndef ALLOW_AUTODIFF_TAMC
322           localTij(sNx+i, 1-j )=localTij(sNx+j,    i    )  C-     Internal exchange for calculations in X
323           localTij(sNx+i,sNy+j)=localTij(sNx+j, sNy+1-i )  #ifdef MULTIDIM_OLD_VERSION
324            IF ( useCubedSphereExchange ) THEN
325    #else
326            IF ( useCubedSphereExchange .AND.
327         &       ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
328    #endif
329             CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )
330            ENDIF
331    #endif
332    
333    C-     Advective flux in X
334            DO j=1-Oly,sNy+Oly
335             DO i=1-Olx,sNx+Olx
336              af(i,j) = 0.
337             ENDDO
338          ENDDO          ENDDO
        ENDDO  
       ENDIF  
339    
340  C-    Advective flux in X  #ifdef ALLOW_AUTODIFF_TAMC
341        DO j=1-Oly,sNy+Oly  #ifndef DISABLE_MULTIDIM_ADVECTION
342         DO i=1-Olx,sNx+Olx  CADJ STORE localTij(:,:)  =
343          af(i,j) = 0.  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
344         ENDDO  #endif
345        ENDDO  #endif /* ALLOW_AUTODIFF_TAMC */
346        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN  
347         CALL GAD_FLUXLIMIT_ADV_X(          IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
348       &      bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
349        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I                              uTrans, uVel, maskLocW, localTij,
350         CALL GAD_DST3_ADV_X(       O                              af, myThid )
351       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
352        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_DST3_ADV_X(      bi,bj,k, dTtracerLev(k),
353         CALL GAD_DST3FL_ADV_X(       I                              uTrans, uVel, maskLocW, localTij,
354       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)       O                              af, myThid )
355        ELSE          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
356         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            CALL GAD_DST3FL_ADV_X(    bi,bj,k, dTtracerLev(k),
357        ENDIF       I                              uTrans, uVel, maskLocW, localTij,
358        DO j=1-Oly,sNy+Oly       O                              af, myThid )
359         DO i=1-Olx,sNx+Olx-1          ELSE
360          localTij(i,j)=localTij(i,j)-deltaTtracer*           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
361       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)          ENDIF
362       &    *recip_rA(i,j,bi,bj)  
363       &    *( af(i+1,j)-af(i,j)  C-     Advective flux in X : done
364       &      -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))         ENDIF
365       &     )  
366         ENDDO  #ifndef ALLOW_AUTODIFF_TAMC
367        ENDDO  C-     Internal exchange for next calculations in Y
368           IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
369             CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )
370           ENDIF
371    #endif
372    
373    C-     Update the local tracer field where needed:
374    
375    C      update in overlap-Only
376           IF ( overlapOnly ) THEN
377            iMinUpd = 1-Olx+1
378            iMaxUpd = sNx+Olx-1
379    C- notes: these 2 lines below have no real effect (because recip_hFac=0
380    C         in corner region) but safer to keep them.
381            IF ( W_edge ) iMinUpd = 1
382            IF ( E_edge ) iMaxUpd = sNx
383    
384            IF ( S_edge ) THEN
385             DO j=1-Oly,0
386              DO i=iMinUpd,iMaxUpd
387               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
388         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
389         &       *recip_rA(i,j,bi,bj)
390         &       *( af(i+1,j)-af(i,j)
391         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
392         &        )
393              ENDDO
394             ENDDO
395            ENDIF
396            IF ( N_edge ) THEN
397             DO j=sNy+1,sNy+Oly
398              DO i=iMinUpd,iMaxUpd
399               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
400         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
401         &       *recip_rA(i,j,bi,bj)
402         &       *( af(i+1,j)-af(i,j)
403         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
404         &        )
405              ENDDO
406             ENDDO
407            ENDIF
408    
409           ELSE
410    C      do not only update the overlap
411            jMinUpd = 1-Oly
412            jMaxUpd = sNy+Oly
413            IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
414            IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
415            DO j=jMinUpd,jMaxUpd
416             DO i=1-Olx+1,sNx+Olx-1
417               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
418         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
419         &       *recip_rA(i,j,bi,bj)
420         &       *( af(i+1,j)-af(i,j)
421         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
422         &        )
423             ENDDO
424            ENDDO
425    C-      keep advective flux (for diagnostics)
426            DO j=1-Oly,sNy+Oly
427             DO i=1-Olx,sNx+Olx
428              afx(i,j) = af(i,j)
429             ENDDO
430            ENDDO
431    
432  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
433  C--   Apply open boundary conditions  C-     Apply open boundary conditions
434        IF (useOBCS) THEN          IF ( useOBCS ) THEN
435         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN           IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
436          CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )            CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
437         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN           ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
438          CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )            CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
439         END IF           ENDIF
440        END IF          ENDIF
441  #endif /* ALLOW_OBCS */  #endif /* ALLOW_OBCS */
442    
443    C-     end if/else update overlap-Only
444           ENDIF
445            
446  C--   End of X direction  C--   End of X direction
447        ENDIF        ENDIF
448    
449    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
450  C--   Y direction  C--   Y direction
451        IF (calc_fluxes_Y) THEN        IF (calc_fluxes_Y) THEN
452    
453  C--   Internal exchange for calculations in Y  C-     Do not compute fluxes if
454        IF (useCubedSphereExchange) THEN  C       a) needed in overlap only
455         DO j=1,Oly  C   and b) the overlap of myTile are not cube-face edges
456          DO i=1,Olx         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
457           localTij( 1-i , 1-j )=localTij(   j   , 1-i )  
458           localTij( 1-i ,sNy+j)=localTij(   j   ,sNy+i)  #ifndef ALLOW_AUTODIFF_TAMC
459           localTij(sNx+i, 1-j )=localTij(sNx+1-j, 1-i )  C-     Internal exchange for calculations in Y
460           localTij(sNx+i,sNy+j)=localTij(sNx+1-j,sNy+i)  #ifdef MULTIDIM_OLD_VERSION
461            IF ( useCubedSphereExchange ) THEN
462    #else
463            IF ( useCubedSphereExchange .AND.
464         &       ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
465    #endif
466             CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )
467            ENDIF
468    #endif
469    
470    C-     Advective flux in Y
471            DO j=1-Oly,sNy+Oly
472             DO i=1-Olx,sNx+Olx
473              af(i,j) = 0.
474             ENDDO
475          ENDDO          ENDDO
        ENDDO  
       ENDIF  
476    
477  C-    Advective flux in Y  #ifdef ALLOW_AUTODIFF_TAMC
478        DO j=1-Oly,sNy+Oly  #ifndef DISABLE_MULTIDIM_ADVECTION
479         DO i=1-Olx,sNx+Olx  CADJ STORE localTij(:,:)  =
480          af(i,j) = 0.  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
481         ENDDO  #endif
482        ENDDO  #endif /* ALLOW_AUTODIFF_TAMC */
483        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN  
484         CALL GAD_FLUXLIMIT_ADV_Y(          IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
485       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
486        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I                              vTrans, vVel, maskLocS, localTij,
487         CALL GAD_DST3_ADV_Y(       O                              af, myThid )
488       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
489        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_DST3_ADV_Y(      bi,bj,k, dTtracerLev(k),
490         CALL GAD_DST3FL_ADV_Y(       I                              vTrans, vVel, maskLocS, localTij,
491       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)       O                              af, myThid )
492        ELSE          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
493         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, dTtracerLev(k),
494        ENDIF       I                              vTrans, vVel, maskLocS, localTij,
495        DO j=1-Oly,sNy+Oly-1       O                              af, myThid )
496         DO i=1-Olx,sNx+Olx          ELSE
497          localTij(i,j)=localTij(i,j)-deltaTtracer*           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
498       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)          ENDIF
      &    *recip_rA(i,j,bi,bj)  
      &    *( af(i,j+1)-af(i,j)  
      &      -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))  
      &     )  
        ENDDO  
       ENDDO  
499    
500  #ifdef ALLOW_OBCS  C-     Advective flux in Y : done
501  C--   Apply open boundary conditions         ENDIF
       IF (useOBCS) THEN  
        IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN  
         CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )  
        ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN  
         CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )  
        END IF  
       END IF  
 #endif /* ALLOW_OBCS */  
502    
503  C--   End of Y direction  #ifndef ALLOW_AUTODIFF_TAMC
504        ENDIF  C-     Internal exchange for next calculations in X
505           IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
506             CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )
507           ENDIF
508    #endif
509    
510        DO j=1-Oly,sNy+Oly  C-     Update the local tracer field where needed:
        DO i=1-Olx,sNx+Olx  
         localTijk(i,j,k)=localTij(i,j)  
        ENDDO  
       ENDDO  
511    
512  C--   End of ipass loop  C      update in overlap-Only
513        ENDDO         IF ( overlapOnly ) THEN
514            jMinUpd = 1-Oly+1
515            jMaxUpd = sNy+Oly-1
516    C- notes: these 2 lines below have no real effect (because recip_hFac=0
517    C         in corner region) but safer to keep them.
518            IF ( S_edge ) jMinUpd = 1
519            IF ( N_edge ) jMaxUpd = sNy
520    
521            IF ( W_edge ) THEN
522             DO j=jMinUpd,jMaxUpd
523              DO i=1-Olx,0
524               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
525         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
526         &       *recip_rA(i,j,bi,bj)
527         &       *( af(i,j+1)-af(i,j)
528         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
529         &        )
530              ENDDO
531             ENDDO
532            ENDIF
533            IF ( E_edge ) THEN
534             DO j=jMinUpd,jMaxUpd
535              DO i=sNx+1,sNx+Olx
536               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
537         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
538         &       *recip_rA(i,j,bi,bj)
539         &       *( af(i,j+1)-af(i,j)
540         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
541         &        )
542              ENDDO
543             ENDDO
544            ENDIF
545    
546  C--   End of K loop for horizontal fluxes         ELSE
547        ENDDO  C      do not only update the overlap
548            iMinUpd = 1-Olx
549            iMaxUpd = sNx+Olx
550            IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
551            IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
552            DO j=1-Oly+1,sNy+Oly-1
553             DO i=iMinUpd,iMaxUpd
554               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
555         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
556         &       *recip_rA(i,j,bi,bj)
557         &       *( af(i,j+1)-af(i,j)
558         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
559         &        )
560             ENDDO
561            ENDDO
562    C-      keep advective flux (for diagnostics)
563            DO j=1-Oly,sNy+Oly
564             DO i=1-Olx,sNx+Olx
565              afy(i,j) = af(i,j)
566             ENDDO
567            ENDDO
568    
569  C--   Start of k loop for vertical flux  #ifdef ALLOW_OBCS
570        DO k=Nr,1,-1  C-     Apply open boundary conditions
571            IF (useOBCS) THEN
572             IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
573              CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
574             ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
575              CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
576             ENDIF
577            ENDIF
578    #endif /* ALLOW_OBCS */
579    
580  C--   kup    Cycles through 1,2 to point to w-layer above  C      end if/else update overlap-Only
581  C--   kDown  Cycles through 2,1 to point to w-layer below         ENDIF
       kup  = 1+MOD(k+1,2)  
       kDown= 1+MOD(k,2)  
582    
583  C--   Get temporary terms used by tendency routines  C--   End of Y direction
584        CALL CALC_COMMON_FACTORS (        ENDIF
      I         bi,bj,iMin,iMax,jMin,jMax,k,  
      O         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         myThid)  
585    
586  C-    Advective flux in R  C--   End of ipass loop
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         af(i,j) = 0.  
        ENDDO  
587        ENDDO        ENDDO
588    
589  C     Note: wVel needs to be masked        IF ( implicitAdvection ) THEN
590        IF (K.GE.2) THEN  C-    explicit advection is done ; store tendency in gTracer:
591  C-    Compute vertical advective flux in the interior:          DO j=1-Oly,sNy+Oly
592         IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           DO i=1-Olx,sNx+Olx
593          CALL GAD_FLUXLIMIT_ADV_R(            gTracer(i,j,k,bi,bj)=
594       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)
595         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN           ENDDO
         CALL GAD_DST3_ADV_R(  
      &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)  
        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN  
         CALL GAD_DST3FL_ADV_R(  
      &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)  
        ELSE  
         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'  
        ENDIF  
 C-    Surface "correction" term at k>1 :  
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          af(i,j) = af(i,j)  
      &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  
      &             rTrans(i,j)*localTijk(i,j,k)  
596          ENDDO          ENDDO
        ENDDO  
597        ELSE        ELSE
598  C-    Surface "correction" term at k=1 :  C-    horizontal advection done; store intermediate result in 3D array:
599         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
600          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
601           af(i,j) = rTrans(i,j)*localTijk(i,j,k)           localTijk(i,j,k)=localTij(i,j)
602          ENDDO          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) = af(i,j)  
603         ENDDO         ENDDO
604        ENDDO        ENDIF
605    
606  C--   Divergence of fluxes  #ifdef ALLOW_DIAGNOSTICS
607        kp1=min(Nr,k+1)          IF ( useDiagnostics ) THEN
608        kp1Msk=1.            kk = -k
609        if (k.EQ.Nr) kp1Msk=0.            diagName = 'ADVx'//diagSufx
610        DO j=1-Oly,sNy+Oly            CALL DIAGNOSTICS_FILL(afx,diagName, kk,1, 2,bi,bj, myThid)
611         DO i=1-Olx,sNx+Olx            diagName = 'ADVy'//diagSufx
612          localTij(i,j)=localTijk(i,j,k)-deltaTtracer*            CALL DIAGNOSTICS_FILL(afy,diagName, kk,1, 2,bi,bj, myThid)
613       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)          ENDIF
614       &    *recip_rA(i,j,bi,bj)  #endif
615       &    *( fVerT(i,j,kUp)-fVerT(i,j,kDown)  
616       &      -tracer(i,j,k,bi,bj)*rA(i,j,bi,bj)*  #ifdef ALLOW_DEBUG
617       &        (wVel(i,j,k,bi,bj)-kp1Msk*wVel(i,j,kp1,bi,bj))        IF ( debugLevel .GE. debLevB
618       &     )*rkFac       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
619          gTracer(i,j,k,bi,bj)=       &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
620       &   (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
621         ENDDO       &   .AND. useCubedSphereExchange ) THEN
622            CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION',
623         &             afx,afy, k, standardMessageUnit,bi,bj,myThid )
624          ENDIF
625    #endif /* ALLOW_DEBUG */
626    
627    C--   End of K loop for horizontal fluxes
628        ENDDO        ENDDO
629    
630    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
631    
632          IF ( .NOT.implicitAdvection ) THEN
633    C--   Start of k loop for vertical flux
634           DO k=Nr,1,-1
635    #ifdef ALLOW_AUTODIFF_TAMC
636             kkey = (igadkey-1)*Nr + k
637    #endif /* ALLOW_AUTODIFF_TAMC */
638    C--   kup    Cycles through 1,2 to point to w-layer above
639    C--   kDown  Cycles through 2,1 to point to w-layer below
640            kup  = 1+MOD(k+1,2)
641            kDown= 1+MOD(k,2)
642    c       kp1=min(Nr,k+1)
643            kp1Msk=1.
644            if (k.EQ.Nr) kp1Msk=0.
645    
646    C-- Compute Vertical transport
647    #ifdef ALLOW_AIM
648    C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
649            IF ( k.EQ.1 .OR.
650         &     (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
651         &              ) THEN
652    #else
653            IF ( k.EQ.1 ) THEN
654    #endif
655    
656    C- Surface interface :
657             DO j=1-Oly,sNy+Oly
658              DO i=1-Olx,sNx+Olx
659               rTransKp1(i,j) = kp1Msk*rTrans(i,j)
660               rTrans(i,j) = 0.
661               fVerT(i,j,kUp) = 0.
662              ENDDO
663             ENDDO
664    
665            ELSE
666    C- Interior interface :
667    
668             DO j=1-Oly,sNy+Oly
669              DO i=1-Olx,sNx+Olx
670               rTransKp1(i,j) = kp1Msk*rTrans(i,j)
671               rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
672         &                 *maskC(i,j,k-1,bi,bj)
673               fVerT(i,j,kUp) = 0.
674              ENDDO
675             ENDDO
676    
677    #ifdef ALLOW_GMREDI
678    C--   Residual transp = Bolus transp + Eulerian transp
679             IF (useGMRedi)
680         &   CALL GMREDI_CALC_WFLOW(
681         &                    rTrans, bi, bj, k, myThid)
682    #endif /* ALLOW_GMREDI */
683    
684    #ifdef ALLOW_AUTODIFF_TAMC
685    CADJ STORE localTijk(:,:,k)  
686    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
687    CADJ STORE rTrans(:,:)  
688    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
689    #endif /* ALLOW_AUTODIFF_TAMC */
690    
691    C-    Compute vertical advective flux in the interior:
692             IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
693    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
694               CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k),
695         I                               rTrans, wVel, localTijk,
696         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
697             ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
698               CALL GAD_DST3_ADV_R(      bi,bj,k, dTtracerLev(k),
699         I                               rTrans, wVel, localTijk,
700         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
701             ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
702               CALL GAD_DST3FL_ADV_R(    bi,bj,k, dTtracerLev(k),
703         I                               rTrans, wVel, localTijk,
704         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
705             ELSE
706              STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
707             ENDIF
708    
709    C- end Surface/Interior if bloc
710            ENDIF
711    
712    #ifdef ALLOW_AUTODIFF_TAMC
713    CADJ STORE rTrans(:,:)  
714    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
715    CADJ STORE rTranskp1(:,:)  
716    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
717    #endif /* ALLOW_AUTODIFF_TAMC */
718    
719    C--   Divergence of vertical fluxes
720            DO j=1-Oly,sNy+Oly
721             DO i=1-Olx,sNx+Olx
722              localTij(i,j)=localTijk(i,j,k)-dTtracerLev(k)*
723         &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
724         &     *recip_rA(i,j,bi,bj)
725         &     *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
726         &       -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
727         &      )*rkFac
728              gTracer(i,j,k,bi,bj)=
729         &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)
730             ENDDO
731            ENDDO
732    
733    #ifdef ALLOW_DIAGNOSTICS
734            IF ( useDiagnostics ) THEN
735              kk = -k
736              diagName = 'ADVr'//diagSufx
737              CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp),
738         &                           diagName, kk,1, 2,bi,bj, myThid)
739            ENDIF
740    #endif
741    
742  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
743        ENDDO         ENDDO
744    C--   end of if not.implicitAdvection block
745          ENDIF
746    
747        RETURN        RETURN
748        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22