/[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.38 by jmc, Sat Nov 5 01:01:51 2005 UTC revision 1.64 by jmc, Tue Mar 29 15:52:57 2011 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  #undef OBCS_MULTIDIM_OLD_VERSION
6    
7  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
8  CBOP  CBOP
# Line 11  C !ROUTINE: GAD_ADVECTION Line 11  C !ROUTINE: GAD_ADVECTION
11  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
12        SUBROUTINE GAD_ADVECTION(        SUBROUTINE GAD_ADVECTION(
13       I     implicitAdvection, advectionScheme, vertAdvecScheme,       I     implicitAdvection, advectionScheme, vertAdvecScheme,
14       I     tracerIdentity,       I     tracerIdentity, deltaTLev,
15       I     uVel, vVel, wVel, tracer,       I     uVel, vVel, wVel, tracer,
16       O     gTracer,       O     gTracer,
17       I     bi,bj, myTime,myIter,myThid)       I     bi,bj, myTime,myIter,myThid)
18    
19  C !DESCRIPTION:  C !DESCRIPTION:
20  C Calculates the tendancy of a tracer due to advection.  C Calculates the tendency of a tracer due to advection.
21  C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}  C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
22  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
23  C direct-space-time method and flux-limiters.  C direct-space-time method and flux-limiters.
24  C  C
25  C The algorithm is as follows:  C The algorithm is as follows:
26  C \begin{itemize}  C \begin{itemize}
# Line 33  C      - \Delta t \partial_r (w\theta^{( Line 33  C      - \Delta t \partial_r (w\theta^{(
33  C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}  C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}
34  C \end{itemize}  C \end{itemize}
35  C  C
36  C The tendancy (output) is over-written by this routine.  C The tendency (output) is over-written by this routine.
37    
38  C !USES: ===============================================================  C !USES: ===============================================================
39        IMPLICIT NONE        IMPLICIT NONE
# Line 50  C !USES: =============================== Line 50  C !USES: ===============================
50  # endif  # endif
51  #endif  #endif
52  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
53    #include "W2_EXCH2_SIZE.h"
54  #include "W2_EXCH2_TOPOLOGY.h"  #include "W2_EXCH2_TOPOLOGY.h"
 #include "W2_EXCH2_PARAMS.h"  
55  #endif /* ALLOW_EXCH2 */  #endif /* ALLOW_EXCH2 */
56    
57  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
# Line 70  C  myThid            :: thread number Line 70  C  myThid            :: thread number
70        LOGICAL implicitAdvection        LOGICAL implicitAdvection
71        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
72        INTEGER tracerIdentity        INTEGER tracerIdentity
73          _RL deltaTLev(Nr)
74        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
75        _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
76        _RL wVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL wVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
# Line 80  C  myThid            :: thread number Line 81  C  myThid            :: thread number
81        INTEGER myThid        INTEGER myThid
82    
83  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
84  C  gTracer           :: tendancy array  C  gTracer           :: tendency array
85        _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)
86    
87  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
88  C  maskUp        :: 2-D array for mask at W points  C  maskUp        :: 2-D array for mask at W points
89  C  maskLocW      :: 2-D array for mask at West points  C  maskLocW      :: 2-D array for mask at West points
90  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  
91  C [iMin,iMax]Upd :: loop range to update tracer field  C [iMin,iMax]Upd :: loop range to update tracer field
92  C [jMin,jMax]Upd :: loop range to update tracer field  C [jMin,jMax]Upd :: loop range to update tracer field
93  C  i,j,k         :: loop indices  C  i,j,k         :: loop indices
94  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
95  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
96  C  kp1           :: =k+1 for k<Nr, =Nr for k=Nr  C  kp1           :: =k+1 for k<Nr, =Nr for k=Nr
97  C  xA,yA         :: areas of X and Y face of tracer cells  C  xA,yA         :: areas of X and Y face of tracer cells
98    C  uFld,vFld     :: 2-D local copy of horizontal velocity, U,V components
99    C  wFld          :: 2-D local copy of vertical velocity
100  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
101  C  rTrans        :: 2-D arrays of volume transports at W points  C  rTrans        :: 2-D arrays of volume transports at W points
102  C  rTransKp1     :: vertical volume transport at interface k+1  C  rTransKp1     :: vertical volume transport at interface k+1
# Line 110  C  calc_fluxes_X :: logical to indicate Line 111  C  calc_fluxes_X :: logical to indicate
111  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
112  C  interiorOnly  :: only update the interior of myTile, but not the edges  C  interiorOnly  :: only update the interior of myTile, but not the edges
113  C  overlapOnly   :: only update the edges of myTile, but not the interior  C  overlapOnly   :: only update the edges of myTile, but not the interior
114  C  nipass        :: number of passes in multi-dimensional method  C  npass         :: number of passes in multi-dimensional method
115  C  ipass         :: number of the current pass being made  C  ipass         :: number of the current pass being made
116  C  myTile        :: variables used to determine which cube face  C  myTile        :: variables used to determine which cube face
117  C  nCFace        :: owns a tile for cube grid runs using  C  nCFace        :: owns a tile for cube grid runs using
118  C                :: multi-dim advection.  C                :: multi-dim advection.
119  C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube  C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
120        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  c     _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121        _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _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  
123        INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd        INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
124        INTEGER i,j,k,kup,kDown        INTEGER i,j,k,kUp,kDown
125        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127          _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128          _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129          _RL wFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 137  C [N,S,E,W]_edge :: true if N,S,E,W edge Line 140  C [N,S,E,W]_edge :: true if N,S,E,W edge
140        _RL kp1Msk        _RL kp1Msk
141        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
142        LOGICAL interiorOnly, overlapOnly        LOGICAL interiorOnly, overlapOnly
143        INTEGER nipass,ipass        INTEGER npass, ipass
144        INTEGER nCFace        INTEGER nCFace
145        LOGICAL N_edge, S_edge, E_edge, W_edge        LOGICAL N_edge, S_edge, E_edge, W_edge
146  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
# Line 145  C [N,S,E,W]_edge :: true if N,S,E,W edge Line 148  C [N,S,E,W]_edge :: true if N,S,E,W edge
148  #endif  #endif
149  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
150        CHARACTER*8 diagName        CHARACTER*8 diagName
151        CHARACTER*4 GAD_DIAG_SUFX, diagSufx        CHARACTER*4 diagSufx
152          LOGICAL     doDiagAdvX, doDiagAdvY, doDiagAdvR
153    C-    Functions:
154          CHARACTER*4 GAD_DIAG_SUFX
155        EXTERNAL    GAD_DIAG_SUFX        EXTERNAL    GAD_DIAG_SUFX
156          LOGICAL  DIAGNOSTICS_IS_ON
157          EXTERNAL DIAGNOSTICS_IS_ON
158  #endif  #endif
159  CEOP  CEOP
160    
161  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
162            act0 = tracerIdentity - 1            act0 = tracerIdentity
163            max0 = maxpass            max0 = maxpass
164            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
165            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
# Line 160  CEOP Line 168  CEOP
168            act3 = myThid - 1            act3 = myThid - 1
169            max3 = nTx*nTy            max3 = nTx*nTy
170            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
171            igadkey = (act0 + 1)            igadkey = act0
172       &                      + act1*max0       &                      + act1*max0
173       &                      + act2*max0*max1       &                      + act2*max0*max1
174       &                      + act3*max0*max1*max2       &                      + act3*max0*max1*max2
175       &                      + act4*max0*max1*max2*max3       &                      + act4*max0*max1*max2*max3
176            if (tracerIdentity.GT.maxpass) then            IF (tracerIdentity.GT.maxpass) THEN
177               print *, 'ph-pass gad_advection ', maxpass, tracerIdentity               print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
178               STOP 'maxpass seems smaller than tracerIdentity'               STOP 'maxpass seems smaller than tracerIdentity'
179            endif            ENDIF
180  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
181    
182  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
183  C--   Set diagnostic suffix for the current tracer  C--   Set diagnostics flags and suffix for the current tracer
184          doDiagAdvX = .FALSE.
185          doDiagAdvY = .FALSE.
186          doDiagAdvR = .FALSE.
187        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
188          diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )          diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
189            diagName = 'ADVx'//diagSufx
190            doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )
191            diagName = 'ADVy'//diagSufx
192            doDiagAdvY = DIAGNOSTICS_IS_ON( diagName, myThid )
193            diagName = 'ADVr'//diagSufx
194            doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )
195        ENDIF        ENDIF
196  #endif  #endif
197    
# Line 193  C     uninitialised but inert locations. Line 210  C     uninitialised but inert locations.
210          fVerT(i,j,1) = 0. _d 0          fVerT(i,j,1) = 0. _d 0
211          fVerT(i,j,2) = 0. _d 0          fVerT(i,j,2) = 0. _d 0
212          rTransKp1(i,j)= 0. _d 0          rTransKp1(i,j)= 0. _d 0
213    #ifdef ALLOW_AUTODIFF_TAMC
214            localTij(i,j) = 0. _d 0
215            wfld(i,j)    = 0. _d 0
216    #endif
217         ENDDO         ENDDO
218        ENDDO        ENDDO
219    
220  C--   Set tile-specific parameters for horizontal fluxes  C--   Set tile-specific parameters for horizontal fluxes
221        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
222         nipass=3         npass  = 3
223  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
224         IF ( nipass.GT.maxcube ) STOP 'maxcube needs to be = 3'         IF ( npass.GT.maxcube ) STOP 'maxcube needs to be = 3'
225  #endif  #endif
226  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
227         myTile = W2_myTileList(bi)         myTile = W2_myTileList(bi,bj)
228         nCFace = exch2_myFace(myTile)         nCFace = exch2_myFace(myTile)
229         N_edge = exch2_isNedge(myTile).EQ.1         N_edge = exch2_isNedge(myTile).EQ.1
230         S_edge = exch2_isSedge(myTile).EQ.1         S_edge = exch2_isSedge(myTile).EQ.1
# Line 217  C--   Set tile-specific parameters for h Line 238  C--   Set tile-specific parameters for h
238         W_edge = .TRUE.         W_edge = .TRUE.
239  #endif  #endif
240        ELSE        ELSE
241         nipass=2         npass  = 2
242         nCFace = bi         nCFace = 0
243         N_edge = .FALSE.         N_edge = .FALSE.
244         S_edge = .FALSE.         S_edge = .FALSE.
245         E_edge = .FALSE.         E_edge = .FALSE.
246         W_edge = .FALSE.         W_edge = .FALSE.
247        ENDIF        ENDIF
248    
       iMin = 1-OLx  
       iMax = sNx+OLx  
       jMin = 1-OLy  
       jMax = sNy+OLy  
   
249  C--   Start of k loop for horizontal fluxes  C--   Start of k loop for horizontal fluxes
250        DO k=1,Nr        DO k=1,Nr
251  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
252           kkey = (igadkey-1)*Nr + k           kkey = (igadkey-1)*Nr + k
253  CADJ STORE tracer(:,:,k,bi,bj) =  CADJ STORE tracer(:,:,k,bi,bj) =
254  CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
255  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
256    
257  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
258        CALL CALC_COMMON_FACTORS (        CALL CALC_COMMON_FACTORS (
259       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         uVel, vVel,
260       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         uFld, vFld, uTrans, vTrans, xA, yA,
261       I         myThid)       I         k,bi,bj, myThid )
262    
263  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
264  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
265        IF (useGMRedi)        IF (useGMRedi)
266       &   CALL GMREDI_CALC_UVFLOW(       &   CALL GMREDI_CALC_UVFLOW(
267       &            uTrans, vTrans, bi, bj, k, myThid)       U                  uFld, vFld, uTrans, vTrans,
268         I                  k, bi, bj, myThid )
269  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
270    
271  C--   Make local copy of tracer array and mask West & South  C--   Make local copy of tracer array and mask West & South
272        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
273         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
274           localTij(i,j)=tracer(i,j,k,bi,bj)           localTij(i,j)=tracer(i,j,k,bi,bj)
275           maskLocW(i,j)=maskW(i,j,k,bi,bj)  #ifdef ALLOW_OBCS
276           maskLocS(i,j)=maskS(i,j,k,bi,bj)           maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
277             maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
278    #else /* ALLOW_OBCS */
279             maskLocW(i,j) = _maskW(i,j,k,bi,bj)
280             maskLocS(i,j) = _maskS(i,j,k,bi,bj)
281    #endif /* ALLOW_OBCS */
282         ENDDO         ENDDO
283        ENDDO        ENDDO
284    
285  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
286        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
287          withSigns = .FALSE.          withSigns = .FALSE.
288          CALL FILL_CS_CORNER_UV_RS(          CALL FILL_CS_CORNER_UV_RS(
289       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )
290        ENDIF        ENDIF
291  #endif  cph-exch2#endif
292    
293  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
294  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.
295        DO ipass=1,nipass        DO ipass=1,npass
296  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
297           passkey = ipass + (k-1)      *maxcube           passkey = ipass
298       &                   + (igadkey-1)*maxcube*Nr       &                   + (k-1)      *maxpass
299           IF (nipass .GT. maxpass) THEN       &                   + (igadkey-1)*maxpass*Nr
300            STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'           IF (npass .GT. maxpass) THEN
301              STOP 'GAD_ADVECTION: npass > maxcube. check tamc.h'
302           ENDIF           ENDIF
303  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
304    
305        interiorOnly = .FALSE.        interiorOnly = .FALSE.
306        overlapOnly  = .FALSE.        overlapOnly  = .FALSE.
307        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 */  
308  C-    CubedSphere : pass 3 times, with partial update of local tracer field  C-    CubedSphere : pass 3 times, with partial update of local tracer field
309         IF (ipass.EQ.1) THEN         IF (ipass.EQ.1) THEN
310          overlapOnly  = MOD(nCFace,3).EQ.0          overlapOnly  = MOD(nCFace,3).EQ.0
# Line 299  C-    CubedSphere : pass 3 times, with p Line 313  C-    CubedSphere : pass 3 times, with p
313          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
314         ELSEIF (ipass.EQ.2) THEN         ELSEIF (ipass.EQ.2) THEN
315          overlapOnly  = MOD(nCFace,3).EQ.2          overlapOnly  = MOD(nCFace,3).EQ.2
316            interiorOnly = MOD(nCFace,3).EQ.1
317          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
318          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 */  
319         ELSE         ELSE
320            interiorOnly = .TRUE.
321          calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6          calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
322          calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3          calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
323         ENDIF         ENDIF
# Line 311  C-    not CubedSphere Line 326  C-    not CubedSphere
326          calc_fluxes_X = MOD(ipass,2).EQ.1          calc_fluxes_X = MOD(ipass,2).EQ.1
327          calc_fluxes_Y = .NOT.calc_fluxes_X          calc_fluxes_Y = .NOT.calc_fluxes_X
328        ENDIF        ENDIF
329    
330  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
331  C--   X direction  C--   X direction
332    C-     Advective flux in X
333            DO j=1-Oly,sNy+Oly
334             DO i=1-Olx,sNx+Olx
335              af(i,j) = 0.
336             ENDDO
337            ENDDO
338    C
339    #ifdef ALLOW_AUTODIFF_TAMC
340    # ifndef DISABLE_MULTIDIM_ADVECTION
341    CADJ STORE localTij(:,:)  =
342    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
343    CADJ STORE af(:,:)  =
344    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
345    # endif
346    #endif /* ALLOW_AUTODIFF_TAMC */
347    C
348        IF (calc_fluxes_X) THEN        IF (calc_fluxes_X) THEN
349    
350  C-     Do not compute fluxes if  C-     Do not compute fluxes if
351  C       a) needed in overlap only  C       a) needed in overlap only
352  C   and b) the overlap of myTile are not cube-face Edges  C   and b) the overlap of myTile are not cube-face Edges
353         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
354    
 #ifndef ALLOW_AUTODIFF_TAMC  
355  C-     Internal exchange for calculations in X  C-     Internal exchange for calculations in X
356  #ifdef MULTIDIM_OLD_VERSION          IF ( overlapOnly ) THEN
357          IF ( useCubedSphereExchange ) THEN           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
358  #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 )  
359          ENDIF          ENDIF
 #endif  
   
 C-     Advective flux in X  
         DO j=1-Oly,sNy+Oly  
          DO i=1-Olx,sNx+Olx  
           af(i,j) = 0.  
          ENDDO  
         ENDDO  
360    
361  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
362  #ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
363  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
364  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
365  #endif  # endif
366  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
367    
368          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
369       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
370            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
371       I                           dTtracerLev(k),uTrans,uVel,localTij,       I                           deltaTLev(k),uTrans,uFld,localTij,
372       O                           af, myThid )       O                           af, myThid )
373          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
374            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
375       I                              uTrans, uVel, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
376       O                              af, myThid )       O                              af, myThid )
377          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
378            CALL GAD_DST3_ADV_X(      bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X(      bi,bj,k, .TRUE., deltaTLev(k),
379       I                              uTrans, uVel, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
380       O                              af, myThid )       O                              af, myThid )
381          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
382            CALL GAD_DST3FL_ADV_X(    bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_X(    bi,bj,k, .TRUE., deltaTLev(k),
383       I                              uTrans, uVel, maskLocW, localTij,       I                              uTrans, uFld, maskLocW, localTij,
384       O                              af, myThid )       O                              af, myThid )
385    #ifndef ALLOW_AUTODIFF_TAMC
386            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
387              CALL GAD_OS7MP_ADV_X(     bi,bj,k, .TRUE., deltaTLev(k),
388         I                              uTrans, uFld, maskLocW, localTij,
389         O                              af, myThid )
390    #endif
391          ELSE          ELSE
392           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
393          ENDIF          ENDIF
394    
 C-     Advective flux in X : done  
        ENDIF  
   
 #ifndef ALLOW_AUTODIFF_TAMC  
395  C-     Internal exchange for next calculations in Y  C-     Internal exchange for next calculations in Y
396         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
397           CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
398         &                              localTij, bi,bj, myThid )
399            ENDIF
400    
401    C-     Advective flux in X : done
402         ENDIF         ENDIF
 #endif  
403    
404  C-     Update the local tracer field where needed:  C-     Update the local tracer field where needed:
405    C      use "maksInC" to prevent updating tracer field in OB regions
406    
407  C      update in overlap-Only  C      update in overlap-Only
408         IF ( overlapOnly ) THEN         IF ( overlapOnly ) THEN
409          iMinUpd = 1-Olx+1          iMinUpd = 1-Olx+1
410          iMaxUpd = sNx+Olx-1          iMaxUpd = sNx+Olx-1
411  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
412  C         in corner region) but safer to keep them.  C         in corner region) but safer to keep them.
413          IF ( W_edge ) iMinUpd = 1          IF ( W_edge ) iMinUpd = 1
414          IF ( E_edge ) iMaxUpd = sNx          IF ( E_edge ) iMaxUpd = sNx
# Line 392  C         in corner region) but safer to Line 416  C         in corner region) but safer to
416          IF ( S_edge ) THEN          IF ( S_edge ) THEN
417           DO j=1-Oly,0           DO j=1-Oly,0
418            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
419             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
420       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
421       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
422         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
423       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
424       &         -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))
425       &        )       &        )*maskInC(i,j,bi,bj)
426            ENDDO            ENDDO
427           ENDDO           ENDDO
428          ENDIF          ENDIF
429          IF ( N_edge ) THEN          IF ( N_edge ) THEN
430           DO j=sNy+1,sNy+Oly           DO j=sNy+1,sNy+Oly
431            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
432             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
433       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
434       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
435         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
436       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
437       &         -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))
438       &        )       &        )*maskInC(i,j,bi,bj)
439            ENDDO            ENDDO
440           ENDDO           ENDDO
441          ENDIF          ENDIF
442    
443         ELSE         ELSE
444  C      do not only update the overlap  C      do not only update the overlap
445          jMinUpd = 1-Oly          jMinUpd = 1-Oly
446          jMaxUpd = sNy+Oly          jMaxUpd = sNy+Oly
447          IF ( interiorOnly .AND. S_edge ) jMinUpd = 1          IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
448          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
449          DO j=jMinUpd,jMaxUpd          DO j=jMinUpd,jMaxUpd
450           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
451             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
452       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
453       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
454         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
455       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
456       &         -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))
457       &        )       &        )*maskInC(i,j,bi,bj)
458           ENDDO           ENDDO
459          ENDDO          ENDDO
460  C-      keep advective flux (for diagnostics)  C-      keep advective flux (for diagnostics)
# Line 438  C-      keep advective flux (for diagnos Line 465  C-      keep advective flux (for diagnos
465          ENDDO          ENDDO
466    
467  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
468    #ifdef OBCS_MULTIDIM_OLD_VERSION
469  C-     Apply open boundary conditions  C-     Apply open boundary conditions
470          IF ( useOBCS ) THEN          IF ( useOBCS ) THEN
471           IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN           IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
# Line 446  C-     Apply open boundary conditions Line 474  C-     Apply open boundary conditions
474            CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )            CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
475  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
476           ELSEIF (tracerIdentity.GE.GAD_TR1) THEN           ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
477            CALL OBCS_APPLY_PTRACER( bi, bj, k,            CALL OBCS_APPLY_PTRACER( bi, bj, k,
478       &         tracerIdentity-GAD_TR1+1, localTij, myThid )       &         tracerIdentity-GAD_TR1+1, localTij, myThid )
479  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
480           ENDIF           ENDIF
481          ENDIF          ENDIF
482    #endif /* OBCS_MULTIDIM_OLD_VERSION */
483  #endif /* ALLOW_OBCS */  #endif /* ALLOW_OBCS */
484    
485  C-     end if/else update overlap-Only  C-     end if/else update overlap-Only
486         ENDIF         ENDIF
487            
488  C--   End of X direction  C--   End of X direction
489        ENDIF        ENDIF
490    
491  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
492  C--   Y direction  C--   Y direction
493    cph-test
494    C-     Advective flux in Y
495            DO j=1-Oly,sNy+Oly
496             DO i=1-Olx,sNx+Olx
497              af(i,j) = 0.
498             ENDDO
499            ENDDO
500    C
501    #ifdef ALLOW_AUTODIFF_TAMC
502    # ifndef DISABLE_MULTIDIM_ADVECTION
503    CADJ STORE localTij(:,:)  =
504    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
505    CADJ STORE af(:,:)  =
506    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
507    # endif
508    #endif /* ALLOW_AUTODIFF_TAMC */
509    C
510        IF (calc_fluxes_Y) THEN        IF (calc_fluxes_Y) THEN
511    
512  C-     Do not compute fluxes if  C-     Do not compute fluxes if
# Line 468  C       a) needed in overlap only Line 514  C       a) needed in overlap only
514  C   and b) the overlap of myTile are not cube-face edges  C   and b) the overlap of myTile are not cube-face edges
515         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
516    
 #ifndef ALLOW_AUTODIFF_TAMC  
517  C-     Internal exchange for calculations in Y  C-     Internal exchange for calculations in Y
518  #ifdef MULTIDIM_OLD_VERSION          IF ( overlapOnly ) THEN
519          IF ( useCubedSphereExchange ) THEN           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
520  #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 )  
521          ENDIF          ENDIF
 #endif  
522    
523  C-     Advective flux in Y  C-     Advective flux in Y
524          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
# Line 487  C-     Advective flux in Y Line 527  C-     Advective flux in Y
527           ENDDO           ENDDO
528          ENDDO          ENDDO
529    
530  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
531  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
532  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
533  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
534  #endif  #endif
535  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
536    
537          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
538       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &     .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
539            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
540       I                           dTtracerLev(k),vTrans,vVel,localTij,       I                           deltaTLev(k),vTrans,vFld,localTij,
541       O                           af, myThid )       O                           af, myThid )
542          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
543            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
544       I                              vTrans, vVel, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
545       O                              af, myThid )       O                              af, myThid )
546          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
547            CALL GAD_DST3_ADV_Y(      bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y(      bi,bj,k, .TRUE., deltaTLev(k),
548       I                              vTrans, vVel, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
549       O                              af, myThid )       O                              af, myThid )
550          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
551            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .TRUE., deltaTLev(k),
552       I                              vTrans, vVel, maskLocS, localTij,       I                              vTrans, vFld, maskLocS, localTij,
553         O                              af, myThid )
554    #ifndef ALLOW_AUTODIFF_TAMC
555            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
556              CALL GAD_OS7MP_ADV_Y(     bi,bj,k, .TRUE., deltaTLev(k),
557         I                              vTrans, vFld, maskLocS, localTij,
558       O                              af, myThid )       O                              af, myThid )
559    #endif
560          ELSE          ELSE
561           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
562          ENDIF          ENDIF
563    
 C-     Advective flux in Y : done  
        ENDIF  
   
 #ifndef ALLOW_AUTODIFF_TAMC  
564  C-     Internal exchange for next calculations in X  C-     Internal exchange for next calculations in X
565         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
566           CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
567         &                              localTij, bi,bj, myThid )
568            ENDIF
569    
570    C-     Advective flux in Y : done
571         ENDIF         ENDIF
 #endif  
572    
573  C-     Update the local tracer field where needed:  C-     Update the local tracer field where needed:
574    C      use "maksInC" to prevent updating tracer field in OB regions
575    
576  C      update in overlap-Only  C      update in overlap-Only
577         IF ( overlapOnly ) THEN         IF ( overlapOnly ) THEN
578          jMinUpd = 1-Oly+1          jMinUpd = 1-Oly+1
579          jMaxUpd = sNy+Oly-1          jMaxUpd = sNy+Oly-1
580  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
581  C         in corner region) but safer to keep them.  C         in corner region) but safer to keep them.
582          IF ( S_edge ) jMinUpd = 1          IF ( S_edge ) jMinUpd = 1
583          IF ( N_edge ) jMaxUpd = sNy          IF ( N_edge ) jMaxUpd = sNy
# Line 539  C         in corner region) but safer to Line 585  C         in corner region) but safer to
585          IF ( W_edge ) THEN          IF ( W_edge ) THEN
586           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
587            DO i=1-Olx,0            DO i=1-Olx,0
588             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
589       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
590       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
591         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
592       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
593       &         -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))
594       &        )       &        )*maskInC(i,j,bi,bj)
595            ENDDO            ENDDO
596           ENDDO           ENDDO
597          ENDIF          ENDIF
598          IF ( E_edge ) THEN          IF ( E_edge ) THEN
599           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
600            DO i=sNx+1,sNx+Olx            DO i=sNx+1,sNx+Olx
601             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
602       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
603       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
604         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
605       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
606       &         -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))
607       &        )       &        )*maskInC(i,j,bi,bj)
608            ENDDO            ENDDO
609           ENDDO           ENDDO
610          ENDIF          ENDIF
# Line 569  C      do not only update the overlap Line 617  C      do not only update the overlap
617          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
618          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
619           DO i=iMinUpd,iMaxUpd           DO i=iMinUpd,iMaxUpd
620             localTij(i,j)=localTij(i,j)-dTtracerLev(k)*             localTij(i,j) = localTij(i,j)
621       &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
622       &       *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
623         &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
624       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
625       &         -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))
626       &        )       &        )*maskInC(i,j,bi,bj)
627           ENDDO           ENDDO
628          ENDDO          ENDDO
629  C-      keep advective flux (for diagnostics)  C-      keep advective flux (for diagnostics)
# Line 585  C-      keep advective flux (for diagnos Line 634  C-      keep advective flux (for diagnos
634          ENDDO          ENDDO
635    
636  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
637    #ifdef OBCS_MULTIDIM_OLD_VERSION
638  C-     Apply open boundary conditions  C-     Apply open boundary conditions
639          IF (useOBCS) THEN          IF (useOBCS) THEN
640           IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN           IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
# Line 593  C-     Apply open boundary conditions Line 643  C-     Apply open boundary conditions
643            CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )            CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
644  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
645           ELSEIF (tracerIdentity.GE.GAD_TR1) THEN           ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
646            CALL OBCS_APPLY_PTRACER( bi, bj, k,            CALL OBCS_APPLY_PTRACER( bi, bj, k,
647       &         tracerIdentity-GAD_TR1+1, localTij, myThid )       &         tracerIdentity-GAD_TR1+1, localTij, myThid )
648  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
649           ENDIF           ENDIF
650          ENDIF          ENDIF
651    #endif /* OBCS_MULTIDIM_OLD_VERSION */
652  #endif /* ALLOW_OBCS */  #endif /* ALLOW_OBCS */
653    
654  C      end if/else update overlap-Only  C      end if/else update overlap-Only
# Line 614  C-    explicit advection is done ; store Line 665  C-    explicit advection is done ; store
665          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
666           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
667            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
668       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
669           ENDDO           ENDDO
670          ENDDO          ENDDO
671        ELSE        ELSE
672  C-    horizontal advection done; store intermediate result in 3D array:  C-    horizontal advection done; store intermediate result in 3D array:
673         DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
674          DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
675           localTijk(i,j,k)=localTij(i,j)            localTijk(i,j,k)=localTij(i,j)
676             ENDDO
677          ENDDO          ENDDO
        ENDDO  
678        ENDIF        ENDIF
679    
680  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
681          IF ( useDiagnostics ) THEN          IF ( doDiagAdvX ) THEN
682            diagName = 'ADVx'//diagSufx            diagName = 'ADVx'//diagSufx
683            CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)
684            ENDIF
685            IF ( doDiagAdvY ) THEN
686            diagName = 'ADVy'//diagSufx            diagName = 'ADVy'//diagSufx
687            CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)
688          ENDIF          ENDIF
# Line 654  C---+----1----+----2----+----3----+----4 Line 707  C---+----1----+----2----+----3----+----4
707        IF ( .NOT.implicitAdvection ) THEN        IF ( .NOT.implicitAdvection ) THEN
708  C--   Start of k loop for vertical flux  C--   Start of k loop for vertical flux
709         DO k=Nr,1,-1         DO k=Nr,1,-1
710  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
711           kkey = (igadkey-1)*Nr + k           kkey = (igadkey-1)*Nr + (Nr-k+1)
712  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
713  C--   kup    Cycles through 1,2 to point to w-layer above  C--   kUp    Cycles through 1,2 to point to w-layer above
714  C--   kDown  Cycles through 2,1 to point to w-layer below  C--   kDown  Cycles through 2,1 to point to w-layer below
715          kup  = 1+MOD(k+1,2)          kUp  = 1+MOD(k+1,2)
716          kDown= 1+MOD(k,2)          kDown= 1+MOD(k,2)
717  c       kp1=min(Nr,k+1)  c       kp1=min(Nr,k+1)
718          kp1Msk=1.          kp1Msk=1.
719          if (k.EQ.Nr) kp1Msk=0.          if (k.EQ.Nr) kp1Msk=0.
720    
721    #ifdef ALLOW_AUTODIFF_TAMC
722    CADJ STORE rtrans(:,:)  =
723    CADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
724    cphCADJ STORE wfld(:,:)  =
725    cphCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
726    #endif
727    
728  C-- Compute Vertical transport  C-- Compute Vertical transport
729  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
730  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
# Line 675  C- a hack to prevent Water-Vapor vert.tr Line 735  C- a hack to prevent Water-Vapor vert.tr
735          IF ( k.EQ.1 ) THEN          IF ( k.EQ.1 ) THEN
736  #endif  #endif
737    
738    #ifdef ALLOW_AUTODIFF_TAMC
739    cphmultiCADJ STORE wfld(:,:)  =
740    cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
741    #endif /* ALLOW_AUTODIFF_TAMC */
742    
743  C- Surface interface :  C- Surface interface :
744           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
745            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
746             rTransKp1(i,j) = kp1Msk*rTrans(i,j)             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
747               wFld(i,j)   = 0.
748             rTrans(i,j) = 0.             rTrans(i,j) = 0.
749             fVerT(i,j,kUp) = 0.             fVerT(i,j,kUp) = 0.
750            ENDDO            ENDDO
751           ENDDO           ENDDO
752    
753          ELSE          ELSE
 C- Interior interface :  
754    
755    #ifdef ALLOW_AUTODIFF_TAMC
756    cphmultiCADJ STORE wfld(:,:)  =
757    cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
758    #endif /* ALLOW_AUTODIFF_TAMC */
759    
760    C- Interior interface :
761           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
762            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
763             rTransKp1(i,j) = kp1Msk*rTrans(i,j)             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
764               wFld(i,j)   = wVel(i,j,k,bi,bj)
765             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
766         &                 *deepFac2F(k)*rhoFacF(k)
767       &                 *maskC(i,j,k-1,bi,bj)       &                 *maskC(i,j,k-1,bi,bj)
768             fVerT(i,j,kUp) = 0.             fVerT(i,j,kUp) = 0.
769            ENDDO            ENDDO
# Line 698  C- Interior interface : Line 771  C- Interior interface :
771    
772  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
773  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
774           IF (useGMRedi)           IF (useGMRedi)
775       &   CALL GMREDI_CALC_WFLOW(       &     CALL GMREDI_CALC_WFLOW(
776       &                    rTrans, bi, bj, k, myThid)       U                 wFld, rTrans,
777         I                 k, bi, bj, myThid )
778  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
779    
780  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
781  CADJ STORE localTijk(:,:,k)    cphmultiCADJ STORE localTijk(:,:,k)
782  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
783  CADJ STORE rTrans(:,:)    cphmultiCADJ STORE rTrans(:,:)
784  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
785  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
786    
787  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
788           IF ( advectionScheme.EQ.ENUM_UPWIND_1RST           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
789       &      .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
790             CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,
791       I                            dTtracerLev(k),rTrans,wVel,localTijk,       I                            deltaTLev(k),rTrans,wFld,localTijk,
792       O                            fVerT(1-Olx,1-Oly,kUp), myThid )       O                            fVerT(1-Olx,1-Oly,kUp), myThid )
793           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
794             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k),             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
795       I                               rTrans, wVel, localTijk,       I                               rTrans, wFld, localTijk,
796       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
797           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
798             CALL GAD_DST3_ADV_R(      bi,bj,k, dTtracerLev(k),             CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTLev(k),
799       I                               rTrans, wVel, localTijk,       I                               rTrans, wFld, localTijk,
800       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
801           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
802             CALL GAD_DST3FL_ADV_R(    bi,bj,k, dTtracerLev(k),             CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTLev(k),
803       I                               rTrans, wVel, localTijk,       I                               rTrans, wFld, localTijk,
804       O                               fVerT(1-Olx,1-Oly,kUp), myThid )       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
805    #ifndef ALLOW_AUTODIFF_TAMC
806             ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
807               CALL GAD_OS7MP_ADV_R(     bi,bj,k, deltaTLev(k),
808         I                               rTrans, wFld, localTijk,
809         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
810    #endif
811           ELSE           ELSE
812            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
813           ENDIF           ENDIF
# Line 735  C-    Compute vertical advective flux in Line 815  C-    Compute vertical advective flux in
815  C- end Surface/Interior if bloc  C- end Surface/Interior if bloc
816          ENDIF          ENDIF
817    
818  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
819  CADJ STORE rTrans(:,:)    cphmultiCADJ STORE rTrans(:,:)
820  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
821  CADJ STORE rTranskp1(:,:)    cphmultiCADJ STORE rTranskp1(:,:)
822  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
823    cph --- following storing of fVerT is critical for correct
824    cph --- gradient with multiDimAdvection
825    cph --- Without it, kDown component is not properly recomputed
826    cph --- This is a TAF bug (and no warning available)
827    CADJ STORE fVerT(:,:,:)
828    CADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
829  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
830    
831  C--   Divergence of vertical fluxes  C--   Divergence of vertical fluxes
832          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
833           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
834            localTij(i,j)=localTijk(i,j,k)-dTtracerLev(k)*            localTij(i,j) = localTijk(i,j,k)
835       &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
836       &     *recip_rA(i,j,bi,bj)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
837       &     *( fVerT(i,j,kDown)-fVerT(i,j,kUp)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
838       &       -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))       &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)
839       &      )*rkSign       &         -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))
840         &        )*rkSign
841            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
842       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
843           ENDDO           ENDDO
844          ENDDO          ENDDO
845    
846  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
847          IF ( useDiagnostics ) THEN          IF ( doDiagAdvR ) THEN
848            diagName = 'ADVr'//diagSufx            diagName = 'ADVr'//diagSufx
849            CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp),            CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp),
850       &                           diagName, k,1, 2,bi,bj, myThid)       &                           diagName, k,1, 2,bi,bj, myThid)
# Line 767  C--   Divergence of vertical fluxes Line 854  C--   Divergence of vertical fluxes
854  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
855         ENDDO         ENDDO
856  C--   end of if not.implicitAdvection block  C--   end of if not.implicitAdvection block
857        ENDIF        ENDIF
858    
859        RETURN        RETURN
860        END        END
Legend:
Removed from v.1.38  
changed lines
  Added in v.1.64
  ViewVC Help
Powered by ViewVC 1.1.22