/[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.25 by heimbach, Wed Jun 30 23:45:35 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    #endif
48    #ifdef ALLOW_EXCH2
49    #include "W2_EXCH2_TOPOLOGY.h"
50    #include "W2_EXCH2_PARAMS.h"
51    #endif /* ALLOW_EXCH2 */
52    
53  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
54  C  bi,bj                :: tile indices  C  implicitAdvection :: implicit vertical advection (later on)
55  C  advectionScheme      :: advection scheme to use  C  advectionScheme   :: advection scheme to use (Horizontal plane)
56  C  tracerIdentity       :: identifier for the tracer (required only for OBCS)  C  vertAdvecScheme   :: advection scheme to use (vertical direction)
57  C  Tracer               :: tracer field  C  tracerIdentity    :: tracer identifier (required only for OBCS)
58  C  myTime               :: current time  C  uVel              :: velocity, zonal component
59  C  myIter               :: iteration number  C  vVel              :: velocity, meridional component
60  C  myThid               :: thread number  C  wVel              :: velocity, vertical component
61        INTEGER bi,bj  C  tracer            :: tracer field
62        INTEGER advectionScheme  C  bi,bj             :: tile indices
63    C  myTime            :: current time
64    C  myIter            :: iteration number
65    C  myThid            :: thread number
66          LOGICAL implicitAdvection
67          INTEGER advectionScheme, vertAdvecScheme
68        INTEGER tracerIdentity        INTEGER tracerIdentity
69        _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)
70          _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
71          _RL wVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
72          _RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
73          INTEGER bi,bj
74        _RL myTime        _RL myTime
75        INTEGER myIter        INTEGER myIter
76        INTEGER myThid        INTEGER myThid
77    
78  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
79  C  gTracer              :: tendancy array  C  gTracer           :: tendancy array
80        _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)
81    
82  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
83  C  maskUp               :: 2-D array for mask at W points  C  maskUp        :: 2-D array for mask at W points
84  C  iMin,iMax,jMin,jMax  :: loop range for called routines  C  iMin,iMax,    :: loop range for called routines
85  C  i,j,k                :: loop indices  C  jMin,jMax     :: loop range for called routines
86  C  kup                  :: index into 2 1/2D array, toggles between 1 and 2  C  i,j,k         :: loop indices
87  C  kdown                :: index into 2 1/2D array, toggles between 2 and 1  C  kup           :: index into 2 1/2D array, toggles between 1 and 2
88  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr  C  kdown         :: index into 2 1/2D array, toggles between 2 and 1
89  C  xA,yA                :: areas of X and Y face of tracer cells  C  kp1           :: =k+1 for k<Nr, =Nr for k=Nr
90  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points  C  xA,yA         :: areas of X and Y face of tracer cells
91  C  af                   :: 2-D array for horizontal advective flux  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
92  C  fVerT                :: 2 1/2D arrays for vertical advective flux  C  rTrans        :: 2-D arrays of volume transports at W points
93  C  localTij             :: 2-D array used as temporary local copy of tracer fld  C  rTransKp1     :: vertical volume transport at interface k+1
94  C  localTijk            :: 3-D array used as temporary local copy of tracer fld  C  af            :: 2-D array for horizontal advective flux
95  C  kp1Msk               :: flag (0,1) to act as over-riding mask for W levels  C  fVerT         :: 2 1/2D arrays for vertical advective flux
96  C  calc_fluxes_X        :: logical to indicate to calculate fluxes in X dir  C  localTij      :: 2-D array, temporary local copy of tracer fld
97  C  calc_fluxes_Y        :: logical to indicate to calculate fluxes in Y dir  C  localTijk     :: 3-D array, temporary local copy of tracer fld
98  C  nipass               :: number of passes to make in multi-dimensional method  C  kp1Msk        :: flag (0,1) for over-riding mask for W levels
99  C  ipass                :: number of the current pass being made  C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
100    C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
101    C  nipass        :: number of passes in multi-dimensional method
102    C  ipass         :: number of the current pass being made
103    C  myTile        :: variables used to determine which cube face
104    C  nCFace        :: owns a tile for cube grid runs using
105    C                :: multi-dim advection.
106        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
108        INTEGER i,j,k,kup,kDown,kp1        INTEGER i,j,k,kup,kDown
109        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
117        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 120  C  ipass                :: number of the Line 119  C  ipass                :: number of the
119        _RL kp1Msk        _RL kp1Msk
120        LOGICAL calc_fluxes_X,calc_fluxes_Y        LOGICAL calc_fluxes_X,calc_fluxes_Y
121        INTEGER nipass,ipass        INTEGER nipass,ipass
122          INTEGER myTile, nCFace
123          LOGICAL southWestCorner
124          LOGICAL southEastCorner
125          LOGICAL northWestCorner
126          LOGICAL northEastCorner
127  CEOP  CEOP
128    
129    #ifdef ALLOW_AUTODIFF_TAMC
130              act0 = tracerIdentity - 1
131              max0 = maxpass
132              act1 = bi - myBxLo(myThid)
133              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
134              act2 = bj - myByLo(myThid)
135              max2 = myByHi(myThid) - myByLo(myThid) + 1
136              act3 = myThid - 1
137              max3 = nTx*nTy
138              act4 = ikey_dynamics - 1
139              igadkey = (act0 + 1)
140         &                      + act1*max0
141         &                      + act2*max0*max1
142         &                      + act3*max0*max1*max2
143         &                      + act4*max0*max1*max2*max3
144              if (tracerIdentity.GT.maxpass) then
145                 print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
146                 STOP 'maxpass seems smaller than tracerIdentity'
147              endif
148    #endif /* ALLOW_AUTODIFF_TAMC */
149    
150  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
151  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
152  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 161  C     uninitialised but inert locations.
161          rTrans(i,j)  = 0. _d 0          rTrans(i,j)  = 0. _d 0
162          fVerT(i,j,1) = 0. _d 0          fVerT(i,j,1) = 0. _d 0
163          fVerT(i,j,2) = 0. _d 0          fVerT(i,j,2) = 0. _d 0
164            rTransKp1(i,j)= 0. _d 0
165         ENDDO         ENDDO
166        ENDDO        ENDDO
167    
# Line 146  C     uninitialised but inert locations. Line 172  C     uninitialised but inert locations.
172    
173  C--   Start of k loop for horizontal fluxes  C--   Start of k loop for horizontal fluxes
174        DO k=1,Nr        DO k=1,Nr
175    #ifdef ALLOW_AUTODIFF_TAMC
176             kkey = (igadkey-1)*Nr + k
177    CADJ STORE tracer(:,:,k,bi,bj) =
178    CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte
179    #endif /* ALLOW_AUTODIFF_TAMC */
180    
181  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
182        CALL CALC_COMMON_FACTORS (        CALL CALC_COMMON_FACTORS (
# Line 153  C--   Get temporary terms used by tenden Line 184  C--   Get temporary terms used by tenden
184       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
185       I         myThid)       I         myThid)
186    
187    #ifdef ALLOW_GMREDI
188    C--   Residual transp = Bolus transp + Eulerian transp
189           IF (useGMRedi)
190         &   CALL GMREDI_CALC_UVFLOW(
191         &            uTrans, vTrans, bi, bj, k, myThid)
192    #endif /* ALLOW_GMREDI */
193    
194  C--   Make local copy of tracer array  C--   Make local copy of tracer array
195        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
196         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
# Line 160  C--   Make local copy of tracer array Line 198  C--   Make local copy of tracer array
198         ENDDO         ENDDO
199        ENDDO        ENDDO
200    
201    cph  The following block is needed for useCubedSphereExchange only,
202    cph  but needs to be set for all cases to avoid spurious
203    cph  TAF dependencies
204           southWestCorner = .TRUE.
205           southEastCorner = .TRUE.
206           northWestCorner = .TRUE.
207           northEastCorner = .TRUE.
208    #ifdef ALLOW_EXCH2
209           myTile = W2_myTileList(bi)
210           nCFace = exch2_myFace(myTile)
211           southWestCorner = exch2_isWedge(myTile).EQ.1
212         &             .AND. exch2_isSedge(myTile).EQ.1
213           southEastCorner = exch2_isEedge(myTile).EQ.1
214         &             .AND. exch2_isSedge(myTile).EQ.1
215           northEastCorner = exch2_isEedge(myTile).EQ.1
216         &             .AND. exch2_isNedge(myTile).EQ.1
217           northWestCorner = exch2_isWedge(myTile).EQ.1
218         &             .AND. exch2_isNedge(myTile).EQ.1
219    #else
220           nCFace = bi
221    #endif
222        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
223    
224         nipass=3         nipass=3
225    #ifdef ALLOW_AUTODIFF_TAMC
226           if ( nipass.GT.maxcube )
227         &      STOP 'maxcube needs to be = 3'
228    #endif
229        ELSE        ELSE
230         nipass=1         nipass=1
231        ENDIF        ENDIF
232         nipass=1  cph       nipass=1
233    
234  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
235    C--   For cube need one pass for each of red, green and blue axes.
236        DO ipass=1,nipass        DO ipass=1,nipass
237    #ifdef ALLOW_AUTODIFF_TAMC
238             passkey = ipass + (k-1)      *maxcube
239         &                   + (igadkey-1)*maxcube*Nr
240             IF (nipass .GT. maxpass) THEN
241              STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
242             ENDIF
243    #endif /* ALLOW_AUTODIFF_TAMC */
244    
245        IF (nipass.EQ.3) THEN        IF (nipass.EQ.3) THEN
246         calc_fluxes_X=.FALSE.         calc_fluxes_X=.FALSE.
247         calc_fluxes_Y=.FALSE.         calc_fluxes_Y=.FALSE.
248         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
249          calc_fluxes_X=.TRUE.          calc_fluxes_X=.TRUE.
250         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
251          calc_fluxes_Y=.TRUE.          calc_fluxes_Y=.TRUE.
252         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
253          calc_fluxes_Y=.TRUE.          calc_fluxes_Y=.TRUE.
254         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
255          calc_fluxes_X=.TRUE.          calc_fluxes_X=.TRUE.
256         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
257          calc_fluxes_Y=.TRUE.          calc_fluxes_Y=.TRUE.
258         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
259          calc_fluxes_X=.TRUE.          calc_fluxes_X=.TRUE.
260         ENDIF         ENDIF
261        ELSE        ELSE
# Line 196  C--   X direction Line 268  C--   X direction
268    
269  C--   Internal exchange for calculations in X  C--   Internal exchange for calculations in X
270        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
271         DO j=1,Oly  C--    For cube face corners we need to duplicate the
272          DO i=1,Olx  C--    i-1 and i+1 values into the null space as follows:
273           localTij( 1-i , 1-j )=localTij( 1-j ,    i    )  C
274           localTij( 1-i ,sNy+j)=localTij( 1-j , sNy+1-i )  C
275           localTij(sNx+i, 1-j )=localTij(sNx+j,    i    )  C      o NW corner: copy T(    0,sNy  ) into T(    0,sNy+1) e.g.
276           localTij(sNx+i,sNy+j)=localTij(sNx+j, sNy+1-i )  C                      |
277          ENDDO  C         x T(0,sNy+1) |
278         ENDDO  C        /\            |
279    C      --||------------|-----------
280    C        ||            |
281    C         x T(0,sNy)   |   x T(1,sNy)
282    C                      |
283    C
284    C      o SW corner: copy T(0,1) into T(0,0) e.g.
285    C                      |
286    C         x T(0,1)     |  x T(1,1)
287    C        ||            |
288    C      --||------------|-----------
289    C        \/            |
290    C         x T(0,0)     |
291    C                      |
292    C
293    C      o NE corner: copy T(sNx+1,sNy  ) into T(sNx+1,sNy+1) e.g.
294    C                      |
295    C                      |   x T(sNx+1,sNy+1)
296    C                      |  /\
297    C      ----------------|--||-------
298    C                      |  ||
299    C         x T(sNx,sNy) |   x T(sNx+1,sNy  )
300    C                      |
301    C      o SE corner: copy T(sNx+1,1    ) into T(sNx+1,0    ) e.g.
302    C                      |
303    C         x T(sNx,1)   |   x T(sNx+1,    1)
304    C                      |  ||
305    C      ----------------|--||-------
306    C                      |  \/
307    C                      |   x T(sNx+1,    0)
308           IF ( southWestCorner ) THEN
309            localTij(0    ,0    )= localTij(0    ,1  )
310           ENDIF
311           IF ( southEastCorner ) THEN
312            localTij(sNx+1,0    )= localTij(sNx+1,1  )
313           ENDIF
314           IF ( northWestCorner ) THEN
315            localTij(0    ,sNy+1)= localTij(0    ,sNy)
316           ENDIF
317           IF ( northEastCorner ) THEN
318            localTij(sNx+1,sNy+1)= localTij(sNx+1,sNy)
319           ENDIF
320        ENDIF        ENDIF
321    
322  C-    Advective flux in X  C-    Advective flux in X
# Line 212  C-    Advective flux in X Line 325  C-    Advective flux in X
325          af(i,j) = 0.          af(i,j) = 0.
326         ENDDO         ENDDO
327        ENDDO        ENDDO
328    
329    #ifdef ALLOW_AUTODIFF_TAMC
330    #ifndef DISABLE_MULTIDIM_ADVECTION
331    CADJ STORE localTij(:,:)  =
332    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
333    #endif
334    #endif /* ALLOW_AUTODIFF_TAMC */
335    
336        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
337         CALL GAD_FLUXLIMIT_ADV_X(         CALL GAD_FLUXLIMIT_ADV_X(
338       &      bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)       &      bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
# Line 222  C-    Advective flux in X Line 343  C-    Advective flux in X
343         CALL GAD_DST3FL_ADV_X(         CALL GAD_DST3FL_ADV_X(
344       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
345        ELSE        ELSE
346         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'         STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
347        ENDIF        ENDIF
348    
349        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
350         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
351          localTij(i,j)=localTij(i,j)-deltaTtracer*          localTij(i,j)=localTij(i,j)-deltaTtracer*
# Line 252  C--   End of X direction Line 374  C--   End of X direction
374  C--   Y direction  C--   Y direction
375        IF (calc_fluxes_Y) THEN        IF (calc_fluxes_Y) THEN
376    
 C--   Internal exchange for calculations in Y  
377        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
378         DO j=1,Oly  C--   Internal exchange for calculations in Y
379          DO i=1,Olx  C--    For cube face corners we need to duplicate the
380           localTij( 1-i , 1-j )=localTij(   j   , 1-i )  C--    j-1 and j+1 values into the null space as follows:
381           localTij( 1-i ,sNy+j)=localTij(   j   ,sNy+i)  C
382           localTij(sNx+i, 1-j )=localTij(sNx+1-j, 1-i )  C      o SW corner: copy T(0,1) into T(0,0) e.g.
383           localTij(sNx+i,sNy+j)=localTij(sNx+1-j,sNy+i)  C                      |
384          ENDDO  C                      |  x T(1,1)
385         ENDDO  C                      |
386    C      ----------------|-----------
387    C                      |
388    C         x T(0,0)<====== x T(1,0)
389    C                      |
390    C
391    C      o NW corner: copy T(    0,sNy  ) into T(    0,sNy+1) e.g.
392    C                      |
393    C         x T(0,sNy+1)<=== x T(1,sNy+1)
394    C                      |
395    C      ----------------|-----------
396    C                      |
397    C                      |   x T(1,sNy)
398    C                      |
399    C
400    C      o NE corner: copy T(sNx+1,sNy  ) into T(sNx+1,sNy+1) e.g.
401    C                      |
402    C      x T(sNx,sNy+1)=====>x T(sNx+1,sNy+1)
403    C                      |    
404    C      ----------------|-----------
405    C                      |    
406    C      x T(sNx,sNy)    |                      
407    C                      |
408    C      o SE corner: copy T(sNx+1,1    ) into T(sNx+1,0    ) e.g.
409    C                      |
410    C         x T(sNx,1)   |                    
411    C                      |    
412    C      ----------------|-----------
413    C                      |    
414    C         x T(sNx,0) =====>x T(sNx+1,    0)
415           IF ( southWestCorner ) THEN
416             localTij(    0,0    ) = localTij(  1,0    )
417           ENDIF
418           IF ( southEastCorner ) THEN
419             localTij(sNx+1,0    ) = localTij(sNx,0    )
420           ENDIF
421           IF ( northWestCorner ) THEN
422             localTij(0    ,sNy+1) = localTij(  1,sNy+1)
423           ENDIF
424           IF ( northEastCorner ) THEN
425             localTij(sNx+1,sNy+1) = localTij(sNx,sNy+1)
426           ENDIF
427        ENDIF        ENDIF
428    
429  C-    Advective flux in Y  C-    Advective flux in Y
# Line 270  C-    Advective flux in Y Line 432  C-    Advective flux in Y
432          af(i,j) = 0.          af(i,j) = 0.
433         ENDDO         ENDDO
434        ENDDO        ENDDO
435    
436    #ifdef ALLOW_AUTODIFF_TAMC
437    #ifndef DISABLE_MULTIDIM_ADVECTION
438    CADJ STORE localTij(:,:)  =
439    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
440    #endif
441    #endif /* ALLOW_AUTODIFF_TAMC */
442    
443        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
444         CALL GAD_FLUXLIMIT_ADV_Y(         CALL GAD_FLUXLIMIT_ADV_Y(
445       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
# Line 282  C-    Advective flux in Y Line 452  C-    Advective flux in Y
452        ELSE        ELSE
453         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
454        ENDIF        ENDIF
455    
456        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
457         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
458          localTij(i,j)=localTij(i,j)-deltaTtracer*          localTij(i,j)=localTij(i,j)-deltaTtracer*
# Line 307  C--   Apply open boundary conditions Line 478  C--   Apply open boundary conditions
478  C--   End of Y direction  C--   End of Y direction
479        ENDIF        ENDIF
480    
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         localTijk(i,j,k)=localTij(i,j)  
        ENDDO  
       ENDDO  
   
481  C--   End of ipass loop  C--   End of ipass loop
482        ENDDO        ENDDO
483    
484          IF ( implicitAdvection ) THEN
485    C-    explicit advection is done ; store tendency in gTracer:
486            DO j=1-Oly,sNy+Oly
487             DO i=1-Olx,sNx+Olx
488              gTracer(i,j,k,bi,bj)=
489         &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
490             ENDDO
491            ENDDO
492          ELSE
493    C-    horizontal advection done; store intermediate result in 3D array:
494           DO j=1-Oly,sNy+Oly
495            DO i=1-Olx,sNx+Olx
496             localTijk(i,j,k)=localTij(i,j)
497            ENDDO
498           ENDDO
499          ENDIF
500    
501  C--   End of K loop for horizontal fluxes  C--   End of K loop for horizontal fluxes
502        ENDDO        ENDDO
503    
504          IF ( .NOT.implicitAdvection ) THEN
505  C--   Start of k loop for vertical flux  C--   Start of k loop for vertical flux
506        DO k=Nr,1,-1         DO k=Nr,1,-1
507    #ifdef ALLOW_AUTODIFF_TAMC
508             kkey = (igadkey-1)*Nr + k
509    #endif /* ALLOW_AUTODIFF_TAMC */
510  C--   kup    Cycles through 1,2 to point to w-layer above  C--   kup    Cycles through 1,2 to point to w-layer above
511  C--   kDown  Cycles through 2,1 to point to w-layer below  C--   kDown  Cycles through 2,1 to point to w-layer below
512        kup  = 1+MOD(k+1,2)          kup  = 1+MOD(k+1,2)
513        kDown= 1+MOD(k,2)          kDown= 1+MOD(k,2)
514    c       kp1=min(Nr,k+1)
515  C--   Get temporary terms used by tendency routines          kp1Msk=1.
516        CALL CALC_COMMON_FACTORS (          if (k.EQ.Nr) kp1Msk=0.
517       I         bi,bj,iMin,iMax,jMin,jMax,k,  
518       O         xA,yA,uTrans,vTrans,rTrans,maskUp,  C-- Compute Vertical transport
519       I         myThid)  #ifdef ALLOW_AIM
520    C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
521            IF ( k.EQ.1 .OR.
522         &     (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
523         &              ) THEN
524    #else
525            IF ( k.EQ.1 ) THEN
526    #endif
527    
528    C- Surface interface :
529             DO j=1-Oly,sNy+Oly
530              DO i=1-Olx,sNx+Olx
531               rTransKp1(i,j) = kp1Msk*rTrans(i,j)
532               rTrans(i,j) = 0.
533               fVerT(i,j,kUp) = 0.
534               af(i,j) = 0.
535              ENDDO
536             ENDDO
537    
538            ELSE
539    C- Interior interface :
540    
541             DO j=1-Oly,sNy+Oly
542              DO i=1-Olx,sNx+Olx
543               rTransKp1(i,j) = kp1Msk*rTrans(i,j)
544               rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
545         &                 *maskC(i,j,k-1,bi,bj)
546               af(i,j) = 0.
547              ENDDO
548             ENDDO
549    
550    #ifdef ALLOW_GMREDI
551    C--   Residual transp = Bolus transp + Eulerian transp
552             IF (useGMRedi)
553         &   CALL GMREDI_CALC_WFLOW(
554         &                    rTrans, bi, bj, k, myThid)
555    #endif /* ALLOW_GMREDI */
556    
557    #ifdef ALLOW_AUTODIFF_TAMC
558    CADJ STORE localTijk(:,:,k)  
559    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
560    CADJ STORE rTrans(:,:)  
561    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
562    #endif /* ALLOW_AUTODIFF_TAMC */
563    
 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  
564  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
565         IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
566          CALL GAD_FLUXLIMIT_ADV_R(            CALL GAD_FLUXLIMIT_ADV_R(
567       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
568         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
569          CALL GAD_DST3_ADV_R(            CALL GAD_DST3_ADV_R(
570       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
571         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
572          CALL GAD_DST3FL_ADV_R(            CALL GAD_DST3FL_ADV_R(
573       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
574         ELSE           ELSE
575          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
576         ENDIF           ENDIF
 C-    Surface "correction" term at k>1 :  
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          af(i,j) = af(i,j)  
      &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  
      &             rTrans(i,j)*localTijk(i,j,k)  
         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  
577  C-    add the advective flux to fVerT  C-    add the advective flux to fVerT
578        DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
579         DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
580          fVerT(i,j,kUp) = af(i,j)             fVerT(i,j,kUp) = af(i,j)
581         ENDDO            ENDDO
582        ENDDO           ENDDO
583    
584  C--   Divergence of fluxes  C- end Surface/Interior if bloc
585        kp1=min(Nr,k+1)          ENDIF
586        kp1Msk=1.  
587        if (k.EQ.Nr) kp1Msk=0.  #ifdef ALLOW_AUTODIFF_TAMC
588        DO j=1-Oly,sNy+Oly  CADJ STORE rTrans(:,:)  
589         DO i=1-Olx,sNx+Olx  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
590          localTij(i,j)=localTijk(i,j,k)-deltaTtracer*  CADJ STORE rTranskp1(:,:)  
591       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
592       &    *recip_rA(i,j,bi,bj)  #endif /* ALLOW_AUTODIFF_TAMC */
593       &    *( fVerT(i,j,kUp)-fVerT(i,j,kDown)  
594       &      -tracer(i,j,k,bi,bj)*rA(i,j,bi,bj)*  C--   Divergence of vertical fluxes
595       &        (wVel(i,j,k,bi,bj)-kp1Msk*wVel(i,j,kp1,bi,bj))          DO j=1-Oly,sNy+Oly
596       &     )*rkFac           DO i=1-Olx,sNx+Olx
597          gTracer(i,j,k,bi,bj)=            localTij(i,j)=localTijk(i,j,k)-deltaTtracer*
598       &   (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer       &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
599         ENDDO       &     *recip_rA(i,j,bi,bj)
600        ENDDO       &     *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
601         &       -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
602         &      )*rkFac
603              gTracer(i,j,k,bi,bj)=
604         &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
605             ENDDO
606            ENDDO
607    
608  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
609        ENDDO         ENDDO
610    C--   end of if not.implicitAdvection block
611          ENDIF
612    
613        RETURN        RETURN
614        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22