/[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.19 by edhill, Sat Mar 27 03:51:51 2004 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-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 36  C !ROUTINE: GAD_ADVECTION Line 9  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         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 70  C !USES: =============================== Line 44  C !USES: ===============================
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  implicitAdvection    :: vertical advection treated implicitly (later on)  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  uVel                 :: velocity, zonal component  C  tracerIdentity    :: tracer identifier (required only for OBCS)
61  C  vVel                 :: velocity, meridional component  C  uVel              :: velocity, zonal component
62  C  wVel                 :: velocity, vertical component  C  vVel              :: velocity, meridional component
63  C  tracer               :: tracer field  C  wVel              :: velocity, vertical component
64  C  bi,bj                :: tile indices  C  tracer            :: tracer field
65  C  myTime               :: current time  C  bi,bj             :: tile indices
66  C  myIter               :: iteration number  C  myTime            :: current time
67  C  myThid               :: thread number  C  myIter            :: iteration number
68    C  myThid            :: thread number
69        LOGICAL implicitAdvection        LOGICAL implicitAdvection
70        INTEGER advectionScheme        INTEGER advectionScheme, vertAdvecScheme
71        INTEGER tracerIdentity        INTEGER tracerIdentity
72        _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)
73        _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 79  C  myThid               :: thread number
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 135  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
# Line 209  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  #ifdef ALLOW_AUTODIFF_TAMC
229         if ( nipass.GT.maxcube )         if ( nipass.GT.maxcube )
# Line 221  C--   Make local copy of tracer array Line 235  C--   Make local copy of tracer array
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)      *maxcube           passkey = ipass + (k-1)      *maxcube
# Line 233  C--   Multiple passes for different dire Line 248  C--   Multiple passes for different dire
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 256  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 321  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 423  c       kp1=min(Nr,k+1) Line 551  c       kp1=min(Nr,k+1)
551          if (k.EQ.Nr) kp1Msk=0.          if (k.EQ.Nr) kp1Msk=0.
552    
553  C-- Compute Vertical transport  C-- Compute Vertical transport
554          IF (k.EQ.1) THEN  #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    
563  C- Surface interface :  C- Surface interface :
564           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
565            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
566             rTransKp1(i,j) = rTrans(i,j)             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
567             rTrans(i,j) = 0.             rTrans(i,j) = 0.
568             fVerT(i,j,kUp) = 0.             fVerT(i,j,kUp) = 0.
569             af(i,j) = 0.             af(i,j) = 0.
# Line 462  CADJ &     = comlev1_bibj_k_gad, key=kke Line 597  CADJ &     = comlev1_bibj_k_gad, key=kke
597  #endif /* ALLOW_AUTODIFF_TAMC */  #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
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22