/[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.3 - (hide annotations) (download)
Fri May 9 21:43:16 2008 UTC (16 years, 1 month ago) by jmc
Branch: MAIN
Changes since 1.2: +1 -4 lines
change option: GAD_ALLOW_SOM_ADVECT to GAD_ALLOW_TS_SOM_ADV which only
applies to files where Temperature & salinity 2nd Order moments are used

1 jmc 1.3 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_som_advect.F,v 1.2 2008/02/12 20:32: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     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_AUTODIFF_TAMC
223     IF ( npass.GT.maxcube ) STOP 'maxcube needs to be = 3'
224     #endif
225     #ifdef ALLOW_EXCH2
226     myTile = W2_myTileList(bi)
227     nCFace = exch2_myFace(myTile)
228     N_edge = exch2_isNedge(myTile).EQ.1
229     S_edge = exch2_isSedge(myTile).EQ.1
230     E_edge = exch2_isEedge(myTile).EQ.1
231     W_edge = exch2_isWedge(myTile).EQ.1
232     #else
233     nCFace = bi
234     N_edge = .TRUE.
235     S_edge = .TRUE.
236     E_edge = .TRUE.
237     W_edge = .TRUE.
238     #endif
239     ELSE
240     npass = 2
241     nCFace = 0
242     N_edge = .FALSE.
243     S_edge = .FALSE.
244     E_edge = .FALSE.
245     W_edge = .FALSE.
246     ENDIF
247    
248 jmc 1.1 limiter = MOD(advectionScheme, 10)
249    
250     C-- Start of k loop for horizontal fluxes
251     DO k=1,Nr
252    
253     C-- Get temporary terms used by tendency routines
254 jmc 1.2 CALL CALC_COMMON_FACTORS (
255 jmc 1.1 I uVel, vVel,
256     O uFld, vFld, uTrans, vTrans, xA, yA,
257     I k,bi,bj, myThid )
258    
259     #ifdef ALLOW_GMREDI
260     C-- Residual transp = Bolus transp + Eulerian transp
261 jmc 1.2 IF (useGMRedi)
262 jmc 1.1 & CALL GMREDI_CALC_UVFLOW(
263     U uFld, vFld, uTrans, vTrans,
264     I k, bi, bj, myThid )
265     #endif /* ALLOW_GMREDI */
266    
267     C-- grid-box volume and tracer content (zero order moment)
268 jmc 1.2 DO j=1-OLy,sNy+OLy
269     DO i=1-OLx,sNx+OLx
270 jmc 1.1 smVol(i,j,k) = rA(i,j,bi,bj)*deepFac2C(k)
271     & *drF(k)*hFacC(i,j,k,bi,bj)
272     & *rhoFacC(k)
273     smTr0(i,j,k) = tracer(i,j,k,bi,bj)*smVol(i,j,k)
274     C- fill empty grid-box:
275     smVol(i,j,k) = smVol(i,j,k)
276     & + (1. _d 0 - maskC(i,j,k,bi,bj))
277 jmc 1.2 ENDDO
278 jmc 1.1 ENDDO
279    
280     C-- Multiple passes for different directions on different tiles
281     C-- For cube need one pass for each of red, green and blue axes.
282 jmc 1.2 DO ipass=1,npass
283 jmc 1.1
284 jmc 1.2 interiorOnly = .FALSE.
285     overlapOnly = .FALSE.
286     IF (useCubedSphereExchange) THEN
287     C- CubedSphere : pass 3 times, with partial update of local tracer field
288     IF (ipass.EQ.1) THEN
289     overlapOnly = MOD(nCFace,3).EQ.0
290     interiorOnly = MOD(nCFace,3).NE.0
291     calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
292     calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
293     ELSEIF (ipass.EQ.2) THEN
294     overlapOnly = MOD(nCFace,3).EQ.2
295     interiorOnly = MOD(nCFace,3).EQ.1
296     calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
297     calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
298     ELSE
299     interiorOnly = .TRUE.
300     calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
301     calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
302     ENDIF
303     ELSE
304 jmc 1.1 C- not CubedSphere
305 jmc 1.2 calc_fluxes_X = MOD(ipass,2).EQ.1
306     calc_fluxes_Y = .NOT.calc_fluxes_X
307     ENDIF
308 jmc 1.1
309     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
310 jmc 1.2
311 jmc 1.1 C-- X direction
312 jmc 1.2 C- Do not compute fluxes if
313     C a) needed in overlap only
314     C and b) the overlap of myTile are not cube-face Edges
315     IF ( calc_fluxes_X .AND.
316     & (.NOT.overlapOnly .OR. N_edge .OR. S_edge)
317     & ) THEN
318    
319     C- Internal exchange for calculations in X
320     IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN
321     CALL GAD_SOM_PREP_CS_CORNER(
322     U smVol, smTr0, smTr, smCorners,
323     I .TRUE., overlapOnly, interiorOnly,
324     I N_edge, S_edge, E_edge, W_edge,
325     I ipass, k, Nr, bi, bj, myThid )
326     ENDIF
327    
328     C- Solve advection in X and update moments
329     IF ( advectionScheme.EQ.ENUM_SOM_PRATHER
330     & .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN
331     CALL GAD_SOM_ADV_X(
332 jmc 1.1 I bi,bj,k, limiter,
333 jmc 1.2 I overlapOnly, interiorOnly,
334     I N_edge, S_edge, E_edge, W_edge,
335 jmc 1.1 I dTtracerLev(k), uTrans,
336     U smVol(1-OLx,1-OLy,k),
337     U smTr0(1-OLx,1-OLy,k),
338     U smTr(1-OLx,1-OLy,k,bi,bj,1),
339     U smTr(1-OLx,1-OLy,k,bi,bj,2),
340     U smTr(1-OLx,1-OLy,k,bi,bj,3),
341     U smTr(1-OLx,1-OLy,k,bi,bj,4),
342     U smTr(1-OLx,1-OLy,k,bi,bj,5),
343     U smTr(1-OLx,1-OLy,k,bi,bj,6),
344     U smTr(1-OLx,1-OLy,k,bi,bj,7),
345     U smTr(1-OLx,1-OLy,k,bi,bj,8),
346     U smTr(1-OLx,1-OLy,k,bi,bj,9),
347     O afx, myThid )
348 jmc 1.2 ELSE
349     STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
350     ENDIF
351 jmc 1.1
352     #ifdef ALLOW_OBCS
353     C- Apply open boundary conditions
354 jmc 1.2 c IF ( useOBCS ) THEN
355 jmc 1.1 ccc localTij(i,j) = smTr0(i,j)/smVol(i,j)
356 jmc 1.2 c IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
357     c CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
358     c ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
359     c CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
360 jmc 1.1 #ifdef ALLOW_PTRACERS
361 jmc 1.2 c ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
362     c CALL OBCS_APPLY_PTRACER( bi, bj, k,
363     c & tracerIdentity-GAD_TR1+1, localTij, myThid )
364 jmc 1.1 #endif /* ALLOW_PTRACERS */
365 jmc 1.2 c ENDIF
366 jmc 1.1 ccc smTr0(i,j) = localTij(i,j)*smVol(i,j)
367 jmc 1.2 c ENDIF
368 jmc 1.1 #endif /* ALLOW_OBCS */
369    
370     C-- End of X direction
371 jmc 1.2 ENDIF
372 jmc 1.1
373     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
374 jmc 1.2
375 jmc 1.1 C-- Y direction
376     C- Do not compute fluxes if
377     C a) needed in overlap only
378     C and b) the overlap of myTile are not cube-face edges
379 jmc 1.2 IF ( calc_fluxes_Y .AND.
380     & (.NOT.overlapOnly .OR. E_edge .OR. W_edge)
381     & ) THEN
382    
383     C- Internal exchange for calculations in Y
384     IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN
385     CALL GAD_SOM_PREP_CS_CORNER(
386     U smVol, smTr0, smTr, smCorners,
387     I .FALSE., overlapOnly, interiorOnly,
388     I N_edge, S_edge, E_edge, W_edge,
389     I iPass, k, Nr, bi, bj, myThid )
390     ENDIF
391    
392     C- Solve advection in Y and update moments
393     IF ( advectionScheme.EQ.ENUM_SOM_PRATHER
394     & .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN
395     CALL GAD_SOM_ADV_Y(
396 jmc 1.1 I bi,bj,k, limiter,
397 jmc 1.2 I overlapOnly, interiorOnly,
398     I N_edge, S_edge, E_edge, W_edge,
399 jmc 1.1 I dTtracerLev(k), vTrans,
400     U smVol(1-OLx,1-OLy,k),
401     U smTr0(1-OLx,1-OLy,k),
402     U smTr(1-OLx,1-OLy,k,bi,bj,1),
403     U smTr(1-OLx,1-OLy,k,bi,bj,2),
404     U smTr(1-OLx,1-OLy,k,bi,bj,3),
405     U smTr(1-OLx,1-OLy,k,bi,bj,4),
406     U smTr(1-OLx,1-OLy,k,bi,bj,5),
407     U smTr(1-OLx,1-OLy,k,bi,bj,6),
408     U smTr(1-OLx,1-OLy,k,bi,bj,7),
409     U smTr(1-OLx,1-OLy,k,bi,bj,8),
410     U smTr(1-OLx,1-OLy,k,bi,bj,9),
411     O afy, myThid )
412 jmc 1.2 ELSE
413     STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
414     ENDIF
415 jmc 1.1
416     #ifdef ALLOW_OBCS
417     C- Apply open boundary conditions
418 jmc 1.2 c IF (useOBCS) THEN
419 jmc 1.1 ccc localTij(i,j) = smTr0(i,j)/smVol(i,j)
420 jmc 1.2 c IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
421     c CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
422     c ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
423     c CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
424 jmc 1.1 #ifdef ALLOW_PTRACERS
425 jmc 1.2 c ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
426     c CALL OBCS_APPLY_PTRACER( bi, bj, k,
427     c & tracerIdentity-GAD_TR1+1, localTij, myThid )
428 jmc 1.1 #endif /* ALLOW_PTRACERS */
429 jmc 1.2 c ENDIF
430 jmc 1.1 ccc smTr0(i,j) = localTij(i,j)*smVol(i,j)
431 jmc 1.2 c ENDIF
432 jmc 1.1 #endif /* ALLOW_OBCS */
433    
434     C-- End of Y direction
435 jmc 1.2 ENDIF
436 jmc 1.1
437     C-- End of ipass loop
438 jmc 1.2 ENDDO
439 jmc 1.1
440 jmc 1.2 IF ( implicitAdvection ) THEN
441 jmc 1.1 C- explicit advection is done ; store tendency in gTracer:
442     DO j=1-OLy,sNy+OLy
443     DO i=1-OLx,sNx+OLx
444     localTr = smTr0(i,j,k)/smVol(i,j,k)
445     gTracer(i,j,k,bi,bj) = ( localTr - tracer(i,j,k,bi,bj) )
446     & / dTtracerLev(k)
447     ENDDO
448     ENDDO
449 jmc 1.2 ELSE
450 jmc 1.1 C- horizontal advection done; store intermediate result in 3D array:
451     c DO j=1-OLy,sNy+OLy
452     c DO i=1-OLx,sNx+OLx
453     c localTijk(i,j,k)=localTij(i,j)
454     c ENDDO
455     c ENDDO
456 jmc 1.2 ENDIF
457 jmc 1.1
458     #ifdef ALLOW_DIAGNOSTICS
459 jmc 1.2 IF ( doDiagAdvX ) THEN
460     diagName = 'ADVx'//diagSufx
461     CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid )
462     ENDIF
463     IF ( doDiagAdvY ) THEN
464     diagName = 'ADVy'//diagSufx
465     CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid )
466     ENDIF
467 jmc 1.1 #endif
468    
469     #ifdef ALLOW_DEBUG
470 jmc 1.2 IF ( debugLevel .GE. debLevB
471 jmc 1.1 & .AND. tracerIdentity.EQ.GAD_TEMPERATURE
472     & .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
473     & .AND. nPx.EQ.1 .AND. nPy.EQ.1
474     & .AND. useCubedSphereExchange ) THEN
475     CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_SOM_ADVECT',
476     & afx,afy, k, standardMessageUnit,bi,bj,myThid )
477 jmc 1.2 ENDIF
478 jmc 1.1 #endif /* ALLOW_DEBUG */
479    
480     C-- End of K loop for horizontal fluxes
481     ENDDO
482    
483     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
484    
485     IF ( .NOT.implicitAdvection ) THEN
486     C-- Apply limiter (if any):
487     CALL GAD_SOM_LIM_R( bi,bj, limiter,
488     U smVol,
489     U smTr0,
490     U smTr(1-OLx,1-OLy,1,bi,bj,1),
491     U smTr(1-OLx,1-OLy,1,bi,bj,2),
492     U smTr(1-OLx,1-OLy,1,bi,bj,3),
493     U smTr(1-OLx,1-OLy,1,bi,bj,4),
494     U smTr(1-OLx,1-OLy,1,bi,bj,5),
495     U smTr(1-OLx,1-OLy,1,bi,bj,6),
496     U smTr(1-OLx,1-OLy,1,bi,bj,7),
497     U smTr(1-OLx,1-OLy,1,bi,bj,8),
498     U smTr(1-OLx,1-OLy,1,bi,bj,9),
499     I myThid )
500    
501     C-- Start of k loop for vertical flux
502     DO k=Nr,1,-1
503     C-- kUp Cycles through 1,2 to point to w-layer above
504     C-- kDown Cycles through 2,1 to point to w-layer below
505     kUp = 1+MOD(Nr-k,2)
506     kDown= 1+MOD(Nr-k+1,2)
507     IF (k.EQ.Nr) THEN
508     C-- Set advective fluxes at the very bottom:
509     DO j=1-OLy,sNy+OLy
510     DO i=1-OLx,sNx+OLx
511     alp (i,j,kDown) = 0. _d 0
512     aln (i,j,kDown) = 0. _d 0
513     fp_v (i,j,kDown) = 0. _d 0
514     fn_v (i,j,kDown) = 0. _d 0
515     fp_o (i,j,kDown) = 0. _d 0
516     fn_o (i,j,kDown) = 0. _d 0
517     fp_x (i,j,kDown) = 0. _d 0
518     fn_x (i,j,kDown) = 0. _d 0
519     fp_y (i,j,kDown) = 0. _d 0
520     fn_y (i,j,kDown) = 0. _d 0
521     fp_z (i,j,kDown) = 0. _d 0
522     fn_z (i,j,kDown) = 0. _d 0
523     fp_xx(i,j,kDown) = 0. _d 0
524     fn_xx(i,j,kDown) = 0. _d 0
525     fp_yy(i,j,kDown) = 0. _d 0
526     fn_yy(i,j,kDown) = 0. _d 0
527     fp_zz(i,j,kDown) = 0. _d 0
528     fn_zz(i,j,kDown) = 0. _d 0
529     fp_xy(i,j,kDown) = 0. _d 0
530     fn_xy(i,j,kDown) = 0. _d 0
531     fp_xz(i,j,kDown) = 0. _d 0
532     fn_xz(i,j,kDown) = 0. _d 0
533     fp_yz(i,j,kDown) = 0. _d 0
534     fn_yz(i,j,kDown) = 0. _d 0
535     ENDDO
536     ENDDO
537     ENDIF
538    
539     C-- Compute Vertical transport
540     #ifdef ALLOW_AIM
541     C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
542     c IF ( k.EQ.1 .OR.
543     c & (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
544     c & ) THEN
545     #else
546     c IF ( k.EQ.1 ) THEN
547     #endif
548     IF ( (rigidLid.OR.nonlinFreeSurf.GE.1) .AND. k.EQ.1 ) THEN
549     C- Surface interface :
550     DO j=1-OLy,sNy+OLy
551     DO i=1-OLx,sNx+OLx
552     wFld(i,j) = 0.
553     rTrans(i,j) = 0.
554     maskUp(i,j) = 0.
555     ENDDO
556     ENDDO
557    
558     ELSEIF ( rigidLid.OR.nonlinFreeSurf.GE.1 ) THEN
559     C- Interior interface :
560     DO j=1-OLy,sNy+OLy
561     DO i=1-OLx,sNx+OLx
562     wFld(i,j) = wVel(i,j,k,bi,bj)
563     rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
564     & *deepFac2F(k)*rhoFacF(k)
565     & *maskC(i,j,k-1,bi,bj)
566     maskUp(i,j) = 1.
567     ENDDO
568     ENDDO
569    
570     ELSE
571     C- Linear Free-Surface: do not mask rTrans :
572     km1= MAX(k-1,1)
573     DO j=1-OLy,sNy+OLy
574     DO i=1-OLx,sNx+OLx
575     wFld(i,j) = wVel(i,j,k,bi,bj)
576     rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
577     & *deepFac2F(k)*rhoFacF(k)
578     maskUp(i,j) = maskC(i,j,km1,bi,bj)*maskC(i,j,k,bi,bj)
579     ENDDO
580     ENDDO
581    
582     C- end Surface/Interior if bloc
583     ENDIF
584    
585     #ifdef ALLOW_GMREDI
586     C-- Residual transp = Bolus transp + Eulerian transp
587 jmc 1.2 IF (useGMRedi .AND. k.GT.1 )
588 jmc 1.1 & CALL GMREDI_CALC_WFLOW(
589     U wFld, rTrans,
590     I k, bi, bj, myThid )
591     #endif /* ALLOW_GMREDI */
592    
593     C- Compute vertical advective flux in the interior:
594     IF ( vertAdvecScheme.EQ.ENUM_SOM_PRATHER
595 jmc 1.2 & .OR. vertAdvecScheme.EQ.ENUM_SOM_LIMITER ) THEN
596     CALL GAD_SOM_ADV_R(
597 jmc 1.1 I bi,bj,k, kUp, kDown,
598     I dTtracerLev(k), rTrans, maskUp,
599     U smVol,
600     U smTr0,
601     U smTr(1-OLx,1-OLy,1,bi,bj,1),
602     U smTr(1-OLx,1-OLy,1,bi,bj,2),
603     U smTr(1-OLx,1-OLy,1,bi,bj,3),
604     U smTr(1-OLx,1-OLy,1,bi,bj,4),
605     U smTr(1-OLx,1-OLy,1,bi,bj,5),
606     U smTr(1-OLx,1-OLy,1,bi,bj,6),
607     U smTr(1-OLx,1-OLy,1,bi,bj,7),
608     U smTr(1-OLx,1-OLy,1,bi,bj,8),
609     U smTr(1-OLx,1-OLy,1,bi,bj,9),
610     U alp, aln, fp_v, fn_v, fp_o, fn_o,
611     U fp_x, fn_x, fp_y, fn_y, fp_z, fn_z,
612     U fp_xx, fn_xx, fp_yy, fn_yy, fp_zz, fn_zz,
613     U fp_xy, fn_xy, fp_xz, fn_xz, fp_yz, fn_yz,
614     O afr, myThid )
615     ELSE
616 jmc 1.2 STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
617 jmc 1.1 ENDIF
618    
619 jmc 1.2 C-- Compute new tracer value and store tracer tendency
620 jmc 1.1 DO j=1-OLy,sNy+OLy
621     DO i=1-OLx,sNx+OLx
622     localTr = smTr0(i,j,k)
623     & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
624     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
625     & *recip_rhoFacC(k)
626     c localTr = smTr0(i,j,k)/smVol(i,j,k)
627     gTracer(i,j,k,bi,bj) = ( localTr - tracer(i,j,k,bi,bj) )
628     & / dTtracerLev(k)
629     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