1 | ! |
---|
2 | ! $Id: vlsplt.F90 5136 2024-07-28 14:17:54Z abarral $ |
---|
3 | ! |
---|
4 | |
---|
5 | SUBROUTINE vlsplt(q, pente_max, masse, w, pbaru, pbarv, pdt, iq) |
---|
6 | USE infotrac, ONLY: nqtot, tracers |
---|
7 | USE lmdz_ssum_scopy, ONLY: scopy |
---|
8 | ! |
---|
9 | ! Auteurs: P.Le Van, F.Hourdin, F.Forget |
---|
10 | ! |
---|
11 | ! ******************************************************************** |
---|
12 | ! Shema d'advection " pseudo amont " . |
---|
13 | ! ******************************************************************** |
---|
14 | ! q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... |
---|
15 | ! |
---|
16 | ! pente_max facteur de limitation des pentes: 2 en general |
---|
17 | ! 0 pour un schema amont |
---|
18 | ! pbaru,pbarv,w flux de masse en u ,v ,w |
---|
19 | ! pdt pas de temps |
---|
20 | ! |
---|
21 | ! -------------------------------------------------------------------- |
---|
22 | IMPLICIT NONE |
---|
23 | ! |
---|
24 | INCLUDE "dimensions.h" |
---|
25 | INCLUDE "paramet.h" |
---|
26 | |
---|
27 | ! |
---|
28 | ! Arguments: |
---|
29 | ! ---------- |
---|
30 | REAL :: masse(ip1jmp1, llm), pente_max |
---|
31 | REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm) |
---|
32 | REAL :: q(ip1jmp1, llm, nqtot) |
---|
33 | REAL :: w(ip1jmp1, llm), pdt |
---|
34 | INTEGER :: iq ! CRisi |
---|
35 | ! |
---|
36 | ! Local |
---|
37 | ! --------- |
---|
38 | ! |
---|
39 | INTEGER :: ij, l |
---|
40 | ! |
---|
41 | REAL :: zm(ip1jmp1, llm, nqtot) |
---|
42 | REAL :: mu(ip1jmp1, llm) |
---|
43 | REAL :: mv(ip1jm, llm) |
---|
44 | REAL :: mw(ip1jmp1, llm + 1) |
---|
45 | REAL :: zq(ip1jmp1, llm, nqtot) |
---|
46 | REAL :: zzpbar, zzw |
---|
47 | INTEGER :: ifils, iq2 ! CRisi |
---|
48 | |
---|
49 | REAL :: qmin, qmax |
---|
50 | DATA qmin, qmax/0., 1.e33/ |
---|
51 | |
---|
52 | zzpbar = 0.5 * pdt |
---|
53 | zzw = pdt |
---|
54 | DO l = 1, llm |
---|
55 | DO ij = iip2, ip1jm |
---|
56 | mu(ij, l) = pbaru(ij, l) * zzpbar |
---|
57 | ENDDO |
---|
58 | DO ij = 1, ip1jm |
---|
59 | mv(ij, l) = pbarv(ij, l) * zzpbar |
---|
60 | ENDDO |
---|
61 | DO ij = 1, ip1jmp1 |
---|
62 | mw(ij, l) = w(ij, l) * zzw |
---|
63 | ENDDO |
---|
64 | ENDDO |
---|
65 | |
---|
66 | DO ij = 1, ip1jmp1 |
---|
67 | mw(ij, llm + 1) = 0. |
---|
68 | ENDDO |
---|
69 | |
---|
70 | CALL SCOPY(ijp1llm, q(1, 1, iq), 1, zq(1, 1, iq), 1) |
---|
71 | CALL SCOPY(ijp1llm, masse, 1, zm(1, 1, iq), 1) |
---|
72 | |
---|
73 | do ifils = 1, tracers(iq)%nqDescen |
---|
74 | iq2 = tracers(iq)%iqDescen(ifils) |
---|
75 | CALL SCOPY(ijp1llm, q(1, 1, iq2), 1, zq(1, 1, iq2), 1) |
---|
76 | enddo |
---|
77 | |
---|
78 | CALL vlx(zq, pente_max, zm, mu, iq) |
---|
79 | CALL vly(zq, pente_max, zm, mv, iq) |
---|
80 | CALL vlz(zq, pente_max, zm, mw, iq) |
---|
81 | CALL vly(zq, pente_max, zm, mv, iq) |
---|
82 | CALL vlx(zq, pente_max, zm, mu, iq) |
---|
83 | |
---|
84 | DO l = 1, llm |
---|
85 | DO ij = 1, ip1jmp1 |
---|
86 | q(ij, l, iq) = zq(ij, l, iq) |
---|
87 | ENDDO |
---|
88 | DO ij = 1, ip1jm + 1, iip1 |
---|
89 | q(ij + iim, l, iq) = q(ij, l, iq) |
---|
90 | ENDDO |
---|
91 | ENDDO |
---|
92 | ! CRisi: aussi pour les fils |
---|
93 | do ifils = 1, tracers(iq)%nqDescen |
---|
94 | iq2 = tracers(iq)%iqDescen(ifils) |
---|
95 | DO l = 1, llm |
---|
96 | DO ij = 1, ip1jmp1 |
---|
97 | q(ij, l, iq2) = zq(ij, l, iq2) |
---|
98 | ENDDO |
---|
99 | DO ij = 1, ip1jm + 1, iip1 |
---|
100 | q(ij + iim, l, iq2) = q(ij, l, iq2) |
---|
101 | ENDDO |
---|
102 | ENDDO |
---|
103 | enddo |
---|
104 | |
---|
105 | END SUBROUTINE vlsplt |
---|
106 | RECURSIVE SUBROUTINE vlx(q, pente_max, masse, u_m, iq) |
---|
107 | USE infotrac, ONLY: nqtot, tracers, & ! CRisi |
---|
108 | min_qParent, min_qMass, min_ratio ! MVals et CRisi |
---|
109 | USE lmdz_iniprint, ONLY: lunout, prt_level |
---|
110 | |
---|
111 | ! Auteurs: P.Le Van, F.Hourdin, F.Forget |
---|
112 | ! |
---|
113 | ! ******************************************************************** |
---|
114 | ! Shema d'advection " pseudo amont " . |
---|
115 | ! ******************************************************************** |
---|
116 | ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... |
---|
117 | ! |
---|
118 | ! |
---|
119 | ! -------------------------------------------------------------------- |
---|
120 | IMPLICIT NONE |
---|
121 | ! |
---|
122 | INCLUDE "dimensions.h" |
---|
123 | INCLUDE "paramet.h" |
---|
124 | ! |
---|
125 | ! |
---|
126 | ! Arguments: |
---|
127 | ! ---------- |
---|
128 | REAL :: masse(ip1jmp1, llm, nqtot), pente_max |
---|
129 | REAL :: u_m(ip1jmp1, llm) |
---|
130 | REAL :: q(ip1jmp1, llm, nqtot) |
---|
131 | INTEGER :: iq ! CRisi |
---|
132 | ! |
---|
133 | ! Local |
---|
134 | ! --------- |
---|
135 | ! |
---|
136 | INTEGER :: ij, l, j, i, iju, ijq, indu(ip1jmp1), niju |
---|
137 | INTEGER :: n0, iadvplus(ip1jmp1, llm), nl(llm) |
---|
138 | ! |
---|
139 | REAL :: new_m, zu_m, zdum(ip1jmp1, llm) |
---|
140 | REAL :: dxq(ip1jmp1, llm), dxqu(ip1jmp1) |
---|
141 | REAL :: zz(ip1jmp1) |
---|
142 | REAL :: adxqu(ip1jmp1), dxqmax(ip1jmp1, llm) |
---|
143 | REAL :: u_mq(ip1jmp1, llm) |
---|
144 | |
---|
145 | ! CRisi |
---|
146 | REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) |
---|
147 | INTEGER :: ifils, iq2 ! CRisi |
---|
148 | |
---|
149 | LOGICAL, SAVE :: first |
---|
150 | DATA first/.TRUE./ |
---|
151 | |
---|
152 | ! calcul de la pente a droite et a gauche de la maille |
---|
153 | |
---|
154 | IF (pente_max>-1.e-5) THEN |
---|
155 | ! IF (pente_max.gt.10) THEN |
---|
156 | |
---|
157 | ! calcul des pentes avec limitation, Van Leer scheme I: |
---|
158 | ! ----------------------------------------------------- |
---|
159 | |
---|
160 | ! calcul de la pente aux points u |
---|
161 | DO l = 1, llm |
---|
162 | DO ij = iip2, ip1jm - 1 |
---|
163 | dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq) |
---|
164 | ENDDO |
---|
165 | DO ij = iip1 + iip1, ip1jm, iip1 |
---|
166 | dxqu(ij) = dxqu(ij - iim) |
---|
167 | ! sigu(ij)=sigu(ij-iim) |
---|
168 | ENDDO |
---|
169 | |
---|
170 | DO ij = iip2, ip1jm |
---|
171 | adxqu(ij) = abs(dxqu(ij)) |
---|
172 | ENDDO |
---|
173 | |
---|
174 | ! calcul de la pente maximum dans la maille en valeur absolue |
---|
175 | |
---|
176 | DO ij = iip2 + 1, ip1jm |
---|
177 | dxqmax(ij, l) = pente_max * & |
---|
178 | min(adxqu(ij - 1), adxqu(ij)) |
---|
179 | ! limitation subtile |
---|
180 | ! , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij))) |
---|
181 | |
---|
182 | ENDDO |
---|
183 | |
---|
184 | DO ij = iip1 + iip1, ip1jm, iip1 |
---|
185 | dxqmax(ij - iim, l) = dxqmax(ij, l) |
---|
186 | ENDDO |
---|
187 | |
---|
188 | DO ij = iip2 + 1, ip1jm |
---|
189 | IF(dxqu(ij - 1) * dxqu(ij)>0) THEN |
---|
190 | dxq(ij, l) = dxqu(ij - 1) + dxqu(ij) |
---|
191 | ELSE |
---|
192 | ! extremum local |
---|
193 | dxq(ij, l) = 0. |
---|
194 | ENDIF |
---|
195 | dxq(ij, l) = 0.5 * dxq(ij, l) |
---|
196 | dxq(ij, l) = & |
---|
197 | sign(min(abs(dxq(ij, l)), dxqmax(ij, l)), dxq(ij, l)) |
---|
198 | ENDDO |
---|
199 | |
---|
200 | ENDDO ! l=1,llm |
---|
201 | !print*,'Ok calcul des pentes' |
---|
202 | |
---|
203 | ELSE ! (pente_max.lt.-1.e-5) |
---|
204 | |
---|
205 | ! Pentes produits: |
---|
206 | ! ---------------- |
---|
207 | |
---|
208 | DO l = 1, llm |
---|
209 | DO ij = iip2, ip1jm - 1 |
---|
210 | dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq) |
---|
211 | ENDDO |
---|
212 | DO ij = iip1 + iip1, ip1jm, iip1 |
---|
213 | dxqu(ij) = dxqu(ij - iim) |
---|
214 | ENDDO |
---|
215 | |
---|
216 | DO ij = iip2 + 1, ip1jm |
---|
217 | zz(ij) = dxqu(ij - 1) * dxqu(ij) |
---|
218 | zz(ij) = zz(ij) + zz(ij) |
---|
219 | IF(zz(ij)>0) THEN |
---|
220 | dxq(ij, l) = zz(ij) / (dxqu(ij - 1) + dxqu(ij)) |
---|
221 | ELSE |
---|
222 | ! extremum local |
---|
223 | dxq(ij, l) = 0. |
---|
224 | ENDIF |
---|
225 | ENDDO |
---|
226 | |
---|
227 | ENDDO |
---|
228 | |
---|
229 | ENDIF ! (pente_max.lt.-1.e-5) |
---|
230 | |
---|
231 | ! bouclage de la pente en iip1: |
---|
232 | ! ----------------------------- |
---|
233 | |
---|
234 | DO l = 1, llm |
---|
235 | DO ij = iip1 + iip1, ip1jm, iip1 |
---|
236 | dxq(ij - iim, l) = dxq(ij, l) |
---|
237 | ENDDO |
---|
238 | DO ij = 1, ip1jmp1 |
---|
239 | iadvplus(ij, l) = 0 |
---|
240 | ENDDO |
---|
241 | |
---|
242 | ENDDO |
---|
243 | ! calcul des flux a gauche et a droite |
---|
244 | |
---|
245 | ! on cumule le flux correspondant a toutes les mailles dont la masse |
---|
246 | ! au travers de la paroi pENDant le pas de temps. |
---|
247 | DO l = 1, llm |
---|
248 | DO ij = iip2, ip1jm - 1 |
---|
249 | IF (u_m(ij, l)>0.) THEN |
---|
250 | zdum(ij, l) = 1. - u_m(ij, l) / masse(ij, l, iq) |
---|
251 | u_mq(ij, l) = u_m(ij, l) * (q(ij, l, iq) + 0.5 * zdum(ij, l) * dxq(ij, l)) |
---|
252 | ELSE |
---|
253 | zdum(ij, l) = 1. + u_m(ij, l) / masse(ij + 1, l, iq) |
---|
254 | u_mq(ij, l) = u_m(ij, l) * (q(ij + 1, l, iq) & |
---|
255 | - 0.5 * zdum(ij, l) * dxq(ij + 1, l)) |
---|
256 | ENDIF |
---|
257 | ENDDO |
---|
258 | ENDDO |
---|
259 | |
---|
260 | ! detection des points ou on advecte plus que la masse de la |
---|
261 | ! maille |
---|
262 | DO l = 1, llm |
---|
263 | DO ij = iip2, ip1jm - 1 |
---|
264 | IF(zdum(ij, l)<0) THEN |
---|
265 | iadvplus(ij, l) = 1 |
---|
266 | u_mq(ij, l) = 0. |
---|
267 | ENDIF |
---|
268 | ENDDO |
---|
269 | ENDDO |
---|
270 | DO l = 1, llm |
---|
271 | DO ij = iip1 + iip1, ip1jm, iip1 |
---|
272 | iadvplus(ij, l) = iadvplus(ij - iim, l) |
---|
273 | ENDDO |
---|
274 | ENDDO |
---|
275 | |
---|
276 | |
---|
277 | ! traitement special pour le cas ou on advecte en longitude plus que le |
---|
278 | ! contenu de la maille. |
---|
279 | ! cette partie est mal vectorisee. |
---|
280 | |
---|
281 | ! calcul du nombre de maille sur lequel on advecte plus que la maille. |
---|
282 | |
---|
283 | n0 = 0 |
---|
284 | DO l = 1, llm |
---|
285 | nl(l) = 0 |
---|
286 | DO ij = iip2, ip1jm |
---|
287 | nl(l) = nl(l) + iadvplus(ij, l) |
---|
288 | ENDDO |
---|
289 | n0 = n0 + nl(l) |
---|
290 | ENDDO |
---|
291 | |
---|
292 | IF(n0>0) THEN |
---|
293 | IF (prt_level > 2) PRINT *, & |
---|
294 | 'Nombre de points pour lesquels on advect plus que le' & |
---|
295 | , 'contenu de la maille : ', n0 |
---|
296 | |
---|
297 | DO l = 1, llm |
---|
298 | IF(nl(l)>0) THEN |
---|
299 | iju = 0 |
---|
300 | ! indicage des mailles concernees par le traitement special |
---|
301 | DO ij = iip2, ip1jm |
---|
302 | IF(iadvplus(ij, l)==1.AND.mod(ij, iip1)/=0) THEN |
---|
303 | iju = iju + 1 |
---|
304 | indu(iju) = ij |
---|
305 | ENDIF |
---|
306 | ENDDO |
---|
307 | niju = iju |
---|
308 | |
---|
309 | ! traitement des mailles |
---|
310 | DO iju = 1, niju |
---|
311 | ij = indu(iju) |
---|
312 | j = (ij - 1) / iip1 + 1 |
---|
313 | zu_m = u_m(ij, l) |
---|
314 | u_mq(ij, l) = 0. |
---|
315 | IF(zu_m>0.) THEN |
---|
316 | ijq = ij |
---|
317 | i = ijq - (j - 1) * iip1 |
---|
318 | ! accumulation pour les mailles completements advectees |
---|
319 | do while(zu_m>masse(ijq, l, iq)) |
---|
320 | u_mq(ij, l) = u_mq(ij, l) + q(ijq, l, iq) & |
---|
321 | * masse(ijq, l, iq) |
---|
322 | zu_m = zu_m - masse(ijq, l, iq) |
---|
323 | i = mod(i - 2 + iim, iim) + 1 |
---|
324 | ijq = (j - 1) * iip1 + i |
---|
325 | ENDDO |
---|
326 | ! ajout de la maille non completement advectee |
---|
327 | u_mq(ij, l) = u_mq(ij, l) + zu_m * & |
---|
328 | (q(ijq, l, iq) + 0.5 * (1. - zu_m / masse(ijq, l, iq)) & |
---|
329 | * dxq(ijq, l)) |
---|
330 | ELSE |
---|
331 | ijq = ij + 1 |
---|
332 | i = ijq - (j - 1) * iip1 |
---|
333 | ! accumulation pour les mailles completements advectees |
---|
334 | do while(-zu_m>masse(ijq, l, iq)) |
---|
335 | u_mq(ij, l) = u_mq(ij, l) - q(ijq, l, iq) & |
---|
336 | * masse(ijq, l, iq) |
---|
337 | zu_m = zu_m + masse(ijq, l, iq) |
---|
338 | i = mod(i, iim) + 1 |
---|
339 | ijq = (j - 1) * iip1 + i |
---|
340 | ENDDO |
---|
341 | ! ajout de la maille non completement advectee |
---|
342 | u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l, iq) - & |
---|
343 | 0.5 * (1. + zu_m / masse(ijq, l, iq)) * dxq(ijq, l)) |
---|
344 | ENDIF |
---|
345 | ENDDO |
---|
346 | ENDIF |
---|
347 | ENDDO |
---|
348 | ENDIF ! n0.gt.0 |
---|
349 | |
---|
350 | |
---|
351 | ! bouclage en latitude |
---|
352 | !print*,'cvant bouclage en latitude' |
---|
353 | DO l = 1, llm |
---|
354 | DO ij = iip1 + iip1, ip1jm, iip1 |
---|
355 | u_mq(ij, l) = u_mq(ij - iim, l) |
---|
356 | ENDDO |
---|
357 | ENDDO |
---|
358 | |
---|
359 | ! CRisi: appel récursif de l'advection sur les fils. |
---|
360 | ! Il faut faire ça avant d'avoir mis à jour q et masse |
---|
361 | !WRITE(*,*) 'vlsplt 326: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen |
---|
362 | |
---|
363 | do ifils = 1, tracers(iq)%nqDescen |
---|
364 | iq2 = tracers(iq)%iqDescen(ifils) |
---|
365 | DO l = 1, llm |
---|
366 | DO ij = iip2, ip1jm |
---|
367 | ! On a besoin de q et masse seulement entre iip2 et ip1jm |
---|
368 | !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) |
---|
369 | ! !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) |
---|
370 | !Mvals: veiller a ce qu'on n'ait pas de denominateur nul |
---|
371 | masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass) |
---|
372 | IF (q(ij, l, iq)>min_qParent) THEN |
---|
373 | Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq) |
---|
374 | else |
---|
375 | Ratio(ij, l, iq2) = min_ratio |
---|
376 | endif |
---|
377 | enddo |
---|
378 | enddo |
---|
379 | enddo |
---|
380 | do ifils = 1, tracers(iq)%nqChildren |
---|
381 | iq2 = tracers(iq)%iqDescen(ifils) |
---|
382 | CALL vlx(Ratio, pente_max, masseq, u_mq, iq2) |
---|
383 | enddo |
---|
384 | ! end CRisi |
---|
385 | |
---|
386 | |
---|
387 | ! calcul des tENDances |
---|
388 | |
---|
389 | DO l = 1, llm |
---|
390 | DO ij = iip2 + 1, ip1jm |
---|
391 | !MVals: veiller a ce qu'on ait pas de denominateur nul |
---|
392 | new_m = max(masse(ij, l, iq) + u_m(ij - 1, l) - u_m(ij, l), min_qMass) |
---|
393 | q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + & |
---|
394 | u_mq(ij - 1, l) - u_mq(ij, l)) & |
---|
395 | / new_m |
---|
396 | masse(ij, l, iq) = new_m |
---|
397 | ENDDO |
---|
398 | DO ij = iip1 + iip1, ip1jm, iip1 |
---|
399 | q(ij - iim, l, iq) = q(ij, l, iq) |
---|
400 | masse(ij - iim, l, iq) = masse(ij, l, iq) |
---|
401 | ENDDO |
---|
402 | ENDDO |
---|
403 | |
---|
404 | ! retablir les fils en rapport de melange par rapport a l'air: |
---|
405 | ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio |
---|
406 | ! puis on boucle en longitude |
---|
407 | do ifils = 1, tracers(iq)%nqDescen |
---|
408 | iq2 = tracers(iq)%iqDescen(ifils) |
---|
409 | DO l = 1, llm |
---|
410 | DO ij = iip2 + 1, ip1jm |
---|
411 | q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2) |
---|
412 | enddo |
---|
413 | DO ij = iip1 + iip1, ip1jm, iip1 |
---|
414 | q(ij - iim, l, iq2) = q(ij, l, iq2) |
---|
415 | enddo ! DO ij=ijb+iip1-1,ije,iip1 |
---|
416 | enddo !DO l=1,llm |
---|
417 | enddo |
---|
418 | |
---|
419 | END SUBROUTINE vlx |
---|
420 | RECURSIVE SUBROUTINE vly(q, pente_max, masse, masse_adv_v, iq) |
---|
421 | USE infotrac, ONLY: nqtot, tracers, & ! CRisi |
---|
422 | min_qParent, min_qMass, min_ratio ! MVals et CRisi |
---|
423 | USE lmdz_ssum_scopy, ONLY: ssum |
---|
424 | USE lmdz_comgeom |
---|
425 | |
---|
426 | ! |
---|
427 | ! Auteurs: P.Le Van, F.Hourdin, F.Forget |
---|
428 | ! |
---|
429 | ! ******************************************************************** |
---|
430 | ! Shema d'advection " pseudo amont " . |
---|
431 | ! ******************************************************************** |
---|
432 | ! q,masse_adv_v,w sont des arguments d'entree pour le s-pg .... |
---|
433 | ! dq sont des arguments de sortie pour le s-pg .... |
---|
434 | ! |
---|
435 | ! |
---|
436 | ! -------------------------------------------------------------------- |
---|
437 | USE comconst_mod, ONLY: pi |
---|
438 | IMPLICIT NONE |
---|
439 | ! |
---|
440 | INCLUDE "dimensions.h" |
---|
441 | INCLUDE "paramet.h" |
---|
442 | ! |
---|
443 | ! |
---|
444 | ! Arguments: |
---|
445 | ! ---------- |
---|
446 | REAL :: masse(ip1jmp1, llm, nqtot), pente_max |
---|
447 | REAL :: masse_adv_v(ip1jm, llm) |
---|
448 | REAL :: q(ip1jmp1, llm, nqtot) |
---|
449 | INTEGER :: iq ! CRisi |
---|
450 | ! |
---|
451 | ! Local |
---|
452 | ! --------- |
---|
453 | ! |
---|
454 | INTEGER :: i, ij, l |
---|
455 | ! |
---|
456 | REAL :: airej2, airejjm, airescb(iim), airesch(iim) |
---|
457 | REAL :: dyq(ip1jmp1, llm), dyqv(ip1jm) |
---|
458 | REAL :: adyqv(ip1jm), dyqmax(ip1jmp1) |
---|
459 | REAL :: qbyv(ip1jm, llm) |
---|
460 | |
---|
461 | REAL :: qpns, qpsn, dyn1, dys1, dyn2, dys2, newmasse, fn, fs |
---|
462 | LOGICAL, SAVE :: first |
---|
463 | |
---|
464 | REAL :: convpn, convps, convmpn, convmps |
---|
465 | REAL :: massepn, masseps, qpn, qps |
---|
466 | REAL :: sinlon(iip1), sinlondlon(iip1) |
---|
467 | REAL :: coslon(iip1), coslondlon(iip1) |
---|
468 | SAVE sinlon, coslon, sinlondlon, coslondlon |
---|
469 | SAVE airej2, airejjm |
---|
470 | |
---|
471 | REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) ! CRisi |
---|
472 | INTEGER :: ifils, iq2 ! CRisi |
---|
473 | |
---|
474 | DATA first/.TRUE./ |
---|
475 | |
---|
476 | ! !WRITE(*,*) 'vly 578: entree, iq=',iq |
---|
477 | |
---|
478 | IF(first) THEN |
---|
479 | PRINT*, 'Shema Amont nouveau appele dans Vanleer ' |
---|
480 | first = .FALSE. |
---|
481 | do i = 2, iip1 |
---|
482 | coslon(i) = cos(rlonv(i)) |
---|
483 | sinlon(i) = sin(rlonv(i)) |
---|
484 | coslondlon(i) = coslon(i) * (rlonu(i) - rlonu(i - 1)) / pi |
---|
485 | sinlondlon(i) = sinlon(i) * (rlonu(i) - rlonu(i - 1)) / pi |
---|
486 | ENDDO |
---|
487 | coslon(1) = coslon(iip1) |
---|
488 | coslondlon(1) = coslondlon(iip1) |
---|
489 | sinlon(1) = sinlon(iip1) |
---|
490 | sinlondlon(1) = sinlondlon(iip1) |
---|
491 | airej2 = SSUM(iim, aire(iip2), 1) |
---|
492 | airejjm = SSUM(iim, aire(ip1jm - iim), 1) |
---|
493 | ENDIF |
---|
494 | |
---|
495 | ! |
---|
496 | !PRINT*,'CALCUL EN LATITUDE' |
---|
497 | |
---|
498 | DO l = 1, llm |
---|
499 | ! |
---|
500 | ! -------------------------------- |
---|
501 | ! CALCUL EN LATITUDE |
---|
502 | ! -------------------------------- |
---|
503 | |
---|
504 | ! On commence par calculer la valeur du traceur moyenne sur le premier cercle |
---|
505 | ! de latitude autour du pole (qpns pour le pole nord et qpsn pour |
---|
506 | ! le pole nord) qui sera utilisee pour evaluer les pentes au pole. |
---|
507 | |
---|
508 | DO i = 1, iim |
---|
509 | airescb(i) = aire(i + iip1) * q(i + iip1, l, iq) |
---|
510 | airesch(i) = aire(i + ip1jm - iip1) * q(i + ip1jm - iip1, l, iq) |
---|
511 | ENDDO |
---|
512 | qpns = SSUM(iim, airescb, 1) / airej2 |
---|
513 | qpsn = SSUM(iim, airesch, 1) / airejjm |
---|
514 | |
---|
515 | ! calcul des pentes aux points v |
---|
516 | |
---|
517 | DO ij = 1, ip1jm |
---|
518 | dyqv(ij) = q(ij, l, iq) - q(ij + iip1, l, iq) |
---|
519 | adyqv(ij) = abs(dyqv(ij)) |
---|
520 | ENDDO |
---|
521 | |
---|
522 | ! calcul des pentes aux points scalaires |
---|
523 | |
---|
524 | DO ij = iip2, ip1jm |
---|
525 | dyq(ij, l) = .5 * (dyqv(ij - iip1) + dyqv(ij)) |
---|
526 | dyqmax(ij) = min(adyqv(ij - iip1), adyqv(ij)) |
---|
527 | dyqmax(ij) = pente_max * dyqmax(ij) |
---|
528 | ENDDO |
---|
529 | |
---|
530 | ! calcul des pentes aux poles |
---|
531 | |
---|
532 | DO ij = 1, iip1 |
---|
533 | dyq(ij, l) = qpns - q(ij + iip1, l, iq) |
---|
534 | dyq(ip1jm + ij, l) = q(ip1jm + ij - iip1, l, iq) - qpsn |
---|
535 | ENDDO |
---|
536 | |
---|
537 | ! filtrage de la derivee |
---|
538 | dyn1 = 0. |
---|
539 | dys1 = 0. |
---|
540 | dyn2 = 0. |
---|
541 | dys2 = 0. |
---|
542 | DO ij = 1, iim |
---|
543 | dyn1 = dyn1 + sinlondlon(ij) * dyq(ij, l) |
---|
544 | dys1 = dys1 + sinlondlon(ij) * dyq(ip1jm + ij, l) |
---|
545 | dyn2 = dyn2 + coslondlon(ij) * dyq(ij, l) |
---|
546 | dys2 = dys2 + coslondlon(ij) * dyq(ip1jm + ij, l) |
---|
547 | ENDDO |
---|
548 | DO ij = 1, iip1 |
---|
549 | dyq(ij, l) = dyn1 * sinlon(ij) + dyn2 * coslon(ij) |
---|
550 | dyq(ip1jm + ij, l) = dys1 * sinlon(ij) + dys2 * coslon(ij) |
---|
551 | ENDDO |
---|
552 | |
---|
553 | ! calcul des pentes limites aux poles |
---|
554 | |
---|
555 | goto 8888 |
---|
556 | fn = 1. |
---|
557 | fs = 1. |
---|
558 | DO ij = 1, iim |
---|
559 | IF(pente_max * adyqv(ij)<abs(dyq(ij, l))) THEN |
---|
560 | fn = min(pente_max * adyqv(ij) / abs(dyq(ij, l)), fn) |
---|
561 | ENDIF |
---|
562 | IF(pente_max * adyqv(ij + ip1jm - iip1)<abs(dyq(ij + ip1jm, l))) THEN |
---|
563 | fs = min(pente_max * adyqv(ij + ip1jm - iip1) / abs(dyq(ij + ip1jm, l)), fs) |
---|
564 | ENDIF |
---|
565 | ENDDO |
---|
566 | DO ij = 1, iip1 |
---|
567 | dyq(ij, l) = fn * dyq(ij, l) |
---|
568 | dyq(ip1jm + ij, l) = fs * dyq(ip1jm + ij, l) |
---|
569 | ENDDO |
---|
570 | 8888 continue |
---|
571 | DO ij = 1, iip1 |
---|
572 | dyq(ij, l) = 0. |
---|
573 | dyq(ip1jm + ij, l) = 0. |
---|
574 | ENDDO |
---|
575 | |
---|
576 | !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
---|
577 | ! En memoire de dIFferents tests sur la |
---|
578 | ! limitation des pentes aux poles. |
---|
579 | !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
---|
580 | ! PRINT*,dyq(1) |
---|
581 | ! PRINT*,dyqv(iip1+1) |
---|
582 | ! appn=abs(dyq(1)/dyqv(iip1+1)) |
---|
583 | ! PRINT*,dyq(ip1jm+1) |
---|
584 | ! PRINT*,dyqv(ip1jm-iip1+1) |
---|
585 | ! apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1)) |
---|
586 | ! DO ij=2,iim |
---|
587 | ! appn=amax1(abs(dyq(ij)/dyqv(ij)),appn) |
---|
588 | ! apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps) |
---|
589 | ! ENDDO |
---|
590 | ! appn=min(pente_max/appn,1.) |
---|
591 | ! apps=min(pente_max/apps,1.) |
---|
592 | ! |
---|
593 | ! |
---|
594 | ! cas ou on a un extremum au pole |
---|
595 | ! |
---|
596 | ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) |
---|
597 | ! & appn=0. |
---|
598 | ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* |
---|
599 | ! & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) |
---|
600 | ! & apps=0. |
---|
601 | ! |
---|
602 | ! limitation des pentes aux poles |
---|
603 | ! DO ij=1,iip1 |
---|
604 | ! dyq(ij)=appn*dyq(ij) |
---|
605 | ! dyq(ip1jm+ij)=apps*dyq(ip1jm+ij) |
---|
606 | ! ENDDO |
---|
607 | ! |
---|
608 | ! test |
---|
609 | ! DO ij=1,iip1 |
---|
610 | ! dyq(iip1+ij)=0. |
---|
611 | ! dyq(ip1jm+ij-iip1)=0. |
---|
612 | ! ENDDO |
---|
613 | ! DO ij=1,ip1jmp1 |
---|
614 | ! dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1)) |
---|
615 | ! ENDDO |
---|
616 | ! |
---|
617 | ! changement 10 07 96 |
---|
618 | ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.) |
---|
619 | ! & THEN |
---|
620 | ! DO ij=1,iip1 |
---|
621 | ! dyqmax(ij)=0. |
---|
622 | ! ENDDO |
---|
623 | ! ELSE |
---|
624 | ! DO ij=1,iip1 |
---|
625 | ! dyqmax(ij)=pente_max*abs(dyqv(ij)) |
---|
626 | ! ENDDO |
---|
627 | ! ENDIF |
---|
628 | ! |
---|
629 | ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* |
---|
630 | ! & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.) |
---|
631 | ! &THEN |
---|
632 | ! DO ij=ip1jm+1,ip1jmp1 |
---|
633 | ! dyqmax(ij)=0. |
---|
634 | ! ENDDO |
---|
635 | ! ELSE |
---|
636 | ! DO ij=ip1jm+1,ip1jmp1 |
---|
637 | ! dyqmax(ij)=pente_max*abs(dyqv(ij-iip1)) |
---|
638 | ! ENDDO |
---|
639 | ! ENDIF |
---|
640 | ! fin changement 10 07 96 |
---|
641 | !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
---|
642 | |
---|
643 | ! calcul des pentes limitees |
---|
644 | |
---|
645 | DO ij = iip2, ip1jm |
---|
646 | IF(dyqv(ij) * dyqv(ij - iip1)>0.) THEN |
---|
647 | dyq(ij, l) = sign(min(abs(dyq(ij, l)), dyqmax(ij)), dyq(ij, l)) |
---|
648 | ELSE |
---|
649 | dyq(ij, l) = 0. |
---|
650 | ENDIF |
---|
651 | ENDDO |
---|
652 | |
---|
653 | ENDDO |
---|
654 | |
---|
655 | ! !WRITE(*,*) 'vly 756' |
---|
656 | DO l = 1, llm |
---|
657 | DO ij = 1, ip1jm |
---|
658 | IF(masse_adv_v(ij, l)>0) THEN |
---|
659 | qbyv(ij, l) = q(ij + iip1, l, iq) + dyq(ij + iip1, l) * & |
---|
660 | 0.5 * (1. - masse_adv_v(ij, l) & |
---|
661 | / masse(ij + iip1, l, iq)) |
---|
662 | ELSE |
---|
663 | qbyv(ij, l) = q(ij, l, iq) - dyq(ij, l) * & |
---|
664 | 0.5 * (1. + masse_adv_v(ij, l) & |
---|
665 | / masse(ij, l, iq)) |
---|
666 | ENDIF |
---|
667 | qbyv(ij, l) = masse_adv_v(ij, l) * qbyv(ij, l) |
---|
668 | ENDDO |
---|
669 | ENDDO |
---|
670 | |
---|
671 | ! CRisi: appel récursif de l'advection sur les fils. |
---|
672 | ! Il faut faire ça avant d'avoir mis à jour q et masse |
---|
673 | ! WRITE(*,*) 'vly 689: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen |
---|
674 | |
---|
675 | do ifils = 1, tracers(iq)%nqDescen |
---|
676 | iq2 = tracers(iq)%iqDescen(ifils) |
---|
677 | DO l = 1, llm |
---|
678 | DO ij = 1, ip1jmp1 |
---|
679 | ! ! attention, chaque fils doit avoir son masseq, sinon, le 1er |
---|
680 | ! ! fils ecrase le masseq de ses freres. |
---|
681 | ! !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) |
---|
682 | ! !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) |
---|
683 | ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul |
---|
684 | masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass) |
---|
685 | IF (q(ij, l, iq)>min_qParent) THEN |
---|
686 | Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq) |
---|
687 | else |
---|
688 | Ratio(ij, l, iq2) = min_ratio |
---|
689 | endif |
---|
690 | enddo |
---|
691 | enddo |
---|
692 | enddo |
---|
693 | |
---|
694 | do ifils = 1, tracers(iq)%nqDescen |
---|
695 | iq2 = tracers(iq)%iqDescen(ifils) |
---|
696 | CALL vly(Ratio, pente_max, masseq, qbyv, iq2) |
---|
697 | enddo |
---|
698 | |
---|
699 | DO l = 1, llm |
---|
700 | DO ij = iip2, ip1jm |
---|
701 | newmasse = masse(ij, l, iq) & |
---|
702 | + masse_adv_v(ij, l) - masse_adv_v(ij - iip1, l) |
---|
703 | q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + qbyv(ij, l) & |
---|
704 | - qbyv(ij - iip1, l)) / newmasse |
---|
705 | masse(ij, l, iq) = newmasse |
---|
706 | ENDDO |
---|
707 | convpn = SSUM(iim, qbyv(1, l), 1) |
---|
708 | convmpn = ssum(iim, masse_adv_v(1, l), 1) |
---|
709 | massepn = ssum(iim, masse(1, l, iq), 1) |
---|
710 | qpn = 0. |
---|
711 | do ij = 1, iim |
---|
712 | qpn = qpn + masse(ij, l, iq) * q(ij, l, iq) |
---|
713 | enddo |
---|
714 | qpn = (qpn + convpn) / (massepn + convmpn) |
---|
715 | do ij = 1, iip1 |
---|
716 | q(ij, l, iq) = qpn |
---|
717 | enddo |
---|
718 | convps = -SSUM(iim, qbyv(ip1jm - iim, l), 1) |
---|
719 | convmps = -ssum(iim, masse_adv_v(ip1jm - iim, l), 1) |
---|
720 | masseps = ssum(iim, masse(ip1jm + 1, l, iq), 1) |
---|
721 | qps = 0. |
---|
722 | do ij = ip1jm + 1, ip1jmp1 - 1 |
---|
723 | qps = qps + masse(ij, l, iq) * q(ij, l, iq) |
---|
724 | enddo |
---|
725 | qps = (qps + convps) / (masseps + convmps) |
---|
726 | do ij = ip1jm + 1, ip1jmp1 |
---|
727 | q(ij, l, iq) = qps |
---|
728 | enddo |
---|
729 | ENDDO |
---|
730 | |
---|
731 | ! retablir les fils en rapport de melange par rapport a l'air: |
---|
732 | do ifils = 1, tracers(iq)%nqDescen |
---|
733 | iq2 = tracers(iq)%iqDescen(ifils) |
---|
734 | DO l = 1, llm |
---|
735 | DO ij = 1, ip1jmp1 |
---|
736 | q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2) |
---|
737 | enddo |
---|
738 | enddo |
---|
739 | enddo |
---|
740 | |
---|
741 | ! !WRITE(*,*) 'vly 853: sortie' |
---|
742 | |
---|
743 | END SUBROUTINE vly |
---|
744 | RECURSIVE SUBROUTINE vlz(q, pente_max, masse, w, iq) |
---|
745 | USE infotrac, ONLY: nqtot, tracers, & ! CRisi |
---|
746 | min_qParent, min_qMass, min_ratio ! MVals et CRisi |
---|
747 | ! |
---|
748 | ! Auteurs: P.Le Van, F.Hourdin, F.Forget |
---|
749 | ! |
---|
750 | ! ******************************************************************** |
---|
751 | ! Shema d'advection " pseudo amont " . |
---|
752 | ! ******************************************************************** |
---|
753 | ! q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg .... |
---|
754 | ! dq sont des arguments de sortie pour le s-pg .... |
---|
755 | ! -------------------------------------------------------------------- |
---|
756 | IMPLICIT NONE |
---|
757 | ! |
---|
758 | INCLUDE "dimensions.h" |
---|
759 | INCLUDE "paramet.h" |
---|
760 | ! |
---|
761 | ! |
---|
762 | ! Arguments: |
---|
763 | ! ---------- |
---|
764 | REAL :: masse(ip1jmp1, llm, nqtot), pente_max |
---|
765 | REAL :: q(ip1jmp1, llm, nqtot) |
---|
766 | REAL :: w(ip1jmp1, llm + 1) |
---|
767 | INTEGER :: iq |
---|
768 | ! |
---|
769 | ! Local |
---|
770 | ! --------- |
---|
771 | ! |
---|
772 | INTEGER :: ij, l |
---|
773 | ! |
---|
774 | REAL :: wq(ip1jmp1, llm + 1), newmasse |
---|
775 | |
---|
776 | REAL :: dzq(ip1jmp1, llm), dzqw(ip1jmp1, llm), adzqw(ip1jmp1, llm), dzqmax |
---|
777 | REAL :: sigw |
---|
778 | |
---|
779 | REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) ! CRisi |
---|
780 | INTEGER :: ifils, iq2 ! CRisi |
---|
781 | |
---|
782 | #ifdef BIDON |
---|
783 | REAL :: temps0,temps1,second |
---|
784 | SAVE temps0,temps1 |
---|
785 | |
---|
786 | DATA temps0,temps1/0.,0./ |
---|
787 | #endif |
---|
788 | |
---|
789 | ! On oriente tout dans le sens de la pression c'est a dire dans le |
---|
790 | ! sens de W |
---|
791 | DO l = 2, llm |
---|
792 | DO ij = 1, ip1jmp1 |
---|
793 | dzqw(ij, l) = q(ij, l - 1, iq) - q(ij, l, iq) |
---|
794 | adzqw(ij, l) = abs(dzqw(ij, l)) |
---|
795 | ENDDO |
---|
796 | ENDDO |
---|
797 | |
---|
798 | DO l = 2, llm - 1 |
---|
799 | DO ij = 1, ip1jmp1 |
---|
800 | IF(dzqw(ij, l) * dzqw(ij, l + 1)>0.) THEN |
---|
801 | dzq(ij, l) = 0.5 * (dzqw(ij, l) + dzqw(ij, l + 1)) |
---|
802 | ELSE |
---|
803 | dzq(ij, l) = 0. |
---|
804 | ENDIF |
---|
805 | dzqmax = pente_max * min(adzqw(ij, l), adzqw(ij, l + 1)) |
---|
806 | dzq(ij, l) = sign(min(abs(dzq(ij, l)), dzqmax), dzq(ij, l)) |
---|
807 | ENDDO |
---|
808 | ENDDO |
---|
809 | |
---|
810 | ! !WRITE(*,*) 'vlz 954' |
---|
811 | DO ij = 1, ip1jmp1 |
---|
812 | dzq(ij, 1) = 0. |
---|
813 | dzq(ij, llm) = 0. |
---|
814 | ENDDO |
---|
815 | |
---|
816 | ! --------------------------------------------------------------- |
---|
817 | ! .... calcul des termes d'advection verticale ....... |
---|
818 | ! --------------------------------------------------------------- |
---|
819 | |
---|
820 | ! calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq |
---|
821 | |
---|
822 | DO l = 1, llm - 1 |
---|
823 | do ij = 1, ip1jmp1 |
---|
824 | IF(w(ij, l + 1)>0.) THEN |
---|
825 | sigw = w(ij, l + 1) / masse(ij, l + 1, iq) |
---|
826 | wq(ij, l + 1) = w(ij, l + 1) * (q(ij, l + 1, iq) & |
---|
827 | + 0.5 * (1. - sigw) * dzq(ij, l + 1)) |
---|
828 | ELSE |
---|
829 | sigw = w(ij, l + 1) / masse(ij, l, iq) |
---|
830 | wq(ij, l + 1) = w(ij, l + 1) * (q(ij, l, iq) - 0.5 * (1. + sigw) * dzq(ij, l)) |
---|
831 | ENDIF |
---|
832 | ENDDO |
---|
833 | ENDDO |
---|
834 | |
---|
835 | DO ij = 1, ip1jmp1 |
---|
836 | wq(ij, llm + 1) = 0. |
---|
837 | wq(ij, 1) = 0. |
---|
838 | ENDDO |
---|
839 | |
---|
840 | ! CRisi: appel récursif de l'advection sur les fils. |
---|
841 | ! Il faut faire ça avant d'avoir mis à jour q et masse |
---|
842 | ! !WRITE(*,*) 'vlsplt 942: iq,nqChildren(iq)=',iq,nqChildren(iq) |
---|
843 | do ifils = 1, tracers(iq)%nqDescen |
---|
844 | iq2 = tracers(iq)%iqDescen(ifils) |
---|
845 | DO l = 1, llm |
---|
846 | DO ij = 1, ip1jmp1 |
---|
847 | ! !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) |
---|
848 | ! !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) |
---|
849 | ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul |
---|
850 | masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass) |
---|
851 | IF (q(ij, l, iq)>min_qParent) THEN |
---|
852 | Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq) |
---|
853 | else |
---|
854 | Ratio(ij, l, iq2) = min_ratio |
---|
855 | endif |
---|
856 | enddo |
---|
857 | enddo |
---|
858 | enddo |
---|
859 | |
---|
860 | do ifils = 1, tracers(iq)%nqChildren |
---|
861 | iq2 = tracers(iq)%iqDescen(ifils) |
---|
862 | CALL vlz(Ratio, pente_max, masseq, wq, iq2) |
---|
863 | enddo |
---|
864 | ! end CRisi |
---|
865 | |
---|
866 | DO l = 1, llm |
---|
867 | DO ij = 1, ip1jmp1 |
---|
868 | newmasse = masse(ij, l, iq) + w(ij, l + 1) - w(ij, l) |
---|
869 | q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + wq(ij, l + 1) - wq(ij, l)) & |
---|
870 | / newmasse |
---|
871 | masse(ij, l, iq) = newmasse |
---|
872 | ENDDO |
---|
873 | ENDDO |
---|
874 | |
---|
875 | ! retablir les fils en rapport de melange par rapport a l'air: |
---|
876 | do ifils = 1, tracers(iq)%nqDescen |
---|
877 | iq2 = tracers(iq)%iqDescen(ifils) |
---|
878 | DO l = 1, llm |
---|
879 | DO ij = 1, ip1jmp1 |
---|
880 | q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2) |
---|
881 | enddo |
---|
882 | enddo |
---|
883 | enddo |
---|
884 | |
---|
885 | END SUBROUTINE vlz |
---|
886 | |
---|
887 | SUBROUTINE minmaxq(zq, qmin, qmax, comment) |
---|
888 | |
---|
889 | INCLUDE "dimensions.h" |
---|
890 | INCLUDE "paramet.h" |
---|
891 | |
---|
892 | CHARACTER(LEN = 20) :: comment |
---|
893 | REAL :: qmin, qmax |
---|
894 | REAL :: zq(ip1jmp1, llm) |
---|
895 | REAL :: zzq(iip1, jjp1, llm) |
---|
896 | |
---|
897 | END SUBROUTINE minmaxq |
---|
898 | |
---|
899 | |
---|
900 | |
---|