1 |
|
2 |
! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/ppm3d.F,v 1.1.1.1 2004/05/19 |
3 |
! 12:53:07 lmdzadmin Exp $ |
4 |
|
5 |
|
6 |
! From lin@explorer.gsfc.nasa.gov Wed Apr 15 17:44:44 1998 |
7 |
! Date: Wed, 15 Apr 1998 11:37:03 -0400 |
8 |
! From: lin@explorer.gsfc.nasa.gov |
9 |
! To: Frederic.Hourdin@lmd.jussieu.fr |
10 |
! Subject: 3D transport module of the GSFC CTM and GEOS GCM |
11 |
|
12 |
|
13 |
! This code is sent to you by S-J Lin, DAO, NASA-GSFC |
14 |
|
15 |
! Note: this version is intended for machines like CRAY |
16 |
! -90. No multitasking directives implemented. |
17 |
|
18 |
|
19 |
! ******************************************************************** |
20 |
|
21 |
! TransPort Core for Goddard Chemistry Transport Model (G-CTM), Goddard |
22 |
! Earth Observing System General Circulation Model (GEOS-GCM), and Data |
23 |
! Assimilation System (GEOS-DAS). |
24 |
|
25 |
! ******************************************************************** |
26 |
|
27 |
! Purpose: given horizontal winds on a hybrid sigma-p surfaces, |
28 |
! one call to tpcore updates the 3-D mixing ratio |
29 |
! fields one time step (NDT). [vertical mass flux is computed |
30 |
! internally consistent with the discretized hydrostatic mass |
31 |
! continuity equation of the C-Grid GEOS-GCM (for IGD=1)]. |
32 |
|
33 |
! Schemes: Multi-dimensional Flux Form Semi-Lagrangian (FFSL) scheme based |
34 |
! on the van Leer or PPM. |
35 |
! (see Lin and Rood 1996). |
36 |
! Version 4.5 |
37 |
! Last modified: Dec. 5, 1996 |
38 |
! Major changes from version 4.0: a more general vertical hybrid sigma- |
39 |
! pressure coordinate. |
40 |
! Subroutines modified: xtp, ytp, fzppm, qckxyz |
41 |
! Subroutines deleted: vanz |
42 |
|
43 |
! Author: Shian-Jiann Lin |
44 |
! mail address: |
45 |
! Shian-Jiann Lin* |
46 |
! Code 910.3, NASA/GSFC, Greenbelt, MD 20771 |
47 |
! Phone: 301-286-9540 |
48 |
! E-mail: lin@dao.gsfc.nasa.gov |
49 |
|
50 |
! *affiliation: |
51 |
! Joint Center for Earth Systems Technology |
52 |
! The University of Maryland Baltimore County |
53 |
! NASA - Goddard Space Flight Center |
54 |
! References: |
55 |
|
56 |
! 1. Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux form semi- |
57 |
! Lagrangian transport schemes. Mon. Wea. Rev., 124, 2046-2070. |
58 |
|
59 |
! 2. Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of |
60 |
! the van Leer-type transport schemes and its applications to the moist- |
61 |
! ure transport in a General Circulation Model. Mon. Wea. Rev., 122, |
62 |
! 1575-1593. |
63 |
|
64 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
65 |
|
66 |
SUBROUTINE ppm3d(igd, q, ps1, ps2, u, v, w, ndt, iord, jord, kord, nc, imr, & |
67 |
jnp, j1, nlay, ap, bp, pt, ae, fill, dum, umax) |
68 |
|
69 |
! implicit none |
70 |
|
71 |
! rajout de déclarations |
72 |
! integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd |
73 |
! integer iu,iiu,j2,jmr,js0,jt |
74 |
! real dtdy,dtdy5,rcap,iml,jn0,imjm,pi,dl,dp |
75 |
! real dt,cr1,maxdt,ztc,d5,sum1,sum2,ru |
76 |
|
77 |
! ******************************************************************** |
78 |
|
79 |
! ============= |
80 |
! INPUT: |
81 |
! ============= |
82 |
|
83 |
! Q(IMR,JNP,NLAY,NC): mixing ratios at current time (t) |
84 |
! NC: total number of constituents |
85 |
! IMR: first dimension (E-W); number of Grid intervals in E-W is IMR |
86 |
! JNP: 2nd dimension (N-S); number of Grid intervals in N-S is JNP-1 |
87 |
! NLAY: 3rd dimension (number of layers); vertical index increases from 1 |
88 |
! at |
89 |
! the model top to NLAY near the surface (see fig. below). |
90 |
! It is assumed that 6 <= NLAY <= JNP (for dynamic memory allocation) |
91 |
|
92 |
! PS1(IMR,JNP): surface pressure at current time (t) |
93 |
! PS2(IMR,JNP): surface pressure at mid-time-level (t+NDT/2) |
94 |
! PS2 is replaced by the predicted PS (at t+NDT) on output. |
95 |
! Note: surface pressure can have any unit or can be multiplied by any |
96 |
! const. |
97 |
|
98 |
! The pressure at layer edges are defined as follows: |
99 |
|
100 |
! p(i,j,k) = AP(k)*PT + BP(k)*PS(i,j) (1) |
101 |
|
102 |
! Where PT is a constant having the same unit as PS. |
103 |
! AP and BP are unitless constants given at layer edges |
104 |
! defining the vertical coordinate. |
105 |
! BP(1) = 0., BP(NLAY+1) = 1. |
106 |
! The pressure at the model top is PTOP = AP(1)*PT |
107 |
|
108 |
! For pure sigma system set AP(k) = 1 for all k, PT = PTOP, |
109 |
! BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP. |
110 |
|
111 |
! Note: the sigma-P coordinate is a subset of Eq. 1, which in turn |
112 |
! is a subset of the following even more general sigma-P-thelta coord. |
113 |
! currently under development. |
114 |
! p(i,j,k) = (AP(k)*PT + BP(k)*PS(i,j))/(D(k)-C(k)*TE**(-1/kapa)) |
115 |
|
116 |
! ///////////////////////////////// |
117 |
! / \ ------------- PTOP -------------- AP(1), BP(1) |
118 |
! | |
119 |
! delp(1) | ........... Q(i,j,1) ............ |
120 |
! | |
121 |
! W(1) \ / --------------------------------- AP(2), BP(2) |
122 |
|
123 |
|
124 |
|
125 |
! W(k-1) / \ --------------------------------- AP(k), BP(k) |
126 |
! | |
127 |
! delp(K) | ........... Q(i,j,k) ............ |
128 |
! | |
129 |
! W(k) \ / --------------------------------- AP(k+1), BP(k+1) |
130 |
|
131 |
|
132 |
|
133 |
! / \ --------------------------------- AP(NLAY), BP(NLAY) |
134 |
! | |
135 |
! delp(NLAY) | ........... Q(i,j,NLAY) ......... |
136 |
! | |
137 |
! W(NLAY)=0 \ / ------------- surface ----------- AP(NLAY+1), BP(NLAY+1) |
138 |
! ////////////////////////////////// |
139 |
|
140 |
! U(IMR,JNP,NLAY) & V(IMR,JNP,NLAY):winds (m/s) at mid-time-level (t+NDT/2) |
141 |
! U and V may need to be polar filtered in advance in some cases. |
142 |
|
143 |
! IGD: grid type on which winds are defined. |
144 |
! IGD = 0: A-Grid [all variables defined at the same point from south |
145 |
! pole (j=1) to north pole (j=JNP) ] |
146 |
|
147 |
! IGD = 1 GEOS-GCM C-Grid |
148 |
! [North] |
149 |
|
150 |
! V(i,j) |
151 |
! | |
152 |
! | |
153 |
! | |
154 |
! U(i-1,j)---Q(i,j)---U(i,j) [EAST] |
155 |
! | |
156 |
! | |
157 |
! | |
158 |
! V(i,j-1) |
159 |
|
160 |
! U(i, 1) is defined at South Pole. |
161 |
! V(i, 1) is half grid north of the South Pole. |
162 |
! V(i,JMR) is half grid south of the North Pole. |
163 |
|
164 |
! V must be defined at j=1 and j=JMR if IGD=1 |
165 |
! V at JNP need not be given. |
166 |
|
167 |
! NDT: time step in seconds (need not be constant during the course of |
168 |
! the integration). Suggested value: 30 min. for 4x5, 15 min. for 2x2.5 |
169 |
! (Lat-Lon) resolution. Smaller values are recommanded if the model |
170 |
! has a well-resolved stratosphere. |
171 |
|
172 |
! J1 defines the size of the polar cap: |
173 |
! South polar cap edge is located at -90 + (j1-1.5)*180/(JNP-1) deg. |
174 |
! North polar cap edge is located at 90 - (j1-1.5)*180/(JNP-1) deg. |
175 |
! There are currently only two choices (j1=2 or 3). |
176 |
! IMR must be an even integer if j1 = 2. Recommended value: J1=3. |
177 |
|
178 |
! IORD, JORD, and KORD are integers controlling various options in E-W, |
179 |
! N-S, |
180 |
! and vertical transport, respectively. Recommended values for positive |
181 |
! definite scalars: IORD=JORD=3, KORD=5. Use KORD=3 for non- |
182 |
! positive definite scalars or when linear correlation between constituents |
183 |
! is to be maintained. |
184 |
|
185 |
! _ORD= |
186 |
! 1: 1st order upstream scheme (too diffusive, not a useful option; it |
187 |
! can be used for debugging purposes; this is THE only known "linear" |
188 |
! monotonic advection scheme.). |
189 |
! 2: 2nd order van Leer (full monotonicity constraint; |
190 |
! see Lin et al 1994, MWR) |
191 |
! 3: monotonic PPM* (slightly improved PPM of Collela & Woodward 1984) |
192 |
! 4: semi-monotonic PPM (same as 3, but overshoots are allowed) |
193 |
! 5: positive-definite PPM (constraint on the subgrid distribution is |
194 |
! only strong enough to prevent generation of negative values; |
195 |
! both overshoots & undershoots are possible). |
196 |
! 6: un-constrained PPM (nearly diffusion free; slightly faster but |
197 |
! positivity not quaranteed. Use this option only when the fields |
198 |
! and winds are very smooth). |
199 |
|
200 |
! *PPM: Piece-wise Parabolic Method |
201 |
|
202 |
! Note that KORD <=2 options are no longer supported. DO not use option 4 |
203 |
! or 5. |
204 |
! for non-positive definite scalars (such as Ertel Potential Vorticity). |
205 |
|
206 |
! The implicit numerical diffusion decreases as _ORD increases. |
207 |
! The last two options (ORDER=5, 6) should only be used when there is |
208 |
! significant explicit diffusion (such as a turbulence parameterization). |
209 |
! You |
210 |
! might get dispersive results otherwise. |
211 |
! No filter of any kind is applied to the constituent fields here. |
212 |
|
213 |
! AE: Radius of the sphere (meters). |
214 |
! Recommended value for the planet earth: 6.371E6 |
215 |
|
216 |
! fill(logical): flag to do filling for negatives (see note below). |
217 |
|
218 |
! Umax: Estimate (upper limit) of the maximum U-wind speed (m/s). |
219 |
! (220 m/s is a good value for troposphere model; 280 m/s otherwise) |
220 |
|
221 |
! ============= |
222 |
! Output |
223 |
! ============= |
224 |
|
225 |
! Q: mixing ratios at future time (t+NDT) (original values are |
226 |
! over-written) |
227 |
! W(NLAY): large-scale vertical mass flux as diagnosed from the hydrostatic |
228 |
! relationship. W will have the same unit as PS1 and PS2 (eg, mb). |
229 |
! W must be divided by NDT to get the correct mass-flux unit. |
230 |
! The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND |
231 |
! is the pressure thickness in the "upwind" direction. For example, |
232 |
! C(k) = W(k)/delp(k) if W(k) > 0; |
233 |
! C(k) = W(k)/delp(k+1) if W(k) < 0. |
234 |
! ( W > 0 is downward, ie, toward surface) |
235 |
! PS2: predicted PS at t+NDT (original values are over-written) |
236 |
|
237 |
! ******************************************************************** |
238 |
! NOTES: |
239 |
! This forward-in-time upstream-biased transport scheme reduces to |
240 |
! the 2nd order center-in-time center-in-space mass continuity eqn. |
241 |
! if Q = 1 (constant fields will remain constant). This also ensures |
242 |
! that the computed vertical velocity to be identical to GEOS-1 GCM |
243 |
! for on-line transport. |
244 |
|
245 |
! A larger polar cap is used if j1=3 (recommended for C-Grid winds or when |
246 |
! winds are noisy near poles). |
247 |
|
248 |
! Flux-Form Semi-Lagrangian transport in the East-West direction is used |
249 |
! when and where Courant number is greater than one. |
250 |
|
251 |
! The user needs to change the parameter Jmax or Kmax if the resolution |
252 |
! is greater than 0.5 deg in N-S or 150 layers in the vertical direction. |
253 |
! (this TransPort Core is otherwise resolution independent and can be used |
254 |
! as a library routine). |
255 |
|
256 |
! PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd |
257 |
! order accurate for non-uniform grid (vertical sigma coord.). |
258 |
|
259 |
! Time step is limitted only by transport in the meridional direction. |
260 |
! (the FFSL scheme is not implemented in the meridional direction). |
261 |
|
262 |
! Since only 1-D limiters are applied, negative values could |
263 |
! potentially be generated when large time step is used and when the |
264 |
! initial fields contain discontinuities. |
265 |
! This does not necessarily imply the integration is unstable. |
266 |
! These negatives are typically very small. A filling algorithm is |
267 |
! activated if the user set "fill" to be true. |
268 |
|
269 |
! The van Leer scheme used here is nearly as accurate as the original PPM |
270 |
! due to the use of a 4th order accurate reference slope. The PPM imple- |
271 |
! mented here is an improvement over the original and is also based on |
272 |
! the 4th order reference slope. |
273 |
|
274 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
275 |
|
276 |
! User modifiable parameters |
277 |
|
278 |
PARAMETER (jmax=361, kmax=150) |
279 |
|
280 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
281 |
|
282 |
! Input-Output arrays |
283 |
|
284 |
|
285 |
REAL q(imr, jnp, nlay, nc), ps1(imr, jnp), ps2(imr, jnp), & |
286 |
u(imr, jnp, nlay), v(imr, jnp, nlay), ap(nlay+1), bp(nlay+1), & |
287 |
w(imr, jnp, nlay), ndt, val(nlay), umax |
288 |
INTEGER igd, iord, jord, kord, nc, imr, jnp, j1, nlay, ae |
289 |
INTEGER imrd2 |
290 |
REAL pt |
291 |
LOGICAL cross, fill, dum |
292 |
|
293 |
! Local dynamic arrays |
294 |
|
295 |
REAL crx(imr, jnp), cry(imr, jnp), xmass(imr, jnp), ymass(imr, jnp), & |
296 |
fx1(imr+1), dpi(imr, jnp, nlay), delp1(imr, jnp, nlay), & |
297 |
wk1(imr, jnp, nlay), pu(imr, jnp), pv(imr, jnp), dc2(imr, jnp), & |
298 |
delp2(imr, jnp, nlay), dq(imr, jnp, nlay, nc), va(imr, jnp), & |
299 |
ua(imr, jnp), qtmp(-imr:2*imr) |
300 |
|
301 |
! Local static arrays |
302 |
|
303 |
REAL dtdx(jmax), dtdx5(jmax), acosp(jmax), cosp(jmax), cose(jmax), & |
304 |
dap(kmax), dbk(kmax) |
305 |
DATA ndt0, nstep/0, 0/ |
306 |
DATA cross/.TRUE./ |
307 |
SAVE dtdy, dtdy5, rcap, js0, jn0, iml, dtdx, dtdx5, acosp, cosp, cose, dap, & |
308 |
dbk |
309 |
|
310 |
|
311 |
jmr = jnp - 1 |
312 |
imjm = imr*jnp |
313 |
j2 = jnp - j1 + 1 |
314 |
nstep = nstep + 1 |
315 |
|
316 |
! *********** Initialization ********************** |
317 |
IF (nstep==1) THEN |
318 |
|
319 |
WRITE (6, *) '------------------------------------ ' |
320 |
WRITE (6, *) 'NASA/GSFC Transport Core Version 4.5' |
321 |
WRITE (6, *) '------------------------------------ ' |
322 |
|
323 |
WRITE (6, *) 'IMR=', imr, ' JNP=', jnp, ' NLAY=', nlay, ' j1=', j1 |
324 |
WRITE (6, *) 'NC=', nc, iord, jord, kord, ndt |
325 |
|
326 |
! controles sur les parametres |
327 |
IF (nlay<6) THEN |
328 |
WRITE (6, *) 'NLAY must be >= 6' |
329 |
STOP |
330 |
END IF |
331 |
IF (jnp<nlay) THEN |
332 |
WRITE (6, *) 'JNP must be >= NLAY' |
333 |
STOP |
334 |
END IF |
335 |
imrd2 = mod(imr, 2) |
336 |
IF (j1==2 .AND. imrd2/=0) THEN |
337 |
WRITE (6, *) 'if j1=2 IMR must be an even integer' |
338 |
STOP |
339 |
END IF |
340 |
|
341 |
|
342 |
IF (jmax<jnp .OR. kmax<nlay) THEN |
343 |
WRITE (6, *) 'Jmax or Kmax is too small' |
344 |
STOP |
345 |
END IF |
346 |
|
347 |
DO k = 1, nlay |
348 |
dap(k) = (ap(k+1)-ap(k))*pt |
349 |
dbk(k) = bp(k+1) - bp(k) |
350 |
END DO |
351 |
|
352 |
pi = 4.*atan(1.) |
353 |
dl = 2.*pi/float(imr) |
354 |
dp = pi/float(jmr) |
355 |
|
356 |
IF (igd==0) THEN |
357 |
! Compute analytic cosine at cell edges |
358 |
CALL cosa(cosp, cose, jnp, pi, dp) |
359 |
ELSE |
360 |
! Define cosine consistent with GEOS-GCM (using dycore2.0 or later) |
361 |
CALL cosc(cosp, cose, jnp, pi, dp) |
362 |
END IF |
363 |
|
364 |
DO j = 2, jmr |
365 |
acosp(j) = 1./cosp(j) |
366 |
END DO |
367 |
|
368 |
! Inverse of the Scaled polar cap area. |
369 |
|
370 |
rcap = dp/(imr*(1.-cos((j1-1.5)*dp))) |
371 |
acosp(1) = rcap |
372 |
acosp(jnp) = rcap |
373 |
END IF |
374 |
|
375 |
IF (ndt0/=ndt) THEN |
376 |
dt = ndt |
377 |
ndt0 = ndt |
378 |
|
379 |
IF (umax<180.) THEN |
380 |
WRITE (6, *) 'Umax may be too small!' |
381 |
END IF |
382 |
cr1 = abs(umax*dt)/(dl*ae) |
383 |
maxdt = dp*ae/abs(umax) + 0.5 |
384 |
WRITE (6, *) 'Largest time step for max(V)=', umax, ' is ', maxdt |
385 |
IF (maxdt<abs(ndt)) THEN |
386 |
WRITE (6, *) 'Warning!!! NDT maybe too large!' |
387 |
END IF |
388 |
|
389 |
IF (cr1>=0.95) THEN |
390 |
js0 = 0 |
391 |
jn0 = 0 |
392 |
iml = imr - 2 |
393 |
ztc = 0. |
394 |
ELSE |
395 |
ztc = acos(cr1)*(180./pi) |
396 |
|
397 |
js0 = float(jmr)*(90.-ztc)/180. + 2 |
398 |
js0 = max(js0, j1+1) |
399 |
iml = min(6*js0/(j1-1)+2, 4*imr/5) |
400 |
jn0 = jnp - js0 + 1 |
401 |
END IF |
402 |
|
403 |
|
404 |
DO j = 2, jmr |
405 |
dtdx(j) = dt/(dl*ae*cosp(j)) |
406 |
|
407 |
dtdx5(j) = 0.5*dtdx(j) |
408 |
END DO |
409 |
|
410 |
|
411 |
dtdy = dt/(ae*dp) |
412 |
dtdy5 = 0.5*dtdy |
413 |
|
414 |
END IF |
415 |
|
416 |
! *********** End Initialization ********************** |
417 |
|
418 |
! delp = pressure thickness: the psudo-density in a hydrostatic system. |
419 |
DO k = 1, nlay |
420 |
DO j = 1, jnp |
421 |
DO i = 1, imr |
422 |
delp1(i, j, k) = dap(k) + dbk(k)*ps1(i, j) |
423 |
delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j) |
424 |
END DO |
425 |
END DO |
426 |
END DO |
427 |
|
428 |
|
429 |
IF (j1/=2) THEN |
430 |
DO ic = 1, nc |
431 |
DO l = 1, nlay |
432 |
DO i = 1, imr |
433 |
q(i, 2, l, ic) = q(i, 1, l, ic) |
434 |
q(i, jmr, l, ic) = q(i, jnp, l, ic) |
435 |
END DO |
436 |
END DO |
437 |
END DO |
438 |
END IF |
439 |
|
440 |
! Compute "tracer density" |
441 |
DO ic = 1, nc |
442 |
DO k = 1, nlay |
443 |
DO j = 1, jnp |
444 |
DO i = 1, imr |
445 |
dq(i, j, k, ic) = q(i, j, k, ic)*delp1(i, j, k) |
446 |
END DO |
447 |
END DO |
448 |
END DO |
449 |
END DO |
450 |
|
451 |
DO k = 1, nlay |
452 |
|
453 |
IF (igd==0) THEN |
454 |
! Convert winds on A-Grid to Courant number on C-Grid. |
455 |
CALL a2c(u(1,1,k), v(1,1,k), imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5) |
456 |
ELSE |
457 |
! Convert winds on C-grid to Courant number |
458 |
DO j = j1, j2 |
459 |
DO i = 2, imr |
460 |
crx(i, j) = dtdx(j)*u(i-1, j, k) |
461 |
END DO |
462 |
END DO |
463 |
|
464 |
|
465 |
DO j = j1, j2 |
466 |
crx(1, j) = dtdx(j)*u(imr, j, k) |
467 |
END DO |
468 |
|
469 |
DO i = 1, imr*jmr |
470 |
cry(i, 2) = dtdy*v(i, 1, k) |
471 |
END DO |
472 |
END IF |
473 |
|
474 |
! Determine JS and JN |
475 |
js = j1 |
476 |
jn = j2 |
477 |
|
478 |
DO j = js0, j1 + 1, -1 |
479 |
DO i = 1, imr |
480 |
IF (abs(crx(i,j))>1.) THEN |
481 |
js = j |
482 |
GO TO 2222 |
483 |
END IF |
484 |
END DO |
485 |
END DO |
486 |
|
487 |
2222 CONTINUE |
488 |
DO j = jn0, j2 - 1 |
489 |
DO i = 1, imr |
490 |
IF (abs(crx(i,j))>1.) THEN |
491 |
jn = j |
492 |
GO TO 2233 |
493 |
END IF |
494 |
END DO |
495 |
END DO |
496 |
2233 CONTINUE |
497 |
|
498 |
IF (j1/=2) THEN ! Enlarged polar cap. |
499 |
DO i = 1, imr |
500 |
dpi(i, 2, k) = 0. |
501 |
dpi(i, jmr, k) = 0. |
502 |
END DO |
503 |
END IF |
504 |
|
505 |
! ******* Compute horizontal mass fluxes ************ |
506 |
|
507 |
! N-S component |
508 |
DO j = j1, j2 + 1 |
509 |
d5 = 0.5*cose(j) |
510 |
DO i = 1, imr |
511 |
ymass(i, j) = cry(i, j)*d5*(delp2(i,j,k)+delp2(i,j-1,k)) |
512 |
END DO |
513 |
END DO |
514 |
|
515 |
DO j = j1, j2 |
516 |
DO i = 1, imr |
517 |
dpi(i, j, k) = (ymass(i,j)-ymass(i,j+1))*acosp(j) |
518 |
END DO |
519 |
END DO |
520 |
|
521 |
! Poles |
522 |
sum1 = ymass(imr, j1) |
523 |
sum2 = ymass(imr, j2+1) |
524 |
DO i = 1, imr - 1 |
525 |
sum1 = sum1 + ymass(i, j1) |
526 |
sum2 = sum2 + ymass(i, j2+1) |
527 |
END DO |
528 |
|
529 |
sum1 = -sum1*rcap |
530 |
sum2 = sum2*rcap |
531 |
DO i = 1, imr |
532 |
dpi(i, 1, k) = sum1 |
533 |
dpi(i, jnp, k) = sum2 |
534 |
END DO |
535 |
|
536 |
! E-W component |
537 |
|
538 |
DO j = j1, j2 |
539 |
DO i = 2, imr |
540 |
pu(i, j) = 0.5*(delp2(i,j,k)+delp2(i-1,j,k)) |
541 |
END DO |
542 |
END DO |
543 |
|
544 |
DO j = j1, j2 |
545 |
pu(1, j) = 0.5*(delp2(1,j,k)+delp2(imr,j,k)) |
546 |
END DO |
547 |
|
548 |
DO j = j1, j2 |
549 |
DO i = 1, imr |
550 |
xmass(i, j) = pu(i, j)*crx(i, j) |
551 |
END DO |
552 |
END DO |
553 |
|
554 |
DO j = j1, j2 |
555 |
DO i = 1, imr - 1 |
556 |
dpi(i, j, k) = dpi(i, j, k) + xmass(i, j) - xmass(i+1, j) |
557 |
END DO |
558 |
END DO |
559 |
|
560 |
DO j = j1, j2 |
561 |
dpi(imr, j, k) = dpi(imr, j, k) + xmass(imr, j) - xmass(1, j) |
562 |
END DO |
563 |
|
564 |
DO j = j1, j2 |
565 |
DO i = 1, imr - 1 |
566 |
ua(i, j) = 0.5*(crx(i,j)+crx(i+1,j)) |
567 |
END DO |
568 |
END DO |
569 |
|
570 |
DO j = j1, j2 |
571 |
ua(imr, j) = 0.5*(crx(imr,j)+crx(1,j)) |
572 |
END DO |
573 |
! cccccccccccccccccccccccccccccccccccccccccccccccccccccc |
574 |
! Rajouts pour LMDZ.3.3 |
575 |
! cccccccccccccccccccccccccccccccccccccccccccccccccccccc |
576 |
DO i = 1, imr |
577 |
DO j = 1, jnp |
578 |
va(i, j) = 0. |
579 |
END DO |
580 |
END DO |
581 |
|
582 |
DO i = 1, imr*(jmr-1) |
583 |
va(i, 2) = 0.5*(cry(i,2)+cry(i,3)) |
584 |
END DO |
585 |
|
586 |
IF (j1==2) THEN |
587 |
imh = imr/2 |
588 |
DO i = 1, imh |
589 |
va(i, 1) = 0.5*(cry(i,2)-cry(i+imh,2)) |
590 |
va(i+imh, 1) = -va(i, 1) |
591 |
va(i, jnp) = 0.5*(cry(i,jnp)-cry(i+imh,jmr)) |
592 |
va(i+imh, jnp) = -va(i, jnp) |
593 |
END DO |
594 |
va(imr, 1) = va(1, 1) |
595 |
va(imr, jnp) = va(1, jnp) |
596 |
END IF |
597 |
|
598 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
599 |
DO ic = 1, nc |
600 |
|
601 |
DO i = 1, imjm |
602 |
wk1(i, 1, 1) = 0. |
603 |
wk1(i, 1, 2) = 0. |
604 |
END DO |
605 |
|
606 |
! E-W advective cross term |
607 |
DO j = j1, j2 |
608 |
IF (j>js .AND. j<jn) GO TO 250 |
609 |
|
610 |
DO i = 1, imr |
611 |
qtmp(i) = q(i, j, k, ic) |
612 |
END DO |
613 |
|
614 |
DO i = -iml, 0 |
615 |
qtmp(i) = q(imr+i, j, k, ic) |
616 |
qtmp(imr+1-i) = q(1-i, j, k, ic) |
617 |
END DO |
618 |
|
619 |
DO i = 1, imr |
620 |
iu = ua(i, j) |
621 |
ru = ua(i, j) - iu |
622 |
iiu = i - iu |
623 |
IF (ua(i,j)>=0.) THEN |
624 |
wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu)) |
625 |
ELSE |
626 |
wk1(i, j, 1) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1)) |
627 |
END IF |
628 |
wk1(i, j, 1) = wk1(i, j, 1) - qtmp(i) |
629 |
END DO |
630 |
250 END DO |
631 |
|
632 |
IF (jn/=0) THEN |
633 |
DO j = js + 1, jn - 1 |
634 |
|
635 |
DO i = 1, imr |
636 |
qtmp(i) = q(i, j, k, ic) |
637 |
END DO |
638 |
|
639 |
qtmp(0) = q(imr, j, k, ic) |
640 |
qtmp(imr+1) = q(1, j, k, ic) |
641 |
|
642 |
DO i = 1, imr |
643 |
iu = i - ua(i, j) |
644 |
wk1(i, j, 1) = ua(i, j)*(qtmp(iu)-qtmp(iu+1)) |
645 |
END DO |
646 |
END DO |
647 |
END IF |
648 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
649 |
! Contribution from the N-S advection |
650 |
DO i = 1, imr*(j2-j1+1) |
651 |
jt = float(j1) - va(i, j1) |
652 |
wk1(i, j1, 2) = va(i, j1)*(q(i,jt,k,ic)-q(i,jt+1,k,ic)) |
653 |
END DO |
654 |
|
655 |
DO i = 1, imjm |
656 |
wk1(i, 1, 1) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 1) |
657 |
wk1(i, 1, 2) = q(i, 1, k, ic) + 0.5*wk1(i, 1, 2) |
658 |
END DO |
659 |
|
660 |
IF (cross) THEN |
661 |
! Add cross terms in the vertical direction. |
662 |
IF (iord>=2) THEN |
663 |
iad = 2 |
664 |
ELSE |
665 |
iad = 1 |
666 |
END IF |
667 |
|
668 |
IF (jord>=2) THEN |
669 |
jad = 2 |
670 |
ELSE |
671 |
jad = 1 |
672 |
END IF |
673 |
CALL xadv(imr, jnp, j1, j2, wk1(1,1,2), ua, js, jn, iml, dc2, iad) |
674 |
CALL yadv(imr, jnp, j1, j2, wk1(1,1,1), va, pv, w, jad) |
675 |
DO j = 1, jnp |
676 |
DO i = 1, imr |
677 |
q(i, j, k, ic) = q(i, j, k, ic) + dc2(i, j) + pv(i, j) |
678 |
END DO |
679 |
END DO |
680 |
END IF |
681 |
|
682 |
CALL xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq(1,1,k,ic), wk1(1,1,2), & |
683 |
crx, fx1, xmass, iord) |
684 |
|
685 |
CALL ytp(imr, jnp, j1, j2, acosp, rcap, dq(1,1,k,ic), wk1(1,1,1), cry, & |
686 |
dc2, ymass, wk1(1,1,3), wk1(1,1,4), wk1(1,1,5), wk1(1,1,6), jord) |
687 |
|
688 |
END DO |
689 |
END DO |
690 |
|
691 |
! ******* Compute vertical mass flux (same unit as PS) *********** |
692 |
|
693 |
! 1st step: compute total column mass CONVERGENCE. |
694 |
|
695 |
DO j = 1, jnp |
696 |
DO i = 1, imr |
697 |
cry(i, j) = dpi(i, j, 1) |
698 |
END DO |
699 |
END DO |
700 |
|
701 |
DO k = 2, nlay |
702 |
DO j = 1, jnp |
703 |
DO i = 1, imr |
704 |
cry(i, j) = cry(i, j) + dpi(i, j, k) |
705 |
END DO |
706 |
END DO |
707 |
END DO |
708 |
|
709 |
DO j = 1, jnp |
710 |
DO i = 1, imr |
711 |
|
712 |
! 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption. |
713 |
! Changes (increases) to surface pressure = total column mass |
714 |
! convergence |
715 |
|
716 |
ps2(i, j) = ps1(i, j) + cry(i, j) |
717 |
|
718 |
! 3rd step: compute vertical mass flux from mass conservation |
719 |
! principle. |
720 |
|
721 |
w(i, j, 1) = dpi(i, j, 1) - dbk(1)*cry(i, j) |
722 |
w(i, j, nlay) = 0. |
723 |
END DO |
724 |
END DO |
725 |
|
726 |
DO k = 2, nlay - 1 |
727 |
DO j = 1, jnp |
728 |
DO i = 1, imr |
729 |
w(i, j, k) = w(i, j, k-1) + dpi(i, j, k) - dbk(k)*cry(i, j) |
730 |
END DO |
731 |
END DO |
732 |
END DO |
733 |
|
734 |
DO k = 1, nlay |
735 |
DO j = 1, jnp |
736 |
DO i = 1, imr |
737 |
delp2(i, j, k) = dap(k) + dbk(k)*ps2(i, j) |
738 |
END DO |
739 |
END DO |
740 |
END DO |
741 |
|
742 |
krd = max(3, kord) |
743 |
DO ic = 1, nc |
744 |
|
745 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
746 |
|
747 |
CALL fzppm(imr, jnp, nlay, j1, dq(1,1,1,ic), w, q(1,1,1,ic), wk1, dpi, & |
748 |
dc2, crx, cry, pu, pv, xmass, ymass, delp1, krd) |
749 |
|
750 |
|
751 |
IF (fill) CALL qckxyz(dq(1,1,1,ic), dc2, imr, jnp, nlay, j1, j2, cosp, & |
752 |
acosp, .FALSE., ic, nstep) |
753 |
|
754 |
! Recover tracer mixing ratio from "density" using predicted |
755 |
! "air density" (pressure thickness) at time-level n+1 |
756 |
|
757 |
DO k = 1, nlay |
758 |
DO j = 1, jnp |
759 |
DO i = 1, imr |
760 |
q(i, j, k, ic) = dq(i, j, k, ic)/delp2(i, j, k) |
761 |
END DO |
762 |
END DO |
763 |
END DO |
764 |
|
765 |
IF (j1/=2) THEN |
766 |
DO k = 1, nlay |
767 |
DO i = 1, imr |
768 |
! j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord |
769 |
q(i, 2, k, ic) = q(i, 1, k, ic) |
770 |
q(i, jmr, k, ic) = q(i, jmp, k, ic) |
771 |
END DO |
772 |
END DO |
773 |
END IF |
774 |
END DO |
775 |
|
776 |
IF (j1/=2) THEN |
777 |
DO k = 1, nlay |
778 |
DO i = 1, imr |
779 |
w(i, 2, k) = w(i, 1, k) |
780 |
w(i, jmr, k) = w(i, jnp, k) |
781 |
END DO |
782 |
END DO |
783 |
END IF |
784 |
|
785 |
RETURN |
786 |
END SUBROUTINE ppm3d |
787 |
|
788 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
789 |
SUBROUTINE fzppm(imr, jnp, nlay, j1, dq, wz, p, dc, dqdt, ar, al, a6, flux, & |
790 |
wk1, wk2, wz2, delp, kord) |
791 |
PARAMETER (kmax=150) |
792 |
PARAMETER (r23=2./3., r3=1./3.) |
793 |
REAL wz(imr, jnp, nlay), p(imr, jnp, nlay), dc(imr, jnp, nlay), & |
794 |
wk1(imr, *), delp(imr, jnp, nlay), dq(imr, jnp, nlay), & |
795 |
dqdt(imr, jnp, nlay) |
796 |
! Assuming JNP >= NLAY |
797 |
REAL ar(imr, *), al(imr, *), a6(imr, *), flux(imr, *), wk2(imr, *), & |
798 |
wz2(imr, *) |
799 |
|
800 |
jmr = jnp - 1 |
801 |
imjm = imr*jnp |
802 |
nlaym1 = nlay - 1 |
803 |
|
804 |
lmt = kord - 3 |
805 |
|
806 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
807 |
! Compute DC for PPM |
808 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
809 |
|
810 |
DO k = 1, nlaym1 |
811 |
DO i = 1, imjm |
812 |
dqdt(i, 1, k) = p(i, 1, k+1) - p(i, 1, k) |
813 |
END DO |
814 |
END DO |
815 |
|
816 |
DO k = 2, nlaym1 |
817 |
DO i = 1, imjm |
818 |
c0 = delp(i, 1, k)/(delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1)) |
819 |
c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k)) |
820 |
c2 = (delp(i,1,k+1)+0.5*delp(i,1,k))/(delp(i,1,k-1)+delp(i,1,k)) |
821 |
tmp = c0*(c1*dqdt(i,1,k)+c2*dqdt(i,1,k-1)) |
822 |
qmax = max(p(i,1,k-1), p(i,1,k), p(i,1,k+1)) - p(i, 1, k) |
823 |
qmin = p(i, 1, k) - min(p(i,1,k-1), p(i,1,k), p(i,1,k+1)) |
824 |
dc(i, 1, k) = sign(min(abs(tmp),qmax,qmin), tmp) |
825 |
END DO |
826 |
END DO |
827 |
|
828 |
|
829 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
830 |
! Loop over latitudes (to save memory) |
831 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
832 |
|
833 |
DO j = 1, jnp |
834 |
IF ((j==2 .OR. j==jmr) .AND. j1/=2) GO TO 2000 |
835 |
|
836 |
DO k = 1, nlay |
837 |
DO i = 1, imr |
838 |
wz2(i, k) = wz(i, j, k) |
839 |
wk1(i, k) = p(i, j, k) |
840 |
wk2(i, k) = delp(i, j, k) |
841 |
flux(i, k) = dc(i, j, k) !this flux is actually the monotone slope |
842 |
END DO |
843 |
END DO |
844 |
|
845 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
846 |
! Compute first guesses at cell interfaces |
847 |
! First guesses are required to be continuous. |
848 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
849 |
|
850 |
! three-cell parabolic subgrid distribution at model top |
851 |
! two-cell parabolic with zero gradient subgrid distribution |
852 |
! at the surface. |
853 |
|
854 |
! First guess top edge value |
855 |
DO i = 1, imr |
856 |
! three-cell PPM |
857 |
! Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp |
858 |
a = 3.*(dqdt(i,j,2)-dqdt(i,j,1)*(wk2(i,2)+wk2(i,3))/(wk2(i,1)+wk2(i, & |
859 |
2)))/((wk2(i,2)+wk2(i,3))*(wk2(i,1)+wk2(i,2)+wk2(i,3))) |
860 |
b = 2.*dqdt(i, j, 1)/(wk2(i,1)+wk2(i,2)) - r23*a*(2.*wk2(i,1)+wk2(i,2)) |
861 |
al(i, 1) = wk1(i, 1) - wk2(i, 1)*(r3*a*wk2(i,1)+0.5*b) |
862 |
al(i, 2) = wk2(i, 1)*(a*wk2(i,1)+b) + al(i, 1) |
863 |
|
864 |
! Check if change sign |
865 |
IF (wk1(i,1)*al(i,1)<=0.) THEN |
866 |
al(i, 1) = 0. |
867 |
flux(i, 1) = 0. |
868 |
ELSE |
869 |
flux(i, 1) = wk1(i, 1) - al(i, 1) |
870 |
END IF |
871 |
END DO |
872 |
|
873 |
! Bottom |
874 |
DO i = 1, imr |
875 |
! 2-cell PPM with zero gradient right at the surface |
876 |
|
877 |
fct = dqdt(i, j, nlaym1)*wk2(i, nlay)**2/((wk2(i,nlay)+wk2(i, & |
878 |
nlaym1))*(2.*wk2(i,nlay)+wk2(i,nlaym1))) |
879 |
ar(i, nlay) = wk1(i, nlay) + fct |
880 |
al(i, nlay) = wk1(i, nlay) - (fct+fct) |
881 |
IF (wk1(i,nlay)*ar(i,nlay)<=0.) ar(i, nlay) = 0. |
882 |
flux(i, nlay) = ar(i, nlay) - wk1(i, nlay) |
883 |
END DO |
884 |
|
885 |
|
886 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
887 |
! 4th order interpolation in the interior. |
888 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
889 |
|
890 |
DO k = 3, nlaym1 |
891 |
DO i = 1, imr |
892 |
c1 = dqdt(i, j, k-1)*wk2(i, k-1)/(wk2(i,k-1)+wk2(i,k)) |
893 |
c2 = 2./(wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1)) |
894 |
a1 = (wk2(i,k-2)+wk2(i,k-1))/(2.*wk2(i,k-1)+wk2(i,k)) |
895 |
a2 = (wk2(i,k)+wk2(i,k+1))/(2.*wk2(i,k)+wk2(i,k-1)) |
896 |
al(i, k) = wk1(i, k-1) + c1 + c2*(wk2(i,k)*(c1*(a1-a2)+a2*flux(i, & |
897 |
k-1))-wk2(i,k-1)*a1*flux(i,k)) |
898 |
END DO |
899 |
END DO |
900 |
|
901 |
DO i = 1, imr*nlaym1 |
902 |
ar(i, 1) = al(i, 2) |
903 |
END DO |
904 |
|
905 |
DO i = 1, imr*nlay |
906 |
a6(i, 1) = 3.*(wk1(i,1)+wk1(i,1)-(al(i,1)+ar(i,1))) |
907 |
END DO |
908 |
|
909 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
910 |
! Top & Bot always monotonic |
911 |
CALL lmtppm(flux(1,1), a6(1,1), ar(1,1), al(1,1), wk1(1,1), imr, 0) |
912 |
CALL lmtppm(flux(1,nlay), a6(1,nlay), ar(1,nlay), al(1,nlay), & |
913 |
wk1(1,nlay), imr, 0) |
914 |
|
915 |
! Interior depending on KORD |
916 |
IF (lmt<=2) CALL lmtppm(flux(1,2), a6(1,2), ar(1,2), al(1,2), wk1(1,2), & |
917 |
imr*(nlay-2), lmt) |
918 |
|
919 |
! ****6***0*********0*********0*********0*********0*********0**********72 |
920 |
|
921 |
DO i = 1, imr*nlaym1 |
922 |
IF (wz2(i,1)>0.) THEN |
923 |
cm = wz2(i, 1)/wk2(i, 1) |
924 |
flux(i, 2) = ar(i, 1) + 0.5*cm*(al(i,1)-ar(i,1)+a6(i,1)*(1.-r23*cm)) |
925 |
ELSE |
926 |
cp = wz2(i, 1)/wk2(i, 2) |
927 |
flux(i, 2) = al(i, 2) + 0.5*cp*(al(i,2)-ar(i,2)-a6(i,2)*(1.+r23*cp)) |
928 |
END IF |
929 |
END DO |
930 |
|
931 |
DO i = 1, imr*nlaym1 |
932 |
flux(i, 2) = wz2(i, 1)*flux(i, 2) |
933 |
END DO |
934 |
|
935 |
DO i = 1, imr |
936 |
dq(i, j, 1) = dq(i, j, 1) - flux(i, 2) |
937 |
dq(i, j, nlay) = dq(i, j, nlay) + flux(i, nlay) |
938 |
END DO |
939 |
|
940 |
DO k = 2, nlaym1 |
941 |
DO i = 1, imr |
942 |
dq(i, j, k) = dq(i, j, k) + flux(i, k) - flux(i, k+1) |
943 |
END DO |
944 |
END DO |
945 |
2000 END DO |
946 |
RETURN |
947 |
END SUBROUTINE fzppm |
948 |
|
949 |
SUBROUTINE xtp(imr, jnp, iml, j1, j2, jn, js, pu, dq, q, uc, fx1, xmass, & |
950 |
iord) |
951 |
DIMENSION uc(imr, *), dc(-iml:imr+iml+1), xmass(imr, jnp), fx1(imr+1), & |
952 |
dq(imr, jnp), qtmp(-iml:imr+1+iml) |
953 |
DIMENSION pu(imr, jnp), q(imr, jnp), isave(imr) |
954 |
|
955 |
imp = imr + 1 |
956 |
|
957 |
! van Leer at high latitudes |
958 |
jvan = max(1, jnp/18) |
959 |
j1vl = j1 + jvan |
960 |
j2vl = j2 - jvan |
961 |
|
962 |
DO j = j1, j2 |
963 |
|
964 |
DO i = 1, imr |
965 |
qtmp(i) = q(i, j) |
966 |
END DO |
967 |
|
968 |
IF (j>=jn .OR. j<=js) GO TO 2222 |
969 |
! ************* Eulerian ********** |
970 |
|
971 |
qtmp(0) = q(imr, j) |
972 |
qtmp(-1) = q(imr-1, j) |
973 |
qtmp(imp) = q(1, j) |
974 |
qtmp(imp+1) = q(2, j) |
975 |
|
976 |
IF (iord==1 .OR. j==j1 .OR. j==j2) THEN |
977 |
DO i = 1, imr |
978 |
iu = float(i) - uc(i, j) |
979 |
fx1(i) = qtmp(iu) |
980 |
END DO |
981 |
ELSE |
982 |
CALL xmist(imr, iml, qtmp, dc) |
983 |
dc(0) = dc(imr) |
984 |
|
985 |
IF (iord==2 .OR. j<=j1vl .OR. j>=j2vl) THEN |
986 |
DO i = 1, imr |
987 |
iu = float(i) - uc(i, j) |
988 |
fx1(i) = qtmp(iu) + dc(iu)*(sign(1.,uc(i,j))-uc(i,j)) |
989 |
END DO |
990 |
ELSE |
991 |
CALL fxppm(imr, iml, uc(1,j), qtmp, dc, fx1, iord) |
992 |
END IF |
993 |
|
994 |
END IF |
995 |
|
996 |
DO i = 1, imr |
997 |
fx1(i) = fx1(i)*xmass(i, j) |
998 |
END DO |
999 |
|
1000 |
GO TO 1309 |
1001 |
|
1002 |
! ***** Conservative (flux-form) Semi-Lagrangian transport ***** |
1003 |
|
1004 |
2222 CONTINUE |
1005 |
|
1006 |
DO i = -iml, 0 |
1007 |
qtmp(i) = q(imr+i, j) |
1008 |
qtmp(imp-i) = q(1-i, j) |
1009 |
END DO |
1010 |
|
1011 |
IF (iord==1 .OR. j==j1 .OR. j==j2) THEN |
1012 |
DO i = 1, imr |
1013 |
itmp = int(uc(i,j)) |
1014 |
isave(i) = i - itmp |
1015 |
iu = i - uc(i, j) |
1016 |
fx1(i) = (uc(i,j)-itmp)*qtmp(iu) |
1017 |
END DO |
1018 |
ELSE |
1019 |
CALL xmist(imr, iml, qtmp, dc) |
1020 |
|
1021 |
DO i = -iml, 0 |
1022 |
dc(i) = dc(imr+i) |
1023 |
dc(imp-i) = dc(1-i) |
1024 |
END DO |
1025 |
|
1026 |
DO i = 1, imr |
1027 |
itmp = int(uc(i,j)) |
1028 |
rut = uc(i, j) - itmp |
1029 |
isave(i) = i - itmp |
1030 |
iu = i - uc(i, j) |
1031 |
fx1(i) = rut*(qtmp(iu)+dc(iu)*(sign(1.,rut)-rut)) |
1032 |
END DO |
1033 |
END IF |
1034 |
|
1035 |
DO i = 1, imr |
1036 |
IF (uc(i,j)>1.) THEN |
1037 |
! DIR$ NOVECTOR |
1038 |
DO ist = isave(i), i - 1 |
1039 |
fx1(i) = fx1(i) + qtmp(ist) |
1040 |
END DO |
1041 |
ELSE IF (uc(i,j)<-1.) THEN |
1042 |
DO ist = i, isave(i) - 1 |
1043 |
fx1(i) = fx1(i) - qtmp(ist) |
1044 |
END DO |
1045 |
! DIR$ VECTOR |
1046 |
END IF |
1047 |
END DO |
1048 |
DO i = 1, imr |
1049 |
fx1(i) = pu(i, j)*fx1(i) |
1050 |
END DO |
1051 |
|
1052 |
! *************************************** |
1053 |
|
1054 |
1309 fx1(imp) = fx1(1) |
1055 |
DO i = 1, imr |
1056 |
dq(i, j) = dq(i, j) + fx1(i) - fx1(i+1) |
1057 |
END DO |
1058 |
|
1059 |
! *************************************** |
1060 |
|
1061 |
END DO |
1062 |
RETURN |
1063 |
END SUBROUTINE xtp |
1064 |
|
1065 |
SUBROUTINE fxppm(imr, iml, ut, p, dc, flux, iord) |
1066 |
PARAMETER (r3=1./3., r23=2./3.) |
1067 |
DIMENSION ut(*), flux(*), p(-iml:imr+iml+1), dc(-iml:imr+iml+1) |
1068 |
DIMENSION ar(0:imr), al(0:imr), a6(0:imr) |
1069 |
INTEGER lmt |
1070 |
! logical first |
1071 |
! data first /.true./ |
1072 |
! SAVE LMT |
1073 |
! if(first) then |
1074 |
|
1075 |
! correction calcul de LMT a chaque passage pour pouvoir choisir |
1076 |
! plusieurs schemas PPM pour differents traceurs |
1077 |
! IF (IORD.LE.0) then |
1078 |
! if(IMR.GE.144) then |
1079 |
! LMT = 0 |
1080 |
! elseif(IMR.GE.72) then |
1081 |
! LMT = 1 |
1082 |
! else |
1083 |
! LMT = 2 |
1084 |
! endif |
1085 |
! else |
1086 |
! LMT = IORD - 3 |
1087 |
! endif |
1088 |
|
1089 |
lmt = iord - 3 |
1090 |
|
1091 |
DO i = 1, imr |
1092 |
al(i) = 0.5*(p(i-1)+p(i)) + (dc(i-1)-dc(i))*r3 |
1093 |
END DO |
1094 |
|
1095 |
DO i = 1, imr - 1 |
1096 |
ar(i) = al(i+1) |
1097 |
END DO |
1098 |
ar(imr) = al(1) |
1099 |
|
1100 |
DO i = 1, imr |
1101 |
a6(i) = 3.*(p(i)+p(i)-(al(i)+ar(i))) |
1102 |
END DO |
1103 |
|
1104 |
IF (lmt<=2) CALL lmtppm(dc(1), a6(1), ar(1), al(1), p(1), imr, lmt) |
1105 |
|
1106 |
al(0) = al(imr) |
1107 |
ar(0) = ar(imr) |
1108 |
a6(0) = a6(imr) |
1109 |
|
1110 |
DO i = 1, imr |
1111 |
IF (ut(i)>0.) THEN |
1112 |
flux(i) = ar(i-1) + 0.5*ut(i)*(al(i-1)-ar(i-1)+a6(i-1)*(1.-r23*ut(i))) |
1113 |
ELSE |
1114 |
flux(i) = al(i) - 0.5*ut(i)*(ar(i)-al(i)+a6(i)*(1.+r23*ut(i))) |
1115 |
END IF |
1116 |
END DO |
1117 |
RETURN |
1118 |
END SUBROUTINE fxppm |
1119 |
|
1120 |
SUBROUTINE xmist(imr, iml, p, dc) |
1121 |
PARAMETER (r24=1./24.) |
1122 |
DIMENSION p(-iml:imr+1+iml), dc(-iml:imr+1+iml) |
1123 |
|
1124 |
DO i = 1, imr |
1125 |
tmp = r24*(8.*(p(i+1)-p(i-1))+p(i-2)-p(i+2)) |
1126 |
pmax = max(p(i-1), p(i), p(i+1)) - p(i) |
1127 |
pmin = p(i) - min(p(i-1), p(i), p(i+1)) |
1128 |
dc(i) = sign(min(abs(tmp),pmax,pmin), tmp) |
1129 |
END DO |
1130 |
RETURN |
1131 |
END SUBROUTINE xmist |
1132 |
|
1133 |
SUBROUTINE ytp(imr, jnp, j1, j2, acosp, rcap, dq, p, vc, dc2, ymass, fx, a6, & |
1134 |
ar, al, jord) |
1135 |
DIMENSION p(imr, jnp), vc(imr, jnp), ymass(imr, jnp), dc2(imr, jnp), & |
1136 |
dq(imr, jnp), acosp(jnp) |
1137 |
! Work array |
1138 |
DIMENSION fx(imr, jnp), ar(imr, jnp), al(imr, jnp), a6(imr, jnp) |
1139 |
|
1140 |
jmr = jnp - 1 |
1141 |
len = imr*(j2-j1+2) |
1142 |
|
1143 |
IF (jord==1) THEN |
1144 |
DO i = 1, len |
1145 |
jt = float(j1) - vc(i, j1) |
1146 |
fx(i, j1) = p(i, jt) |
1147 |
END DO |
1148 |
ELSE |
1149 |
|
1150 |
CALL ymist(imr, jnp, j1, p, dc2, 4) |
1151 |
|
1152 |
IF (jord<=0 .OR. jord>=3) THEN |
1153 |
|
1154 |
CALL fyppm(vc, p, dc2, fx, imr, jnp, j1, j2, a6, ar, al, jord) |
1155 |
|
1156 |
ELSE |
1157 |
DO i = 1, len |
1158 |
jt = float(j1) - vc(i, j1) |
1159 |
fx(i, j1) = p(i, jt) + (sign(1.,vc(i,j1))-vc(i,j1))*dc2(i, jt) |
1160 |
END DO |
1161 |
END IF |
1162 |
END IF |
1163 |
|
1164 |
DO i = 1, len |
1165 |
fx(i, j1) = fx(i, j1)*ymass(i, j1) |
1166 |
END DO |
1167 |
|
1168 |
DO j = j1, j2 |
1169 |
DO i = 1, imr |
1170 |
dq(i, j) = dq(i, j) + (fx(i,j)-fx(i,j+1))*acosp(j) |
1171 |
END DO |
1172 |
END DO |
1173 |
|
1174 |
! Poles |
1175 |
sum1 = fx(imr, j1) |
1176 |
sum2 = fx(imr, j2+1) |
1177 |
DO i = 1, imr - 1 |
1178 |
sum1 = sum1 + fx(i, j1) |
1179 |
sum2 = sum2 + fx(i, j2+1) |
1180 |
END DO |
1181 |
|
1182 |
sum1 = dq(1, 1) - sum1*rcap |
1183 |
sum2 = dq(1, jnp) + sum2*rcap |
1184 |
DO i = 1, imr |
1185 |
dq(i, 1) = sum1 |
1186 |
dq(i, jnp) = sum2 |
1187 |
END DO |
1188 |
|
1189 |
IF (j1/=2) THEN |
1190 |
DO i = 1, imr |
1191 |
dq(i, 2) = sum1 |
1192 |
dq(i, jmr) = sum2 |
1193 |
END DO |
1194 |
END IF |
1195 |
|
1196 |
RETURN |
1197 |
END SUBROUTINE ytp |
1198 |
|
1199 |
SUBROUTINE ymist(imr, jnp, j1, p, dc, id) |
1200 |
PARAMETER (r24=1./24.) |
1201 |
DIMENSION p(imr, jnp), dc(imr, jnp) |
1202 |
|
1203 |
imh = imr/2 |
1204 |
jmr = jnp - 1 |
1205 |
ijm3 = imr*(jmr-3) |
1206 |
|
1207 |
IF (id==2) THEN |
1208 |
DO i = 1, imr*(jmr-1) |
1209 |
tmp = 0.25*(p(i,3)-p(i,1)) |
1210 |
pmax = max(p(i,1), p(i,2), p(i,3)) - p(i, 2) |
1211 |
pmin = p(i, 2) - min(p(i,1), p(i,2), p(i,3)) |
1212 |
dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp) |
1213 |
END DO |
1214 |
ELSE |
1215 |
DO i = 1, imh |
1216 |
! J=2 |
1217 |
tmp = (8.*(p(i,3)-p(i,1))+p(i+imh,2)-p(i,4))*r24 |
1218 |
pmax = max(p(i,1), p(i,2), p(i,3)) - p(i, 2) |
1219 |
pmin = p(i, 2) - min(p(i,1), p(i,2), p(i,3)) |
1220 |
dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp) |
1221 |
! J=JMR |
1222 |
tmp = (8.*(p(i,jnp)-p(i,jmr-1))+p(i,jmr-2)-p(i+imh,jmr))*r24 |
1223 |
pmax = max(p(i,jmr-1), p(i,jmr), p(i,jnp)) - p(i, jmr) |
1224 |
pmin = p(i, jmr) - min(p(i,jmr-1), p(i,jmr), p(i,jnp)) |
1225 |
dc(i, jmr) = sign(min(abs(tmp),pmin,pmax), tmp) |
1226 |
END DO |
1227 |
DO i = imh + 1, imr |
1228 |
! J=2 |
1229 |
tmp = (8.*(p(i,3)-p(i,1))+p(i-imh,2)-p(i,4))*r24 |
1230 |
pmax = max(p(i,1), p(i,2), p(i,3)) - p(i, 2) |
1231 |
pmin = p(i, 2) - min(p(i,1), p(i,2), p(i,3)) |
1232 |
dc(i, 2) = sign(min(abs(tmp),pmin,pmax), tmp) |
1233 |
! J=JMR |
1234 |
tmp = (8.*(p(i,jnp)-p(i,jmr-1))+p(i,jmr-2)-p(i-imh,jmr))*r24 |
1235 |
pmax = max(p(i,jmr-1), p(i,jmr), p(i,jnp)) - p(i, jmr) |
1236 |
pmin = p(i, jmr) - min(p(i,jmr-1), p(i,jmr), p(i,jnp)) |
1237 |
dc(i, jmr) = sign(min(abs(tmp),pmin,pmax), tmp) |
1238 |
END DO |
1239 |
|
1240 |
DO i = 1, ijm3 |
1241 |
tmp = (8.*(p(i,4)-p(i,2))+p(i,1)-p(i,5))*r24 |
1242 |
pmax = max(p(i,2), p(i,3), p(i,4)) - p(i, 3) |
1243 |
pmin = p(i, 3) - min(p(i,2), p(i,3), p(i,4)) |
1244 |
dc(i, 3) = sign(min(abs(tmp),pmin,pmax), tmp) |
1245 |
END DO |
1246 |
END IF |
1247 |
|
1248 |
IF (j1/=2) THEN |
1249 |
DO i = 1, imr |
1250 |
dc(i, 1) = 0. |
1251 |
dc(i, jnp) = 0. |
1252 |
END DO |
1253 |
ELSE |
1254 |
! Determine slopes in polar caps for scalars! |
1255 |
|
1256 |
DO i = 1, imh |
1257 |
! South |
1258 |
tmp = 0.25*(p(i,2)-p(i+imh,2)) |
1259 |
pmax = max(p(i,2), p(i,1), p(i+imh,2)) - p(i, 1) |
1260 |
pmin = p(i, 1) - min(p(i,2), p(i,1), p(i+imh,2)) |
1261 |
dc(i, 1) = sign(min(abs(tmp),pmax,pmin), tmp) |
1262 |
! North. |
1263 |
tmp = 0.25*(p(i+imh,jmr)-p(i,jmr)) |
1264 |
pmax = max(p(i+imh,jmr), p(i,jnp), p(i,jmr)) - p(i, jnp) |
1265 |
pmin = p(i, jnp) - min(p(i+imh,jmr), p(i,jnp), p(i,jmr)) |
1266 |
dc(i, jnp) = sign(min(abs(tmp),pmax,pmin), tmp) |
1267 |
END DO |
1268 |
|
1269 |
DO i = imh + 1, imr |
1270 |
dc(i, 1) = -dc(i-imh, 1) |
1271 |
dc(i, jnp) = -dc(i-imh, jnp) |
1272 |
END DO |
1273 |
END IF |
1274 |
RETURN |
1275 |
END SUBROUTINE ymist |
1276 |
|
1277 |
SUBROUTINE fyppm(vc, p, dc, flux, imr, jnp, j1, j2, a6, ar, al, jord) |
1278 |
PARAMETER (r3=1./3., r23=2./3.) |
1279 |
REAL vc(imr, *), flux(imr, *), p(imr, *), dc(imr, *) |
1280 |
! Local work arrays. |
1281 |
REAL ar(imr, jnp), al(imr, jnp), a6(imr, jnp) |
1282 |
INTEGER lmt |
1283 |
! logical first |
1284 |
! data first /.true./ |
1285 |
! SAVE LMT |
1286 |
|
1287 |
imh = imr/2 |
1288 |
jmr = jnp - 1 |
1289 |
j11 = j1 - 1 |
1290 |
imjm1 = imr*(j2-j1+2) |
1291 |
len = imr*(j2-j1+3) |
1292 |
! if(first) then |
1293 |
! IF(JORD.LE.0) then |
1294 |
! if(JMR.GE.90) then |
1295 |
! LMT = 0 |
1296 |
! elseif(JMR.GE.45) then |
1297 |
! LMT = 1 |
1298 |
! else |
1299 |
! LMT = 2 |
1300 |
! endif |
1301 |
! else |
1302 |
! LMT = JORD - 3 |
1303 |
! endif |
1304 |
|
1305 |
! first = .false. |
1306 |
! endif |
1307 |
|
1308 |
! modifs pour pouvoir choisir plusieurs schemas PPM |
1309 |
lmt = jord - 3 |
1310 |
|
1311 |
DO i = 1, imr*jmr |
1312 |
al(i, 2) = 0.5*(p(i,1)+p(i,2)) + (dc(i,1)-dc(i,2))*r3 |
1313 |
ar(i, 1) = al(i, 2) |
1314 |
END DO |
1315 |
|
1316 |
! Poles: |
1317 |
|
1318 |
DO i = 1, imh |
1319 |
al(i, 1) = al(i+imh, 2) |
1320 |
al(i+imh, 1) = al(i, 2) |
1321 |
|
1322 |
ar(i, jnp) = ar(i+imh, jmr) |
1323 |
ar(i+imh, jnp) = ar(i, jmr) |
1324 |
END DO |
1325 |
|
1326 |
! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
1327 |
! Rajout pour LMDZ.3.3 |
1328 |
! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
1329 |
ar(imr, 1) = al(1, 1) |
1330 |
ar(imr, jnp) = al(1, jnp) |
1331 |
! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
1332 |
|
1333 |
|
1334 |
DO i = 1, len |
1335 |
a6(i, j11) = 3.*(p(i,j11)+p(i,j11)-(al(i,j11)+ar(i,j11))) |
1336 |
END DO |
1337 |
|
1338 |
IF (lmt<=2) CALL lmtppm(dc(1,j11), a6(1,j11), ar(1,j11), al(1,j11), & |
1339 |
p(1,j11), len, lmt) |
1340 |
|
1341 |
|
1342 |
DO i = 1, imjm1 |
1343 |
IF (vc(i,j1)>0.) THEN |
1344 |
flux(i, j1) = ar(i, j11) + 0.5*vc(i, j1)*(al(i,j11)-ar(i,j11)+a6(i,j11) & |
1345 |
*(1.-r23*vc(i,j1))) |
1346 |
ELSE |
1347 |
flux(i, j1) = al(i, j1) - 0.5*vc(i, j1)*(ar(i,j1)-al(i,j1)+a6(i,j1)*(1. & |
1348 |
+r23*vc(i,j1))) |
1349 |
END IF |
1350 |
END DO |
1351 |
RETURN |
1352 |
END SUBROUTINE fyppm |
1353 |
|
1354 |
SUBROUTINE yadv(imr, jnp, j1, j2, p, va, ady, wk, iad) |
1355 |
REAL p(imr, jnp), ady(imr, jnp), va(imr, jnp) |
1356 |
REAL wk(imr, -1:jnp+2) |
1357 |
|
1358 |
jmr = jnp - 1 |
1359 |
imh = imr/2 |
1360 |
DO j = 1, jnp |
1361 |
DO i = 1, imr |
1362 |
wk(i, j) = p(i, j) |
1363 |
END DO |
1364 |
END DO |
1365 |
! Poles: |
1366 |
DO i = 1, imh |
1367 |
wk(i, -1) = p(i+imh, 3) |
1368 |
wk(i+imh, -1) = p(i, 3) |
1369 |
wk(i, 0) = p(i+imh, 2) |
1370 |
wk(i+imh, 0) = p(i, 2) |
1371 |
wk(i, jnp+1) = p(i+imh, jmr) |
1372 |
wk(i+imh, jnp+1) = p(i, jmr) |
1373 |
wk(i, jnp+2) = p(i+imh, jnp-2) |
1374 |
wk(i+imh, jnp+2) = p(i, jnp-2) |
1375 |
END DO |
1376 |
|
1377 |
IF (iad==2) THEN |
1378 |
DO j = j1 - 1, j2 + 1 |
1379 |
DO i = 1, imr |
1380 |
jp = nint(va(i,j)) |
1381 |
rv = jp - va(i, j) |
1382 |
jp = j - jp |
1383 |
a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i, jp) |
1384 |
b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1)) |
1385 |
ady(i, j) = wk(i, jp) + rv*(a1*rv+b1) - wk(i, j) |
1386 |
END DO |
1387 |
END DO |
1388 |
|
1389 |
ELSE IF (iad==1) THEN |
1390 |
DO j = j1 - 1, j2 + 1 |
1391 |
DO i = 1, imr |
1392 |
jp = float(j) - va(i, j) |
1393 |
ady(i, j) = va(i, j)*(wk(i,jp)-wk(i,jp+1)) |
1394 |
END DO |
1395 |
END DO |
1396 |
END IF |
1397 |
|
1398 |
IF (j1/=2) THEN |
1399 |
sum1 = 0. |
1400 |
sum2 = 0. |
1401 |
DO i = 1, imr |
1402 |
sum1 = sum1 + ady(i, 2) |
1403 |
sum2 = sum2 + ady(i, jmr) |
1404 |
END DO |
1405 |
sum1 = sum1/imr |
1406 |
sum2 = sum2/imr |
1407 |
|
1408 |
DO i = 1, imr |
1409 |
ady(i, 2) = sum1 |
1410 |
ady(i, jmr) = sum2 |
1411 |
ady(i, 1) = sum1 |
1412 |
ady(i, jnp) = sum2 |
1413 |
END DO |
1414 |
ELSE |
1415 |
! Poles: |
1416 |
sum1 = 0. |
1417 |
sum2 = 0. |
1418 |
DO i = 1, imr |
1419 |
sum1 = sum1 + ady(i, 1) |
1420 |
sum2 = sum2 + ady(i, jnp) |
1421 |
END DO |
1422 |
sum1 = sum1/imr |
1423 |
sum2 = sum2/imr |
1424 |
|
1425 |
DO i = 1, imr |
1426 |
ady(i, 1) = sum1 |
1427 |
ady(i, jnp) = sum2 |
1428 |
END DO |
1429 |
END IF |
1430 |
|
1431 |
RETURN |
1432 |
END SUBROUTINE yadv |
1433 |
|
1434 |
SUBROUTINE xadv(imr, jnp, j1, j2, p, ua, js, jn, iml, adx, iad) |
1435 |
REAL p(imr, jnp), adx(imr, jnp), qtmp(-imr:imr+imr), ua(imr, jnp) |
1436 |
|
1437 |
jmr = jnp - 1 |
1438 |
DO j = j1, j2 |
1439 |
IF (j>js .AND. j<jn) GO TO 1309 |
1440 |
|
1441 |
DO i = 1, imr |
1442 |
qtmp(i) = p(i, j) |
1443 |
END DO |
1444 |
|
1445 |
DO i = -iml, 0 |
1446 |
qtmp(i) = p(imr+i, j) |
1447 |
qtmp(imr+1-i) = p(1-i, j) |
1448 |
END DO |
1449 |
|
1450 |
IF (iad==2) THEN |
1451 |
DO i = 1, imr |
1452 |
ip = nint(ua(i,j)) |
1453 |
ru = ip - ua(i, j) |
1454 |
ip = i - ip |
1455 |
a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip) |
1456 |
b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1)) |
1457 |
adx(i, j) = qtmp(ip) + ru*(a1*ru+b1) |
1458 |
END DO |
1459 |
ELSE IF (iad==1) THEN |
1460 |
DO i = 1, imr |
1461 |
iu = ua(i, j) |
1462 |
ru = ua(i, j) - iu |
1463 |
iiu = i - iu |
1464 |
IF (ua(i,j)>=0.) THEN |
1465 |
adx(i, j) = qtmp(iiu) + ru*(qtmp(iiu-1)-qtmp(iiu)) |
1466 |
ELSE |
1467 |
adx(i, j) = qtmp(iiu) + ru*(qtmp(iiu)-qtmp(iiu+1)) |
1468 |
END IF |
1469 |
END DO |
1470 |
END IF |
1471 |
|
1472 |
DO i = 1, imr |
1473 |
adx(i, j) = adx(i, j) - p(i, j) |
1474 |
END DO |
1475 |
1309 END DO |
1476 |
|
1477 |
! Eulerian upwind |
1478 |
|
1479 |
DO j = js + 1, jn - 1 |
1480 |
|
1481 |
DO i = 1, imr |
1482 |
qtmp(i) = p(i, j) |
1483 |
END DO |
1484 |
|
1485 |
qtmp(0) = p(imr, j) |
1486 |
qtmp(imr+1) = p(1, j) |
1487 |
|
1488 |
IF (iad==2) THEN |
1489 |
qtmp(-1) = p(imr-1, j) |
1490 |
qtmp(imr+2) = p(2, j) |
1491 |
DO i = 1, imr |
1492 |
ip = nint(ua(i,j)) |
1493 |
ru = ip - ua(i, j) |
1494 |
ip = i - ip |
1495 |
a1 = 0.5*(qtmp(ip+1)+qtmp(ip-1)) - qtmp(ip) |
1496 |
b1 = 0.5*(qtmp(ip+1)-qtmp(ip-1)) |
1497 |
adx(i, j) = qtmp(ip) - p(i, j) + ru*(a1*ru+b1) |
1498 |
END DO |
1499 |
ELSE IF (iad==1) THEN |
1500 |
! 1st order |
1501 |
DO i = 1, imr |
1502 |
ip = i - ua(i, j) |
1503 |
adx(i, j) = ua(i, j)*(qtmp(ip)-qtmp(ip+1)) |
1504 |
END DO |
1505 |
END IF |
1506 |
END DO |
1507 |
|
1508 |
IF (j1/=2) THEN |
1509 |
DO i = 1, imr |
1510 |
adx(i, 2) = 0. |
1511 |
adx(i, jmr) = 0. |
1512 |
END DO |
1513 |
END IF |
1514 |
! set cross term due to x-adv at the poles to zero. |
1515 |
DO i = 1, imr |
1516 |
adx(i, 1) = 0. |
1517 |
adx(i, jnp) = 0. |
1518 |
END DO |
1519 |
RETURN |
1520 |
END SUBROUTINE xadv |
1521 |
|
1522 |
SUBROUTINE lmtppm(dc, a6, ar, al, p, im, lmt) |
1523 |
|
1524 |
! A6 = CURVATURE OF THE TEST PARABOLA |
1525 |
! AR = RIGHT EDGE VALUE OF THE TEST PARABOLA |
1526 |
! AL = LEFT EDGE VALUE OF THE TEST PARABOLA |
1527 |
! DC = 0.5 * MISMATCH |
1528 |
! P = CELL-AVERAGED VALUE |
1529 |
! IM = VECTOR LENGTH |
1530 |
|
1531 |
! OPTIONS: |
1532 |
|
1533 |
! LMT = 0: FULL MONOTONICITY |
1534 |
! LMT = 1: SEMI-MONOTONIC CONSTRAINT (NO UNDERSHOOTS) |
1535 |
! LMT = 2: POSITIVE-DEFINITE CONSTRAINT |
1536 |
|
1537 |
PARAMETER (r12=1./12.) |
1538 |
DIMENSION a6(im), ar(im), al(im), p(im), dc(im) |
1539 |
|
1540 |
IF (lmt==0) THEN |
1541 |
! Full constraint |
1542 |
DO i = 1, im |
1543 |
IF (dc(i)==0.) THEN |
1544 |
ar(i) = p(i) |
1545 |
al(i) = p(i) |
1546 |
a6(i) = 0. |
1547 |
ELSE |
1548 |
da1 = ar(i) - al(i) |
1549 |
da2 = da1**2 |
1550 |
a6da = a6(i)*da1 |
1551 |
IF (a6da<-da2) THEN |
1552 |
a6(i) = 3.*(al(i)-p(i)) |
1553 |
ar(i) = al(i) - a6(i) |
1554 |
ELSE IF (a6da>da2) THEN |
1555 |
a6(i) = 3.*(ar(i)-p(i)) |
1556 |
al(i) = ar(i) - a6(i) |
1557 |
END IF |
1558 |
END IF |
1559 |
END DO |
1560 |
ELSE IF (lmt==1) THEN |
1561 |
! Semi-monotonic constraint |
1562 |
DO i = 1, im |
1563 |
IF (abs(ar(i)-al(i))>=-a6(i)) GO TO 150 |
1564 |
IF (p(i)<ar(i) .AND. p(i)<al(i)) THEN |
1565 |
ar(i) = p(i) |
1566 |
al(i) = p(i) |
1567 |
a6(i) = 0. |
1568 |
ELSE IF (ar(i)>al(i)) THEN |
1569 |
a6(i) = 3.*(al(i)-p(i)) |
1570 |
ar(i) = al(i) - a6(i) |
1571 |
ELSE |
1572 |
a6(i) = 3.*(ar(i)-p(i)) |
1573 |
al(i) = ar(i) - a6(i) |
1574 |
END IF |
1575 |
150 END DO |
1576 |
ELSE IF (lmt==2) THEN |
1577 |
DO i = 1, im |
1578 |
IF (abs(ar(i)-al(i))>=-a6(i)) GO TO 250 |
1579 |
fmin = p(i) + 0.25*(ar(i)-al(i))**2/a6(i) + a6(i)*r12 |
1580 |
IF (fmin>=0.) GO TO 250 |
1581 |
IF (p(i)<ar(i) .AND. p(i)<al(i)) THEN |
1582 |
ar(i) = p(i) |
1583 |
al(i) = p(i) |
1584 |
a6(i) = 0. |
1585 |
ELSE IF (ar(i)>al(i)) THEN |
1586 |
a6(i) = 3.*(al(i)-p(i)) |
1587 |
ar(i) = al(i) - a6(i) |
1588 |
ELSE |
1589 |
a6(i) = 3.*(ar(i)-p(i)) |
1590 |
al(i) = ar(i) - a6(i) |
1591 |
END IF |
1592 |
250 END DO |
1593 |
END IF |
1594 |
RETURN |
1595 |
END SUBROUTINE lmtppm |
1596 |
|
1597 |
SUBROUTINE a2c(u, v, imr, jmr, j1, j2, crx, cry, dtdx5, dtdy5) |
1598 |
DIMENSION u(imr, *), v(imr, *), crx(imr, *), cry(imr, *), dtdx5(*) |
1599 |
|
1600 |
DO j = j1, j2 |
1601 |
DO i = 2, imr |
1602 |
crx(i, j) = dtdx5(j)*(u(i,j)+u(i-1,j)) |
1603 |
END DO |
1604 |
END DO |
1605 |
|
1606 |
DO j = j1, j2 |
1607 |
crx(1, j) = dtdx5(j)*(u(1,j)+u(imr,j)) |
1608 |
END DO |
1609 |
|
1610 |
DO i = 1, imr*jmr |
1611 |
cry(i, 2) = dtdy5*(v(i,2)+v(i,1)) |
1612 |
END DO |
1613 |
RETURN |
1614 |
END SUBROUTINE a2c |
1615 |
|
1616 |
SUBROUTINE cosa(cosp, cose, jnp, pi, dp) |
1617 |
DIMENSION cosp(*), cose(*) |
1618 |
|
1619 |
jmr = jnp - 1 |
1620 |
DO j = 2, jnp |
1621 |
ph5 = -0.5*pi + (float(j-1)-0.5)*dp |
1622 |
cose(j) = cos(ph5) |
1623 |
END DO |
1624 |
|
1625 |
jeq = (jnp+1)/2 |
1626 |
IF (jmr==2*(jmr/2)) THEN |
1627 |
DO j = jnp, jeq + 1, -1 |
1628 |
cose(j) = cose(jnp+2-j) |
1629 |
END DO |
1630 |
ELSE |
1631 |
! cell edge at equator. |
1632 |
cose(jeq+1) = 1. |
1633 |
DO j = jnp, jeq + 2, -1 |
1634 |
cose(j) = cose(jnp+2-j) |
1635 |
END DO |
1636 |
END IF |
1637 |
|
1638 |
DO j = 2, jmr |
1639 |
cosp(j) = 0.5*(cose(j)+cose(j+1)) |
1640 |
END DO |
1641 |
cosp(1) = 0. |
1642 |
cosp(jnp) = 0. |
1643 |
RETURN |
1644 |
END SUBROUTINE cosa |
1645 |
|
1646 |
SUBROUTINE cosc(cosp, cose, jnp, pi, dp) |
1647 |
DIMENSION cosp(*), cose(*) |
1648 |
|
1649 |
phi = -0.5*pi |
1650 |
DO j = 2, jnp - 1 |
1651 |
phi = phi + dp |
1652 |
cosp(j) = cos(phi) |
1653 |
END DO |
1654 |
cosp(1) = 0. |
1655 |
cosp(jnp) = 0. |
1656 |
|
1657 |
DO j = 2, jnp |
1658 |
cose(j) = 0.5*(cosp(j)+cosp(j-1)) |
1659 |
END DO |
1660 |
|
1661 |
DO j = 2, jnp - 1 |
1662 |
cosp(j) = 0.5*(cose(j)+cose(j+1)) |
1663 |
END DO |
1664 |
RETURN |
1665 |
END SUBROUTINE cosc |
1666 |
|
1667 |
SUBROUTINE qckxyz(q, qtmp, imr, jnp, nlay, j1, j2, cosp, acosp, cross, ic, & |
1668 |
nstep) |
1669 |
|
1670 |
PARAMETER (tiny=1.E-60) |
1671 |
DIMENSION q(imr, jnp, nlay), qtmp(imr, jnp), cosp(*), acosp(*) |
1672 |
LOGICAL cross |
1673 |
|
1674 |
nlaym1 = nlay - 1 |
1675 |
len = imr*(j2-j1+1) |
1676 |
ip = 0 |
1677 |
|
1678 |
! Top layer |
1679 |
l = 1 |
1680 |
icr = 1 |
1681 |
CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, tiny) |
1682 |
IF (ipy==0) GO TO 50 |
1683 |
CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny) |
1684 |
IF (ipx==0) GO TO 50 |
1685 |
|
1686 |
IF (cross) THEN |
1687 |
CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny) |
1688 |
END IF |
1689 |
IF (icr==0) GO TO 50 |
1690 |
|
1691 |
! Vertical filling... |
1692 |
DO i = 1, len |
1693 |
IF (q(i,j1,1)<0.) THEN |
1694 |
ip = ip + 1 |
1695 |
q(i, j1, 2) = q(i, j1, 2) + q(i, j1, 1) |
1696 |
q(i, j1, 1) = 0. |
1697 |
END IF |
1698 |
END DO |
1699 |
|
1700 |
50 CONTINUE |
1701 |
DO l = 2, nlaym1 |
1702 |
icr = 1 |
1703 |
|
1704 |
CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, tiny) |
1705 |
IF (ipy==0) GO TO 225 |
1706 |
CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny) |
1707 |
IF (ipx==0) GO TO 225 |
1708 |
IF (cross) THEN |
1709 |
CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny) |
1710 |
END IF |
1711 |
IF (icr==0) GO TO 225 |
1712 |
|
1713 |
DO i = 1, len |
1714 |
IF (q(i,j1,l)<0.) THEN |
1715 |
|
1716 |
ip = ip + 1 |
1717 |
! From above |
1718 |
qup = q(i, j1, l-1) |
1719 |
qly = -q(i, j1, l) |
1720 |
dup = min(qly, qup) |
1721 |
q(i, j1, l-1) = qup - dup |
1722 |
q(i, j1, l) = dup - qly |
1723 |
! Below |
1724 |
q(i, j1, l+1) = q(i, j1, l+1) + q(i, j1, l) |
1725 |
q(i, j1, l) = 0. |
1726 |
END IF |
1727 |
END DO |
1728 |
225 END DO |
1729 |
|
1730 |
! BOTTOM LAYER |
1731 |
sum = 0. |
1732 |
l = nlay |
1733 |
|
1734 |
CALL filns(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, ipy, tiny) |
1735 |
IF (ipy==0) GO TO 911 |
1736 |
CALL filew(q(1,1,l), qtmp, imr, jnp, j1, j2, ipx, tiny) |
1737 |
IF (ipx==0) GO TO 911 |
1738 |
|
1739 |
CALL filcr(q(1,1,l), imr, jnp, j1, j2, cosp, acosp, icr, tiny) |
1740 |
IF (icr==0) GO TO 911 |
1741 |
|
1742 |
DO i = 1, len |
1743 |
IF (q(i,j1,l)<0.) THEN |
1744 |
ip = ip + 1 |
1745 |
|
1746 |
! From above |
1747 |
|
1748 |
qup = q(i, j1, nlaym1) |
1749 |
qly = -q(i, j1, l) |
1750 |
dup = min(qly, qup) |
1751 |
q(i, j1, nlaym1) = qup - dup |
1752 |
! From "below" the surface. |
1753 |
sum = sum + qly - dup |
1754 |
q(i, j1, l) = 0. |
1755 |
END IF |
1756 |
END DO |
1757 |
|
1758 |
911 CONTINUE |
1759 |
|
1760 |
IF (ip>imr) THEN |
1761 |
WRITE (6, *) 'IC=', ic, ' STEP=', nstep, ' Vertical filling pts=', ip |
1762 |
END IF |
1763 |
|
1764 |
IF (sum>1.E-25) THEN |
1765 |
WRITE (6, *) ic, nstep, ' Mass source from the ground=', sum |
1766 |
END IF |
1767 |
RETURN |
1768 |
END SUBROUTINE qckxyz |
1769 |
|
1770 |
SUBROUTINE filcr(q, imr, jnp, j1, j2, cosp, acosp, icr, tiny) |
1771 |
DIMENSION q(imr, *), cosp(*), acosp(*) |
1772 |
|
1773 |
icr = 0 |
1774 |
DO j = j1 + 1, j2 - 1 |
1775 |
DO i = 1, imr - 1 |
1776 |
IF (q(i,j)<0.) THEN |
1777 |
icr = 1 |
1778 |
dq = -q(i, j)*cosp(j) |
1779 |
! N-E |
1780 |
dn = q(i+1, j+1)*cosp(j+1) |
1781 |
d0 = max(0., dn) |
1782 |
d1 = min(dq, d0) |
1783 |
q(i+1, j+1) = (dn-d1)*acosp(j+1) |
1784 |
dq = dq - d1 |
1785 |
! S-E |
1786 |
ds = q(i+1, j-1)*cosp(j-1) |
1787 |
d0 = max(0., ds) |
1788 |
d2 = min(dq, d0) |
1789 |
q(i+1, j-1) = (ds-d2)*acosp(j-1) |
1790 |
q(i, j) = (d2-dq)*acosp(j) + tiny |
1791 |
END IF |
1792 |
END DO |
1793 |
IF (icr==0 .AND. q(imr,j)>=0.) GO TO 65 |
1794 |
DO i = 2, imr |
1795 |
IF (q(i,j)<0.) THEN |
1796 |
icr = 1 |
1797 |
dq = -q(i, j)*cosp(j) |
1798 |
! N-W |
1799 |
dn = q(i-1, j+1)*cosp(j+1) |
1800 |
d0 = max(0., dn) |
1801 |
d1 = min(dq, d0) |
1802 |
q(i-1, j+1) = (dn-d1)*acosp(j+1) |
1803 |
dq = dq - d1 |
1804 |
! S-W |
1805 |
ds = q(i-1, j-1)*cosp(j-1) |
1806 |
d0 = max(0., ds) |
1807 |
d2 = min(dq, d0) |
1808 |
q(i-1, j-1) = (ds-d2)*acosp(j-1) |
1809 |
q(i, j) = (d2-dq)*acosp(j) + tiny |
1810 |
END IF |
1811 |
END DO |
1812 |
! ***************************************** |
1813 |
! i=1 |
1814 |
i = 1 |
1815 |
IF (q(i,j)<0.) THEN |
1816 |
icr = 1 |
1817 |
dq = -q(i, j)*cosp(j) |
1818 |
! N-W |
1819 |
dn = q(imr, j+1)*cosp(j+1) |
1820 |
d0 = max(0., dn) |
1821 |
d1 = min(dq, d0) |
1822 |
q(imr, j+1) = (dn-d1)*acosp(j+1) |
1823 |
dq = dq - d1 |
1824 |
! S-W |
1825 |
ds = q(imr, j-1)*cosp(j-1) |
1826 |
d0 = max(0., ds) |
1827 |
d2 = min(dq, d0) |
1828 |
q(imr, j-1) = (ds-d2)*acosp(j-1) |
1829 |
q(i, j) = (d2-dq)*acosp(j) + tiny |
1830 |
END IF |
1831 |
! ***************************************** |
1832 |
! i=IMR |
1833 |
i = imr |
1834 |
IF (q(i,j)<0.) THEN |
1835 |
icr = 1 |
1836 |
dq = -q(i, j)*cosp(j) |
1837 |
! N-E |
1838 |
dn = q(1, j+1)*cosp(j+1) |
1839 |
d0 = max(0., dn) |
1840 |
d1 = min(dq, d0) |
1841 |
q(1, j+1) = (dn-d1)*acosp(j+1) |
1842 |
dq = dq - d1 |
1843 |
! S-E |
1844 |
ds = q(1, j-1)*cosp(j-1) |
1845 |
d0 = max(0., ds) |
1846 |
d2 = min(dq, d0) |
1847 |
q(1, j-1) = (ds-d2)*acosp(j-1) |
1848 |
q(i, j) = (d2-dq)*acosp(j) + tiny |
1849 |
END IF |
1850 |
! ***************************************** |
1851 |
65 END DO |
1852 |
|
1853 |
DO i = 1, imr |
1854 |
IF (q(i,j1)<0. .OR. q(i,j2)<0.) THEN |
1855 |
icr = 1 |
1856 |
GO TO 80 |
1857 |
END IF |
1858 |
END DO |
1859 |
|
1860 |
80 CONTINUE |
1861 |
|
1862 |
IF (q(1,1)<0. .OR. q(1,jnp)<0.) THEN |
1863 |
icr = 1 |
1864 |
END IF |
1865 |
|
1866 |
RETURN |
1867 |
END SUBROUTINE filcr |
1868 |
|
1869 |
SUBROUTINE filns(q, imr, jnp, j1, j2, cosp, acosp, ipy, tiny) |
1870 |
DIMENSION q(imr, *), cosp(*), acosp(*) |
1871 |
! logical first |
1872 |
! data first /.true./ |
1873 |
! save cap1 |
1874 |
|
1875 |
! if(first) then |
1876 |
dp = 4.*atan(1.)/float(jnp-1) |
1877 |
cap1 = imr*(1.-cos((j1-1.5)*dp))/dp |
1878 |
! first = .false. |
1879 |
! endif |
1880 |
|
1881 |
ipy = 0 |
1882 |
DO j = j1 + 1, j2 - 1 |
1883 |
DO i = 1, imr |
1884 |
IF (q(i,j)<0.) THEN |
1885 |
ipy = 1 |
1886 |
dq = -q(i, j)*cosp(j) |
1887 |
! North |
1888 |
dn = q(i, j+1)*cosp(j+1) |
1889 |
d0 = max(0., dn) |
1890 |
d1 = min(dq, d0) |
1891 |
q(i, j+1) = (dn-d1)*acosp(j+1) |
1892 |
dq = dq - d1 |
1893 |
! South |
1894 |
ds = q(i, j-1)*cosp(j-1) |
1895 |
d0 = max(0., ds) |
1896 |
d2 = min(dq, d0) |
1897 |
q(i, j-1) = (ds-d2)*acosp(j-1) |
1898 |
q(i, j) = (d2-dq)*acosp(j) + tiny |
1899 |
END IF |
1900 |
END DO |
1901 |
END DO |
1902 |
|
1903 |
DO i = 1, imr |
1904 |
IF (q(i,j1)<0.) THEN |
1905 |
ipy = 1 |
1906 |
dq = -q(i, j1)*cosp(j1) |
1907 |
! North |
1908 |
dn = q(i, j1+1)*cosp(j1+1) |
1909 |
d0 = max(0., dn) |
1910 |
d1 = min(dq, d0) |
1911 |
q(i, j1+1) = (dn-d1)*acosp(j1+1) |
1912 |
q(i, j1) = (d1-dq)*acosp(j1) + tiny |
1913 |
END IF |
1914 |
END DO |
1915 |
|
1916 |
j = j2 |
1917 |
DO i = 1, imr |
1918 |
IF (q(i,j)<0.) THEN |
1919 |
ipy = 1 |
1920 |
dq = -q(i, j)*cosp(j) |
1921 |
! South |
1922 |
ds = q(i, j-1)*cosp(j-1) |
1923 |
d0 = max(0., ds) |
1924 |
d2 = min(dq, d0) |
1925 |
q(i, j-1) = (ds-d2)*acosp(j-1) |
1926 |
q(i, j) = (d2-dq)*acosp(j) + tiny |
1927 |
END IF |
1928 |
END DO |
1929 |
|
1930 |
! Check Poles. |
1931 |
IF (q(1,1)<0.) THEN |
1932 |
dq = q(1, 1)*cap1/float(imr)*acosp(j1) |
1933 |
DO i = 1, imr |
1934 |
q(i, 1) = 0. |
1935 |
q(i, j1) = q(i, j1) + dq |
1936 |
IF (q(i,j1)<0.) ipy = 1 |
1937 |
END DO |
1938 |
END IF |
1939 |
|
1940 |
IF (q(1,jnp)<0.) THEN |
1941 |
dq = q(1, jnp)*cap1/float(imr)*acosp(j2) |
1942 |
DO i = 1, imr |
1943 |
q(i, jnp) = 0. |
1944 |
q(i, j2) = q(i, j2) + dq |
1945 |
IF (q(i,j2)<0.) ipy = 1 |
1946 |
END DO |
1947 |
END IF |
1948 |
|
1949 |
RETURN |
1950 |
END SUBROUTINE filns |
1951 |
|
1952 |
SUBROUTINE filew(q, qtmp, imr, jnp, j1, j2, ipx, tiny) |
1953 |
DIMENSION q(imr, *), qtmp(jnp, imr) |
1954 |
|
1955 |
ipx = 0 |
1956 |
! Copy & swap direction for vectorization. |
1957 |
DO i = 1, imr |
1958 |
DO j = j1, j2 |
1959 |
qtmp(j, i) = q(i, j) |
1960 |
END DO |
1961 |
END DO |
1962 |
|
1963 |
DO i = 2, imr - 1 |
1964 |
DO j = j1, j2 |
1965 |
IF (qtmp(j,i)<0.) THEN |
1966 |
ipx = 1 |
1967 |
! west |
1968 |
d0 = max(0., qtmp(j,i-1)) |
1969 |
d1 = min(-qtmp(j,i), d0) |
1970 |
qtmp(j, i-1) = qtmp(j, i-1) - d1 |
1971 |
qtmp(j, i) = qtmp(j, i) + d1 |
1972 |
! east |
1973 |
d0 = max(0., qtmp(j,i+1)) |
1974 |
d2 = min(-qtmp(j,i), d0) |
1975 |
qtmp(j, i+1) = qtmp(j, i+1) - d2 |
1976 |
qtmp(j, i) = qtmp(j, i) + d2 + tiny |
1977 |
END IF |
1978 |
END DO |
1979 |
END DO |
1980 |
|
1981 |
i = 1 |
1982 |
DO j = j1, j2 |
1983 |
IF (qtmp(j,i)<0.) THEN |
1984 |
ipx = 1 |
1985 |
! west |
1986 |
d0 = max(0., qtmp(j,imr)) |
1987 |
d1 = min(-qtmp(j,i), d0) |
1988 |
qtmp(j, imr) = qtmp(j, imr) - d1 |
1989 |
qtmp(j, i) = qtmp(j, i) + d1 |
1990 |
! east |
1991 |
d0 = max(0., qtmp(j,i+1)) |
1992 |
d2 = min(-qtmp(j,i), d0) |
1993 |
qtmp(j, i+1) = qtmp(j, i+1) - d2 |
1994 |
|
1995 |
qtmp(j, i) = qtmp(j, i) + d2 + tiny |
1996 |
END IF |
1997 |
END DO |
1998 |
i = imr |
1999 |
DO j = j1, j2 |
2000 |
IF (qtmp(j,i)<0.) THEN |
2001 |
ipx = 1 |
2002 |
! west |
2003 |
d0 = max(0., qtmp(j,i-1)) |
2004 |
d1 = min(-qtmp(j,i), d0) |
2005 |
qtmp(j, i-1) = qtmp(j, i-1) - d1 |
2006 |
qtmp(j, i) = qtmp(j, i) + d1 |
2007 |
! east |
2008 |
d0 = max(0., qtmp(j,1)) |
2009 |
d2 = min(-qtmp(j,i), d0) |
2010 |
qtmp(j, 1) = qtmp(j, 1) - d2 |
2011 |
|
2012 |
qtmp(j, i) = qtmp(j, i) + d2 + tiny |
2013 |
END IF |
2014 |
END DO |
2015 |
|
2016 |
IF (ipx/=0) THEN |
2017 |
DO j = j1, j2 |
2018 |
DO i = 1, imr |
2019 |
q(i, j) = qtmp(j, i) |
2020 |
END DO |
2021 |
END DO |
2022 |
ELSE |
2023 |
|
2024 |
! Poles. |
2025 |
IF (q(1,1)<0 .OR. q(1,jnp)<0.) ipx = 1 |
2026 |
END IF |
2027 |
RETURN |
2028 |
END SUBROUTINE filew |