/[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.41 by jmc, Sun Jun 18 23:27:44 2006 UTC revision 1.74 by jmc, Fri Apr 4 20:28:13 2014 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5  #undef MULTIDIM_OLD_VERSION  #ifdef ALLOW_AUTODIFF
6    # include "AUTODIFF_OPTIONS.h"
7    #endif
8    
9  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
10  CBOP  CBOP
# Line 11  C !ROUTINE: GAD_ADVECTION Line 13  C !ROUTINE: GAD_ADVECTION
13  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
14        SUBROUTINE GAD_ADVECTION(        SUBROUTINE GAD_ADVECTION(
15       I     implicitAdvection, advectionScheme, vertAdvecScheme,       I     implicitAdvection, advectionScheme, vertAdvecScheme,
16       I     tracerIdentity,       I     trIdentity, deltaTLev,
17       I     uVel, vVel, wVel, tracer,       I     uFld, vFld, wFld, tracer,
18       O     gTracer,       O     gTracer,
19       I     bi,bj, myTime,myIter,myThid)       I     bi,bj, myTime,myIter,myThid)
20    
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.
# Line 33  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 42  C !USES: =============================== Line 44  C !USES: ===============================
44  #include "PARAMS.h"  #include "PARAMS.h"
45  #include "GRID.h"  #include "GRID.h"
46  #include "GAD.h"  #include "GAD.h"
47  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
48  # include "tamc.h"  # include "tamc.h"
49  # include "tamc_keys.h"  # include "tamc_keys.h"
50  # ifdef ALLOW_PTRACERS  # ifdef ALLOW_PTRACERS
# Line 50  C !USES: =============================== Line 52  C !USES: ===============================
52  # endif  # endif
53  #endif  #endif
54  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
55    #include "W2_EXCH2_SIZE.h"
56  #include "W2_EXCH2_TOPOLOGY.h"  #include "W2_EXCH2_TOPOLOGY.h"
 #include "W2_EXCH2_PARAMS.h"  
57  #endif /* ALLOW_EXCH2 */  #endif /* ALLOW_EXCH2 */
58    
59  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
60  C  implicitAdvection :: implicit vertical advection (later on)  C  implicitAdvection :: implicit vertical advection (later on)
61  C  advectionScheme   :: advection scheme to use (Horizontal plane)  C  advectionScheme   :: advection scheme to use (Horizontal plane)
62  C  vertAdvecScheme   :: advection scheme to use (vertical direction)  C  vertAdvecScheme   :: advection scheme to use (vertical direction)
63  C  tracerIdentity    :: tracer identifier (required only for OBCS)  C  trIdentity        :: tracer identifier
64  C  uVel              :: velocity, zonal component  C  uFld              :: Advection velocity field, zonal component
65  C  vVel              :: velocity, meridional component  C  vFld              :: Advection velocity field, meridional component
66  C  wVel              :: velocity, vertical component  C  wFld              :: Advection velocity field, vertical component
67  C  tracer            :: tracer field  C  tracer            :: tracer field
68  C  bi,bj             :: tile indices  C  bi,bj             :: tile indices
69  C  myTime            :: current time  C  myTime            :: current time
# Line 69  C  myIter            :: iteration number Line 71  C  myIter            :: iteration number
71  C  myThid            :: thread number  C  myThid            :: thread number
72        LOGICAL implicitAdvection        LOGICAL implicitAdvection
73        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
74        INTEGER tracerIdentity        INTEGER trIdentity
75        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL deltaTLev(Nr)
76        _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
77        _RL wVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
78        _RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
79          _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
80        INTEGER bi,bj        INTEGER bi,bj
81        _RL myTime        _RL myTime
82        INTEGER myIter        INTEGER myIter
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
 C  iMin,iMax,    :: loop range for called routines  
 C  jMin,jMax     :: loop range for called routines  
101  C [iMin,iMax]Upd :: loop range to update tracer field  C [iMin,iMax]Upd :: loop range to update tracer field
102  C [jMin,jMax]Upd :: loop range to update tracer field  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
 C  uFld,vFld     :: 2-D local copy of horizontal velocity, U,V components  
 C  wFld          :: 2-D local copy of vertical velocity  
107  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
108  C  rTrans        :: 2-D arrays of volume transports at W points  C  rTrans        :: 2-D arrays of volume transports at W points
109  C  rTransKp1     :: vertical volume transport at interface k+1  C  rTransKp      :: vertical volume transport at interface k+1
110  C  af            :: 2-D array for horizontal advective flux  C  af            :: 2-D array for horizontal advective flux
111  C  afx           :: 2-D array for horizontal advective flux, x direction  C  afx           :: 2-D array for horizontal advective flux, x direction
112  C  afy           :: 2-D array for horizontal advective flux, y direction  C  afy           :: 2-D array for horizontal advective flux, y direction
113  C  fVerT         :: 2 1/2D arrays for vertical advective flux  C  fVerT         :: 2 1/2D arrays for vertical advective flux
114  C  localTij      :: 2-D array, temporary local copy of tracer fld  C  localTij      :: 2-D array, temporary local copy of tracer field
115  C  localTijk     :: 3-D array, temporary local copy of tracer fld  C  localT3d      :: 3-D array, temporary local copy of tracer field
116  C  kp1Msk        :: flag (0,1) for over-riding mask for W levels  C  kp1Msk        :: flag (0,1) for over-riding mask for W levels
117  C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir  C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
118  C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir  C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
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
# Line 121  C [N,S,E,W]_edge :: true if N,S,E,W edge Line 127  C [N,S,E,W]_edge :: true if N,S,E,W edge
127  c     _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)
       _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL wFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
134        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139        _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140        _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
142        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
143        _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL localT3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
144    #ifdef GAD_MULTIDIM_COMPRESSIBLE
145          _RL tmpTrac
146          _RL localVol(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147          _RL locVol3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
148    #endif
149        _RL kp1Msk        _RL kp1Msk
150        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
151        LOGICAL interiorOnly, overlapOnly        LOGICAL interiorOnly, overlapOnly
152        INTEGER nipass,ipass        INTEGER npass, ipass
153        INTEGER nCFace        INTEGER nCFace
154        LOGICAL N_edge, S_edge, E_edge, W_edge        LOGICAL N_edge, S_edge, E_edge, W_edge
155    #ifdef ALLOW_AUTODIFF_TAMC
156    C     msgBuf     :: Informational/error message buffer
157          CHARACTER*(MAX_LEN_MBUF) msgBuf
158    #endif
159  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
160        INTEGER myTile        INTEGER myTile
161  #endif  #endif
162  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
163        CHARACTER*8 diagName        CHARACTER*8 diagName
164        CHARACTER*4 GAD_DIAG_SUFX, diagSufx        CHARACTER*4 diagSufx
165        EXTERNAL    GAD_DIAG_SUFX        LOGICAL     doDiagAdvX, doDiagAdvY, doDiagAdvR
166  #endif  #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 165  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  #ifdef ALLOW_DIAGNOSTICS
194  C--   Set diagnostic suffix for the current tracer  C--   Set diagnostics flags and suffix for the current tracer
195          doDiagAdvX = .FALSE.
196          doDiagAdvY = .FALSE.
197          doDiagAdvR = .FALSE.
198        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
199          diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )          diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
200            diagName = 'ADVx'//diagSufx
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        ENDIF
207  #endif  #endif /* ALLOW_DIAGNOSTICS */
208    
209  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
210  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
# Line 190  C     point numbers. This prevents spuri Line 213  C     point numbers. This prevents spuri
213  C     uninitialised but inert locations.  C     uninitialised but inert locations.
214        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
215         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
216          xA(i,j)      = 0. _d 0  C-    xA,yA,vFld,uTrans,vTrans are set over the full domain
217          yA(i,j)      = 0. _d 0  C      => no need for extra initialisation
218          uTrans(i,j)  = 0. _d 0  c       xA(i,j)      = 0. _d 0
219          vTrans(i,j)  = 0. _d 0  c       yA(i,j)      = 0. _d 0
220    c       uTrans(i,j)  = 0. _d 0
221    c       vTrans(i,j)  = 0. _d 0
222    C-    rTransKp is set over the full domain: no need for extra initialisation
223    c       rTransKp(i,j)= 0. _d 0
224    C-    rTrans and fVerT need to be initialised to zero:
225          rTrans(i,j)  = 0. _d 0          rTrans(i,j)  = 0. _d 0
226          fVerT(i,j,1) = 0. _d 0          fVerT(i,j,1) = 0. _d 0
227          fVerT(i,j,2) = 0. _d 0          fVerT(i,j,2) = 0. _d 0
228          rTransKp1(i,j)= 0. _d 0  #ifdef ALLOW_AUTODIFF
229  #ifdef ALLOW_AUTODIFF_TAMC  # ifdef GAD_MULTIDIM_COMPRESSIBLE
230            localVol(i,j) = 0. _d 0
231    # endif
232          localTij(i,j) = 0. _d 0          localTij(i,j) = 0. _d 0
233  #endif  #endif /* ALLOW_AUTODIFF */
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 225  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 = bi         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 (        DO j=1-OLy,sNy+OLy
276       I         uVel, vVel,         DO i=1-OLx,sNx+OLx
277       O         uFld, vFld, uTrans, vTrans, xA, yA,           xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
278       I         k,bi,bj, myThid )       &           *drF(k)*_hFacW(i,j,k,bi,bj)
279             yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
280  #ifdef ALLOW_GMREDI       &           *drF(k)*_hFacS(i,j,k,bi,bj)
281  C--   Residual transp = Bolus transp + Eulerian transp         ENDDO
282        IF (useGMRedi)        ENDDO
283       &   CALL GMREDI_CALC_UVFLOW(  C--   Calculate "volume transports" through tracer cell faces.
284       U                  uFld, vFld, uTrans, vTrans,  C     anelastic: scaled by rhoFacC (~ mass transport)
285       I                  k, bi, bj, myThid )        DO j=1-OLy,sNy+OLy
286  #endif /* ALLOW_GMREDI */         DO i=1-OLx,sNx+OLx
287             uTrans(i,j) = uFld(i,j,k)*xA(i,j)*rhoFacC(k)
288             vTrans(i,j) = vFld(i,j,k)*yA(i,j)*rhoFacC(k)
289           ENDDO
290          ENDDO
291    
292  C--   Make local copy of tracer array and mask West & South  C--   Make local copy of tracer array and mask West & South
293        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
294         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
295           localTij(i,j)=tracer(i,j,k,bi,bj)           localTij(i,j) = tracer(i,j,k,bi,bj)
296           maskLocW(i,j)=maskW(i,j,k,bi,bj)  #ifdef GAD_MULTIDIM_COMPRESSIBLE
297           maskLocS(i,j)=maskS(i,j,k,bi,bj)           localVol(i,j) = rA(i,j,bi,bj)*deepFac2C(k)
298         &                  *rhoFacC(k)*drF(k)*hFacC(i,j,k,bi,bj)
299         &                 + ( oneRS - maskC(i,j,k,bi,bj) )
300    #endif
301    #ifdef ALLOW_OBCS
302             maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
303             maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
304    #else /* ALLOW_OBCS */
305             maskLocW(i,j) = _maskW(i,j,k,bi,bj)
306             maskLocS(i,j) = _maskS(i,j,k,bi,bj)
307    #endif /* ALLOW_OBCS */
308         ENDDO         ENDDO
309        ENDDO        ENDDO
310    
 #ifndef ALLOW_AUTODIFF_TAMC  
311        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
312          withSigns = .FALSE.          withSigns = .FALSE.
313          CALL FILL_CS_CORNER_UV_RS(          CALL FILL_CS_CORNER_UV_RS(
314       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )
315        ENDIF        ENDIF
 #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        interiorOnly = .FALSE.        interiorOnly = .FALSE.
330        overlapOnly  = .FALSE.        overlapOnly  = .FALSE.
331        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
 #ifdef MULTIDIM_OLD_VERSION  
 C-    CubedSphere : pass 3 times, with full update of local tracer field  
        IF (ipass.EQ.1) THEN  
         calc_fluxes_X = nCFace.EQ.1 .OR. nCFace.EQ.2  
         calc_fluxes_Y = nCFace.EQ.4 .OR. nCFace.EQ.5  
        ELSEIF (ipass.EQ.2) THEN  
         calc_fluxes_X = nCFace.EQ.3 .OR. nCFace.EQ.4  
         calc_fluxes_Y = nCFace.EQ.6 .OR. nCFace.EQ.1  
 #else /* MULTIDIM_OLD_VERSION */  
332  C-    CubedSphere : pass 3 times, with partial update of local tracer field  C-    CubedSphere : pass 3 times, with partial update of local tracer field
333         IF (ipass.EQ.1) THEN         IF (ipass.EQ.1) THEN
334          overlapOnly  = MOD(nCFace,3).EQ.0          overlapOnly  = MOD(nCFace,3).EQ.0
# Line 308  C-    CubedSphere : pass 3 times, with p Line 337  C-    CubedSphere : pass 3 times, with p
337          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
338         ELSEIF (ipass.EQ.2) THEN         ELSEIF (ipass.EQ.2) THEN
339          overlapOnly  = MOD(nCFace,3).EQ.2          overlapOnly  = MOD(nCFace,3).EQ.2
340            interiorOnly = MOD(nCFace,3).EQ.1
341          calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4          calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
342          calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1          calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
 #endif /* MULTIDIM_OLD_VERSION */  
343         ELSE         ELSE
344            interiorOnly = .TRUE.
345          calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6          calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
346          calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3          calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
347         ENDIF         ENDIF
# Line 323  C-    not CubedSphere Line 353  C-    not CubedSphere
353    
354  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
355  C--   X direction  C--   X direction
356  C-     Advective flux in X  
357          DO j=1-Oly,sNy+Oly  #ifdef ALLOW_AUTODIFF
358           DO i=1-Olx,sNx+Olx  C-     Always reset advective flux in X
359            DO j=1-OLy,sNy+OLy
360             DO i=1-OLx,sNx+OLx
361            af(i,j) = 0.            af(i,j) = 0.
362           ENDDO           ENDDO
363          ENDDO          ENDDO
 C  
 #ifdef ALLOW_AUTODIFF_TAMC  
364  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
365  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
366  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
367  CADJ STORE af(:,:)  =  CADJ STORE af(:,:)  =
368  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
369  # endif  # endif
370  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
371  C  
372        IF (calc_fluxes_X) THEN        IF (calc_fluxes_X) THEN
373    
374  C-     Do not compute fluxes if  C-     Do not compute fluxes if
# Line 346  C       a) needed in overlap only Line 376  C       a) needed in overlap only
376  C   and b) the overlap of myTile are not cube-face Edges  C   and b) the overlap of myTile are not cube-face Edges
377         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
378    
 #ifndef ALLOW_AUTODIFF_TAMC  
379  C-     Internal exchange for calculations in X  C-     Internal exchange for calculations in X
380  #ifdef MULTIDIM_OLD_VERSION          IF ( overlapOnly ) THEN
381          IF ( useCubedSphereExchange ) THEN           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
382  #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 )  
383          ENDIF          ENDIF
 #endif  
384    
385  #ifdef ALLOW_AUTODIFF_TAMC  C-     Advective flux in X
386    #ifndef ALLOW_AUTODIFF
387            DO j=1-OLy,sNy+OLy
388             DO i=1-OLx,sNx+OLx
389              af(i,j) = 0.
390             ENDDO
391            ENDDO
392    #else /* ALLOW_AUTODIFF */
393  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
394  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
395  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
396  # endif  # endif
397  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
398    
399          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
400       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
401            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_X(    bi,bj,k, advectionScheme, .TRUE.,
402       I                           dTtracerLev(k),uTrans,uFld,localTij,       I             deltaTLev(k),uTrans,uFld(1-OLx,1-OLy,k), localTij,
403       O                           af, myThid )       O             af, myThid )
404          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
405            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
406       I                              uTrans, uFld, maskLocW, localTij,       I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
407       O                              af, myThid )       O             af, myThid )
408          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
409            CALL GAD_DST3_ADV_X(      bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X(      bi,bj,k, .TRUE., deltaTLev(k),
410       I                              uTrans, uFld, maskLocW, localTij,       I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
411       O                              af, myThid )       O             af, myThid )
412          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
413            CALL GAD_DST3FL_ADV_X(    bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_X(    bi,bj,k, .TRUE., deltaTLev(k),
414       I                              uTrans, uFld, maskLocW, localTij,       I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
415       O                              af, myThid )       O             af, myThid )
416    #ifndef ALLOW_AUTODIFF
417            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
418              CALL GAD_OS7MP_ADV_X(     bi,bj,k, .TRUE., deltaTLev(k),
419         I             uTrans, uFld(1-OLx,1-OLy,k), maskLocW, localTij,
420         O             af, myThid )
421    #endif
422          ELSE          ELSE
423           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
424          ENDIF          ENDIF
425    
426  C-     Advective flux in X : done  #ifdef ALLOW_OBCS
427         ENDIF          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    
 #ifndef ALLOW_AUTODIFF_TAMC  
435  C-     Internal exchange for next calculations in Y  C-     Internal exchange for next calculations in Y
436         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
437           CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )           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         ENDIF
 #endif  
443    
444  C-     Update the local tracer field where needed:  C-     Update the local tracer field where needed:
445    C      use "maksInC" to prevent updating tracer field in OB regions
446    #ifdef ALLOW_AUTODIFF_TAMC
447    # ifdef GAD_MULTIDIM_COMPRESSIBLE
448    CADJ STORE localVol(:,:)  =
449    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
450    CADJ STORE localTij(:,:)  =
451    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
452    # endif
453    #endif /* ALLOW_AUTODIFF_TAMC */
454    
455  C      update in overlap-Only  C      update in overlap-Only
456         IF ( overlapOnly ) THEN         IF ( overlapOnly ) THEN
457          iMinUpd = 1-Olx+1          iMinUpd = 1-OLx+1
458          iMaxUpd = sNx+Olx-1          iMaxUpd = sNx+OLx-1
459  C- notes: these 2 lines below have no real effect (because recip_hFac=0  C- notes: these 2 lines below have no real effect (because recip_hFac=0
460  C         in corner region) but safer to keep them.  C         in corner region) but safer to keep them.
461          IF ( W_edge ) iMinUpd = 1          IF ( W_edge ) iMinUpd = 1
462          IF ( E_edge ) iMaxUpd = sNx          IF ( E_edge ) iMaxUpd = sNx
463    
464          IF ( S_edge ) THEN          IF ( S_edge ) THEN
465           DO j=1-Oly,0           DO j=1-OLy,0
466            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
467             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*  #ifdef GAD_MULTIDIM_COMPRESSIBLE
468       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)             tmpTrac = localTij(i,j)*localVol(i,j)
469       &       *recip_rA(i,j,bi,bj)       &      -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)       &       *( af(i+1,j)-af(i,j)
481       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
482       &        )       &        )*maskInC(i,j,bi,bj)
483    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
484            ENDDO            ENDDO
485           ENDDO           ENDDO
486          ENDIF          ENDIF
487          IF ( N_edge ) THEN          IF ( N_edge ) THEN
488           DO j=sNy+1,sNy+Oly           DO j=sNy+1,sNy+OLy
489            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
490             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*  #ifdef GAD_MULTIDIM_COMPRESSIBLE
491       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)             tmpTrac = localTij(i,j)*localVol(i,j)
492       &       *recip_rA(i,j,bi,bj)       &      -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)       &       *( af(i+1,j)-af(i,j)
504       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
505       &        )       &        )*maskInC(i,j,bi,bj)
506    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
507            ENDDO            ENDDO
508           ENDDO           ENDDO
509          ENDIF          ENDIF
510    
511         ELSE         ELSE
512  C      do not only update the overlap  C      do not only update the overlap
513          jMinUpd = 1-Oly          jMinUpd = 1-OLy
514          jMaxUpd = sNy+Oly          jMaxUpd = sNy+OLy
515          IF ( interiorOnly .AND. S_edge ) jMinUpd = 1          IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
516          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
517          DO j=jMinUpd,jMaxUpd          DO j=jMinUpd,jMaxUpd
518           DO i=1-Olx+1,sNx+Olx-1           DO i=1-OLx+1,sNx+OLx-1
519             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*  #ifdef GAD_MULTIDIM_COMPRESSIBLE
520       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)             tmpTrac = localTij(i,j)*localVol(i,j)
521       &       *recip_rA(i,j,bi,bj)       &      -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)       &       *( af(i+1,j)-af(i,j)
533       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
534       &        )       &        )*maskInC(i,j,bi,bj)
535    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
536           ENDDO           ENDDO
537          ENDDO          ENDDO
538  C-      keep advective flux (for diagnostics)  C-      keep advective flux (for diagnostics)
539          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
540           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
541            afx(i,j) = af(i,j)            afx(i,j) = af(i,j)
542           ENDDO           ENDDO
543          ENDDO          ENDDO
544    
 #ifdef ALLOW_OBCS  
 C-     Apply open boundary conditions  
         IF ( useOBCS ) THEN  
          IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN  
           CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )  
          ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN  
           CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )  
 #ifdef ALLOW_PTRACERS  
          ELSEIF (tracerIdentity.GE.GAD_TR1) THEN  
           CALL OBCS_APPLY_PTRACER( bi, bj, k,  
      &         tracerIdentity-GAD_TR1+1, localTij, myThid )  
 #endif /* ALLOW_PTRACERS */  
          ENDIF  
         ENDIF  
 #endif /* ALLOW_OBCS */  
   
545  C-     end if/else update overlap-Only  C-     end if/else update overlap-Only
546         ENDIF         ENDIF
547    
# Line 479  C--   End of X direction Line 550  C--   End of X direction
550    
551  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
552  C--   Y direction  C--   Y direction
553  cph-test  
554  C-     Advective flux in Y  #ifdef ALLOW_AUTODIFF
555          DO j=1-Oly,sNy+Oly  C-     Always reset advective flux in Y
556           DO i=1-Olx,sNx+Olx          DO j=1-OLy,sNy+OLy
557             DO i=1-OLx,sNx+OLx
558            af(i,j) = 0.            af(i,j) = 0.
559           ENDDO           ENDDO
560          ENDDO          ENDDO
 C  
 #ifdef ALLOW_AUTODIFF_TAMC  
561  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
562  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
563  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
564  CADJ STORE af(:,:)  =  CADJ STORE af(:,:)  =
565  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
566  # endif  # endif
567  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
568  C  
569        IF (calc_fluxes_Y) THEN        IF (calc_fluxes_Y) THEN
570    
571  C-     Do not compute fluxes if  C-     Do not compute fluxes if
# Line 503  C       a) needed in overlap only Line 573  C       a) needed in overlap only
573  C   and b) the overlap of myTile are not cube-face edges  C   and b) the overlap of myTile are not cube-face edges
574         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
575    
 #ifndef ALLOW_AUTODIFF_TAMC  
576  C-     Internal exchange for calculations in Y  C-     Internal exchange for calculations in Y
577  #ifdef MULTIDIM_OLD_VERSION          IF ( overlapOnly ) THEN
578          IF ( useCubedSphereExchange ) THEN           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
579  #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 )  
580          ENDIF          ENDIF
 #endif  
581    
582  C-     Advective flux in Y  C-     Advective flux in Y
583          DO j=1-Oly,sNy+Oly  #ifndef ALLOW_AUTODIFF
584           DO i=1-Olx,sNx+Olx          DO j=1-OLy,sNy+OLy
585             DO i=1-OLx,sNx+OLx
586            af(i,j) = 0.            af(i,j) = 0.
587           ENDDO           ENDDO
588          ENDDO          ENDDO
589    #else /* ALLOW_AUTODIFF */
 #ifdef ALLOW_AUTODIFF_TAMC  
590  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
591  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
592  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
593  #endif  #endif
594  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
595    
596          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
597       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
598            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_Y(    bi,bj,k, advectionScheme, .TRUE.,
599       I                           dTtracerLev(k),vTrans,vFld,localTij,       I             deltaTLev(k),vTrans,vFld(1-OLx,1-OLy,k), localTij,
600       O                           af, myThid )       O             af, myThid )
601          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
602            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
603       I                              vTrans, vFld, maskLocS, localTij,       I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
604       O                              af, myThid )       O             af, myThid )
605          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
606            CALL GAD_DST3_ADV_Y(      bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y(      bi,bj,k, .TRUE., deltaTLev(k),
607       I                              vTrans, vFld, maskLocS, localTij,       I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
608       O                              af, myThid )       O             af, myThid )
609          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
610            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .TRUE., deltaTLev(k),
611       I                              vTrans, vFld, maskLocS, localTij,       I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
612       O                              af, myThid )       O             af, myThid )
613    #ifndef ALLOW_AUTODIFF
614            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
615              CALL GAD_OS7MP_ADV_Y(     bi,bj,k, .TRUE., deltaTLev(k),
616         I             vTrans, vFld(1-OLx,1-OLy,k), maskLocS, localTij,
617         O             af, myThid )
618    #endif
619          ELSE          ELSE
620           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
621          ENDIF          ENDIF
622    
623  C-     Advective flux in Y : done  #ifdef ALLOW_OBCS
624         ENDIF          IF ( useOBCS ) THEN
625    C-      replace advective flux with 1st order upwind scheme estimate
626              CALL OBCS_U1_ADV_TRACER( .FALSE., trIdentity, bi, bj, k,
627         I                             maskLocS, vTrans, localTij,
628         U                             af, myThid )
629            ENDIF
630    #endif /* ALLOW_OBCS */
631    
 #ifndef ALLOW_AUTODIFF_TAMC  
632  C-     Internal exchange for next calculations in X  C-     Internal exchange for next calculations in X
633         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
634           CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )           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         ENDIF
 #endif  
640    
641  C-     Update the local tracer field where needed:  C-     Update the local tracer field where needed:
642    C      use "maksInC" to prevent updating tracer field in OB regions
643    #ifdef ALLOW_AUTODIFF_TAMC
644    # ifdef GAD_MULTIDIM_COMPRESSIBLE
645    CADJ STORE localVol(:,:)  =
646    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
647    CADJ STORE localTij(:,:)  =
648    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
649    # endif
650    #endif /* ALLOW_AUTODIFF_TAMC */
651    
652  C      update in overlap-Only  C      update in overlap-Only
653         IF ( overlapOnly ) THEN         IF ( overlapOnly ) THEN
654          jMinUpd = 1-Oly+1          jMinUpd = 1-OLy+1
655          jMaxUpd = sNy+Oly-1          jMaxUpd = sNy+OLy-1
656  C- notes: these 2 lines below have no real effect (because recip_hFac=0  C- notes: these 2 lines below have no real effect (because recip_hFac=0
657  C         in corner region) but safer to keep them.  C         in corner region) but safer to keep them.
658          IF ( S_edge ) jMinUpd = 1          IF ( S_edge ) jMinUpd = 1
# Line 573  C         in corner region) but safer to Line 660  C         in corner region) but safer to
660    
661          IF ( W_edge ) THEN          IF ( W_edge ) THEN
662           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
663            DO i=1-Olx,0            DO i=1-OLx,0
664             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*  #ifdef GAD_MULTIDIM_COMPRESSIBLE
665       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)             tmpTrac = localTij(i,j)*localVol(i,j)
666       &       *recip_rA(i,j,bi,bj)       &      -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)       &       *( af(i,j+1)-af(i,j)
678       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
679       &        )       &        )*maskInC(i,j,bi,bj)
680    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
681            ENDDO            ENDDO
682           ENDDO           ENDDO
683          ENDIF          ENDIF
684          IF ( E_edge ) THEN          IF ( E_edge ) THEN
685           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
686            DO i=sNx+1,sNx+Olx            DO i=sNx+1,sNx+OLx
687             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*  #ifdef GAD_MULTIDIM_COMPRESSIBLE
688       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)             tmpTrac = localTij(i,j)*localVol(i,j)
689       &       *recip_rA(i,j,bi,bj)       &      -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)       &       *( af(i,j+1)-af(i,j)
701       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
702       &        )       &        )*maskInC(i,j,bi,bj)
703    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
704            ENDDO            ENDDO
705           ENDDO           ENDDO
706          ENDIF          ENDIF
707    
708         ELSE         ELSE
709  C      do not only update the overlap  C      do not only update the overlap
710          iMinUpd = 1-Olx          iMinUpd = 1-OLx
711          iMaxUpd = sNx+Olx          iMaxUpd = sNx+OLx
712          IF ( interiorOnly .AND. W_edge ) iMinUpd = 1          IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
713          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
714          DO j=1-Oly+1,sNy+Oly-1          DO j=1-OLy+1,sNy+OLy-1
715           DO i=iMinUpd,iMaxUpd           DO i=iMinUpd,iMaxUpd
716             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*  #ifdef GAD_MULTIDIM_COMPRESSIBLE
717       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)             tmpTrac = localTij(i,j)*localVol(i,j)
718       &       *recip_rA(i,j,bi,bj)       &      -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)       &       *( af(i,j+1)-af(i,j)
730       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
731       &        )       &        )*maskInC(i,j,bi,bj)
732    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
733           ENDDO           ENDDO
734          ENDDO          ENDDO
735  C-      keep advective flux (for diagnostics)  C-      keep advective flux (for diagnostics)
736          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
737           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
738            afy(i,j) = af(i,j)            afy(i,j) = af(i,j)
739           ENDDO           ENDDO
740          ENDDO          ENDDO
741    
 #ifdef ALLOW_OBCS  
 C-     Apply open boundary conditions  
         IF (useOBCS) THEN  
          IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN  
           CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )  
          ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN  
           CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )  
 #ifdef ALLOW_PTRACERS  
          ELSEIF (tracerIdentity.GE.GAD_TR1) THEN  
           CALL OBCS_APPLY_PTRACER( bi, bj, k,  
      &         tracerIdentity-GAD_TR1+1, localTij, myThid )  
 #endif /* ALLOW_PTRACERS */  
          ENDIF  
         ENDIF  
 #endif /* ALLOW_OBCS */  
   
742  C      end if/else update overlap-Only  C      end if/else update overlap-Only
743         ENDIF         ENDIF
744    
# Line 646  C--   End of ipass loop Line 750  C--   End of ipass loop
750    
751        IF ( implicitAdvection ) THEN        IF ( implicitAdvection ) THEN
752  C-    explicit advection is done ; store tendency in gTracer:  C-    explicit advection is done ; store tendency in gTracer:
753          DO j=1-Oly,sNy+Oly  #ifdef GAD_MULTIDIM_COMPRESSIBLE
754           DO i=1-Olx,sNx+Olx          STOP 'GAD_ADVECTION: missing code for implicitAdvection'
755    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
756            DO j=1-OLy,sNy+OLy
757             DO i=1-OLx,sNx+OLx
758            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
759       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
760           ENDDO           ENDDO
761          ENDDO          ENDDO
762        ELSE        ELSE
763  C-    horizontal advection done; store intermediate result in 3D array:  C-    horizontal advection done; store intermediate result in 3D array:
764         DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
765          DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
766           localTijk(i,j,k)=localTij(i,j)  #ifdef GAD_MULTIDIM_COMPRESSIBLE
767              locVol3d(i,j,k) = localVol(i,j)
768    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
769              localT3d(i,j,k) = localTij(i,j)
770             ENDDO
771          ENDDO          ENDDO
        ENDDO  
772        ENDIF        ENDIF
773    
774  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
775          IF ( useDiagnostics ) THEN          IF ( doDiagAdvX ) THEN
776            diagName = 'ADVx'//diagSufx            diagName = 'ADVx'//diagSufx
777            CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( afx, diagName, k,1, 2,bi,bj, myThid )
778            ENDIF
779            IF ( doDiagAdvY ) THEN
780            diagName = 'ADVy'//diagSufx            diagName = 'ADVy'//diagSufx
781            CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( afy, diagName, k,1, 2,bi,bj, myThid )
782          ENDIF          ENDIF
783  #endif  #endif /* ALLOW_DIAGNOSTICS */
784    
785  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
786        IF ( debugLevel .GE. debLevB        IF ( debugLevel .GE. debLevC
787       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE       &   .AND. trIdentity.EQ.GAD_TEMPERATURE
788       &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0       &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
789       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
790       &   .AND. useCubedSphereExchange ) THEN       &   .AND. useCubedSphereExchange ) THEN
# Line 689  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)
            wFld(i,j)   = 0.  
833             rTrans(i,j) = 0.             rTrans(i,j) = 0.
834             fVerT(i,j,kUp) = 0.             fVerT(i,j,kUp) = 0.
835            ENDDO            ENDDO
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             wFld(i,j)   = wVel(i,j,k,bi,bj)             rTransKp(i,j) = kp1Msk*rTrans(i,j)
844             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)             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       U                 wFld, rTrans,  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
      I                 k, bi, bj, myThid )  
 #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 ( advectionScheme.EQ.ENUM_UPWIND_1RST           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
860       &      .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
861             CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_R(    bi,bj,k, advectionScheme,
862       I                            dTtracerLev(k),rTrans,wFld,localTijk,       I              deltaTLev(k),rTrans,wFld(1-OLx,1-OLy,k),localT3d,
863       O                            fVerT(1-Olx,1-Oly,kUp), myThid )       O              fVerT(1-OLx,1-OLy,kUp), myThid )
864           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
865             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k),             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
866       I                               rTrans, wFld, localTijk,       I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
867       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O              fVerT(1-OLx,1-OLy,kUp), myThid )
868           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
869             CALL GAD_DST3_ADV_R(      bi,bj,k, dTtracerLev(k),             CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTLev(k),
870       I                               rTrans, wFld, localTijk,       I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
871       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O              fVerT(1-OLx,1-OLy,kUp), myThid )
872           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
873             CALL GAD_DST3FL_ADV_R(    bi,bj,k, dTtracerLev(k),             CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTLev(k),
874       I                               rTrans, wFld, localTijk,       I              rTrans, wFld(1-OLx,1-OLy,k), localT3d,
875       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O              fVerT(1-OLx,1-OLy,kUp), myThid )
876    #ifndef ALLOW_AUTODIFF
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 773  C-    Compute vertical advective flux in Line 886  C-    Compute vertical advective flux in
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)-dTtracerLev(k)*           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,kDown)-fVerT(i,j,kUp)       &                   *rkSign*maskInC(i,j,bi,bj)
909       &       -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))            localVol(i,j) = locVol3d(i,j,k)
910       &      )*rkSign       &      -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))/dTtracerLev(k)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
941           ENDDO           ENDDO
942          ENDDO          ENDDO
943    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
944    
945  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
946          IF ( useDiagnostics ) THEN          IF ( doDiagAdvR ) THEN
947            diagName = 'ADVr'//diagSufx            diagName = 'ADVr'//diagSufx
948            CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp),            CALL DIAGNOSTICS_FILL( fVerT(1-OLx,1-OLy,kUp),
949       &                           diagName, k,1, 2,bi,bj, myThid)       &                           diagName, k,1, 2,bi,bj, myThid )
950          ENDIF          ENDIF
951  #endif  #endif /* ALLOW_DIAGNOSTICS */
952    
953  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
954         ENDDO         ENDDO
Legend:
Removed from v.1.41  
changed lines
  Added in v.1.74
  ViewVC Help
Powered by ViewVC 1.1.22