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