/[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.16 by heimbach, Tue Nov 25 23:31:44 2003 UTC revision 1.26 by cnh, Wed Jul 7 20:09:42 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  
   
 #include "PACKAGES_CONFIG.h"  
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
 #ifdef ALLOW_AUTODIFF  
 # include "CPP_OPTIONS.h"  
 #endif  
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 66  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  #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  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 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)
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 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 129  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 203  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 215  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 227  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 250  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    C         x T(0,sNy+1) |
278    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            DO J=1,OLy
310             DO I=1,OLx
311              localTij(1-I, 1-J   )= localTij(1-J  ,1  )
312             ENDDO
313          ENDDO          ENDDO
314         ENDDO         ENDIF
315           IF ( southEastCorner ) THEN
316            DO J=1,OLy
317             DO I=1,OLx
318              localTij(sNx+I, 1-J )=localTij(sNx+J, I  )
319             ENDDO
320            ENDDO
321           ENDIF
322           IF ( northWestCorner ) THEN
323            DO J=1,OLy
324             DO I=1,OLx
325              localTij( 1-I ,sNy+J)=localTij( 1-J , sNy+1-I )
326             ENDDO
327            ENDDO
328           ENDIF
329           IF ( northEastCorner ) THEN
330            DO J=1,OLy
331             DO I=1,OLx
332              localTij(sNx+I,sNy+J)=localTij(sNx+J, sNy+1-I )
333             ENDDO
334            ENDDO
335           ENDIF
336        ENDIF        ENDIF
337    
338  C-    Advective flux in X  C-    Advective flux in X
# Line 315  C--   End of X direction Line 390  C--   End of X direction
390  C--   Y direction  C--   Y direction
391        IF (calc_fluxes_Y) THEN        IF (calc_fluxes_Y) THEN
392    
 C--   Internal exchange for calculations in Y  
393        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
394         DO j=1,Oly  C--   Internal exchange for calculations in Y
395          DO i=1,Olx  C--    For cube face corners we need to duplicate the
396           localTij( 1-i , 1-j )=localTij(   j   , 1-i )  C--    j-1 and j+1 values into the null space as follows:
397           localTij( 1-i ,sNy+j)=localTij(   j   ,sNy+i)  C
398           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.
399           localTij(sNx+i,sNy+j)=localTij(sNx+1-j,sNy+i)  C                      |
400    C                      |  x T(1,1)
401    C                      |
402    C      ----------------|-----------
403    C                      |
404    C         x T(0,0)<====== x T(1,0)
405    C                      |
406    C
407    C      o NW corner: copy T(    0,sNy  ) into T(    0,sNy+1) e.g.
408    C                      |
409    C         x T(0,sNy+1)<=== x T(1,sNy+1)
410    C                      |
411    C      ----------------|-----------
412    C                      |
413    C                      |   x T(1,sNy)
414    C                      |
415    C
416    C      o NE corner: copy T(sNx+1,sNy  ) into T(sNx+1,sNy+1) e.g.
417    C                      |
418    C      x T(sNx,sNy+1)=====>x T(sNx+1,sNy+1)
419    C                      |    
420    C      ----------------|-----------
421    C                      |    
422    C      x T(sNx,sNy)    |                      
423    C                      |
424    C      o SE corner: copy T(sNx+1,1    ) into T(sNx+1,0    ) e.g.
425    C                      |
426    C         x T(sNx,1)   |                    
427    C                      |    
428    C      ----------------|-----------
429    C                      |    
430    C         x T(sNx,0) =====>x T(sNx+1,    0)
431           IF ( southWestCorner ) THEN
432            DO J=1,Oly
433             DO I=1,Olx
434              localTij( 1-i , 1-j ) = localTij(j   , 1-i )
435             ENDDO
436          ENDDO          ENDDO
437         ENDDO         ENDIF
438           IF ( southEastCorner ) THEN
439            DO J=1,Oly
440             DO I=1,Olx
441              localTij(sNx+i, 1-j ) = localTij(sNx+1-j, 1-i )
442             ENDDO
443            ENDDO
444           ENDIF
445           IF ( northWestCorner ) THEN
446            DO J=1,Oly
447             DO I=1,Olx
448              localTij( 1-i ,sNy+j) = localTij(j   ,sNy+i)
449             ENDDO
450            ENDDO
451           ENDIF
452           IF ( northEastCorner ) THEN
453            DO J=1,Oly
454             DO I=1,Olx
455              localTij(sNx+i,sNy+j) = localTij(sNx+1-j,sNy+i)
456             ENDDO
457            ENDDO
458           ENDIF
459        ENDIF        ENDIF
460    
461  C-    Advective flux in Y  C-    Advective flux in Y
# Line 379  C--   Apply open boundary conditions Line 510  C--   Apply open boundary conditions
510  C--   End of Y direction  C--   End of Y direction
511        ENDIF        ENDIF
512    
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         localTijk(i,j,k)=localTij(i,j)  
        ENDDO  
       ENDDO  
   
513  C--   End of ipass loop  C--   End of ipass loop
514        ENDDO        ENDDO
515    
516          IF ( implicitAdvection ) THEN
517    C-    explicit advection is done ; store tendency in gTracer:
518            DO j=1-Oly,sNy+Oly
519             DO i=1-Olx,sNx+Olx
520              gTracer(i,j,k,bi,bj)=
521         &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
522             ENDDO
523            ENDDO
524          ELSE
525    C-    horizontal advection done; store intermediate result in 3D array:
526           DO j=1-Oly,sNy+Oly
527            DO i=1-Olx,sNx+Olx
528             localTijk(i,j,k)=localTij(i,j)
529            ENDDO
530           ENDDO
531          ENDIF
532    
533  C--   End of K loop for horizontal fluxes  C--   End of K loop for horizontal fluxes
534        ENDDO        ENDDO
535    
536          IF ( .NOT.implicitAdvection ) THEN
537  C--   Start of k loop for vertical flux  C--   Start of k loop for vertical flux
538        DO k=Nr,1,-1         DO k=Nr,1,-1
539  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
540           kkey = (igadkey-1)*Nr + k           kkey = (igadkey-1)*Nr + k
541  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
542  C--   kup    Cycles through 1,2 to point to w-layer above  C--   kup    Cycles through 1,2 to point to w-layer above
543  C--   kDown  Cycles through 2,1 to point to w-layer below  C--   kDown  Cycles through 2,1 to point to w-layer below
544        kup  = 1+MOD(k+1,2)          kup  = 1+MOD(k+1,2)
545        kDown= 1+MOD(k,2)          kDown= 1+MOD(k,2)
546  c     kp1=min(Nr,k+1)  c       kp1=min(Nr,k+1)
547        kp1Msk=1.          kp1Msk=1.
548        if (k.EQ.Nr) kp1Msk=0.          if (k.EQ.Nr) kp1Msk=0.
549    
550  C-- Compute Vertical transport  C-- Compute Vertical transport
551  C     Note: wVel needs to be masked  #ifdef ALLOW_AIM
552    C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
553            IF ( k.EQ.1 .OR.
554         &     (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
555         &              ) THEN
556    #else
557            IF ( k.EQ.1 ) THEN
558    #endif
559    
       IF (k.EQ.1) THEN  
560  C- Surface interface :  C- Surface interface :
561             DO j=1-Oly,sNy+Oly
562              DO i=1-Olx,sNx+Olx
563               rTransKp1(i,j) = kp1Msk*rTrans(i,j)
564               rTrans(i,j) = 0.
565               fVerT(i,j,kUp) = 0.
566               af(i,j) = 0.
567              ENDDO
568             ENDDO
569    
570         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.  
          af(i,j) = 0.  
         ENDDO  
        ENDDO  
   
       ELSE  
571  C- Interior interface :  C- Interior interface :
572         DO j=1-Oly,sNy+Oly  
573          DO i=1-Olx,sNx+Olx           DO j=1-Oly,sNy+Oly
574           rTransKp1(i,j) = kp1Msk*rTrans(i,j)            DO i=1-Olx,sNx+Olx
575           rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
576       &               *maskC(i,j,k-1,bi,bj)             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
577           af(i,j) = 0.       &                 *maskC(i,j,k-1,bi,bj)
578          ENDDO             af(i,j) = 0.
579         ENDDO            ENDDO
580             ENDDO
581    
582  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
583  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
584         IF (useGMRedi)           IF (useGMRedi)
585       &   CALL GMREDI_CALC_WFLOW(       &   CALL GMREDI_CALC_WFLOW(
586       &                    rTrans, bi, bj, k, myThid)       &                    rTrans, bi, bj, k, myThid)
587  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
# Line 446  CADJ &     = comlev1_bibj_k_gad, key=kke Line 594  CADJ &     = comlev1_bibj_k_gad, key=kke
594  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
595    
596  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
597         IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
598          CALL GAD_FLUXLIMIT_ADV_R(            CALL GAD_FLUXLIMIT_ADV_R(
599       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
600         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
601          CALL GAD_DST3_ADV_R(            CALL GAD_DST3_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_FLUX_LIMIT ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
604          CALL GAD_DST3FL_ADV_R(            CALL GAD_DST3FL_ADV_R(
605       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
606         ELSE           ELSE
607          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
608         ENDIF           ENDIF
609  C-    add the advective flux to fVerT  C-    add the advective flux to fVerT
610         DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
611          DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
612           fVerT(i,j,kUp) = af(i,j)             fVerT(i,j,kUp) = af(i,j)
613          ENDDO            ENDDO
614         ENDDO           ENDDO
615    
616  C- end Surface/Interior if bloc  C- end Surface/Interior if bloc
617        ENDIF          ENDIF
618    
619  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
620  CADJ STORE rTrans(:,:)    CADJ STORE rTrans(:,:)  
# Line 475  CADJ STORE rTranskp1(:,:) Line 623  CADJ STORE rTranskp1(:,:)
623  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte  CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
624  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
625    
626  C--   Divergence of fluxes  C--   Divergence of vertical fluxes
627        DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
628         DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
629          localTij(i,j)=localTijk(i,j,k)-deltaTtracer*            localTij(i,j)=localTijk(i,j,k)-deltaTtracer*
630       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
631       &    *recip_rA(i,j,bi,bj)       &     *recip_rA(i,j,bi,bj)
632       &    *( fVerT(i,j,kUp)-fVerT(i,j,kDown)       &     *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
633       &      -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))       &       -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
634       &     )*rkFac       &      )*rkFac
635          gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
636       &   (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
637         ENDDO           ENDDO
638        ENDDO          ENDDO
639    
640  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
641        ENDDO         ENDDO
642    C--   end of if not.implicitAdvection block
643          ENDIF
644    
645        RETURN        RETURN
646        END        END
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.22