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

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

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


Revision 1.2 - (hide annotations) (download)
Tue Feb 12 20:32:34 2008 UTC (16 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59q, checkpoint59p, checkpoint59o
Changes since 1.1: +196 -166 lines
prather advection scheme (SOM) coded for Cubed-Sphere grid

1 jmc 1.2 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_som_advect.F,v 1.1 2007/01/16 04:38:34 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "GAD_OPTIONS.h"
5    
6     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7     CBOP
8     C !ROUTINE: GAD_SOM_ADVECT
9    
10     C !INTERFACE: ==========================================================
11     SUBROUTINE GAD_SOM_ADVECT(
12     I implicitAdvection, advectionScheme, vertAdvecScheme,
13     I tracerIdentity,
14     I uVel, vVel, wVel, tracer,
15     U smTr,
16     O gTracer,
17     I bi,bj, myTime,myIter,myThid)
18    
19     C !DESCRIPTION:
20     C Calculates the tendency of a tracer due to advection.
21     C It uses the 2nd-Order moment advection scheme with multi-dimensional method
22     C see Prather, 1986, JGR, v.91, D-6, pp.6671-6681.
23     C
24     C The tendency (output) is over-written by this routine.
25    
26     C !USES: ===============================================================
27     IMPLICIT NONE
28     #include "SIZE.h"
29     #include "EEPARAMS.h"
30     #include "PARAMS.h"
31     #include "GRID.h"
32     #include "GAD.h"
33 jmc 1.2 #ifdef ALLOW_EXCH2
34     #include "W2_EXCH2_TOPOLOGY.h"
35     #include "W2_EXCH2_PARAMS.h"
36     #endif /* ALLOW_EXCH2 */
37 jmc 1.1
38     C !INPUT PARAMETERS: ===================================================
39     C implicitAdvection :: implicit vertical advection (later on)
40     C advectionScheme :: advection scheme to use (Horizontal plane)
41     C vertAdvecScheme :: advection scheme to use (vertical direction)
42     C tracerIdentity :: tracer identifier (required only for OBCS)
43     C uVel :: velocity, zonal component
44     C vVel :: velocity, meridional component
45     C wVel :: velocity, vertical component
46     C tracer :: tracer field
47     C bi,bj :: tile indices
48     C myTime :: current time
49     C myIter :: iteration number
50     C myThid :: thread number
51     LOGICAL implicitAdvection
52     INTEGER advectionScheme, vertAdvecScheme
53     INTEGER tracerIdentity
54     _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
55     _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
56     _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
57     _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
58     INTEGER bi,bj
59     _RL myTime
60     INTEGER myIter
61     INTEGER myThid
62    
63     C !OUTPUT PARAMETERS: ==================================================
64     C smTr :: tracer 1rst & 2nd Order moments
65     C gTracer :: tendency array
66     _RL smTr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy,nSOM)
67     _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
68    
69     #ifdef GAD_ALLOW_SOM_ADVECT
70     C !LOCAL VARIABLES: ====================================================
71     C maskUp :: 2-D array mask for W points
72     C i,j,k :: loop indices
73     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
75     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
77     C wFld :: 2-D local copy of vertical velocity
78     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
80     C afx :: 2-D array for horizontal advective flux, x direction
81     C afy :: 2-D array for horizontal advective flux, y direction
82     C afr :: 2-D array for vertical advective flux
83     C fVerT :: 2 1/2D arrays for vertical advective flux
84     C localTij :: 2-D array, temporary local copy of tracer fld
85     C localTijk :: 3-D array, temporary local copy of tracer fld
86     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
88     C interiorOnly :: only update the interior of myTile, but not the edges
89     C overlapOnly :: only update the edges of myTile, but not the interior
90     C npass :: number of passes in multi-dimensional method
91     C ipass :: number of the current pass being made
92     C myTile :: variables used to determine which cube face
93     C nCFace :: owns a tile for cube grid runs using
94     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
96     C msgBuf :: Informational/error meesage buffer
97     _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98     INTEGER i,j,k,km1,kUp,kDown
99     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101     _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102     _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103     _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104     _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105     _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106     _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107     _RL afx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108     _RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109     _RL afr (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110     ccc _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111     c _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
112     _RL smVol (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
113     _RL smTr0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
114     _RL alp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
115     _RL aln (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
116     _RL fp_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
117     _RL fn_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
118     _RL fp_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
119     _RL fn_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
120     _RL fp_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
121     _RL fn_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
122     _RL fp_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
123     _RL fn_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
124     _RL fp_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
125     _RL fn_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
126     _RL fp_xx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
127     _RL fn_xx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
128     _RL fp_yy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
129     _RL fn_yy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
130     _RL fp_zz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
131     _RL fn_zz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
132     _RL fp_xy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
133     _RL fn_xy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
134     _RL fp_xz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
135     _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)
137     _RL fn_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
138 jmc 1.2 _RL smCorners(OLx,OLy,4,-1:nSOM)
139 jmc 1.1 _RL localTr
140     LOGICAL calc_fluxes_X, calc_fluxes_Y
141 jmc 1.2 LOGICAL interiorOnly, overlapOnly
142 jmc 1.1 INTEGER limiter
143     INTEGER npass, ipass
144 jmc 1.2 INTEGER nCFace, n
145     LOGICAL N_edge, S_edge, E_edge, W_edge
146 jmc 1.1 CHARACTER*(MAX_LEN_MBUF) msgBuf
147 jmc 1.2 #ifdef ALLOW_EXCH2
148     INTEGER myTile
149     #endif
150 jmc 1.1 #ifdef ALLOW_DIAGNOSTICS
151     CHARACTER*8 diagName
152 jmc 1.2 CHARACTER*4 diagSufx
153     LOGICAL doDiagAdvX, doDiagAdvY, doDiagAdvR
154     C- Functions:
155     CHARACTER*4 GAD_DIAG_SUFX
156 jmc 1.1 EXTERNAL GAD_DIAG_SUFX
157 jmc 1.2 LOGICAL DIAGNOSTICS_IS_ON
158     EXTERNAL DIAGNOSTICS_IS_ON
159 jmc 1.1 #endif
160     CEOP
161    
162     #ifdef ALLOW_DIAGNOSTICS
163 jmc 1.2 C-- Set diagnostics flags and suffix for the current tracer
164     doDiagAdvX = .FALSE.
165     doDiagAdvY = .FALSE.
166     doDiagAdvR = .FALSE.
167 jmc 1.1 IF ( useDiagnostics ) THEN
168     diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
169 jmc 1.2 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 jmc 1.1 ENDIF
176     #endif
177    
178     C-- Set up work arrays with valid (i.e. not NaN) values
179     C These inital values do not alter the numerical results. They
180     C just ensure that all memory references are to valid floating
181     C point numbers. This prevents spurious hardware signals due to
182     C uninitialised but inert locations.
183 jmc 1.2 DO j=1-OLy,sNy+OLy
184     DO i=1-OLx,sNx+OLx
185     afx(i,j) = 0.
186     afy(i,j) = 0.
187 jmc 1.1 C- xA,yA,uFld,vFld,uTrans,vTrans are set over the full domain
188     C in CALC_COMMON_FACTORS: no need for extra initialisation
189     c xA(i,j) = 0. _d 0
190     c yA(i,j) = 0. _d 0
191     c uTrans(i,j) = 0. _d 0
192     c vTrans(i,j) = 0. _d 0
193     C- rTrans is set over the full domain: no need for extra initialisation
194     c rTrans(i,j) = 0. _d 0
195 jmc 1.2 ENDDO
196     ENDDO
197     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 jmc 1.1
207     IF ( implicitAdvection ) THEN
208     WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ',
209     & 'not coded for implicit-vertical Advection'
210     CALL PRINT_ERROR( msgBuf, myThid )
211     STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'
212     ENDIF
213     IF ( vertAdvecScheme .NE. advectionScheme ) THEN
214     WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ',
215     & 'not coded for different vertAdvecScheme'
216     CALL PRINT_ERROR( msgBuf, myThid )
217     STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'
218     ENDIF
219    
220 jmc 1.2 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 jmc 1.1 limiter = MOD(advectionScheme, 10)
250    
251     C-- Start of k loop for horizontal fluxes
252     DO k=1,Nr
253    
254     C-- Get temporary terms used by tendency routines
255 jmc 1.2 CALL CALC_COMMON_FACTORS (
256 jmc 1.1 I uVel, vVel,
257     O uFld, vFld, uTrans, vTrans, xA, yA,
258     I k,bi,bj, myThid )
259    
260     #ifdef ALLOW_GMREDI
261     C-- Residual transp = Bolus transp + Eulerian transp
262 jmc 1.2 IF (useGMRedi)
263 jmc 1.1 & CALL GMREDI_CALC_UVFLOW(
264     U uFld, vFld, uTrans, vTrans,
265     I k, bi, bj, myThid )
266     #endif /* ALLOW_GMREDI */
267    
268     C-- grid-box volume and tracer content (zero order moment)
269 jmc 1.2 DO j=1-OLy,sNy+OLy
270     DO i=1-OLx,sNx+OLx
271 jmc 1.1 smVol(i,j,k) = rA(i,j,bi,bj)*deepFac2C(k)
272     & *drF(k)*hFacC(i,j,k,bi,bj)
273     & *rhoFacC(k)
274     smTr0(i,j,k) = tracer(i,j,k,bi,bj)*smVol(i,j,k)
275     C- fill empty grid-box:
276     smVol(i,j,k) = smVol(i,j,k)
277     & + (1. _d 0 - maskC(i,j,k,bi,bj))
278 jmc 1.2 ENDDO
279 jmc 1.1 ENDDO
280    
281     C-- Multiple passes for different directions on different tiles
282     C-- For cube need one pass for each of red, green and blue axes.
283 jmc 1.2 DO ipass=1,npass
284 jmc 1.1
285 jmc 1.2 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 jmc 1.1 C- not CubedSphere
306 jmc 1.2 calc_fluxes_X = MOD(ipass,2).EQ.1
307     calc_fluxes_Y = .NOT.calc_fluxes_X
308     ENDIF
309 jmc 1.1
310     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
311 jmc 1.2
312 jmc 1.1 C-- X direction
313 jmc 1.2 C- Do not compute fluxes if
314     C a) needed in overlap only
315     C and b) the overlap of myTile are not cube-face Edges
316     IF ( calc_fluxes_X .AND.
317     & (.NOT.overlapOnly .OR. N_edge .OR. S_edge)
318     & ) THEN
319    
320     C- Internal exchange for calculations in X
321     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 jmc 1.1 I bi,bj,k, limiter,
334 jmc 1.2 I overlapOnly, interiorOnly,
335     I N_edge, S_edge, E_edge, W_edge,
336 jmc 1.1 I dTtracerLev(k), uTrans,
337     U smVol(1-OLx,1-OLy,k),
338     U smTr0(1-OLx,1-OLy,k),
339     U smTr(1-OLx,1-OLy,k,bi,bj,1),
340     U smTr(1-OLx,1-OLy,k,bi,bj,2),
341     U smTr(1-OLx,1-OLy,k,bi,bj,3),
342     U smTr(1-OLx,1-OLy,k,bi,bj,4),
343     U smTr(1-OLx,1-OLy,k,bi,bj,5),
344     U smTr(1-OLx,1-OLy,k,bi,bj,6),
345     U smTr(1-OLx,1-OLy,k,bi,bj,7),
346     U smTr(1-OLx,1-OLy,k,bi,bj,8),
347     U smTr(1-OLx,1-OLy,k,bi,bj,9),
348     O afx, myThid )
349 jmc 1.2 ELSE
350     STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
351     ENDIF
352 jmc 1.1
353     #ifdef ALLOW_OBCS
354     C- Apply open boundary conditions
355 jmc 1.2 c IF ( useOBCS ) THEN
356 jmc 1.1 ccc localTij(i,j) = smTr0(i,j)/smVol(i,j)
357 jmc 1.2 c IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
358     c CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
359     c ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
360     c CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
361 jmc 1.1 #ifdef ALLOW_PTRACERS
362 jmc 1.2 c ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
363     c CALL OBCS_APPLY_PTRACER( bi, bj, k,
364     c & tracerIdentity-GAD_TR1+1, localTij, myThid )
365 jmc 1.1 #endif /* ALLOW_PTRACERS */
366 jmc 1.2 c ENDIF
367 jmc 1.1 ccc smTr0(i,j) = localTij(i,j)*smVol(i,j)
368 jmc 1.2 c ENDIF
369 jmc 1.1 #endif /* ALLOW_OBCS */
370    
371     C-- End of X direction
372 jmc 1.2 ENDIF
373 jmc 1.1
374     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
375 jmc 1.2
376 jmc 1.1 C-- Y direction
377     C- Do not compute fluxes if
378     C a) needed in overlap only
379     C and b) the overlap of myTile are not cube-face edges
380 jmc 1.2 IF ( calc_fluxes_Y .AND.
381     & (.NOT.overlapOnly .OR. E_edge .OR. W_edge)
382     & ) THEN
383    
384     C- Internal exchange for calculations in Y
385     IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN
386     CALL GAD_SOM_PREP_CS_CORNER(
387     U smVol, smTr0, smTr, smCorners,
388     I .FALSE., overlapOnly, interiorOnly,
389     I N_edge, S_edge, E_edge, W_edge,
390     I iPass, k, Nr, bi, bj, myThid )
391     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 jmc 1.1 I bi,bj,k, limiter,
398 jmc 1.2 I overlapOnly, interiorOnly,
399     I N_edge, S_edge, E_edge, W_edge,
400 jmc 1.1 I dTtracerLev(k), vTrans,
401     U smVol(1-OLx,1-OLy,k),
402     U smTr0(1-OLx,1-OLy,k),
403     U smTr(1-OLx,1-OLy,k,bi,bj,1),
404     U smTr(1-OLx,1-OLy,k,bi,bj,2),
405     U smTr(1-OLx,1-OLy,k,bi,bj,3),
406     U smTr(1-OLx,1-OLy,k,bi,bj,4),
407     U smTr(1-OLx,1-OLy,k,bi,bj,5),
408     U smTr(1-OLx,1-OLy,k,bi,bj,6),
409     U smTr(1-OLx,1-OLy,k,bi,bj,7),
410     U smTr(1-OLx,1-OLy,k,bi,bj,8),
411     U smTr(1-OLx,1-OLy,k,bi,bj,9),
412     O afy, myThid )
413 jmc 1.2 ELSE
414     STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
415     ENDIF
416 jmc 1.1
417     #ifdef ALLOW_OBCS
418     C- Apply open boundary conditions
419 jmc 1.2 c IF (useOBCS) THEN
420 jmc 1.1 ccc localTij(i,j) = smTr0(i,j)/smVol(i,j)
421 jmc 1.2 c IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
422     c CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
423     c ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
424     c CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
425 jmc 1.1 #ifdef ALLOW_PTRACERS
426 jmc 1.2 c ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
427     c CALL OBCS_APPLY_PTRACER( bi, bj, k,
428     c & tracerIdentity-GAD_TR1+1, localTij, myThid )
429 jmc 1.1 #endif /* ALLOW_PTRACERS */
430 jmc 1.2 c ENDIF
431 jmc 1.1 ccc smTr0(i,j) = localTij(i,j)*smVol(i,j)
432 jmc 1.2 c ENDIF
433 jmc 1.1 #endif /* ALLOW_OBCS */
434    
435     C-- End of Y direction
436 jmc 1.2 ENDIF
437 jmc 1.1
438     C-- End of ipass loop
439 jmc 1.2 ENDDO
440 jmc 1.1
441 jmc 1.2 IF ( implicitAdvection ) THEN
442 jmc 1.1 C- explicit advection is done ; store tendency in gTracer:
443     DO j=1-OLy,sNy+OLy
444     DO i=1-OLx,sNx+OLx
445     localTr = smTr0(i,j,k)/smVol(i,j,k)
446     gTracer(i,j,k,bi,bj) = ( localTr - tracer(i,j,k,bi,bj) )
447     & / dTtracerLev(k)
448     ENDDO
449     ENDDO
450 jmc 1.2 ELSE
451 jmc 1.1 C- horizontal advection done; store intermediate result in 3D array:
452     c DO j=1-OLy,sNy+OLy
453     c DO i=1-OLx,sNx+OLx
454     c localTijk(i,j,k)=localTij(i,j)
455     c ENDDO
456     c ENDDO
457 jmc 1.2 ENDIF
458 jmc 1.1
459     #ifdef ALLOW_DIAGNOSTICS
460 jmc 1.2 IF ( doDiagAdvX ) THEN
461     diagName = 'ADVx'//diagSufx
462     CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid )
463     ENDIF
464     IF ( doDiagAdvY ) THEN
465     diagName = 'ADVy'//diagSufx
466     CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid )
467     ENDIF
468 jmc 1.1 #endif
469    
470     #ifdef ALLOW_DEBUG
471 jmc 1.2 IF ( debugLevel .GE. debLevB
472 jmc 1.1 & .AND. tracerIdentity.EQ.GAD_TEMPERATURE
473     & .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
474     & .AND. nPx.EQ.1 .AND. nPy.EQ.1
475     & .AND. useCubedSphereExchange ) THEN
476     CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_SOM_ADVECT',
477     & afx,afy, k, standardMessageUnit,bi,bj,myThid )
478 jmc 1.2 ENDIF
479 jmc 1.1 #endif /* ALLOW_DEBUG */
480    
481     C-- End of K loop for horizontal fluxes
482     ENDDO
483    
484     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
485    
486     IF ( .NOT.implicitAdvection ) THEN
487     C-- Apply limiter (if any):
488     CALL GAD_SOM_LIM_R( bi,bj, limiter,
489     U smVol,
490     U smTr0,
491     U smTr(1-OLx,1-OLy,1,bi,bj,1),
492     U smTr(1-OLx,1-OLy,1,bi,bj,2),
493     U smTr(1-OLx,1-OLy,1,bi,bj,3),
494     U smTr(1-OLx,1-OLy,1,bi,bj,4),
495     U smTr(1-OLx,1-OLy,1,bi,bj,5),
496     U smTr(1-OLx,1-OLy,1,bi,bj,6),
497     U smTr(1-OLx,1-OLy,1,bi,bj,7),
498     U smTr(1-OLx,1-OLy,1,bi,bj,8),
499     U smTr(1-OLx,1-OLy,1,bi,bj,9),
500     I myThid )
501    
502     C-- Start of k loop for vertical flux
503     DO k=Nr,1,-1
504     C-- kUp Cycles through 1,2 to point to w-layer above
505     C-- kDown Cycles through 2,1 to point to w-layer below
506     kUp = 1+MOD(Nr-k,2)
507     kDown= 1+MOD(Nr-k+1,2)
508     IF (k.EQ.Nr) THEN
509     C-- Set advective fluxes at the very bottom:
510     DO j=1-OLy,sNy+OLy
511     DO i=1-OLx,sNx+OLx
512     alp (i,j,kDown) = 0. _d 0
513     aln (i,j,kDown) = 0. _d 0
514     fp_v (i,j,kDown) = 0. _d 0
515     fn_v (i,j,kDown) = 0. _d 0
516     fp_o (i,j,kDown) = 0. _d 0
517     fn_o (i,j,kDown) = 0. _d 0
518     fp_x (i,j,kDown) = 0. _d 0
519     fn_x (i,j,kDown) = 0. _d 0
520     fp_y (i,j,kDown) = 0. _d 0
521     fn_y (i,j,kDown) = 0. _d 0
522     fp_z (i,j,kDown) = 0. _d 0
523     fn_z (i,j,kDown) = 0. _d 0
524     fp_xx(i,j,kDown) = 0. _d 0
525     fn_xx(i,j,kDown) = 0. _d 0
526     fp_yy(i,j,kDown) = 0. _d 0
527     fn_yy(i,j,kDown) = 0. _d 0
528     fp_zz(i,j,kDown) = 0. _d 0
529     fn_zz(i,j,kDown) = 0. _d 0
530     fp_xy(i,j,kDown) = 0. _d 0
531     fn_xy(i,j,kDown) = 0. _d 0
532     fp_xz(i,j,kDown) = 0. _d 0
533     fn_xz(i,j,kDown) = 0. _d 0
534     fp_yz(i,j,kDown) = 0. _d 0
535     fn_yz(i,j,kDown) = 0. _d 0
536     ENDDO
537     ENDDO
538     ENDIF
539    
540     C-- Compute Vertical transport
541     #ifdef ALLOW_AIM
542     C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
543     c IF ( k.EQ.1 .OR.
544     c & (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
545     c & ) THEN
546     #else
547     c IF ( k.EQ.1 ) THEN
548     #endif
549     IF ( (rigidLid.OR.nonlinFreeSurf.GE.1) .AND. k.EQ.1 ) THEN
550     C- Surface interface :
551     DO j=1-OLy,sNy+OLy
552     DO i=1-OLx,sNx+OLx
553     wFld(i,j) = 0.
554     rTrans(i,j) = 0.
555     maskUp(i,j) = 0.
556     ENDDO
557     ENDDO
558    
559     ELSEIF ( rigidLid.OR.nonlinFreeSurf.GE.1 ) THEN
560     C- Interior interface :
561     DO j=1-OLy,sNy+OLy
562     DO i=1-OLx,sNx+OLx
563     wFld(i,j) = wVel(i,j,k,bi,bj)
564     rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
565     & *deepFac2F(k)*rhoFacF(k)
566     & *maskC(i,j,k-1,bi,bj)
567     maskUp(i,j) = 1.
568     ENDDO
569     ENDDO
570    
571     ELSE
572     C- Linear Free-Surface: do not mask rTrans :
573     km1= MAX(k-1,1)
574     DO j=1-OLy,sNy+OLy
575     DO i=1-OLx,sNx+OLx
576     wFld(i,j) = wVel(i,j,k,bi,bj)
577     rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
578     & *deepFac2F(k)*rhoFacF(k)
579     maskUp(i,j) = maskC(i,j,km1,bi,bj)*maskC(i,j,k,bi,bj)
580     ENDDO
581     ENDDO
582    
583     C- end Surface/Interior if bloc
584     ENDIF
585    
586     #ifdef ALLOW_GMREDI
587     C-- Residual transp = Bolus transp + Eulerian transp
588 jmc 1.2 IF (useGMRedi .AND. k.GT.1 )
589 jmc 1.1 & CALL GMREDI_CALC_WFLOW(
590     U wFld, rTrans,
591     I k, bi, bj, myThid )
592     #endif /* ALLOW_GMREDI */
593    
594     C- Compute vertical advective flux in the interior:
595     IF ( vertAdvecScheme.EQ.ENUM_SOM_PRATHER
596 jmc 1.2 & .OR. vertAdvecScheme.EQ.ENUM_SOM_LIMITER ) THEN
597     CALL GAD_SOM_ADV_R(
598 jmc 1.1 I bi,bj,k, kUp, kDown,
599     I dTtracerLev(k), rTrans, maskUp,
600     U smVol,
601     U smTr0,
602     U smTr(1-OLx,1-OLy,1,bi,bj,1),
603     U smTr(1-OLx,1-OLy,1,bi,bj,2),
604     U smTr(1-OLx,1-OLy,1,bi,bj,3),
605     U smTr(1-OLx,1-OLy,1,bi,bj,4),
606     U smTr(1-OLx,1-OLy,1,bi,bj,5),
607     U smTr(1-OLx,1-OLy,1,bi,bj,6),
608     U smTr(1-OLx,1-OLy,1,bi,bj,7),
609     U smTr(1-OLx,1-OLy,1,bi,bj,8),
610     U smTr(1-OLx,1-OLy,1,bi,bj,9),
611     U alp, aln, fp_v, fn_v, fp_o, fn_o,
612     U fp_x, fn_x, fp_y, fn_y, fp_z, fn_z,
613     U fp_xx, fn_xx, fp_yy, fn_yy, fp_zz, fn_zz,
614     U fp_xy, fn_xy, fp_xz, fn_xz, fp_yz, fn_yz,
615     O afr, myThid )
616     ELSE
617 jmc 1.2 STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
618 jmc 1.1 ENDIF
619    
620 jmc 1.2 C-- Compute new tracer value and store tracer tendency
621 jmc 1.1 DO j=1-OLy,sNy+OLy
622     DO i=1-OLx,sNx+OLx
623     localTr = smTr0(i,j,k)
624     & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
625     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
626     & *recip_rhoFacC(k)
627     c localTr = smTr0(i,j,k)/smVol(i,j,k)
628     gTracer(i,j,k,bi,bj) = ( localTr - tracer(i,j,k,bi,bj) )
629     & / dTtracerLev(k)
630     ENDDO
631     ENDDO
632    
633     #ifdef ALLOW_DIAGNOSTICS
634 jmc 1.2 IF ( doDiagAdvR ) THEN
635     diagName = 'ADVr'//diagSufx
636     CALL DIAGNOSTICS_FILL( afr,
637     & diagName, k,1, 2,bi,bj, myThid )
638 jmc 1.1 ENDIF
639     #endif
640    
641 jmc 1.2 C-- End of k loop for vertical flux
642 jmc 1.1 ENDDO
643     C-- end of if not.implicitAdvection block
644     ENDIF
645    
646     #endif /* GAD_ALLOW_SOM_ADVECT */
647    
648     RETURN
649     END

  ViewVC Help
Powered by ViewVC 1.1.22