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

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

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

Revision 1.8 - (show annotations) (download)
Sun Jun 28 01:05:41 2009 UTC (14 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62y, checkpoint62x, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61s, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.7: +2 -2 lines 
add bj in exch2 arrays and S/R
1 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_som_advect.F,v 1.7 2009/06/26 23:10:09 jahn Exp $
2 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, deltaTLev,
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 #ifdef ALLOW_EXCH2
34 #include "W2_EXCH2_SIZE.h"
35 #include "W2_EXCH2_TOPOLOGY.h"
36 #endif /* ALLOW_EXCH2 */
37
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 deltaTLev(Nr)
55 _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 _RL smCorners(OLx,OLy,4,-1:nSOM)
137 c _RL localTr
138 LOGICAL calc_fluxes_X, calc_fluxes_Y
139 LOGICAL interiorOnly, overlapOnly
140 INTEGER limiter
141 INTEGER npass, ipass
142 INTEGER nCFace, n
143 LOGICAL N_edge, S_edge, E_edge, W_edge
144 CHARACTER*(MAX_LEN_MBUF) msgBuf
145 #ifdef ALLOW_EXCH2
146 INTEGER myTile
147 #endif
148 #ifdef ALLOW_DIAGNOSTICS
149 CHARACTER*8 diagName
150 CHARACTER*4 diagSufx
151 LOGICAL doDiagAdvX, doDiagAdvY, doDiagAdvR
152 C- Functions:
153 CHARACTER*4 GAD_DIAG_SUFX
154 EXTERNAL GAD_DIAG_SUFX
155 LOGICAL DIAGNOSTICS_IS_ON
156 EXTERNAL DIAGNOSTICS_IS_ON
157 #endif
158 CEOP
159
160 #ifdef ALLOW_DIAGNOSTICS
161 C-- Set diagnostics flags and suffix for the current tracer
162 doDiagAdvX = .FALSE.
163 doDiagAdvY = .FALSE.
164 doDiagAdvR = .FALSE.
165 IF ( useDiagnostics ) THEN
166 diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
167 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 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 DO j=1-OLy,sNy+OLy
182 DO i=1-OLx,sNx+OLx
183 afx(i,j) = 0.
184 afy(i,j) = 0.
185 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 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
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 C-- Set tile-specific parameters for horizontal fluxes
219 IF (useCubedSphereExchange) THEN
220 npass = 3
221 #ifdef ALLOW_EXCH2
222 myTile = W2_myTileList(bi,bj)
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 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 CALL CALC_COMMON_FACTORS (
251 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 IF (useGMRedi)
258 & 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 DO j=1-OLy,sNy+OLy
265 DO i=1-OLx,sNx+OLx
266 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 ENDDO
274 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 DO ipass=1,npass
279
280 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 C- not CubedSphere
301 calc_fluxes_X = MOD(ipass,2).EQ.1
302 calc_fluxes_Y = .NOT.calc_fluxes_X
303 ENDIF
304
305 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
306
307 C-- X direction
308 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 I bi,bj,k, limiter,
329 I overlapOnly, interiorOnly,
330 I N_edge, S_edge, E_edge, W_edge,
331 I deltaTLev(k), uTrans,
332 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 ELSE
345 STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
346 ENDIF
347
348 #ifdef ALLOW_OBCS
349 C- Apply open boundary conditions
350 c IF ( useOBCS ) THEN
351 ccc localTij(i,j) = smTr0(i,j)/smVol(i,j)
352 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 #ifdef ALLOW_PTRACERS
357 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 #endif /* ALLOW_PTRACERS */
361 c ENDIF
362 ccc smTr0(i,j) = localTij(i,j)*smVol(i,j)
363 c ENDIF
364 #endif /* ALLOW_OBCS */
365
366 C-- End of X direction
367 ENDIF
368
369 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
370
371 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 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 I bi,bj,k, limiter,
393 I overlapOnly, interiorOnly,
394 I N_edge, S_edge, E_edge, W_edge,
395 I deltaTLev(k), vTrans,
396 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 ELSE
409 STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
410 ENDIF
411
412 #ifdef ALLOW_OBCS
413 C- Apply open boundary conditions
414 c IF (useOBCS) THEN
415 ccc localTij(i,j) = smTr0(i,j)/smVol(i,j)
416 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 #ifdef ALLOW_PTRACERS
421 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 #endif /* ALLOW_PTRACERS */
425 c ENDIF
426 ccc smTr0(i,j) = localTij(i,j)*smVol(i,j)
427 c ENDIF
428 #endif /* ALLOW_OBCS */
429
430 C-- End of Y direction
431 ENDIF
432
433 C-- End of ipass loop
434 ENDDO
435
436 IF ( implicitAdvection ) THEN
437 C- explicit advection is done ; store tendency in gTracer:
438 DO j=1-OLy,sNy+OLy
439 DO i=1-OLx,sNx+OLx
440 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 c & / deltaTLev(k)
444 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 & /deltaTLev(k)
451 ENDDO
452 ENDDO
453 ENDIF
454
455 #ifdef ALLOW_DIAGNOSTICS
456 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 #endif
465
466 #ifdef ALLOW_DEBUG
467 IF ( debugLevel .GE. debLevB
468 & .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 ENDIF
475 #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 IF (useGMRedi .AND. k.GT.1 )
585 & 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 & .OR. vertAdvecScheme.EQ.ENUM_SOM_LIMITER ) THEN
593 CALL GAD_SOM_ADV_R(
594 I bi,bj,k, kUp, kDown,
595 I deltaTLev(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 STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM'
614 ENDIF
615
616 C-- Compute new tracer value and store tracer tendency
617 DO j=1-OLy,sNy+OLy
618 DO i=1-OLx,sNx+OLx
619 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 c & / deltaTLev(k)
623 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 & *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 & /deltaTLev(k)
630 ENDDO
631 ENDDO
632
633 #ifdef ALLOW_DIAGNOSTICS
634 IF ( doDiagAdvR ) THEN
635 diagName = 'ADVr'//diagSufx
636 CALL DIAGNOSTICS_FILL( afr,
637 & diagName, k,1, 2,bi,bj, myThid )
638 ENDIF
639 #endif
640
641 C-- End of k loop for vertical flux
642 ENDDO
643 C-- end of if not.implicitAdvection block
644 ENDIF
645
646 RETURN
647 END
  ViewVC Help
Powered by ViewVC 1.1.22