/[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.29 by jmc, Fri Sep 24 16:52:44 2004 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    #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 10  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    
21  C !DESCRIPTION:  C !DESCRIPTION:
22  C Calculates the tendancy of a tracer due to advection.  C Calculates the tendency of a tracer due to advection.
23  C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}  C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
24  C and can only be used for the non-linear advection schemes such as the  C and can only be used for the non-linear advection schemes such as the
25  C direct-space-time method and flux-limiters.  C direct-space-time method and flux-limiters.
26  C  C
27  C The algorithm is as follows:  C The algorithm is as follows:
28  C \begin{itemize}  C \begin{itemize}
# Line 32  C      - \Delta t \partial_r (w\theta^{( Line 35  C      - \Delta t \partial_r (w\theta^{(
35  C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}  C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}
36  C \end{itemize}  C \end{itemize}
37  C  C
38  C The tendancy (output) is over-written by this routine.  C The tendency (output) is over-written by this routine.
39    
40  C !USES: ===============================================================  C !USES: ===============================================================
41        IMPLICIT NONE        IMPLICIT NONE
# Line 41  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 49  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 68  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
83        INTEGER myThid        INTEGER myThid
84    
85  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
86  C  gTracer           :: tendancy 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
99  C  maskLocW      :: 2-D array for mask at West points  C  maskLocW      :: 2-D array for mask at West points
100  C  maskLocS      :: 2-D array for mask at South points  C  maskLocS      :: 2-D array for mask at South points
101  C  iMin,iMax,    :: loop range for called routines  C [iMin,iMax]Upd :: loop range to update tracer field
102  C  jMin,jMax     :: loop range for called routines  C [jMin,jMax]Upd :: loop range to update tracer field
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
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
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  nipass        :: number of passes in multi-dimensional method  C  interiorOnly  :: only update the interior of myTile, but not the edges
120    C  overlapOnly   :: only update the edges of myTile, but not the interior
121    C  npass         :: number of passes in multi-dimensional method
122  C  ipass         :: number of the current pass being made  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        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
127    c     _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128        _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _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)
130        INTEGER iMin,iMax,jMin,jMax        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 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)
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        INTEGER nipass,ipass        LOGICAL interiorOnly, overlapOnly
152        INTEGER myTile, nCFace        INTEGER npass, ipass
153          INTEGER nCFace
154          LOGICAL N_edge, S_edge, E_edge, W_edge
155    #ifdef ALLOW_AUTODIFF_TAMC
156    C     msgBuf     :: Informational/error message buffer
157          CHARACTER*(MAX_LEN_MBUF) msgBuf
158    #endif
159    #ifdef ALLOW_EXCH2
160          INTEGER myTile
161    #endif
162    #ifdef ALLOW_DIAGNOSTICS
163          CHARACTER*8 diagName
164          CHARACTER*4 diagSufx
165          LOGICAL     doDiagAdvX, doDiagAdvY, doDiagAdvR
166    #endif /* ALLOW_DIAGNOSTICS */
167  CEOP  CEOP
168    
169  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
170            act0 = tracerIdentity - 1            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 141  CEOP Line 176  CEOP
176            act3 = myThid - 1            act3 = myThid - 1
177            max3 = nTx*nTy            max3 = nTx*nTy
178            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
179            igadkey = (act0 + 1)            igadkey = act0
180       &                      + act1*max0       &                      + act1*max0
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            endif       &      maxpass, trIdentity
188               CALL PRINT_ERROR( msgBuf, myThid )
189               STOP 'ABNORMAL END: S/R GAD_ADVECTION'
190              ENDIF
191  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
192    
193    #ifdef ALLOW_DIAGNOSTICS
194    C--   Set diagnostics flags and suffix for the current tracer
195          doDiagAdvX = .FALSE.
196          doDiagAdvY = .FALSE.
197          doDiagAdvR = .FALSE.
198          IF ( useDiagnostics ) THEN
199            diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
200            diagName = 'ADVx'//diagSufx
201            doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )
202            diagName = 'ADVy'//diagSufx
203            doDiagAdvY = DIAGNOSTICS_IS_ON( diagName, myThid )
204            diagName = 'ADVr'//diagSufx
205            doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )
206          ENDIF
207    #endif /* ALLOW_DIAGNOSTICS */
208    
209  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
210  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
211  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
# Line 159  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 GAD_MULTIDIM_COMPRESSIBLE
230            localVol(i,j) = 0. _d 0
231    # endif
232            localTij(i,j) = 0. _d 0
233    #endif /* ALLOW_AUTODIFF */
234         ENDDO         ENDDO
235        ENDDO        ENDDO
236    
237        iMin = 1-OLx  C--   Set tile-specific parameters for horizontal fluxes
238        iMax = sNx+OLx        IF (useCubedSphereExchange) THEN
239        jMin = 1-OLy         npass  = 3
240        jMax = sNy+OLy  #ifdef ALLOW_AUTODIFF_TAMC
241           IF ( npass.GT.maxcube ) STOP 'maxcube needs to be = 3'
242    #endif
243    #ifdef ALLOW_EXCH2
244           myTile = W2_myTileList(bi,bj)
245           nCFace = exch2_myFace(myTile)
246           N_edge = exch2_isNedge(myTile).EQ.1
247           S_edge = exch2_isSedge(myTile).EQ.1
248           E_edge = exch2_isEedge(myTile).EQ.1
249           W_edge = exch2_isWedge(myTile).EQ.1
250    #else
251           nCFace = bi
252           N_edge = .TRUE.
253           S_edge = .TRUE.
254           E_edge = .TRUE.
255           W_edge = .TRUE.
256    #endif
257          ELSE
258           npass  = 2
259           nCFace = 0
260           N_edge = .FALSE.
261           S_edge = .FALSE.
262           E_edge = .FALSE.
263           W_edge = .FALSE.
264          ENDIF
265    
266  C--   Start of k loop for horizontal fluxes  C--   Start of k loop for horizontal fluxes
267        DO k=1,Nr        DO k=1,Nr
268  #ifdef ALLOW_AUTODIFF_TAMC  #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         bi,bj,iMin,iMax,jMin,jMax,k,         DO i=1-OLx,sNx+OLx
277       O         xA,yA,uTrans,vTrans,rTrans,maskUp,           xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
278       I         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       &            uTrans, vTrans, bi, bj, k, myThid)  C     anelastic: scaled by rhoFacC (~ mass transport)
285  #endif /* ALLOW_GMREDI */        DO j=1-OLy,sNy+OLy
286           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    
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
 #ifdef ALLOW_EXCH2  
        myTile = W2_myTileList(bi)  
        nCFace = exch2_myFace(myTile)  
 #else  
        nCFace = bi  
 #endif  
       IF (useCubedSphereExchange) THEN  
   
        nipass=3  
 #ifdef ALLOW_AUTODIFF_TAMC  
        if ( nipass.GT.maxcube )  
      &      STOP 'maxcube needs to be = 3'  
 #endif  
       ELSE  
        nipass=1  
       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.
319        DO ipass=1,nipass        DO ipass=1,npass
320  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
321           passkey = ipass + (k-1)      *maxcube           passkey = ipass
322       &                   + (igadkey-1)*maxcube*Nr       &                   + (k-1)      *maxpass
323           IF (nipass .GT. maxpass) THEN       &                   + (igadkey-1)*maxpass*Nr
324            STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'           IF (npass .GT. maxpass) THEN
325              STOP 'GAD_ADVECTION: npass > maxcube. check tamc.h'
326           ENDIF           ENDIF
327  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
328    
329        IF (nipass.EQ.3) THEN        interiorOnly = .FALSE.
330         calc_fluxes_X=.FALSE.        overlapOnly  = .FALSE.
331         calc_fluxes_Y=.FALSE.        IF (useCubedSphereExchange) THEN
332         IF (ipass.EQ.1 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.2) ) THEN  C-    CubedSphere : pass 3 times, with partial update of local tracer field
333          calc_fluxes_X=.TRUE.         IF (ipass.EQ.1) THEN
334         ELSEIF (ipass.EQ.1 .AND. (nCFace.EQ.4 .OR. nCFace.EQ.5) ) THEN          overlapOnly  = MOD(nCFace,3).EQ.0
335          calc_fluxes_Y=.TRUE.          interiorOnly = MOD(nCFace,3).NE.0
336         ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.6) ) THEN          calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
337          calc_fluxes_Y=.TRUE.          calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
338         ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.3 .OR. nCFace.EQ.4) ) THEN         ELSEIF (ipass.EQ.2) THEN
339          calc_fluxes_X=.TRUE.          overlapOnly  = MOD(nCFace,3).EQ.2
340         ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.2 .OR. nCFace.EQ.3) ) THEN          interiorOnly = MOD(nCFace,3).EQ.1
341          calc_fluxes_Y=.TRUE.          calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
342         ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.5 .OR. nCFace.EQ.6) ) THEN          calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
343          calc_fluxes_X=.TRUE.         ELSE
344            interiorOnly = .TRUE.
345            calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
346            calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
347         ENDIF         ENDIF
348        ELSE        ELSE
349         calc_fluxes_X=.TRUE.  C-    not CubedSphere
350         calc_fluxes_Y=.TRUE.          calc_fluxes_X = MOD(ipass,2).EQ.1
351            calc_fluxes_Y = .NOT.calc_fluxes_X
352        ENDIF        ENDIF
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    
357    #ifdef ALLOW_AUTODIFF
358    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.
362             ENDDO
363            ENDDO
364    # ifndef DISABLE_MULTIDIM_ADVECTION
365    CADJ STORE localTij(:,:)  =
366    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
367    CADJ STORE af(:,:)  =
368    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
369    # endif
370    #endif /* ALLOW_AUTODIFF */
371    
372        IF (calc_fluxes_X) THEN        IF (calc_fluxes_X) THEN
373    
374  C--   Internal exchange for calculations in X  C-     Do not compute fluxes if
375        IF (useCubedSphereExchange) THEN  C       a) needed in overlap only
376           CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )  C   and b) the overlap of myTile are not cube-face Edges
377        ENDIF         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
378    
379    C-     Internal exchange for calculations in X
380            IF ( overlapOnly ) THEN
381             CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
382         &                              localTij, bi,bj, myThid )
383            ENDIF
384    
385  C-    Advective flux in X  C-     Advective flux in X
386        DO j=1-Oly,sNy+Oly  #ifndef ALLOW_AUTODIFF
387         DO i=1-Olx,sNx+Olx          DO j=1-OLy,sNy+OLy
388          afx(i,j) = 0.           DO i=1-OLx,sNx+OLx
389         ENDDO            af(i,j) = 0.
390        ENDDO           ENDDO
391            ENDDO
392    #else /* ALLOW_AUTODIFF */
393    # ifndef DISABLE_MULTIDIM_ADVECTION
394    CADJ STORE localTij(:,:)  =
395    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
396    # endif
397    #endif /* ALLOW_AUTODIFF */
398    
399  #ifdef ALLOW_AUTODIFF_TAMC          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
400  #ifndef DISABLE_MULTIDIM_ADVECTION       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
401  CADJ STORE localTij(:,:)  =            CALL GAD_DST2U1_ADV_X(    bi,bj,k, advectionScheme, .TRUE.,
402  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte       I             deltaTLev(k),uTrans,uFld(1-OLx,1-OLy,k), localTij,
403         O             af, myThid )
404            ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
405              CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
406         I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
407         O             af, myThid )
408            ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
409              CALL GAD_DST3_ADV_X(      bi,bj,k, .TRUE., deltaTLev(k),
410         I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
411         O             af, myThid )
412            ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
413              CALL GAD_DST3FL_ADV_X(    bi,bj,k, .TRUE., deltaTLev(k),
414         I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
415         O             af, myThid )
416    #ifndef ALLOW_AUTODIFF
417            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
418              CALL GAD_OS7MP_ADV_X(     bi,bj,k, .TRUE., deltaTLev(k),
419         I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
420         O             af, myThid )
421  #endif  #endif
422            ELSE
423             STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
424            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
436            IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
437             CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
438         &                              localTij, bi,bj, myThid )
439            ENDIF
440    
441    C-     Advective flux in X : done
442           ENDIF
443    
444    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 */  #endif /* ALLOW_AUTODIFF_TAMC */
454    
455        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN  C      update in overlap-Only
456          CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, deltaTtracer,         IF ( overlapOnly ) THEN
457       I                            uTrans, uVel, maskLocW, localTij,          iMinUpd = 1-OLx+1
458       O                            afx, myThid )          iMaxUpd = sNx+OLx-1
459        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN  C- notes: these 2 lines below have no real effect (because recip_hFac=0
460          CALL GAD_DST3_ADV_X(      bi,bj,k, deltaTtracer,  C         in corner region) but safer to keep them.
461       I                            uTrans, uVel, maskLocW, localTij,          IF ( W_edge ) iMinUpd = 1
462       O                            afx, myThid )          IF ( E_edge ) iMaxUpd = sNx
463        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN  
464          CALL GAD_DST3FL_ADV_X(    bi,bj,k, deltaTtracer,          IF ( S_edge ) THEN
465       I                            uTrans, uVel, maskLocW, localTij,           DO j=1-OLy,0
466       O                            afx, myThid )            DO i=iMinUpd,iMaxUpd
467        ELSE  #ifdef GAD_MULTIDIM_COMPRESSIBLE
468         STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'             tmpTrac = localTij(i,j)*localVol(i,j)
469        ENDIF       &      -deltaTLev(k)*( af(i+1,j) - af(i,j) )
470         &                   *maskInC(i,j,bi,bj)
471               localVol(i,j) = localVol(i,j)
472         &      -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) )
473         &                   *maskInC(i,j,bi,bj)
474               localTij(i,j) = tmpTrac/localVol(i,j)
475    #else /* GAD_MULTIDIM_COMPRESSIBLE */
476               localTij(i,j) = localTij(i,j)
477         &      -deltaTLev(k)*recip_rhoFacC(k)
478         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
479         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
480         &       *( af(i+1,j)-af(i,j)
481         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
482         &        )*maskInC(i,j,bi,bj)
483    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
484              ENDDO
485             ENDDO
486            ENDIF
487            IF ( N_edge ) THEN
488             DO j=sNy+1,sNy+OLy
489              DO i=iMinUpd,iMaxUpd
490    #ifdef GAD_MULTIDIM_COMPRESSIBLE
491               tmpTrac = localTij(i,j)*localVol(i,j)
492         &      -deltaTLev(k)*( af(i+1,j) - af(i,j) )
493         &                   *maskInC(i,j,bi,bj)
494               localVol(i,j) = localVol(i,j)
495         &      -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) )
496         &                   *maskInC(i,j,bi,bj)
497               localTij(i,j) = tmpTrac/localVol(i,j)
498    #else /* GAD_MULTIDIM_COMPRESSIBLE */
499               localTij(i,j) = localTij(i,j)
500         &      -deltaTLev(k)*recip_rhoFacC(k)
501         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
502         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
503         &       *( af(i+1,j)-af(i,j)
504         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
505         &        )*maskInC(i,j,bi,bj)
506    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
507              ENDDO
508             ENDDO
509            ENDIF
510    
511        DO j=1-Oly,sNy+Oly         ELSE
512         DO i=1-Olx,sNx+Olx-1  C      do not only update the overlap
513          localTij(i,j)=localTij(i,j)-deltaTtracer*          jMinUpd = 1-OLy
514       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)          jMaxUpd = sNy+OLy
515       &    *recip_rA(i,j,bi,bj)          IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
516       &    *( afx(i+1,j)-afx(i,j)          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
517       &      -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))          DO j=jMinUpd,jMaxUpd
518       &     )           DO i=1-OLx+1,sNx+OLx-1
519         ENDDO  #ifdef GAD_MULTIDIM_COMPRESSIBLE
520        ENDDO             tmpTrac = localTij(i,j)*localVol(i,j)
521         &      -deltaTLev(k)*( af(i+1,j) - af(i,j) )
522         &                   *maskInC(i,j,bi,bj)
523               localVol(i,j) = localVol(i,j)
524         &      -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) )
525         &                   *maskInC(i,j,bi,bj)
526               localTij(i,j) = tmpTrac/localVol(i,j)
527    #else /* GAD_MULTIDIM_COMPRESSIBLE */
528               localTij(i,j) = localTij(i,j)
529         &      -deltaTLev(k)*recip_rhoFacC(k)
530         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
531         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
532         &       *( af(i+1,j)-af(i,j)
533         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
534         &        )*maskInC(i,j,bi,bj)
535    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
536             ENDDO
537            ENDDO
538    C-      keep advective flux (for diagnostics)
539            DO j=1-OLy,sNy+OLy
540             DO i=1-OLx,sNx+OLx
541              afx(i,j) = af(i,j)
542             ENDDO
543            ENDDO
544    
545  #ifdef ALLOW_OBCS  C-     end if/else update overlap-Only
546  C--   Apply open boundary conditions         ENDIF
       IF (useOBCS) THEN  
        IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN  
         CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )  
        ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN  
         CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )  
        END IF  
       END IF  
 #endif /* ALLOW_OBCS */  
547    
548  C--   End of X direction  C--   End of X direction
549        ENDIF        ENDIF
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
       IF (calc_fluxes_Y) THEN  
553    
554  C--   Internal exchange for calculations in Y  #ifdef ALLOW_AUTODIFF
555        IF (useCubedSphereExchange) THEN  C-     Always reset advective flux in Y
556           CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )          DO j=1-OLy,sNy+OLy
557        ENDIF           DO i=1-OLx,sNx+OLx
558              af(i,j) = 0.
559             ENDDO
560            ENDDO
561    # ifndef DISABLE_MULTIDIM_ADVECTION
562    CADJ STORE localTij(:,:)  =
563    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
564    CADJ STORE af(:,:)  =
565    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
566    # endif
567    #endif /* ALLOW_AUTODIFF */
568    
569  C-    Advective flux in Y        IF (calc_fluxes_Y) THEN
570        DO j=1-Oly,sNy+Oly  
571         DO i=1-Olx,sNx+Olx  C-     Do not compute fluxes if
572          afy(i,j) = 0.  C       a) needed in overlap only
573         ENDDO  C   and b) the overlap of myTile are not cube-face edges
574        ENDDO         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
575    
576    C-     Internal exchange for calculations in Y
577            IF ( overlapOnly ) THEN
578             CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
579         &                              localTij, bi,bj, myThid )
580            ENDIF
581    
582  #ifdef ALLOW_AUTODIFF_TAMC  C-     Advective flux in Y
583    #ifndef ALLOW_AUTODIFF
584            DO j=1-OLy,sNy+OLy
585             DO i=1-OLx,sNx+OLx
586              af(i,j) = 0.
587             ENDDO
588            ENDDO
589    #else /* ALLOW_AUTODIFF */
590  #ifndef DISABLE_MULTIDIM_ADVECTION  #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_FLUX_LIMIT) THEN          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
597          CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, deltaTtracer,       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
598       I                            vTrans, vVel, maskLocS, localTij,            CALL GAD_DST2U1_ADV_Y(    bi,bj,k, advectionScheme, .TRUE.,
599       O                            afy, myThid )       I             deltaTLev(k),vTrans,vFld(1-OLx,1-OLy,k), localTij,
600        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       O             af, myThid )
601          CALL GAD_DST3_ADV_Y(      bi,bj,k, deltaTtracer,          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
602       I                            vTrans, vVel, maskLocS, localTij,            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
603       O                            afy, myThid )       I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
604        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       O             af, myThid )
605          CALL GAD_DST3FL_ADV_Y(    bi,bj,k, deltaTtracer,          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
606       I                            vTrans, vVel, maskLocS, localTij,            CALL GAD_DST3_ADV_Y(      bi,bj,k, .TRUE., deltaTLev(k),
607       O                            afy, myThid )       I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
608        ELSE       O             af, myThid )
609         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
610        ENDIF            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .TRUE., deltaTLev(k),
611         I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
612        DO j=1-Oly,sNy+Oly-1       O             af, myThid )
613         DO i=1-Olx,sNx+Olx  #ifndef ALLOW_AUTODIFF
614          localTij(i,j)=localTij(i,j)-deltaTtracer*          ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
615       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)            CALL GAD_OS7MP_ADV_Y(     bi,bj,k, .TRUE., deltaTLev(k),
616       &    *recip_rA(i,j,bi,bj)       I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
617       &    *( afy(i,j+1)-afy(i,j)       O             af, myThid )
618       &      -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))  #endif
619       &     )          ELSE
620         ENDDO           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
621        ENDDO          ENDIF
622    
623  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
624  C--   Apply open boundary conditions          IF ( useOBCS ) THEN
625        IF (useOBCS) THEN  C-      replace advective flux with 1st order upwind scheme estimate
626         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN            CALL OBCS_U1_ADV_TRACER( .FALSE., trIdentity, bi, bj, k,
627          CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )       I                             maskLocS, vTrans, localTij,
628         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN       U                             af, myThid )
629          CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )          ENDIF
        END IF  
       END IF  
630  #endif /* ALLOW_OBCS */  #endif /* ALLOW_OBCS */
631    
632    C-     Internal exchange for next calculations in X
633            IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
634             CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
635         &                              localTij, bi,bj, myThid )
636            ENDIF
637    
638    C-     Advective flux in Y : done
639           ENDIF
640    
641    C-     Update the local tracer field where needed:
642    C      use "maksInC" to prevent updating tracer field in OB regions
643    #ifdef ALLOW_AUTODIFF_TAMC
644    # ifdef GAD_MULTIDIM_COMPRESSIBLE
645    CADJ STORE localVol(:,:)  =
646    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
647    CADJ STORE localTij(:,:)  =
648    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
649    # endif
650    #endif /* ALLOW_AUTODIFF_TAMC */
651    
652    C      update in overlap-Only
653           IF ( overlapOnly ) THEN
654            jMinUpd = 1-OLy+1
655            jMaxUpd = sNy+OLy-1
656    C- notes: these 2 lines below have no real effect (because recip_hFac=0
657    C         in corner region) but safer to keep them.
658            IF ( S_edge ) jMinUpd = 1
659            IF ( N_edge ) jMaxUpd = sNy
660    
661            IF ( W_edge ) THEN
662             DO j=jMinUpd,jMaxUpd
663              DO i=1-OLx,0
664    #ifdef GAD_MULTIDIM_COMPRESSIBLE
665               tmpTrac = localTij(i,j)*localVol(i,j)
666         &      -deltaTLev(k)*( af(i,j+1) - af(i,j) )
667         &                   *maskInC(i,j,bi,bj)
668               localVol(i,j) = localVol(i,j)
669         &      -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) )
670         &                   *maskInC(i,j,bi,bj)
671               localTij(i,j) = tmpTrac/localVol(i,j)
672    #else /* GAD_MULTIDIM_COMPRESSIBLE */
673               localTij(i,j) = localTij(i,j)
674         &      -deltaTLev(k)*recip_rhoFacC(k)
675         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
676         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
677         &       *( af(i,j+1)-af(i,j)
678         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
679         &        )*maskInC(i,j,bi,bj)
680    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
681              ENDDO
682             ENDDO
683            ENDIF
684            IF ( E_edge ) THEN
685             DO j=jMinUpd,jMaxUpd
686              DO i=sNx+1,sNx+OLx
687    #ifdef GAD_MULTIDIM_COMPRESSIBLE
688               tmpTrac = localTij(i,j)*localVol(i,j)
689         &      -deltaTLev(k)*( af(i,j+1) - af(i,j) )
690         &                   *maskInC(i,j,bi,bj)
691               localVol(i,j) = localVol(i,j)
692         &      -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) )
693         &                   *maskInC(i,j,bi,bj)
694               localTij(i,j) = tmpTrac/localVol(i,j)
695    #else /* GAD_MULTIDIM_COMPRESSIBLE */
696               localTij(i,j) = localTij(i,j)
697         &      -deltaTLev(k)*recip_rhoFacC(k)
698         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
699         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
700         &       *( af(i,j+1)-af(i,j)
701         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
702         &        )*maskInC(i,j,bi,bj)
703    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
704              ENDDO
705             ENDDO
706            ENDIF
707    
708           ELSE
709    C      do not only update the overlap
710            iMinUpd = 1-OLx
711            iMaxUpd = sNx+OLx
712            IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
713            IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
714            DO j=1-OLy+1,sNy+OLy-1
715             DO i=iMinUpd,iMaxUpd
716    #ifdef GAD_MULTIDIM_COMPRESSIBLE
717               tmpTrac = localTij(i,j)*localVol(i,j)
718         &      -deltaTLev(k)*( af(i,j+1) - af(i,j) )
719         &                   *maskInC(i,j,bi,bj)
720               localVol(i,j) = localVol(i,j)
721         &      -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) )
722         &                   *maskInC(i,j,bi,bj)
723               localTij(i,j) = tmpTrac/localVol(i,j)
724    #else /* GAD_MULTIDIM_COMPRESSIBLE */
725               localTij(i,j) = localTij(i,j)
726         &      -deltaTLev(k)*recip_rhoFacC(k)
727         &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
728         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
729         &       *( af(i,j+1)-af(i,j)
730         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
731         &        )*maskInC(i,j,bi,bj)
732    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
733             ENDDO
734            ENDDO
735    C-      keep advective flux (for diagnostics)
736            DO j=1-OLy,sNy+OLy
737             DO i=1-OLx,sNx+OLx
738              afy(i,j) = af(i,j)
739             ENDDO
740            ENDDO
741    
742    C      end if/else update overlap-Only
743           ENDIF
744    
745  C--   End of Y direction  C--   End of Y direction
746        ENDIF        ENDIF
747    
# Line 392  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))/deltaTtracer       &     (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
771          ENDDO          ENDDO
        ENDDO  
772        ENDIF        ENDIF
773    
774    #ifdef ALLOW_DIAGNOSTICS
775            IF ( doDiagAdvX ) THEN
776              diagName = 'ADVx'//diagSufx
777              CALL DIAGNOSTICS_FILL( afx, diagName, k,1, 2,bi,bj, myThid )
778            ENDIF
779            IF ( doDiagAdvY ) THEN
780              diagName = 'ADVy'//diagSufx
781              CALL DIAGNOSTICS_FILL( afy, diagName, k,1, 2,bi,bj, myThid )
782            ENDIF
783    #endif /* ALLOW_DIAGNOSTICS */
784    
785  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
786        IF ( debugLevel .GE. debLevB        IF ( debugLevel .GE. debLevC
787       &   .AND. k.EQ.3 .AND. myIter.EQ.1+nIter0       &   .AND. trIdentity.EQ.GAD_TEMPERATURE
788         &   .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
791          CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION',          CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION',
# Line 425  C---+----1----+----2----+----3----+----4 Line 801  C---+----1----+----2----+----3----+----4
801        IF ( .NOT.implicitAdvection ) THEN        IF ( .NOT.implicitAdvection ) THEN
802  C--   Start of k loop for vertical flux  C--   Start of k loop for vertical flux
803         DO k=Nr,1,-1         DO k=Nr,1,-1
804  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
805           kkey = (igadkey-1)*Nr + k           kkey = (igadkey-1)*Nr + (Nr-k+1)
806  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
807  C--   kup    Cycles through 1,2 to point to w-layer above  C--   kUp    Cycles through 1,2 to point to w-layer above
808  C--   kDown  Cycles through 2,1 to point to w-layer below  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
815    CADJ STORE rtrans(:,:)  =
816    CADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
817    #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    
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)
833             rTrans(i,j) = 0.             rTrans(i,j) = 0.
834             fVerT(i,j,kUp) = 0.             fVerT(i,j,kUp) = 0.
835            ENDDO            ENDDO
836           ENDDO           ENDDO
837    
838          ELSE          ELSE
 C- Interior interface :  
839    
840           DO j=1-Oly,sNy+Oly  C- Interior interface :
841            DO i=1-Olx,sNx+Olx           DO j=1-OLy,sNy+OLy
842             rTransKp1(i,j) = kp1Msk*rTrans(i,j)            DO i=1-OLx,sNx+OLx
843             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)             rTransKp(i,j) = kp1Msk*rTrans(i,j)
844               rTrans(i,j) = wFld(i,j,k)*rA(i,j,bi,bj)
845         &                 *deepFac2F(k)*rhoFacF(k)
846       &                 *maskC(i,j,k-1,bi,bj)       &                 *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    
851  #ifdef ALLOW_GMREDI  #ifdef ALLOW_AUTODIFF_TAMC
852  C--   Residual transp = Bolus transp + Eulerian transp  cphmultiCADJ STORE localT3d(:,:,k)
853           IF (useGMRedi)  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
854       &   CALL GMREDI_CALC_WFLOW(  cphmultiCADJ STORE rTrans(:,:)
855       &                    rTrans, bi, bj, k, myThid)  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
 #endif /* ALLOW_GMREDI */  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE localTijk(:,:,k)    
 CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  
 CADJ STORE rTrans(:,:)    
 CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=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_FLUX_LIMIT) THEN           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
860  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
861             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTtracer,             CALL GAD_DST2U1_ADV_R(    bi,bj,k, advectionScheme,
862       I                               rTrans, wVel, 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_DST3 ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
865             CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTtracer,             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
866       I                               rTrans, wVel, 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_FLUX_LIMIT ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
869             CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTtracer,             CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTLev(k),
870       I                               rTrans, wVel, 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
873               CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTLev(k),
874         I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
875         O              fVerT(1-OLx,1-OLy,kUp), myThid )
876    #ifndef ALLOW_AUTODIFF
877             ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
878               CALL GAD_OS7MP_ADV_R(     bi,bj,k, deltaTLev(k),
879         I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
880         O              fVerT(1-OLx,1-OLy,kUp), myThid )
881    #endif
882           ELSE           ELSE
883            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
884           ENDIF           ENDIF
# Line 502  C---+----1----+----2----+----3----+----4 Line 886  C---+----1----+----2----+----3----+----4
886  C- end Surface/Interior if bloc  C- end Surface/Interior if bloc
887          ENDIF          ENDIF
888    
889  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
890  CADJ STORE rTrans(:,:)    cphmultiCADJ STORE rTrans(:,:)
891  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
892  CADJ STORE rTranskp1(:,:)    cphmultiCADJ STORE rTranskp(:,:)
893  CADJ &     = 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
895    cph --- gradient with multiDimAdvection
896    cph --- Without it, kDown component is not properly recomputed
897    cph --- This is a TAF bug (and no warning available)
898    CADJ STORE fVerT(:,:,:)
899    CADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
900  #endif /* ALLOW_AUTODIFF_TAMC */  #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)-deltaTtracer*           DO i=1-OLx,sNx+OLx
906       &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)            tmpTrac = localT3d(i,j,k)*locVol3d(i,j,k)
907       &     *recip_rA(i,j,bi,bj)       &      -deltaTLev(k)*( fVerT(i,j,kDown)-fVerT(i,j,kUp) )
908       &     *( fVerT(i,j,kUp)-fVerT(i,j,kDown)       &                   *rkSign*maskInC(i,j,bi,bj)
909       &       -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))            localVol(i,j) = locVol3d(i,j,k)
910       &      )*rkFac       &      -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)
935         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
936         &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)
937         &         -tracer(i,j,k,bi,bj)*(rTransKp(i,j)-rTrans(i,j))
938         &        )*rkSign*maskInC(i,j,bi,bj)
939            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
940       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer       &     (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
946            IF ( doDiagAdvR ) THEN
947              diagName = 'ADVr'//diagSufx
948              CALL DIAGNOSTICS_FILL( fVerT(1-OLx,1-OLy,kUp),
949         &                           diagName, k,1, 2,bi,bj, myThid )
950            ENDIF
951    #endif /* ALLOW_DIAGNOSTICS */
952    
953  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
954         ENDDO         ENDDO
955  C--   end of if not.implicitAdvection block  C--   end of if not.implicitAdvection block
956        ENDIF        ENDIF
957    
958        RETURN        RETURN
959        END        END
Legend:
Removed from v.1.29  
changed lines
  Added in v.1.74
  ViewVC Help
Powered by ViewVC 1.1.22