/[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.4 - (hide annotations) (download)
Fri May 9 22:32:35 2008 UTC (16 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint60, checkpoint61, checkpoint59r, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a
Changes since 1.3: +1 -4 lines 
clean-up
1 jmc 1.4 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_som_advect.F,v 1.3 2008/05/09 21:43:16 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     C !LOCAL VARIABLES: ====================================================
70     C maskUp :: 2-D array mask for W points
71     C i,j,k :: loop indices
72     C kUp :: index into 2 1/2D array, toggles between 1 and 2
73     C kDown :: index into 2 1/2D array, toggles between 2 and 1
74     C xA,yA :: areas of X and Y face of tracer cells
75     C uFld,vFld :: 2-D local copy of horizontal velocity, U,V components
76     C wFld :: 2-D local copy of vertical velocity
77     C uTrans,vTrans :: 2-D arrays of volume transports at U,V points
78     C rTrans :: 2-D arrays of volume transports at W points
79     C afx :: 2-D array for horizontal advective flux, x direction
80     C afy :: 2-D array for horizontal advective flux, y direction
81     C afr :: 2-D array for vertical advective flux
82     C fVerT :: 2 1/2D arrays for vertical advective flux
83     C localTij :: 2-D array, temporary local copy of tracer fld
84     C localTijk :: 3-D array, temporary local copy of tracer fld
85     C calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
86     C calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
87     C interiorOnly :: only update the interior of myTile, but not the edges
88     C overlapOnly :: only update the edges of myTile, but not the interior
89     C npass :: number of passes in multi-dimensional method
90     C ipass :: number of the current pass being made
91     C myTile :: variables used to determine which cube face
92     C nCFace :: owns a tile for cube grid runs using
93     C :: multi-dim advection.
94     C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
95     C msgBuf :: Informational/error meesage buffer
96     _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97     INTEGER i,j,k,km1,kUp,kDown
98     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100     _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101     _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102     _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103     _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104     _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105     _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106     _RL afx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107     _RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108     _RL afr (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109     ccc _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110     c _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
111     _RL smVol (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
112     _RL smTr0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
113     _RL alp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
114     _RL aln (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
115     _RL fp_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
116     _RL fn_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
117     _RL fp_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
118     _RL fn_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
119     _RL fp_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
120     _RL fn_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
121     _RL fp_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
122     _RL fn_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
123     _RL fp_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
124     _RL fn_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
125     _RL fp_xx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
126     _RL fn_xx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
127     _RL fp_yy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
128     _RL fn_yy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
129     _RL fp_zz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
130     _RL fn_zz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
131     _RL fp_xy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
132     _RL fn_xy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
133     _RL fp_xz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
134     _RL fn_xz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
135     _RL fp_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
136     _RL fn_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
137 jmc 1.2 _RL smCorners(OLx,OLy,4,-1:nSOM)
138 jmc 1.1 _RL localTr
139     LOGICAL calc_fluxes_X, calc_fluxes_Y
140 jmc 1.2 LOGICAL interiorOnly, overlapOnly
141 jmc 1.1 INTEGER limiter
142     INTEGER npass, ipass
143 jmc 1.2 INTEGER nCFace, n
144     LOGICAL N_edge, S_edge, E_edge, W_edge
145 jmc 1.1 CHARACTER*(MAX_LEN_MBUF) msgBuf
146 jmc 1.2 #ifdef ALLOW_EXCH2
147     INTEGER myTile
148     #endif
149 jmc 1.1 #ifdef ALLOW_DIAGNOSTICS
150     CHARACTER*8 diagName
151 jmc 1.2 CHARACTER*4 diagSufx
152     LOGICAL doDiagAdvX, doDiagAdvY, doDiagAdvR
153     C- Functions:
154     CHARACTER*4 GAD_DIAG_SUFX
155 jmc 1.1 EXTERNAL GAD_DIAG_SUFX
156 jmc 1.2 LOGICAL DIAGNOSTICS_IS_ON
157     EXTERNAL DIAGNOSTICS_IS_ON
158 jmc 1.1 #endif
159     CEOP
160    
161     #ifdef ALLOW_DIAGNOSTICS
162 jmc 1.2 C-- Set diagnostics flags and suffix for the current tracer
163     doDiagAdvX = .FALSE.
164     doDiagAdvY = .FALSE.
165     doDiagAdvR = .FALSE.
166 jmc 1.1 IF ( useDiagnostics ) THEN
167     diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
168 jmc 1.2 diagName = 'ADVx'//diagSufx
169     doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )
170     diagName = 'ADVy'//diagSufx
171     doDiagAdvY = DIAGNOSTICS_IS_ON( diagName, myThid )
172     diagName = 'ADVr'//diagSufx
173     doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )
174 jmc 1.1 ENDIF
175     #endif
176    
177     C-- Set up work arrays with valid (i.e. not NaN) values
178     C These inital values do not alter the numerical results. They
179     C just ensure that all memory references are to valid floating
180     C point numbers. This prevents spurious hardware signals due to
181     C uninitialised but inert locations.
182 jmc 1.2 DO j=1-OLy,sNy+OLy
183     DO i=1-OLx,sNx+OLx
184     afx(i,j) = 0.
185     afy(i,j) = 0.
186 jmc 1.1 C- xA,yA,uFld,vFld,uTrans,vTrans are set over the full domain
187     C in CALC_COMMON_FACTORS: no need for extra initialisation
188     c xA(i,j) = 0. _d 0
189     c yA(i,j) = 0. _d 0
190     c uTrans(i,j) = 0. _d 0
191     c vTrans(i,j) = 0. _d 0
192     C- rTrans is set over the full domain: no need for extra initialisation
193     c rTrans(i,j) = 0. _d 0
194 jmc 1.2 ENDDO
195     ENDDO
196     DO n=-1,nSOM
197     DO k=1,4
198     DO j=1,OLy
199     DO i=1,OLx
200     smCorners(i,j,k,n) = 0.
201     ENDDO
202     ENDDO
203     ENDDO
204     ENDDO
205 jmc 1.1
206     IF ( implicitAdvection ) THEN
207     WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ',
208     & 'not coded for implicit-vertical Advection'
209     CALL PRINT_ERROR( msgBuf, myThid )
210     STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'
211     ENDIF
212     IF ( vertAdvecScheme .NE. advectionScheme ) THEN
213     WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ',
214     & 'not coded for different vertAdvecScheme'
215     CALL PRINT_ERROR( msgBuf, myThid )
216     STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'
217     ENDIF
218    
219 jmc 1.2 C-- Set tile-specific parameters for horizontal fluxes
220     IF (useCubedSphereExchange) THEN
221     npass = 3
222     #ifdef ALLOW_EXCH2
223     myTile = W2_myTileList(bi)
224     nCFace = exch2_myFace(myTile)
225     N_edge = exch2_isNedge(myTile).EQ.1
226     S_edge = exch2_isSedge(myTile).EQ.1
227     E_edge = exch2_isEedge(myTile).EQ.1
228     W_edge = exch2_isWedge(myTile).EQ.1
229     #else
230     nCFace = bi
231     N_edge = .TRUE.
232     S_edge = .TRUE.
233     E_edge = .TRUE.
234     W_edge = .TRUE.
235     #endif
236     ELSE
237     npass = 2
238     nCFace = 0
239     N_edge = .FALSE.
240     S_edge = .FALSE.
241     E_edge = .FALSE.
242     W_edge = .FALSE.
243     ENDIF
244    
245 jmc 1.1 limiter = MOD(advectionScheme, 10)
246    
247     C-- Start of k loop for horizontal fluxes
248     DO k=1,Nr
249    
250     C-- Get temporary terms used by tendency routines
251 jmc 1.2 CALL CALC_COMMON_FACTORS (
252 jmc 1.1 I uVel, vVel,
253     O uFld, vFld, uTrans, vTrans, xA, yA,
254     I k,bi,bj, myThid )
255    
256     #ifdef ALLOW_GMREDI
257     C-- Residual transp = Bolus transp + Eulerian transp
258 jmc 1.2 IF (useGMRedi)
259 jmc 1.1 & CALL GMREDI_CALC_UVFLOW(
260     U uFld, vFld, uTrans, vTrans,
261     I k, bi, bj, myThid )
262     #endif /* ALLOW_GMREDI */
263    
264     C-- grid-box volume and tracer content (zero order moment)
265 jmc 1.2 DO j=1-OLy,sNy+OLy
266     DO i=1-OLx,sNx+OLx
267 jmc 1.1 smVol(i,j,k) = rA(i,j,bi,bj)*deepFac2C(k)
268     & *drF(k)*hFacC(i,j,k,bi,bj)
269     & *rhoFacC(k)
270     smTr0(i,j,k) = tracer(i,j,k,bi,bj)*smVol(i,j,k)
271     C- fill empty grid-box:
272     smVol(i,j,k) = smVol(i,j,k)
273     & + (1. _d 0 - maskC(i,j,k,bi,bj))
274 jmc 1.2 ENDDO
275 jmc 1.1 ENDDO
276    
277     C-- Multiple passes for different directions on different tiles
278     C-- For cube need one pass for each of red, green and blue axes.
279 jmc 1.2 DO ipass=1,npass
280 jmc 1.1
281 jmc 1.2 interiorOnly = .FALSE.
282     overlapOnly = .FALSE.
283     IF (useCubedSphereExchange) THEN
284     C- CubedSphere : pass 3 times, with partial update of local tracer field
285     IF (ipass.EQ.1) THEN
286     overlapOnly = MOD(nCFace,3).EQ.0
287     interiorOnly = MOD(nCFace,3).NE.0
288     calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
289     calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
290     ELSEIF (ipass.EQ.2) THEN
291     overlapOnly = MOD(nCFace,3).EQ.2
292     interiorOnly = MOD(nCFace,3).EQ.1
293     calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
294     calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
295     ELSE
296     interiorOnly = .TRUE.
297     calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
298     calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
299     ENDIF
300     ELSE
301 jmc 1.1 C- not CubedSphere
302 jmc 1.2 calc_fluxes_X = MOD(ipass,2).EQ.1
303     calc_fluxes_Y = .NOT.calc_fluxes_X
304     ENDIF
305 jmc 1.1
306     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
307 jmc 1.2
308 jmc 1.1 C-- X direction
309 jmc 1.2 C- Do not compute fluxes if
310     C a) needed in overlap only
311     C and b) the overlap of myTile are not cube-face Edges
312     IF ( calc_fluxes_X .AND.
313     & (.NOT.overlapOnly .OR. N_edge .OR. S_edge)
314     & ) THEN
315    
316     C- Internal exchange for calculations in X
317     IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN
318     CALL GAD_SOM_PREP_CS_CORNER(
319     U smVol, smTr0, smTr, smCorners,
320     I .TRUE., overlapOnly, interiorOnly,
321     I N_edge, S_edge, E_edge, W_edge,
322     I ipass, k, Nr, bi, bj, myThid )
323     ENDIF
324    
325     C- Solve advection in X and update moments
326     IF ( advectionScheme.EQ.ENUM_SOM_PRATHER
327     & .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN
328     CALL GAD_SOM_ADV_X(
329 jmc 1.1 I bi,bj,k, limiter,
330 jmc 1.2 I overlapOnly, interiorOnly,
331     I N_edge, S_edge, E_edge, W_edge,
332 jmc 1.1 I dTtracerLev(k), uTrans,
333     U smVol(1-OLx,1-OLy,k),
334     U smTr0(1-OLx,1-OLy,k),
335     U smTr(1-OLx,1-OLy,k,bi,bj,1),
336     U smTr(1-OLx,1-OLy,k,bi,bj,2),
337     U smTr(1-OLx,1-OLy,k,bi,bj,3),
338     U smTr(1-OLx,1-OLy,k,bi,bj,4),
339     U smTr(1-OLx,1-OLy,k,bi,bj,5),
340     U smTr(1-OLx,1-OLy,k,bi,bj,6),
341     U smTr(1-OLx,1-OLy,k,bi,bj,7),
342     U smTr(1-OLx,1-OLy,k,bi,bj,8),
343     U smTr(1-OLx,1-OLy,k,bi,bj,9),
344     O afx, myThid )
345 jmc 1.2 ELSE
346     STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
347     ENDIF
348 jmc 1.1
349     #ifdef ALLOW_OBCS
350     C- Apply open boundary conditions
351 jmc 1.2 c IF ( useOBCS ) THEN
352 jmc 1.1 ccc localTij(i,j) = smTr0(i,j)/smVol(i,j)
353 jmc 1.2 c IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
354     c CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
355     c ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
356     c CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
357 jmc 1.1 #ifdef ALLOW_PTRACERS
358 jmc 1.2 c ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
359     c CALL OBCS_APPLY_PTRACER( bi, bj, k,
360     c & tracerIdentity-GAD_TR1+1, localTij, myThid )
361 jmc 1.1 #endif /* ALLOW_PTRACERS */
362 jmc 1.2 c ENDIF
363 jmc 1.1 ccc smTr0(i,j) = localTij(i,j)*smVol(i,j)
364 jmc 1.2 c ENDIF
365 jmc 1.1 #endif /* ALLOW_OBCS */
366    
367     C-- End of X direction
368 jmc 1.2 ENDIF
369 jmc 1.1
370     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
371 jmc 1.2
372 jmc 1.1 C-- Y direction
373     C- Do not compute fluxes if
374     C a) needed in overlap only
375     C and b) the overlap of myTile are not cube-face edges
376 jmc 1.2 IF ( calc_fluxes_Y .AND.
377     & (.NOT.overlapOnly .OR. E_edge .OR. W_edge)
378     & ) THEN
379    
380     C- Internal exchange for calculations in Y
381     IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN
382     CALL GAD_SOM_PREP_CS_CORNER(
383     U smVol, smTr0, smTr, smCorners,
384     I .FALSE., overlapOnly, interiorOnly,
385     I N_edge, S_edge, E_edge, W_edge,
386     I iPass, k, Nr, bi, bj, myThid )
387     ENDIF
388    
389     C- Solve advection in Y and update moments
390     IF ( advectionScheme.EQ.ENUM_SOM_PRATHER
391     & .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN
392     CALL GAD_SOM_ADV_Y(
393 jmc 1.1 I bi,bj,k, limiter,
394 jmc 1.2 I overlapOnly, interiorOnly,
395     I N_edge, S_edge, E_edge, W_edge,
396 jmc 1.1 I dTtracerLev(k), vTrans,
397     U smVol(1-OLx,1-OLy,k),
398     U smTr0(1-OLx,1-OLy,k),
399     U smTr(1-OLx,1-OLy,k,bi,bj,1),
400     U smTr(1-OLx,1-OLy,k,bi,bj,2),
401     U smTr(1-OLx,1-OLy,k,bi,bj,3),
402     U smTr(1-OLx,1-OLy,k,bi,bj,4),
403     U smTr(1-OLx,1-OLy,k,bi,bj,5),
404     U smTr(1-OLx,1-OLy,k,bi,bj,6),
405     U smTr(1-OLx,1-OLy,k,bi,bj,7),
406     U smTr(1-OLx,1-OLy,k,bi,bj,8),
407     U smTr(1-OLx,1-OLy,k,bi,bj,9),
408     O afy, myThid )
409 jmc 1.2 ELSE
410     STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
411     ENDIF
412 jmc 1.1
413     #ifdef ALLOW_OBCS
414     C- Apply open boundary conditions
415 jmc 1.2 c IF (useOBCS) THEN
416 jmc 1.1 ccc localTij(i,j) = smTr0(i,j)/smVol(i,j)
417 jmc 1.2 c IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
418     c CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
419     c ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
420     c CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
421 jmc 1.1 #ifdef ALLOW_PTRACERS
422 jmc 1.2 c ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
423     c CALL OBCS_APPLY_PTRACER( bi, bj, k,
424     c & tracerIdentity-GAD_TR1+1, localTij, myThid )
425 jmc 1.1 #endif /* ALLOW_PTRACERS */
426 jmc 1.2 c ENDIF
427 jmc 1.1 ccc smTr0(i,j) = localTij(i,j)*smVol(i,j)
428 jmc 1.2 c ENDIF
429 jmc 1.1 #endif /* ALLOW_OBCS */
430    
431     C-- End of Y direction
432 jmc 1.2 ENDIF
433 jmc 1.1
434     C-- End of ipass loop
435 jmc 1.2 ENDDO
436 jmc 1.1
437 jmc 1.2 IF ( implicitAdvection ) THEN
438 jmc 1.1 C- explicit advection is done ; store tendency in gTracer:
439     DO j=1-OLy,sNy+OLy
440     DO i=1-OLx,sNx+OLx
441     localTr = smTr0(i,j,k)/smVol(i,j,k)
442     gTracer(i,j,k,bi,bj) = ( localTr - tracer(i,j,k,bi,bj) )
443     & / dTtracerLev(k)
444     ENDDO
445     ENDDO
446 jmc 1.2 ELSE
447 jmc 1.1 C- horizontal advection done; store intermediate result in 3D array:
448     c DO j=1-OLy,sNy+OLy
449     c DO i=1-OLx,sNx+OLx
450     c localTijk(i,j,k)=localTij(i,j)
451     c ENDDO
452     c ENDDO
453 jmc 1.2 ENDIF
454 jmc 1.1
455     #ifdef ALLOW_DIAGNOSTICS
456 jmc 1.2 IF ( doDiagAdvX ) THEN
457     diagName = 'ADVx'//diagSufx
458     CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid )
459     ENDIF
460     IF ( doDiagAdvY ) THEN
461     diagName = 'ADVy'//diagSufx
462     CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid )
463     ENDIF
464 jmc 1.1 #endif
465    
466     #ifdef ALLOW_DEBUG
467 jmc 1.2 IF ( debugLevel .GE. debLevB
468 jmc 1.1 & .AND. tracerIdentity.EQ.GAD_TEMPERATURE
469     & .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
470     & .AND. nPx.EQ.1 .AND. nPy.EQ.1
471     & .AND. useCubedSphereExchange ) THEN
472     CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_SOM_ADVECT',
473     & afx,afy, k, standardMessageUnit,bi,bj,myThid )
474 jmc 1.2 ENDIF
475 jmc 1.1 #endif /* ALLOW_DEBUG */
476    
477     C-- End of K loop for horizontal fluxes
478     ENDDO
479    
480     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
481    
482     IF ( .NOT.implicitAdvection ) THEN
483     C-- Apply limiter (if any):
484     CALL GAD_SOM_LIM_R( bi,bj, limiter,
485     U smVol,
486     U smTr0,
487     U smTr(1-OLx,1-OLy,1,bi,bj,1),
488     U smTr(1-OLx,1-OLy,1,bi,bj,2),
489     U smTr(1-OLx,1-OLy,1,bi,bj,3),
490     U smTr(1-OLx,1-OLy,1,bi,bj,4),
491     U smTr(1-OLx,1-OLy,1,bi,bj,5),
492     U smTr(1-OLx,1-OLy,1,bi,bj,6),
493     U smTr(1-OLx,1-OLy,1,bi,bj,7),
494     U smTr(1-OLx,1-OLy,1,bi,bj,8),
495     U smTr(1-OLx,1-OLy,1,bi,bj,9),
496     I myThid )
497    
498     C-- Start of k loop for vertical flux
499     DO k=Nr,1,-1
500     C-- kUp Cycles through 1,2 to point to w-layer above
501     C-- kDown Cycles through 2,1 to point to w-layer below
502     kUp = 1+MOD(Nr-k,2)
503     kDown= 1+MOD(Nr-k+1,2)
504     IF (k.EQ.Nr) THEN
505     C-- Set advective fluxes at the very bottom:
506     DO j=1-OLy,sNy+OLy
507     DO i=1-OLx,sNx+OLx
508     alp (i,j,kDown) = 0. _d 0
509     aln (i,j,kDown) = 0. _d 0
510     fp_v (i,j,kDown) = 0. _d 0
511     fn_v (i,j,kDown) = 0. _d 0
512     fp_o (i,j,kDown) = 0. _d 0
513     fn_o (i,j,kDown) = 0. _d 0
514     fp_x (i,j,kDown) = 0. _d 0
515     fn_x (i,j,kDown) = 0. _d 0
516     fp_y (i,j,kDown) = 0. _d 0
517     fn_y (i,j,kDown) = 0. _d 0
518     fp_z (i,j,kDown) = 0. _d 0
519     fn_z (i,j,kDown) = 0. _d 0
520     fp_xx(i,j,kDown) = 0. _d 0
521     fn_xx(i,j,kDown) = 0. _d 0
522     fp_yy(i,j,kDown) = 0. _d 0
523     fn_yy(i,j,kDown) = 0. _d 0
524     fp_zz(i,j,kDown) = 0. _d 0
525     fn_zz(i,j,kDown) = 0. _d 0
526     fp_xy(i,j,kDown) = 0. _d 0
527     fn_xy(i,j,kDown) = 0. _d 0
528     fp_xz(i,j,kDown) = 0. _d 0
529     fn_xz(i,j,kDown) = 0. _d 0
530     fp_yz(i,j,kDown) = 0. _d 0
531     fn_yz(i,j,kDown) = 0. _d 0
532     ENDDO
533     ENDDO
534     ENDIF
535    
536     C-- Compute Vertical transport
537     #ifdef ALLOW_AIM
538     C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
539     c IF ( k.EQ.1 .OR.
540     c & (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
541     c & ) THEN
542     #else
543     c IF ( k.EQ.1 ) THEN
544     #endif
545     IF ( (rigidLid.OR.nonlinFreeSurf.GE.1) .AND. k.EQ.1 ) THEN
546     C- Surface interface :
547     DO j=1-OLy,sNy+OLy
548     DO i=1-OLx,sNx+OLx
549     wFld(i,j) = 0.
550     rTrans(i,j) = 0.
551     maskUp(i,j) = 0.
552     ENDDO
553     ENDDO
554    
555     ELSEIF ( rigidLid.OR.nonlinFreeSurf.GE.1 ) THEN
556     C- Interior interface :
557     DO j=1-OLy,sNy+OLy
558     DO i=1-OLx,sNx+OLx
559     wFld(i,j) = wVel(i,j,k,bi,bj)
560     rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
561     & *deepFac2F(k)*rhoFacF(k)
562     & *maskC(i,j,k-1,bi,bj)
563     maskUp(i,j) = 1.
564     ENDDO
565     ENDDO
566    
567     ELSE
568     C- Linear Free-Surface: do not mask rTrans :
569     km1= MAX(k-1,1)
570     DO j=1-OLy,sNy+OLy
571     DO i=1-OLx,sNx+OLx
572     wFld(i,j) = wVel(i,j,k,bi,bj)
573     rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
574     & *deepFac2F(k)*rhoFacF(k)
575     maskUp(i,j) = maskC(i,j,km1,bi,bj)*maskC(i,j,k,bi,bj)
576     ENDDO
577     ENDDO
578    
579     C- end Surface/Interior if bloc
580     ENDIF
581    
582     #ifdef ALLOW_GMREDI
583     C-- Residual transp = Bolus transp + Eulerian transp
584 jmc 1.2 IF (useGMRedi .AND. k.GT.1 )
585 jmc 1.1 & CALL GMREDI_CALC_WFLOW(
586     U wFld, rTrans,
587     I k, bi, bj, myThid )
588     #endif /* ALLOW_GMREDI */
589    
590     C- Compute vertical advective flux in the interior:
591     IF ( vertAdvecScheme.EQ.ENUM_SOM_PRATHER
592 jmc 1.2 & .OR. vertAdvecScheme.EQ.ENUM_SOM_LIMITER ) THEN
593     CALL GAD_SOM_ADV_R(
594 jmc 1.1 I bi,bj,k, kUp, kDown,
595     I dTtracerLev(k), rTrans, maskUp,
596     U smVol,
597     U smTr0,
598     U smTr(1-OLx,1-OLy,1,bi,bj,1),
599     U smTr(1-OLx,1-OLy,1,bi,bj,2),
600     U smTr(1-OLx,1-OLy,1,bi,bj,3),
601     U smTr(1-OLx,1-OLy,1,bi,bj,4),
602     U smTr(1-OLx,1-OLy,1,bi,bj,5),
603     U smTr(1-OLx,1-OLy,1,bi,bj,6),
604     U smTr(1-OLx,1-OLy,1,bi,bj,7),
605     U smTr(1-OLx,1-OLy,1,bi,bj,8),
606     U smTr(1-OLx,1-OLy,1,bi,bj,9),
607     U alp, aln, fp_v, fn_v, fp_o, fn_o,
608     U fp_x, fn_x, fp_y, fn_y, fp_z, fn_z,
609     U fp_xx, fn_xx, fp_yy, fn_yy, fp_zz, fn_zz,
610     U fp_xy, fn_xy, fp_xz, fn_xz, fp_yz, fn_yz,
611     O afr, myThid )
612     ELSE
613 jmc 1.2 STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
614 jmc 1.1 ENDIF
615    
616 jmc 1.2 C-- Compute new tracer value and store tracer tendency
617 jmc 1.1 DO j=1-OLy,sNy+OLy
618     DO i=1-OLx,sNx+OLx
619     localTr = smTr0(i,j,k)
620     & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
621     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
622     & *recip_rhoFacC(k)
623     c localTr = smTr0(i,j,k)/smVol(i,j,k)
624     gTracer(i,j,k,bi,bj) = ( localTr - tracer(i,j,k,bi,bj) )
625     & / dTtracerLev(k)
626     ENDDO
627     ENDDO
628    
629     #ifdef ALLOW_DIAGNOSTICS
630 jmc 1.2 IF ( doDiagAdvR ) THEN
631     diagName = 'ADVr'//diagSufx
632     CALL DIAGNOSTICS_FILL( afr,
633     & diagName, k,1, 2,bi,bj, myThid )
634 jmc 1.1 ENDIF
635     #endif
636    
637 jmc 1.2 C-- End of k loop for vertical flux
638 jmc 1.1 ENDDO
639     C-- end of if not.implicitAdvection block
640     ENDIF
641    
642     RETURN
643     END
  ViewVC Help
Powered by ViewVC 1.1.22