/[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.12 by heimbach, Sun Mar 24 02:12:50 2002 UTC revision 1.27 by heimbach, Fri Sep 17 23:02:01 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: Generic 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 62  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  #ifdef ALLOW_AUTODIFF_TAMC
45  # include "tamc.h"  # include "tamc.h"
46  # include "tamc_keys.h"  # include "tamc_keys.h"
47    # ifdef ALLOW_PTRACERS
48    #  include "PTRACERS_SIZE.h"
49    # endif
50  #endif  #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 Tracer(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 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)
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  iMin,iMax,    :: loop range for called routines
88  C  i,j,k                :: loop indices  C  jMin,jMax     :: loop range for called routines
89  C  kup                  :: index into 2 1/2D array, toggles between 1 and 2  C  i,j,k         :: loop indices
90  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
91  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr  C  kdown         :: index into 2 1/2D array, toggles between 2 and 1
92  C  xA,yA                :: areas of X and Y face of tracer cells  C  kp1           :: =k+1 for k<Nr, =Nr for k=Nr
93  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
94  C  rTransKp1            :: vertical volume transport at interface k+1  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
95  C  af                   :: 2-D array for horizontal advective flux  C  rTrans        :: 2-D arrays of volume transports at W points
96  C  fVerT                :: 2 1/2D arrays for vertical advective flux  C  rTransKp1     :: vertical volume transport at interface k+1
97  C  localTij             :: 2-D array used as temporary local copy of tracer fld  C  af            :: 2-D array for horizontal advective flux
98  C  localTijk            :: 3-D array used as temporary local copy of tracer fld  C  fVerT         :: 2 1/2D arrays for vertical advective flux
99  C  kp1Msk               :: flag (0,1) to act as over-riding mask for W levels  C  localTij      :: 2-D array, temporary local copy of tracer fld
100  C  calc_fluxes_X        :: logical to indicate to calculate fluxes in X dir  C  localTijk     :: 3-D array, temporary local copy of tracer fld
101  C  calc_fluxes_Y        :: logical to indicate to calculate fluxes in Y dir  C  kp1Msk        :: flag (0,1) for over-riding mask for W levels
102  C  nipass               :: number of passes to make in multi-dimensional method  C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
103  C  ipass                :: number of the current pass being made  C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
104    C  nipass        :: number of passes in multi-dimensional method
105    C  ipass         :: number of the current pass being made
106    C  myTile        :: variables used to determine which cube face
107    C  nCFace        :: owns a tile for cube grid runs using
108    C                :: multi-dim advection.
109        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
111        INTEGER i,j,k,kup,kDown        INTEGER i,j,k,kup,kDown
# Line 125  C  ipass                :: number of the Line 122  C  ipass                :: number of the
122        _RL kp1Msk        _RL kp1Msk
123        LOGICAL calc_fluxes_X,calc_fluxes_Y        LOGICAL calc_fluxes_X,calc_fluxes_Y
124        INTEGER nipass,ipass        INTEGER nipass,ipass
125          INTEGER myTile, nCFace
126          LOGICAL southWestCorner
127          LOGICAL southEastCorner
128          LOGICAL northWestCorner
129          LOGICAL northEastCorner
130  CEOP  CEOP
131    
132  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
133              act0 = tracerIdentity - 1
134              max0 = maxpass
135            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
136            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
137            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
# Line 135  CEOP Line 139  CEOP
139            act3 = myThid - 1            act3 = myThid - 1
140            max3 = nTx*nTy            max3 = nTx*nTy
141            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
142            ikey = (act1 + 1) + act2*max1            igadkey = (act0 + 1)
143       &                      + act3*max1*max2       &                      + act1*max0
144       &                      + act4*max1*max2*max3       &                      + act2*max0*max1
145         &                      + act3*max0*max1*max2
146         &                      + act4*max0*max1*max2*max3
147              if (tracerIdentity.GT.maxpass) then
148                 print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
149                 STOP 'maxpass seems smaller than tracerIdentity'
150              endif
151  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
152    
153  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
# Line 166  C     uninitialised but inert locations. Line 176  C     uninitialised but inert locations.
176  C--   Start of k loop for horizontal fluxes  C--   Start of k loop for horizontal fluxes
177        DO k=1,Nr        DO k=1,Nr
178  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
179           kkey = (ikey-1)*Nr + k           kkey = (igadkey-1)*Nr + k
180  CADJ STORE tracer(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE tracer(:,:,k,bi,bj) =
181    CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte
182  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
183    
184  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
# Line 190  C--   Make local copy of tracer array Line 201  C--   Make local copy of tracer array
201         ENDDO         ENDDO
202        ENDDO        ENDDO
203    
204    cph  The following block is needed for useCubedSphereExchange only,
205    cph  but needs to be set for all cases to avoid spurious
206    cph  TAF dependencies
207           southWestCorner = .TRUE.
208           southEastCorner = .TRUE.
209           northWestCorner = .TRUE.
210           northEastCorner = .TRUE.
211    #ifdef ALLOW_EXCH2
212           myTile = W2_myTileList(bi)
213           nCFace = exch2_myFace(myTile)
214           southWestCorner = exch2_isWedge(myTile).EQ.1
215         &             .AND. exch2_isSedge(myTile).EQ.1
216           southEastCorner = exch2_isEedge(myTile).EQ.1
217         &             .AND. exch2_isSedge(myTile).EQ.1
218           northEastCorner = exch2_isEedge(myTile).EQ.1
219         &             .AND. exch2_isNedge(myTile).EQ.1
220           northWestCorner = exch2_isWedge(myTile).EQ.1
221         &             .AND. exch2_isNedge(myTile).EQ.1
222    #else
223           nCFace = bi
224    #endif
225        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
226    
227         nipass=3         nipass=3
228    #ifdef ALLOW_AUTODIFF_TAMC
229           if ( nipass.GT.maxcube )
230         &      STOP 'maxcube needs to be = 3'
231    #endif
232        ELSE        ELSE
233         nipass=1         nipass=1
234        ENDIF        ENDIF
235  cph       nipass=1  cph       nipass=1
236    
237  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
238    C--   For cube need one pass for each of red, green and blue axes.
239        DO ipass=1,nipass        DO ipass=1,nipass
240  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
241           passkey = ipass + (k-1)   *maxpass           passkey = ipass + (k-1)      *maxcube
242       &                   + (ikey-1)*maxpass*Nr       &                   + (igadkey-1)*maxcube*Nr
243           IF (nipass .GT. maxpass) THEN           IF (nipass .GT. maxpass) THEN
244            STOP 'GAD_ADVECTION: nipass > maxpass. check tamc.h'            STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
245           ENDIF           ENDIF
246  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
247    
248        IF (nipass.EQ.3) THEN        IF (nipass.EQ.3) THEN
249         calc_fluxes_X=.FALSE.         calc_fluxes_X=.FALSE.
250         calc_fluxes_Y=.FALSE.         calc_fluxes_Y=.FALSE.
251         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
252          calc_fluxes_X=.TRUE.          calc_fluxes_X=.TRUE.
253         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
254          calc_fluxes_Y=.TRUE.          calc_fluxes_Y=.TRUE.
255         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
256          calc_fluxes_Y=.TRUE.          calc_fluxes_Y=.TRUE.
257         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
258          calc_fluxes_X=.TRUE.          calc_fluxes_X=.TRUE.
259         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
260          calc_fluxes_Y=.TRUE.          calc_fluxes_Y=.TRUE.
261         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
262          calc_fluxes_X=.TRUE.          calc_fluxes_X=.TRUE.
263         ENDIF         ENDIF
264        ELSE        ELSE
# Line 233  C--   X direction Line 271  C--   X direction
271    
272  C--   Internal exchange for calculations in X  C--   Internal exchange for calculations in X
273        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
274         DO j=1,Oly  C--    For cube face corners we need to duplicate the
275          DO i=1,Olx  C--    i-1 and i+1 values into the null space as follows:
276           localTij( 1-i , 1-j )=localTij( 1-j ,    i    )  C
277           localTij( 1-i ,sNy+j)=localTij( 1-j , sNy+1-i )  C
278           localTij(sNx+i, 1-j )=localTij(sNx+j,    i    )  C      o NW corner: copy T(    0,sNy  ) into T(    0,sNy+1) e.g.
279           localTij(sNx+i,sNy+j)=localTij(sNx+j, sNy+1-i )  C                      |
280    C         x T(0,sNy+1) |
281    C        /\            |
282    C      --||------------|-----------
283    C        ||            |
284    C         x T(0,sNy)   |   x T(1,sNy)
285    C                      |
286    C
287    C      o SW corner: copy T(0,1) into T(0,0) e.g.
288    C                      |
289    C         x T(0,1)     |  x T(1,1)
290    C        ||            |
291    C      --||------------|-----------
292    C        \/            |
293    C         x T(0,0)     |
294    C                      |
295    C
296    C      o NE corner: copy T(sNx+1,sNy  ) into T(sNx+1,sNy+1) e.g.
297    C                      |
298    C                      |   x T(sNx+1,sNy+1)
299    C                      |  /\
300    C      ----------------|--||-------
301    C                      |  ||
302    C         x T(sNx,sNy) |   x T(sNx+1,sNy  )
303    C                      |
304    C      o SE corner: copy T(sNx+1,1    ) into T(sNx+1,0    ) e.g.
305    C                      |
306    C         x T(sNx,1)   |   x T(sNx+1,    1)
307    C                      |  ||
308    C      ----------------|--||-------
309    C                      |  \/
310    C                      |   x T(sNx+1,    0)
311           IF ( southWestCorner ) THEN
312            DO J=1,OLy
313             DO I=1,OLx
314              localTij(1-I, 1-J   )= localTij(1-J  ,1  )
315             ENDDO
316          ENDDO          ENDDO
317         ENDDO         ENDIF
318           IF ( southEastCorner ) THEN
319            DO J=1,OLy
320             DO I=1,OLx
321              localTij(sNx+I, 1-J )=localTij(sNx+J, I  )
322             ENDDO
323            ENDDO
324           ENDIF
325           IF ( northWestCorner ) THEN
326            DO J=1,OLy
327             DO I=1,OLx
328              localTij( 1-I ,sNy+J)=localTij( 1-J , sNy+1-I )
329             ENDDO
330            ENDDO
331           ENDIF
332           IF ( northEastCorner ) THEN
333            DO J=1,OLy
334             DO I=1,OLx
335              localTij(sNx+I,sNy+J)=localTij(sNx+J, sNy+1-I )
336             ENDDO
337            ENDDO
338           ENDIF
339        ENDIF        ENDIF
340    
341  C-    Advective flux in X  C-    Advective flux in X
# Line 252  C-    Advective flux in X Line 347  C-    Advective flux in X
347    
348  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
349  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
350  CADJ STORE localTij(:,:)  = comlev1_bibj_pass, key=passkey, byte=isbyte  CADJ STORE localTij(:,:)  =
351    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
352  #endif  #endif
353  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
354    
# Line 297  C--   End of X direction Line 393  C--   End of X direction
393  C--   Y direction  C--   Y direction
394        IF (calc_fluxes_Y) THEN        IF (calc_fluxes_Y) THEN
395    
 C--   Internal exchange for calculations in Y  
396        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
397         DO j=1,Oly  C--   Internal exchange for calculations in Y
398          DO i=1,Olx  C--    For cube face corners we need to duplicate the
399           localTij( 1-i , 1-j )=localTij(   j   , 1-i )  C--    j-1 and j+1 values into the null space as follows:
400           localTij( 1-i ,sNy+j)=localTij(   j   ,sNy+i)  C
401           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.
402           localTij(sNx+i,sNy+j)=localTij(sNx+1-j,sNy+i)  C                      |
403    C                      |  x T(1,1)
404    C                      |
405    C      ----------------|-----------
406    C                      |
407    C         x T(0,0)<====== x T(1,0)
408    C                      |
409    C
410    C      o NW corner: copy T(    0,sNy  ) into T(    0,sNy+1) e.g.
411    C                      |
412    C         x T(0,sNy+1)<=== x T(1,sNy+1)
413    C                      |
414    C      ----------------|-----------
415    C                      |
416    C                      |   x T(1,sNy)
417    C                      |
418    C
419    C      o NE corner: copy T(sNx+1,sNy  ) into T(sNx+1,sNy+1) e.g.
420    C                      |
421    C      x T(sNx,sNy+1)=====>x T(sNx+1,sNy+1)
422    C                      |    
423    C      ----------------|-----------
424    C                      |    
425    C      x T(sNx,sNy)    |                      
426    C                      |
427    C      o SE corner: copy T(sNx+1,1    ) into T(sNx+1,0    ) e.g.
428    C                      |
429    C         x T(sNx,1)   |                    
430    C                      |    
431    C      ----------------|-----------
432    C                      |    
433    C         x T(sNx,0) =====>x T(sNx+1,    0)
434           IF ( southWestCorner ) THEN
435            DO J=1,Oly
436             DO I=1,Olx
437              localTij( 1-i , 1-j ) = localTij(j   , 1-i )
438             ENDDO
439          ENDDO          ENDDO
440         ENDDO         ENDIF
441           IF ( southEastCorner ) THEN
442            DO J=1,Oly
443             DO I=1,Olx
444              localTij(sNx+i, 1-j ) = localTij(sNx+1-j, 1-i )
445             ENDDO
446            ENDDO
447           ENDIF
448           IF ( northWestCorner ) THEN
449            DO J=1,Oly
450             DO I=1,Olx
451              localTij( 1-i ,sNy+j) = localTij(j   ,sNy+i)
452             ENDDO
453            ENDDO
454           ENDIF
455           IF ( northEastCorner ) THEN
456            DO J=1,Oly
457             DO I=1,Olx
458              localTij(sNx+i,sNy+j) = localTij(sNx+1-j,sNy+i)
459             ENDDO
460            ENDDO
461           ENDIF
462        ENDIF        ENDIF
463    
464  C-    Advective flux in Y  C-    Advective flux in Y
# Line 318  C-    Advective flux in Y Line 470  C-    Advective flux in Y
470    
471  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
472  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
473  CADJ STORE localTij(:,:)  = comlev1_bibj_pass, key=passkey, byte=isbyte  CADJ STORE localTij(:,:)  =
474    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
475  #endif  #endif
476  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
477    
# Line 360  C--   Apply open boundary conditions Line 513  C--   Apply open boundary conditions
513  C--   End of Y direction  C--   End of Y direction
514        ENDIF        ENDIF
515    
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         localTijk(i,j,k)=localTij(i,j)  
        ENDDO  
       ENDDO  
   
516  C--   End of ipass loop  C--   End of ipass loop
517        ENDDO        ENDDO
518    
519          IF ( implicitAdvection ) THEN
520    C-    explicit advection is done ; store tendency in gTracer:
521            DO j=1-Oly,sNy+Oly
522             DO i=1-Olx,sNx+Olx
523              gTracer(i,j,k,bi,bj)=
524         &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
525             ENDDO
526            ENDDO
527          ELSE
528    C-    horizontal advection done; store intermediate result in 3D array:
529           DO j=1-Oly,sNy+Oly
530            DO i=1-Olx,sNx+Olx
531             localTijk(i,j,k)=localTij(i,j)
532            ENDDO
533           ENDDO
534          ENDIF
535    
536  C--   End of K loop for horizontal fluxes  C--   End of K loop for horizontal fluxes
537        ENDDO        ENDDO
538    
539          IF ( .NOT.implicitAdvection ) THEN
540  C--   Start of k loop for vertical flux  C--   Start of k loop for vertical flux
541        DO k=Nr,1,-1         DO k=Nr,1,-1
542  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
543           kkey = (ikey-1)*Nr + k           kkey = (igadkey-1)*Nr + k
544  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
545  C--   kup    Cycles through 1,2 to point to w-layer above  C--   kup    Cycles through 1,2 to point to w-layer above
546  C--   kDown  Cycles through 2,1 to point to w-layer below  C--   kDown  Cycles through 2,1 to point to w-layer below
547        kup  = 1+MOD(k+1,2)          kup  = 1+MOD(k+1,2)
548        kDown= 1+MOD(k,2)          kDown= 1+MOD(k,2)
549  c     kp1=min(Nr,k+1)  c       kp1=min(Nr,k+1)
550        kp1Msk=1.          kp1Msk=1.
551        if (k.EQ.Nr) kp1Msk=0.          if (k.EQ.Nr) kp1Msk=0.
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE localTijk(:,:,k)    
 CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE rTrans(:,:)    
 CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
552    
553  C-- Compute Vertical transport  C-- Compute Vertical transport
554  C     Note: wVel needs to be masked  #ifdef ALLOW_AIM
555    C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
556            IF ( k.EQ.1 .OR.
557         &     (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
558         &              ) THEN
559    #else
560            IF ( k.EQ.1 ) THEN
561    #endif
562    
       IF (k.EQ.1) THEN  
563  C- Surface interface :  C- Surface interface :
564             DO j=1-Oly,sNy+Oly
565              DO i=1-Olx,sNx+Olx
566               rTransKp1(i,j) = kp1Msk*rTrans(i,j)
567               rTrans(i,j) = 0.
568               fVerT(i,j,kUp) = 0.
569               af(i,j) = 0.
570              ENDDO
571             ENDDO
572    
573         DO j=1-Oly,sNy+Oly          ELSE
         DO i=1-Olx,sNx+Olx  
          rTransKp1(i,j) = rTrans(i,j)  
          rTrans(i,j) = 0.  
          fVerT(i,j,kUp) = 0.  
         ENDDO  
        ENDDO  
   
       ELSE  
574  C- Interior interface :  C- Interior interface :
575         DO j=1-Oly,sNy+Oly  
576          DO i=1-Olx,sNx+Olx           DO j=1-Oly,sNy+Oly
577           rTransKp1(i,j) = kp1Msk*rTrans(i,j)            DO i=1-Olx,sNx+Olx
578           rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
579       &               *maskC(i,j,k-1,bi,bj)             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
580           af(i,j) = 0.       &                 *maskC(i,j,k-1,bi,bj)
581          ENDDO             af(i,j) = 0.
582         ENDDO            ENDDO
583             ENDDO
584    
585  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
586  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
587         IF (useGMRedi)           IF (useGMRedi)
588       &   CALL GMREDI_CALC_WFLOW(       &   CALL GMREDI_CALC_WFLOW(
589       &                    rTrans, bi, bj, k, myThid)       &                    rTrans, bi, bj, k, myThid)
590  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
591    
592    #ifdef ALLOW_AUTODIFF_TAMC
593    CADJ STORE localTijk(:,:,k)  
594    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
595    CADJ STORE rTrans(:,:)  
596    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
597    #endif /* ALLOW_AUTODIFF_TAMC */
598    
599  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
600         IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
601          CALL GAD_FLUXLIMIT_ADV_R(            CALL GAD_FLUXLIMIT_ADV_R(
602       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
603         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
604          CALL GAD_DST3_ADV_R(            CALL GAD_DST3_ADV_R(
605       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
606         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
607          CALL GAD_DST3FL_ADV_R(            CALL GAD_DST3FL_ADV_R(
608       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
609         ELSE           ELSE
610          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
611         ENDIF           ENDIF
612  C-    add the advective flux to fVerT  C-    add the advective flux to fVerT
613         DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
614          DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
615           fVerT(i,j,kUp) = af(i,j)             fVerT(i,j,kUp) = af(i,j)
616          ENDDO            ENDDO
617         ENDDO           ENDDO
618    
619  C- end Surface/Interior if bloc  C- end Surface/Interior if bloc
620        ENDIF          ENDIF
621    
622  C--   Divergence of fluxes  #ifdef ALLOW_AUTODIFF_TAMC
623        DO j=1-Oly,sNy+Oly  CADJ STORE rTrans(:,:)  
624         DO i=1-Olx,sNx+Olx  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
625          localTij(i,j)=localTijk(i,j,k)-deltaTtracer*  CADJ STORE rTranskp1(:,:)  
626       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
627       &    *recip_rA(i,j,bi,bj)  #endif /* ALLOW_AUTODIFF_TAMC */
628       &    *( fVerT(i,j,kUp)-fVerT(i,j,kDown)  
629       &      -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))  C--   Divergence of vertical fluxes
630       &     )*rkFac          DO j=1-Oly,sNy+Oly
631          gTracer(i,j,k,bi,bj)=           DO i=1-Olx,sNx+Olx
632       &   (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer            localTij(i,j)=localTijk(i,j,k)-deltaTtracer*
633         ENDDO       &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
634        ENDDO       &     *recip_rA(i,j,bi,bj)
635         &     *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
636         &       -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
637         &      )*rkFac
638              gTracer(i,j,k,bi,bj)=
639         &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
640             ENDDO
641            ENDDO
642    
643  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
644        ENDDO         ENDDO
645    C--   end of if not.implicitAdvection block
646          ENDIF
647    
648        RETURN        RETURN
649        END        END
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22