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