/[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.24 by dimitri, Mon Jun 28 21:10:55 2004 UTC revision 1.33 by jmc, Thu Dec 16 22:28:43 2004 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    #undef MULTIDIM_OLD_VERSION
6    
7  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
8  CBOP  CBOP
# Line 44  C !USES: =============================== Line 45  C !USES: ===============================
45  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
46  # include "tamc.h"  # include "tamc.h"
47  # include "tamc_keys.h"  # include "tamc_keys.h"
48    # ifdef ALLOW_PTRACERS
49    #  include "PTRACERS_SIZE.h"
50    # endif
51  #endif  #endif
52  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
53  #include "W2_EXCH2_TOPOLOGY.h"  #include "W2_EXCH2_TOPOLOGY.h"
# Line 81  C  gTracer           :: tendancy array Line 85  C  gTracer           :: tendancy array
85    
86  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
87  C  maskUp        :: 2-D array for mask at W points  C  maskUp        :: 2-D array for mask at W points
88    C  maskLocW      :: 2-D array for mask at West points
89    C  maskLocS      :: 2-D array for mask at South points
90  C  iMin,iMax,    :: loop range for called routines  C  iMin,iMax,    :: loop range for called routines
91  C  jMin,jMax     :: loop range for called routines  C  jMin,jMax     :: loop range for called routines
92    C [iMin,iMax]Upd :: loop range to update tracer field
93    C [jMin,jMax]Upd :: loop range to update tracer field
94  C  i,j,k         :: loop indices  C  i,j,k         :: loop indices
95  C  kup           :: index into 2 1/2D array, toggles between 1 and 2  C  kup           :: index into 2 1/2D array, toggles between 1 and 2
96  C  kdown         :: index into 2 1/2D array, toggles between 2 and 1  C  kdown         :: index into 2 1/2D array, toggles between 2 and 1
# Line 92  C  uTrans,vTrans :: 2-D arrays of volume Line 100  C  uTrans,vTrans :: 2-D arrays of volume
100  C  rTrans        :: 2-D arrays of volume transports at W points  C  rTrans        :: 2-D arrays of volume transports at W points
101  C  rTransKp1     :: vertical volume transport at interface k+1  C  rTransKp1     :: vertical volume transport at interface k+1
102  C  af            :: 2-D array for horizontal advective flux  C  af            :: 2-D array for horizontal advective flux
103    C  afx           :: 2-D array for horizontal advective flux, x direction
104    C  afy           :: 2-D array for horizontal advective flux, y direction
105  C  fVerT         :: 2 1/2D arrays for vertical advective flux  C  fVerT         :: 2 1/2D arrays for vertical advective flux
106  C  localTij      :: 2-D array, temporary local copy of tracer fld  C  localTij      :: 2-D array, temporary local copy of tracer fld
107  C  localTijk     :: 3-D array, temporary local copy of tracer fld  C  localTijk     :: 3-D array, temporary local copy of tracer fld
108  C  kp1Msk        :: flag (0,1) for over-riding mask for W levels  C  kp1Msk        :: flag (0,1) for over-riding mask for W levels
109  C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir  C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
110  C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir  C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
111    C  interiorOnly  :: only update the interior of myTile, but not the edges
112    C  overlapOnly   :: only update the edges of myTile, but not the interior
113  C  nipass        :: number of passes in multi-dimensional method  C  nipass        :: number of passes in multi-dimensional method
114  C  ipass         :: number of the current pass being made  C  ipass         :: number of the current pass being made
115  C  myTile        :: variables used to determine which cube face  C  myTile        :: variables used to determine which cube face
116  C  nCFace        :: owns a tile for cube grid runs using  C  nCFace        :: owns a tile for cube grid runs using
117  C                :: multi-dim advection.  C                :: multi-dim advection.
118    C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
119        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120          _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121          _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
123          INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
124        INTEGER i,j,k,kup,kDown        INTEGER i,j,k,kup,kDown
125        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 113  C                :: multi-dim advection. Line 129  C                :: multi-dim advection.
129        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132          _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133          _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
135        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136        _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
137        _RL kp1Msk        _RL kp1Msk
138        LOGICAL calc_fluxes_X,calc_fluxes_Y        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
139          LOGICAL interiorOnly, overlapOnly
140        INTEGER nipass,ipass        INTEGER nipass,ipass
141        INTEGER myTile, nCFace        INTEGER myTile, nCFace
142        LOGICAL southWestCorner        LOGICAL N_edge, S_edge, E_edge, W_edge
143        LOGICAL southEastCorner  #ifdef ALLOW_DIAGNOSTICS
144        LOGICAL northWestCorner        INTEGER     kk
145        LOGICAL northEastCorner        CHARACTER*8 diagName
146          CHARACTER*4 GAD_DIAG_SUFX, diagSufx
147          EXTERNAL    GAD_DIAG_SUFX
148    #endif
149  CEOP  CEOP
150    
151  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 147  CEOP Line 169  CEOP
169            endif            endif
170  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
171    
172    #ifdef ALLOW_DIAGNOSTICS
173    C--   Set diagnostic suffix for the current tracer
174          IF ( useDiagnostics ) THEN
175            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
176          ENDIF
177    #endif
178    
179  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
180  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
181  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
# Line 165  C     uninitialised but inert locations. Line 194  C     uninitialised but inert locations.
194         ENDDO         ENDDO
195        ENDDO        ENDDO
196    
197    C--   Set tile-specific parameters for horizontal fluxes
198          IF (useCubedSphereExchange) THEN
199           nipass=3
200    #ifdef ALLOW_AUTODIFF_TAMC
201           IF ( nipass.GT.maxcube ) STOP 'maxcube needs to be = 3'
202    #endif
203    #ifdef ALLOW_EXCH2
204           myTile = W2_myTileList(bi)
205           nCFace = exch2_myFace(myTile)
206           N_edge = exch2_isNedge(myTile).EQ.1
207           S_edge = exch2_isSedge(myTile).EQ.1
208           E_edge = exch2_isEedge(myTile).EQ.1
209           W_edge = exch2_isWedge(myTile).EQ.1
210    #else
211           nCFace = bi
212           N_edge = .TRUE.
213           S_edge = .TRUE.
214           E_edge = .TRUE.
215           W_edge = .TRUE.
216    #endif
217          ELSE
218           nipass=2
219           N_edge = .FALSE.
220           S_edge = .FALSE.
221           E_edge = .FALSE.
222           W_edge = .FALSE.
223          ENDIF
224    
225        iMin = 1-OLx        iMin = 1-OLx
226        iMax = sNx+OLx        iMax = sNx+OLx
227        jMin = 1-OLy        jMin = 1-OLy
# Line 186  C--   Get temporary terms used by tenden Line 243  C--   Get temporary terms used by tenden
243    
244  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
245  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
246         IF (useGMRedi)        IF (useGMRedi)
247       &   CALL GMREDI_CALC_UVFLOW(       &   CALL GMREDI_CALC_UVFLOW(
248       &            uTrans, vTrans, bi, bj, k, myThid)       &            uTrans, vTrans, bi, bj, k, myThid)
249  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
250    
251  C--   Make local copy of tracer array  C--   Make local copy of tracer array and mask West & South
252        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
253         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
254          localTij(i,j)=tracer(i,j,k,bi,bj)           localTij(i,j)=tracer(i,j,k,bi,bj)
255             maskLocW(i,j)=maskW(i,j,k,bi,bj)
256             maskLocS(i,j)=maskS(i,j,k,bi,bj)
257         ENDDO         ENDDO
258        ENDDO        ENDDO
259    
260    #ifndef ALLOW_AUTODIFF_TAMC
261        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
262         southWestCorner = .TRUE.          withSigns = .FALSE.
263         southEastCorner = .TRUE.          CALL FILL_CS_CORNER_UV_RS(
264         northWestCorner = .TRUE.       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )
        northEastCorner = .TRUE.  
 #ifdef ALLOW_EXCH2  
        myTile = W2_myTileList(bi)  
        nCFace = exch2_myFace(myTile)  
        southWestCorner = exch2_isWedge(myTile).EQ.1  
      &             .AND. exch2_isSedge(myTile).EQ.1  
        southEastCorner = exch2_isEedge(myTile).EQ.1  
      &             .AND. exch2_isSedge(myTile).EQ.1  
        northEastCorner = exch2_isEedge(myTile).EQ.1  
      &             .AND. exch2_isNedge(myTile).EQ.1  
        northWestCorner = exch2_isWedge(myTile).EQ.1  
      &             .AND. exch2_isNedge(myTile).EQ.1  
 #else  
        nCFace = bi  
 #endif  
   
        nipass=3  
 #ifdef ALLOW_AUTODIFF_TAMC  
        if ( nipass.GT.maxcube )  
      &      STOP 'maxcube needs to be = 3'  
 #endif  
       ELSE  
        nipass=1  
265        ENDIF        ENDIF
266  cph       nipass=1  #endif
267    
268  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
269  C--   For cube need one pass for each of red, green and blue axes.  C--   For cube need one pass for each of red, green and blue axes.
# Line 239  C--   For cube need one pass for each of Line 276  C--   For cube need one pass for each of
276           ENDIF           ENDIF
277  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
278    
279        IF (nipass.EQ.3) THEN        interiorOnly = .FALSE.
280         calc_fluxes_X=.FALSE.        overlapOnly  = .FALSE.
281         calc_fluxes_Y=.FALSE.        IF (useCubedSphereExchange) THEN
282         IF (ipass.EQ.1 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.2) ) THEN  #ifdef MULTIDIM_OLD_VERSION
283          calc_fluxes_X=.TRUE.  C-    CubedSphere : pass 3 times, with full update of local tracer field
284         ELSEIF (ipass.EQ.1 .AND. (nCFace.EQ.4 .OR. nCFace.EQ.5) ) THEN         IF (ipass.EQ.1) THEN
285          calc_fluxes_Y=.TRUE.          calc_fluxes_X = nCFace.EQ.1 .OR. nCFace.EQ.2
286         ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.1 .OR. nCFace.EQ.6) ) THEN          calc_fluxes_Y = nCFace.EQ.4 .OR. nCFace.EQ.5
287          calc_fluxes_Y=.TRUE.         ELSEIF (ipass.EQ.2) THEN
288         ELSEIF (ipass.EQ.2 .AND. (nCFace.EQ.3 .OR. nCFace.EQ.4) ) THEN          calc_fluxes_X = nCFace.EQ.3 .OR. nCFace.EQ.4
289          calc_fluxes_X=.TRUE.          calc_fluxes_Y = nCFace.EQ.6 .OR. nCFace.EQ.1
290         ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.2 .OR. nCFace.EQ.3) ) THEN  #else /* MULTIDIM_OLD_VERSION */
291          calc_fluxes_Y=.TRUE.  C-    CubedSphere : pass 3 times, with partial update of local tracer field
292         ELSEIF (ipass.EQ.3 .AND. (nCFace.EQ.5 .OR. nCFace.EQ.6) ) THEN         IF (ipass.EQ.1) THEN
293          calc_fluxes_X=.TRUE.          overlapOnly  = MOD(nCFace,3).EQ.0
294            interiorOnly = MOD(nCFace,3).NE.0
295            calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
296            calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
297           ELSEIF (ipass.EQ.2) THEN
298            overlapOnly  = MOD(nCFace,3).EQ.2
299            calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
300            calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
301    #endif /* MULTIDIM_OLD_VERSION */
302           ELSE
303            calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
304            calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
305         ENDIF         ENDIF
306        ELSE        ELSE
307         calc_fluxes_X=.TRUE.  C-    not CubedSphere
308         calc_fluxes_Y=.TRUE.          calc_fluxes_X = MOD(ipass,2).EQ.1
309            calc_fluxes_Y = .NOT.calc_fluxes_X
310        ENDIF        ENDIF
311    
312    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
313  C--   X direction  C--   X direction
314        IF (calc_fluxes_X) THEN        IF (calc_fluxes_X) THEN
315    
316  C--   Internal exchange for calculations in X  C-     Do not compute fluxes if
317        IF (useCubedSphereExchange) THEN  C       a) needed in overlap only
318  C--    For cube face corners we need to duplicate the  C   and b) the overlap of myTile are not cube-face Edges
319  C--    i-1 and i+1 values into the null space as follows:         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
320  C  
321  C  #ifndef ALLOW_AUTODIFF_TAMC
322  C      o NW corner: copy T(    0,sNy  ) into T(    0,sNy+1) e.g.  C-     Internal exchange for calculations in X
323  C                      |  #ifdef MULTIDIM_OLD_VERSION
324  C         x T(0,sNy+1) |          IF ( useCubedSphereExchange ) THEN
325  C        /\            |  #else
326  C      --||------------|-----------          IF ( useCubedSphereExchange .AND.
327  C        ||            |       &       ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
328  C         x T(0,sNy)   |   x T(1,sNy)  #endif
329  C                      |           CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )
330  C          ENDIF
331  C      o SW corner: copy T(0,1) into T(0,0) e.g.  #endif
 C                      |  
 C         x T(0,1)     |  x T(1,1)  
 C        ||            |  
 C      --||------------|-----------  
 C        \/            |  
 C         x T(0,0)     |  
 C                      |  
 C  
 C      o NE corner: copy T(sNx+1,sNy  ) into T(sNx+1,sNy+1) e.g.  
 C                      |  
 C                      |   x T(sNx+1,sNy+1)  
 C                      |  /\  
 C      ----------------|--||-------  
 C                      |  ||  
 C         x T(sNx,sNy) |   x T(sNx+1,sNy  )  
 C                      |  
 C      o SE corner: copy T(sNx+1,1    ) into T(sNx+1,0    ) e.g.  
 C                      |  
 C         x T(sNx,1)   |   x T(sNx+1,    1)  
 C                      |  ||  
 C      ----------------|--||-------  
 C                      |  \/  
 C                      |   x T(sNx+1,    0)  
        IF ( southWestCorner ) THEN  
         localTij(0    ,0    )= localTij(0    ,1  )  
        ENDIF  
        IF ( southEastCorner ) THEN  
         localTij(sNx+1,0    )= localTij(sNx+1,1  )  
        ENDIF  
        IF ( northWestCorner ) THEN  
         localTij(0    ,sNy+1)= localTij(0    ,sNy)  
        ENDIF  
        IF ( northEastCorner ) THEN  
         localTij(sNx+1,sNy+1)= localTij(sNx+1,sNy)  
        ENDIF  
       ENDIF  
332    
333  C-    Advective flux in X  C-     Advective flux in X
334        DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
335         DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
336          af(i,j) = 0.            af(i,j) = 0.
337         ENDDO           ENDDO
338        ENDDO          ENDDO
339    
340  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
341  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
# Line 330  CADJ &     comlev1_bibj_k_gad_pass, key= Line 344  CADJ &     comlev1_bibj_k_gad_pass, key=
344  #endif  #endif
345  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
346    
347        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
348         CALL GAD_FLUXLIMIT_ADV_X(            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
349       &      bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)       I                              uTrans, uVel, maskLocW, localTij,
350        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       O                              af, myThid )
351         CALL GAD_DST3_ADV_X(          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
352       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)            CALL GAD_DST3_ADV_X(      bi,bj,k, dTtracerLev(k),
353        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       I                              uTrans, uVel, maskLocW, localTij,
354         CALL GAD_DST3FL_ADV_X(       O                              af, myThid )
355       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
356        ELSE            CALL GAD_DST3FL_ADV_X(    bi,bj,k, dTtracerLev(k),
357         STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'       I                              uTrans, uVel, maskLocW, localTij,
358        ENDIF       O                              af, myThid )
359            ELSE
360             STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
361            ENDIF
362    
363        DO j=1-Oly,sNy+Oly  C-     Advective flux in X : done
364         DO i=1-Olx,sNx+Olx-1         ENDIF
365          localTij(i,j)=localTij(i,j)-deltaTtracer*  
366       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  #ifndef ALLOW_AUTODIFF_TAMC
367       &    *recip_rA(i,j,bi,bj)  C-     Internal exchange for next calculations in Y
368       &    *( af(i+1,j)-af(i,j)         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
369       &      -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))           CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )
370       &     )         ENDIF
371         ENDDO  #endif
372        ENDDO  
373    C-     Update the local tracer field where needed:
374    
375    C      update in overlap-Only
376           IF ( overlapOnly ) THEN
377            iMinUpd = 1-Olx+1
378            iMaxUpd = sNx+Olx-1
379    C- notes: these 2 lines below have no real effect (because recip_hFac=0
380    C         in corner region) but safer to keep them.
381            IF ( W_edge ) iMinUpd = 1
382            IF ( E_edge ) iMaxUpd = sNx
383    
384            IF ( S_edge ) THEN
385             DO j=1-Oly,0
386              DO i=iMinUpd,iMaxUpd
387               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
388         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
389         &       *recip_rA(i,j,bi,bj)
390         &       *( af(i+1,j)-af(i,j)
391         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
392         &        )
393              ENDDO
394             ENDDO
395            ENDIF
396            IF ( N_edge ) THEN
397             DO j=sNy+1,sNy+Oly
398              DO i=iMinUpd,iMaxUpd
399               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
400         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
401         &       *recip_rA(i,j,bi,bj)
402         &       *( af(i+1,j)-af(i,j)
403         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
404         &        )
405              ENDDO
406             ENDDO
407            ENDIF
408    
409           ELSE
410    C      do not only update the overlap
411            jMinUpd = 1-Oly
412            jMaxUpd = sNy+Oly
413            IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
414            IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
415            DO j=jMinUpd,jMaxUpd
416             DO i=1-Olx+1,sNx+Olx-1
417               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
418         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
419         &       *recip_rA(i,j,bi,bj)
420         &       *( af(i+1,j)-af(i,j)
421         &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
422         &        )
423             ENDDO
424            ENDDO
425    C-      keep advective flux (for diagnostics)
426            DO j=1-Oly,sNy+Oly
427             DO i=1-Olx,sNx+Olx
428              afx(i,j) = af(i,j)
429             ENDDO
430            ENDDO
431    
432  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
433  C--   Apply open boundary conditions  C-     Apply open boundary conditions
434        IF (useOBCS) THEN          IF ( useOBCS ) THEN
435         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN           IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
436          CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )            CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
437         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN           ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
438          CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )            CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
439         END IF           ENDIF
440        END IF          ENDIF
441  #endif /* ALLOW_OBCS */  #endif /* ALLOW_OBCS */
442    
443    C-     end if/else update overlap-Only
444           ENDIF
445            
446  C--   End of X direction  C--   End of X direction
447        ENDIF        ENDIF
448    
449    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
450  C--   Y direction  C--   Y direction
451        IF (calc_fluxes_Y) THEN        IF (calc_fluxes_Y) THEN
452    
453        IF (useCubedSphereExchange) THEN  C-     Do not compute fluxes if
454  C--   Internal exchange for calculations in Y  C       a) needed in overlap only
455  C--    For cube face corners we need to duplicate the  C   and b) the overlap of myTile are not cube-face edges
456  C--    j-1 and j+1 values into the null space as follows:         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
457  C  
458  C      o SW corner: copy T(0,1) into T(0,0) e.g.  #ifndef ALLOW_AUTODIFF_TAMC
459  C                      |  C-     Internal exchange for calculations in Y
460  C                      |  x T(1,1)  #ifdef MULTIDIM_OLD_VERSION
461  C                      |          IF ( useCubedSphereExchange ) THEN
462  C      ----------------|-----------  #else
463  C                      |          IF ( useCubedSphereExchange .AND.
464  C         x T(0,0)<====== x T(1,0)       &       ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
465  C                      |  #endif
466  C           CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )
467  C      o NW corner: copy T(    0,sNy  ) into T(    0,sNy+1) e.g.          ENDIF
468  C                      |  #endif
 C         x T(0,sNy+1)<=== x T(1,sNy+1)  
 C                      |  
 C      ----------------|-----------  
 C                      |  
 C                      |   x T(1,sNy)  
 C                      |  
 C  
 C      o NE corner: copy T(sNx+1,sNy  ) into T(sNx+1,sNy+1) e.g.  
 C                      |  
 C      x T(sNx,sNy+1)=====>x T(sNx+1,sNy+1)  
 C                      |      
 C      ----------------|-----------  
 C                      |      
 C      x T(sNx,sNy)    |                        
 C                      |  
 C      o SE corner: copy T(sNx+1,1    ) into T(sNx+1,0    ) e.g.  
 C                      |  
 C         x T(sNx,1)   |                      
 C                      |      
 C      ----------------|-----------  
 C                      |      
 C         x T(sNx,0) =====>x T(sNx+1,    0)  
        IF ( southWestCorner ) THEN  
          localTij(    0,0    ) = localTij(  1,0    )  
        ENDIF  
        IF ( southEastCorner ) THEN  
          localTij(sNx+1,0    ) = localTij(sNx,0    )  
        ENDIF  
        IF ( northWestCorner ) THEN  
          localTij(0    ,sNy+1) = localTij(  1,sNy+1)  
        ENDIF  
        IF ( northEastCorner ) THEN  
          localTij(sNx+1,sNy+1) = localTij(sNx,sNy+1)  
        ENDIF  
       ENDIF  
469    
470  C-    Advective flux in Y  C-     Advective flux in Y
471        DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
472         DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
473          af(i,j) = 0.            af(i,j) = 0.
474         ENDDO           ENDDO
475        ENDDO          ENDDO
476    
477  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
478  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
# Line 437  CADJ &     comlev1_bibj_k_gad_pass, key= Line 481  CADJ &     comlev1_bibj_k_gad_pass, key=
481  #endif  #endif
482  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
483    
484        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
485         CALL GAD_FLUXLIMIT_ADV_Y(            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
486       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)       I                              vTrans, vVel, maskLocS, localTij,
487        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       O                              af, myThid )
488         CALL GAD_DST3_ADV_Y(          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
489       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)            CALL GAD_DST3_ADV_Y(      bi,bj,k, dTtracerLev(k),
490        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       I                              vTrans, vVel, maskLocS, localTij,
491         CALL GAD_DST3FL_ADV_Y(       O                              af, myThid )
492       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
493        ELSE            CALL GAD_DST3FL_ADV_Y(    bi,bj,k, dTtracerLev(k),
494         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'       I                              vTrans, vVel, maskLocS, localTij,
495        ENDIF       O                              af, myThid )
496            ELSE
497             STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
498            ENDIF
499    
500        DO j=1-Oly,sNy+Oly-1  C-     Advective flux in Y : done
501         DO i=1-Olx,sNx+Olx         ENDIF
502          localTij(i,j)=localTij(i,j)-deltaTtracer*  
503       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  #ifndef ALLOW_AUTODIFF_TAMC
504       &    *recip_rA(i,j,bi,bj)  C-     Internal exchange for next calculations in X
505       &    *( af(i,j+1)-af(i,j)         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
506       &      -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))           CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )
507       &     )         ENDIF
508         ENDDO  #endif
509        ENDDO  
510    C-     Update the local tracer field where needed:
511    
512    C      update in overlap-Only
513           IF ( overlapOnly ) THEN
514            jMinUpd = 1-Oly+1
515            jMaxUpd = sNy+Oly-1
516    C- notes: these 2 lines below have no real effect (because recip_hFac=0
517    C         in corner region) but safer to keep them.
518            IF ( S_edge ) jMinUpd = 1
519            IF ( N_edge ) jMaxUpd = sNy
520    
521            IF ( W_edge ) THEN
522             DO j=jMinUpd,jMaxUpd
523              DO i=1-Olx,0
524               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
525         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
526         &       *recip_rA(i,j,bi,bj)
527         &       *( af(i,j+1)-af(i,j)
528         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
529         &        )
530              ENDDO
531             ENDDO
532            ENDIF
533            IF ( E_edge ) THEN
534             DO j=jMinUpd,jMaxUpd
535              DO i=sNx+1,sNx+Olx
536               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
537         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
538         &       *recip_rA(i,j,bi,bj)
539         &       *( af(i,j+1)-af(i,j)
540         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
541         &        )
542              ENDDO
543             ENDDO
544            ENDIF
545    
546           ELSE
547    C      do not only update the overlap
548            iMinUpd = 1-Olx
549            iMaxUpd = sNx+Olx
550            IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
551            IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
552            DO j=1-Oly+1,sNy+Oly-1
553             DO i=iMinUpd,iMaxUpd
554               localTij(i,j)=localTij(i,j)-dTtracerLev(k)*
555         &       _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
556         &       *recip_rA(i,j,bi,bj)
557         &       *( af(i,j+1)-af(i,j)
558         &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
559         &        )
560             ENDDO
561            ENDDO
562    C-      keep advective flux (for diagnostics)
563            DO j=1-Oly,sNy+Oly
564             DO i=1-Olx,sNx+Olx
565              afy(i,j) = af(i,j)
566             ENDDO
567            ENDDO
568    
569  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
570  C--   Apply open boundary conditions  C-     Apply open boundary conditions
571        IF (useOBCS) THEN          IF (useOBCS) THEN
572         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN           IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
573          CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )            CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
574         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN           ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
575          CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )            CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
576         END IF           ENDIF
577        END IF          ENDIF
578  #endif /* ALLOW_OBCS */  #endif /* ALLOW_OBCS */
579    
580    C      end if/else update overlap-Only
581           ENDIF
582    
583  C--   End of Y direction  C--   End of Y direction
584        ENDIF        ENDIF
585    
# Line 483  C-    explicit advection is done ; store Line 591  C-    explicit advection is done ; store
591          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
592           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
593            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
594       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)
595           ENDDO           ENDDO
596          ENDDO          ENDDO
597        ELSE        ELSE
# Line 495  C-    horizontal advection done; store i Line 603  C-    horizontal advection done; store i
603         ENDDO         ENDDO
604        ENDIF        ENDIF
605    
606    #ifdef ALLOW_DIAGNOSTICS
607            IF ( useDiagnostics ) THEN
608              kk = -k
609              diagName = 'ADVx'//diagSufx
610              CALL DIAGNOSTICS_FILL(afx,diagName, kk,1, 2,bi,bj, myThid)
611              diagName = 'ADVy'//diagSufx
612              CALL DIAGNOSTICS_FILL(afy,diagName, kk,1, 2,bi,bj, myThid)
613            ENDIF
614    #endif
615    
616    #ifdef ALLOW_DEBUG
617          IF ( debugLevel .GE. debLevB
618         &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
619         &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
620         &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
621         &   .AND. useCubedSphereExchange ) THEN
622            CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION',
623         &             afx,afy, k, standardMessageUnit,bi,bj,myThid )
624          ENDIF
625    #endif /* ALLOW_DEBUG */
626    
627  C--   End of K loop for horizontal fluxes  C--   End of K loop for horizontal fluxes
628        ENDDO        ENDDO
629    
630    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
631    
632        IF ( .NOT.implicitAdvection ) THEN        IF ( .NOT.implicitAdvection ) THEN
633  C--   Start of k loop for vertical flux  C--   Start of k loop for vertical flux
634         DO k=Nr,1,-1         DO k=Nr,1,-1
# Line 528  C- Surface interface : Line 659  C- Surface interface :
659             rTransKp1(i,j) = kp1Msk*rTrans(i,j)             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
660             rTrans(i,j) = 0.             rTrans(i,j) = 0.
661             fVerT(i,j,kUp) = 0.             fVerT(i,j,kUp) = 0.
            af(i,j) = 0.  
662            ENDDO            ENDDO
663           ENDDO           ENDDO
664    
# Line 540  C- Interior interface : Line 670  C- Interior interface :
670             rTransKp1(i,j) = kp1Msk*rTrans(i,j)             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
671             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
672       &                 *maskC(i,j,k-1,bi,bj)       &                 *maskC(i,j,k-1,bi,bj)
673             af(i,j) = 0.             fVerT(i,j,kUp) = 0.
674            ENDDO            ENDDO
675           ENDDO           ENDDO
676    
# Line 560  CADJ &     = comlev1_bibj_k_gad, key=kke Line 690  CADJ &     = comlev1_bibj_k_gad, key=kke
690    
691  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
692           IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN           IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
693            CALL GAD_FLUXLIMIT_ADV_R(  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
694       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k),
695         I                               rTrans, wVel, localTijk,
696         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
697           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
698            CALL GAD_DST3_ADV_R(             CALL GAD_DST3_ADV_R(      bi,bj,k, dTtracerLev(k),
699       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       I                               rTrans, wVel, localTijk,
700         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
701           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
702            CALL GAD_DST3FL_ADV_R(             CALL GAD_DST3FL_ADV_R(    bi,bj,k, dTtracerLev(k),
703       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       I                               rTrans, wVel, localTijk,
704         O                               fVerT(1-Olx,1-Oly,kUp), myThid )
705           ELSE           ELSE
706            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'            STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
707           ENDIF           ENDIF
 C-    add the advective flux to fVerT  
          DO j=1-Oly,sNy+Oly  
           DO i=1-Olx,sNx+Olx  
            fVerT(i,j,kUp) = af(i,j)  
           ENDDO  
          ENDDO  
708    
709  C- end Surface/Interior if bloc  C- end Surface/Interior if bloc
710          ENDIF          ENDIF
# Line 591  CADJ &     = comlev1_bibj_k_gad, key=kke Line 719  CADJ &     = comlev1_bibj_k_gad, key=kke
719  C--   Divergence of vertical fluxes  C--   Divergence of vertical fluxes
720          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
721           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
722            localTij(i,j)=localTijk(i,j,k)-deltaTtracer*            localTij(i,j)=localTijk(i,j,k)-dTtracerLev(k)*
723       &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &     _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
724       &     *recip_rA(i,j,bi,bj)       &     *recip_rA(i,j,bi,bj)
725       &     *( fVerT(i,j,kUp)-fVerT(i,j,kDown)       &     *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
726       &       -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))       &       -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
727       &      )*rkFac       &      )*rkFac
728            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
729       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)
730           ENDDO           ENDDO
731          ENDDO          ENDDO
732    
733    #ifdef ALLOW_DIAGNOSTICS
734            IF ( useDiagnostics ) THEN
735              kk = -k
736              diagName = 'ADVr'//diagSufx
737              CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp),
738         &                           diagName, kk,1, 2,bi,bj, myThid)
739            ENDIF
740    #endif
741    
742  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
743         ENDDO         ENDDO
744  C--   end of if not.implicitAdvection block  C--   end of if not.implicitAdvection block
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22