/[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.63 - (show annotations) (download)
Sun Oct 31 15:20:48 2010 UTC (13 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62u, checkpoint62t, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62o, checkpoint62n
Changes since 1.62: +14 -26 lines 
-remove MULTIDIM_OLD_VERSION code (was disabled 6.years ago).
-remove calls to OBCS_APPLY (applying future time-step OB values);
-use instead "maskInC" to keep current-time-step OB values in OB regions.
1 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_advection.F,v 1.62 2009/06/28 01:05:41 jmc Exp $
2 C $Name: $
3
4 #include "GAD_OPTIONS.h"
5 #undef OBCS_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, deltaTLev,
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_SIZE.h"
54 #include "W2_EXCH2_TOPOLOGY.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 deltaTLev(Nr)
74 _RL uVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
75 _RL vVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
76 _RL wVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
77 _RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
78 INTEGER bi,bj
79 _RL myTime
80 INTEGER myIter
81 INTEGER myThid
82
83 C !OUTPUT PARAMETERS: ==================================================
84 C gTracer :: tendency array
85 _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
86
87 C !LOCAL VARIABLES: ====================================================
88 C maskUp :: 2-D array for mask at W points
89 C maskLocW :: 2-D array for mask at West points
90 C maskLocS :: 2-D array for mask at South points
91 C [iMin,iMax]Upd :: loop range to update tracer field
92 C [jMin,jMax]Upd :: loop range to update tracer field
93 C i,j,k :: loop indices
94 C kUp :: index into 2 1/2D array, toggles between 1 and 2
95 C kDown :: index into 2 1/2D array, toggles between 2 and 1
96 C kp1 :: =k+1 for k<Nr, =Nr for k=Nr
97 C xA,yA :: areas of X and Y face of tracer cells
98 C uFld,vFld :: 2-D local copy of horizontal velocity, U,V components
99 C wFld :: 2-D local copy of vertical velocity
100 C uTrans,vTrans :: 2-D arrays of volume transports at U,V points
101 C rTrans :: 2-D arrays of volume transports at W points
102 C rTransKp1 :: vertical volume transport at interface k+1
103 C af :: 2-D array for horizontal advective flux
104 C afx :: 2-D array for horizontal advective flux, x direction
105 C afy :: 2-D array for horizontal advective flux, y direction
106 C fVerT :: 2 1/2D arrays for vertical advective flux
107 C localTij :: 2-D array, temporary local copy of tracer fld
108 C localTijk :: 3-D array, temporary local copy of tracer fld
109 C kp1Msk :: flag (0,1) for over-riding mask for W levels
110 C calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
111 C calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
112 C interiorOnly :: only update the interior of myTile, but not the edges
113 C overlapOnly :: only update the edges of myTile, but not the interior
114 C npass :: number of passes in multi-dimensional method
115 C ipass :: number of the current pass being made
116 C myTile :: variables used to determine which cube face
117 C nCFace :: owns a tile for cube grid runs using
118 C :: multi-dim advection.
119 C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
120 c _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121 _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122 _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123 INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
124 INTEGER i,j,k,kUp,kDown
125 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127 _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128 _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129 _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133 _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134 _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135 _RL afx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136 _RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
138 _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139 _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
140 _RL kp1Msk
141 LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
142 LOGICAL interiorOnly, overlapOnly
143 INTEGER npass, ipass
144 INTEGER nCFace
145 LOGICAL N_edge, S_edge, E_edge, W_edge
146 #ifdef ALLOW_EXCH2
147 INTEGER myTile
148 #endif
149 #ifdef ALLOW_DIAGNOSTICS
150 CHARACTER*8 diagName
151 CHARACTER*4 diagSufx
152 LOGICAL doDiagAdvX, doDiagAdvY, doDiagAdvR
153 C- Functions:
154 CHARACTER*4 GAD_DIAG_SUFX
155 EXTERNAL GAD_DIAG_SUFX
156 LOGICAL DIAGNOSTICS_IS_ON
157 EXTERNAL DIAGNOSTICS_IS_ON
158 #endif
159 CEOP
160
161 #ifdef ALLOW_AUTODIFF_TAMC
162 act0 = tracerIdentity
163 max0 = maxpass
164 act1 = bi - myBxLo(myThid)
165 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
166 act2 = bj - myByLo(myThid)
167 max2 = myByHi(myThid) - myByLo(myThid) + 1
168 act3 = myThid - 1
169 max3 = nTx*nTy
170 act4 = ikey_dynamics - 1
171 igadkey = act0
172 & + act1*max0
173 & + act2*max0*max1
174 & + act3*max0*max1*max2
175 & + act4*max0*max1*max2*max3
176 IF (tracerIdentity.GT.maxpass) THEN
177 print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
178 STOP 'maxpass seems smaller than tracerIdentity'
179 ENDIF
180 #endif /* ALLOW_AUTODIFF_TAMC */
181
182 #ifdef ALLOW_DIAGNOSTICS
183 C-- Set diagnostics flags and suffix for the current tracer
184 doDiagAdvX = .FALSE.
185 doDiagAdvY = .FALSE.
186 doDiagAdvR = .FALSE.
187 IF ( useDiagnostics ) THEN
188 diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
189 diagName = 'ADVx'//diagSufx
190 doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid )
191 diagName = 'ADVy'//diagSufx
192 doDiagAdvY = DIAGNOSTICS_IS_ON( diagName, myThid )
193 diagName = 'ADVr'//diagSufx
194 doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )
195 ENDIF
196 #endif
197
198 C-- Set up work arrays with valid (i.e. not NaN) values
199 C These inital values do not alter the numerical results. They
200 C just ensure that all memory references are to valid floating
201 C point numbers. This prevents spurious hardware signals due to
202 C uninitialised but inert locations.
203 DO j=1-OLy,sNy+OLy
204 DO i=1-OLx,sNx+OLx
205 xA(i,j) = 0. _d 0
206 yA(i,j) = 0. _d 0
207 uTrans(i,j) = 0. _d 0
208 vTrans(i,j) = 0. _d 0
209 rTrans(i,j) = 0. _d 0
210 fVerT(i,j,1) = 0. _d 0
211 fVerT(i,j,2) = 0. _d 0
212 rTransKp1(i,j)= 0. _d 0
213 #ifdef ALLOW_AUTODIFF_TAMC
214 localTij(i,j) = 0. _d 0
215 wfld(i,j) = 0. _d 0
216 #endif
217 ENDDO
218 ENDDO
219
220 C-- Set tile-specific parameters for horizontal fluxes
221 IF (useCubedSphereExchange) THEN
222 npass = 3
223 #ifdef ALLOW_AUTODIFF_TAMC
224 IF ( npass.GT.maxcube ) STOP 'maxcube needs to be = 3'
225 #endif
226 #ifdef ALLOW_EXCH2
227 myTile = W2_myTileList(bi,bj)
228 nCFace = exch2_myFace(myTile)
229 N_edge = exch2_isNedge(myTile).EQ.1
230 S_edge = exch2_isSedge(myTile).EQ.1
231 E_edge = exch2_isEedge(myTile).EQ.1
232 W_edge = exch2_isWedge(myTile).EQ.1
233 #else
234 nCFace = bi
235 N_edge = .TRUE.
236 S_edge = .TRUE.
237 E_edge = .TRUE.
238 W_edge = .TRUE.
239 #endif
240 ELSE
241 npass = 2
242 nCFace = 0
243 N_edge = .FALSE.
244 S_edge = .FALSE.
245 E_edge = .FALSE.
246 W_edge = .FALSE.
247 ENDIF
248
249 C-- Start of k loop for horizontal fluxes
250 DO k=1,Nr
251 #ifdef ALLOW_AUTODIFF_TAMC
252 kkey = (igadkey-1)*Nr + k
253 CADJ STORE tracer(:,:,k,bi,bj) =
254 CADJ & comlev1_bibj_k_gad, key=kkey, kind=isbyte
255 #endif /* ALLOW_AUTODIFF_TAMC */
256
257 C-- Get temporary terms used by tendency routines
258 CALL CALC_COMMON_FACTORS (
259 I uVel, vVel,
260 O uFld, vFld, uTrans, vTrans, xA, yA,
261 I k,bi,bj, myThid )
262
263 #ifdef ALLOW_GMREDI
264 C-- Residual transp = Bolus transp + Eulerian transp
265 IF (useGMRedi)
266 & CALL GMREDI_CALC_UVFLOW(
267 U uFld, vFld, uTrans, vTrans,
268 I k, bi, bj, myThid )
269 #endif /* ALLOW_GMREDI */
270
271 C-- Make local copy of tracer array and mask West & South
272 DO j=1-OLy,sNy+OLy
273 DO i=1-OLx,sNx+OLx
274 localTij(i,j)=tracer(i,j,k,bi,bj)
275 maskLocW(i,j)=_maskW(i,j,k,bi,bj)
276 maskLocS(i,j)=_maskS(i,j,k,bi,bj)
277 ENDDO
278 ENDDO
279
280 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
281 IF (useCubedSphereExchange) THEN
282 withSigns = .FALSE.
283 CALL FILL_CS_CORNER_UV_RS(
284 & withSigns, maskLocW,maskLocS, bi,bj, myThid )
285 ENDIF
286 cph-exch2#endif
287
288 C-- Multiple passes for different directions on different tiles
289 C-- For cube need one pass for each of red, green and blue axes.
290 DO ipass=1,npass
291 #ifdef ALLOW_AUTODIFF_TAMC
292 passkey = ipass
293 & + (k-1) *maxpass
294 & + (igadkey-1)*maxpass*Nr
295 IF (npass .GT. maxpass) THEN
296 STOP 'GAD_ADVECTION: npass > maxcube. check tamc.h'
297 ENDIF
298 #endif /* ALLOW_AUTODIFF_TAMC */
299
300 interiorOnly = .FALSE.
301 overlapOnly = .FALSE.
302 IF (useCubedSphereExchange) THEN
303 C- CubedSphere : pass 3 times, with partial update of local tracer field
304 IF (ipass.EQ.1) THEN
305 overlapOnly = MOD(nCFace,3).EQ.0
306 interiorOnly = MOD(nCFace,3).NE.0
307 calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
308 calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
309 ELSEIF (ipass.EQ.2) THEN
310 overlapOnly = MOD(nCFace,3).EQ.2
311 interiorOnly = MOD(nCFace,3).EQ.1
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 ELSE
315 interiorOnly = .TRUE.
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, kind=isbyte
338 CADJ STORE af(:,:) =
339 CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=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 C- Internal exchange for calculations in X
351 IF ( overlapOnly ) THEN
352 CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
353 & localTij, bi,bj, myThid )
354 ENDIF
355
356 #ifdef ALLOW_AUTODIFF_TAMC
357 # ifndef DISABLE_MULTIDIM_ADVECTION
358 CADJ STORE localTij(:,:) =
359 CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
360 # endif
361 #endif /* ALLOW_AUTODIFF_TAMC */
362
363 IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
364 & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
365 CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
366 I deltaTLev(k),uTrans,uFld,localTij,
367 O af, myThid )
368 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
369 CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
370 I uTrans, uFld, maskLocW, localTij,
371 O af, myThid )
372 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
373 CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
374 I uTrans, uFld, maskLocW, localTij,
375 O af, myThid )
376 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
377 CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
378 I uTrans, uFld, maskLocW, localTij,
379 O af, myThid )
380 #ifndef ALLOW_AUTODIFF_TAMC
381 ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
382 CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
383 I uTrans, uFld, maskLocW, localTij,
384 O af, myThid )
385 #endif
386 ELSE
387 STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
388 ENDIF
389
390 C- Internal exchange for next calculations in Y
391 IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
392 CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
393 & localTij, bi,bj, myThid )
394 ENDIF
395
396 C- Advective flux in X : done
397 ENDIF
398
399 C- Update the local tracer field where needed:
400 C use "maksInC" to prevent updating tracer field in OB regions
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 & -deltaTLev(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 & )*maskInC(i,j,bi,bj)
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 & -deltaTLev(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 & )*maskInC(i,j,bi,bj)
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 & -deltaTLev(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 & )*maskInC(i,j,bi,bj)
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 #ifdef OBCS_MULTIDIM_OLD_VERSION
464 C- Apply open boundary conditions
465 IF ( useOBCS ) THEN
466 IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
467 CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
468 ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
469 CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
470 #ifdef ALLOW_PTRACERS
471 ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
472 CALL OBCS_APPLY_PTRACER( bi, bj, k,
473 & tracerIdentity-GAD_TR1+1, localTij, myThid )
474 #endif /* ALLOW_PTRACERS */
475 ENDIF
476 ENDIF
477 #endif /* OBCS_MULTIDIM_OLD_VERSION */
478 #endif /* ALLOW_OBCS */
479
480 C- end if/else update overlap-Only
481 ENDIF
482
483 C-- End of X direction
484 ENDIF
485
486 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
487 C-- Y direction
488 cph-test
489 C- Advective flux in Y
490 DO j=1-Oly,sNy+Oly
491 DO i=1-Olx,sNx+Olx
492 af(i,j) = 0.
493 ENDDO
494 ENDDO
495 C
496 #ifdef ALLOW_AUTODIFF_TAMC
497 # ifndef DISABLE_MULTIDIM_ADVECTION
498 CADJ STORE localTij(:,:) =
499 CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
500 CADJ STORE af(:,:) =
501 CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
502 # endif
503 #endif /* ALLOW_AUTODIFF_TAMC */
504 C
505 IF (calc_fluxes_Y) THEN
506
507 C- Do not compute fluxes if
508 C a) needed in overlap only
509 C and b) the overlap of myTile are not cube-face edges
510 IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
511
512 C- Internal exchange for calculations in Y
513 IF ( overlapOnly ) THEN
514 CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
515 & localTij, bi,bj, myThid )
516 ENDIF
517
518 C- Advective flux in Y
519 DO j=1-Oly,sNy+Oly
520 DO i=1-Olx,sNx+Olx
521 af(i,j) = 0.
522 ENDDO
523 ENDDO
524
525 #ifdef ALLOW_AUTODIFF_TAMC
526 #ifndef DISABLE_MULTIDIM_ADVECTION
527 CADJ STORE localTij(:,:) =
528 CADJ & comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
529 #endif
530 #endif /* ALLOW_AUTODIFF_TAMC */
531
532 IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
533 & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
534 CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
535 I deltaTLev(k),vTrans,vFld,localTij,
536 O af, myThid )
537 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
538 CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
539 I vTrans, vFld, maskLocS, localTij,
540 O af, myThid )
541 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
542 CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
543 I vTrans, vFld, maskLocS, localTij,
544 O af, myThid )
545 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
546 CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
547 I vTrans, vFld, maskLocS, localTij,
548 O af, myThid )
549 #ifndef ALLOW_AUTODIFF_TAMC
550 ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
551 CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
552 I vTrans, vFld, maskLocS, localTij,
553 O af, myThid )
554 #endif
555 ELSE
556 STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
557 ENDIF
558
559 C- Internal exchange for next calculations in X
560 IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
561 CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
562 & localTij, bi,bj, myThid )
563 ENDIF
564
565 C- Advective flux in Y : done
566 ENDIF
567
568 C- Update the local tracer field where needed:
569 C use "maksInC" to prevent updating tracer field in OB regions
570
571 C update in overlap-Only
572 IF ( overlapOnly ) THEN
573 jMinUpd = 1-Oly+1
574 jMaxUpd = sNy+Oly-1
575 C- notes: these 2 lines below have no real effect (because recip_hFac=0
576 C in corner region) but safer to keep them.
577 IF ( S_edge ) jMinUpd = 1
578 IF ( N_edge ) jMaxUpd = sNy
579
580 IF ( W_edge ) THEN
581 DO j=jMinUpd,jMaxUpd
582 DO i=1-Olx,0
583 localTij(i,j) = localTij(i,j)
584 & -deltaTLev(k)*recip_rhoFacC(k)
585 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
586 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
587 & *( af(i,j+1)-af(i,j)
588 & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
589 & )*maskInC(i,j,bi,bj)
590 ENDDO
591 ENDDO
592 ENDIF
593 IF ( E_edge ) THEN
594 DO j=jMinUpd,jMaxUpd
595 DO i=sNx+1,sNx+Olx
596 localTij(i,j) = localTij(i,j)
597 & -deltaTLev(k)*recip_rhoFacC(k)
598 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
599 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
600 & *( af(i,j+1)-af(i,j)
601 & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
602 & )*maskInC(i,j,bi,bj)
603 ENDDO
604 ENDDO
605 ENDIF
606
607 ELSE
608 C do not only update the overlap
609 iMinUpd = 1-Olx
610 iMaxUpd = sNx+Olx
611 IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
612 IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
613 DO j=1-Oly+1,sNy+Oly-1
614 DO i=iMinUpd,iMaxUpd
615 localTij(i,j) = localTij(i,j)
616 & -deltaTLev(k)*recip_rhoFacC(k)
617 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
618 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
619 & *( af(i,j+1)-af(i,j)
620 & -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
621 & )*maskInC(i,j,bi,bj)
622 ENDDO
623 ENDDO
624 C- keep advective flux (for diagnostics)
625 DO j=1-Oly,sNy+Oly
626 DO i=1-Olx,sNx+Olx
627 afy(i,j) = af(i,j)
628 ENDDO
629 ENDDO
630
631 #ifdef ALLOW_OBCS
632 #ifdef OBCS_MULTIDIM_OLD_VERSION
633 C- Apply open boundary conditions
634 IF (useOBCS) THEN
635 IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
636 CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid )
637 ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
638 CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid )
639 #ifdef ALLOW_PTRACERS
640 ELSEIF (tracerIdentity.GE.GAD_TR1) THEN
641 CALL OBCS_APPLY_PTRACER( bi, bj, k,
642 & tracerIdentity-GAD_TR1+1, localTij, myThid )
643 #endif /* ALLOW_PTRACERS */
644 ENDIF
645 ENDIF
646 #endif /* OBCS_MULTIDIM_OLD_VERSION */
647 #endif /* ALLOW_OBCS */
648
649 C end if/else update overlap-Only
650 ENDIF
651
652 C-- End of Y direction
653 ENDIF
654
655 C-- End of ipass loop
656 ENDDO
657
658 IF ( implicitAdvection ) THEN
659 C- explicit advection is done ; store tendency in gTracer:
660 DO j=1-Oly,sNy+Oly
661 DO i=1-Olx,sNx+Olx
662 gTracer(i,j,k,bi,bj)=
663 & (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
664 ENDDO
665 ENDDO
666 ELSE
667 C- horizontal advection done; store intermediate result in 3D array:
668 DO j=1-Oly,sNy+Oly
669 DO i=1-Olx,sNx+Olx
670 localTijk(i,j,k)=localTij(i,j)
671 ENDDO
672 ENDDO
673 ENDIF
674
675 #ifdef ALLOW_DIAGNOSTICS
676 IF ( doDiagAdvX ) THEN
677 diagName = 'ADVx'//diagSufx
678 CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)
679 ENDIF
680 IF ( doDiagAdvY ) THEN
681 diagName = 'ADVy'//diagSufx
682 CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)
683 ENDIF
684 #endif
685
686 #ifdef ALLOW_DEBUG
687 IF ( debugLevel .GE. debLevB
688 & .AND. tracerIdentity.EQ.GAD_TEMPERATURE
689 & .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
690 & .AND. nPx.EQ.1 .AND. nPy.EQ.1
691 & .AND. useCubedSphereExchange ) THEN
692 CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_ADVECTION',
693 & afx,afy, k, standardMessageUnit,bi,bj,myThid )
694 ENDIF
695 #endif /* ALLOW_DEBUG */
696
697 C-- End of K loop for horizontal fluxes
698 ENDDO
699
700 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
701
702 IF ( .NOT.implicitAdvection ) THEN
703 C-- Start of k loop for vertical flux
704 DO k=Nr,1,-1
705 #ifdef ALLOW_AUTODIFF_TAMC
706 kkey = (igadkey-1)*Nr + (Nr-k+1)
707 #endif /* ALLOW_AUTODIFF_TAMC */
708 C-- kUp Cycles through 1,2 to point to w-layer above
709 C-- kDown Cycles through 2,1 to point to w-layer below
710 kUp = 1+MOD(k+1,2)
711 kDown= 1+MOD(k,2)
712 c kp1=min(Nr,k+1)
713 kp1Msk=1.
714 if (k.EQ.Nr) kp1Msk=0.
715
716 #ifdef ALLOW_AUTODIFF_TAMC
717 CADJ STORE rtrans(:,:) =
718 CADJ & comlev1_bibj_k_gad, key=kkey, kind=isbyte
719 cphCADJ STORE wfld(:,:) =
720 cphCADJ & comlev1_bibj_k_gad, key=kkey, kind=isbyte
721 #endif
722
723 C-- Compute Vertical transport
724 #ifdef ALLOW_AIM
725 C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
726 IF ( k.EQ.1 .OR.
727 & (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr)
728 & ) THEN
729 #else
730 IF ( k.EQ.1 ) THEN
731 #endif
732
733 #ifdef ALLOW_AUTODIFF_TAMC
734 cphmultiCADJ STORE wfld(:,:) =
735 cphmultiCADJ & comlev1_bibj_k_gad, key=kkey, kind=isbyte
736 #endif /* ALLOW_AUTODIFF_TAMC */
737
738 C- Surface interface :
739 DO j=1-Oly,sNy+Oly
740 DO i=1-Olx,sNx+Olx
741 rTransKp1(i,j) = kp1Msk*rTrans(i,j)
742 wFld(i,j) = 0.
743 rTrans(i,j) = 0.
744 fVerT(i,j,kUp) = 0.
745 ENDDO
746 ENDDO
747
748 ELSE
749
750 #ifdef ALLOW_AUTODIFF_TAMC
751 cphmultiCADJ STORE wfld(:,:) =
752 cphmultiCADJ & comlev1_bibj_k_gad, key=kkey, kind=isbyte
753 #endif /* ALLOW_AUTODIFF_TAMC */
754
755 C- Interior interface :
756 DO j=1-Oly,sNy+Oly
757 DO i=1-Olx,sNx+Olx
758 rTransKp1(i,j) = kp1Msk*rTrans(i,j)
759 wFld(i,j) = wVel(i,j,k,bi,bj)
760 rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
761 & *deepFac2F(k)*rhoFacF(k)
762 & *maskC(i,j,k-1,bi,bj)
763 fVerT(i,j,kUp) = 0.
764 ENDDO
765 ENDDO
766
767 #ifdef ALLOW_GMREDI
768 C-- Residual transp = Bolus transp + Eulerian transp
769 IF (useGMRedi)
770 & CALL GMREDI_CALC_WFLOW(
771 U wFld, rTrans,
772 I k, bi, bj, myThid )
773 #endif /* ALLOW_GMREDI */
774
775 #ifdef ALLOW_AUTODIFF_TAMC
776 cphmultiCADJ STORE localTijk(:,:,k)
777 cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, kind=isbyte
778 cphmultiCADJ STORE rTrans(:,:)
779 cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, kind=isbyte
780 #endif /* ALLOW_AUTODIFF_TAMC */
781
782 C- Compute vertical advective flux in the interior:
783 IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
784 & .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
785 CALL GAD_DST2U1_ADV_R( bi,bj,k, advectionScheme,
786 I deltaTLev(k),rTrans,wFld,localTijk,
787 O fVerT(1-Olx,1-Oly,kUp), myThid )
788 ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
789 CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
790 I rTrans, wFld, localTijk,
791 O fVerT(1-Olx,1-Oly,kUp), myThid )
792 ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
793 CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
794 I rTrans, wFld, localTijk,
795 O fVerT(1-Olx,1-Oly,kUp), myThid )
796 ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
797 CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
798 I rTrans, wFld, localTijk,
799 O fVerT(1-Olx,1-Oly,kUp), myThid )
800 #ifndef ALLOW_AUTODIFF_TAMC
801 ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
802 CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
803 I rTrans, wFld, localTijk,
804 O fVerT(1-Olx,1-Oly,kUp), myThid )
805 #endif
806 ELSE
807 STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
808 ENDIF
809
810 C- end Surface/Interior if bloc
811 ENDIF
812
813 #ifdef ALLOW_AUTODIFF_TAMC
814 cphmultiCADJ STORE rTrans(:,:)
815 cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, kind=isbyte
816 cphmultiCADJ STORE rTranskp1(:,:)
817 cphmultiCADJ & = comlev1_bibj_k_gad, key=kkey, kind=isbyte
818 cph --- following storing of fVerT is critical for correct
819 cph --- gradient with multiDimAdvection
820 cph --- Without it, kDown component is not properly recomputed
821 cph --- This is a TAF bug (and no warning available)
822 CADJ STORE fVerT(:,:,:)
823 CADJ & = comlev1_bibj_k_gad, key=kkey, kind=isbyte
824 #endif /* ALLOW_AUTODIFF_TAMC */
825
826 C-- Divergence of vertical fluxes
827 DO j=1-Oly,sNy+Oly
828 DO i=1-Olx,sNx+Olx
829 localTij(i,j) = localTijk(i,j,k)
830 & -deltaTLev(k)*recip_rhoFacC(k)
831 & *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
832 & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
833 & *( fVerT(i,j,kDown)-fVerT(i,j,kUp)
834 & -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))
835 & )*rkSign
836 gTracer(i,j,k,bi,bj)=
837 & (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
838 ENDDO
839 ENDDO
840
841 #ifdef ALLOW_DIAGNOSTICS
842 IF ( doDiagAdvR ) THEN
843 diagName = 'ADVr'//diagSufx
844 CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp),
845 & diagName, k,1, 2,bi,bj, myThid)
846 ENDIF
847 #endif
848
849 C-- End of K loop for vertical flux
850 ENDDO
851 C-- end of if not.implicitAdvection block
852 ENDIF
853
854 RETURN
855 END
  ViewVC Help
Powered by ViewVC 1.1.22