/[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.6 - (hide annotations) (download)
Tue May 12 19:56:35 2009 UTC (15 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint61o, checkpoint61r, checkpoint61p, checkpoint61q
Changes since 1.5: +2 -2 lines
new header file "W2_EXCH2_SIZE.h" coming with new W2-Exch2 topology code

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

  ViewVC Help
Powered by ViewVC 1.1.22