/[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.49 - (show annotations) (download)
Fri Jul 27 22:18:57 2007 UTC (16 years, 10 months ago) by heimbach
Branch: MAIN
Changes since 1.48: +11 -11 lines 
Comment all relevant #ifndef ALLOW_AUTODIFF_TAMC
that used to hide exch2 or cubed-sphere specific code
(commented via 'cph-exch2')
1 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.48 2007/04/04 01:39:06 jmc 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 tendency 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 tendency (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 :: tendency 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 cph-exch2#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 cph-exch2#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 cph-exch2#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 cph-exch2#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, .TRUE.,
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, .TRUE., 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, .TRUE., 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, .TRUE., dTtracerLev(k),
384 I uTrans, uFld, maskLocW, localTij,
385 O af, myThid )
386 #ifndef ALLOW_AUTODIFF_TAMC
387 ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
388 CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
389 I uTrans, uFld, maskLocW, localTij,
390 O af, myThid )
391 #endif
392 ELSE
393 STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
394 ENDIF
395
396 C- Advective flux in X : done
397 ENDIF
398
399 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
400 C- Internal exchange for next calculations in Y
401 IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
402 CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )
403 ENDIF
404 cph-exch2#endif
405
406 C- Update the local tracer field where needed:
407
408 C update in overlap-Only
409 IF ( overlapOnly ) THEN
410 iMinUpd = 1-Olx+1
411 iMaxUpd = sNx+Olx-1
412 C- notes: these 2 lines below have no real effect (because recip_hFac=0
413 C in corner region) but safer to keep them.
414 IF ( W_edge ) iMinUpd = 1
415 IF ( E_edge ) iMaxUpd = sNx
416
417 IF ( S_edge ) THEN
418 DO j=1-Oly,0
419 DO i=iMinUpd,iMaxUpd
420 localTij(i,j) = localTij(i,j)
421 & -dTtracerLev(k)*recip_rhoFacC(k)
422 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
423 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
424 & *( af(i+1,j)-af(i,j)
425 & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
426 & )
427 ENDDO
428 ENDDO
429 ENDIF
430 IF ( N_edge ) THEN
431 DO j=sNy+1,sNy+Oly
432 DO i=iMinUpd,iMaxUpd
433 localTij(i,j) = localTij(i,j)
434 & -dTtracerLev(k)*recip_rhoFacC(k)
435 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
436 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
437 & *( af(i+1,j)-af(i,j)
438 & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
439 & )
440 ENDDO
441 ENDDO
442 ENDIF
443
444 ELSE
445 C do not only update the overlap
446 jMinUpd = 1-Oly
447 jMaxUpd = sNy+Oly
448 IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
449 IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
450 DO j=jMinUpd,jMaxUpd
451 DO i=1-Olx+1,sNx+Olx-1
452 localTij(i,j) = localTij(i,j)
453 & -dTtracerLev(k)*recip_rhoFacC(k)
454 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
455 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
456 & *( af(i+1,j)-af(i,j)
457 & -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
458 & )
459 ENDDO
460 ENDDO
461 C- keep advective flux (for diagnostics)
462 DO j=1-Oly,sNy+Oly
463 DO i=1-Olx,sNx+Olx
464 afx(i,j) = af(i,j)
465 ENDDO
466 ENDDO
467
468 #ifdef ALLOW_OBCS
469 C- Apply open boundary conditions
470 IF ( useOBCS ) THEN
471 IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
472 CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
473 ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
474 CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
475 #ifdef ALLOW_PTRACERS
476 ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
477 CALL OBCS_APPLY_PTRACER( bi, bj, k,
478 & tracerIdentity-GAD_TR1+1, localTij, myThid )
479 #endif /* ALLOW_PTRACERS */
480 ENDIF
481 ENDIF
482 #endif /* ALLOW_OBCS */
483
484 C- end if/else update overlap-Only
485 ENDIF
486
487 C-- End of X direction
488 ENDIF
489
490 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
491 C-- Y direction
492 cph-test
493 C- Advective flux in Y
494 DO j=1-Oly,sNy+Oly
495 DO i=1-Olx,sNx+Olx
496 af(i,j) = 0.
497 ENDDO
498 ENDDO
499 C
500 #ifdef ALLOW_AUTODIFF_TAMC
501 # ifndef DISABLE_MULTIDIM_ADVECTION
502 CADJ STORE localTij(:,:) =
503 CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
504 CADJ STORE af(:,:) =
505 CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
506 # endif
507 #endif /* ALLOW_AUTODIFF_TAMC */
508 C
509 IF (calc_fluxes_Y) THEN
510
511 C- Do not compute fluxes if
512 C a) needed in overlap only
513 C and b) the overlap of myTile are not cube-face edges
514 IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
515
516 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
517 C- Internal exchange for calculations in Y
518 #ifdef MULTIDIM_OLD_VERSION
519 IF ( useCubedSphereExchange ) THEN
520 #else
521 IF ( useCubedSphereExchange .AND.
522 & ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
523 #endif
524 CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid )
525 ENDIF
526 cph-exch2#endif
527
528 C- Advective flux in Y
529 DO j=1-Oly,sNy+Oly
530 DO i=1-Olx,sNx+Olx
531 af(i,j) = 0.
532 ENDDO
533 ENDDO
534
535 #ifdef ALLOW_AUTODIFF_TAMC
536 #ifndef DISABLE_MULTIDIM_ADVECTION
537 CADJ STORE localTij(:,:) =
538 CADJ & comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
539 #endif
540 #endif /* ALLOW_AUTODIFF_TAMC */
541
542 IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
543 & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
544 CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
545 I dTtracerLev(k),vTrans,vFld,localTij,
546 O af, myThid )
547 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
548 CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
549 I vTrans, vFld, maskLocS, localTij,
550 O af, myThid )
551 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
552 CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
553 I vTrans, vFld, maskLocS, localTij,
554 O af, myThid )
555 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
556 CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
557 I vTrans, vFld, maskLocS, localTij,
558 O af, myThid )
559 #ifndef ALLOW_AUTODIFF_TAMC
560 ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
561 CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
562 I vTrans, vFld, maskLocS, localTij,
563 O af, myThid )
564 #endif
565 ELSE
566 STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
567 ENDIF
568
569 C- Advective flux in Y : done
570 ENDIF
571
572 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
573 C- Internal exchange for next calculations in X
574 IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
575 CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid )
576 ENDIF
577 cph-exch2#endif
578
579 C- Update the local tracer field where needed:
580
581 C update in overlap-Only
582 IF ( overlapOnly ) THEN
583 jMinUpd = 1-Oly+1
584 jMaxUpd = sNy+Oly-1
585 C- notes: these 2 lines below have no real effect (because recip_hFac=0
586 C in corner region) but safer to keep them.
587 IF ( S_edge ) jMinUpd = 1
588 IF ( N_edge ) jMaxUpd = sNy
589
590 IF ( W_edge ) THEN
591 DO j=jMinUpd,jMaxUpd
592 DO i=1-Olx,0
593 localTij(i,j) = localTij(i,j)
594 & -dTtracerLev(k)*recip_rhoFacC(k)
595 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
596 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
597 & *( af(i,j+1)-af(i,j)
598 & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
599 & )
600 ENDDO
601 ENDDO
602 ENDIF
603 IF ( E_edge ) THEN
604 DO j=jMinUpd,jMaxUpd
605 DO i=sNx+1,sNx+Olx
606 localTij(i,j) = localTij(i,j)
607 & -dTtracerLev(k)*recip_rhoFacC(k)
608 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
609 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
610 & *( af(i,j+1)-af(i,j)
611 & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
612 & )
613 ENDDO
614 ENDDO
615 ENDIF
616
617 ELSE
618 C do not only update the overlap
619 iMinUpd = 1-Olx
620 iMaxUpd = sNx+Olx
621 IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
622 IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
623 DO j=1-Oly+1,sNy+Oly-1
624 DO i=iMinUpd,iMaxUpd
625 localTij(i,j) = localTij(i,j)
626 & -dTtracerLev(k)*recip_rhoFacC(k)
627 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
628 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
629 & *( af(i,j+1)-af(i,j)
630 & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
631 & )
632 ENDDO
633 ENDDO
634 C- keep advective flux (for diagnostics)
635 DO j=1-Oly,sNy+Oly
636 DO i=1-Olx,sNx+Olx
637 afy(i,j) = af(i,j)
638 ENDDO
639 ENDDO
640
641 #ifdef ALLOW_OBCS
642 C- Apply open boundary conditions
643 IF (useOBCS) THEN
644 IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
645 CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
646 ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
647 CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
648 #ifdef ALLOW_PTRACERS
649 ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
650 CALL OBCS_APPLY_PTRACER( bi, bj, k,
651 & tracerIdentity-GAD_TR1+1, localTij, myThid )
652 #endif /* ALLOW_PTRACERS */
653 ENDIF
654 ENDIF
655 #endif /* ALLOW_OBCS */
656
657 C end if/else update overlap-Only
658 ENDIF
659
660 C-- End of Y direction
661 ENDIF
662
663 C-- End of ipass loop
664 ENDDO
665
666 IF ( implicitAdvection ) THEN
667 C- explicit advection is done ; store tendency in gTracer:
668 DO j=1-Oly,sNy+Oly
669 DO i=1-Olx,sNx+Olx
670 gTracer(i,j,k,bi,bj)=
671 & (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)
672 ENDDO
673 ENDDO
674 ELSE
675 C- horizontal advection done; store intermediate result in 3D array:
676 DO j=1-Oly,sNy+Oly
677 DO i=1-Olx,sNx+Olx
678 localTijk(i,j,k)=localTij(i,j)
679 ENDDO
680 ENDDO
681 ENDIF
682
683 #ifdef ALLOW_DIAGNOSTICS
684 IF ( useDiagnostics ) THEN
685 diagName = 'ADVx'//diagSufx
686 CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)
687 diagName = 'ADVy'//diagSufx
688 CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)
689 ENDIF
690 #endif
691
692 #ifdef ALLOW_DEBUG
693 IF ( debugLevel .GE. debLevB
694 & .AND. tracerIdentity.EQ.GAD_TEMPERATURE
695 & .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
696 & .AND. nPx.EQ.1 .AND. nPy.EQ.1
697 & .AND. useCubedSphereExchange ) THEN
698 CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION',
699 & afx,afy, k, standardMessageUnit,bi,bj,myThid )
700 ENDIF
701 #endif /* ALLOW_DEBUG */
702
703 C-- End of K loop for horizontal fluxes
704 ENDDO
705
706 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
707
708 IF ( .NOT.implicitAdvection ) THEN
709 C-- Start of k loop for vertical flux
710 DO k=Nr,1,-1
711 #ifdef ALLOW_AUTODIFF_TAMC
712 kkey = (igadkey-1)*Nr + k
713 #endif /* ALLOW_AUTODIFF_TAMC */
714 C-- kUp Cycles through 1,2 to point to w-layer above
715 C-- kDown Cycles through 2,1 to point to w-layer below
716 kUp = 1+MOD(k+1,2)
717 kDown= 1+MOD(k,2)
718 c kp1=min(Nr,k+1)
719 kp1Msk=1.
720 if (k.EQ.Nr) kp1Msk=0.
721
722 C-- Compute Vertical transport
723 #ifdef ALLOW_AIM
724 C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
725 IF ( k.EQ.1 .OR.
726 & (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
727 & ) THEN
728 #else
729 IF ( k.EQ.1 ) THEN
730 #endif
731
732 #ifdef ALLOW_AUTODIFF_TAMC
733 CADJ STORE rtrans(:,:) =
734 CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
735 CADJ STORE wfld(:,:) =
736 CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
737 #endif /* ALLOW_AUTODIFF_TAMC */
738
739 C- Surface interface :
740 DO j=1-Oly,sNy+Oly
741 DO i=1-Olx,sNx+Olx
742 rTransKp1(i,j) = kp1Msk*rTrans(i,j)
743 wFld(i,j) = 0.
744 rTrans(i,j) = 0.
745 fVerT(i,j,kUp) = 0.
746 ENDDO
747 ENDDO
748
749 ELSE
750
751 #ifdef ALLOW_AUTODIFF_TAMC
752 CADJ STORE rtrans(:,:) =
753 CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
754 CADJ STORE wfld(:,:) =
755 CADJ & comlev1_bibj_k_gad, key=kkey, byte=isbyte
756 #endif /* ALLOW_AUTODIFF_TAMC */
757
758 C- Interior interface :
759 DO j=1-Oly,sNy+Oly
760 DO i=1-Olx,sNx+Olx
761 rTransKp1(i,j) = kp1Msk*rTrans(i,j)
762 wFld(i,j) = wVel(i,j,k,bi,bj)
763 rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
764 & *deepFac2F(k)*rhoFacF(k)
765 & *maskC(i,j,k-1,bi,bj)
766 fVerT(i,j,kUp) = 0.
767 ENDDO
768 ENDDO
769
770 #ifdef ALLOW_GMREDI
771 C-- Residual transp = Bolus transp + Eulerian transp
772 IF (useGMRedi)
773 & CALL GMREDI_CALC_WFLOW(
774 U wFld, rTrans,
775 I k, bi, bj, myThid )
776 #endif /* ALLOW_GMREDI */
777
778 #ifdef ALLOW_AUTODIFF_TAMC
779 CADJ STORE localTijk(:,:,k)
780 CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
781 CADJ STORE rTrans(:,:)
782 CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
783 #endif /* ALLOW_AUTODIFF_TAMC */
784
785 C- Compute vertical advective flux in the interior:
786 IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
787 & .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
788 CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,
789 I dTtracerLev(k),rTrans,wFld,localTijk,
790 O fVerT(1-Olx,1-Oly,kUp), myThid )
791 ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
792 CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k),
793 I rTrans, wFld, localTijk,
794 O fVerT(1-Olx,1-Oly,kUp), myThid )
795 ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
796 CALL GAD_DST3_ADV_R( bi,bj,k, dTtracerLev(k),
797 I rTrans, wFld, localTijk,
798 O fVerT(1-Olx,1-Oly,kUp), myThid )
799 ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
800 CALL GAD_DST3FL_ADV_R( bi,bj,k, dTtracerLev(k),
801 I rTrans, wFld, localTijk,
802 O fVerT(1-Olx,1-Oly,kUp), myThid )
803 #ifndef ALLOW_AUTODIFF_TAMC
804 ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
805 CALL GAD_OS7MP_ADV_R( bi,bj,k, dTtracerLev(k),
806 I rTrans, wFld, localTijk,
807 O fVerT(1-Olx,1-Oly,kUp), myThid )
808 #endif
809 ELSE
810 STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
811 ENDIF
812
813 C- end Surface/Interior if bloc
814 ENDIF
815
816 #ifdef ALLOW_AUTODIFF_TAMC
817 CADJ STORE rTrans(:,:)
818 CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
819 CADJ STORE rTranskp1(:,:)
820 CADJ & = comlev1_bibj_k_gad, key=kkey, byte=isbyte
821 #endif /* ALLOW_AUTODIFF_TAMC */
822
823 C-- Divergence of vertical fluxes
824 DO j=1-Oly,sNy+Oly
825 DO i=1-Olx,sNx+Olx
826 localTij(i,j) = localTijk(i,j,k)
827 & -dTtracerLev(k)*recip_rhoFacC(k)
828 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
829 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
830 & *( fVerT(i,j,kDown)-fVerT(i,j,kUp)
831 & -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))
832 & )*rkSign
833 gTracer(i,j,k,bi,bj)=
834 & (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k)
835 ENDDO
836 ENDDO
837
838 #ifdef ALLOW_DIAGNOSTICS
839 IF ( useDiagnostics ) THEN
840 diagName = 'ADVr'//diagSufx
841 CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp),
842 & diagName, k,1, 2,bi,bj, myThid)
843 ENDIF
844 #endif
845
846 C-- End of K loop for vertical flux
847 ENDDO
848 C-- end of if not.implicitAdvection block
849 ENDIF
850
851 RETURN
852 END
  ViewVC Help
Powered by ViewVC 1.1.22