/[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.58 by jmc, Wed Oct 22 00:28:47 2008 UTC revision 1.74 by jmc, Fri Apr 4 20:28:13 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  #undef MULTIDIM_OLD_VERSION  #ifdef ALLOW_AUTODIFF
6    # include "AUTODIFF_OPTIONS.h"
7    #endif
8    
9  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
10  CBOP  CBOP
# Line 11  C !ROUTINE: GAD_ADVECTION Line 13  C !ROUTINE: GAD_ADVECTION
13  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
14        SUBROUTINE GAD_ADVECTION(        SUBROUTINE GAD_ADVECTION(
15       I     implicitAdvection, advectionScheme, vertAdvecScheme,       I     implicitAdvection, advectionScheme, vertAdvecScheme,
16       I     tracerIdentity,       I     trIdentity, deltaTLev,
17       I     uVel, vVel, wVel, tracer,       I     uFld, vFld, wFld, tracer,
18       O     gTracer,       O     gTracer,
19       I     bi,bj, myTime,myIter,myThid)       I     bi,bj, myTime,myIter,myThid)
20    
# Line 42  C !USES: =============================== Line 44  C !USES: ===============================
44  #include "PARAMS.h"  #include "PARAMS.h"
45  #include "GRID.h"  #include "GRID.h"
46  #include "GAD.h"  #include "GAD.h"
47  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
48  # include "tamc.h"  # include "tamc.h"
49  # include "tamc_keys.h"  # include "tamc_keys.h"
50  # ifdef ALLOW_PTRACERS  # ifdef ALLOW_PTRACERS
# Line 50  C !USES: =============================== Line 52  C !USES: ===============================
52  # endif  # endif
53  #endif  #endif
54  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
55    #include "W2_EXCH2_SIZE.h"
56  #include "W2_EXCH2_TOPOLOGY.h"  #include "W2_EXCH2_TOPOLOGY.h"
 #include "W2_EXCH2_PARAMS.h"  
57  #endif /* ALLOW_EXCH2 */  #endif /* ALLOW_EXCH2 */
58    
59  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
60  C  implicitAdvection :: implicit vertical advection (later on)  C  implicitAdvection :: implicit vertical advection (later on)
61  C  advectionScheme   :: advection scheme to use (Horizontal plane)  C  advectionScheme   :: advection scheme to use (Horizontal plane)
62  C  vertAdvecScheme   :: advection scheme to use (vertical direction)  C  vertAdvecScheme   :: advection scheme to use (vertical direction)
63  C  tracerIdentity    :: tracer identifier (required only for OBCS)  C  trIdentity        :: tracer identifier
64  C  uVel              :: velocity, zonal component  C  uFld              :: Advection velocity field, zonal component
65  C  vVel              :: velocity, meridional component  C  vFld              :: Advection velocity field, meridional component
66  C  wVel              :: velocity, vertical component  C  wFld              :: Advection velocity field, vertical component
67  C  tracer            :: tracer field  C  tracer            :: tracer field
68  C  bi,bj             :: tile indices  C  bi,bj             :: tile indices
69  C  myTime            :: current time  C  myTime            :: current time
# Line 69  C  myIter            :: iteration number Line 71  C  myIter            :: iteration number
71  C  myThid            :: thread number  C  myThid            :: thread number
72        LOGICAL implicitAdvection        LOGICAL implicitAdvection
73        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
74        INTEGER tracerIdentity        INTEGER trIdentity
75        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL deltaTLev(Nr)
76        _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
77        _RL wVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
78        _RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _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
81        _RL myTime        _RL myTime
82        INTEGER myIter        INTEGER myIter
# Line 81  C  myThid            :: thread number Line 84  C  myThid            :: thread number
84    
85  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
86  C  gTracer           :: tendency array  C  gTracer           :: tendency array
87        _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
88    
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: ====================================================  C !LOCAL VARIABLES: ====================================================
98  C  maskUp        :: 2-D array for mask at W points  C  maskUp        :: 2-D array for mask at W points
# Line 92  C [jMin,jMax]Upd :: loop range to update Line 103  C [jMin,jMax]Upd :: loop range to update
103  C  i,j,k         :: loop indices  C  i,j,k         :: loop indices
104  C  kUp           :: index into 2 1/2D array, toggles between 1 and 2  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  C  kDown         :: index into 2 1/2D array, toggles between 2 and 1
 C  kp1           :: =k+1 for k<Nr, =Nr for k=Nr  
106  C  xA,yA         :: areas of X and Y face of tracer cells  C  xA,yA         :: areas of X and Y face of tracer cells
 C  uFld,vFld     :: 2-D local copy of horizontal velocity, U,V components  
 C  wFld          :: 2-D local copy of vertical velocity  
107  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points  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  C  rTrans        :: 2-D arrays of volume transports at W points
109  C  rTransKp1     :: vertical volume transport at interface k+1  C  rTransKp      :: vertical volume transport at interface k+1
110  C  af            :: 2-D array for horizontal advective flux  C  af            :: 2-D array for horizontal advective flux
111  C  afx           :: 2-D array for horizontal advective flux, x direction  C  afx           :: 2-D array for horizontal advective flux, x direction
112  C  afy           :: 2-D array for horizontal advective flux, y direction  C  afy           :: 2-D array for horizontal advective flux, y direction
113  C  fVerT         :: 2 1/2D arrays for vertical advective flux  C  fVerT         :: 2 1/2D arrays for vertical advective flux
114  C  localTij      :: 2-D array, temporary local copy of tracer fld  C  localTij      :: 2-D array, temporary local copy of tracer field
115  C  localTijk     :: 3-D array, temporary local copy of tracer fld  C  localT3d      :: 3-D array, temporary local copy of tracer field
116  C  kp1Msk        :: flag (0,1) for over-riding mask for W levels  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  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  C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
# Line 123  c     _RS maskUp  (1-OLx:sNx+OLx,1-OLy:s Line 131  c     _RS maskUp  (1-OLx:sNx+OLx,1-OLy:s
131        INTEGER i,j,k,kUp,kDown        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)
       _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL wFld    (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 rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _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)        _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140        _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _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        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
151        LOGICAL interiorOnly, overlapOnly        LOGICAL interiorOnly, overlapOnly
152        INTEGER npass, ipass        INTEGER npass, ipass
153        INTEGER nCFace        INTEGER nCFace
154        LOGICAL N_edge, S_edge, E_edge, W_edge        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  #ifdef ALLOW_EXCH2
160        INTEGER myTile        INTEGER myTile
161  #endif  #endif
# Line 149  c     _RS maskUp  (1-OLx:sNx+OLx,1-OLy:s Line 163  c     _RS maskUp  (1-OLx:sNx+OLx,1-OLy:s
163        CHARACTER*8 diagName        CHARACTER*8 diagName
164        CHARACTER*4 diagSufx        CHARACTER*4 diagSufx
165        LOGICAL     doDiagAdvX, doDiagAdvY, doDiagAdvR        LOGICAL     doDiagAdvX, doDiagAdvY, doDiagAdvR
166  C-    Functions:  #endif /* ALLOW_DIAGNOSTICS */
       CHARACTER*4 GAD_DIAG_SUFX  
       EXTERNAL    GAD_DIAG_SUFX  
       LOGICAL  DIAGNOSTICS_IS_ON  
       EXTERNAL DIAGNOSTICS_IS_ON  
 #endif  
167  CEOP  CEOP
168    
169  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
170            act0 = tracerIdentity            act0 = trIdentity
171            max0 = maxpass            max0 = maxpass
172            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
173            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
# Line 172  CEOP Line 181  CEOP
181       &                      + act2*max0*max1       &                      + act2*max0*max1
182       &                      + act3*max0*max1*max2       &                      + act3*max0*max1*max2
183       &                      + act4*max0*max1*max2*max3       &                      + act4*max0*max1*max2*max3
184            IF (tracerIdentity.GT.maxpass) THEN            IF (trIdentity.GT.maxpass) THEN
185               print *, 'ph-pass gad_advection ', maxpass, tracerIdentity             WRITE(msgBuf,'(A,2I3)')
186               STOP 'maxpass seems smaller than tracerIdentity'       &      'GAD_ADVECTION: maxpass < trIdentity ',
187         &      maxpass, trIdentity
188               CALL PRINT_ERROR( msgBuf, myThid )
189               STOP 'ABNORMAL END: S/R GAD_ADVECTION'
190            ENDIF            ENDIF
191  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
192    
# Line 184  C--   Set diagnostics flags and suffix f Line 196  C--   Set diagnostics flags and suffix f
196        doDiagAdvY = .FALSE.        doDiagAdvY = .FALSE.
197        doDiagAdvR = .FALSE.        doDiagAdvR = .FALSE.
198        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
199          diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )          diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
200          diagName = 'ADVx'//diagSufx          diagName = 'ADVx'//diagSufx
201          doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )          doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )
202          diagName = 'ADVy'//diagSufx          diagName = 'ADVy'//diagSufx
# Line 192  C--   Set diagnostics flags and suffix f Line 204  C--   Set diagnostics flags and suffix f
204          diagName = 'ADVr'//diagSufx          diagName = 'ADVr'//diagSufx
205          doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )          doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )
206        ENDIF        ENDIF
207  #endif  #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 201  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          rTransKp1(i,j)= 0. _d 0  #ifdef ALLOW_AUTODIFF
229  #ifdef ALLOW_AUTODIFF_TAMC  # ifdef GAD_MULTIDIM_COMPRESSIBLE
230            localVol(i,j) = 0. _d 0
231    # endif
232          localTij(i,j) = 0. _d 0          localTij(i,j) = 0. _d 0
233          wfld(i,j)    = 0. _d 0  #endif /* ALLOW_AUTODIFF */
 #endif  
234         ENDDO         ENDDO
235        ENDDO        ENDDO
236    
# Line 223  C--   Set tile-specific parameters for h Line 241  C--   Set tile-specific parameters for h
241         IF ( npass.GT.maxcube ) STOP 'maxcube needs to be = 3'         IF ( npass.GT.maxcube ) STOP 'maxcube needs to be = 3'
242  #endif  #endif
243  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
244         myTile = W2_myTileList(bi)         myTile = W2_myTileList(bi,bj)
245         nCFace = exch2_myFace(myTile)         nCFace = exch2_myFace(myTile)
246         N_edge = exch2_isNedge(myTile).EQ.1         N_edge = exch2_isNedge(myTile).EQ.1
247         S_edge = exch2_isSedge(myTile).EQ.1         S_edge = exch2_isSedge(myTile).EQ.1
# Line 250  C--   Start of k loop for horizontal flu Line 268  C--   Start of k loop for horizontal flu
268  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
269           kkey = (igadkey-1)*Nr + k           kkey = (igadkey-1)*Nr + k
270  CADJ STORE tracer(:,:,k,bi,bj) =  CADJ STORE tracer(:,:,k,bi,bj) =
271  CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
272  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
273    
274  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
275        CALL CALC_COMMON_FACTORS (        DO j=1-OLy,sNy+OLy
276       I         uVel, vVel,         DO i=1-OLx,sNx+OLx
277       O         uFld, vFld, uTrans, vTrans, xA, yA,           xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
278       I         k,bi,bj, myThid )       &           *drF(k)*_hFacW(i,j,k,bi,bj)
279             yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
280  #ifdef ALLOW_GMREDI       &           *drF(k)*_hFacS(i,j,k,bi,bj)
281  C--   Residual transp = Bolus transp + Eulerian transp         ENDDO
282        IF (useGMRedi)        ENDDO
283       &   CALL GMREDI_CALC_UVFLOW(  C--   Calculate "volume transports" through tracer cell faces.
284       U                  uFld, vFld, uTrans, vTrans,  C     anelastic: scaled by rhoFacC (~ mass transport)
285       I                  k, bi, bj, myThid )        DO j=1-OLy,sNy+OLy
286  #endif /* ALLOW_GMREDI */         DO i=1-OLx,sNx+OLx
287             uTrans(i,j) = uFld(i,j,k)*xA(i,j)*rhoFacC(k)
288             vTrans(i,j) = vFld(i,j,k)*yA(i,j)*rhoFacC(k)
289           ENDDO
290          ENDDO
291    
292  C--   Make local copy of tracer array and mask West & South  C--   Make local copy of tracer array and mask West & South
293        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
294         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
295           localTij(i,j)=tracer(i,j,k,bi,bj)           localTij(i,j) = tracer(i,j,k,bi,bj)
296           maskLocW(i,j)=_maskW(i,j,k,bi,bj)  #ifdef GAD_MULTIDIM_COMPRESSIBLE
297           maskLocS(i,j)=_maskS(i,j,k,bi,bj)           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
302             maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
303             maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
304    #else /* ALLOW_OBCS */
305             maskLocW(i,j) = _maskW(i,j,k,bi,bj)
306             maskLocS(i,j) = _maskS(i,j,k,bi,bj)
307    #endif /* ALLOW_OBCS */
308         ENDDO         ENDDO
309        ENDDO        ENDDO
310    
 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC  
311        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
312          withSigns = .FALSE.          withSigns = .FALSE.
313          CALL FILL_CS_CORNER_UV_RS(          CALL FILL_CS_CORNER_UV_RS(
314       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )
315        ENDIF        ENDIF
 cph-exch2#endif  
316    
317  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
318  C--   For cube need one pass for each of red, green and blue axes.  C--   For cube need one pass for each of red, green and blue axes.
# Line 299  C--   For cube need one pass for each of Line 329  C--   For cube need one pass for each of
329        interiorOnly = .FALSE.        interiorOnly = .FALSE.
330        overlapOnly  = .FALSE.        overlapOnly  = .FALSE.
331        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
 #ifdef MULTIDIM_OLD_VERSION  
 C-    CubedSphere : pass 3 times, with full update of local tracer field  
        IF (ipass.EQ.1) THEN  
         calc_fluxes_X = nCFace.EQ.1 .OR. nCFace.EQ.2  
         calc_fluxes_Y = nCFace.EQ.4 .OR. nCFace.EQ.5  
        ELSEIF (ipass.EQ.2) THEN  
         calc_fluxes_X = nCFace.EQ.3 .OR. nCFace.EQ.4  
         calc_fluxes_Y = nCFace.EQ.6 .OR. nCFace.EQ.1  
 #else /* MULTIDIM_OLD_VERSION */  
332  C-    CubedSphere : pass 3 times, with partial update of local tracer field  C-    CubedSphere : pass 3 times, with partial update of local tracer field
333         IF (ipass.EQ.1) THEN         IF (ipass.EQ.1) THEN
334          overlapOnly  = MOD(nCFace,3).EQ.0          overlapOnly  = MOD(nCFace,3).EQ.0
# Line 319  C-    CubedSphere : pass 3 times, with p Line 340  C-    CubedSphere : pass 3 times, with p
340          interiorOnly = MOD(nCFace,3).EQ.1          interiorOnly = MOD(nCFace,3).EQ.1
341          calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4          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          calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
 #endif /* MULTIDIM_OLD_VERSION */  
343         ELSE         ELSE
344          interiorOnly = .TRUE.          interiorOnly = .TRUE.
345          calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6          calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
# Line 333  C-    not CubedSphere Line 353  C-    not CubedSphere
353    
354  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
355  C--   X direction  C--   X direction
356  C-     Advective flux in X  
357          DO j=1-Oly,sNy+Oly  #ifdef ALLOW_AUTODIFF
358           DO i=1-Olx,sNx+Olx  C-     Always reset advective flux in X
359            DO j=1-OLy,sNy+OLy
360             DO i=1-OLx,sNx+OLx
361            af(i,j) = 0.            af(i,j) = 0.
362           ENDDO           ENDDO
363          ENDDO          ENDDO
 C  
 #ifdef ALLOW_AUTODIFF_TAMC  
364  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
365  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
366  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
367  CADJ STORE af(:,:)  =  CADJ STORE af(:,:)  =
368  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
369  # endif  # endif
370  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
371  C  
372        IF (calc_fluxes_X) THEN        IF (calc_fluxes_X) THEN
373    
374  C-     Do not compute fluxes if  C-     Do not compute fluxes if
# Line 357  C   and b) the overlap of myTile are not Line 377  C   and b) the overlap of myTile are not
377         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
378    
379  C-     Internal exchange for calculations in X  C-     Internal exchange for calculations in X
 #ifdef MULTIDIM_OLD_VERSION  
         IF ( useCubedSphereExchange ) THEN  
 #else  
380          IF ( overlapOnly ) THEN          IF ( overlapOnly ) THEN
 #endif  
381           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
382       &                              localTij, bi,bj, myThid )       &                              localTij, bi,bj, myThid )
383          ENDIF          ENDIF
384    
385  #ifdef ALLOW_AUTODIFF_TAMC  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  # ifndef DISABLE_MULTIDIM_ADVECTION
394  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
395  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
396  # endif  # endif
397  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
398    
399          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
400       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
401            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,            CALL GAD_DST2U1_ADV_X(    bi,bj,k, advectionScheme, .TRUE.,
402       I                           dTtracerLev(k),uTrans,uFld,localTij,       I             deltaTLev(k),uTrans,uFld(1-OLx,1-OLy,k), localTij,
403       O                           af, myThid )       O             af, myThid )
404          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
405            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
406       I                              uTrans, uFld, maskLocW, localTij,       I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
407       O                              af, myThid )       O             af, myThid )
408          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
409            CALL GAD_DST3_ADV_X(      bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_DST3_ADV_X(      bi,bj,k, .TRUE., deltaTLev(k),
410       I                              uTrans, uFld, maskLocW, localTij,       I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
411       O                              af, myThid )       O             af, myThid )
412          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
413            CALL GAD_DST3FL_ADV_X(    bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_DST3FL_ADV_X(    bi,bj,k, .TRUE., deltaTLev(k),
414       I                              uTrans, uFld, maskLocW, localTij,       I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
415       O                              af, myThid )       O             af, myThid )
416  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF
417          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
418            CALL GAD_OS7MP_ADV_X(     bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_OS7MP_ADV_X(     bi,bj,k, .TRUE., deltaTLev(k),
419       I                              uTrans, uFld, maskLocW, localTij,       I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
420       O                              af, myThid )       O             af, myThid )
421  #endif  #endif
422          ELSE          ELSE
423           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
424          ENDIF          ENDIF
425    
426    #ifdef ALLOW_OBCS
427            IF ( useOBCS ) THEN
428    C-      replace advective flux with 1st order upwind scheme estimate
429              CALL OBCS_U1_ADV_TRACER( .TRUE., trIdentity, bi, bj, k,
430         I                             maskLocW, uTrans, localTij,
431         U                             af, myThid )
432            ENDIF
433    #endif /* ALLOW_OBCS */
434    
435  C-     Internal exchange for next calculations in Y  C-     Internal exchange for next calculations in Y
436          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
437           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
# Line 410  C-     Advective flux in X : done Line 442  C-     Advective flux in X : done
442         ENDIF         ENDIF
443    
444  C-     Update the local tracer field where needed:  C-     Update the local tracer field where needed:
445    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  C      update in overlap-Only
456         IF ( overlapOnly ) THEN         IF ( overlapOnly ) THEN
457          iMinUpd = 1-Olx+1          iMinUpd = 1-OLx+1
458          iMaxUpd = sNx+Olx-1          iMaxUpd = sNx+OLx-1
459  C- notes: these 2 lines below have no real effect (because recip_hFac=0  C- notes: these 2 lines below have no real effect (because recip_hFac=0
460  C         in corner region) but safer to keep them.  C         in corner region) but safer to keep them.
461          IF ( W_edge ) iMinUpd = 1          IF ( W_edge ) iMinUpd = 1
462          IF ( E_edge ) iMaxUpd = sNx          IF ( E_edge ) iMaxUpd = sNx
463    
464          IF ( S_edge ) THEN          IF ( S_edge ) THEN
465           DO j=1-Oly,0           DO j=1-OLy,0
466            DO i=iMinUpd,iMaxUpd            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)             localTij(i,j) = localTij(i,j)
477       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
478       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
479       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
480       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
481       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -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            ENDDO
485           ENDDO           ENDDO
486          ENDIF          ENDIF
487          IF ( N_edge ) THEN          IF ( N_edge ) THEN
488           DO j=sNy+1,sNy+Oly           DO j=sNy+1,sNy+OLy
489            DO i=iMinUpd,iMaxUpd            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)             localTij(i,j) = localTij(i,j)
500       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
501       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
502       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
503       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
504       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -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            ENDDO
508           ENDDO           ENDDO
509          ENDIF          ENDIF
510    
511         ELSE         ELSE
512  C      do not only update the overlap  C      do not only update the overlap
513          jMinUpd = 1-Oly          jMinUpd = 1-OLy
514          jMaxUpd = sNy+Oly          jMaxUpd = sNy+OLy
515          IF ( interiorOnly .AND. S_edge ) jMinUpd = 1          IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
516          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
517          DO j=jMinUpd,jMaxUpd          DO j=jMinUpd,jMaxUpd
518           DO i=1-Olx+1,sNx+Olx-1           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)             localTij(i,j) = localTij(i,j)
529       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
530       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
531       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
532       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
533       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -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           ENDDO
537          ENDDO          ENDDO
538  C-      keep advective flux (for diagnostics)  C-      keep advective flux (for diagnostics)
539          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
540           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
541            afx(i,j) = af(i,j)            afx(i,j) = af(i,j)
542           ENDDO           ENDDO
543          ENDDO          ENDDO
544    
 #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 )  
 #ifdef ALLOW_PTRACERS  
          ELSEIF (tracerIdentity.GE.GAD_TR1) THEN  
           CALL OBCS_APPLY_PTRACER( bi, bj, k,  
      &         tracerIdentity-GAD_TR1+1, localTij, myThid )  
 #endif /* ALLOW_PTRACERS */  
          ENDIF  
         ENDIF  
 #endif /* ALLOW_OBCS */  
   
545  C-     end if/else update overlap-Only  C-     end if/else update overlap-Only
546         ENDIF         ENDIF
547    
# Line 495  C--   End of X direction Line 550  C--   End of X direction
550    
551  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
552  C--   Y direction  C--   Y direction
553  cph-test  
554  C-     Advective flux in Y  #ifdef ALLOW_AUTODIFF
555          DO j=1-Oly,sNy+Oly  C-     Always reset advective flux in Y
556           DO i=1-Olx,sNx+Olx          DO j=1-OLy,sNy+OLy
557             DO i=1-OLx,sNx+OLx
558            af(i,j) = 0.            af(i,j) = 0.
559           ENDDO           ENDDO
560          ENDDO          ENDDO
 C  
 #ifdef ALLOW_AUTODIFF_TAMC  
561  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
562  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
563  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
564  CADJ STORE af(:,:)  =  CADJ STORE af(:,:)  =
565  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
566  # endif  # endif
567  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
568  C  
569        IF (calc_fluxes_Y) THEN        IF (calc_fluxes_Y) THEN
570    
571  C-     Do not compute fluxes if  C-     Do not compute fluxes if
# Line 520  C   and b) the overlap of myTile are not Line 574  C   and b) the overlap of myTile are not
574         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
575    
576  C-     Internal exchange for calculations in Y  C-     Internal exchange for calculations in Y
 #ifdef MULTIDIM_OLD_VERSION  
         IF ( useCubedSphereExchange ) THEN  
 #else  
577          IF ( overlapOnly ) THEN          IF ( overlapOnly ) THEN
 #endif  
578           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
579       &                              localTij, bi,bj, myThid )       &                              localTij, bi,bj, myThid )
580          ENDIF          ENDIF
581    
582  C-     Advective flux in Y  C-     Advective flux in Y
583          DO j=1-Oly,sNy+Oly  #ifndef ALLOW_AUTODIFF
584           DO i=1-Olx,sNx+Olx          DO j=1-OLy,sNy+OLy
585             DO i=1-OLx,sNx+OLx
586            af(i,j) = 0.            af(i,j) = 0.
587           ENDDO           ENDDO
588          ENDDO          ENDDO
589    #else /* ALLOW_AUTODIFF */
 #ifdef ALLOW_AUTODIFF_TAMC  
590  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
591  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
592  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
593  #endif  #endif
594  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
595    
596          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
597       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
598            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,            CALL GAD_DST2U1_ADV_Y(    bi,bj,k, advectionScheme, .TRUE.,
599       I                           dTtracerLev(k),vTrans,vFld,localTij,       I             deltaTLev(k),vTrans,vFld(1-OLx,1-OLy,k), localTij,
600       O                           af, myThid )       O             af, myThid )
601          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
602            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
603       I                              vTrans, vFld, maskLocS, localTij,       I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
604       O                              af, myThid )       O             af, myThid )
605          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
606            CALL GAD_DST3_ADV_Y(      bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_DST3_ADV_Y(      bi,bj,k, .TRUE., deltaTLev(k),
607       I                              vTrans, vFld, maskLocS, localTij,       I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
608       O                              af, myThid )       O             af, myThid )
609          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
610            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .TRUE., deltaTLev(k),
611       I                              vTrans, vFld, maskLocS, localTij,       I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
612       O                              af, myThid )       O             af, myThid )
613  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF
614          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
615            CALL GAD_OS7MP_ADV_Y(     bi,bj,k, .TRUE., dTtracerLev(k),            CALL GAD_OS7MP_ADV_Y(     bi,bj,k, .TRUE., deltaTLev(k),
616       I                              vTrans, vFld, maskLocS, localTij,       I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
617       O                              af, myThid )       O             af, myThid )
618  #endif  #endif
619          ELSE          ELSE
620           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
621          ENDIF          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  C-     Internal exchange for next calculations in X
633          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
634           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
# Line 580  C-     Advective flux in Y : done Line 639  C-     Advective flux in Y : done
639         ENDIF         ENDIF
640    
641  C-     Update the local tracer field where needed:  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  C      update in overlap-Only
653         IF ( overlapOnly ) THEN         IF ( overlapOnly ) THEN
654          jMinUpd = 1-Oly+1          jMinUpd = 1-OLy+1
655          jMaxUpd = sNy+Oly-1          jMaxUpd = sNy+OLy-1
656  C- notes: these 2 lines below have no real effect (because recip_hFac=0  C- notes: these 2 lines below have no real effect (because recip_hFac=0
657  C         in corner region) but safer to keep them.  C         in corner region) but safer to keep them.
658          IF ( S_edge ) jMinUpd = 1          IF ( S_edge ) jMinUpd = 1
# Line 592  C         in corner region) but safer to Line 660  C         in corner region) but safer to
660    
661          IF ( W_edge ) THEN          IF ( W_edge ) THEN
662           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
663            DO i=1-Olx,0            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)             localTij(i,j) = localTij(i,j)
674       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
675       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
676       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
677       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
678       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -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            ENDDO
682           ENDDO           ENDDO
683          ENDIF          ENDIF
684          IF ( E_edge ) THEN          IF ( E_edge ) THEN
685           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
686            DO i=sNx+1,sNx+Olx            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)             localTij(i,j) = localTij(i,j)
697       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
698       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
699       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
700       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
701       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -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            ENDDO
705           ENDDO           ENDDO
706          ENDIF          ENDIF
707    
708         ELSE         ELSE
709  C      do not only update the overlap  C      do not only update the overlap
710          iMinUpd = 1-Olx          iMinUpd = 1-OLx
711          iMaxUpd = sNx+Olx          iMaxUpd = sNx+OLx
712          IF ( interiorOnly .AND. W_edge ) iMinUpd = 1          IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
713          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
714          DO j=1-Oly+1,sNy+Oly-1          DO j=1-OLy+1,sNy+OLy-1
715           DO i=iMinUpd,iMaxUpd           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)             localTij(i,j) = localTij(i,j)
726       &      -dTtracerLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
727       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
728       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
729       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
730       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -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           ENDDO
734          ENDDO          ENDDO
735  C-      keep advective flux (for diagnostics)  C-      keep advective flux (for diagnostics)
736          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
737           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
738            afy(i,j) = af(i,j)            afy(i,j) = af(i,j)
739           ENDDO           ENDDO
740          ENDDO          ENDDO
741    
 #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 )  
 #ifdef ALLOW_PTRACERS  
          ELSEIF (tracerIdentity.GE.GAD_TR1) THEN  
           CALL OBCS_APPLY_PTRACER( bi, bj, k,  
      &         tracerIdentity-GAD_TR1+1, localTij, myThid )  
 #endif /* ALLOW_PTRACERS */  
          ENDIF  
         ENDIF  
 #endif /* ALLOW_OBCS */  
   
742  C      end if/else update overlap-Only  C      end if/else update overlap-Only
743         ENDIF         ENDIF
744    
# Line 668  C--   End of ipass loop Line 750  C--   End of ipass loop
750    
751        IF ( implicitAdvection ) THEN        IF ( implicitAdvection ) THEN
752  C-    explicit advection is done ; store tendency in gTracer:  C-    explicit advection is done ; store tendency in gTracer:
753          DO j=1-Oly,sNy+Oly  #ifdef GAD_MULTIDIM_COMPRESSIBLE
754           DO i=1-Olx,sNx+Olx          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,bi,bj)=            gTracer(i,j,k,bi,bj)=
759       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
760           ENDDO           ENDDO
761          ENDDO          ENDDO
762        ELSE        ELSE
763  C-    horizontal advection done; store intermediate result in 3D array:  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            localTijk(i,j,k)=localTij(i,j)  #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           ENDDO
771          ENDDO          ENDDO
772        ENDIF        ENDIF
# Line 686  C-    horizontal advection done; store i Line 774  C-    horizontal advection done; store i
774  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
775          IF ( doDiagAdvX ) THEN          IF ( doDiagAdvX ) THEN
776            diagName = 'ADVx'//diagSufx            diagName = 'ADVx'//diagSufx
777            CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( afx, diagName, k,1, 2,bi,bj, myThid )
778          ENDIF          ENDIF
779          IF ( doDiagAdvY ) THEN          IF ( doDiagAdvY ) THEN
780            diagName = 'ADVy'//diagSufx            diagName = 'ADVy'//diagSufx
781            CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( afy, diagName, k,1, 2,bi,bj, myThid )
782          ENDIF          ENDIF
783  #endif  #endif /* ALLOW_DIAGNOSTICS */
784    
785  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
786        IF ( debugLevel .GE. debLevB        IF ( debugLevel .GE. debLevC
787       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE       &   .AND. trIdentity.EQ.GAD_TEMPERATURE
788       &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0       &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
789       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
790       &   .AND. useCubedSphereExchange ) THEN       &   .AND. useCubedSphereExchange ) THEN
# Line 720  C--   kUp    Cycles through 1,2 to point Line 808  C--   kUp    Cycles through 1,2 to point
808  C--   kDown  Cycles through 2,1 to point to w-layer below  C--   kDown  Cycles through 2,1 to point to w-layer below
809          kUp  = 1+MOD(k+1,2)          kUp  = 1+MOD(k+1,2)
810          kDown= 1+MOD(k,2)          kDown= 1+MOD(k,2)
 c       kp1=min(Nr,k+1)  
811          kp1Msk=1.          kp1Msk=1.
812          if (k.EQ.Nr) kp1Msk=0.          IF (k.EQ.Nr) kp1Msk=0.
813    
814  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
815  CADJ STORE rtrans(:,:)  =  CADJ STORE rtrans(:,:)  =
816  CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
 cphCADJ STORE wfld(:,:)  =  
 cphCADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte  
817  #endif  #endif
818    
819  C-- Compute Vertical transport  C-- Compute Vertical transport
820  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
821  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
822          IF ( k.EQ.1 .OR.          IF ( k.EQ.1 .OR.
823       &     (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)       &     (useAIM .AND. trIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
824       &              ) THEN       &              ) THEN
825  #else  #else
826          IF ( k.EQ.1 ) THEN          IF ( k.EQ.1 ) THEN
827  #endif  #endif
828    
 #ifdef ALLOW_AUTODIFF_TAMC  
 cphmultiCADJ STORE wfld(:,:)  =  
 cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
829  C- Surface interface :  C- Surface interface :
830           DO j=1-Oly,sNy+Oly           DO j=1-OLy,sNy+OLy
831            DO i=1-Olx,sNx+Olx            DO i=1-OLx,sNx+OLx
832             rTransKp1(i,j) = kp1Msk*rTrans(i,j)             rTransKp(i,j) = kp1Msk*rTrans(i,j)
            wFld(i,j)   = 0.  
833             rTrans(i,j) = 0.             rTrans(i,j) = 0.
834             fVerT(i,j,kUp) = 0.             fVerT(i,j,kUp) = 0.
835            ENDDO            ENDDO
# Line 758  C- Surface interface : Line 837  C- Surface interface :
837    
838          ELSE          ELSE
839    
 #ifdef ALLOW_AUTODIFF_TAMC  
 cphmultiCADJ STORE wfld(:,:)  =  
 cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
840  C- Interior interface :  C- Interior interface :
841           DO j=1-Oly,sNy+Oly           DO j=1-OLy,sNy+OLy
842            DO i=1-Olx,sNx+Olx            DO i=1-OLx,sNx+OLx
843             rTransKp1(i,j) = kp1Msk*rTrans(i,j)             rTransKp(i,j) = kp1Msk*rTrans(i,j)
844             wFld(i,j)   = wVel(i,j,k,bi,bj)             rTrans(i,j) = wFld(i,j,k)*rA(i,j,bi,bj)
            rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)  
845       &                 *deepFac2F(k)*rhoFacF(k)       &                 *deepFac2F(k)*rhoFacF(k)
846       &                 *maskC(i,j,k-1,bi,bj)       &                 *maskC(i,j,k-1,bi,bj)
847             fVerT(i,j,kUp) = 0.             fVerT(i,j,kUp) = 0.
848            ENDDO            ENDDO
849           ENDDO           ENDDO
850    
 #ifdef ALLOW_GMREDI  
 C--   Residual transp = Bolus transp + Eulerian transp  
          IF (useGMRedi)  
      &     CALL GMREDI_CALC_WFLOW(  
      U                 wFld, rTrans,  
      I                 k, bi, bj, myThid )  
 #endif /* ALLOW_GMREDI */  
   
851  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
852  cphmultiCADJ STORE localTijk(:,:,k)  cphmultiCADJ STORE localT3d(:,:,k)
853  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
854  cphmultiCADJ STORE rTrans(:,:)  cphmultiCADJ STORE rTrans(:,:)
855  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
856  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
857    
858  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
859           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
860       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
861             CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_R(    bi,bj,k, advectionScheme,
862       I                            dTtracerLev(k),rTrans,wFld,localTijk,       I              deltaTLev(k),rTrans,wFld(1-OLx,1-OLy,k),localT3d,
863       O                            fVerT(1-Olx,1-Oly,kUp), myThid )       O              fVerT(1-OLx,1-OLy,kUp), myThid )
864           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
865             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k),             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
866       I                               rTrans, wFld, localTijk,       I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
867       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O              fVerT(1-OLx,1-OLy,kUp), myThid )
868           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
869             CALL GAD_DST3_ADV_R(      bi,bj,k, dTtracerLev(k),             CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTLev(k),
870       I                               rTrans, wFld, localTijk,       I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
871       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O              fVerT(1-OLx,1-OLy,kUp), myThid )
872           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
873             CALL GAD_DST3FL_ADV_R(    bi,bj,k, dTtracerLev(k),             CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTLev(k),
874       I                               rTrans, wFld, localTijk,       I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
875       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O              fVerT(1-OLx,1-OLy,kUp), myThid )
876  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF
877           ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
878             CALL GAD_OS7MP_ADV_R(     bi,bj,k, dTtracerLev(k),             CALL GAD_OS7MP_ADV_R(     bi,bj,k, deltaTLev(k),
879       I                               rTrans, wFld, localTijk,       I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
880       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O              fVerT(1-OLx,1-OLy,kUp), myThid )
881  #endif  #endif
882           ELSE           ELSE
883            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
# Line 823  C- end Surface/Interior if bloc Line 888  C- end Surface/Interior if bloc
888    
889  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
890  cphmultiCADJ STORE rTrans(:,:)  cphmultiCADJ STORE rTrans(:,:)
891  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
892  cphmultiCADJ STORE rTranskp1(:,:)  cphmultiCADJ STORE rTranskp(:,:)
893  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
894  cph --- following storing of fVerT is critical for correct  cph --- following storing of fVerT is critical for correct
895  cph --- gradient with multiDimAdvection  cph --- gradient with multiDimAdvection
896  cph --- Without it, kDown component is not properly recomputed  cph --- Without it, kDown component is not properly recomputed
897  cph --- This is a TAF bug (and no warning available)  cph --- This is a TAF bug (and no warning available)
898  CADJ STORE fVerT(:,:,:)  CADJ STORE fVerT(:,:,:)
899  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  CADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
900  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
901    
902  C--   Divergence of vertical fluxes  C--   Divergence of vertical fluxes
903          DO j=1-Oly,sNy+Oly  #ifdef GAD_MULTIDIM_COMPRESSIBLE
904           DO i=1-Olx,sNx+Olx          DO j=1-OLy,sNy+OLy
905            localTij(i,j) = localTijk(i,j,k)           DO i=1-OLx,sNx+OLx
906       &      -dTtracerLev(k)*recip_rhoFacC(k)            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,bi,bj)=
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,bi,bj) =
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)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
935       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
936       &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)       &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)
937       &         -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(rTransKp(i,j)-rTrans(i,j))
938       &        )*rkSign       &        )*rkSign*maskInC(i,j,bi,bj)
939            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
940       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
941           ENDDO           ENDDO
942          ENDDO          ENDDO
943    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
944    
945  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
946          IF ( doDiagAdvR ) THEN          IF ( doDiagAdvR ) THEN
947            diagName = 'ADVr'//diagSufx            diagName = 'ADVr'//diagSufx
948            CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp),            CALL DIAGNOSTICS_FILL( fVerT(1-OLx,1-OLy,kUp),
949       &                           diagName, k,1, 2,bi,bj, myThid)       &                           diagName, k,1, 2,bi,bj, myThid )
950          ENDIF          ENDIF
951  #endif  #endif /* ALLOW_DIAGNOSTICS */
952    
953  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
954         ENDDO         ENDDO
Legend:
Removed from v.1.58  
changed lines
  Added in v.1.74
  ViewVC Help
Powered by ViewVC 1.1.22