/[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.18 by jmc, Wed Jan 7 21:35:00 2004 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: 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(        SUBROUTINE GAD_ADVECTION(
12       I           implicitAdvection, advectionScheme, tracerIdentity,       I     implicitAdvection, advectionScheme, vertAdvecScheme,
13       I           uVel, vVel, wVel, tracer,       I     tracerIdentity,
14       O           gTracer,       I     uVel, vVel, wVel, tracer,
15       I           bi,bj, myTime,myIter,myThid)       O     gTracer,
16  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|       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 71  C !USES: =============================== Line 45  C !USES: ===============================
45  # include "tamc.h"  # include "tamc.h"
46  # include "tamc_keys.h"  # include "tamc_keys.h"
47  #endif  #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  implicitAdvection    :: vertical advection treated implicitly (later on)  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  uVel                 :: velocity, zonal component  C  tracerIdentity    :: tracer identifier (required only for OBCS)
58  C  vVel                 :: velocity, meridional component  C  uVel              :: velocity, zonal component
59  C  wVel                 :: velocity, vertical component  C  vVel              :: velocity, meridional component
60  C  tracer               :: tracer field  C  wVel              :: velocity, vertical component
61  C  bi,bj                :: tile indices  C  tracer            :: tracer field
62  C  myTime               :: current time  C  bi,bj             :: tile indices
63  C  myIter               :: iteration number  C  myTime            :: current time
64  C  myThid               :: thread number  C  myIter            :: iteration number
65    C  myThid            :: thread number
66        LOGICAL implicitAdvection        LOGICAL implicitAdvection
67        INTEGER advectionScheme        INTEGER advectionScheme, vertAdvecScheme
68        INTEGER tracerIdentity        INTEGER tracerIdentity
69        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
70        _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
# Line 97  C  myThid               :: thread number Line 76  C  myThid               :: thread number
76        INTEGER myThid        INTEGER myThid
77    
78  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
79  C  gTracer              :: tendancy array  C  gTracer           :: tendancy array
80        _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)
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  rTransKp1            :: vertical volume transport at interface k+1  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
92  C  af                   :: 2-D array for horizontal advective flux  C  rTrans        :: 2-D arrays of volume transports at W points
93  C  fVerT                :: 2 1/2D arrays for vertical advective flux  C  rTransKp1     :: vertical volume transport at interface k+1
94  C  localTij             :: 2-D array used as temporary local copy of tracer fld  C  af            :: 2-D array for horizontal advective flux
95  C  localTijk            :: 3-D array used as temporary local copy of tracer fld  C  fVerT         :: 2 1/2D arrays for vertical advective flux
96  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
97  C  calc_fluxes_X        :: logical to indicate to calculate fluxes in X dir  C  localTijk     :: 3-D array, temporary local copy of tracer fld
98  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
99  C  nipass               :: number of passes to make in multi-dimensional method  C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
100  C  ipass                :: number of the current pass being made  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        INTEGER i,j,k,kup,kDown
# Line 135  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  #ifdef ALLOW_AUTODIFF_TAMC
# Line 209  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  #ifdef ALLOW_AUTODIFF_TAMC
226         if ( nipass.GT.maxcube )         if ( nipass.GT.maxcube )
# Line 221  C--   Make local copy of tracer array Line 232  C--   Make local copy of tracer array
232  cph       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  #ifdef ALLOW_AUTODIFF_TAMC
238           passkey = ipass + (k-1)      *maxcube           passkey = ipass + (k-1)      *maxcube
# Line 233  C--   Multiple passes for different dire Line 245  C--   Multiple passes for different dire
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 256  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 321  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 423  c       kp1=min(Nr,k+1) Line 516  c       kp1=min(Nr,k+1)
516          if (k.EQ.Nr) kp1Msk=0.          if (k.EQ.Nr) kp1Msk=0.
517    
518  C-- Compute Vertical transport  C-- Compute Vertical transport
519          IF (k.EQ.1) THEN  #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 :  C- Surface interface :
529           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
530            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
531             rTransKp1(i,j) = rTrans(i,j)             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
532             rTrans(i,j) = 0.             rTrans(i,j) = 0.
533             fVerT(i,j,kUp) = 0.             fVerT(i,j,kUp) = 0.
534             af(i,j) = 0.             af(i,j) = 0.
# Line 462  CADJ &     = comlev1_bibj_k_gad, key=kke Line 562  CADJ &     = comlev1_bibj_k_gad, key=kke
562  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
563    
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
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22