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