/[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.1 by adcroft, Mon Sep 10 01:22:48 2001 UTC revision 1.68 by jmc, Tue Jun 21 16:13:52 2011 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6        SUBROUTINE GAD_ADVECTION(bi,bj,advectionScheme,tracerIdentity,  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7       U                         Tracer,Gtracer,  CBOP
8       I                         myTime,myIter,myThid)  C !ROUTINE: GAD_ADVECTION
9  C     /==========================================================\  
10  C     | SUBROUTINE GAD_ADVECTION                                 |  C !INTERFACE: ==========================================================
11  C     | o Solves the pure advection tracer equation.             |        SUBROUTINE GAD_ADVECTION(
12  C     |==========================================================|       I     implicitAdvection, advectionScheme, vertAdvecScheme,
13  C     \==========================================================/       I     tracerIdentity, deltaTLev,
14        IMPLICIT NONE       I     uVel, vVel, wVel, tracer,
15         O     gTracer,
16         I     bi,bj, myTime,myIter,myThid)
17    
18    C !DESCRIPTION:
19    C Calculates the tendency of a tracer due to advection.
20    C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
21    C and can only be used for the non-linear advection schemes such as the
22    C direct-space-time method and flux-limiters.
23    C
24    C The algorithm is as follows:
25    C \begin{itemize}
26    C \item{$\theta^{(n+1/3)} = \theta^{(n)}
27    C      - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$}
28    C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)}
29    C      - \Delta t \partial_y (v\theta^{(n+1/3)}) + \theta^{(n)} \partial_y v$}
30    C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)}
31    C      - \Delta t \partial_r (w\theta^{(n+2/3)}) + \theta^{(n)} \partial_r w$}
32    C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}
33    C \end{itemize}
34    C
35    C The tendency (output) is over-written by this routine.
36    
37  C     == Global variables ===  C !USES: ===============================================================
38          IMPLICIT NONE
39  #include "SIZE.h"  #include "SIZE.h"
40  #include "EEPARAMS.h"  #include "EEPARAMS.h"
41  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
42  #include "GRID.h"  #include "GRID.h"
43  #include "GAD.h"  #include "GAD.h"
44    #ifdef ALLOW_AUTODIFF_TAMC
45  C     == Routine arguments ==  # include "tamc.h"
46        INTEGER bi,bj  # include "tamc_keys.h"
47        INTEGER advectionScheme  # ifdef ALLOW_PTRACERS
48    #  include "PTRACERS_SIZE.h"
49    # endif
50    #endif
51    #ifdef ALLOW_EXCH2
52    #include "W2_EXCH2_SIZE.h"
53    #include "W2_EXCH2_TOPOLOGY.h"
54    #endif /* ALLOW_EXCH2 */
55    
56    C !INPUT PARAMETERS: ===================================================
57    C  implicitAdvection :: implicit vertical advection (later on)
58    C  advectionScheme   :: advection scheme to use (Horizontal plane)
59    C  vertAdvecScheme   :: advection scheme to use (vertical direction)
60    C  tracerIdentity    :: tracer identifier (required only for OBCS)
61    C  uVel              :: velocity, zonal component
62    C  vVel              :: velocity, meridional component
63    C  wVel              :: velocity, vertical component
64    C  tracer            :: tracer field
65    C  bi,bj             :: tile indices
66    C  myTime            :: current time
67    C  myIter            :: iteration number
68    C  myThid            :: thread number
69          LOGICAL implicitAdvection
70          INTEGER advectionScheme, vertAdvecScheme
71        INTEGER tracerIdentity        INTEGER tracerIdentity
72        _RL Tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL deltaTLev(Nr)
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     == Local variables  C !OUTPUT PARAMETERS: ==================================================
83        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C  gTracer           :: tendency array
84        INTEGER iMin,iMax,jMin,jMax        _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
85        INTEGER i,j,k,kup,kDown,kp1  
86    C !LOCAL VARIABLES: ====================================================
87    C  maskUp        :: 2-D array for mask at W points
88    C  maskLocW      :: 2-D array for mask at West points
89    C  maskLocS      :: 2-D array for mask at South points
90    C [iMin,iMax]Upd :: loop range to update tracer field
91    C [jMin,jMax]Upd :: loop range to update tracer field
92    C  i,j,k         :: loop indices
93    C  kUp           :: index into 2 1/2D array, toggles between 1 and 2
94    C  kDown         :: index into 2 1/2D array, toggles between 2 and 1
95    C  kp1           :: =k+1 for k<Nr, =Nr for k=Nr
96    C  xA,yA         :: areas of X and Y face of tracer cells
97    C  uFld,vFld     :: 2-D local copy of horizontal velocity, U,V components
98    C  wFld          :: 2-D local copy of vertical velocity
99    C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
100    C  rTrans        :: 2-D arrays of volume transports at W points
101    C  rTransKp1     :: vertical volume transport at interface k+1
102    C  af            :: 2-D array for horizontal advective flux
103    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  npass         :: 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    c     _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 iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
123          INTEGER i,j,k,kUp,kDown
124        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126          _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127          _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128          _RL wFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134          _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135          _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
137        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138        _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
139        _RL kp1Msk        _RL kp1Msk
140          LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
141          LOGICAL interiorOnly, overlapOnly
142          INTEGER npass, ipass
143          INTEGER nCFace
144          LOGICAL N_edge, S_edge, E_edge, W_edge
145    C     msgBuf     :: Informational/error message buffer
146          CHARACTER*(MAX_LEN_MBUF) msgBuf
147    #ifdef ALLOW_EXCH2
148          INTEGER myTile
149    #endif
150    #ifdef ALLOW_DIAGNOSTICS
151          CHARACTER*8 diagName
152          CHARACTER*4 diagSufx
153          LOGICAL     doDiagAdvX, doDiagAdvY, doDiagAdvR
154    C-    Functions:
155          CHARACTER*4 GAD_DIAG_SUFX
156          EXTERNAL    GAD_DIAG_SUFX
157          LOGICAL  DIAGNOSTICS_IS_ON
158          EXTERNAL DIAGNOSTICS_IS_ON
159    #endif
160    CEOP
161    
162    #ifdef ALLOW_AUTODIFF_TAMC
163              act0 = tracerIdentity
164              max0 = maxpass
165              act1 = bi - myBxLo(myThid)
166              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
167              act2 = bj - myByLo(myThid)
168              max2 = myByHi(myThid) - myByLo(myThid) + 1
169              act3 = myThid - 1
170              max3 = nTx*nTy
171              act4 = ikey_dynamics - 1
172              igadkey = act0
173         &                      + act1*max0
174         &                      + act2*max0*max1
175         &                      + act3*max0*max1*max2
176         &                      + act4*max0*max1*max2*max3
177              IF (tracerIdentity.GT.maxpass) THEN
178               WRITE(msgBuf,'(A,2I3)')
179         &      'GAD_ADVECTION: maxpass < tracerIdentity ',
180         &      maxpass, tracerIdentity
181               CALL PRINT_ERROR( msgBuf, myThid )
182               STOP 'ABNORMAL END: S/R GAD_ADVECTION'
183              ENDIF
184    #endif /* ALLOW_AUTODIFF_TAMC */
185    
186    #ifdef ALLOW_DIAGNOSTICS
187    C--   Set diagnostics flags and suffix for the current tracer
188          doDiagAdvX = .FALSE.
189          doDiagAdvY = .FALSE.
190          doDiagAdvR = .FALSE.
191          IF ( useDiagnostics ) THEN
192            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
193            diagName = 'ADVx'//diagSufx
194            doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )
195            diagName = 'ADVy'//diagSufx
196            doDiagAdvY = DIAGNOSTICS_IS_ON( diagName, myThid )
197            diagName = 'ADVr'//diagSufx
198            doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )
199          ENDIF
200    #endif
201    
202  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
203  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
# Line 60  C     uninitialised but inert locations. Line 213  C     uninitialised but inert locations.
213          rTrans(i,j)  = 0. _d 0          rTrans(i,j)  = 0. _d 0
214          fVerT(i,j,1) = 0. _d 0          fVerT(i,j,1) = 0. _d 0
215          fVerT(i,j,2) = 0. _d 0          fVerT(i,j,2) = 0. _d 0
216            rTransKp1(i,j)= 0. _d 0
217    #ifdef ALLOW_AUTODIFF_TAMC
218            localTij(i,j) = 0. _d 0
219            wfld(i,j)    = 0. _d 0
220    #endif
221         ENDDO         ENDDO
222        ENDDO        ENDDO
223    
224        iMin = 1-OLx  C--   Set tile-specific parameters for horizontal fluxes
225        iMax = sNx+OLx        IF (useCubedSphereExchange) THEN
226        jMin = 1-OLy         npass  = 3
227        jMax = sNy+OLy  #ifdef ALLOW_AUTODIFF_TAMC
228           IF ( npass.GT.maxcube ) STOP 'maxcube needs to be = 3'
229    #endif
230    #ifdef ALLOW_EXCH2
231           myTile = W2_myTileList(bi,bj)
232           nCFace = exch2_myFace(myTile)
233           N_edge = exch2_isNedge(myTile).EQ.1
234           S_edge = exch2_isSedge(myTile).EQ.1
235           E_edge = exch2_isEedge(myTile).EQ.1
236           W_edge = exch2_isWedge(myTile).EQ.1
237    #else
238           nCFace = bi
239           N_edge = .TRUE.
240           S_edge = .TRUE.
241           E_edge = .TRUE.
242           W_edge = .TRUE.
243    #endif
244          ELSE
245           npass  = 2
246           nCFace = 0
247           N_edge = .FALSE.
248           S_edge = .FALSE.
249           E_edge = .FALSE.
250           W_edge = .FALSE.
251          ENDIF
252    
253  C--   Start of k loop for horizontal fluxes  C--   Start of k loop for horizontal fluxes
254        DO k=1,Nr        DO k=1,Nr
255    #ifdef ALLOW_AUTODIFF_TAMC
256             kkey = (igadkey-1)*Nr + k
257    CADJ STORE tracer(:,:,k,bi,bj) =
258    CADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
259    #endif /* ALLOW_AUTODIFF_TAMC */
260    
261  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
262        CALL CALC_COMMON_FACTORS (        CALL CALC_COMMON_FACTORS (
263       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         uVel, vVel,
264       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         uFld, vFld, uTrans, vTrans, xA, yA,
265       I         myThid)       I         k,bi,bj, myThid )
266    
267    #ifdef ALLOW_GMREDI
268    C--   Residual transp = Bolus transp + Eulerian transp
269          IF (useGMRedi)
270         &   CALL GMREDI_CALC_UVFLOW(
271         U                  uFld, vFld, uTrans, vTrans,
272         I                  k, bi, bj, myThid )
273    #endif /* ALLOW_GMREDI */
274    
275  C--   Make local copy of tracer array  C--   Make local copy of tracer array and mask West & South
276        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
277         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
278          localTij(i,j)=tracer(i,j,k,bi,bj)           localTij(i,j)=tracer(i,j,k,bi,bj)
279    #ifdef ALLOW_OBCS
280             maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
281             maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
282    #else /* ALLOW_OBCS */
283             maskLocW(i,j) = _maskW(i,j,k,bi,bj)
284             maskLocS(i,j) = _maskS(i,j,k,bi,bj)
285    #endif /* ALLOW_OBCS */
286         ENDDO         ENDDO
287        ENDDO        ENDDO
288    
289  C-    Advective flux in X  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
290        DO j=1-Oly,sNy+Oly        IF (useCubedSphereExchange) THEN
291         DO i=1-Olx,sNx+Olx          withSigns = .FALSE.
292          af(i,j) = 0.          CALL FILL_CS_CORNER_UV_RS(
293         ENDDO       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )
294        ENDDO        ENDIF
295        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN  cph-exch2#endif
296         CALL GAD_FLUXLIMIT_ADV_X(  
297       &      bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)  C--   Multiple passes for different directions on different tiles
298        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN  C--   For cube need one pass for each of red, green and blue axes.
299         CALL GAD_DST3_ADV_X(        DO ipass=1,npass
300       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)  #ifdef ALLOW_AUTODIFF_TAMC
301        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN           passkey = ipass
302         CALL GAD_DST3FL_ADV_X(       &                   + (k-1)      *maxpass
303       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)       &                   + (igadkey-1)*maxpass*Nr
304             IF (npass .GT. maxpass) THEN
305              STOP 'GAD_ADVECTION: npass > maxcube. check tamc.h'
306             ENDIF
307    #endif /* ALLOW_AUTODIFF_TAMC */
308    
309          interiorOnly = .FALSE.
310          overlapOnly  = .FALSE.
311          IF (useCubedSphereExchange) THEN
312    C-    CubedSphere : pass 3 times, with partial update of local tracer field
313           IF (ipass.EQ.1) THEN
314            overlapOnly  = MOD(nCFace,3).EQ.0
315            interiorOnly = MOD(nCFace,3).NE.0
316            calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
317            calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
318           ELSEIF (ipass.EQ.2) THEN
319            overlapOnly  = MOD(nCFace,3).EQ.2
320            interiorOnly = MOD(nCFace,3).EQ.1
321            calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
322            calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
323           ELSE
324            interiorOnly = .TRUE.
325            calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
326            calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
327           ENDIF
328        ELSE        ELSE
329         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'  C-    not CubedSphere
330            calc_fluxes_X = MOD(ipass,2).EQ.1
331            calc_fluxes_Y = .NOT.calc_fluxes_X
332        ENDIF        ENDIF
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx-1  
         localTij(i,j)=localTij(i,j)-deltaTtracer*  
      &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  
      &    *recip_rA(i,j,bi,bj)  
      &    *( af(i+1,j)-af(i,j)  
      &      -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))  
      &     )  
        ENDDO  
       ENDDO  
333    
334  #ifdef ALLOW_OBCS  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
335  C--   Apply open boundary conditions  C--   X direction
336        IF (useOBCS) THEN  C-     Advective flux in X
337         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN          DO j=1-Oly,sNy+Oly
338          CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )           DO i=1-Olx,sNx+Olx
339         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN            af(i,j) = 0.
340          CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )           ENDDO
341         END IF          ENDDO
342        END IF  C
343  #endif /* ALLOW_OBCS */  #ifdef ALLOW_AUTODIFF_TAMC
344    # ifndef DISABLE_MULTIDIM_ADVECTION
345    CADJ STORE localTij(:,:)  =
346    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
347    CADJ STORE af(:,:)  =
348    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
349    # endif
350    #endif /* ALLOW_AUTODIFF_TAMC */
351    C
352          IF (calc_fluxes_X) THEN
353    
354    C-     Do not compute fluxes if
355    C       a) needed in overlap only
356    C   and b) the overlap of myTile are not cube-face Edges
357           IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
358    
359    C-     Internal exchange for calculations in X
360            IF ( overlapOnly ) THEN
361             CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
362         &                              localTij, bi,bj, myThid )
363            ENDIF
364    
365    #ifdef ALLOW_AUTODIFF_TAMC
366    # ifndef DISABLE_MULTIDIM_ADVECTION
367    CADJ STORE localTij(:,:)  =
368    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
369    # endif
370    #endif /* ALLOW_AUTODIFF_TAMC */
371    
372            IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
373         &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
374              CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
375         I                           deltaTLev(k),uTrans,uFld,localTij,
376         O                           af, myThid )
377            ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
378              CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
379         I                              uTrans, uFld, maskLocW, localTij,
380         O                              af, myThid )
381            ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
382              CALL GAD_DST3_ADV_X(      bi,bj,k, .TRUE., deltaTLev(k),
383         I                              uTrans, uFld, maskLocW, localTij,
384         O                              af, myThid )
385            ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
386              CALL GAD_DST3FL_ADV_X(    bi,bj,k, .TRUE., deltaTLev(k),
387         I                              uTrans, uFld, maskLocW, localTij,
388         O                              af, myThid )
389    #ifndef ALLOW_AUTODIFF_TAMC
390            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
391              CALL GAD_OS7MP_ADV_X(     bi,bj,k, .TRUE., deltaTLev(k),
392         I                              uTrans, uFld, maskLocW, localTij,
393         O                              af, myThid )
394    #endif
395            ELSE
396             STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
397            ENDIF
398    
399    C-     Internal exchange for next calculations in Y
400            IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
401             CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
402         &                              localTij, bi,bj, myThid )
403            ENDIF
404    
405  C-    Advective flux in Y  C-     Advective flux in X : done
406        DO j=1-Oly,sNy+Oly         ENDIF
407         DO i=1-Olx,sNx+Olx  
408          af(i,j) = 0.  C-     Update the local tracer field where needed:
409         ENDDO  C      use "maksInC" to prevent updating tracer field in OB regions
410        ENDDO  
411        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN  C      update in overlap-Only
412         CALL GAD_FLUXLIMIT_ADV_Y(         IF ( overlapOnly ) THEN
413       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)          iMinUpd = 1-Olx+1
414        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          iMaxUpd = sNx+Olx-1
415         CALL GAD_DST3_ADV_Y(  C- notes: these 2 lines below have no real effect (because recip_hFac=0
416       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)  C         in corner region) but safer to keep them.
417        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          IF ( W_edge ) iMinUpd = 1
418         CALL GAD_DST3FL_ADV_Y(          IF ( E_edge ) iMaxUpd = sNx
419       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)  
420        ELSE          IF ( S_edge ) THEN
421         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'           DO j=1-Oly,0
422              DO i=iMinUpd,iMaxUpd
423               localTij(i,j) = localTij(i,j)
424         &      -deltaTLev(k)*recip_rhoFacC(k)
425         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
426         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
427         &       *( af(i+1,j)-af(i,j)
428         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
429         &        )*maskInC(i,j,bi,bj)
430              ENDDO
431             ENDDO
432            ENDIF
433            IF ( N_edge ) THEN
434             DO j=sNy+1,sNy+Oly
435              DO i=iMinUpd,iMaxUpd
436               localTij(i,j) = localTij(i,j)
437         &      -deltaTLev(k)*recip_rhoFacC(k)
438         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
439         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
440         &       *( af(i+1,j)-af(i,j)
441         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
442         &        )*maskInC(i,j,bi,bj)
443              ENDDO
444             ENDDO
445            ENDIF
446    
447           ELSE
448    C      do not only update the overlap
449            jMinUpd = 1-Oly
450            jMaxUpd = sNy+Oly
451            IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
452            IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
453            DO j=jMinUpd,jMaxUpd
454             DO i=1-Olx+1,sNx+Olx-1
455               localTij(i,j) = localTij(i,j)
456         &      -deltaTLev(k)*recip_rhoFacC(k)
457         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
458         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
459         &       *( af(i+1,j)-af(i,j)
460         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
461         &        )*maskInC(i,j,bi,bj)
462             ENDDO
463            ENDDO
464    C-      keep advective flux (for diagnostics)
465            DO j=1-Oly,sNy+Oly
466             DO i=1-Olx,sNx+Olx
467              afx(i,j) = af(i,j)
468             ENDDO
469            ENDDO
470    
471    C-     end if/else update overlap-Only
472           ENDIF
473    
474    C--   End of X direction
475        ENDIF        ENDIF
       DO j=1-Oly,sNy+Oly-1  
        DO i=1-Olx,sNx+Olx  
         localTij(i,j)=localTij(i,j)-deltaTtracer*  
      &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  
      &    *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  
 #ifdef ALLOW_OBCS  
 C--   Apply open boundary conditions  
       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 */  
       DO j=1-Oly,sNy+Oly-1  
        DO i=1-Olx,sNx+Olx  
         localTijk(i,j,k)=localTij(i,j)  
        ENDDO  
       ENDDO  
476    
477    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
478    C--   Y direction
479    cph-test
480    C-     Advective flux in Y
481            DO j=1-Oly,sNy+Oly
482             DO i=1-Olx,sNx+Olx
483              af(i,j) = 0.
484             ENDDO
485            ENDDO
486    C
487    #ifdef ALLOW_AUTODIFF_TAMC
488    # ifndef DISABLE_MULTIDIM_ADVECTION
489    CADJ STORE localTij(:,:)  =
490    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
491    CADJ STORE af(:,:)  =
492    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
493    # endif
494    #endif /* ALLOW_AUTODIFF_TAMC */
495    C
496          IF (calc_fluxes_Y) THEN
497    
498    C-     Do not compute fluxes if
499    C       a) needed in overlap only
500    C   and b) the overlap of myTile are not cube-face edges
501           IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
502    
503    C-     Internal exchange for calculations in Y
504            IF ( overlapOnly ) THEN
505             CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
506         &                              localTij, bi,bj, myThid )
507            ENDIF
508    
509    C-     Advective flux in Y
510            DO j=1-Oly,sNy+Oly
511             DO i=1-Olx,sNx+Olx
512              af(i,j) = 0.
513             ENDDO
514            ENDDO
515    
516  C--   End of K loop for horizontal fluxes  #ifdef ALLOW_AUTODIFF_TAMC
517        ENDDO  #ifndef DISABLE_MULTIDIM_ADVECTION
518    CADJ STORE localTij(:,:)  =
519    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
520    #endif
521    #endif /* ALLOW_AUTODIFF_TAMC */
522    
523            IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
524         &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
525              CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
526         I                           deltaTLev(k),vTrans,vFld,localTij,
527         O                           af, myThid )
528            ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
529              CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
530         I                              vTrans, vFld, maskLocS, localTij,
531         O                              af, myThid )
532            ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
533              CALL GAD_DST3_ADV_Y(      bi,bj,k, .TRUE., deltaTLev(k),
534         I                              vTrans, vFld, maskLocS, localTij,
535         O                              af, myThid )
536            ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
537              CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .TRUE., deltaTLev(k),
538         I                              vTrans, vFld, maskLocS, localTij,
539         O                              af, myThid )
540    #ifndef ALLOW_AUTODIFF_TAMC
541            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
542              CALL GAD_OS7MP_ADV_Y(     bi,bj,k, .TRUE., deltaTLev(k),
543         I                              vTrans, vFld, maskLocS, localTij,
544         O                              af, myThid )
545    #endif
546            ELSE
547             STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
548            ENDIF
549    
550    C-     Internal exchange for next calculations in X
551            IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
552             CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
553         &                              localTij, bi,bj, myThid )
554            ENDIF
555    
556  C--   Start of k loop for vertical flux  C-     Advective flux in Y : done
557        DO k=Nr,1,-1         ENDIF
558    
559  C--   kup    Cycles through 1,2 to point to w-layer above  C-     Update the local tracer field where needed:
560  C--   kDown  Cycles through 2,1 to point to w-layer below  C      use "maksInC" to prevent updating tracer field in OB regions
       kup  = 1+MOD(k+1,2)  
       kDown= 1+MOD(k,2)  
561    
562  C--   Get temporary terms used by tendency routines  C      update in overlap-Only
563        CALL CALC_COMMON_FACTORS (         IF ( overlapOnly ) THEN
564       I         bi,bj,iMin,iMax,jMin,jMax,k,          jMinUpd = 1-Oly+1
565       O         xA,yA,uTrans,vTrans,rTrans,maskUp,          jMaxUpd = sNy+Oly-1
566       I         myThid)  C- notes: these 2 lines below have no real effect (because recip_hFac=0
567    C         in corner region) but safer to keep them.
568  C-    Advective flux in R          IF ( S_edge ) jMinUpd = 1
569        DO j=1-Oly,sNy+Oly          IF ( N_edge ) jMaxUpd = sNy
570         DO i=1-Olx,sNx+Olx  
571          af(i,j) = 0.          IF ( W_edge ) THEN
572         ENDDO           DO j=jMinUpd,jMaxUpd
573        ENDDO            DO i=1-Olx,0
574               localTij(i,j) = localTij(i,j)
575         &      -deltaTLev(k)*recip_rhoFacC(k)
576         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
577         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
578         &       *( af(i,j+1)-af(i,j)
579         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
580         &        )*maskInC(i,j,bi,bj)
581              ENDDO
582             ENDDO
583            ENDIF
584            IF ( E_edge ) THEN
585             DO j=jMinUpd,jMaxUpd
586              DO i=sNx+1,sNx+Olx
587               localTij(i,j) = localTij(i,j)
588         &      -deltaTLev(k)*recip_rhoFacC(k)
589         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
590         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
591         &       *( af(i,j+1)-af(i,j)
592         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
593         &        )*maskInC(i,j,bi,bj)
594              ENDDO
595             ENDDO
596            ENDIF
597    
 C     Note: wVel needs to be masked  
       IF (K.GE.2) THEN  
 C-    Compute vertical advective flux in the interior:  
        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN  
         CALL GAD_FLUXLIMIT_ADV_R(  
      &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)  
        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN  
 c       CALL GAD_DST3_ADV_R(  
 c    &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)  
         STOP 'GAD_ADVECTION: adv. scheme not avail. yet'  
        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN  
 c       CALL GAD_DST3FL_ADV_R(  
 c    &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)  
         STOP 'GAD_ADVECTION: adv. scheme not avail. yet'  
598         ELSE         ELSE
599          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'  C      do not only update the overlap
600            iMinUpd = 1-Olx
601            iMaxUpd = sNx+Olx
602            IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
603            IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
604            DO j=1-Oly+1,sNy+Oly-1
605             DO i=iMinUpd,iMaxUpd
606               localTij(i,j) = localTij(i,j)
607         &      -deltaTLev(k)*recip_rhoFacC(k)
608         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
609         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
610         &       *( af(i,j+1)-af(i,j)
611         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
612         &        )*maskInC(i,j,bi,bj)
613             ENDDO
614            ENDDO
615    C-      keep advective flux (for diagnostics)
616            DO j=1-Oly,sNy+Oly
617             DO i=1-Olx,sNx+Olx
618              afy(i,j) = af(i,j)
619             ENDDO
620            ENDDO
621    
622    C      end if/else update overlap-Only
623         ENDIF         ENDIF
624  C-    Surface "correction" term at k>1 :  
625         DO j=1-Oly,sNy+Oly  C--   End of Y direction
626          DO i=1-Olx,sNx+Olx        ENDIF
627           af(i,j) = af(i,j)  
628       &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  C--   End of ipass loop
629       &             rTrans(i,j)*localTijk(i,j,k)        ENDDO
630    
631          IF ( implicitAdvection ) THEN
632    C-    explicit advection is done ; store tendency in gTracer:
633            DO j=1-Oly,sNy+Oly
634             DO i=1-Olx,sNx+Olx
635              gTracer(i,j,k,bi,bj)=
636         &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
637             ENDDO
638          ENDDO          ENDDO
        ENDDO  
639        ELSE        ELSE
640  C-    Surface "correction" term at k=1 :  C-    horizontal advection done; store intermediate result in 3D array:
641         DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
642          DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
643           af(i,j) = rTrans(i,j)*localTijk(i,j,k)            localTijk(i,j,k)=localTij(i,j)
644             ENDDO
645          ENDDO          ENDDO
        ENDDO  
646        ENDIF        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)  
        ENDDO  
       ENDDO  
647    
648  C--   Divergence of fluxes  #ifdef ALLOW_DIAGNOSTICS
649        kp1=min(Nr,k+1)          IF ( doDiagAdvX ) THEN
650        kp1Msk=1.            diagName = 'ADVx'//diagSufx
651        if (k.EQ.Nr) kp1Msk=0.            CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)
652        DO j=1-Oly,sNy+Oly          ENDIF
653         DO i=1-Olx,sNx+Olx          IF ( doDiagAdvY ) THEN
654          localTij(i,j)=localTijk(i,j,k)-deltaTtracer*            diagName = 'ADVy'//diagSufx
655       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)            CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)
656       &    *recip_rA(i,j,bi,bj)          ENDIF
657       &    *( fVerT(i,j,kUp)-fVerT(i,j,kDown)  #endif
658       &      -tracer(i,j,k,bi,bj)*rA(i,j,bi,bj)*  
659       &        (wVel(i,j,k,bi,bj)-kp1Msk*wVel(i,j,kp1,bi,bj))  #ifdef ALLOW_DEBUG
660       &     )*rkFac        IF ( debugLevel .GE. debLevC
661          gTracer(i,j,k,bi,bj)=       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
662       &   (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer       &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
663         ENDDO       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
664         &   .AND. useCubedSphereExchange ) THEN
665            CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION',
666         &             afx,afy, k, standardMessageUnit,bi,bj,myThid )
667          ENDIF
668    #endif /* ALLOW_DEBUG */
669    
670    C--   End of K loop for horizontal fluxes
671        ENDDO        ENDDO
672    
673    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
674    
675          IF ( .NOT.implicitAdvection ) THEN
676    C--   Start of k loop for vertical flux
677           DO k=Nr,1,-1
678    #ifdef ALLOW_AUTODIFF_TAMC
679             kkey = (igadkey-1)*Nr + (Nr-k+1)
680    #endif /* ALLOW_AUTODIFF_TAMC */
681    C--   kUp    Cycles through 1,2 to point to w-layer above
682    C--   kDown  Cycles through 2,1 to point to w-layer below
683            kUp  = 1+MOD(k+1,2)
684            kDown= 1+MOD(k,2)
685    c       kp1=min(Nr,k+1)
686            kp1Msk=1.
687            if (k.EQ.Nr) kp1Msk=0.
688    
689    #ifdef ALLOW_AUTODIFF_TAMC
690    CADJ STORE rtrans(:,:)  =
691    CADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
692    cphCADJ STORE wfld(:,:)  =
693    cphCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
694    #endif
695    
696    C-- Compute Vertical transport
697    #ifdef ALLOW_AIM
698    C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
699            IF ( k.EQ.1 .OR.
700         &     (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
701         &              ) THEN
702    #else
703            IF ( k.EQ.1 ) THEN
704    #endif
705    
706    #ifdef ALLOW_AUTODIFF_TAMC
707    cphmultiCADJ STORE wfld(:,:)  =
708    cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
709    #endif /* ALLOW_AUTODIFF_TAMC */
710    
711    C- Surface interface :
712             DO j=1-Oly,sNy+Oly
713              DO i=1-Olx,sNx+Olx
714               rTransKp1(i,j) = kp1Msk*rTrans(i,j)
715               wFld(i,j)   = 0.
716               rTrans(i,j) = 0.
717               fVerT(i,j,kUp) = 0.
718              ENDDO
719             ENDDO
720    
721            ELSE
722    
723    #ifdef ALLOW_AUTODIFF_TAMC
724    cphmultiCADJ STORE wfld(:,:)  =
725    cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
726    #endif /* ALLOW_AUTODIFF_TAMC */
727    
728    C- Interior interface :
729             DO j=1-Oly,sNy+Oly
730              DO i=1-Olx,sNx+Olx
731               rTransKp1(i,j) = kp1Msk*rTrans(i,j)
732               wFld(i,j)   = wVel(i,j,k,bi,bj)
733               rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
734         &                 *deepFac2F(k)*rhoFacF(k)
735         &                 *maskC(i,j,k-1,bi,bj)
736               fVerT(i,j,kUp) = 0.
737              ENDDO
738             ENDDO
739    
740    #ifdef ALLOW_GMREDI
741    C--   Residual transp = Bolus transp + Eulerian transp
742             IF (useGMRedi)
743         &     CALL GMREDI_CALC_WFLOW(
744         U                 wFld, rTrans,
745         I                 k, bi, bj, myThid )
746    #endif /* ALLOW_GMREDI */
747    
748    #ifdef ALLOW_AUTODIFF_TAMC
749    cphmultiCADJ STORE localTijk(:,:,k)
750    cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
751    cphmultiCADJ STORE rTrans(:,:)
752    cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
753    #endif /* ALLOW_AUTODIFF_TAMC */
754    
755    C-    Compute vertical advective flux in the interior:
756             IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
757         &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
758               CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,
759         I                            deltaTLev(k),rTrans,wFld,localTijk,
760         O                            fVerT(1-Olx,1-Oly,kUp), myThid )
761             ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
762               CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
763         I                               rTrans, wFld, localTijk,
764         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
765             ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
766               CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTLev(k),
767         I                               rTrans, wFld, localTijk,
768         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
769             ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
770               CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTLev(k),
771         I                               rTrans, wFld, localTijk,
772         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
773    #ifndef ALLOW_AUTODIFF_TAMC
774             ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
775               CALL GAD_OS7MP_ADV_R(     bi,bj,k, deltaTLev(k),
776         I                               rTrans, wFld, localTijk,
777         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
778    #endif
779             ELSE
780              STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
781             ENDIF
782    
783    C- end Surface/Interior if bloc
784            ENDIF
785    
786    #ifdef ALLOW_AUTODIFF_TAMC
787    cphmultiCADJ STORE rTrans(:,:)
788    cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
789    cphmultiCADJ STORE rTranskp1(:,:)
790    cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
791    cph --- following storing of fVerT is critical for correct
792    cph --- gradient with multiDimAdvection
793    cph --- Without it, kDown component is not properly recomputed
794    cph --- This is a TAF bug (and no warning available)
795    CADJ STORE fVerT(:,:,:)
796    CADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
797    #endif /* ALLOW_AUTODIFF_TAMC */
798    
799    C--   Divergence of vertical fluxes
800            DO j=1-Oly,sNy+Oly
801             DO i=1-Olx,sNx+Olx
802              localTij(i,j) = localTijk(i,j,k)
803         &      -deltaTLev(k)*recip_rhoFacC(k)
804         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
805         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
806         &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)
807         &         -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))
808         &        )*rkSign
809              gTracer(i,j,k,bi,bj)=
810         &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
811             ENDDO
812            ENDDO
813    
814    #ifdef ALLOW_DIAGNOSTICS
815            IF ( doDiagAdvR ) THEN
816              diagName = 'ADVr'//diagSufx
817              CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp),
818         &                           diagName, k,1, 2,bi,bj, myThid)
819            ENDIF
820    #endif
821    
822  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
823        ENDDO         ENDDO
824    C--   end of if not.implicitAdvection block
825          ENDIF
826    
827        RETURN        RETURN
828        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.68
  ViewVC Help
Powered by ViewVC 1.1.22