/[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.5 by adcroft, Fri Sep 21 13:11:43 2001 UTC revision 1.29 by jmc, Fri Sep 24 16:52:44 2004 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
 CBOI  
 C !TITLE: pkg/generic\_advdiff  
 C !AUTHORS: adcroft@mit.edu  
 C !INTRODUCTION:  
 C \section{Generica Advection Diffusion Package}  
 C  
 C Package "generic\_advdiff" provides a common set of routines for calculating  
 C advective/diffusive fluxes for tracers (cell centered quantities on a C-grid).  
 C  
 C Many different advection schemes are available: the standard centered  
 C second order, centered fourth order and upwind biased third order schemes  
 C are known as linear methods and require some stable time-stepping method  
 C such as Adams-Bashforth. Alternatives such as flux-limited schemes are  
 C stable in the forward sense and are best combined with the multi-dimensional  
 C method provided in gad\_advection.  
 C  
 C There are two high-level routines:  
 C  \begin{itemize}  
 C  \item{GAD\_CALC\_RHS} calculates all fluxes at time level "n" and is used  
 C  for the standard linear schemes. This must be used in conjuction with  
 C  Adams-Bashforth time-stepping. Diffusive and parameterized fluxes are  
 C  always calculated here.  
 C  \item{GAD\_ADVECTION} calculates just the advective fluxes using the  
 C  non-linear schemes and can not be used in conjuction with Adams-Bashforth  
 C  time-stepping.  
 C  \end{itemize}  
 CEOI  
   
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7  CBOP  CBOP
8  C !ROUTINE: GAD_ADVECTION  C !ROUTINE: GAD_ADVECTION
9    
10  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
11        SUBROUTINE GAD_ADVECTION(bi,bj,advectionScheme,tracerIdentity,        SUBROUTINE GAD_ADVECTION(
12       U                         Tracer,Gtracer,       I     implicitAdvection, advectionScheme, vertAdvecScheme,
13       I                         myTime,myIter,myThid)       I     tracerIdentity,
14         I     uVel, vVel, wVel, tracer,
15         O     gTracer,
16         I     bi,bj, myTime,myIter,myThid)
17    
18  C !DESCRIPTION:  C !DESCRIPTION:
19  C Calculates the tendancy of a tracer due to advection.  C Calculates the tendancy of a tracer due to advection.
# Line 63  C !USES: =============================== Line 39  C !USES: ===============================
39  #include "SIZE.h"  #include "SIZE.h"
40  #include "EEPARAMS.h"  #include "EEPARAMS.h"
41  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
42  #include "GRID.h"  #include "GRID.h"
43  #include "GAD.h"  #include "GAD.h"
44    #ifdef ALLOW_AUTODIFF_TAMC
45    # include "tamc.h"
46    # include "tamc_keys.h"
47    # ifdef ALLOW_PTRACERS
48    #  include "PTRACERS_SIZE.h"
49    # endif
50    #endif
51    #ifdef ALLOW_EXCH2
52    #include "W2_EXCH2_TOPOLOGY.h"
53    #include "W2_EXCH2_PARAMS.h"
54    #endif /* ALLOW_EXCH2 */
55    
56  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
57  C  bi,bj                :: tile indices  C  implicitAdvection :: implicit vertical advection (later on)
58  C  advectionScheme      :: advection scheme to use  C  advectionScheme   :: advection scheme to use (Horizontal plane)
59  C  tracerIdentity       :: identifier for the tracer (required only for OBCS)  C  vertAdvecScheme   :: advection scheme to use (vertical direction)
60  C  Tracer               :: tracer field  C  tracerIdentity    :: tracer identifier (required only for OBCS)
61  C  myTime               :: current time  C  uVel              :: velocity, zonal component
62  C  myIter               :: iteration number  C  vVel              :: velocity, meridional component
63  C  myThid               :: thread number  C  wVel              :: velocity, vertical component
64        INTEGER bi,bj  C  tracer            :: tracer field
65        INTEGER advectionScheme  C  bi,bj             :: tile indices
66    C  myTime            :: current time
67    C  myIter            :: iteration number
68    C  myThid            :: thread number
69          LOGICAL implicitAdvection
70          INTEGER advectionScheme, vertAdvecScheme
71        INTEGER tracerIdentity        INTEGER tracerIdentity
72        _RL Gtracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
73          _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
74          _RL wVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
75          _RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
76          INTEGER bi,bj
77        _RL myTime        _RL myTime
78        INTEGER myIter        INTEGER myIter
79        INTEGER myThid        INTEGER myThid
80    
81  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
82  C  gTracer              :: tendancy array  C  gTracer           :: tendancy array
83        _RL Tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
84    
85  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
86  C  maskUp               :: 2-D array for mask at W points  C  maskUp        :: 2-D array for mask at W points
87  C  iMin,iMax,jMin,jMax  :: loop range for called routines  C  maskLocW      :: 2-D array for mask at West points
88  C  i,j,k                :: loop indices  C  maskLocS      :: 2-D array for mask at South points
89  C  kup                  :: index into 2 1/2D array, toggles between 1 and 2  C  iMin,iMax,    :: loop range for called routines
90  C  kdown                :: index into 2 1/2D array, toggles between 2 and 1  C  jMin,jMax     :: loop range for called routines
91  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr  C  i,j,k         :: loop indices
92  C  xA,yA                :: areas of X and Y face of tracer cells  C  kup           :: index into 2 1/2D array, toggles between 1 and 2
93  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points  C  kdown         :: index into 2 1/2D array, toggles between 2 and 1
94  C  af                   :: 2-D array for horizontal advective flux  C  kp1           :: =k+1 for k<Nr, =Nr for k=Nr
95  C  fVerT                :: 2 1/2D arrays for vertical advective flux  C  xA,yA         :: areas of X and Y face of tracer cells
96  C  localTij             :: 2-D array used as temporary local copy of tracer fld  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
97  C  localTijk            :: 3-D array used as temporary local copy of tracer fld  C  rTrans        :: 2-D arrays of volume transports at W points
98  C  kp1Msk               :: flag (0,1) to act as over-riding mask for W levels  C  rTransKp1     :: vertical volume transport at interface k+1
99  C  calc_fluxes_X        :: logical to indicate to calculate fluxes in X dir  C  afx           :: 2-D array for horizontal advective flux, x direction
100  C  calc_fluxes_Y        :: logical to indicate to calculate fluxes in Y dir  C  afy           :: 2-D array for horizontal advective flux, y direction
101  C  nipass               :: number of passes to make in multi-dimensional method  C  fVerT         :: 2 1/2D arrays for vertical advective flux
102  C  ipass                :: number of the current pass being made  C  localTij      :: 2-D array, temporary local copy of tracer fld
103    C  localTijk     :: 3-D array, temporary local copy of tracer fld
104    C  kp1Msk        :: flag (0,1) for over-riding mask for W levels
105    C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
106    C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
107    C  nipass        :: number of passes in multi-dimensional method
108    C  ipass         :: number of the current pass being made
109    C  myTile        :: variables used to determine which cube face
110    C  nCFace        :: owns a tile for cube grid runs using
111    C                :: multi-dim advection.
112        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113          _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114          _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
116        INTEGER i,j,k,kup,kDown,kp1        INTEGER i,j,k,kup,kDown
117        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123          _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124          _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
126        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127        _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
128        _RL kp1Msk        _RL kp1Msk
129        LOGICAL calc_fluxes_X,calc_fluxes_Y        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
130        INTEGER nipass,ipass        INTEGER nipass,ipass
131          INTEGER myTile, nCFace
132  CEOP  CEOP
133    
134    #ifdef ALLOW_AUTODIFF_TAMC
135              act0 = tracerIdentity - 1
136              max0 = maxpass
137              act1 = bi - myBxLo(myThid)
138              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
139              act2 = bj - myByLo(myThid)
140              max2 = myByHi(myThid) - myByLo(myThid) + 1
141              act3 = myThid - 1
142              max3 = nTx*nTy
143              act4 = ikey_dynamics - 1
144              igadkey = (act0 + 1)
145         &                      + act1*max0
146         &                      + act2*max0*max1
147         &                      + act3*max0*max1*max2
148         &                      + act4*max0*max1*max2*max3
149              if (tracerIdentity.GT.maxpass) then
150                 print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
151                 STOP 'maxpass seems smaller than tracerIdentity'
152              endif
153    #endif /* ALLOW_AUTODIFF_TAMC */
154    
155  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
156  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
157  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
# Line 136  C     uninitialised but inert locations. Line 166  C     uninitialised but inert locations.
166          rTrans(i,j)  = 0. _d 0          rTrans(i,j)  = 0. _d 0
167          fVerT(i,j,1) = 0. _d 0          fVerT(i,j,1) = 0. _d 0
168          fVerT(i,j,2) = 0. _d 0          fVerT(i,j,2) = 0. _d 0
169            rTransKp1(i,j)= 0. _d 0
170         ENDDO         ENDDO
171        ENDDO        ENDDO
172    
# Line 146  C     uninitialised but inert locations. Line 177  C     uninitialised but inert locations.
177    
178  C--   Start of k loop for horizontal fluxes  C--   Start of k loop for horizontal fluxes
179        DO k=1,Nr        DO k=1,Nr
180    #ifdef ALLOW_AUTODIFF_TAMC
181             kkey = (igadkey-1)*Nr + k
182    CADJ STORE tracer(:,:,k,bi,bj) =
183    CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte
184    #endif /* ALLOW_AUTODIFF_TAMC */
185    
186  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
187        CALL CALC_COMMON_FACTORS (        CALL CALC_COMMON_FACTORS (
# Line 153  C--   Get temporary terms used by tenden Line 189  C--   Get temporary terms used by tenden
189       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
190       I         myThid)       I         myThid)
191    
192  C--   Make local copy of tracer array  #ifdef ALLOW_GMREDI
193    C--   Residual transp = Bolus transp + Eulerian transp
194           IF (useGMRedi)
195         &   CALL GMREDI_CALC_UVFLOW(
196         &            uTrans, vTrans, bi, bj, k, myThid)
197    #endif /* ALLOW_GMREDI */
198    
199    C--   Make local copy of tracer array and mask West & South
200        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
201         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
202          localTij(i,j)=tracer(i,j,k,bi,bj)          localTij(i,j)=tracer(i,j,k,bi,bj)
203            maskLocW(i,j)=maskW(i,j,k,bi,bj)
204            maskLocS(i,j)=maskS(i,j,k,bi,bj)
205         ENDDO         ENDDO
206        ENDDO        ENDDO
207    
208        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
209            withSigns = .FALSE.
210            CALL FILL_CS_CORNER_UV_RS(
211         &            withSigns, maskLocW,maskLocS, bi,bj, myThid )
212          ENDIF
213    #ifdef ALLOW_EXCH2
214           myTile = W2_myTileList(bi)
215           nCFace = exch2_myFace(myTile)
216    #else
217           nCFace = bi
218    #endif
219          IF (useCubedSphereExchange) THEN
220    
221         nipass=3         nipass=3
222    #ifdef ALLOW_AUTODIFF_TAMC
223           if ( nipass.GT.maxcube )
224         &      STOP 'maxcube needs to be = 3'
225    #endif
226        ELSE        ELSE
227         nipass=1         nipass=1
228        ENDIF        ENDIF
        nipass=1  
229    
230  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
231    C--   For cube need one pass for each of red, green and blue axes.
232        DO ipass=1,nipass        DO ipass=1,nipass
233    #ifdef ALLOW_AUTODIFF_TAMC
234             passkey = ipass + (k-1)      *maxcube
235         &                   + (igadkey-1)*maxcube*Nr
236             IF (nipass .GT. maxpass) THEN
237              STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
238             ENDIF
239    #endif /* ALLOW_AUTODIFF_TAMC */
240    
241        IF (nipass.EQ.3) THEN        IF (nipass.EQ.3) THEN
242         calc_fluxes_X=.FALSE.         calc_fluxes_X=.FALSE.
243         calc_fluxes_Y=.FALSE.         calc_fluxes_Y=.FALSE.
244         IF (ipass.EQ.1 .AND. (bi.EQ.1 .OR. bi.EQ.2) ) THEN         IF (ipass.EQ.1 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.2) ) THEN
245          calc_fluxes_X=.TRUE.          calc_fluxes_X=.TRUE.
246         ELSEIF (ipass.EQ.1 .AND. (bi.EQ.4 .OR. bi.EQ.5) ) THEN         ELSEIF (ipass.EQ.1 .AND. (nCFace.EQ.4 .OR. nCFace.EQ.5) ) THEN
247          calc_fluxes_Y=.TRUE.          calc_fluxes_Y=.TRUE.
248         ELSEIF (ipass.EQ.2 .AND. (bi.EQ.1 .OR. bi.EQ.6) ) THEN         ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.6) ) THEN
249          calc_fluxes_Y=.TRUE.          calc_fluxes_Y=.TRUE.
250         ELSEIF (ipass.EQ.2 .AND. (bi.EQ.3 .OR. bi.EQ.4) ) THEN         ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.3 .OR. nCFace.EQ.4) ) THEN
251          calc_fluxes_X=.TRUE.          calc_fluxes_X=.TRUE.
252         ELSEIF (ipass.EQ.3 .AND. (bi.EQ.2 .OR. bi.EQ.3) ) THEN         ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.2 .OR. nCFace.EQ.3) ) THEN
253          calc_fluxes_Y=.TRUE.          calc_fluxes_Y=.TRUE.
254         ELSEIF (ipass.EQ.3 .AND. (bi.EQ.5 .OR. bi.EQ.6) ) THEN         ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.5 .OR. nCFace.EQ.6) ) THEN
255          calc_fluxes_X=.TRUE.          calc_fluxes_X=.TRUE.
256         ENDIF         ENDIF
257        ELSE        ELSE
# Line 191  C--   Multiple passes for different dire Line 259  C--   Multiple passes for different dire
259         calc_fluxes_Y=.TRUE.         calc_fluxes_Y=.TRUE.
260        ENDIF        ENDIF
261    
262    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
263  C--   X direction  C--   X direction
264        IF (calc_fluxes_X) THEN        IF (calc_fluxes_X) THEN
265    
266  C--   Internal exchange for calculations in X  C--   Internal exchange for calculations in X
267        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
268         DO j=1,Oly           CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )
         DO i=1,Olx  
          localTij( 1-i , 1-j )=localTij( 1-j ,    i    )  
          localTij( 1-i ,sNy+j)=localTij( 1-j , sNy+1-i )  
          localTij(sNx+i, 1-j )=localTij(sNx+j,    i    )  
          localTij(sNx+i,sNy+j)=localTij(sNx+j, sNy+1-i )  
         ENDDO  
        ENDDO  
269        ENDIF        ENDIF
270    
271  C-    Advective flux in X  C-    Advective flux in X
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          af(i,j) = 0.          afx(i,j) = 0.
275         ENDDO         ENDDO
276        ENDDO        ENDDO
277    
278    #ifdef ALLOW_AUTODIFF_TAMC
279    #ifndef DISABLE_MULTIDIM_ADVECTION
280    CADJ STORE localTij(:,:)  =
281    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
282    #endif
283    #endif /* ALLOW_AUTODIFF_TAMC */
284    
285        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
286         CALL GAD_FLUXLIMIT_ADV_X(          CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, deltaTtracer,
287       &      bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)       I                            uTrans, uVel, maskLocW, localTij,
288         O                            afx, myThid )
289        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
290         CALL GAD_DST3_ADV_X(          CALL GAD_DST3_ADV_X(      bi,bj,k, deltaTtracer,
291       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)       I                            uTrans, uVel, maskLocW, localTij,
292         O                            afx, myThid )
293        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
294         CALL GAD_DST3FL_ADV_X(          CALL GAD_DST3FL_ADV_X(    bi,bj,k, deltaTtracer,
295       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)       I                            uTrans, uVel, maskLocW, localTij,
296         O                            afx, myThid )
297        ELSE        ELSE
298         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'         STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
299        ENDIF        ENDIF
300    
301        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
302         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
303          localTij(i,j)=localTij(i,j)-deltaTtracer*          localTij(i,j)=localTij(i,j)-deltaTtracer*
304       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
305       &    *recip_rA(i,j,bi,bj)       &    *recip_rA(i,j,bi,bj)
306       &    *( af(i+1,j)-af(i,j)       &    *( afx(i+1,j)-afx(i,j)
307       &      -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))
308       &     )       &     )
309         ENDDO         ENDDO
# Line 249  C--   Apply open boundary conditions Line 323  C--   Apply open boundary conditions
323  C--   End of X direction  C--   End of X direction
324        ENDIF        ENDIF
325    
326    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
327  C--   Y direction  C--   Y direction
328        IF (calc_fluxes_Y) THEN        IF (calc_fluxes_Y) THEN
329    
330  C--   Internal exchange for calculations in Y  C--   Internal exchange for calculations in Y
331        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
332         DO j=1,Oly           CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )
         DO i=1,Olx  
          localTij( 1-i , 1-j )=localTij(   j   , 1-i )  
          localTij( 1-i ,sNy+j)=localTij(   j   ,sNy+i)  
          localTij(sNx+i, 1-j )=localTij(sNx+1-j, 1-i )  
          localTij(sNx+i,sNy+j)=localTij(sNx+1-j,sNy+i)  
         ENDDO  
        ENDDO  
333        ENDIF        ENDIF
334    
335  C-    Advective flux in Y  C-    Advective flux in Y
336        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
337         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
338          af(i,j) = 0.          afy(i,j) = 0.
339         ENDDO         ENDDO
340        ENDDO        ENDDO
341    
342    #ifdef ALLOW_AUTODIFF_TAMC
343    #ifndef DISABLE_MULTIDIM_ADVECTION
344    CADJ STORE localTij(:,:)  =
345    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
346    #endif
347    #endif /* ALLOW_AUTODIFF_TAMC */
348    
349        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
350         CALL GAD_FLUXLIMIT_ADV_Y(          CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, deltaTtracer,
351       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)       I                            vTrans, vVel, maskLocS, localTij,
352         O                            afy, myThid )
353        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
354         CALL GAD_DST3_ADV_Y(          CALL GAD_DST3_ADV_Y(      bi,bj,k, deltaTtracer,
355       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)       I                            vTrans, vVel, maskLocS, localTij,
356         O                            afy, myThid )
357        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
358         CALL GAD_DST3FL_ADV_Y(          CALL GAD_DST3FL_ADV_Y(    bi,bj,k, deltaTtracer,
359       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)       I                            vTrans, vVel, maskLocS, localTij,
360         O                            afy, myThid )
361        ELSE        ELSE
362         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
363        ENDIF        ENDIF
364    
365        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
366         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
367          localTij(i,j)=localTij(i,j)-deltaTtracer*          localTij(i,j)=localTij(i,j)-deltaTtracer*
368       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
369       &    *recip_rA(i,j,bi,bj)       &    *recip_rA(i,j,bi,bj)
370       &    *( af(i,j+1)-af(i,j)       &    *( afy(i,j+1)-afy(i,j)
371       &      -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))
372       &     )       &     )
373         ENDDO         ENDDO
# Line 307  C--   Apply open boundary conditions Line 387  C--   Apply open boundary conditions
387  C--   End of Y direction  C--   End of Y direction
388        ENDIF        ENDIF
389    
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         localTijk(i,j,k)=localTij(i,j)  
        ENDDO  
       ENDDO  
   
390  C--   End of ipass loop  C--   End of ipass loop
391        ENDDO        ENDDO
392    
393          IF ( implicitAdvection ) THEN
394    C-    explicit advection is done ; store tendency in gTracer:
395            DO j=1-Oly,sNy+Oly
396             DO i=1-Olx,sNx+Olx
397              gTracer(i,j,k,bi,bj)=
398         &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
399             ENDDO
400            ENDDO
401          ELSE
402    C-    horizontal advection done; store intermediate result in 3D array:
403           DO j=1-Oly,sNy+Oly
404            DO i=1-Olx,sNx+Olx
405             localTijk(i,j,k)=localTij(i,j)
406            ENDDO
407           ENDDO
408          ENDIF
409    
410    #ifdef ALLOW_DEBUG
411          IF ( debugLevel .GE. debLevB
412         &   .AND. k.EQ.3 .AND. myIter.EQ.1+nIter0
413         &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
414         &   .AND. useCubedSphereExchange ) THEN
415            CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION',
416         &             afx,afy, k, standardMessageUnit,bi,bj,myThid )
417          ENDIF
418    #endif /* ALLOW_DEBUG */
419    
420  C--   End of K loop for horizontal fluxes  C--   End of K loop for horizontal fluxes
421        ENDDO        ENDDO
422    
423  C--   Start of k loop for vertical flux  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
       DO k=Nr,1,-1  
424    
425          IF ( .NOT.implicitAdvection ) THEN
426    C--   Start of k loop for vertical flux
427           DO k=Nr,1,-1
428    #ifdef ALLOW_AUTODIFF_TAMC
429             kkey = (igadkey-1)*Nr + k
430    #endif /* ALLOW_AUTODIFF_TAMC */
431  C--   kup    Cycles through 1,2 to point to w-layer above  C--   kup    Cycles through 1,2 to point to w-layer above
432  C--   kDown  Cycles through 2,1 to point to w-layer below  C--   kDown  Cycles through 2,1 to point to w-layer below
433        kup  = 1+MOD(k+1,2)          kup  = 1+MOD(k+1,2)
434        kDown= 1+MOD(k,2)          kDown= 1+MOD(k,2)
435    c       kp1=min(Nr,k+1)
436  C--   Get temporary terms used by tendency routines          kp1Msk=1.
437        CALL CALC_COMMON_FACTORS (          if (k.EQ.Nr) kp1Msk=0.
438       I         bi,bj,iMin,iMax,jMin,jMax,k,  
439       O         xA,yA,uTrans,vTrans,rTrans,maskUp,  C-- Compute Vertical transport
440       I         myThid)  #ifdef ALLOW_AIM
441    C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
442            IF ( k.EQ.1 .OR.
443         &     (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
444         &              ) THEN
445    #else
446            IF ( k.EQ.1 ) THEN
447    #endif
448    
449    C- Surface interface :
450             DO j=1-Oly,sNy+Oly
451              DO i=1-Olx,sNx+Olx
452               rTransKp1(i,j) = kp1Msk*rTrans(i,j)
453               rTrans(i,j) = 0.
454               fVerT(i,j,kUp) = 0.
455              ENDDO
456             ENDDO
457    
458            ELSE
459    C- Interior interface :
460    
461             DO j=1-Oly,sNy+Oly
462              DO i=1-Olx,sNx+Olx
463               rTransKp1(i,j) = kp1Msk*rTrans(i,j)
464               rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
465         &                 *maskC(i,j,k-1,bi,bj)
466               fVerT(i,j,kUp) = 0.
467              ENDDO
468             ENDDO
469    
470    #ifdef ALLOW_GMREDI
471    C--   Residual transp = Bolus transp + Eulerian transp
472             IF (useGMRedi)
473         &   CALL GMREDI_CALC_WFLOW(
474         &                    rTrans, bi, bj, k, myThid)
475    #endif /* ALLOW_GMREDI */
476    
477    #ifdef ALLOW_AUTODIFF_TAMC
478    CADJ STORE localTijk(:,:,k)  
479    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
480    CADJ STORE rTrans(:,:)  
481    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
482    #endif /* ALLOW_AUTODIFF_TAMC */
483    
 C-    Advective flux in R  
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         af(i,j) = 0.  
        ENDDO  
       ENDDO  
   
 C     Note: wVel needs to be masked  
       IF (K.GE.2) THEN  
484  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
485         IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
486          CALL GAD_FLUXLIMIT_ADV_R(  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
487       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTtracer,
488         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I                               rTrans, wVel, localTijk,
489          CALL GAD_DST3_ADV_R(       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
490       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
491         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN             CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTtracer,
492          CALL GAD_DST3FL_ADV_R(       I                               rTrans, wVel, localTijk,
493       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
494         ELSE           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
495          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'             CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTtracer,
496         ENDIF       I                               rTrans, wVel, localTijk,
497  C-    Surface "correction" term at k>1 :       O                               fVerT(1-Olx,1-Oly,kUp), myThid )
498         DO j=1-Oly,sNy+Oly           ELSE
499          DO i=1-Olx,sNx+Olx            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
500           af(i,j) = af(i,j)           ENDIF
501       &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  
502       &             rTrans(i,j)*localTijk(i,j,k)  C- end Surface/Interior if bloc
503            ENDIF
504    
505    #ifdef ALLOW_AUTODIFF_TAMC
506    CADJ STORE rTrans(:,:)  
507    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
508    CADJ STORE rTranskp1(:,:)  
509    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
510    #endif /* ALLOW_AUTODIFF_TAMC */
511    
512    C--   Divergence of vertical fluxes
513            DO j=1-Oly,sNy+Oly
514             DO i=1-Olx,sNx+Olx
515              localTij(i,j)=localTijk(i,j,k)-deltaTtracer*
516         &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
517         &     *recip_rA(i,j,bi,bj)
518         &     *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
519         &       -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
520         &      )*rkFac
521              gTracer(i,j,k,bi,bj)=
522         &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
523             ENDDO
524          ENDDO          ENDDO
        ENDDO  
       ELSE  
 C-    Surface "correction" term at k=1 :  
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          af(i,j) = rTrans(i,j)*localTijk(i,j,k)  
         ENDDO  
        ENDDO  
       ENDIF  
 C-    add the advective flux to fVerT  
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fVerT(i,j,kUp) = af(i,j)  
        ENDDO  
       ENDDO  
   
 C--   Divergence of fluxes  
       kp1=min(Nr,k+1)  
       kp1Msk=1.  
       if (k.EQ.Nr) kp1Msk=0.  
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         localTij(i,j)=localTijk(i,j,k)-deltaTtracer*  
      &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  
      &    *recip_rA(i,j,bi,bj)  
      &    *( fVerT(i,j,kUp)-fVerT(i,j,kDown)  
      &      -tracer(i,j,k,bi,bj)*rA(i,j,bi,bj)*  
      &        (wVel(i,j,k,bi,bj)-kp1Msk*wVel(i,j,kp1,bi,bj))  
      &     )*rkFac  
         gTracer(i,j,k,bi,bj)=  
      &   (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer  
        ENDDO  
       ENDDO  
525    
526  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
527        ENDDO         ENDDO
528    C--   end of if not.implicitAdvection block
529          ENDIF
530    
531        RETURN        RETURN
532        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.22