/[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.7 - (hide annotations) (download)
Fri Jun 26 23:10:09 2009 UTC (14 years, 11 months ago) by jahn
Branch: MAIN
Changes since 1.6: +10 -9 lines
add package longstep

1 jahn 1.7 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_som_advect.F,v 1.6 2009/05/12 19:56:35 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 jahn 1.7 I tracerIdentity, deltaTLev,
14 jmc 1.1 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 jahn 1.7 _RL deltaTLev(Nr)
55 jmc 1.1 _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
56     _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
57     _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
58     _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
59     INTEGER bi,bj
60     _RL myTime
61     INTEGER myIter
62     INTEGER myThid
63    
64     C !OUTPUT PARAMETERS: ==================================================
65     C smTr :: tracer 1rst & 2nd Order moments
66     C gTracer :: tendency array
67     _RL smTr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy,nSOM)
68     _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
69    
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 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     _RL smVol (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
111     _RL smTr0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
112     _RL alp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
113     _RL aln (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
114     _RL fp_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
115     _RL fn_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
116     _RL fp_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
117     _RL fn_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
118     _RL fp_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
119     _RL fn_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
120     _RL fp_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
121     _RL fn_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
122     _RL fp_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
123     _RL fn_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
124     _RL fp_xx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
125     _RL fn_xx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
126     _RL fp_yy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
127     _RL fn_yy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
128     _RL fp_zz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
129     _RL fn_zz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
130     _RL fp_xy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
131     _RL fn_xy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
132     _RL fp_xz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
133     _RL fn_xz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
134     _RL fp_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
135     _RL fn_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
136 jmc 1.2 _RL smCorners(OLx,OLy,4,-1:nSOM)
137 jmc 1.5 c _RL localTr
138 jmc 1.1 LOGICAL calc_fluxes_X, calc_fluxes_Y
139 jmc 1.2 LOGICAL interiorOnly, overlapOnly
140 jmc 1.1 INTEGER limiter
141     INTEGER npass, ipass
142 jmc 1.2 INTEGER nCFace, n
143     LOGICAL N_edge, S_edge, E_edge, W_edge
144 jmc 1.1 CHARACTER*(MAX_LEN_MBUF) msgBuf
145 jmc 1.2 #ifdef ALLOW_EXCH2
146     INTEGER myTile
147     #endif
148 jmc 1.1 #ifdef ALLOW_DIAGNOSTICS
149     CHARACTER*8 diagName
150 jmc 1.2 CHARACTER*4 diagSufx
151     LOGICAL doDiagAdvX, doDiagAdvY, doDiagAdvR
152     C- Functions:
153     CHARACTER*4 GAD_DIAG_SUFX
154 jmc 1.1 EXTERNAL GAD_DIAG_SUFX
155 jmc 1.2 LOGICAL DIAGNOSTICS_IS_ON
156     EXTERNAL DIAGNOSTICS_IS_ON
157 jmc 1.1 #endif
158     CEOP
159    
160     #ifdef ALLOW_DIAGNOSTICS
161 jmc 1.2 C-- Set diagnostics flags and suffix for the current tracer
162     doDiagAdvX = .FALSE.
163     doDiagAdvY = .FALSE.
164     doDiagAdvR = .FALSE.
165 jmc 1.1 IF ( useDiagnostics ) THEN
166     diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
167 jmc 1.2 diagName = 'ADVx'//diagSufx
168     doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )
169     diagName = 'ADVy'//diagSufx
170     doDiagAdvY = DIAGNOSTICS_IS_ON( diagName, myThid )
171     diagName = 'ADVr'//diagSufx
172     doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )
173 jmc 1.1 ENDIF
174     #endif
175    
176     C-- Set up work arrays with valid (i.e. not NaN) values
177     C These inital values do not alter the numerical results. They
178     C just ensure that all memory references are to valid floating
179     C point numbers. This prevents spurious hardware signals due to
180     C uninitialised but inert locations.
181 jmc 1.2 DO j=1-OLy,sNy+OLy
182     DO i=1-OLx,sNx+OLx
183     afx(i,j) = 0.
184     afy(i,j) = 0.
185 jmc 1.1 C- xA,yA,uFld,vFld,uTrans,vTrans are set over the full domain
186     C in CALC_COMMON_FACTORS: no need for extra initialisation
187     c xA(i,j) = 0. _d 0
188     c yA(i,j) = 0. _d 0
189     c uTrans(i,j) = 0. _d 0
190     c vTrans(i,j) = 0. _d 0
191     C- rTrans is set over the full domain: no need for extra initialisation
192     c rTrans(i,j) = 0. _d 0
193 jmc 1.2 ENDDO
194     ENDDO
195     DO n=-1,nSOM
196     DO k=1,4
197     DO j=1,OLy
198     DO i=1,OLx
199     smCorners(i,j,k,n) = 0.
200     ENDDO
201     ENDDO
202     ENDDO
203     ENDDO
204 jmc 1.1
205     IF ( implicitAdvection ) THEN
206     WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ',
207     & 'not coded for implicit-vertical Advection'
208     CALL PRINT_ERROR( msgBuf, myThid )
209     STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'
210     ENDIF
211     IF ( vertAdvecScheme .NE. advectionScheme ) THEN
212     WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ',
213     & 'not coded for different vertAdvecScheme'
214     CALL PRINT_ERROR( msgBuf, myThid )
215     STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT'
216     ENDIF
217    
218 jmc 1.2 C-- Set tile-specific parameters for horizontal fluxes
219     IF (useCubedSphereExchange) THEN
220     npass = 3
221     #ifdef ALLOW_EXCH2
222     myTile = W2_myTileList(bi)
223     nCFace = exch2_myFace(myTile)
224     N_edge = exch2_isNedge(myTile).EQ.1
225     S_edge = exch2_isSedge(myTile).EQ.1
226     E_edge = exch2_isEedge(myTile).EQ.1
227     W_edge = exch2_isWedge(myTile).EQ.1
228     #else
229     nCFace = bi
230     N_edge = .TRUE.
231     S_edge = .TRUE.
232     E_edge = .TRUE.
233     W_edge = .TRUE.
234     #endif
235     ELSE
236     npass = 2
237     nCFace = 0
238     N_edge = .FALSE.
239     S_edge = .FALSE.
240     E_edge = .FALSE.
241     W_edge = .FALSE.
242     ENDIF
243    
244 jmc 1.1 limiter = MOD(advectionScheme, 10)
245    
246     C-- Start of k loop for horizontal fluxes
247     DO k=1,Nr
248    
249     C-- Get temporary terms used by tendency routines
250 jmc 1.2 CALL CALC_COMMON_FACTORS (
251 jmc 1.1 I uVel, vVel,
252     O uFld, vFld, uTrans, vTrans, xA, yA,
253     I k,bi,bj, myThid )
254    
255     #ifdef ALLOW_GMREDI
256     C-- Residual transp = Bolus transp + Eulerian transp
257 jmc 1.2 IF (useGMRedi)
258 jmc 1.1 & CALL GMREDI_CALC_UVFLOW(
259     U uFld, vFld, uTrans, vTrans,
260     I k, bi, bj, myThid )
261     #endif /* ALLOW_GMREDI */
262    
263     C-- grid-box volume and tracer content (zero order moment)
264 jmc 1.2 DO j=1-OLy,sNy+OLy
265     DO i=1-OLx,sNx+OLx
266 jmc 1.1 smVol(i,j,k) = rA(i,j,bi,bj)*deepFac2C(k)
267     & *drF(k)*hFacC(i,j,k,bi,bj)
268     & *rhoFacC(k)
269     smTr0(i,j,k) = tracer(i,j,k,bi,bj)*smVol(i,j,k)
270     C- fill empty grid-box:
271     smVol(i,j,k) = smVol(i,j,k)
272     & + (1. _d 0 - maskC(i,j,k,bi,bj))
273 jmc 1.2 ENDDO
274 jmc 1.1 ENDDO
275    
276     C-- Multiple passes for different directions on different tiles
277     C-- For cube need one pass for each of red, green and blue axes.
278 jmc 1.2 DO ipass=1,npass
279 jmc 1.1
280 jmc 1.2 interiorOnly = .FALSE.
281     overlapOnly = .FALSE.
282     IF (useCubedSphereExchange) THEN
283     C- CubedSphere : pass 3 times, with partial update of local tracer field
284     IF (ipass.EQ.1) THEN
285     overlapOnly = MOD(nCFace,3).EQ.0
286     interiorOnly = MOD(nCFace,3).NE.0
287     calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
288     calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
289     ELSEIF (ipass.EQ.2) THEN
290     overlapOnly = MOD(nCFace,3).EQ.2
291     interiorOnly = MOD(nCFace,3).EQ.1
292     calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
293     calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
294     ELSE
295     interiorOnly = .TRUE.
296     calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
297     calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
298     ENDIF
299     ELSE
300 jmc 1.1 C- not CubedSphere
301 jmc 1.2 calc_fluxes_X = MOD(ipass,2).EQ.1
302     calc_fluxes_Y = .NOT.calc_fluxes_X
303     ENDIF
304 jmc 1.1
305     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
306 jmc 1.2
307 jmc 1.1 C-- X direction
308 jmc 1.2 C- Do not compute fluxes if
309     C a) needed in overlap only
310     C and b) the overlap of myTile are not cube-face Edges
311     IF ( calc_fluxes_X .AND.
312     & (.NOT.overlapOnly .OR. N_edge .OR. S_edge)
313     & ) THEN
314    
315     C- Internal exchange for calculations in X
316     IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN
317     CALL GAD_SOM_PREP_CS_CORNER(
318     U smVol, smTr0, smTr, smCorners,
319     I .TRUE., overlapOnly, interiorOnly,
320     I N_edge, S_edge, E_edge, W_edge,
321     I ipass, k, Nr, bi, bj, myThid )
322     ENDIF
323    
324     C- Solve advection in X and update moments
325     IF ( advectionScheme.EQ.ENUM_SOM_PRATHER
326     & .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN
327     CALL GAD_SOM_ADV_X(
328 jmc 1.1 I bi,bj,k, limiter,
329 jmc 1.2 I overlapOnly, interiorOnly,
330     I N_edge, S_edge, E_edge, W_edge,
331 jahn 1.7 I deltaTLev(k), uTrans,
332 jmc 1.1 U smVol(1-OLx,1-OLy,k),
333     U smTr0(1-OLx,1-OLy,k),
334     U smTr(1-OLx,1-OLy,k,bi,bj,1),
335     U smTr(1-OLx,1-OLy,k,bi,bj,2),
336     U smTr(1-OLx,1-OLy,k,bi,bj,3),
337     U smTr(1-OLx,1-OLy,k,bi,bj,4),
338     U smTr(1-OLx,1-OLy,k,bi,bj,5),
339     U smTr(1-OLx,1-OLy,k,bi,bj,6),
340     U smTr(1-OLx,1-OLy,k,bi,bj,7),
341     U smTr(1-OLx,1-OLy,k,bi,bj,8),
342     U smTr(1-OLx,1-OLy,k,bi,bj,9),
343     O afx, myThid )
344 jmc 1.2 ELSE
345     STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
346     ENDIF
347 jmc 1.1
348     #ifdef ALLOW_OBCS
349     C- Apply open boundary conditions
350 jmc 1.2 c IF ( useOBCS ) THEN
351 jmc 1.1 ccc localTij(i,j) = smTr0(i,j)/smVol(i,j)
352 jmc 1.2 c IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
353     c CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
354     c ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
355     c CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
356 jmc 1.1 #ifdef ALLOW_PTRACERS
357 jmc 1.2 c ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
358     c CALL OBCS_APPLY_PTRACER( bi, bj, k,
359     c & tracerIdentity-GAD_TR1+1, localTij, myThid )
360 jmc 1.1 #endif /* ALLOW_PTRACERS */
361 jmc 1.2 c ENDIF
362 jmc 1.1 ccc smTr0(i,j) = localTij(i,j)*smVol(i,j)
363 jmc 1.2 c ENDIF
364 jmc 1.1 #endif /* ALLOW_OBCS */
365    
366     C-- End of X direction
367 jmc 1.2 ENDIF
368 jmc 1.1
369     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
370 jmc 1.2
371 jmc 1.1 C-- Y direction
372     C- Do not compute fluxes if
373     C a) needed in overlap only
374     C and b) the overlap of myTile are not cube-face edges
375 jmc 1.2 IF ( calc_fluxes_Y .AND.
376     & (.NOT.overlapOnly .OR. E_edge .OR. W_edge)
377     & ) THEN
378    
379     C- Internal exchange for calculations in Y
380     IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN
381     CALL GAD_SOM_PREP_CS_CORNER(
382     U smVol, smTr0, smTr, smCorners,
383     I .FALSE., overlapOnly, interiorOnly,
384     I N_edge, S_edge, E_edge, W_edge,
385     I iPass, k, Nr, bi, bj, myThid )
386     ENDIF
387    
388     C- Solve advection in Y and update moments
389     IF ( advectionScheme.EQ.ENUM_SOM_PRATHER
390     & .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN
391     CALL GAD_SOM_ADV_Y(
392 jmc 1.1 I bi,bj,k, limiter,
393 jmc 1.2 I overlapOnly, interiorOnly,
394     I N_edge, S_edge, E_edge, W_edge,
395 jahn 1.7 I deltaTLev(k), vTrans,
396 jmc 1.1 U smVol(1-OLx,1-OLy,k),
397     U smTr0(1-OLx,1-OLy,k),
398     U smTr(1-OLx,1-OLy,k,bi,bj,1),
399     U smTr(1-OLx,1-OLy,k,bi,bj,2),
400     U smTr(1-OLx,1-OLy,k,bi,bj,3),
401     U smTr(1-OLx,1-OLy,k,bi,bj,4),
402     U smTr(1-OLx,1-OLy,k,bi,bj,5),
403     U smTr(1-OLx,1-OLy,k,bi,bj,6),
404     U smTr(1-OLx,1-OLy,k,bi,bj,7),
405     U smTr(1-OLx,1-OLy,k,bi,bj,8),
406     U smTr(1-OLx,1-OLy,k,bi,bj,9),
407     O afy, myThid )
408 jmc 1.2 ELSE
409     STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
410     ENDIF
411 jmc 1.1
412     #ifdef ALLOW_OBCS
413     C- Apply open boundary conditions
414 jmc 1.2 c IF (useOBCS) THEN
415 jmc 1.1 ccc localTij(i,j) = smTr0(i,j)/smVol(i,j)
416 jmc 1.2 c IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
417     c CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
418     c ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
419     c CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
420 jmc 1.1 #ifdef ALLOW_PTRACERS
421 jmc 1.2 c ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
422     c CALL OBCS_APPLY_PTRACER( bi, bj, k,
423     c & tracerIdentity-GAD_TR1+1, localTij, myThid )
424 jmc 1.1 #endif /* ALLOW_PTRACERS */
425 jmc 1.2 c ENDIF
426 jmc 1.1 ccc smTr0(i,j) = localTij(i,j)*smVol(i,j)
427 jmc 1.2 c ENDIF
428 jmc 1.1 #endif /* ALLOW_OBCS */
429    
430     C-- End of Y direction
431 jmc 1.2 ENDIF
432 jmc 1.1
433     C-- End of ipass loop
434 jmc 1.2 ENDDO
435 jmc 1.1
436 jmc 1.2 IF ( implicitAdvection ) THEN
437 jmc 1.1 C- explicit advection is done ; store tendency in gTracer:
438     DO j=1-OLy,sNy+OLy
439     DO i=1-OLx,sNx+OLx
440 jmc 1.5 C-- without rescaling of tendencies:
441     c localTr = smTr0(i,j,k)/smVol(i,j,k)
442     c gTracer(i,j,k,bi,bj) = ( localTr - tracer(i,j,k,bi,bj) )
443 jahn 1.7 c & / deltaTLev(k)
444 jmc 1.5 C-- consistent with rescaling of tendencies (in FREESURF_RESCALE_G):
445     gTracer(i,j,k,bi,bj) =
446     & ( smTr0(i,j,k) - tracer(i,j,k,bi,bj)*smVol(i,j,k) )
447     & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
448     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
449     & *recip_rhoFacC(k)
450 jahn 1.7 & /deltaTLev(k)
451 jmc 1.1 ENDDO
452     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 jahn 1.7 I deltaTLev(k), rTrans, maskUp,
596 jmc 1.1 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 jmc 1.5 C-- without rescaling of tendencies:
620     c localTr = smTr0(i,j,k)/smVol(i,j,k)
621     c gTracer(i,j,k,bi,bj) = ( localTr - tracer(i,j,k,bi,bj) )
622 jahn 1.7 c & / deltaTLev(k)
623 jmc 1.5 C-- consistent with rescaling of tendencies (in FREESURF_RESCALE_G):
624     gTracer(i,j,k,bi,bj) =
625     & ( smTr0(i,j,k) - tracer(i,j,k,bi,bj)*smVol(i,j,k) )
626 jmc 1.1 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
627     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
628     & *recip_rhoFacC(k)
629 jahn 1.7 & /deltaTLev(k)
630 jmc 1.1 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     RETURN
647     END

  ViewVC Help
Powered by ViewVC 1.1.22