/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_som_advect.F
ViewVC logotype

Diff of /MITgcm/pkg/generic_advdiff/gad_som_advect.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by jmc, Tue Jan 16 04:38:34 2007 UTC revision 1.2 by jmc, Tue Feb 12 20:32:34 2008 UTC
# Line 30  C !USES: =============================== Line 30  C !USES: ===============================
30  #include "PARAMS.h"  #include "PARAMS.h"
31  #include "GRID.h"  #include "GRID.h"
32  #include "GAD.h"  #include "GAD.h"
33    #ifdef ALLOW_EXCH2
34    #include "W2_EXCH2_TOPOLOGY.h"
35    #include "W2_EXCH2_PARAMS.h"
36    #endif /* ALLOW_EXCH2 */
37    
38  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
39  C  implicitAdvection :: implicit vertical advection (later on)  C  implicitAdvection :: implicit vertical advection (later on)
# Line 64  C  gTracer           :: tendency array Line 68  C  gTracer           :: tendency array
68    
69  #ifdef GAD_ALLOW_SOM_ADVECT  #ifdef GAD_ALLOW_SOM_ADVECT
70  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
 C  maskLocW      :: 2-D array for mask at West points  
 C  maskLocS      :: 2-D array for mask at South points  
71  C  maskUp        :: 2-D array mask for W points  C  maskUp        :: 2-D array mask for W points
 C  iMin,iMax,    :: loop range for called routines  
 C  jMin,jMax     :: loop range for called routines  
 C [iMin,iMax]Upd :: loop range to update tracer field  
 C [jMin,jMax]Upd :: loop range to update tracer field  
72  C  i,j,k         :: loop indices  C  i,j,k         :: loop indices
73  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
74  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
 C  kp1           :: =k+1 for k<Nr, =Nr for k=Nr  
75  C  xA,yA         :: areas of X and Y face of tracer cells  C  xA,yA         :: areas of X and Y face of tracer cells
76  C  uFld,vFld     :: 2-D local copy of horizontal velocity, U,V components  C  uFld,vFld     :: 2-D local copy of horizontal velocity, U,V components
77  C  wFld          :: 2-D local copy of vertical velocity  C  wFld          :: 2-D local copy of vertical velocity
78  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
79  C  rTrans        :: 2-D arrays of volume transports at W points  C  rTrans        :: 2-D arrays of volume transports at W points
 C  rTransKp1     :: vertical volume transport at interface k+1  
 C  af            :: 2-D array for horizontal advective flux  
80  C  afx           :: 2-D array for horizontal advective flux, x direction  C  afx           :: 2-D array for horizontal advective flux, x direction
81  C  afy           :: 2-D array for horizontal advective flux, y direction  C  afy           :: 2-D array for horizontal advective flux, y direction
82  C  afr           :: 2-D array for vertical advective flux  C  afr           :: 2-D array for vertical advective flux
83  C  fVerT         :: 2 1/2D arrays for vertical advective flux  C  fVerT         :: 2 1/2D arrays for vertical advective flux
84  C  localTij      :: 2-D array, temporary local copy of tracer fld  C  localTij      :: 2-D array, temporary local copy of tracer fld
85  C  localTijk     :: 3-D array, temporary local copy of tracer fld  C  localTijk     :: 3-D array, temporary local copy of tracer fld
 C  kp1Msk        :: flag (0,1) for over-riding mask for W levels  
86  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
87  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
88  C  interiorOnly  :: only update the interior of myTile, but not the edges  C  interiorOnly  :: only update the interior of myTile, but not the edges
# Line 100  C  nCFace        :: owns a tile for cube Line 94  C  nCFace        :: owns a tile for cube
94  C                :: multi-dim advection.  C                :: multi-dim advection.
95  C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube  C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
96  C  msgBuf        :: Informational/error meesage buffer  C  msgBuf        :: Informational/error meesage buffer
 c     _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
 c     _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
97        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       INTEGER iMin,iMax,jMin,jMax  
 c     INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd  
98        INTEGER i,j,k,km1,kUp,kDown        INTEGER i,j,k,km1,kUp,kDown
99        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 114  c     INTEGER iMinUpd,iMaxUpd,jMinUpd,jM Line 104  c     INTEGER iMinUpd,iMaxUpd,jMinUpd,jM
104        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 c     _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
 c     _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
107        _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108        _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109        _RL afr     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL afr     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 147  c     _RL localTijk(1-OLx:sNx+OLx,1-OLy: Line 135  c     _RL localTijk(1-OLx:sNx+OLx,1-OLy:
135        _RL  fn_xz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL  fn_xz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
136        _RL  fp_yz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL  fp_yz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
137        _RL  fn_yz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL  fn_yz  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
138          _RL  smCorners(OLx,OLy,4,-1:nSOM)
139        _RL  localTr        _RL  localTr
 c     _RL kp1Msk  
140        LOGICAL calc_fluxes_X, calc_fluxes_Y        LOGICAL calc_fluxes_X, calc_fluxes_Y
141  c     LOGICAL interiorOnly, overlapOnly        LOGICAL interiorOnly, overlapOnly
142        INTEGER limiter        INTEGER limiter
143        INTEGER npass, ipass        INTEGER npass, ipass
144  c     INTEGER nCFace        INTEGER nCFace, n
145  c     LOGICAL N_edge, S_edge, E_edge, W_edge        LOGICAL N_edge, S_edge, E_edge, W_edge
146        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
147    #ifdef ALLOW_EXCH2
148          INTEGER myTile
149    #endif
150  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
151        CHARACTER*8 diagName        CHARACTER*8 diagName
152        CHARACTER*4 GAD_DIAG_SUFX, diagSufx        CHARACTER*4 diagSufx
153          LOGICAL     doDiagAdvX, doDiagAdvY, doDiagAdvR
154    C-    Functions:
155          CHARACTER*4 GAD_DIAG_SUFX
156        EXTERNAL    GAD_DIAG_SUFX        EXTERNAL    GAD_DIAG_SUFX
157          LOGICAL  DIAGNOSTICS_IS_ON
158          EXTERNAL DIAGNOSTICS_IS_ON
159  #endif  #endif
160  CEOP  CEOP
161    
162  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
163  C--   Set diagnostic suffix for the current tracer  C--   Set diagnostics flags and suffix for the current tracer
164          doDiagAdvX = .FALSE.
165          doDiagAdvY = .FALSE.
166          doDiagAdvR = .FALSE.
167        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
168          diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )          diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
169            diagName = 'ADVx'//diagSufx
170            doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )
171            diagName = 'ADVy'//diagSufx
172            doDiagAdvY = DIAGNOSTICS_IS_ON( diagName, myThid )
173            diagName = 'ADVr'//diagSufx
174            doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )
175        ENDIF        ENDIF
176  #endif  #endif
177    
# Line 175  C     These inital values do not alter t Line 180  C     These inital values do not alter t
180  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
181  C     point numbers. This prevents spurious hardware signals due to  C     point numbers. This prevents spurious hardware signals due to
182  C     uninitialised but inert locations.  C     uninitialised but inert locations.
183  c     DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
184  c      DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
185             afx(i,j) = 0.
186             afy(i,j) = 0.
187  C-    xA,yA,uFld,vFld,uTrans,vTrans are set over the full domain  C-    xA,yA,uFld,vFld,uTrans,vTrans are set over the full domain
188  C     in CALC_COMMON_FACTORS: no need for extra initialisation  C     in CALC_COMMON_FACTORS: no need for extra initialisation
189  c       xA(i,j)      = 0. _d 0  c       xA(i,j)      = 0. _d 0
# Line 185  c       uTrans(i,j)  = 0. _d 0 Line 192  c       uTrans(i,j)  = 0. _d 0
192  c       vTrans(i,j)  = 0. _d 0  c       vTrans(i,j)  = 0. _d 0
193  C-    rTrans is set over the full domain: no need for extra initialisation  C-    rTrans is set over the full domain: no need for extra initialisation
194  c       rTrans(i,j)  = 0. _d 0  c       rTrans(i,j)  = 0. _d 0
195  c       rTransKp1(i,j)= 0. _d 0         ENDDO
196  c      ENDDO        ENDDO
197  c     ENDDO        DO n=-1,nSOM
198           DO k=1,4
199            DO j=1,OLy
200             DO i=1,OLx
201              smCorners(i,j,k,n) = 0.
202             ENDDO
203            ENDDO
204           ENDDO
205          ENDDO
206    
 C--   Set tile-specific parameters for horizontal fluxes  
       IF (useCubedSphereExchange) THEN  
         WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ',  
      &     'not coded for CubedSphere (useCubedSphereExchange=T)'  
         CALL PRINT_ERROR( msgBuf, myThid )  
         STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'  
       ENDIF  
207        IF ( implicitAdvection ) THEN        IF ( implicitAdvection ) THEN
208          WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ',          WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ',
209       &     'not coded for implicit-vertical Advection'       &     'not coded for implicit-vertical Advection'
# Line 209  C--   Set tile-specific parameters for h Line 217  C--   Set tile-specific parameters for h
217          STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'          STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'
218        ENDIF        ENDIF
219    
220    C--   Set tile-specific parameters for horizontal fluxes
221          IF (useCubedSphereExchange) THEN
222           npass  = 3
223    #ifdef ALLOW_AUTODIFF_TAMC
224           IF ( npass.GT.maxcube ) STOP 'maxcube needs to be = 3'
225    #endif
226    #ifdef ALLOW_EXCH2
227           myTile = W2_myTileList(bi)
228           nCFace = exch2_myFace(myTile)
229           N_edge = exch2_isNedge(myTile).EQ.1
230           S_edge = exch2_isSedge(myTile).EQ.1
231           E_edge = exch2_isEedge(myTile).EQ.1
232           W_edge = exch2_isWedge(myTile).EQ.1
233    #else
234           nCFace = bi
235           N_edge = .TRUE.
236           S_edge = .TRUE.
237           E_edge = .TRUE.
238           W_edge = .TRUE.
239    #endif
240          ELSE
241           npass  = 2
242           nCFace = 0
243           N_edge = .FALSE.
244           S_edge = .FALSE.
245           E_edge = .FALSE.
246           W_edge = .FALSE.
247          ENDIF
248    
249        limiter = MOD(advectionScheme, 10)        limiter = MOD(advectionScheme, 10)
       npass= 2  
       iMin = 1-OLx  
       iMax = sNx+OLx  
       jMin = 1-OLy  
       jMax = sNy+OLy  
250    
251  C--   Start of k loop for horizontal fluxes  C--   Start of k loop for horizontal fluxes
252        DO k=1,Nr        DO k=1,Nr
253    
254  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
255        CALL CALC_COMMON_FACTORS (         CALL CALC_COMMON_FACTORS (
256       I         uVel, vVel,       I         uVel, vVel,
257       O         uFld, vFld, uTrans, vTrans, xA, yA,       O         uFld, vFld, uTrans, vTrans, xA, yA,
258       I         k,bi,bj, myThid )       I         k,bi,bj, myThid )
259    
260  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
261  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
262        IF (useGMRedi)         IF (useGMRedi)
263       &   CALL GMREDI_CALC_UVFLOW(       &   CALL GMREDI_CALC_UVFLOW(
264       U                  uFld, vFld, uTrans, vTrans,       U                  uFld, vFld, uTrans, vTrans,
265       I                  k, bi, bj, myThid )       I                  k, bi, bj, myThid )
266  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
267    
268  C--   grid-box volume and tracer content (zero order moment)  C--   grid-box volume and tracer content (zero order moment)
269        DO j=1-OLy,sNy+OLy         DO j=1-OLy,sNy+OLy
270         DO i=1-OLx,sNx+OLx          DO i=1-OLx,sNx+OLx
271           smVol(i,j,k) = rA(i,j,bi,bj)*deepFac2C(k)           smVol(i,j,k) = rA(i,j,bi,bj)*deepFac2C(k)
272       &                *drF(k)*hFacC(i,j,k,bi,bj)       &                *drF(k)*hFacC(i,j,k,bi,bj)
273       &                *rhoFacC(k)       &                *rhoFacC(k)
# Line 243  C--   grid-box volume and tracer content Line 275  C--   grid-box volume and tracer content
275  C-    fill empty grid-box:  C-    fill empty grid-box:
276           smVol(i,j,k) = smVol(i,j,k)           smVol(i,j,k) = smVol(i,j,k)
277       &                + (1. _d 0 - maskC(i,j,k,bi,bj))       &                + (1. _d 0 - maskC(i,j,k,bi,bj))
278            ENDDO
279         ENDDO         ENDDO
       ENDDO  
   
 c     IF (useCubedSphereExchange) THEN  
 c       CALL FILL_CS_CORNER_UV_RS(  
 c    &            .FALSE., maskLocW,maskLocS, bi,bj, myThid )  
 c     ENDIF  
280    
281  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
282  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.
283        DO ipass=1,npass         DO ipass=1,npass
284    
285            interiorOnly = .FALSE.
286            overlapOnly  = .FALSE.
287            IF (useCubedSphereExchange) THEN
288    C-    CubedSphere : pass 3 times, with partial update of local tracer field
289             IF (ipass.EQ.1) THEN
290              overlapOnly   = MOD(nCFace,3).EQ.0
291              interiorOnly  = MOD(nCFace,3).NE.0
292              calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
293              calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
294             ELSEIF (ipass.EQ.2) THEN
295              overlapOnly   = MOD(nCFace,3).EQ.2
296              interiorOnly  = MOD(nCFace,3).EQ.1
297              calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
298              calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
299             ELSE
300              interiorOnly  = .TRUE.
301              calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
302              calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
303             ENDIF
304            ELSE
305  C-    not CubedSphere  C-    not CubedSphere
306          calc_fluxes_X = MOD(ipass,2).EQ.1            calc_fluxes_X = MOD(ipass,2).EQ.1
307          calc_fluxes_Y = .NOT.calc_fluxes_X            calc_fluxes_Y = .NOT.calc_fluxes_X
308            ENDIF
309    
310  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
 C--   X direction  
       IF (calc_fluxes_X) THEN  
311    
312  C-     Advective flux in X  C--   X direction
313          DO j=1-OLy,sNy+OLy  C-     Do not compute fluxes if
314           DO i=1-OLx,sNx+OLx  C       a) needed in overlap only
315            afx(i,j) = 0.  C   and b) the overlap of myTile are not cube-face Edges
316           ENDDO          IF ( calc_fluxes_X .AND.
317          ENDDO       &      (.NOT.overlapOnly .OR. N_edge .OR. S_edge)
318         &     ) THEN
319          IF ( advectionScheme.EQ.ENUM_SOM_PRATHER  
320       &     .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN  C-     Internal exchange for calculations in X
321            CALL GAD_SOM_ADV_X(            IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN
322               CALL GAD_SOM_PREP_CS_CORNER(
323         U                     smVol, smTr0, smTr, smCorners,
324         I                     .TRUE., overlapOnly, interiorOnly,
325         I                     N_edge, S_edge, E_edge, W_edge,
326         I                     ipass, k, Nr, bi, bj, myThid )
327              ENDIF
328    
329    C-     Solve advection in X and update moments
330              IF ( advectionScheme.EQ.ENUM_SOM_PRATHER
331         &       .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN
332               CALL GAD_SOM_ADV_X(
333       I                     bi,bj,k, limiter,       I                     bi,bj,k, limiter,
334         I                     overlapOnly, interiorOnly,
335         I                     N_edge, S_edge, E_edge, W_edge,
336       I                     dTtracerLev(k), uTrans,       I                     dTtracerLev(k), uTrans,
337       U                     smVol(1-OLx,1-OLy,k),       U                     smVol(1-OLx,1-OLy,k),
338       U                     smTr0(1-OLx,1-OLy,k),       U                     smTr0(1-OLx,1-OLy,k),
# Line 287  C-     Advective flux in X Line 346  C-     Advective flux in X
346       U                     smTr(1-OLx,1-OLy,k,bi,bj,8),       U                     smTr(1-OLx,1-OLy,k,bi,bj,8),
347       U                     smTr(1-OLx,1-OLy,k,bi,bj,9),       U                     smTr(1-OLx,1-OLy,k,bi,bj,9),
348       O                     afx, myThid )       O                     afx, myThid )
349          ELSE            ELSE
350           STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'             STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
351          ENDIF            ENDIF
   
 C-     Update the local tracer field where needed:  
   
 C-      keep advective flux (for diagnostics)  
 c       DO j=1-OLy,sNy+OLy  
 c        DO i=1-OLx,sNx+OLx  
 c         afx(i,j) = af(i,j)  
 c        ENDDO  
 c       ENDDO  
352    
353  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
354  C-     Apply open boundary conditions  C-     Apply open boundary conditions
355  c       IF ( useOBCS ) THEN  c         IF ( useOBCS ) THEN
356  ccc        localTij(i,j) = smTr0(i,j)/smVol(i,j)  ccc        localTij(i,j) = smTr0(i,j)/smVol(i,j)
357  c        IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN  c          IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
358  c         CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )  c           CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
359  c        ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN  c          ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
360  c         CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )  c           CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
361  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
362  c        ELSEIF (tracerIdentity.GE.GAD_TR1) THEN  c          ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
363  c         CALL OBCS_APPLY_PTRACER( bi, bj, k,  c           CALL OBCS_APPLY_PTRACER( bi, bj, k,
364  c    &         tracerIdentity-GAD_TR1+1, localTij, myThid )  c    &             tracerIdentity-GAD_TR1+1, localTij, myThid )
365  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
366  c        ENDIF  c          ENDIF
367  ccc        smTr0(i,j) = localTij(i,j)*smVol(i,j)  ccc        smTr0(i,j) = localTij(i,j)*smVol(i,j)
368  c       ENDIF  c         ENDIF
369  #endif /* ALLOW_OBCS */  #endif /* ALLOW_OBCS */
370    
371  C--   End of X direction  C--   End of X direction
372        ENDIF          ENDIF
373    
374  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
 C--   Y direction  
 C  
       IF (calc_fluxes_Y) THEN  
375    
376    C--   Y direction
377  C-     Do not compute fluxes if  C-     Do not compute fluxes if
378  C       a) needed in overlap only  C       a) needed in overlap only
379  C   and b) the overlap of myTile are not cube-face edges  C   and b) the overlap of myTile are not cube-face edges
380  c      IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN          IF ( calc_fluxes_Y .AND.
381         &      (.NOT.overlapOnly .OR. E_edge .OR. W_edge)
382  C-     Advective flux in Y       &     ) THEN
383          DO j=1-OLy,sNy+OLy  
384           DO i=1-OLx,sNx+OLx  C-     Internal exchange for calculations in Y
385            afy(i,j) = 0.            IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN
386           ENDDO             CALL GAD_SOM_PREP_CS_CORNER(
387          ENDDO       U                     smVol, smTr0, smTr, smCorners,
388         I                     .FALSE., overlapOnly, interiorOnly,
389          IF ( advectionScheme.EQ.ENUM_SOM_PRATHER       I                     N_edge, S_edge, E_edge, W_edge,
390       &     .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN       I                     iPass, k, Nr, bi, bj, myThid )
391            CALL GAD_SOM_ADV_Y(            ENDIF
392    
393    C-     Solve advection in Y and update moments
394              IF ( advectionScheme.EQ.ENUM_SOM_PRATHER
395         &       .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN
396               CALL GAD_SOM_ADV_Y(
397       I                     bi,bj,k, limiter,       I                     bi,bj,k, limiter,
398         I                     overlapOnly, interiorOnly,
399         I                     N_edge, S_edge, E_edge, W_edge,
400       I                     dTtracerLev(k), vTrans,       I                     dTtracerLev(k), vTrans,
401       U                     smVol(1-OLx,1-OLy,k),       U                     smVol(1-OLx,1-OLy,k),
402       U                     smTr0(1-OLx,1-OLy,k),       U                     smTr0(1-OLx,1-OLy,k),
# Line 355  C-     Advective flux in Y Line 410  C-     Advective flux in Y
410       U                     smTr(1-OLx,1-OLy,k,bi,bj,8),       U                     smTr(1-OLx,1-OLy,k,bi,bj,8),
411       U                     smTr(1-OLx,1-OLy,k,bi,bj,9),       U                     smTr(1-OLx,1-OLy,k,bi,bj,9),
412       O                     afy, myThid )       O                     afy, myThid )
413          ELSE            ELSE
414           STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'             STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
415          ENDIF            ENDIF
   
 C-     Advective flux in Y : done  
 c      ENDIF  
   
 C-     Update the local tracer field where needed:  
   
 C-      keep advective flux (for diagnostics)  
 c       DO j=1-OLy,sNy+OLy  
 c        DO i=1-OLx,sNx+OLx  
 c         afy(i,j) = af(i,j)  
 c        ENDDO  
 c       ENDDO  
416    
417  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
418  C-     Apply open boundary conditions  C-     Apply open boundary conditions
419  c       IF (useOBCS) THEN  c         IF (useOBCS) THEN
420  ccc        localTij(i,j) = smTr0(i,j)/smVol(i,j)  ccc        localTij(i,j) = smTr0(i,j)/smVol(i,j)
421  c        IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN  c          IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
422  c         CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )  c           CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
423  c        ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN  c          ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
424  c         CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )  c           CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
425  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
426  c        ELSEIF (tracerIdentity.GE.GAD_TR1) THEN  c          ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
427  c         CALL OBCS_APPLY_PTRACER( bi, bj, k,  c           CALL OBCS_APPLY_PTRACER( bi, bj, k,
428  c    &         tracerIdentity-GAD_TR1+1, localTij, myThid )  c    &             tracerIdentity-GAD_TR1+1, localTij, myThid )
429  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
430  c        ENDIF  c          ENDIF
431  ccc        smTr0(i,j) = localTij(i,j)*smVol(i,j)  ccc        smTr0(i,j) = localTij(i,j)*smVol(i,j)
432  c       ENDIF  c         ENDIF
433  #endif /* ALLOW_OBCS */  #endif /* ALLOW_OBCS */
434    
 C      end if/else update overlap-Only  
 c      ENDIF  
   
435  C--   End of Y direction  C--   End of Y direction
436        ENDIF          ENDIF
437    
438  C--   End of ipass loop  C--   End of ipass loop
439        ENDDO         ENDDO
440    
441        IF ( implicitAdvection ) THEN         IF ( implicitAdvection ) THEN
442  C-    explicit advection is done ; store tendency in gTracer:  C-    explicit advection is done ; store tendency in gTracer:
443          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
444           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
# Line 407  C-    explicit advection is done ; store Line 447  C-    explicit advection is done ; store
447       &                         / dTtracerLev(k)       &                         / dTtracerLev(k)
448           ENDDO           ENDDO
449          ENDDO          ENDDO
450        ELSE         ELSE
451  C-    horizontal advection done; store intermediate result in 3D array:  C-    horizontal advection done; store intermediate result in 3D array:
452  c       DO j=1-OLy,sNy+OLy  c       DO j=1-OLy,sNy+OLy
453  c        DO i=1-OLx,sNx+OLx  c        DO i=1-OLx,sNx+OLx
454  c         localTijk(i,j,k)=localTij(i,j)  c         localTijk(i,j,k)=localTij(i,j)
455  c        ENDDO  c        ENDDO
456  c       ENDDO  c       ENDDO
457        ENDIF         ENDIF
458    
459  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
460          IF ( useDiagnostics ) THEN         IF ( doDiagAdvX ) THEN
461            diagName = 'ADVx'//diagSufx           diagName = 'ADVx'//diagSufx
462            CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)           CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid )
463            diagName = 'ADVy'//diagSufx         ENDIF
464            CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)         IF ( doDiagAdvY ) THEN
465          ENDIF           diagName = 'ADVy'//diagSufx
466             CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid )
467           ENDIF
468  #endif  #endif
469    
470  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
471        IF ( debugLevel .GE. debLevB         IF ( debugLevel .GE. debLevB
472       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
473       &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0       &   .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
474       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
475       &   .AND. useCubedSphereExchange ) THEN       &   .AND. useCubedSphereExchange ) THEN
476          CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_SOM_ADVECT',          CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_SOM_ADVECT',
477       &             afx,afy, k, standardMessageUnit,bi,bj,myThid )       &             afx,afy, k, standardMessageUnit,bi,bj,myThid )
478        ENDIF         ENDIF
479  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
480    
481  C--   End of K loop for horizontal fluxes  C--   End of K loop for horizontal fluxes
# Line 463  C--   kUp    Cycles through 1,2 to point Line 505  C--   kUp    Cycles through 1,2 to point
505  C--   kDown  Cycles through 2,1 to point to w-layer below  C--   kDown  Cycles through 2,1 to point to w-layer below
506          kUp  = 1+MOD(Nr-k,2)          kUp  = 1+MOD(Nr-k,2)
507          kDown= 1+MOD(Nr-k+1,2)          kDown= 1+MOD(Nr-k+1,2)
 c       kp1=min(Nr,k+1)  
 c       kp1Msk = 1. _d 0  
508          IF (k.EQ.Nr) THEN          IF (k.EQ.Nr) THEN
 c        kp1Msk = 0. _d 0  
509  C--   Set advective fluxes at the very bottom:  C--   Set advective fluxes at the very bottom:
510           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
511            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
# Line 546  C- end Surface/Interior if bloc Line 585  C- end Surface/Interior if bloc
585    
586  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
587  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
588           IF (useGMRedi .AND. k.GT.1 )          IF (useGMRedi .AND. k.GT.1 )
589       &     CALL GMREDI_CALC_WFLOW(       &     CALL GMREDI_CALC_WFLOW(
590       U                 wFld, rTrans,       U                 wFld, rTrans,
591       I                 k, bi, bj, myThid )       I                 k, bi, bj, myThid )
# Line 554  C--   Residual transp = Bolus transp + E Line 593  C--   Residual transp = Bolus transp + E
593    
594  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
595          IF ( vertAdvecScheme.EQ.ENUM_SOM_PRATHER          IF ( vertAdvecScheme.EQ.ENUM_SOM_PRATHER
596       &      .OR. vertAdvecScheme.EQ.ENUM_SOM_LIMITER ) THEN       &     .OR. vertAdvecScheme.EQ.ENUM_SOM_LIMITER ) THEN
597            CALL GAD_SOM_ADV_R(             CALL GAD_SOM_ADV_R(
598       I                     bi,bj,k, kUp, kDown,       I                     bi,bj,k, kUp, kDown,
599       I                     dTtracerLev(k), rTrans, maskUp,       I                     dTtracerLev(k), rTrans, maskUp,
600       U                     smVol,       U                     smVol,
# Line 575  C-    Compute vertical advective flux in Line 614  C-    Compute vertical advective flux in
614       U                     fp_xy, fn_xy, fp_xz, fn_xz, fp_yz, fn_yz,       U                     fp_xy, fn_xy, fp_xz, fn_xz, fp_yz, fn_yz,
615       O                     afr, myThid )       O                     afr, myThid )
616          ELSE          ELSE
617            STOP 'GAD_SOM_ADVECT: adv.scheme incompatibale with SOM'             STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
618          ENDIF          ENDIF
619    
620  C--   Divergence of vertical fluxes  C--   Compute new tracer value and store tracer tendency
621          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
622           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
 c         localTij(i,j) = localTijk(i,j,k)  
 c    &      -dTtracerLev(k)*recip_rhoFacC(k)  
 c    &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  
 c    &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)  
 c    &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)  
 c    &         -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))  
 c    &        )*rkSign  
 c         gTracer(i,j,k,bi,bj)=  
 c    &     (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)  
623            localTr = smTr0(i,j,k)            localTr = smTr0(i,j,k)
624       &            *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &            *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
625       &            *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)       &            *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
# Line 601  c         localTr = smTr0(i,j,k)/smVol(i Line 631  c         localTr = smTr0(i,j,k)/smVol(i
631          ENDDO          ENDDO
632    
633  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
634          IF ( useDiagnostics ) THEN          IF ( doDiagAdvR ) THEN
635            diagName = 'ADVr'//diagSufx           diagName = 'ADVr'//diagSufx
636            CALL DIAGNOSTICS_FILL( afr,           CALL DIAGNOSTICS_FILL( afr,
637       &                           diagName, k,1, 2,bi,bj, myThid)       &                          diagName, k,1, 2,bi,bj, myThid )
638          ENDIF          ENDIF
639  #endif  #endif
640    
641  C--   End of K loop for vertical flux  C--   End of k loop for vertical flux
642         ENDDO         ENDDO
643  C--   end of if not.implicitAdvection block  C--   end of if not.implicitAdvection block
644        ENDIF        ENDIF
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22