1 | |
---|
2 | ! $Id: thermcell.F90 1992 2014-03-05 13:19:12Z abarral $ |
---|
3 | |
---|
4 | SUBROUTINE calcul_sec(ngrid, nlay, ptimestep, pplay, pplev, pphi, zlev, pu, & |
---|
5 | pv, pt, po, zmax, wmax, zw2, lmix & ! s |
---|
6 | ! ,pu_therm,pv_therm |
---|
7 | , r_aspect, l_mix, w2di, tho) |
---|
8 | |
---|
9 | USE dimphy |
---|
10 | IMPLICIT NONE |
---|
11 | |
---|
12 | ! ======================================================================= |
---|
13 | |
---|
14 | ! Calcul du transport verticale dans la couche limite en presence |
---|
15 | ! de "thermiques" explicitement representes |
---|
16 | |
---|
17 | ! Réécriture à partir d'un listing papier à Habas, le 14/02/00 |
---|
18 | |
---|
19 | ! le thermique est supposé homogène et dissipé par mélange avec |
---|
20 | ! son environnement. la longueur l_mix contrôle l'efficacité du |
---|
21 | ! mélange |
---|
22 | |
---|
23 | ! Le calcul du transport des différentes espèces se fait en prenant |
---|
24 | ! en compte: |
---|
25 | ! 1. un flux de masse montant |
---|
26 | ! 2. un flux de masse descendant |
---|
27 | ! 3. un entrainement |
---|
28 | ! 4. un detrainement |
---|
29 | |
---|
30 | ! ======================================================================= |
---|
31 | |
---|
32 | ! ----------------------------------------------------------------------- |
---|
33 | ! declarations: |
---|
34 | ! ------------- |
---|
35 | |
---|
36 | include "dimensions.h" |
---|
37 | ! ccc#include "dimphy.h" |
---|
38 | include "YOMCST.h" |
---|
39 | |
---|
40 | ! arguments: |
---|
41 | ! ---------- |
---|
42 | |
---|
43 | INTEGER ngrid, nlay, w2di |
---|
44 | REAL tho |
---|
45 | REAL ptimestep, l_mix, r_aspect |
---|
46 | REAL pt(ngrid, nlay), pdtadj(ngrid, nlay) |
---|
47 | REAL pu(ngrid, nlay), pduadj(ngrid, nlay) |
---|
48 | REAL pv(ngrid, nlay), pdvadj(ngrid, nlay) |
---|
49 | REAL po(ngrid, nlay), pdoadj(ngrid, nlay) |
---|
50 | REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1) |
---|
51 | REAL pphi(ngrid, nlay) |
---|
52 | |
---|
53 | INTEGER idetr |
---|
54 | SAVE idetr |
---|
55 | DATA idetr/3/ |
---|
56 | !$OMP THREADPRIVATE(idetr) |
---|
57 | ! local: |
---|
58 | ! ------ |
---|
59 | |
---|
60 | INTEGER ig, k, l, lmaxa(klon), lmix(klon) |
---|
61 | REAL zsortie1d(klon) |
---|
62 | ! CR: on remplace lmax(klon,klev+1) |
---|
63 | INTEGER lmax(klon), lmin(klon), lentr(klon) |
---|
64 | REAL linter(klon) |
---|
65 | REAL zmix(klon), fracazmix(klon) |
---|
66 | ! RC |
---|
67 | REAL zmax(klon), zw, zw2(klon, klev+1), ztva(klon, klev) |
---|
68 | |
---|
69 | REAL zlev(klon, klev+1), zlay(klon, klev) |
---|
70 | REAL zh(klon, klev), zdhadj(klon, klev) |
---|
71 | REAL ztv(klon, klev) |
---|
72 | REAL zu(klon, klev), zv(klon, klev), zo(klon, klev) |
---|
73 | REAL wh(klon, klev+1) |
---|
74 | REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1) |
---|
75 | REAL zla(klon, klev+1) |
---|
76 | REAL zwa(klon, klev+1) |
---|
77 | REAL zld(klon, klev+1) |
---|
78 | ! real zwd(klon,klev+1) |
---|
79 | REAL zsortie(klon, klev) |
---|
80 | REAL zva(klon, klev) |
---|
81 | REAL zua(klon, klev) |
---|
82 | REAL zoa(klon, klev) |
---|
83 | |
---|
84 | REAL zha(klon, klev) |
---|
85 | REAL wa_moy(klon, klev+1) |
---|
86 | REAL fraca(klon, klev+1) |
---|
87 | REAL fracc(klon, klev+1) |
---|
88 | REAL zf, zf2 |
---|
89 | REAL thetath2(klon, klev), wth2(klon, klev) |
---|
90 | ! common/comtherm/thetath2,wth2 |
---|
91 | |
---|
92 | REAL count_time |
---|
93 | ! integer isplit,nsplit |
---|
94 | INTEGER isplit, nsplit, ialt |
---|
95 | PARAMETER (nsplit=10) |
---|
96 | DATA isplit/0/ |
---|
97 | SAVE isplit |
---|
98 | !$OMP THREADPRIVATE(isplit) |
---|
99 | |
---|
100 | LOGICAL sorties |
---|
101 | REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev) |
---|
102 | REAL zpspsk(klon, klev) |
---|
103 | |
---|
104 | ! real wmax(klon,klev),wmaxa(klon) |
---|
105 | REAL wmax(klon), wmaxa(klon) |
---|
106 | REAL wa(klon, klev, klev+1) |
---|
107 | REAL wd(klon, klev+1) |
---|
108 | REAL larg_part(klon, klev, klev+1) |
---|
109 | REAL fracd(klon, klev+1) |
---|
110 | REAL xxx(klon, klev+1) |
---|
111 | REAL larg_cons(klon, klev+1) |
---|
112 | REAL larg_detr(klon, klev+1) |
---|
113 | REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev) |
---|
114 | REAL pu_therm(klon, klev), pv_therm(klon, klev) |
---|
115 | REAL fm(klon, klev+1), entr(klon, klev) |
---|
116 | REAL fmc(klon, klev+1) |
---|
117 | |
---|
118 | ! CR:nouvelles variables |
---|
119 | REAL f_star(klon, klev+1), entr_star(klon, klev) |
---|
120 | REAL entr_star_tot(klon), entr_star2(klon) |
---|
121 | REAL zalim(klon) |
---|
122 | INTEGER lalim(klon) |
---|
123 | REAL norme(klon) |
---|
124 | REAL f(klon), f0(klon) |
---|
125 | REAL zlevinter(klon) |
---|
126 | LOGICAL therm |
---|
127 | LOGICAL first |
---|
128 | DATA first/.FALSE./ |
---|
129 | SAVE first |
---|
130 | !$OMP THREADPRIVATE(first) |
---|
131 | ! RC |
---|
132 | |
---|
133 | CHARACTER *2 str2 |
---|
134 | CHARACTER *10 str10 |
---|
135 | |
---|
136 | CHARACTER (LEN=20) :: modname = 'calcul_sec' |
---|
137 | CHARACTER (LEN=80) :: abort_message |
---|
138 | |
---|
139 | |
---|
140 | ! LOGICAL vtest(klon),down |
---|
141 | |
---|
142 | EXTERNAL scopy |
---|
143 | |
---|
144 | INTEGER ncorrec |
---|
145 | SAVE ncorrec |
---|
146 | DATA ncorrec/0/ |
---|
147 | !$OMP THREADPRIVATE(ncorrec) |
---|
148 | |
---|
149 | |
---|
150 | ! ----------------------------------------------------------------------- |
---|
151 | ! initialisation: |
---|
152 | ! --------------- |
---|
153 | |
---|
154 | sorties = .TRUE. |
---|
155 | IF (ngrid/=klon) THEN |
---|
156 | PRINT * |
---|
157 | PRINT *, 'STOP dans convadj' |
---|
158 | PRINT *, 'ngrid =', ngrid |
---|
159 | PRINT *, 'klon =', klon |
---|
160 | END IF |
---|
161 | |
---|
162 | ! ----------------------------------------------------------------------- |
---|
163 | ! incrementation eventuelle de tendances precedentes: |
---|
164 | ! --------------------------------------------------- |
---|
165 | |
---|
166 | ! print*,'0 OK convect8' |
---|
167 | |
---|
168 | DO l = 1, nlay |
---|
169 | DO ig = 1, ngrid |
---|
170 | zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa |
---|
171 | zh(ig, l) = pt(ig, l)/zpspsk(ig, l) |
---|
172 | zu(ig, l) = pu(ig, l) |
---|
173 | zv(ig, l) = pv(ig, l) |
---|
174 | zo(ig, l) = po(ig, l) |
---|
175 | ztv(ig, l) = zh(ig, l)*(1.+0.61*zo(ig,l)) |
---|
176 | END DO |
---|
177 | END DO |
---|
178 | |
---|
179 | ! print*,'1 OK convect8' |
---|
180 | ! -------------------- |
---|
181 | |
---|
182 | |
---|
183 | ! + + + + + + + + + + + |
---|
184 | |
---|
185 | |
---|
186 | ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz |
---|
187 | ! wh,wt,wo ... |
---|
188 | |
---|
189 | ! + + + + + + + + + + + zh,zu,zv,zo,rho |
---|
190 | |
---|
191 | |
---|
192 | ! -------------------- zlev(1) |
---|
193 | ! \\\\\\\\\\\\\\\\\\\\ |
---|
194 | |
---|
195 | |
---|
196 | |
---|
197 | ! ----------------------------------------------------------------------- |
---|
198 | ! Calcul des altitudes des couches |
---|
199 | ! ----------------------------------------------------------------------- |
---|
200 | |
---|
201 | DO l = 2, nlay |
---|
202 | DO ig = 1, ngrid |
---|
203 | zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg |
---|
204 | END DO |
---|
205 | END DO |
---|
206 | DO ig = 1, ngrid |
---|
207 | zlev(ig, 1) = 0. |
---|
208 | zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg |
---|
209 | END DO |
---|
210 | DO l = 1, nlay |
---|
211 | DO ig = 1, ngrid |
---|
212 | zlay(ig, l) = pphi(ig, l)/rg |
---|
213 | END DO |
---|
214 | END DO |
---|
215 | |
---|
216 | ! print*,'2 OK convect8' |
---|
217 | ! ----------------------------------------------------------------------- |
---|
218 | ! Calcul des densites |
---|
219 | ! ----------------------------------------------------------------------- |
---|
220 | |
---|
221 | DO l = 1, nlay |
---|
222 | DO ig = 1, ngrid |
---|
223 | rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*zh(ig,l)) |
---|
224 | END DO |
---|
225 | END DO |
---|
226 | |
---|
227 | DO l = 2, nlay |
---|
228 | DO ig = 1, ngrid |
---|
229 | rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1)) |
---|
230 | END DO |
---|
231 | END DO |
---|
232 | |
---|
233 | DO k = 1, nlay |
---|
234 | DO l = 1, nlay + 1 |
---|
235 | DO ig = 1, ngrid |
---|
236 | wa(ig, k, l) = 0. |
---|
237 | END DO |
---|
238 | END DO |
---|
239 | END DO |
---|
240 | |
---|
241 | ! print*,'3 OK convect8' |
---|
242 | ! ------------------------------------------------------------------ |
---|
243 | ! Calcul de w2, quarre de w a partir de la cape |
---|
244 | ! a partir de w2, on calcule wa, vitesse de l'ascendance |
---|
245 | |
---|
246 | ! ATTENTION: Dans cette version, pour cause d'economie de memoire, |
---|
247 | ! w2 est stoke dans wa |
---|
248 | |
---|
249 | ! ATTENTION: dans convect8, on n'utilise le calcule des wa |
---|
250 | ! independants par couches que pour calculer l'entrainement |
---|
251 | ! a la base et la hauteur max de l'ascendance. |
---|
252 | |
---|
253 | ! Indicages: |
---|
254 | ! l'ascendance provenant du niveau k traverse l'interface l avec |
---|
255 | ! une vitesse wa(k,l). |
---|
256 | |
---|
257 | ! -------------------- |
---|
258 | |
---|
259 | ! + + + + + + + + + + |
---|
260 | |
---|
261 | ! wa(k,l) ---- -------------------- l |
---|
262 | ! /\ |
---|
263 | ! /||\ + + + + + + + + + + |
---|
264 | ! || |
---|
265 | ! || -------------------- |
---|
266 | ! || |
---|
267 | ! || + + + + + + + + + + |
---|
268 | ! || |
---|
269 | ! || -------------------- |
---|
270 | ! ||__ |
---|
271 | ! |___ + + + + + + + + + + k |
---|
272 | |
---|
273 | ! -------------------- |
---|
274 | |
---|
275 | |
---|
276 | |
---|
277 | ! ------------------------------------------------------------------ |
---|
278 | |
---|
279 | ! CR: ponderation entrainement des couches instables |
---|
280 | ! def des entr_star tels que entr=f*entr_star |
---|
281 | DO l = 1, klev |
---|
282 | DO ig = 1, ngrid |
---|
283 | entr_star(ig, l) = 0. |
---|
284 | END DO |
---|
285 | END DO |
---|
286 | ! determination de la longueur de la couche d entrainement |
---|
287 | DO ig = 1, ngrid |
---|
288 | lentr(ig) = 1 |
---|
289 | END DO |
---|
290 | |
---|
291 | ! on ne considere que les premieres couches instables |
---|
292 | therm = .FALSE. |
---|
293 | DO k = nlay - 2, 1, -1 |
---|
294 | DO ig = 1, ngrid |
---|
295 | IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN |
---|
296 | lentr(ig) = k + 1 |
---|
297 | therm = .TRUE. |
---|
298 | END IF |
---|
299 | END DO |
---|
300 | END DO |
---|
301 | ! limitation de la valeur du lentr |
---|
302 | ! do ig=1,ngrid |
---|
303 | ! lentr(ig)=min(5,lentr(ig)) |
---|
304 | ! enddo |
---|
305 | ! determination du lmin: couche d ou provient le thermique |
---|
306 | DO ig = 1, ngrid |
---|
307 | lmin(ig) = 1 |
---|
308 | END DO |
---|
309 | DO ig = 1, ngrid |
---|
310 | DO l = nlay, 2, -1 |
---|
311 | IF (ztv(ig,l-1)>ztv(ig,l)) THEN |
---|
312 | lmin(ig) = l - 1 |
---|
313 | END IF |
---|
314 | END DO |
---|
315 | END DO |
---|
316 | ! initialisations |
---|
317 | DO ig = 1, ngrid |
---|
318 | zalim(ig) = 0. |
---|
319 | norme(ig) = 0. |
---|
320 | lalim(ig) = 1 |
---|
321 | END DO |
---|
322 | DO k = 1, klev - 1 |
---|
323 | DO ig = 1, ngrid |
---|
324 | zalim(ig) = zalim(ig) + zlev(ig, k)*max(0., (ztv(ig,k)-ztv(ig, & |
---|
325 | k+1))/(zlev(ig,k+1)-zlev(ig,k))) |
---|
326 | ! s *(zlev(ig,k+1)-zlev(ig,k)) |
---|
327 | norme(ig) = norme(ig) + max(0., (ztv(ig,k)-ztv(ig,k+1))/(zlev(ig, & |
---|
328 | k+1)-zlev(ig,k))) |
---|
329 | ! s *(zlev(ig,k+1)-zlev(ig,k)) |
---|
330 | END DO |
---|
331 | END DO |
---|
332 | DO ig = 1, ngrid |
---|
333 | IF (norme(ig)>1.E-10) THEN |
---|
334 | zalim(ig) = max(10.*zalim(ig)/norme(ig), zlev(ig,2)) |
---|
335 | ! zalim(ig)=min(zalim(ig),zlev(ig,lentr(ig))) |
---|
336 | END IF |
---|
337 | END DO |
---|
338 | ! détermination du lalim correspondant |
---|
339 | DO k = 1, klev - 1 |
---|
340 | DO ig = 1, ngrid |
---|
341 | IF ((zalim(ig)>zlev(ig,k)) .AND. (zalim(ig)<=zlev(ig,k+1))) THEN |
---|
342 | lalim(ig) = k |
---|
343 | END IF |
---|
344 | END DO |
---|
345 | END DO |
---|
346 | |
---|
347 | ! definition de l'entrainement des couches |
---|
348 | DO l = 1, klev - 1 |
---|
349 | DO ig = 1, ngrid |
---|
350 | IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<lentr(ig)) THEN |
---|
351 | entr_star(ig, l) = max((ztv(ig,l)-ztv(ig,l+1)), 0.) & ! s |
---|
352 | ! *(zlev(ig,l+1)-zlev(ig,l)) |
---|
353 | *sqrt(zlev(ig,l+1)) |
---|
354 | ! autre def |
---|
355 | ! entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1) |
---|
356 | ! s /zlev(ig,lentr(ig)+2)))**(3./2.) |
---|
357 | END IF |
---|
358 | END DO |
---|
359 | END DO |
---|
360 | ! nouveau test |
---|
361 | ! if (therm) then |
---|
362 | DO l = 1, klev - 1 |
---|
363 | DO ig = 1, ngrid |
---|
364 | IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<=lalim(ig) .AND. & |
---|
365 | zalim(ig)>1.E-10) THEN |
---|
366 | ! if (l.le.lentr(ig)) then |
---|
367 | ! entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1) |
---|
368 | ! s /zalim(ig)))**(3./2.) |
---|
369 | ! write(10,*)zlev(ig,l),entr_star(ig,l) |
---|
370 | END IF |
---|
371 | END DO |
---|
372 | END DO |
---|
373 | ! endif |
---|
374 | ! pas de thermique si couche 1 stable |
---|
375 | DO ig = 1, ngrid |
---|
376 | IF (lmin(ig)>5) THEN |
---|
377 | DO l = 1, klev |
---|
378 | entr_star(ig, l) = 0. |
---|
379 | END DO |
---|
380 | END IF |
---|
381 | END DO |
---|
382 | ! calcul de l entrainement total |
---|
383 | DO ig = 1, ngrid |
---|
384 | entr_star_tot(ig) = 0. |
---|
385 | END DO |
---|
386 | DO ig = 1, ngrid |
---|
387 | DO k = 1, klev |
---|
388 | entr_star_tot(ig) = entr_star_tot(ig) + entr_star(ig, k) |
---|
389 | END DO |
---|
390 | END DO |
---|
391 | ! Calcul entrainement normalise |
---|
392 | DO ig = 1, ngrid |
---|
393 | IF (entr_star_tot(ig)>1.E-10) THEN |
---|
394 | ! do l=1,lentr(ig) |
---|
395 | DO l = 1, klev |
---|
396 | ! def possibles pour entr_star: zdthetadz, dthetadz, zdtheta |
---|
397 | entr_star(ig, l) = entr_star(ig, l)/entr_star_tot(ig) |
---|
398 | END DO |
---|
399 | END IF |
---|
400 | END DO |
---|
401 | |
---|
402 | ! print*,'fin calcul entr_star' |
---|
403 | DO k = 1, klev |
---|
404 | DO ig = 1, ngrid |
---|
405 | ztva(ig, k) = ztv(ig, k) |
---|
406 | END DO |
---|
407 | END DO |
---|
408 | ! RC |
---|
409 | ! print*,'7 OK convect8' |
---|
410 | DO k = 1, klev + 1 |
---|
411 | DO ig = 1, ngrid |
---|
412 | zw2(ig, k) = 0. |
---|
413 | fmc(ig, k) = 0. |
---|
414 | ! CR |
---|
415 | f_star(ig, k) = 0. |
---|
416 | ! RC |
---|
417 | larg_cons(ig, k) = 0. |
---|
418 | larg_detr(ig, k) = 0. |
---|
419 | wa_moy(ig, k) = 0. |
---|
420 | END DO |
---|
421 | END DO |
---|
422 | |
---|
423 | ! print*,'8 OK convect8' |
---|
424 | DO ig = 1, ngrid |
---|
425 | linter(ig) = 1. |
---|
426 | lmaxa(ig) = 1 |
---|
427 | lmix(ig) = 1 |
---|
428 | wmaxa(ig) = 0. |
---|
429 | END DO |
---|
430 | |
---|
431 | ! CR: |
---|
432 | DO l = 1, nlay - 2 |
---|
433 | DO ig = 1, ngrid |
---|
434 | IF (ztv(ig,l)>ztv(ig,l+1) .AND. entr_star(ig,l)>1.E-10 .AND. & |
---|
435 | zw2(ig,l)<1E-10) THEN |
---|
436 | f_star(ig, l+1) = entr_star(ig, l) |
---|
437 | ! test:calcul de dteta |
---|
438 | zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* & |
---|
439 | (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l)) |
---|
440 | larg_detr(ig, l) = 0. |
---|
441 | ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+entr_star(ig, & |
---|
442 | l)>1.E-10)) THEN |
---|
443 | f_star(ig, l+1) = f_star(ig, l) + entr_star(ig, l) |
---|
444 | ztva(ig, l) = (f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)*ztv(ig,l))/ & |
---|
445 | f_star(ig, l+1) |
---|
446 | zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/f_star(ig,l+1))**2 + & |
---|
447 | 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l)) |
---|
448 | END IF |
---|
449 | ! determination de zmax continu par interpolation lineaire |
---|
450 | IF (zw2(ig,l+1)<0.) THEN |
---|
451 | ! test |
---|
452 | IF (abs(zw2(ig,l+1)-zw2(ig,l))<1E-10) THEN |
---|
453 | ! print*,'pb linter' |
---|
454 | END IF |
---|
455 | linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( & |
---|
456 | ig,l)) |
---|
457 | zw2(ig, l+1) = 0. |
---|
458 | lmaxa(ig) = l |
---|
459 | ELSE |
---|
460 | IF (zw2(ig,l+1)<0.) THEN |
---|
461 | ! print*,'pb1 zw2<0' |
---|
462 | END IF |
---|
463 | wa_moy(ig, l+1) = sqrt(zw2(ig,l+1)) |
---|
464 | END IF |
---|
465 | IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN |
---|
466 | ! lmix est le niveau de la couche ou w (wa_moy) est maximum |
---|
467 | lmix(ig) = l + 1 |
---|
468 | wmaxa(ig) = wa_moy(ig, l+1) |
---|
469 | END IF |
---|
470 | END DO |
---|
471 | END DO |
---|
472 | ! print*,'fin calcul zw2' |
---|
473 | |
---|
474 | ! Calcul de la couche correspondant a la hauteur du thermique |
---|
475 | DO ig = 1, ngrid |
---|
476 | lmax(ig) = lentr(ig) |
---|
477 | ! lmax(ig)=lalim(ig) |
---|
478 | END DO |
---|
479 | DO ig = 1, ngrid |
---|
480 | DO l = nlay, lentr(ig) + 1, -1 |
---|
481 | ! do l=nlay,lalim(ig)+1,-1 |
---|
482 | IF (zw2(ig,l)<=1.E-10) THEN |
---|
483 | lmax(ig) = l - 1 |
---|
484 | END IF |
---|
485 | END DO |
---|
486 | END DO |
---|
487 | ! pas de thermique si couche 1 stable |
---|
488 | DO ig = 1, ngrid |
---|
489 | IF (lmin(ig)>5) THEN |
---|
490 | lmax(ig) = 1 |
---|
491 | lmin(ig) = 1 |
---|
492 | lentr(ig) = 1 |
---|
493 | lalim(ig) = 1 |
---|
494 | END IF |
---|
495 | END DO |
---|
496 | |
---|
497 | ! Determination de zw2 max |
---|
498 | DO ig = 1, ngrid |
---|
499 | wmax(ig) = 0. |
---|
500 | END DO |
---|
501 | |
---|
502 | DO l = 1, nlay |
---|
503 | DO ig = 1, ngrid |
---|
504 | IF (l<=lmax(ig)) THEN |
---|
505 | IF (zw2(ig,l)<0.) THEN |
---|
506 | ! print*,'pb2 zw2<0' |
---|
507 | END IF |
---|
508 | zw2(ig, l) = sqrt(zw2(ig,l)) |
---|
509 | wmax(ig) = max(wmax(ig), zw2(ig,l)) |
---|
510 | ELSE |
---|
511 | zw2(ig, l) = 0. |
---|
512 | END IF |
---|
513 | END DO |
---|
514 | END DO |
---|
515 | |
---|
516 | ! Longueur caracteristique correspondant a la hauteur des thermiques. |
---|
517 | DO ig = 1, ngrid |
---|
518 | zmax(ig) = 0. |
---|
519 | zlevinter(ig) = zlev(ig, 1) |
---|
520 | END DO |
---|
521 | DO ig = 1, ngrid |
---|
522 | ! calcul de zlevinter |
---|
523 | zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + & |
---|
524 | zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig))) |
---|
525 | zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,lmin(ig))) |
---|
526 | END DO |
---|
527 | DO ig = 1, ngrid |
---|
528 | ! write(8,*)zmax(ig),lmax(ig),lentr(ig),lmin(ig) |
---|
529 | END DO |
---|
530 | ! on stope après les calculs de zmax et wmax |
---|
531 | RETURN |
---|
532 | |
---|
533 | ! print*,'avant fermeture' |
---|
534 | ! Fermeture,determination de f |
---|
535 | ! Attention! entrainement normalisé ou pas? |
---|
536 | DO ig = 1, ngrid |
---|
537 | entr_star2(ig) = 0. |
---|
538 | END DO |
---|
539 | DO ig = 1, ngrid |
---|
540 | IF (entr_star_tot(ig)<1.E-10) THEN |
---|
541 | f(ig) = 0. |
---|
542 | ELSE |
---|
543 | DO k = lmin(ig), lentr(ig) |
---|
544 | ! do k=lmin(ig),lalim(ig) |
---|
545 | entr_star2(ig) = entr_star2(ig) + entr_star(ig, k)**2/(rho(ig,k)*( & |
---|
546 | zlev(ig,k+1)-zlev(ig,k))) |
---|
547 | END DO |
---|
548 | ! Nouvelle fermeture |
---|
549 | f(ig) = wmax(ig)/(max(500.,zmax(ig))*r_aspect*entr_star2(ig)) |
---|
550 | ! s *entr_star_tot(ig) |
---|
551 | ! test |
---|
552 | ! if (first) then |
---|
553 | f(ig) = f(ig) + (f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)*wmax(ig)) |
---|
554 | ! endif |
---|
555 | END IF |
---|
556 | f0(ig) = f(ig) |
---|
557 | ! first=.true. |
---|
558 | END DO |
---|
559 | ! print*,'apres fermeture' |
---|
560 | ! on stoppe après la fermeture |
---|
561 | RETURN |
---|
562 | ! Calcul de l'entrainement |
---|
563 | DO k = 1, klev |
---|
564 | DO ig = 1, ngrid |
---|
565 | entr(ig, k) = f(ig)*entr_star(ig, k) |
---|
566 | END DO |
---|
567 | END DO |
---|
568 | ! on stoppe après le calcul de entr |
---|
569 | ! RETURN |
---|
570 | ! CR:test pour entrainer moins que la masse |
---|
571 | ! do ig=1,ngrid |
---|
572 | ! do l=1,lentr(ig) |
---|
573 | ! if ((entr(ig,l)*ptimestep).gt.(0.9*masse(ig,l))) then |
---|
574 | ! entr(ig,l+1)=entr(ig,l+1)+entr(ig,l) |
---|
575 | ! s -0.9*masse(ig,l)/ptimestep |
---|
576 | ! entr(ig,l)=0.9*masse(ig,l)/ptimestep |
---|
577 | ! endif |
---|
578 | ! enddo |
---|
579 | ! enddo |
---|
580 | ! CR: fin test |
---|
581 | ! Calcul des flux |
---|
582 | DO ig = 1, ngrid |
---|
583 | DO l = 1, lmax(ig) - 1 |
---|
584 | fmc(ig, l+1) = fmc(ig, l) + entr(ig, l) |
---|
585 | END DO |
---|
586 | END DO |
---|
587 | |
---|
588 | ! RC |
---|
589 | |
---|
590 | |
---|
591 | ! print*,'9 OK convect8' |
---|
592 | ! print*,'WA1 ',wa_moy |
---|
593 | |
---|
594 | ! determination de l'indice du debut de la mixed layer ou w decroit |
---|
595 | |
---|
596 | ! calcul de la largeur de chaque ascendance dans le cas conservatif. |
---|
597 | ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant |
---|
598 | ! d'une couche est égale à la hauteur de la couche alimentante. |
---|
599 | ! La vitesse maximale dans l'ascendance est aussi prise comme estimation |
---|
600 | ! de la vitesse d'entrainement horizontal dans la couche alimentante. |
---|
601 | |
---|
602 | DO l = 2, nlay |
---|
603 | DO ig = 1, ngrid |
---|
604 | IF (l<=lmaxa(ig)) THEN |
---|
605 | zw = max(wa_moy(ig,l), 1.E-10) |
---|
606 | larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw) |
---|
607 | END IF |
---|
608 | END DO |
---|
609 | END DO |
---|
610 | |
---|
611 | DO l = 2, nlay |
---|
612 | DO ig = 1, ngrid |
---|
613 | IF (l<=lmaxa(ig)) THEN |
---|
614 | ! if (idetr.eq.0) then |
---|
615 | ! cette option est finalement en dur. |
---|
616 | IF ((l_mix*zlev(ig,l))<0.) THEN |
---|
617 | ! print*,'pb l_mix*zlev<0' |
---|
618 | END IF |
---|
619 | ! CR: test: nouvelle def de lambda |
---|
620 | ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) |
---|
621 | IF (zw2(ig,l)>1.E-10) THEN |
---|
622 | larg_detr(ig, l) = sqrt((l_mix/zw2(ig,l))*zlev(ig,l)) |
---|
623 | ELSE |
---|
624 | larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l)) |
---|
625 | END IF |
---|
626 | ! RC |
---|
627 | ! else if (idetr.eq.1) then |
---|
628 | ! larg_detr(ig,l)=larg_cons(ig,l) |
---|
629 | ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig)) |
---|
630 | ! else if (idetr.eq.2) then |
---|
631 | ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) |
---|
632 | ! s *sqrt(wa_moy(ig,l)) |
---|
633 | ! else if (idetr.eq.4) then |
---|
634 | ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) |
---|
635 | ! s *wa_moy(ig,l) |
---|
636 | ! endif |
---|
637 | END IF |
---|
638 | END DO |
---|
639 | END DO |
---|
640 | |
---|
641 | ! print*,'10 OK convect8' |
---|
642 | ! print*,'WA2 ',wa_moy |
---|
643 | ! calcul de la fraction de la maille concernée par l'ascendance en tenant |
---|
644 | ! compte de l'epluchage du thermique. |
---|
645 | |
---|
646 | ! CR def de zmix continu (profil parabolique des vitesses) |
---|
647 | DO ig = 1, ngrid |
---|
648 | IF (lmix(ig)>1.) THEN |
---|
649 | ! test |
---|
650 | IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- & |
---|
651 | (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- & |
---|
652 | zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- & |
---|
653 | (zlev(ig,lmix(ig)))))>1E-10) THEN |
---|
654 | |
---|
655 | zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) & |
---|
656 | )**2-(zlev(ig,lmix(ig)+1))**2)-(zw2(ig,lmix(ig))-zw2(ig, & |
---|
657 | lmix(ig)+1))*((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))/ & |
---|
658 | (2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- & |
---|
659 | (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- & |
---|
660 | zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig)))))) |
---|
661 | ELSE |
---|
662 | zmix(ig) = zlev(ig, lmix(ig)) |
---|
663 | ! print*,'pb zmix' |
---|
664 | END IF |
---|
665 | ELSE |
---|
666 | zmix(ig) = 0. |
---|
667 | END IF |
---|
668 | ! test |
---|
669 | IF ((zmax(ig)-zmix(ig))<0.) THEN |
---|
670 | zmix(ig) = 0.99*zmax(ig) |
---|
671 | ! print*,'pb zmix>zmax' |
---|
672 | END IF |
---|
673 | END DO |
---|
674 | |
---|
675 | ! calcul du nouveau lmix correspondant |
---|
676 | DO ig = 1, ngrid |
---|
677 | DO l = 1, klev |
---|
678 | IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN |
---|
679 | lmix(ig) = l |
---|
680 | END IF |
---|
681 | END DO |
---|
682 | END DO |
---|
683 | |
---|
684 | DO l = 2, nlay |
---|
685 | DO ig = 1, ngrid |
---|
686 | IF (larg_cons(ig,l)>1.) THEN |
---|
687 | ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK' |
---|
688 | fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig)) |
---|
689 | ! test |
---|
690 | fraca(ig, l) = max(fraca(ig,l), 0.) |
---|
691 | fraca(ig, l) = min(fraca(ig,l), 0.5) |
---|
692 | fracd(ig, l) = 1. - fraca(ig, l) |
---|
693 | fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig)) |
---|
694 | ELSE |
---|
695 | ! wa_moy(ig,l)=0. |
---|
696 | fraca(ig, l) = 0. |
---|
697 | fracc(ig, l) = 0. |
---|
698 | fracd(ig, l) = 1. |
---|
699 | END IF |
---|
700 | END DO |
---|
701 | END DO |
---|
702 | ! CR: calcul de fracazmix |
---|
703 | DO ig = 1, ngrid |
---|
704 | fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ & |
---|
705 | (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + & |
---|
706 | fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca(ig & |
---|
707 | ,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig))) |
---|
708 | END DO |
---|
709 | |
---|
710 | DO l = 2, nlay |
---|
711 | DO ig = 1, ngrid |
---|
712 | IF (larg_cons(ig,l)>1.) THEN |
---|
713 | IF (l>lmix(ig)) THEN |
---|
714 | ! test |
---|
715 | IF (zmax(ig)-zmix(ig)<1.E-10) THEN |
---|
716 | ! print*,'pb xxx' |
---|
717 | xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig)) |
---|
718 | ELSE |
---|
719 | xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig)) |
---|
720 | END IF |
---|
721 | IF (idetr==0) THEN |
---|
722 | fraca(ig, l) = fracazmix(ig) |
---|
723 | ELSE IF (idetr==1) THEN |
---|
724 | fraca(ig, l) = fracazmix(ig)*xxx(ig, l) |
---|
725 | ELSE IF (idetr==2) THEN |
---|
726 | fraca(ig, l) = fracazmix(ig)*(1.-(1.-xxx(ig,l))**2) |
---|
727 | ELSE |
---|
728 | fraca(ig, l) = fracazmix(ig)*xxx(ig, l)**2 |
---|
729 | END IF |
---|
730 | ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL' |
---|
731 | fraca(ig, l) = max(fraca(ig,l), 0.) |
---|
732 | fraca(ig, l) = min(fraca(ig,l), 0.5) |
---|
733 | fracd(ig, l) = 1. - fraca(ig, l) |
---|
734 | fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig)) |
---|
735 | END IF |
---|
736 | END IF |
---|
737 | END DO |
---|
738 | END DO |
---|
739 | |
---|
740 | ! print*,'fin calcul fraca' |
---|
741 | ! print*,'11 OK convect8' |
---|
742 | ! print*,'Ea3 ',wa_moy |
---|
743 | ! ------------------------------------------------------------------ |
---|
744 | ! Calcul de fracd, wd |
---|
745 | ! somme wa - wd = 0 |
---|
746 | ! ------------------------------------------------------------------ |
---|
747 | |
---|
748 | |
---|
749 | DO ig = 1, ngrid |
---|
750 | fm(ig, 1) = 0. |
---|
751 | fm(ig, nlay+1) = 0. |
---|
752 | END DO |
---|
753 | |
---|
754 | DO l = 2, nlay |
---|
755 | DO ig = 1, ngrid |
---|
756 | fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l) |
---|
757 | ! CR:test |
---|
758 | IF (entr(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) THEN |
---|
759 | fm(ig, l) = fm(ig, l-1) |
---|
760 | ! write(1,*)'ajustement fm, l',l |
---|
761 | END IF |
---|
762 | ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l) |
---|
763 | ! RC |
---|
764 | END DO |
---|
765 | DO ig = 1, ngrid |
---|
766 | IF (fracd(ig,l)<0.1) THEN |
---|
767 | abort_message = 'fracd trop petit' |
---|
768 | CALL abort_gcm(modname, abort_message, 1) |
---|
769 | |
---|
770 | ELSE |
---|
771 | ! vitesse descendante "diagnostique" |
---|
772 | wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l)) |
---|
773 | END IF |
---|
774 | END DO |
---|
775 | END DO |
---|
776 | |
---|
777 | DO l = 1, nlay |
---|
778 | DO ig = 1, ngrid |
---|
779 | ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) |
---|
780 | masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg |
---|
781 | END DO |
---|
782 | END DO |
---|
783 | |
---|
784 | ! print*,'12 OK convect8' |
---|
785 | ! print*,'WA4 ',wa_moy |
---|
786 | ! c------------------------------------------------------------------ |
---|
787 | ! calcul du transport vertical |
---|
788 | ! ------------------------------------------------------------------ |
---|
789 | |
---|
790 | GO TO 4444 |
---|
791 | ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep |
---|
792 | DO l = 2, nlay - 1 |
---|
793 | DO ig = 1, ngrid |
---|
794 | IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( & |
---|
795 | ig,l+1)) THEN |
---|
796 | ! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM=' |
---|
797 | ! s ,fm(ig,l+1)*ptimestep |
---|
798 | ! s ,' M=',masse(ig,l),masse(ig,l+1) |
---|
799 | END IF |
---|
800 | END DO |
---|
801 | END DO |
---|
802 | |
---|
803 | DO l = 1, nlay |
---|
804 | DO ig = 1, ngrid |
---|
805 | IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN |
---|
806 | ! print*,'WARN!!! E>M ig=',ig,' l=',l,' E==' |
---|
807 | ! s ,entr(ig,l)*ptimestep |
---|
808 | ! s ,' M=',masse(ig,l) |
---|
809 | END IF |
---|
810 | END DO |
---|
811 | END DO |
---|
812 | |
---|
813 | DO l = 1, nlay |
---|
814 | DO ig = 1, ngrid |
---|
815 | IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN |
---|
816 | ! print*,'WARN!!! fm exagere ig=',ig,' l=',l |
---|
817 | ! s ,' FM=',fm(ig,l) |
---|
818 | END IF |
---|
819 | IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN |
---|
820 | ! print*,'WARN!!! masse exagere ig=',ig,' l=',l |
---|
821 | ! s ,' M=',masse(ig,l) |
---|
822 | ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)', |
---|
823 | ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l) |
---|
824 | ! print*,'zlev(ig,l+1),zlev(ig,l)' |
---|
825 | ! s ,zlev(ig,l+1),zlev(ig,l) |
---|
826 | ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)' |
---|
827 | ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1) |
---|
828 | END IF |
---|
829 | IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN |
---|
830 | ! print*,'WARN!!! entr exagere ig=',ig,' l=',l |
---|
831 | ! s ,' E=',entr(ig,l) |
---|
832 | END IF |
---|
833 | END DO |
---|
834 | END DO |
---|
835 | |
---|
836 | 4444 CONTINUE |
---|
837 | |
---|
838 | ! CR:redefinition du entr |
---|
839 | DO l = 1, nlay |
---|
840 | DO ig = 1, ngrid |
---|
841 | detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1) |
---|
842 | IF (detr(ig,l)<0.) THEN |
---|
843 | ! entr(ig,l)=entr(ig,l)-detr(ig,l) |
---|
844 | fm(ig, l+1) = fm(ig, l) + entr(ig, l) |
---|
845 | detr(ig, l) = 0. |
---|
846 | ! print*,'WARNING !!! detrainement negatif ',ig,l |
---|
847 | END IF |
---|
848 | END DO |
---|
849 | END DO |
---|
850 | ! RC |
---|
851 | IF (w2di==1) THEN |
---|
852 | fm0 = fm0 + ptimestep*(fm-fm0)/tho |
---|
853 | entr0 = entr0 + ptimestep*(entr-entr0)/tho |
---|
854 | ELSE |
---|
855 | fm0 = fm |
---|
856 | entr0 = entr |
---|
857 | END IF |
---|
858 | |
---|
859 | IF (1==1) THEN |
---|
860 | CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, & |
---|
861 | zha) |
---|
862 | CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, & |
---|
863 | zoa) |
---|
864 | ELSE |
---|
865 | CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, & |
---|
866 | zdhadj, zha) |
---|
867 | CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, & |
---|
868 | pdoadj, zoa) |
---|
869 | END IF |
---|
870 | |
---|
871 | IF (1==0) THEN |
---|
872 | CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, & |
---|
873 | zu, zv, pduadj, pdvadj, zua, zva) |
---|
874 | ELSE |
---|
875 | CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, & |
---|
876 | zua) |
---|
877 | CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, & |
---|
878 | zva) |
---|
879 | END IF |
---|
880 | |
---|
881 | DO l = 1, nlay |
---|
882 | DO ig = 1, ngrid |
---|
883 | zf = 0.5*(fracc(ig,l)+fracc(ig,l+1)) |
---|
884 | zf2 = zf/(1.-zf) |
---|
885 | thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l))**2 |
---|
886 | wth2(ig, l) = zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2 |
---|
887 | END DO |
---|
888 | END DO |
---|
889 | |
---|
890 | |
---|
891 | |
---|
892 | ! print*,'13 OK convect8' |
---|
893 | ! print*,'WA5 ',wa_moy |
---|
894 | DO l = 1, nlay |
---|
895 | DO ig = 1, ngrid |
---|
896 | pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l) |
---|
897 | END DO |
---|
898 | END DO |
---|
899 | |
---|
900 | |
---|
901 | ! do l=1,nlay |
---|
902 | ! do ig=1,ngrid |
---|
903 | ! if(abs(pdtadj(ig,l))*86400..gt.500.) then |
---|
904 | ! print*,'WARN!!! ig=',ig,' l=',l |
---|
905 | ! s ,' pdtadj=',pdtadj(ig,l) |
---|
906 | ! endif |
---|
907 | ! if(abs(pdoadj(ig,l))*86400..gt.1.) then |
---|
908 | ! print*,'WARN!!! ig=',ig,' l=',l |
---|
909 | ! s ,' pdoadj=',pdoadj(ig,l) |
---|
910 | ! endif |
---|
911 | ! enddo |
---|
912 | ! enddo |
---|
913 | |
---|
914 | ! print*,'14 OK convect8' |
---|
915 | ! ------------------------------------------------------------------ |
---|
916 | ! Calculs pour les sorties |
---|
917 | ! ------------------------------------------------------------------ |
---|
918 | |
---|
919 | IF (sorties) THEN |
---|
920 | DO l = 1, nlay |
---|
921 | DO ig = 1, ngrid |
---|
922 | zla(ig, l) = (1.-fracd(ig,l))*zmax(ig) |
---|
923 | zld(ig, l) = fracd(ig, l)*zmax(ig) |
---|
924 | IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ & |
---|
925 | (1.-fracd(ig,l)) |
---|
926 | END DO |
---|
927 | END DO |
---|
928 | |
---|
929 | ! deja fait |
---|
930 | ! do l=1,nlay |
---|
931 | ! do ig=1,ngrid |
---|
932 | ! detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1) |
---|
933 | ! if (detr(ig,l).lt.0.) then |
---|
934 | ! entr(ig,l)=entr(ig,l)-detr(ig,l) |
---|
935 | ! detr(ig,l)=0. |
---|
936 | ! print*,'WARNING !!! detrainement negatif ',ig,l |
---|
937 | ! endif |
---|
938 | ! enddo |
---|
939 | ! enddo |
---|
940 | |
---|
941 | ! print*,'15 OK convect8' |
---|
942 | |
---|
943 | isplit = isplit + 1 |
---|
944 | |
---|
945 | |
---|
946 | ! #define und |
---|
947 | GO TO 123 |
---|
948 | #ifdef und |
---|
949 | CALL writeg1d(1, nlay, wd, 'wd ', 'wd ') |
---|
950 | CALL writeg1d(1, nlay, zwa, 'wa ', 'wa ') |
---|
951 | CALL writeg1d(1, nlay, fracd, 'fracd ', 'fracd ') |
---|
952 | CALL writeg1d(1, nlay, fraca, 'fraca ', 'fraca ') |
---|
953 | CALL writeg1d(1, nlay, wa_moy, 'wam ', 'wam ') |
---|
954 | CALL writeg1d(1, nlay, zla, 'la ', 'la ') |
---|
955 | CALL writeg1d(1, nlay, zld, 'ld ', 'ld ') |
---|
956 | CALL writeg1d(1, nlay, pt, 'pt ', 'pt ') |
---|
957 | CALL writeg1d(1, nlay, zh, 'zh ', 'zh ') |
---|
958 | CALL writeg1d(1, nlay, zha, 'zha ', 'zha ') |
---|
959 | CALL writeg1d(1, nlay, zu, 'zu ', 'zu ') |
---|
960 | CALL writeg1d(1, nlay, zv, 'zv ', 'zv ') |
---|
961 | CALL writeg1d(1, nlay, zo, 'zo ', 'zo ') |
---|
962 | CALL writeg1d(1, nlay, wh, 'wh ', 'wh ') |
---|
963 | CALL writeg1d(1, nlay, wu, 'wu ', 'wu ') |
---|
964 | CALL writeg1d(1, nlay, wv, 'wv ', 'wv ') |
---|
965 | CALL writeg1d(1, nlay, wo, 'w15uo ', 'wXo ') |
---|
966 | CALL writeg1d(1, nlay, zdhadj, 'zdhadj ', 'zdhadj ') |
---|
967 | CALL writeg1d(1, nlay, pduadj, 'pduadj ', 'pduadj ') |
---|
968 | CALL writeg1d(1, nlay, pdvadj, 'pdvadj ', 'pdvadj ') |
---|
969 | CALL writeg1d(1, nlay, pdoadj, 'pdoadj ', 'pdoadj ') |
---|
970 | CALL writeg1d(1, nlay, entr, 'entr ', 'entr ') |
---|
971 | CALL writeg1d(1, nlay, detr, 'detr ', 'detr ') |
---|
972 | CALL writeg1d(1, nlay, fm, 'fm ', 'fm ') |
---|
973 | |
---|
974 | CALL writeg1d(1, nlay, pdtadj, 'pdtadj ', 'pdtadj ') |
---|
975 | CALL writeg1d(1, nlay, pplay, 'pplay ', 'pplay ') |
---|
976 | CALL writeg1d(1, nlay, pplev, 'pplev ', 'pplev ') |
---|
977 | |
---|
978 | ! recalcul des flux en diagnostique... |
---|
979 | ! print*,'PAS DE TEMPS ',ptimestep |
---|
980 | CALL dt2f(pplev, pplay, pt, pdtadj, wh) |
---|
981 | CALL writeg1d(1, nlay, wh, 'wh2 ', 'wh2 ') |
---|
982 | #endif |
---|
983 | 123 CONTINUE |
---|
984 | |
---|
985 | END IF |
---|
986 | |
---|
987 | ! if(wa_moy(1,4).gt.1.e-10) stop |
---|
988 | |
---|
989 | ! print*,'19 OK convect8' |
---|
990 | RETURN |
---|
991 | END SUBROUTINE calcul_sec |
---|
992 | |
---|
993 | SUBROUTINE fermeture_seche(ngrid, nlay, pplay, pplev, pphi, zlev, rhobarz, & |
---|
994 | f0, zpspsk, alim_star, zh, zo, lentr, lmin, nu_min, nu_max, r_aspect, & |
---|
995 | zmax, wmax) |
---|
996 | |
---|
997 | USE dimphy |
---|
998 | IMPLICIT NONE |
---|
999 | |
---|
1000 | include "dimensions.h" |
---|
1001 | ! ccc#include "dimphy.h" |
---|
1002 | include "YOMCST.h" |
---|
1003 | |
---|
1004 | INTEGER ngrid, nlay |
---|
1005 | REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1) |
---|
1006 | REAL pphi(ngrid, nlay) |
---|
1007 | REAL zlev(klon, klev+1) |
---|
1008 | REAL alim_star(klon, klev) |
---|
1009 | REAL f0(klon) |
---|
1010 | INTEGER lentr(klon) |
---|
1011 | INTEGER lmin(klon) |
---|
1012 | REAL zmax(klon) |
---|
1013 | REAL wmax(klon) |
---|
1014 | REAL nu_min |
---|
1015 | REAL nu_max |
---|
1016 | REAL r_aspect |
---|
1017 | REAL rhobarz(klon, klev+1) |
---|
1018 | REAL zh(klon, klev) |
---|
1019 | REAL zo(klon, klev) |
---|
1020 | REAL zpspsk(klon, klev) |
---|
1021 | |
---|
1022 | INTEGER ig, l |
---|
1023 | |
---|
1024 | REAL f_star(klon, klev+1) |
---|
1025 | REAL detr_star(klon, klev) |
---|
1026 | REAL entr_star(klon, klev) |
---|
1027 | REAL zw2(klon, klev+1) |
---|
1028 | REAL linter(klon) |
---|
1029 | INTEGER lmix(klon) |
---|
1030 | INTEGER lmax(klon) |
---|
1031 | REAL zlevinter(klon) |
---|
1032 | REAL wa_moy(klon, klev+1) |
---|
1033 | REAL wmaxa(klon) |
---|
1034 | REAL ztv(klon, klev) |
---|
1035 | REAL ztva(klon, klev) |
---|
1036 | REAL nu(klon, klev) |
---|
1037 | ! real zmax0_sec(klon) |
---|
1038 | ! save zmax0_sec |
---|
1039 | REAL, SAVE, ALLOCATABLE :: zmax0_sec(:) |
---|
1040 | !$OMP THREADPRIVATE(zmax0_sec) |
---|
1041 | LOGICAL, SAVE :: first = .TRUE. |
---|
1042 | !$OMP THREADPRIVATE(first) |
---|
1043 | |
---|
1044 | IF (first) THEN |
---|
1045 | ALLOCATE (zmax0_sec(klon)) |
---|
1046 | first = .FALSE. |
---|
1047 | END IF |
---|
1048 | |
---|
1049 | DO l = 1, nlay |
---|
1050 | DO ig = 1, ngrid |
---|
1051 | ztv(ig, l) = zh(ig, l)/zpspsk(ig, l) |
---|
1052 | ztv(ig, l) = ztv(ig, l)*(1.+retv*zo(ig,l)) |
---|
1053 | END DO |
---|
1054 | END DO |
---|
1055 | DO l = 1, nlay - 2 |
---|
1056 | DO ig = 1, ngrid |
---|
1057 | IF (ztv(ig,l)>ztv(ig,l+1) .AND. alim_star(ig,l)>1.E-10 .AND. & |
---|
1058 | zw2(ig,l)<1E-10) THEN |
---|
1059 | f_star(ig, l+1) = alim_star(ig, l) |
---|
1060 | ! test:calcul de dteta |
---|
1061 | zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* & |
---|
1062 | (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l)) |
---|
1063 | ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+alim_star(ig, & |
---|
1064 | l))>1.E-10) THEN |
---|
1065 | ! estimation du detrainement a partir de la geometrie du pas |
---|
1066 | ! precedent |
---|
1067 | ! tests sur la definition du detr |
---|
1068 | nu(ig, l) = (nu_min+nu_max)/2.*(1.-(nu_max-nu_min)/(nu_max+nu_min)* & |
---|
1069 | tanh((((ztva(ig,l-1)-ztv(ig,l))/ztv(ig,l))/0.0005))) |
---|
1070 | |
---|
1071 | detr_star(ig, l) = rhobarz(ig, l)*sqrt(zw2(ig,l))/ & |
---|
1072 | (r_aspect*zmax0_sec(ig))* & ! s |
---|
1073 | ! /(r_aspect*zmax0(ig))* |
---|
1074 | (sqrt(nu(ig,l)*zlev(ig,l+1)/sqrt(zw2(ig,l)))-sqrt(nu(ig,l)*zlev(ig, & |
---|
1075 | l)/sqrt(zw2(ig,l)))) |
---|
1076 | detr_star(ig, l) = detr_star(ig, l)/f0(ig) |
---|
1077 | IF ((detr_star(ig,l))>f_star(ig,l)) THEN |
---|
1078 | detr_star(ig, l) = f_star(ig, l) |
---|
1079 | END IF |
---|
1080 | entr_star(ig, l) = 0.9*detr_star(ig, l) |
---|
1081 | IF ((l<lentr(ig))) THEN |
---|
1082 | entr_star(ig, l) = 0. |
---|
1083 | ! detr_star(ig,l)=0. |
---|
1084 | END IF |
---|
1085 | ! print*,'ok detr_star' |
---|
1086 | ! prise en compte du detrainement dans le calcul du flux |
---|
1087 | f_star(ig, l+1) = f_star(ig, l) + alim_star(ig, l) + & |
---|
1088 | entr_star(ig, l) - detr_star(ig, l) |
---|
1089 | ! test sur le signe de f_star |
---|
1090 | IF ((f_star(ig,l+1)+detr_star(ig,l))>1.E-10) THEN |
---|
1091 | ! AM on melange Tl et qt du thermique |
---|
1092 | ztva(ig, l) = (f_star(ig,l)*ztva(ig,l-1)+(entr_star(ig, & |
---|
1093 | l)+alim_star(ig,l))*ztv(ig,l))/(f_star(ig,l+1)+detr_star(ig,l)) |
---|
1094 | zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/(f_star(ig, & |
---|
1095 | l+1)+detr_star(ig,l)))**2 + 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, & |
---|
1096 | l)*(zlev(ig,l+1)-zlev(ig,l)) |
---|
1097 | END IF |
---|
1098 | END IF |
---|
1099 | |
---|
1100 | IF (zw2(ig,l+1)<0.) THEN |
---|
1101 | linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( & |
---|
1102 | ig,l)) |
---|
1103 | zw2(ig, l+1) = 0. |
---|
1104 | ! print*,'linter=',linter(ig) |
---|
1105 | ELSE |
---|
1106 | wa_moy(ig, l+1) = sqrt(zw2(ig,l+1)) |
---|
1107 | END IF |
---|
1108 | IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN |
---|
1109 | ! lmix est le niveau de la couche ou w (wa_moy) est maximum |
---|
1110 | lmix(ig) = l + 1 |
---|
1111 | wmaxa(ig) = wa_moy(ig, l+1) |
---|
1112 | END IF |
---|
1113 | END DO |
---|
1114 | END DO |
---|
1115 | ! print*,'fin calcul zw2' |
---|
1116 | |
---|
1117 | ! Calcul de la couche correspondant a la hauteur du thermique |
---|
1118 | DO ig = 1, ngrid |
---|
1119 | lmax(ig) = lentr(ig) |
---|
1120 | END DO |
---|
1121 | DO ig = 1, ngrid |
---|
1122 | DO l = nlay, lentr(ig) + 1, -1 |
---|
1123 | IF (zw2(ig,l)<=1.E-10) THEN |
---|
1124 | lmax(ig) = l - 1 |
---|
1125 | END IF |
---|
1126 | END DO |
---|
1127 | END DO |
---|
1128 | ! pas de thermique si couche 1 stable |
---|
1129 | DO ig = 1, ngrid |
---|
1130 | IF (lmin(ig)>1) THEN |
---|
1131 | lmax(ig) = 1 |
---|
1132 | lmin(ig) = 1 |
---|
1133 | lentr(ig) = 1 |
---|
1134 | END IF |
---|
1135 | END DO |
---|
1136 | |
---|
1137 | ! Determination de zw2 max |
---|
1138 | DO ig = 1, ngrid |
---|
1139 | wmax(ig) = 0. |
---|
1140 | END DO |
---|
1141 | |
---|
1142 | DO l = 1, nlay |
---|
1143 | DO ig = 1, ngrid |
---|
1144 | IF (l<=lmax(ig)) THEN |
---|
1145 | IF (zw2(ig,l)<0.) THEN |
---|
1146 | ! print*,'pb2 zw2<0' |
---|
1147 | END IF |
---|
1148 | zw2(ig, l) = sqrt(zw2(ig,l)) |
---|
1149 | wmax(ig) = max(wmax(ig), zw2(ig,l)) |
---|
1150 | ELSE |
---|
1151 | zw2(ig, l) = 0. |
---|
1152 | END IF |
---|
1153 | END DO |
---|
1154 | END DO |
---|
1155 | |
---|
1156 | ! Longueur caracteristique correspondant a la hauteur des thermiques. |
---|
1157 | DO ig = 1, ngrid |
---|
1158 | zmax(ig) = 0. |
---|
1159 | zlevinter(ig) = zlev(ig, 1) |
---|
1160 | END DO |
---|
1161 | DO ig = 1, ngrid |
---|
1162 | ! calcul de zlevinter |
---|
1163 | zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + & |
---|
1164 | zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig))) |
---|
1165 | ! pour le cas ou on prend tjs lmin=1 |
---|
1166 | ! zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig))) |
---|
1167 | zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,1)) |
---|
1168 | zmax0_sec(ig) = zmax(ig) |
---|
1169 | END DO |
---|
1170 | |
---|
1171 | RETURN |
---|
1172 | END SUBROUTINE fermeture_seche |
---|