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

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

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


Revision 1.17 - (show annotations) (download)
Tue Sep 18 21:29:36 2012 UTC (11 years, 8 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64f, checkpoint64
Changes since 1.16: +5 -4 lines
- fix ADVr diagnostics in case of implicit vertical advection.

1 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_implicit_r.F,v 1.16 2012/09/11 01:32:02 heimbach Exp $
2 C $Name: $
3
4 #include "GAD_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: GAD_IMPLICIT_R
8 C !INTERFACE:
9 SUBROUTINE GAD_IMPLICIT_R(
10 I implicitAdvection, advectionScheme, tracerIdentity,
11 I deltaTLev,
12 I kappaRX, recip_hFac, wVel, tracer,
13 U gTracer,
14 I bi, bj, myTime, myIter, myThid )
15 C !DESCRIPTION:
16 C Solve implicitly vertical advection and diffusion for one tracer.
17
18 C !USES:
19 IMPLICIT NONE
20 C == Global data ==
21 #include "SIZE.h"
22 #include "EEPARAMS.h"
23 #include "PARAMS.h"
24 #include "GRID.h"
25 #include "SURFACE.h"
26 #include "GAD.h"
27
28 C !INPUT/OUTPUT PARAMETERS:
29 C == Routine Arguments ==
30 C implicitAdvection :: if True, treat vertical advection implicitly
31 C advectionScheme :: advection scheme to use
32 C tracerIdentity :: Identifier for the tracer
33 C kappaRX :: 3-D array for vertical diffusion coefficient
34 C recip_hFac :: inverse of cell open-depth factor
35 C wVel :: vertical component of the velcity field
36 C tracer :: tracer field at current time step
37 C gTracer :: future tracer field
38 C bi,bj :: tile indices
39 C myTime :: current time
40 C myIter :: current iteration number
41 C myThid :: thread number
42 LOGICAL implicitAdvection
43 INTEGER advectionScheme
44 INTEGER tracerIdentity
45 _RL deltaTLev(Nr)
46 _RL kappaRX (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
47 _RS recip_hFac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
48 _RL wVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
49 _RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
50 _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
51 INTEGER bi, bj
52 _RL myTime
53 INTEGER myIter, myThid
54
55 #ifdef ALLOW_DIAGNOSTICS
56 C !FUNCTIONS:
57 CHARACTER*4 GAD_DIAG_SUFX
58 EXTERNAL GAD_DIAG_SUFX
59 LOGICAL DIAGNOSTICS_IS_ON
60 EXTERNAL DIAGNOSTICS_IS_ON
61 #endif
62
63 C !LOCAL VARIABLES:
64 C == Local variables ==
65 C iMin,iMax,jMin,jMax :: computational domain
66 C i,j,k :: loop indices
67 C a5d :: 2nd lower diagonal of the pentadiagonal matrix
68 C b5d :: 1rst lower diagonal of the pentadiagonal matrix
69 C c5d :: main diagonal of the pentadiagonal matrix
70 C d5d :: 1rst upper diagonal of the pentadiagonal matrix
71 C e5d :: 2nd upper diagonal of the pentadiagonal matrix
72 C rTrans :: vertical volume transport at interface k
73 C localTijk :: local copy of tracer (for Non-Lin Adv.Scheme)
74 C diagonalNumber :: number of non-zero diagonals in the matrix
75 C errCode :: > 0 if singular matrix
76 INTEGER iMin,iMax,jMin,jMax
77 PARAMETER( iMin = 1, iMax = sNx )
78 PARAMETER( jMin = 1, jMax = sNy )
79 INTEGER i,j,k
80 INTEGER diagonalNumber, errCode
81 _RL a5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
82 _RL b5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
83 _RL c5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
84 _RL d5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
85 _RL e5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
86 _RL wFld (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
87 _RL rTrans (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
88 _RL localTijk(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
89 #ifdef ALLOW_DIAGNOSTICS
90 CHARACTER*8 diagName
91 CHARACTER*4 diagSufx
92 LOGICAL diagDif, diagAdv
93 INTEGER km1, km2, kp1, kp2
94 _RL df (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
95 _RL af (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
96 _RL div(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
97 _RL flx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
98 #endif
99 CEOP
100
101 C-- no need to solve anything with only 1 level:
102 IF (Nr.GT.1) THEN
103
104 C-- Initialise
105 DO k=1,Nr
106 DO j=1-Oly,sNy+Oly
107 DO i=1-Olx,sNx+Olx
108 a5d(i,j,k) = 0. _d 0
109 b5d(i,j,k) = 0. _d 0
110 c5d(i,j,k) = 1. _d 0
111 d5d(i,j,k) = 0. _d 0
112 e5d(i,j,k) = 0. _d 0
113 ENDDO
114 ENDDO
115 ENDDO
116 diagonalNumber = 1
117
118 C-- Non-Linear Advection scheme: keep a local copy of tracer field
119 IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.
120 & advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
121 IF ( multiDimAdvection ) THEN
122 DO k=1,Nr
123 DO j=1-Oly,sNy+Oly
124 DO i=1-Olx,sNx+Olx
125 localTijk(i,j,k) = gTracer(i,j,k,bi,bj)
126 ENDDO
127 ENDDO
128 ENDDO
129 ELSE
130 DO k=1,Nr
131 DO j=1-Oly,sNy+Oly
132 DO i=1-Olx,sNx+Olx
133 localTijk(i,j,k) = tracer(i,j,k,bi,bj)
134 ENDDO
135 ENDDO
136 ENDDO
137 ENDIF
138 ENDIF
139
140 IF (implicitDiffusion) THEN
141 C-- set the tri-diagonal matrix to solve the implicit diffusion problem
142 diagonalNumber = 3
143 C- 1rst lower diagonal :
144 DO k=2,Nr
145 DO j=jMin,jMax
146 DO i=iMin,iMax
147 b5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k-1,bi,bj)
148 & *recip_hFac(i,j,k)*recip_drF(k)
149 & *kappaRX(i,j, k )*recip_drC( k )
150 ENDDO
151 ENDDO
152 ENDDO
153 C- 1rst upper diagonal :
154 DO k=1,Nr-1
155 DO j=jMin,jMax
156 DO i=iMin,iMax
157 d5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k+1,bi,bj)
158 & *recip_hFac(i,j,k)*recip_drF(k)
159 & *KappaRX(i,j,k+1)*recip_drC(k+1)
160 ENDDO
161 ENDDO
162 ENDDO
163 C- Main diagonal :
164 DO k=1,Nr
165 DO j=jMin,jMax
166 DO i=iMin,iMax
167 c5d(i,j,k) = 1. _d 0 - b5d(i,j,k) - d5d(i,j,k)
168 ENDDO
169 ENDDO
170 ENDDO
171
172 C-- end if implicitDiffusion
173 ENDIF
174
175 IF (implicitAdvection) THEN
176
177 DO k=Nr,1,-1
178
179 C-- Compute transport
180 IF (k.EQ.1) THEN
181 DO j=1-Oly,sNy+Oly
182 DO i=1-Olx,sNx+Olx
183 wFld(i,j) = 0. _d 0
184 rTrans(i,j) = 0. _d 0
185 ENDDO
186 ENDDO
187 ELSE
188 DO j=1-Oly,sNy+Oly
189 DO i=1-Olx,sNx+Olx
190 wFld(i,j) = wVel(i,j,k,bi,bj)
191 #ifdef ALLOW_AUTODIFF_OPENAD
192 rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)%v*maskC(i,j,k-1,bi,bj)
193 #else
194 rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)*maskC(i,j,k-1,bi,bj)
195 #endif
196 ENDDO
197 ENDDO
198 #ifdef ALLOW_GMREDI
199 C-- Residual transp = Bolus transp + Eulerian transp
200 IF (useGMRedi)
201 & CALL GMREDI_CALC_WFLOW(
202 U wFld, rTrans,
203 I k, bi, bj, myThid )
204 #endif /* ALLOW_GMREDI */
205 ENDIF
206
207 #ifdef ALLOW_AIM
208 C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
209 IF ( k.GE.2 .AND.
210 & (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.k.LT.Nr)
211 & ) THEN
212 #else
213 IF ( k.GE.2 ) THEN
214 #endif
215
216 IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
217 diagonalNumber = 3
218 CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
219 I deltaTLev, rTrans, recip_hFac,
220 U b5d, c5d, d5d,
221 I myThid )
222 ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
223 & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
224 diagonalNumber = 3
225 CALL GAD_DST2U1_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
226 I advectionScheme, deltaTLev,
227 I rTrans, recip_hFac,
228 U b5d, c5d, d5d,
229 I myThid )
230 ELSEIF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
231 diagonalNumber = 3
232 CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
233 I deltaTLev, rTrans, recip_hFac, localTijk,
234 U b5d, c5d, d5d,
235 I myThid )
236 ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_3RD
237 & .OR. advectionScheme.EQ.ENUM_CENTERED_4TH
238 & .OR. advectionScheme.EQ.ENUM_DST3 ) THEN
239 diagonalNumber = 5
240 CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
241 I advectionScheme, deltaTLev,
242 I rTrans, recip_hFac,
243 U a5d, b5d, c5d, d5d, e5d,
244 I myThid )
245 ELSEIF ( advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
246 diagonalNumber = 5
247 CALL GAD_DST3FL_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
248 I deltaTLev, rTrans, recip_hFac, localTijk,
249 U a5d, b5d, c5d, d5d, e5d,
250 I myThid )
251 ELSE
252 STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'
253 ENDIF
254
255 ENDIF
256
257 C-- end k loop
258 ENDDO
259
260 C-- end if implicitAdvection
261 ENDIF
262
263 IF ( diagonalNumber .EQ. 3 ) THEN
264 C-- Solve tri-diagonal system :
265 CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
266 I b5d, c5d, d5d,
267 U gTracer,
268 O errCode,
269 I bi, bj, myThid )
270 IF (errCode.GE.1) THEN
271 STOP 'GAD_IMPLICIT_R: error when solving 3-Diag problem'
272 ENDIF
273 ELSEIF ( diagonalNumber .EQ. 5 ) THEN
274 C-- Solve penta-diagonal system :
275 CALL SOLVE_PENTADIAGONAL( iMin,iMax, jMin,jMax,
276 I a5d, b5d, c5d, d5d, e5d,
277 U gTracer,
278 O errCode,
279 I bi, bj, myThid )
280 IF (errCode.GE.1) THEN
281 STOP 'GAD_IMPLICIT_R: error when solving 5-Diag problem'
282 ENDIF
283 ELSEIF ( diagonalNumber .NE. 1 ) THEN
284 STOP 'GAD_IMPLICIT_R: no solver available'
285 ENDIF
286
287 #ifdef ALLOW_DIAGNOSTICS
288 C-- Set diagnostic suffix for the current tracer
289 IF ( useDiagnostics ) THEN
290 diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
291 diagName = 'DFrI'//diagSufx
292 diagDif = implicitDiffusion
293 IF ( diagDif ) diagDif = DIAGNOSTICS_IS_ON(diagName,myThid)
294 diagName = 'ADVr'//diagSufx
295 diagAdv = implicitAdvection
296 IF ( diagAdv ) diagAdv = DIAGNOSTICS_IS_ON(diagName,myThid)
297
298 IF ( diagDif .OR. diagAdv ) THEN
299 DO j=1-OLy,sNy+OLy
300 DO i=1-OLx,sNx+OLx
301 flx(i,j) = 0. _d 0
302 ENDDO
303 ENDDO
304 DO k= Nr,1,-1
305 IF ( implicitDiffusion .AND. k.GE.2 ) THEN
306 DO j=jMin,jMax
307 DO i=iMin,iMax
308 df(i,j) =
309 #ifdef ALLOW_AUTODIFF_OPENAD
310 & -rA(i,j,bi,bj)%v
311 #else
312 & -rA(i,j,bi,bj)
313 #endif
314 & * KappaRX(i,j,k)*recip_drC(k)*rkSign
315 & * (gTracer(i,j,k,bi,bj) - gTracer(i,j,k-1,bi,bj))
316 & * maskC(i,j,k,bi,bj)
317 & * maskC(i,j,k-1,bi,bj)
318 ENDDO
319 ENDDO
320 ELSE
321 DO j=1-OLy,sNy+OLy
322 DO i=1-OLx,sNx+OLx
323 df(i,j) = 0. _d 0
324 ENDDO
325 ENDDO
326 ENDIF
327 C- Note: Needs to explicitly increment counter (call DIAGNOSTICS_COUNT)
328 C since skipping k=1 DIAGNOSTICS_FILL call.
329 IF ( diagDif .AND. k.GE.2 ) THEN
330 diagName = 'DFrI'//diagSufx
331 CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
332 IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
333 ENDIF
334 IF ( diagAdv ) THEN
335 km1=MAX(1,k-1)
336 km2=MAX(1,k-2)
337 kp1=MIN(Nr,k+1)
338 kp2=MIN(Nr,k+2)
339 C-- Flux_divergence*deltaT = Tr^n - Tr^n+1 = [A-I](Tr^n+1)
340 C = deltaT*rkSign*[ Flx_k+1 - Flx_k ]/dz
341 DO j=jMin,jMax
342 DO i=iMin,iMax
343 div(i,j) = gTracer(i,j,k,bi,bj)*( c5d(i,j,k) - 1. _d 0 )
344 & + gTracer(i,j,km1,bi,bj)*b5d(i,j,k)
345 & + gTracer(i,j,kp1,bi,bj)*d5d(i,j,k)
346 ENDDO
347 ENDDO
348 IF ( diagonalNumber .EQ. 5 ) THEN
349 DO j=jMin,jMax
350 DO i=iMin,iMax
351 div(i,j) = div(i,j)
352 & + gTracer(i,j,km2,bi,bj)*a5d(i,j,k)
353 & + gTracer(i,j,kp2,bi,bj)*e5d(i,j,k)
354 ENDDO
355 ENDDO
356 ENDIF
357 #ifdef NONLIN_FRSURF
358 IF ( nonlinFreeSurf.GT.0 ) THEN
359 C-- use future hFac to stay consistent with solver matrix
360 IF ( select_rStar.GT.0 ) THEN
361 DO j=jMin,jMax
362 DO i=iMin,iMax
363 div(i,j) = div(i,j)*h0FacC(i,j,k,bi,bj)*drF(k)
364 & *rStarFacC(i,j,bi,bj)
365 ENDDO
366 ENDDO
367 ELSEIF ( selectSigmaCoord.NE.0 ) THEN
368 DO j=jMin,jMax
369 DO i=iMin,iMax
370 div(i,j) = div(i,j)*(
371 & _hFacC(i,j,k,bi,bj)*drF(k)
372 & + dBHybSigF(k)*dEtaHdt(i,j,bi,bj)*deltaTfreesurf
373 & )
374 ENDDO
375 ENDDO
376 ELSE
377 DO j=jMin,jMax
378 DO i=iMin,iMax
379 IF ( k.EQ.kSurfC(i,j,bi,bj) ) THEN
380 div(i,j) = div(i,j)*hFac_surfC(i,j,bi,bj)*drF(k)
381 ELSE
382 div(i,j) = div(i,j)*_hFacC(i,j,k,bi,bj)*drF(k)
383 ENDIF
384 ENDDO
385 ENDDO
386 ENDIF
387 ELSE
388 #else /* NONLIN_FRSURF */
389 IF ( .TRUE. ) THEN
390 #endif /* NONLIN_FRSURF */
391 C-- use current hFac (consistent with solver matrix)
392 DO j=jMin,jMax
393 DO i=iMin,iMax
394 div(i,j) = div(i,j)*_hFacC(i,j,k,bi,bj)*drF(k)
395 ENDDO
396 ENDDO
397 ENDIF
398 DO j=jMin,jMax
399 DO i=iMin,iMax
400 flx(i,j) = flx(i,j)
401 #ifdef ALLOW_AUTODIFF_OPENAD
402 & - rkSign*div(i,j)*rA(i,j,bi,bj)%v/deltaTLev(k)
403 #else
404 & - rkSign*div(i,j)*rA(i,j,bi,bj)/deltaTLev(k)
405 #endif
406 af(i,j) = flx(i,j) - df(i,j)
407 ENDDO
408 ENDDO
409 diagName = 'ADVr'//diagSufx
410 CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
411 ENDIF
412 ENDDO
413 ENDIF
414 ENDIF
415 #endif /* ALLOW_DIAGNOSTICS */
416
417 C-- end if Nr > 1
418 ENDIF
419
420 RETURN
421 END

  ViewVC Help
Powered by ViewVC 1.1.22