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

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

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

Revision 1.43 - (show annotations) (download)
Tue Dec 5 05:26:46 2006 UTC (17 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint58t_post, mitgcm_mapl_00, checkpoint58s_post
Changes since 1.42: +37 -29 lines 
start to implement deep-atmosphere and/or anelastic formulation
1 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.42 2006/07/18 03:23:50 heimbach Exp $
2 C $Name: $
3
4 #include "GAD_OPTIONS.h"
5 #undef MULTIDIM_OLD_VERSION
6
7 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
8 CBOP
9 C !ROUTINE: GAD_ADVECTION
10
11 C !INTERFACE: ==========================================================
12 SUBROUTINE GAD_ADVECTION(
13 I implicitAdvection, advectionScheme, vertAdvecScheme,
14 I tracerIdentity,
15 I uVel, vVel, wVel, tracer,
16 O gTracer,
17 I bi,bj, myTime,myIter,myThid)
18
19 C !DESCRIPTION:
20 C Calculates the tendancy of a tracer due to advection.
21 C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
22 C and can only be used for the non-linear advection schemes such as the
23 C direct-space-time method and flux-limiters.
24 C
25 C The algorithm is as follows:
26 C \begin{itemize}
27 C \item{$\theta^{(n+1/3)} = \theta^{(n)}
28 C - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$}
29 C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)}
30 C - \Delta t \partial_y (v\theta^{(n+1/3)}) + \theta^{(n)} \partial_y v$}
31 C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)}
32 C - \Delta t \partial_r (w\theta^{(n+2/3)}) + \theta^{(n)} \partial_r w$}
33 C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}
34 C \end{itemize}
35 C
36 C The tendancy (output) is over-written by this routine.
37
38 C !USES: ===============================================================
39 IMPLICIT NONE
40 #include "SIZE.h"
41 #include "EEPARAMS.h"
42 #include "PARAMS.h"
43 #include "GRID.h"
44 #include "GAD.h"
45 #ifdef ALLOW_AUTODIFF_TAMC
46 # include "tamc.h"
47 # include "tamc_keys.h"
48 # ifdef ALLOW_PTRACERS
49 # include "PTRACERS_SIZE.h"
50 # endif
51 #endif
52 #ifdef ALLOW_EXCH2
53 #include "W2_EXCH2_TOPOLOGY.h"
54 #include "W2_EXCH2_PARAMS.h"
55 #endif /* ALLOW_EXCH2 */
56
57 C !INPUT PARAMETERS: ===================================================
58 C implicitAdvection :: implicit vertical advection (later on)
59 C advectionScheme :: advection scheme to use (Horizontal plane)
60 C vertAdvecScheme :: advection scheme to use (vertical direction)
61 C tracerIdentity :: tracer identifier (required only for OBCS)
62 C uVel :: velocity, zonal component
63 C vVel :: velocity, meridional component
64 C wVel :: velocity, vertical component
65 C tracer :: tracer field
66 C bi,bj :: tile indices
67 C myTime :: current time
68 C myIter :: iteration number
69 C myThid :: thread number
70 LOGICAL implicitAdvection
71 INTEGER advectionScheme, vertAdvecScheme
72 INTEGER tracerIdentity
73 _RL uVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
74 _RL vVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
75 _RL wVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
76 _RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
77 INTEGER bi,bj
78 _RL myTime
79 INTEGER myIter
80 INTEGER myThid
81
82 C !OUTPUT PARAMETERS: ==================================================
83 C gTracer :: tendancy array
84 _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
85
86 C !LOCAL VARIABLES: ====================================================
87 C maskUp :: 2-D array for mask at W points
88 C maskLocW :: 2-D array for mask at West points
89 C maskLocS :: 2-D array for mask at South points
90 C iMin,iMax, :: loop range for called routines
91 C jMin,jMax :: loop range for called routines
92 C [iMin,iMax]Upd :: loop range to update tracer field
93 C [jMin,jMax]Upd :: loop range to update tracer field
94 C i,j,k :: loop indices
95 C kUp :: index into 2 1/2D array, toggles between 1 and 2
96 C kDown :: index into 2 1/2D array, toggles between 2 and 1
97 C kp1 :: =k+1 for k<Nr, =Nr for k=Nr
98 C xA,yA :: areas of X and Y face of tracer cells
99 C uFld,vFld :: 2-D local copy of horizontal velocity, U,V components
100 C wFld :: 2-D local copy of vertical velocity
101 C uTrans,vTrans :: 2-D arrays of volume transports at U,V points
102 C rTrans :: 2-D arrays of volume transports at W points
103 C rTransKp1 :: vertical volume transport at interface k+1
104 C af :: 2-D array for horizontal advective flux
105 C afx :: 2-D array for horizontal advective flux, x direction
106 C afy :: 2-D array for horizontal advective flux, y direction
107 C fVerT :: 2 1/2D arrays for vertical advective flux
108 C localTij :: 2-D array, temporary local copy of tracer fld
109 C localTijk :: 3-D array, temporary local copy of tracer fld
110 C kp1Msk :: flag (0,1) for over-riding mask for W levels
111 C calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
112 C calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
113 C interiorOnly :: only update the interior of myTile, but not the edges
114 C overlapOnly :: only update the edges of myTile, but not the interior
115 C nipass :: number of passes in multi-dimensional method
116 C ipass :: number of the current pass being made
117 C myTile :: variables used to determine which cube face
118 C nCFace :: owns a tile for cube grid runs using
119 C :: multi-dim advection.
120 C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
121 c _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122 _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123 _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124 INTEGER iMin,iMax,jMin,jMax
125 INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
126 INTEGER i,j,k,kUp,kDown
127 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129 _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130 _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131 _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135 _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136 _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137 _RL afx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138 _RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
140 _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141 _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
142 _RL kp1Msk
143 LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
144 LOGICAL interiorOnly, overlapOnly
145 INTEGER nipass,ipass
146 INTEGER nCFace
147 LOGICAL N_edge, S_edge, E_edge, W_edge
148 #ifdef ALLOW_EXCH2
149 INTEGER myTile
150 #endif
151 #ifdef ALLOW_DIAGNOSTICS
152 CHARACTER*8 diagName
153 CHARACTER*4 GAD_DIAG_SUFX, diagSufx
154 EXTERNAL GAD_DIAG_SUFX
155 #endif
156 CEOP
157
158 #ifdef ALLOW_AUTODIFF_TAMC
159 act0 = tracerIdentity - 1
160 max0 = maxpass
161 act1 = bi - myBxLo(myThid)
162 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
163 act2 = bj - myByLo(myThid)
164 max2 = myByHi(myThid) - myByLo(myThid) + 1
165 act3 = myThid - 1
166 max3 = nTx*nTy
167 act4 = ikey_dynamics - 1
168 igadkey = (act0 + 1)
169 & + act1*max0
170 & + act2*max0*max1
171 & + act3*max0*max1*max2
172 & + act4*max0*max1*max2*max3
173 if (tracerIdentity.GT.maxpass) then
174 print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
175 STOP 'maxpass seems smaller than tracerIdentity'
176 endif
177 #endif /* ALLOW_AUTODIFF_TAMC */
178
179 #ifdef ALLOW_DIAGNOSTICS
180 C-- Set diagnostic suffix for the current tracer
181 IF ( useDiagnostics ) THEN
182 diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
183 ENDIF
184 #endif
185
186 C-- Set up work arrays with valid (i.e. not NaN) values
187 C These inital values do not alter the numerical results. They
188 C just ensure that all memory references are to valid floating
189 C point numbers. This prevents spurious hardware signals due to
190 C uninitialised but inert locations.
191 DO j=1-OLy,sNy+OLy
192 DO i=1-OLx,sNx+OLx
193 xA(i,j) = 0. _d 0
194 yA(i,j) = 0. _d 0
195 uTrans(i,j) = 0. _d 0
196 vTrans(i,j) = 0. _d 0
197 rTrans(i,j) = 0. _d 0
198 fVerT(i,j,1) = 0. _d 0
199 fVerT(i,j,2) = 0. _d 0
200 rTransKp1(i,j)= 0. _d 0
201 #ifdef ALLOW_AUTODIFF_TAMC
202 localTij(i,j) = 0. _d 0
203 wfld(i,j) = 0. _d 0
204 #endif
205 ENDDO
206 ENDDO
207
208 C-- Set tile-specific parameters for horizontal fluxes
209 IF (useCubedSphereExchange) THEN
210 nipass=3
211 #ifdef ALLOW_AUTODIFF_TAMC
212 IF ( nipass.GT.maxcube ) STOP 'maxcube needs to be = 3'
213 #endif
214 #ifdef ALLOW_EXCH2
215 myTile = W2_myTileList(bi)
216 nCFace = exch2_myFace(myTile)
217 N_edge = exch2_isNedge(myTile).EQ.1
218 S_edge = exch2_isSedge(myTile).EQ.1
219 E_edge = exch2_isEedge(myTile).EQ.1
220 W_edge = exch2_isWedge(myTile).EQ.1
221 #else
222 nCFace = bi
223 N_edge = .TRUE.
224 S_edge = .TRUE.
225 E_edge = .TRUE.
226 W_edge = .TRUE.
227 #endif
228 ELSE
229 nipass=2
230 nCFace = bi
231 N_edge = .FALSE.
232 S_edge = .FALSE.
233 E_edge = .FALSE.
234 W_edge = .FALSE.
235 ENDIF
236
237 iMin = 1-OLx
238 iMax = sNx+OLx
239 jMin = 1-OLy
240 jMax = sNy+OLy
241
242 C-- Start of k loop for horizontal fluxes
243 DO k=1,Nr
244 #ifdef ALLOW_AUTODIFF_TAMC
245 kkey = (igadkey-1)*Nr + k
246 CADJ STORE tracer(:,:,k,bi,bj) =
247 CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
248 #endif /* ALLOW_AUTODIFF_TAMC */
249
250 C-- Get temporary terms used by tendency routines
251 CALL CALC_COMMON_FACTORS (
252 I uVel, vVel,
253 O uFld, vFld, uTrans, vTrans, xA, yA,
254 I k,bi,bj, myThid )
255
256 #ifdef ALLOW_GMREDI
257 C-- Residual transp = Bolus transp + Eulerian transp
258 IF (useGMRedi)
259 & CALL GMREDI_CALC_UVFLOW(
260 U uFld, vFld, uTrans, vTrans,
261 I k, bi, bj, myThid )
262 #endif /* ALLOW_GMREDI */
263
264 C-- Make local copy of tracer array and mask West & South
265 DO j=1-OLy,sNy+OLy
266 DO i=1-OLx,sNx+OLx
267 localTij(i,j)=tracer(i,j,k,bi,bj)
268 maskLocW(i,j)=maskW(i,j,k,bi,bj)
269 maskLocS(i,j)=maskS(i,j,k,bi,bj)
270 ENDDO
271 ENDDO
272
273 #ifndef ALLOW_AUTODIFF_TAMC
274 IF (useCubedSphereExchange) THEN
275 withSigns = .FALSE.
276 CALL FILL_CS_CORNER_UV_RS(
277 & withSigns, maskLocW,maskLocS, bi,bj, myThid )
278 ENDIF
279 #endif
280
281 C-- Multiple passes for different directions on different tiles
282 C-- For cube need one pass for each of red, green and blue axes.
283 DO ipass=1,nipass
284 #ifdef ALLOW_AUTODIFF_TAMC
285 passkey = ipass + (k-1) *maxcube
286 & + (igadkey-1)*maxcube*Nr
287 IF (nipass .GT. maxpass) THEN
288 STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
289 ENDIF
290 #endif /* ALLOW_AUTODIFF_TAMC */
291
292 interiorOnly = .FALSE.
293 overlapOnly = .FALSE.
294 IF (useCubedSphereExchange) THEN
295 #ifdef MULTIDIM_OLD_VERSION
296 C- CubedSphere : pass 3 times, with full update of local tracer field
297 IF (ipass.EQ.1) THEN
298 calc_fluxes_X = nCFace.EQ.1 .OR. nCFace.EQ.2
299 calc_fluxes_Y = nCFace.EQ.4 .OR. nCFace.EQ.5
300 ELSEIF (ipass.EQ.2) THEN
301 calc_fluxes_X = nCFace.EQ.3 .OR. nCFace.EQ.4
302 calc_fluxes_Y = nCFace.EQ.6 .OR. nCFace.EQ.1
303 #else /* MULTIDIM_OLD_VERSION */
304 C- CubedSphere : pass 3 times, with partial update of local tracer field
305 IF (ipass.EQ.1) THEN
306 overlapOnly = MOD(nCFace,3).EQ.0
307 interiorOnly = MOD(nCFace,3).NE.0
308 calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
309 calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
310 ELSEIF (ipass.EQ.2) THEN
311 overlapOnly = MOD(nCFace,3).EQ.2
312 calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
313 calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
314 #endif /* MULTIDIM_OLD_VERSION */
315 ELSE
316 calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
317 calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
318 ENDIF
319 ELSE
320 C- not CubedSphere
321 calc_fluxes_X = MOD(ipass,2).EQ.1
322 calc_fluxes_Y = .NOT.calc_fluxes_X
323 ENDIF
324
325 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
326 C-- X direction
327 C- Advective flux in X
328 DO j=1-Oly,sNy+Oly
329 DO i=1-Olx,sNx+Olx
330 af(i,j) = 0.
331 ENDDO
332 ENDDO
333 C
334 #ifdef ALLOW_AUTODIFF_TAMC
335 # ifndef DISABLE_MULTIDIM_ADVECTION
336 CADJ STORE localTij(:,:) =
337 CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
338 CADJ STORE af(:,:) =
339 CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
340 # endif
341 #endif /* ALLOW_AUTODIFF_TAMC */
342 C
343 IF (calc_fluxes_X) THEN
344
345 C- Do not compute fluxes if
346 C a) needed in overlap only
347 C and b) the overlap of myTile are not cube-face Edges
348 IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
349
350 #ifndef ALLOW_AUTODIFF_TAMC
351 C- Internal exchange for calculations in X
352 #ifdef MULTIDIM_OLD_VERSION
353 IF ( useCubedSphereExchange ) THEN
354 #else
355 IF ( useCubedSphereExchange .AND.
356 & ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
357 #endif
358 CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )
359 ENDIF
360 #endif
361
362 #ifdef ALLOW_AUTODIFF_TAMC
363 # ifndef DISABLE_MULTIDIM_ADVECTION
364 CADJ STORE localTij(:,:) =
365 CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
366 # endif
367 #endif /* ALLOW_AUTODIFF_TAMC */
368
369 IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
370 & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
371 CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,
372 I dTtracerLev(k),uTrans,uFld,localTij,
373 O af, myThid )
374 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
375 CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
376 I uTrans, uFld, maskLocW, localTij,
377 O af, myThid )
378 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
379 CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
380 I uTrans, uFld, maskLocW, localTij,
381 O af, myThid )
382 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
383 CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),
384 I uTrans, uFld, maskLocW, localTij,
385 O af, myThid )
386 ELSE
387 STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
388 ENDIF
389
390 C- Advective flux in X : done
391 ENDIF
392
393 #ifndef ALLOW_AUTODIFF_TAMC
394 C- Internal exchange for next calculations in Y
395 IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
396 CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )
397 ENDIF
398 #endif
399
400 C- Update the local tracer field where needed:
401
402 C update in overlap-Only
403 IF ( overlapOnly ) THEN
404 iMinUpd = 1-Olx+1
405 iMaxUpd = sNx+Olx-1
406 C- notes: these 2 lines below have no real effect (because recip_hFac=0
407 C in corner region) but safer to keep them.
408 IF ( W_edge ) iMinUpd = 1
409 IF ( E_edge ) iMaxUpd = sNx
410
411 IF ( S_edge ) THEN
412 DO j=1-Oly,0
413 DO i=iMinUpd,iMaxUpd
414 localTij(i,j) = localTij(i,j)
415 & -dTtracerLev(k)*recip_rhoFacC(k)
416 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
417 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
418 & *( af(i+1,j)-af(i,j)
419 & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
420 & )
421 ENDDO
422 ENDDO
423 ENDIF
424 IF ( N_edge ) THEN
425 DO j=sNy+1,sNy+Oly
426 DO i=iMinUpd,iMaxUpd
427 localTij(i,j) = localTij(i,j)
428 & -dTtracerLev(k)*recip_rhoFacC(k)
429 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
430 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
431 & *( af(i+1,j)-af(i,j)
432 & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
433 & )
434 ENDDO
435 ENDDO
436 ENDIF
437
438 ELSE
439 C do not only update the overlap
440 jMinUpd = 1-Oly
441 jMaxUpd = sNy+Oly
442 IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
443 IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
444 DO j=jMinUpd,jMaxUpd
445 DO i=1-Olx+1,sNx+Olx-1
446 localTij(i,j) = localTij(i,j)
447 & -dTtracerLev(k)*recip_rhoFacC(k)
448 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
449 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
450 & *( af(i+1,j)-af(i,j)
451 & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
452 & )
453 ENDDO
454 ENDDO
455 C- keep advective flux (for diagnostics)
456 DO j=1-Oly,sNy+Oly
457 DO i=1-Olx,sNx+Olx
458 afx(i,j) = af(i,j)
459 ENDDO
460 ENDDO
461
462 #ifdef ALLOW_OBCS
463 C- Apply open boundary conditions
464 IF ( useOBCS ) THEN
465 IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
466 CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
467 ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
468 CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
469 #ifdef ALLOW_PTRACERS
470 ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
471 CALL OBCS_APPLY_PTRACER( bi, bj, k,
472 & tracerIdentity-GAD_TR1+1, localTij, myThid )
473 #endif /* ALLOW_PTRACERS */
474 ENDIF
475 ENDIF
476 #endif /* ALLOW_OBCS */
477
478 C- end if/else update overlap-Only
479 ENDIF
480
481 C-- End of X direction
482 ENDIF
483
484 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
485 C-- Y direction
486 cph-test
487 C- Advective flux in Y
488 DO j=1-Oly,sNy+Oly
489 DO i=1-Olx,sNx+Olx
490 af(i,j) = 0.
491 ENDDO
492 ENDDO
493 C
494 #ifdef ALLOW_AUTODIFF_TAMC
495 # ifndef DISABLE_MULTIDIM_ADVECTION
496 CADJ STORE localTij(:,:) =
497 CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
498 CADJ STORE af(:,:) =
499 CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
500 # endif
501 #endif /* ALLOW_AUTODIFF_TAMC */
502 C
503 IF (calc_fluxes_Y) THEN
504
505 C- Do not compute fluxes if
506 C a) needed in overlap only
507 C and b) the overlap of myTile are not cube-face edges
508 IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
509
510 #ifndef ALLOW_AUTODIFF_TAMC
511 C- Internal exchange for calculations in Y
512 #ifdef MULTIDIM_OLD_VERSION
513 IF ( useCubedSphereExchange ) THEN
514 #else
515 IF ( useCubedSphereExchange .AND.
516 & ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
517 #endif
518 CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )
519 ENDIF
520 #endif
521
522 C- Advective flux in Y
523 DO j=1-Oly,sNy+Oly
524 DO i=1-Olx,sNx+Olx
525 af(i,j) = 0.
526 ENDDO
527 ENDDO
528
529 #ifdef ALLOW_AUTODIFF_TAMC
530 #ifndef DISABLE_MULTIDIM_ADVECTION
531 CADJ STORE localTij(:,:) =
532 CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
533 #endif
534 #endif /* ALLOW_AUTODIFF_TAMC */
535
536 IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
537 & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
538 CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,
539 I dTtracerLev(k),vTrans,vFld,localTij,
540 O af, myThid )
541 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
542 CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
543 I vTrans, vFld, maskLocS, localTij,
544 O af, myThid )
545 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
546 CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
547 I vTrans, vFld, maskLocS, localTij,
548 O af, myThid )
549 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
550 CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),
551 I vTrans, vFld, maskLocS, localTij,
552 O af, myThid )
553 ELSE
554 STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
555 ENDIF
556
557 C- Advective flux in Y : done
558 ENDIF
559
560 #ifndef ALLOW_AUTODIFF_TAMC
561 C- Internal exchange for next calculations in X
562 IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
563 CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )
564 ENDIF
565 #endif
566
567 C- Update the local tracer field where needed:
568
569 C update in overlap-Only
570 IF ( overlapOnly ) THEN
571 jMinUpd = 1-Oly+1
572 jMaxUpd = sNy+Oly-1
573 C- notes: these 2 lines below have no real effect (because recip_hFac=0
574 C in corner region) but safer to keep them.
575 IF ( S_edge ) jMinUpd = 1
576 IF ( N_edge ) jMaxUpd = sNy
577
578 IF ( W_edge ) THEN
579 DO j=jMinUpd,jMaxUpd
580 DO i=1-Olx,0
581 localTij(i,j) = localTij(i,j)
582 & -dTtracerLev(k)*recip_rhoFacC(k)
583 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
584 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
585 & *( af(i,j+1)-af(i,j)
586 & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
587 & )
588 ENDDO
589 ENDDO
590 ENDIF
591 IF ( E_edge ) THEN
592 DO j=jMinUpd,jMaxUpd
593 DO i=sNx+1,sNx+Olx
594 localTij(i,j) = localTij(i,j)
595 & -dTtracerLev(k)*recip_rhoFacC(k)
596 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
597 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
598 & *( af(i,j+1)-af(i,j)
599 & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
600 & )
601 ENDDO
602 ENDDO
603 ENDIF
604
605 ELSE
606 C do not only update the overlap
607 iMinUpd = 1-Olx
608 iMaxUpd = sNx+Olx
609 IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
610 IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
611 DO j=1-Oly+1,sNy+Oly-1
612 DO i=iMinUpd,iMaxUpd
613 localTij(i,j) = localTij(i,j)
614 & -dTtracerLev(k)*recip_rhoFacC(k)
615 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
616 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
617 & *( af(i,j+1)-af(i,j)
618 & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
619 & )
620 ENDDO
621 ENDDO
622 C- keep advective flux (for diagnostics)
623 DO j=1-Oly,sNy+Oly
624 DO i=1-Olx,sNx+Olx
625 afy(i,j) = af(i,j)
626 ENDDO
627 ENDDO
628
629 #ifdef ALLOW_OBCS
630 C- Apply open boundary conditions
631 IF (useOBCS) THEN
632 IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
633 CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
634 ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
635 CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
636 #ifdef ALLOW_PTRACERS
637 ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
638 CALL OBCS_APPLY_PTRACER( bi, bj, k,
639 & tracerIdentity-GAD_TR1+1, localTij, myThid )
640 #endif /* ALLOW_PTRACERS */
641 ENDIF
642 ENDIF
643 #endif /* ALLOW_OBCS */
644
645 C end if/else update overlap-Only
646 ENDIF
647
648 C-- End of Y direction
649 ENDIF
650
651 C-- End of ipass loop
652 ENDDO
653
654 IF ( implicitAdvection ) THEN
655 C- explicit advection is done ; store tendency in gTracer:
656 DO j=1-Oly,sNy+Oly
657 DO i=1-Olx,sNx+Olx
658 gTracer(i,j,k,bi,bj)=
659 & (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)
660 ENDDO
661 ENDDO
662 ELSE
663 C- horizontal advection done; store intermediate result in 3D array:
664 DO j=1-Oly,sNy+Oly
665 DO i=1-Olx,sNx+Olx
666 localTijk(i,j,k)=localTij(i,j)
667 ENDDO
668 ENDDO
669 ENDIF
670
671 #ifdef ALLOW_DIAGNOSTICS
672 IF ( useDiagnostics ) THEN
673 diagName = 'ADVx'//diagSufx
674 CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)
675 diagName = 'ADVy'//diagSufx
676 CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)
677 ENDIF
678 #endif
679
680 #ifdef ALLOW_DEBUG
681 IF ( debugLevel .GE. debLevB
682 & .AND. tracerIdentity.EQ.GAD_TEMPERATURE
683 & .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
684 & .AND. nPx.EQ.1 .AND. nPy.EQ.1
685 & .AND. useCubedSphereExchange ) THEN
686 CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION',
687 & afx,afy, k, standardMessageUnit,bi,bj,myThid )
688 ENDIF
689 #endif /* ALLOW_DEBUG */
690
691 C-- End of K loop for horizontal fluxes
692 ENDDO
693
694 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
695
696 IF ( .NOT.implicitAdvection ) THEN
697 C-- Start of k loop for vertical flux
698 DO k=Nr,1,-1
699 #ifdef ALLOW_AUTODIFF_TAMC
700 kkey = (igadkey-1)*Nr + k
701 #endif /* ALLOW_AUTODIFF_TAMC */
702 C-- kUp Cycles through 1,2 to point to w-layer above
703 C-- kDown Cycles through 2,1 to point to w-layer below
704 kUp = 1+MOD(k+1,2)
705 kDown= 1+MOD(k,2)
706 c kp1=min(Nr,k+1)
707 kp1Msk=1.
708 if (k.EQ.Nr) kp1Msk=0.
709
710 C-- Compute Vertical transport
711 #ifdef ALLOW_AIM
712 C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
713 IF ( k.EQ.1 .OR.
714 & (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
715 & ) THEN
716 #else
717 IF ( k.EQ.1 ) THEN
718 #endif
719
720 #ifdef ALLOW_AUTODIFF_TAMC
721 CADJ STORE rtrans(:,:) =
722 CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
723 CADJ STORE wfld(:,:) =
724 CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
725 #endif /* ALLOW_AUTODIFF_TAMC */
726
727 C- Surface interface :
728 DO j=1-Oly,sNy+Oly
729 DO i=1-Olx,sNx+Olx
730 rTransKp1(i,j) = kp1Msk*rTrans(i,j)
731 wFld(i,j) = 0.
732 rTrans(i,j) = 0.
733 fVerT(i,j,kUp) = 0.
734 ENDDO
735 ENDDO
736
737 ELSE
738
739 #ifdef ALLOW_AUTODIFF_TAMC
740 CADJ STORE rtrans(:,:) =
741 CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
742 CADJ STORE wfld(:,:) =
743 CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
744 #endif /* ALLOW_AUTODIFF_TAMC */
745
746 C- Interior interface :
747 DO j=1-Oly,sNy+Oly
748 DO i=1-Olx,sNx+Olx
749 rTransKp1(i,j) = kp1Msk*rTrans(i,j)
750 wFld(i,j) = wVel(i,j,k,bi,bj)
751 rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
752 & *deepFac2F(k)*rhoFacF(k)
753 & *maskC(i,j,k-1,bi,bj)
754 fVerT(i,j,kUp) = 0.
755 ENDDO
756 ENDDO
757
758 #ifdef ALLOW_GMREDI
759 C-- Residual transp = Bolus transp + Eulerian transp
760 IF (useGMRedi)
761 & CALL GMREDI_CALC_WFLOW(
762 U wFld, rTrans,
763 I k, bi, bj, myThid )
764 #endif /* ALLOW_GMREDI */
765
766 #ifdef ALLOW_AUTODIFF_TAMC
767 CADJ STORE localTijk(:,:,k)
768 CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
769 CADJ STORE rTrans(:,:)
770 CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
771 #endif /* ALLOW_AUTODIFF_TAMC */
772
773 C- Compute vertical advective flux in the interior:
774 IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
775 & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
776 CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,
777 I dTtracerLev(k),rTrans,wFld,localTijk,
778 O fVerT(1-Olx,1-Oly,kUp), myThid )
779 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
780 CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k),
781 I rTrans, wFld, localTijk,
782 O fVerT(1-Olx,1-Oly,kUp), myThid )
783 ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
784 CALL GAD_DST3_ADV_R( bi,bj,k, dTtracerLev(k),
785 I rTrans, wFld, localTijk,
786 O fVerT(1-Olx,1-Oly,kUp), myThid )
787 ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
788 CALL GAD_DST3FL_ADV_R( bi,bj,k, dTtracerLev(k),
789 I rTrans, wFld, localTijk,
790 O fVerT(1-Olx,1-Oly,kUp), myThid )
791 ELSE
792 STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
793 ENDIF
794
795 C- end Surface/Interior if bloc
796 ENDIF
797
798 #ifdef ALLOW_AUTODIFF_TAMC
799 CADJ STORE rTrans(:,:)
800 CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
801 CADJ STORE rTranskp1(:,:)
802 CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
803 #endif /* ALLOW_AUTODIFF_TAMC */
804
805 C-- Divergence of vertical fluxes
806 DO j=1-Oly,sNy+Oly
807 DO i=1-Olx,sNx+Olx
808 localTij(i,j) = localTijk(i,j,k)
809 & -dTtracerLev(k)*recip_rhoFacC(k)
810 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
811 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
812 & *( fVerT(i,j,kDown)-fVerT(i,j,kUp)
813 & -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))
814 & )*rkSign
815 gTracer(i,j,k,bi,bj)=
816 & (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)
817 ENDDO
818 ENDDO
819
820 #ifdef ALLOW_DIAGNOSTICS
821 IF ( useDiagnostics ) THEN
822 diagName = 'ADVr'//diagSufx
823 CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp),
824 & diagName, k,1, 2,bi,bj, myThid)
825 ENDIF
826 #endif
827
828 C-- End of K loop for vertical flux
829 ENDDO
830 C-- end of if not.implicitAdvection block
831 ENDIF
832
833 RETURN
834 END
  ViewVC Help
Powered by ViewVC 1.1.22