GCC Code Coverage Report


Directory: ./
File: phys/thermcell_old.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 1771 0.0%
Branches: 0 1518 0.0%

Line Branch Exec Source
1 SUBROUTINE thermcell_2002(ngrid, nlay, ptimestep, iflag_thermals, pplay, &
2 pplev, pphi, pu, pv, pt, po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0, &
3 fraca, wa_moy, r_aspect, l_mix, w2di, tho)
4
5 USE dimphy
6 USE write_field_phy
7 IMPLICIT NONE
8
9 ! =======================================================================
10
11 ! Calcul du transport verticale dans la couche limite en presence
12 ! de "thermiques" explicitement representes
13
14 ! R��criture � partir d'un listing papier � Habas, le 14/02/00
15
16 ! le thermique est suppos� homog�ne et dissip� par m�lange avec
17 ! son environnement. la longueur l_mix contr�le l'efficacit� du
18 ! m�lange
19
20 ! Le calcul du transport des diff�rentes esp�ces se fait en prenant
21 ! en compte:
22 ! 1. un flux de masse montant
23 ! 2. un flux de masse descendant
24 ! 3. un entrainement
25 ! 4. un detrainement
26
27 ! =======================================================================
28
29 ! -----------------------------------------------------------------------
30 ! declarations:
31 ! -------------
32
33 include "YOMCST.h"
34
35 ! arguments:
36 ! ----------
37
38 INTEGER ngrid, nlay, w2di, iflag_thermals
39 REAL tho
40 REAL ptimestep, l_mix, r_aspect
41 REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
42 REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
43 REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
44 REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
45 REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
46 REAL pphi(ngrid, nlay)
47 REAL fraca(ngrid, nlay+1), zw2(ngrid, nlay+1)
48
49 INTEGER, SAVE :: idetr = 3, lev_out = 1
50 !$OMP THREADPRIVATE(idetr,lev_out)
51
52 ! local:
53 ! ------
54
55 INTEGER, SAVE :: dvdq = 0, flagdq = 0, dqimpl = 1
56 LOGICAL, SAVE :: debut = .TRUE.
57 !$OMP THREADPRIVATE(dvdq,flagdq,debut,dqimpl)
58
59 INTEGER ig, k, l, lmax(klon, klev+1), lmaxa(klon), lmix(klon)
60 REAL zmax(klon), zw, zz, ztva(klon, klev), zzz
61
62 REAL zlev(klon, klev+1), zlay(klon, klev)
63 REAL zh(klon, klev), zdhadj(klon, klev)
64 REAL ztv(klon, klev)
65 REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
66 REAL wh(klon, klev+1)
67 REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
68 REAL zla(klon, klev+1)
69 REAL zwa(klon, klev+1)
70 REAL zld(klon, klev+1)
71 REAL zwd(klon, klev+1)
72 REAL zsortie(klon, klev)
73 REAL zva(klon, klev)
74 REAL zua(klon, klev)
75 REAL zoa(klon, klev)
76
77 REAL zha(klon, klev)
78 REAL wa_moy(klon, klev+1)
79 REAL fracc(klon, klev+1)
80 REAL zf, zf2
81 REAL thetath2(klon, klev), wth2(klon, klev)
82 ! common/comtherm/thetath2,wth2
83
84 REAL count_time
85
86 LOGICAL sorties
87 REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
88 REAL zpspsk(klon, klev)
89
90 REAL wmax(klon, klev), wmaxa(klon)
91
92 REAL wa(klon, klev, klev+1)
93 REAL wd(klon, klev+1)
94 REAL larg_part(klon, klev, klev+1)
95 REAL fracd(klon, klev+1)
96 REAL xxx(klon, klev+1)
97 REAL larg_cons(klon, klev+1)
98 REAL larg_detr(klon, klev+1)
99 REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
100 REAL pu_therm(klon, klev), pv_therm(klon, klev)
101 REAL fm(klon, klev+1), entr(klon, klev)
102 REAL fmc(klon, klev+1)
103
104 CHARACTER (LEN=2) :: str2
105 CHARACTER (LEN=10) :: str10
106
107 CHARACTER (LEN=20) :: modname = 'thermcell2002'
108 CHARACTER (LEN=80) :: abort_message
109
110 LOGICAL vtest(klon), down
111
112 EXTERNAL scopy
113
114 INTEGER ncorrec, ll
115 SAVE ncorrec
116 DATA ncorrec/0/
117 !$OMP THREADPRIVATE(ncorrec)
118
119
120 ! -----------------------------------------------------------------------
121 ! initialisation:
122 ! ---------------
123
124 sorties = .TRUE.
125 IF (ngrid/=klon) THEN
126 PRINT *
127 PRINT *, 'STOP dans convadj'
128 PRINT *, 'ngrid =', ngrid
129 PRINT *, 'klon =', klon
130 END IF
131
132 ! -----------------------------------------------------------------------
133 ! incrementation eventuelle de tendances precedentes:
134 ! ---------------------------------------------------
135
136 ! print*,'0 OK convect8'
137
138 DO l = 1, nlay
139 DO ig = 1, ngrid
140 zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa
141 zh(ig, l) = pt(ig, l)/zpspsk(ig, l)
142 zu(ig, l) = pu(ig, l)
143 zv(ig, l) = pv(ig, l)
144 zo(ig, l) = po(ig, l)
145 ztv(ig, l) = zh(ig, l)*(1.+0.61*zo(ig,l))
146 END DO
147 END DO
148
149 ! print*,'1 OK convect8'
150 ! --------------------
151
152
153 ! + + + + + + + + + + +
154
155
156 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
157 ! wh,wt,wo ...
158
159 ! + + + + + + + + + + + zh,zu,zv,zo,rho
160
161
162 ! -------------------- zlev(1)
163 ! \\\\\\\\\\\\\\\\\\\
164
165
166 ! -----------------------------------------------------------------------
167 ! Calcul des altitudes des couches
168 ! -----------------------------------------------------------------------
169
170 IF (debut) THEN
171 flagdq = (iflag_thermals-1000)/100
172 dvdq = (iflag_thermals-(1000+flagdq*100))/10
173 IF (flagdq==2) dqimpl = -1
174 IF (flagdq==3) dqimpl = 1
175 debut = .FALSE.
176 END IF
177 PRINT *, 'TH flag th ', iflag_thermals, flagdq, dvdq, dqimpl
178
179 DO l = 2, nlay
180 DO ig = 1, ngrid
181 zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
182 END DO
183 END DO
184 DO ig = 1, ngrid
185 zlev(ig, 1) = 0.
186 zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
187 END DO
188 DO l = 1, nlay
189 DO ig = 1, ngrid
190 zlay(ig, l) = pphi(ig, l)/rg
191 END DO
192 END DO
193
194 ! print*,'2 OK convect8'
195 ! -----------------------------------------------------------------------
196 ! Calcul des densites
197 ! -----------------------------------------------------------------------
198
199 DO l = 1, nlay
200 DO ig = 1, ngrid
201 rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*zh(ig,l))
202 END DO
203 END DO
204
205 DO l = 2, nlay
206 DO ig = 1, ngrid
207 rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
208 END DO
209 END DO
210
211 DO k = 1, nlay
212 DO l = 1, nlay + 1
213 DO ig = 1, ngrid
214 wa(ig, k, l) = 0.
215 END DO
216 END DO
217 END DO
218
219 ! print*,'3 OK convect8'
220 ! ------------------------------------------------------------------
221 ! Calcul de w2, quarre de w a partir de la cape
222 ! a partir de w2, on calcule wa, vitesse de l'ascendance
223
224 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
225 ! w2 est stoke dans wa
226
227 ! ATTENTION: dans convect8, on n'utilise le calcule des wa
228 ! independants par couches que pour calculer l'entrainement
229 ! a la base et la hauteur max de l'ascendance.
230
231 ! Indicages:
232 ! l'ascendance provenant du niveau k traverse l'interface l avec
233 ! une vitesse wa(k,l).
234
235 ! --------------------
236
237 ! + + + + + + + + + +
238
239 ! wa(k,l) ---- -------------------- l
240 ! / ! /||\ + + + + + + + + + +
241 ! ||
242 ! || --------------------
243 ! ||
244 ! || + + + + + + + + + +
245 ! ||
246 ! || --------------------
247 ! ||__
248 ! |___ + + + + + + + + + + k
249
250 ! --------------------
251
252
253
254 ! ------------------------------------------------------------------
255
256
257 DO k = 1, nlay - 1
258 DO ig = 1, ngrid
259 wa(ig, k, k) = 0.
260 wa(ig, k, k+1) = 2.*rg*(ztv(ig,k)-ztv(ig,k+1))/ztv(ig, k+1)* &
261 (zlev(ig,k+1)-zlev(ig,k))
262 END DO
263 DO l = k + 1, nlay - 1
264 DO ig = 1, ngrid
265 wa(ig, k, l+1) = wa(ig, k, l) + 2.*rg*(ztv(ig,k)-ztv(ig,l))/ztv(ig, l &
266 )*(zlev(ig,l+1)-zlev(ig,l))
267 END DO
268 END DO
269 DO ig = 1, ngrid
270 wa(ig, k, nlay+1) = 0.
271 END DO
272 END DO
273
274 ! print*,'4 OK convect8'
275 ! Calcul de la couche correspondant a la hauteur du thermique
276 DO k = 1, nlay - 1
277 DO ig = 1, ngrid
278 lmax(ig, k) = k
279 END DO
280 DO l = nlay, k + 1, -1
281 DO ig = 1, ngrid
282 IF (wa(ig,k,l)<=1.E-10) lmax(ig, k) = l - 1
283 END DO
284 END DO
285 END DO
286
287 ! print*,'5 OK convect8'
288 ! Calcule du w max du thermique
289 DO k = 1, nlay
290 DO ig = 1, ngrid
291 wmax(ig, k) = 0.
292 END DO
293 END DO
294
295 DO k = 1, nlay - 1
296 DO l = k, nlay
297 DO ig = 1, ngrid
298 IF (l<=lmax(ig,k)) THEN
299 wa(ig, k, l) = sqrt(wa(ig,k,l))
300 wmax(ig, k) = max(wmax(ig,k), wa(ig,k,l))
301 ELSE
302 wa(ig, k, l) = 0.
303 END IF
304 END DO
305 END DO
306 END DO
307
308 DO k = 1, nlay - 1
309 DO ig = 1, ngrid
310 pu_therm(ig, k) = sqrt(wmax(ig,k))
311 pv_therm(ig, k) = sqrt(wmax(ig,k))
312 END DO
313 END DO
314
315 ! print*,'6 OK convect8'
316 ! Longueur caracteristique correspondant a la hauteur des thermiques.
317 DO ig = 1, ngrid
318 zmax(ig) = 500.
319 END DO
320 ! print*,'LMAX LMAX LMAX '
321 DO k = 1, nlay - 1
322 DO ig = 1, ngrid
323 zmax(ig) = max(zmax(ig), zlev(ig,lmax(ig,k))-zlev(ig,k))
324 END DO
325 ! print*,k,lmax(1,k)
326 END DO
327 ! print*,'ZMAX ZMAX ZMAX ',zmax
328 ! call dump2d(iim,jjm-1,zmax(2:ngrid-1),'ZMAX ')
329
330 ! print*,'OKl336'
331 ! Calcul de l'entrainement.
332 ! Le rapport d'aspect relie la largeur de l'ascendance a l'epaisseur
333 ! de la couche d'alimentation en partant du principe que la vitesse
334 ! maximum dans l'ascendance est la vitesse d'entrainement horizontale.
335 DO k = 1, nlay
336 DO ig = 1, ngrid
337 zzz = rho(ig, k)*wmax(ig, k)*(zlev(ig,k+1)-zlev(ig,k))/ &
338 (zmax(ig)*r_aspect)
339 IF (w2di==2) THEN
340 entr(ig, k) = entr(ig, k) + ptimestep*(zzz-entr(ig,k))/tho
341 ELSE
342 entr(ig, k) = zzz
343 END IF
344 ztva(ig, k) = ztv(ig, k)
345 END DO
346 END DO
347
348
349 ! print*,'7 OK convect8'
350 DO k = 1, klev + 1
351 DO ig = 1, ngrid
352 zw2(ig, k) = 0.
353 fmc(ig, k) = 0.
354 larg_cons(ig, k) = 0.
355 larg_detr(ig, k) = 0.
356 wa_moy(ig, k) = 0.
357 END DO
358 END DO
359
360 ! print*,'8 OK convect8'
361 DO ig = 1, ngrid
362 lmaxa(ig) = 1
363 lmix(ig) = 1
364 wmaxa(ig) = 0.
365 END DO
366
367
368 ! print*,'OKl372'
369 DO l = 1, nlay - 2
370 DO ig = 1, ngrid
371 ! if (zw2(ig,l).lt.1.e-10.and.ztv(ig,l).gt.ztv(ig,l+1)) then
372 ! print*,'COUCOU ',l,zw2(ig,l),ztv(ig,l),ztv(ig,l+1)
373 IF (zw2(ig,l)<1.E-10 .AND. ztv(ig,l)>ztv(ig,l+1) .AND. &
374 entr(ig,l)>1.E-10) THEN
375 ! print*,'COUCOU cas 1'
376 ! Initialisation de l'ascendance
377 ! lmix(ig)=1
378 ztva(ig, l) = ztv(ig, l)
379 fmc(ig, l) = 0.
380 fmc(ig, l+1) = entr(ig, l)
381 zw2(ig, l) = 0.
382 ! if (.not.ztv(ig,l+1).gt.150.) then
383 ! print*,'ig,l+1,ztv(ig,l+1)'
384 ! print*, ig,l+1,ztv(ig,l+1)
385 ! endif
386 zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
387 (zlev(ig,l+1)-zlev(ig,l))
388 larg_detr(ig, l) = 0.
389 ELSE IF (zw2(ig,l)>=1.E-10 .AND. fmc(ig,l)+entr(ig,l)>1.E-10) THEN
390 ! Incrementation...
391 fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
392 ! if (.not.fmc(ig,l+1).gt.1.e-15) then
393 ! print*,'ig,l+1,fmc(ig,l+1)'
394 ! print*, ig,l+1,fmc(ig,l+1)
395 ! print*,'Fmc ',(fmc(ig,ll),ll=1,klev+1)
396 ! print*,'W2 ',(zw2(ig,ll),ll=1,klev+1)
397 ! print*,'Tv ',(ztv(ig,ll),ll=1,klev)
398 ! print*,'Entr ',(entr(ig,ll),ll=1,klev)
399 ! endif
400 ztva(ig, l) = (fmc(ig,l)*ztva(ig,l-1)+entr(ig,l)*ztv(ig,l))/ &
401 fmc(ig, l+1)
402 ! mise a jour de la vitesse ascendante (l'air entraine de la couche
403 ! consideree commence avec une vitesse nulle).
404 zw2(ig, l+1) = zw2(ig, l)*(fmc(ig,l)/fmc(ig,l+1))**2 + &
405 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
406 END IF
407 IF (zw2(ig,l+1)<0.) THEN
408 zw2(ig, l+1) = 0.
409 lmaxa(ig) = l
410 ELSE
411 wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
412 END IF
413 IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
414 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
415 lmix(ig) = l + 1
416 wmaxa(ig) = wa_moy(ig, l+1)
417 END IF
418 ! print*,'COUCOU cas 2 LMIX=',lmix(ig),wa_moy(ig,l+1),wmaxa(ig)
419 END DO
420 END DO
421
422 ! print*,'9 OK convect8'
423 ! print*,'WA1 ',wa_moy
424
425 ! determination de l'indice du debut de la mixed layer ou w decroit
426
427 ! calcul de la largeur de chaque ascendance dans le cas conservatif.
428 ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
429 ! d'une couche est �gale � la hauteur de la couche alimentante.
430 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
431 ! de la vitesse d'entrainement horizontal dans la couche alimentante.
432
433 ! print*,'OKl439'
434 DO l = 2, nlay
435 DO ig = 1, ngrid
436 IF (l<=lmaxa(ig)) THEN
437 zw = max(wa_moy(ig,l), 1.E-10)
438 larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw)
439 END IF
440 END DO
441 END DO
442
443 DO l = 2, nlay
444 DO ig = 1, ngrid
445 IF (l<=lmaxa(ig)) THEN
446 ! if (idetr.eq.0) then
447 ! cette option est finalement en dur.
448 larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
449 ! else if (idetr.eq.1) then
450 ! larg_detr(ig,l)=larg_cons(ig,l)
451 ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
452 ! else if (idetr.eq.2) then
453 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
454 ! s *sqrt(wa_moy(ig,l))
455 ! else if (idetr.eq.4) then
456 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
457 ! s *wa_moy(ig,l)
458 ! endif
459 END IF
460 END DO
461 END DO
462
463 ! print*,'10 OK convect8'
464 ! print*,'WA2 ',wa_moy
465 ! calcul de la fraction de la maille concern�e par l'ascendance en tenant
466 ! compte de l'epluchage du thermique.
467
468 DO l = 2, nlay
469 DO ig = 1, ngrid
470 IF (larg_cons(ig,l)>1.) THEN
471 ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
472 fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
473 IF (l>lmix(ig)) THEN
474 xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
475 IF (idetr==0) THEN
476 fraca(ig, l) = fraca(ig, lmix(ig))
477 ELSE IF (idetr==1) THEN
478 fraca(ig, l) = fraca(ig, lmix(ig))*xxx(ig, l)
479 ELSE IF (idetr==2) THEN
480 fraca(ig, l) = fraca(ig, lmix(ig))*(1.-(1.-xxx(ig,l))**2)
481 ELSE
482 fraca(ig, l) = fraca(ig, lmix(ig))*xxx(ig, l)**2
483 END IF
484 END IF
485 ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
486 fraca(ig, l) = max(fraca(ig,l), 0.)
487 fraca(ig, l) = min(fraca(ig,l), 0.5)
488 fracd(ig, l) = 1. - fraca(ig, l)
489 fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
490 ELSE
491 ! wa_moy(ig,l)=0.
492 fraca(ig, l) = 0.
493 fracc(ig, l) = 0.
494 fracd(ig, l) = 1.
495 END IF
496 END DO
497 END DO
498
499 ! print*,'11 OK convect8'
500 ! print*,'Ea3 ',wa_moy
501 ! ------------------------------------------------------------------
502 ! Calcul de fracd, wd
503 ! somme wa - wd = 0
504 ! ------------------------------------------------------------------
505
506
507 DO ig = 1, ngrid
508 fm(ig, 1) = 0.
509 fm(ig, nlay+1) = 0.
510 END DO
511
512 DO l = 2, nlay
513 DO ig = 1, ngrid
514 fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
515 END DO
516 DO ig = 1, ngrid
517 IF (fracd(ig,l)<0.1) THEN
518 abort_message = 'fracd trop petit'
519 CALL abort_physic(modname, abort_message, 1)
520 ELSE
521 ! vitesse descendante "diagnostique"
522 wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
523 END IF
524 END DO
525 END DO
526
527 DO l = 1, nlay
528 DO ig = 1, ngrid
529 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
530 masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
531 END DO
532 END DO
533
534 ! print*,'12 OK convect8'
535 ! print*,'WA4 ',wa_moy
536 ! c------------------------------------------------------------------
537 ! calcul du transport vertical
538 ! ------------------------------------------------------------------
539
540 GO TO 4444
541 ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
542 DO l = 2, nlay - 1
543 DO ig = 1, ngrid
544 IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( &
545 ig,l+1)) THEN
546 ! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
547 ! s ,fm(ig,l+1)*ptimestep
548 ! s ,' M=',masse(ig,l),masse(ig,l+1)
549 END IF
550 END DO
551 END DO
552
553 DO l = 1, nlay
554 DO ig = 1, ngrid
555 IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
556 ! print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
557 ! s ,entr(ig,l)*ptimestep
558 ! s ,' M=',masse(ig,l)
559 END IF
560 END DO
561 END DO
562
563 DO l = 1, nlay
564 DO ig = 1, ngrid
565 IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
566 ! print*,'WARN!!! fm exagere ig=',ig,' l=',l
567 ! s ,' FM=',fm(ig,l)
568 END IF
569 IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
570 ! print*,'WARN!!! masse exagere ig=',ig,' l=',l
571 ! s ,' M=',masse(ig,l)
572 ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
573 ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
574 ! print*,'zlev(ig,l+1),zlev(ig,l)'
575 ! s ,zlev(ig,l+1),zlev(ig,l)
576 ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
577 ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
578 END IF
579 IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
580 ! print*,'WARN!!! entr exagere ig=',ig,' l=',l
581 ! s ,' E=',entr(ig,l)
582 END IF
583 END DO
584 END DO
585
586 4444 CONTINUE
587 ! print*,'OK 444 '
588
589 IF (w2di==1) THEN
590 fm0 = fm0 + ptimestep*(fm-fm0)/tho
591 entr0 = entr0 + ptimestep*(entr-entr0)/tho
592 ELSE
593 fm0 = fm
594 entr0 = entr
595 END IF
596
597 IF (flagdq==0) THEN
598 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, &
599 zha)
600 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, &
601 zoa)
602 PRINT *, 'THERMALS OPT 1'
603 ELSE IF (flagdq==1) THEN
604 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
605 zdhadj, zha)
606 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
607 pdoadj, zoa)
608 PRINT *, 'THERMALS OPT 2'
609 ELSE
610 CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zh, &
611 zdhadj, zha, lev_out)
612 CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zo, &
613 pdoadj, zoa, lev_out)
614 PRINT *, 'THERMALS OPT 3', dqimpl
615 END IF
616
617 PRINT *, 'TH VENT ', dvdq
618 IF (dvdq==0) THEN
619 ! print*,'TH VENT OK ',dvdq
620 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
621 zua)
622 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
623 zva)
624 ELSE IF (dvdq==1) THEN
625 CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
626 zu, zv, pduadj, pdvadj, zua, zva)
627 ELSE IF (dvdq==2) THEN
628 CALL thermcell_dv2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, &
629 zmax, zu, zv, pduadj, pdvadj, zua, zva, lev_out)
630 ELSE IF (dvdq==3) THEN
631 CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zu, &
632 pduadj, zua, lev_out)
633 CALL thermcell_dq(ngrid, nlay, dqimpl, ptimestep, fm0, entr0, masse, zv, &
634 pdvadj, zva, lev_out)
635 END IF
636
637 ! CALL writefield_phy('duadj',pduadj,klev)
638
639 DO l = 1, nlay
640 DO ig = 1, ngrid
641 zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
642 zf2 = zf/(1.-zf)
643 thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l))**2
644 wth2(ig, l) = zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
645 END DO
646 END DO
647
648
649
650 ! print*,'13 OK convect8'
651 ! print*,'WA5 ',wa_moy
652 DO l = 1, nlay
653 DO ig = 1, ngrid
654 pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l)
655 END DO
656 END DO
657
658
659 ! do l=1,nlay
660 ! do ig=1,ngrid
661 ! if(abs(pdtadj(ig,l))*86400..gt.500.) then
662 ! print*,'WARN!!! ig=',ig,' l=',l
663 ! s ,' pdtadj=',pdtadj(ig,l)
664 ! endif
665 ! if(abs(pdoadj(ig,l))*86400..gt.1.) then
666 ! print*,'WARN!!! ig=',ig,' l=',l
667 ! s ,' pdoadj=',pdoadj(ig,l)
668 ! endif
669 ! enddo
670 ! enddo
671
672 ! print*,'14 OK convect8'
673 ! ------------------------------------------------------------------
674 ! Calculs pour les sorties
675 ! ------------------------------------------------------------------
676
677 IF (sorties) THEN
678 DO l = 1, nlay
679 DO ig = 1, ngrid
680 zla(ig, l) = (1.-fracd(ig,l))*zmax(ig)
681 zld(ig, l) = fracd(ig, l)*zmax(ig)
682 IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ &
683 (1.-fracd(ig,l))
684 END DO
685 END DO
686
687 DO l = 1, nlay
688 DO ig = 1, ngrid
689 detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1)
690 IF (detr(ig,l)<0.) THEN
691 entr(ig, l) = entr(ig, l) - detr(ig, l)
692 detr(ig, l) = 0.
693 ! print*,'WARNING !!! detrainement negatif ',ig,l
694 END IF
695 END DO
696 END DO
697 END IF
698
699 ! print*,'15 OK convect8'
700
701
702 ! if(wa_moy(1,4).gt.1.e-10) stop
703
704 ! print*,'19 OK convect8'
705 RETURN
706 END SUBROUTINE thermcell_2002
707
708 SUBROUTINE thermcell_cld(ngrid, nlay, ptimestep, pplay, pplev, pphi, zlev, &
709 debut, pu, pv, pt, po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0, zqla, &
710 lmax, zmax_sec, wmax_sec, zw_sec, lmix_sec, ratqscth, ratqsdiff & ! s
711 ! ,pu_therm,pv_therm
712 , r_aspect, l_mix, w2di, tho)
713
714 USE dimphy
715 IMPLICIT NONE
716
717 ! =======================================================================
718
719 ! Calcul du transport verticale dans la couche limite en presence
720 ! de "thermiques" explicitement representes
721
722 ! R��criture � partir d'un listing papier � Habas, le 14/02/00
723
724 ! le thermique est suppos� homog�ne et dissip� par m�lange avec
725 ! son environnement. la longueur l_mix contr�le l'efficacit� du
726 ! m�lange
727
728 ! Le calcul du transport des diff�rentes esp�ces se fait en prenant
729 ! en compte:
730 ! 1. un flux de masse montant
731 ! 2. un flux de masse descendant
732 ! 3. un entrainement
733 ! 4. un detrainement
734
735 ! =======================================================================
736
737 ! -----------------------------------------------------------------------
738 ! declarations:
739 ! -------------
740
741 include "YOMCST.h"
742 include "YOETHF.h"
743 include "FCTTRE.h"
744
745 ! arguments:
746 ! ----------
747
748 INTEGER ngrid, nlay, w2di
749 REAL tho
750 REAL ptimestep, l_mix, r_aspect
751 REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
752 REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
753 REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
754 REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
755 REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
756 REAL pphi(ngrid, nlay)
757
758 INTEGER idetr
759 SAVE idetr
760 DATA idetr/3/
761 !$OMP THREADPRIVATE(idetr)
762
763 ! local:
764 ! ------
765
766 INTEGER ig, k, l, lmaxa(klon), lmix(klon)
767 REAL zsortie1d(klon)
768 ! CR: on remplace lmax(klon,klev+1)
769 INTEGER lmax(klon), lmin(klon), lentr(klon)
770 REAL linter(klon)
771 REAL zmix(klon), fracazmix(klon)
772 REAL alpha
773 SAVE alpha
774 DATA alpha/1./
775 !$OMP THREADPRIVATE(alpha)
776
777 ! RC
778 REAL zmax(klon), zw, zz, zw2(klon, klev+1), ztva(klon, klev), zzz
779 REAL zmax_sec(klon)
780 REAL zmax_sec2(klon)
781 REAL zw_sec(klon, klev+1)
782 INTEGER lmix_sec(klon)
783 REAL w_est(klon, klev+1)
784 ! on garde le zmax du pas de temps precedent
785 ! real zmax0(klon)
786 ! save zmax0
787 ! real zmix0(klon)
788 ! save zmix0
789 REAL, SAVE, ALLOCATABLE :: zmax0(:), zmix0(:)
790 !$OMP THREADPRIVATE(zmax0, zmix0)
791
792 REAL zlev(klon, klev+1), zlay(klon, klev)
793 REAL deltaz(klon, klev)
794 REAL zh(klon, klev), zdhadj(klon, klev)
795 REAL zthl(klon, klev), zdthladj(klon, klev)
796 REAL ztv(klon, klev)
797 REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
798 REAL zl(klon, klev)
799 REAL wh(klon, klev+1)
800 REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
801 REAL zla(klon, klev+1)
802 REAL zwa(klon, klev+1)
803 REAL zld(klon, klev+1)
804 REAL zwd(klon, klev+1)
805 REAL zsortie(klon, klev)
806 REAL zva(klon, klev)
807 REAL zua(klon, klev)
808 REAL zoa(klon, klev)
809
810 REAL zta(klon, klev)
811 REAL zha(klon, klev)
812 REAL wa_moy(klon, klev+1)
813 REAL fraca(klon, klev+1)
814 REAL fracc(klon, klev+1)
815 REAL zf, zf2
816 REAL thetath2(klon, klev), wth2(klon, klev), wth3(klon, klev)
817 REAL q2(klon, klev)
818 REAL dtheta(klon, klev)
819 ! common/comtherm/thetath2,wth2
820
821 REAL ratqscth(klon, klev)
822 REAL sum
823 REAL sumdiff
824 REAL ratqsdiff(klon, klev)
825 REAL count_time
826 INTEGER ialt
827
828 LOGICAL sorties
829 REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
830 REAL zpspsk(klon, klev)
831
832 ! real wmax(klon,klev),wmaxa(klon)
833 REAL wmax(klon), wmaxa(klon)
834 REAL wmax_sec(klon)
835 REAL wmax_sec2(klon)
836 REAL wa(klon, klev, klev+1)
837 REAL wd(klon, klev+1)
838 REAL larg_part(klon, klev, klev+1)
839 REAL fracd(klon, klev+1)
840 REAL xxx(klon, klev+1)
841 REAL larg_cons(klon, klev+1)
842 REAL larg_detr(klon, klev+1)
843 REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
844 REAL massetot(klon, klev)
845 REAL detr0(klon, klev)
846 REAL alim0(klon, klev)
847 REAL pu_therm(klon, klev), pv_therm(klon, klev)
848 REAL fm(klon, klev+1), entr(klon, klev)
849 REAL fmc(klon, klev+1)
850
851 REAL zcor, zdelta, zcvm5, qlbef
852 REAL tbef(klon), qsatbef(klon)
853 REAL dqsat_dt, dt, num, denom
854 REAL reps, rlvcp, ddt0
855 REAL ztla(klon, klev), zqla(klon, klev), zqta(klon, klev)
856 ! CR niveau de condensation
857 REAL nivcon(klon)
858 REAL zcon(klon)
859 REAL zqsat(klon, klev)
860 REAL zqsatth(klon, klev)
861 PARAMETER (ddt0=.01)
862
863
864 ! CR:nouvelles variables
865 REAL f_star(klon, klev+1), entr_star(klon, klev)
866 REAL detr_star(klon, klev)
867 REAL alim_star_tot(klon), alim_star2(klon)
868 REAL entr_star_tot(klon)
869 REAL detr_star_tot(klon)
870 REAL alim_star(klon, klev)
871 REAL alim(klon, klev)
872 REAL nu(klon, klev)
873 REAL nu_e(klon, klev)
874 REAL nu_min
875 REAL nu_max
876 REAL nu_r
877 REAL f(klon)
878 ! real f(klon), f0(klon)
879 ! save f0
880 REAL, SAVE, ALLOCATABLE :: f0(:)
881 !$OMP THREADPRIVATE(f0)
882
883 REAL f_old
884 REAL zlevinter(klon)
885 LOGICAL, SAVE :: first = .TRUE.
886 !$OMP THREADPRIVATE(first)
887 ! data first /.false./
888 ! save first
889 LOGICAL nuage
890 ! save nuage
891 LOGICAL boucle
892 LOGICAL therm
893 LOGICAL debut
894 LOGICAL rale
895 INTEGER test(klon)
896 INTEGER signe_zw2
897 ! RC
898
899 CHARACTER *2 str2
900 CHARACTER *10 str10
901
902 CHARACTER (LEN=20) :: modname = 'thermcell_cld'
903 CHARACTER (LEN=80) :: abort_message
904
905 LOGICAL vtest(klon), down
906 LOGICAL zsat(klon)
907
908 EXTERNAL scopy
909
910 INTEGER ncorrec, ll
911 SAVE ncorrec
912 DATA ncorrec/0/
913 !$OMP THREADPRIVATE(ncorrec)
914
915
916
917 ! -----------------------------------------------------------------------
918 ! initialisation:
919 ! ---------------
920
921 IF (first) THEN
922 ALLOCATE (zmix0(klon))
923 ALLOCATE (zmax0(klon))
924 ALLOCATE (f0(klon))
925 first = .FALSE.
926 END IF
927
928 sorties = .FALSE.
929 ! print*,'NOUVEAU DETR PLUIE '
930 IF (ngrid/=klon) THEN
931 PRINT *
932 PRINT *, 'STOP dans convadj'
933 PRINT *, 'ngrid =', ngrid
934 PRINT *, 'klon =', klon
935 END IF
936
937 ! Initialisation
938 rlvcp = rlvtt/rcpd
939 reps = rd/rv
940 ! initialisations de zqsat
941 DO ll = 1, nlay
942 DO ig = 1, ngrid
943 zqsat(ig, ll) = 0.
944 zqsatth(ig, ll) = 0.
945 END DO
946 END DO
947
948 ! on met le first a true pour le premier passage de la journ�e
949 DO ig = 1, klon
950 test(ig) = 0
951 END DO
952 IF (debut) THEN
953 DO ig = 1, klon
954 test(ig) = 1
955 f0(ig) = 0.
956 zmax0(ig) = 0.
957 END DO
958 END IF
959 DO ig = 1, klon
960 IF ((.NOT. debut) .AND. (f0(ig)<1.E-10)) THEN
961 test(ig) = 1
962 END IF
963 END DO
964 ! do ig=1,klon
965 ! print*,'test(ig)',test(ig),zmax0(ig)
966 ! enddo
967 nuage = .FALSE.
968 ! -----------------------------------------------------------------------
969 ! AM Calcul de T,q,ql a partir de Tl et qT
970 ! ---------------------------------------------------
971
972 ! Pr Tprec=Tl calcul de qsat
973 ! Si qsat>qT T=Tl, q=qT
974 ! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
975 ! On cherche DDT < DDT0
976
977 ! defaut
978 DO ll = 1, nlay
979 DO ig = 1, ngrid
980 zo(ig, ll) = po(ig, ll)
981 zl(ig, ll) = 0.
982 zh(ig, ll) = pt(ig, ll)
983 END DO
984 END DO
985 DO ig = 1, ngrid
986 zsat(ig) = .FALSE.
987 END DO
988
989
990 DO ll = 1, nlay
991 ! les points insatures sont definitifs
992 DO ig = 1, ngrid
993 tbef(ig) = pt(ig, ll)
994 zdelta = max(0., sign(1.,rtt-tbef(ig)))
995 qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, ll)
996 qsatbef(ig) = min(0.5, qsatbef(ig))
997 zcor = 1./(1.-retv*qsatbef(ig))
998 qsatbef(ig) = qsatbef(ig)*zcor
999 zsat(ig) = (max(0.,po(ig,ll)-qsatbef(ig))>1.E-10)
1000 END DO
1001
1002 DO ig = 1, ngrid
1003 IF (zsat(ig) .AND. (1==1)) THEN
1004 qlbef = max(0., po(ig,ll)-qsatbef(ig))
1005 ! si sature: ql est surestime, d'ou la sous-relax
1006 dt = 0.5*rlvcp*qlbef
1007 ! write(18,*),'DT0=',DT
1008 ! on pourra enchainer 2 ou 3 calculs sans Do while
1009 DO WHILE (abs(dt)>ddt0)
1010 ! il faut verifier si c,a conserve quand on repasse en insature ...
1011 tbef(ig) = tbef(ig) + dt
1012 zdelta = max(0., sign(1.,rtt-tbef(ig)))
1013 qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, ll)
1014 qsatbef(ig) = min(0.5, qsatbef(ig))
1015 zcor = 1./(1.-retv*qsatbef(ig))
1016 qsatbef(ig) = qsatbef(ig)*zcor
1017 ! on veut le signe de qlbef
1018 qlbef = po(ig, ll) - qsatbef(ig)
1019 zdelta = max(0., sign(1.,rtt-tbef(ig)))
1020 zcvm5 = r5les*(1.-zdelta) + r5ies*zdelta
1021 zcor = 1./(1.-retv*qsatbef(ig))
1022 dqsat_dt = foede(tbef(ig), zdelta, zcvm5, qsatbef(ig), zcor)
1023 num = -tbef(ig) + pt(ig, ll) + rlvcp*qlbef
1024 denom = 1. + rlvcp*dqsat_dt
1025 IF (denom<1.E-10) THEN
1026 PRINT *, 'pb denom'
1027 END IF
1028 dt = num/denom
1029 END DO
1030 ! on ecrit de maniere conservative (sat ou non)
1031 zl(ig, ll) = max(0., qlbef)
1032 ! T = Tl +Lv/Cp ql
1033 zh(ig, ll) = pt(ig, ll) + rlvcp*zl(ig, ll)
1034 zo(ig, ll) = po(ig, ll) - zl(ig, ll)
1035 END IF
1036 ! on ecrit zqsat
1037 zqsat(ig, ll) = qsatbef(ig)
1038 END DO
1039 END DO
1040 ! AM fin
1041
1042 ! -----------------------------------------------------------------------
1043 ! incrementation eventuelle de tendances precedentes:
1044 ! ---------------------------------------------------
1045
1046 ! print*,'0 OK convect8'
1047
1048 DO l = 1, nlay
1049 DO ig = 1, ngrid
1050 zpspsk(ig, l) = (pplay(ig,l)/100000.)**rkappa
1051 ! zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
1052 ! zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
1053 zu(ig, l) = pu(ig, l)
1054 zv(ig, l) = pv(ig, l)
1055 ! zo(ig,l)=po(ig,l)
1056 ! ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
1057 ! AM attention zh est maintenant le profil de T et plus le profil de
1058 ! theta !
1059
1060 ! T-> Theta
1061 ztv(ig, l) = zh(ig, l)/zpspsk(ig, l)
1062 ! AM Theta_v
1063 ztv(ig, l) = ztv(ig, l)*(1.+retv*(zo(ig,l))-zl(ig,l))
1064 ! AM Thetal
1065 zthl(ig, l) = pt(ig, l)/zpspsk(ig, l)
1066
1067 END DO
1068 END DO
1069
1070 ! print*,'1 OK convect8'
1071 ! --------------------
1072
1073
1074 ! + + + + + + + + + + +
1075
1076
1077 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
1078 ! wh,wt,wo ...
1079
1080 ! + + + + + + + + + + + zh,zu,zv,zo,rho
1081
1082
1083 ! -------------------- zlev(1)
1084 ! \\\\\\\\\\\\\\\\\\\
1085
1086
1087 ! -----------------------------------------------------------------------
1088 ! Calcul des altitudes des couches
1089 ! -----------------------------------------------------------------------
1090
1091 DO l = 2, nlay
1092 DO ig = 1, ngrid
1093 zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
1094 END DO
1095 END DO
1096 DO ig = 1, ngrid
1097 zlev(ig, 1) = 0.
1098 zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
1099 END DO
1100 DO l = 1, nlay
1101 DO ig = 1, ngrid
1102 zlay(ig, l) = pphi(ig, l)/rg
1103 END DO
1104 END DO
1105 ! calcul de deltaz
1106 DO l = 1, nlay
1107 DO ig = 1, ngrid
1108 deltaz(ig, l) = zlev(ig, l+1) - zlev(ig, l)
1109 END DO
1110 END DO
1111
1112 ! print*,'2 OK convect8'
1113 ! -----------------------------------------------------------------------
1114 ! Calcul des densites
1115 ! -----------------------------------------------------------------------
1116
1117 DO l = 1, nlay
1118 DO ig = 1, ngrid
1119 ! rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
1120 rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*ztv(ig,l))
1121 END DO
1122 END DO
1123
1124 DO l = 2, nlay
1125 DO ig = 1, ngrid
1126 rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
1127 END DO
1128 END DO
1129
1130 DO k = 1, nlay
1131 DO l = 1, nlay + 1
1132 DO ig = 1, ngrid
1133 wa(ig, k, l) = 0.
1134 END DO
1135 END DO
1136 END DO
1137 ! Cr:ajout:calcul de la masse
1138 DO l = 1, nlay
1139 DO ig = 1, ngrid
1140 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1141 masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
1142 END DO
1143 END DO
1144 ! print*,'3 OK convect8'
1145 ! ------------------------------------------------------------------
1146 ! Calcul de w2, quarre de w a partir de la cape
1147 ! a partir de w2, on calcule wa, vitesse de l'ascendance
1148
1149 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
1150 ! w2 est stoke dans wa
1151
1152 ! ATTENTION: dans convect8, on n'utilise le calcule des wa
1153 ! independants par couches que pour calculer l'entrainement
1154 ! a la base et la hauteur max de l'ascendance.
1155
1156 ! Indicages:
1157 ! l'ascendance provenant du niveau k traverse l'interface l avec
1158 ! une vitesse wa(k,l).
1159
1160 ! --------------------
1161
1162 ! + + + + + + + + + +
1163
1164 ! wa(k,l) ---- -------------------- l
1165 ! / ! /||\ + + + + + + + + + +
1166 ! ||
1167 ! || --------------------
1168 ! ||
1169 ! || + + + + + + + + + +
1170 ! ||
1171 ! || --------------------
1172 ! ||__
1173 ! |___ + + + + + + + + + + k
1174
1175 ! --------------------
1176
1177
1178
1179 ! ------------------------------------------------------------------
1180
1181 ! CR: ponderation entrainement des couches instables
1182 ! def des alim_star tels que alim=f*alim_star
1183 DO l = 1, klev
1184 DO ig = 1, ngrid
1185 alim_star(ig, l) = 0.
1186 alim(ig, l) = 0.
1187 END DO
1188 END DO
1189 ! determination de la longueur de la couche d entrainement
1190 DO ig = 1, ngrid
1191 lentr(ig) = 1
1192 END DO
1193
1194 ! on ne considere que les premieres couches instables
1195 therm = .FALSE.
1196 DO k = nlay - 2, 1, -1
1197 DO ig = 1, ngrid
1198 IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN
1199 lentr(ig) = k + 1
1200 therm = .TRUE.
1201 END IF
1202 END DO
1203 END DO
1204
1205 ! determination du lmin: couche d ou provient le thermique
1206 DO ig = 1, ngrid
1207 lmin(ig) = 1
1208 END DO
1209 DO ig = 1, ngrid
1210 DO l = nlay, 2, -1
1211 IF (ztv(ig,l-1)>ztv(ig,l)) THEN
1212 lmin(ig) = l - 1
1213 END IF
1214 END DO
1215 END DO
1216
1217 ! definition de l'entrainement des couches
1218 DO l = 1, klev - 1
1219 DO ig = 1, ngrid
1220 IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<lentr(ig)) THEN
1221 ! def possibles pour alim_star: zdthetadz, dthetadz, zdtheta
1222 alim_star(ig, l) = max((ztv(ig,l)-ztv(ig,l+1)), 0.) & ! s
1223 ! *(zlev(ig,l+1)-zlev(ig,l))
1224 *sqrt(zlev(ig,l+1))
1225 ! alim_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
1226 ! s /zlev(ig,lentr(ig)+2)))**(3./2.)
1227 END IF
1228 END DO
1229 END DO
1230
1231 ! pas de thermique si couche 1 stable
1232 DO ig = 1, ngrid
1233 ! if (lmin(ig).gt.1) then
1234 ! CRnouveau test
1235 IF (alim_star(ig,1)<1.E-10) THEN
1236 DO l = 1, klev
1237 alim_star(ig, l) = 0.
1238 END DO
1239 END IF
1240 END DO
1241 ! calcul de l entrainement total
1242 DO ig = 1, ngrid
1243 alim_star_tot(ig) = 0.
1244 entr_star_tot(ig) = 0.
1245 detr_star_tot(ig) = 0.
1246 END DO
1247 DO ig = 1, ngrid
1248 DO k = 1, klev
1249 alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, k)
1250 END DO
1251 END DO
1252
1253 ! Calcul entrainement normalise
1254 DO ig = 1, ngrid
1255 IF (alim_star_tot(ig)>1.E-10) THEN
1256 ! do l=1,lentr(ig)
1257 DO l = 1, klev
1258 ! def possibles pour entr_star: zdthetadz, dthetadz, zdtheta
1259 alim_star(ig, l) = alim_star(ig, l)/alim_star_tot(ig)
1260 END DO
1261 END IF
1262 END DO
1263
1264 ! print*,'fin calcul alim_star'
1265
1266 ! AM:initialisations
1267 DO k = 1, nlay
1268 DO ig = 1, ngrid
1269 ztva(ig, k) = ztv(ig, k)
1270 ztla(ig, k) = zthl(ig, k)
1271 zqla(ig, k) = 0.
1272 zqta(ig, k) = po(ig, k)
1273 zsat(ig) = .FALSE.
1274 END DO
1275 END DO
1276 DO k = 1, klev
1277 DO ig = 1, ngrid
1278 detr_star(ig, k) = 0.
1279 entr_star(ig, k) = 0.
1280 detr(ig, k) = 0.
1281 entr(ig, k) = 0.
1282 END DO
1283 END DO
1284 ! print*,'7 OK convect8'
1285 DO k = 1, klev + 1
1286 DO ig = 1, ngrid
1287 zw2(ig, k) = 0.
1288 fmc(ig, k) = 0.
1289 ! CR
1290 f_star(ig, k) = 0.
1291 ! RC
1292 larg_cons(ig, k) = 0.
1293 larg_detr(ig, k) = 0.
1294 wa_moy(ig, k) = 0.
1295 END DO
1296 END DO
1297
1298 ! n print*,'8 OK convect8'
1299 DO ig = 1, ngrid
1300 linter(ig) = 1.
1301 lmaxa(ig) = 1
1302 lmix(ig) = 1
1303 wmaxa(ig) = 0.
1304 END DO
1305
1306 nu_min = l_mix
1307 nu_max = 1000.
1308 ! do ig=1,ngrid
1309 ! nu_max=wmax_sec(ig)
1310 ! enddo
1311 DO ig = 1, ngrid
1312 DO k = 1, klev
1313 nu(ig, k) = 0.
1314 nu_e(ig, k) = 0.
1315 END DO
1316 END DO
1317 ! Calcul de l'exc�s de temp�rature du � la diffusion turbulente
1318 DO ig = 1, ngrid
1319 DO l = 1, klev
1320 dtheta(ig, l) = 0.
1321 END DO
1322 END DO
1323 DO ig = 1, ngrid
1324 DO l = 1, lentr(ig) - 1
1325 dtheta(ig, l) = sqrt(10.*0.4*zlev(ig,l+1)**2*1.*((ztv(ig,l+1)- &
1326 ztv(ig,l))/(zlev(ig,l+1)-zlev(ig,l)))**2)
1327 END DO
1328 END DO
1329 ! do l=1,nlay-2
1330 DO l = 1, klev - 1
1331 DO ig = 1, ngrid
1332 IF (ztv(ig,l)>ztv(ig,l+1) .AND. alim_star(ig,l)>1.E-10 .AND. &
1333 zw2(ig,l)<1E-10) THEN
1334 ! AM
1335 ! test:on rajoute un exc�s de T dans couche alim
1336 ! ztla(ig,l)=zthl(ig,l)+dtheta(ig,l)
1337 ztla(ig, l) = zthl(ig, l)
1338 ! test: on rajoute un exc�s de q dans la couche alim
1339 ! zqta(ig,l)=po(ig,l)+0.001
1340 zqta(ig, l) = po(ig, l)
1341 zqla(ig, l) = zl(ig, l)
1342 ! AM
1343 f_star(ig, l+1) = alim_star(ig, l)
1344 ! test:calcul de dteta
1345 zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
1346 (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
1347 w_est(ig, l+1) = zw2(ig, l+1)
1348 larg_detr(ig, l) = 0.
1349 ! print*,'coucou boucle 1'
1350 ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+alim_star(ig, &
1351 l))>1.E-10) THEN
1352 ! print*,'coucou boucle 2'
1353 ! estimation du detrainement a partir de la geometrie du pas
1354 ! precedent
1355 IF ((test(ig)==1) .OR. ((.NOT. debut) .AND. (f0(ig)<1.E-10))) THEN
1356 detr_star(ig, l) = 0.
1357 entr_star(ig, l) = 0.
1358 ! print*,'coucou test(ig)',test(ig),f0(ig),zmax0(ig)
1359 ELSE
1360 ! print*,'coucou debut detr'
1361 ! tests sur la definition du detr
1362 IF (zqla(ig,l-1)>1.E-10) THEN
1363 nuage = .TRUE.
1364 END IF
1365
1366 w_est(ig, l+1) = zw2(ig, l)*((f_star(ig,l))**2)/(f_star(ig,l)+ &
1367 alim_star(ig,l))**2 + 2.*rg*(ztva(ig,l-1)-ztv(ig,l))/ztv(ig, l)*( &
1368 zlev(ig,l+1)-zlev(ig,l))
1369 IF (w_est(ig,l+1)<0.) THEN
1370 w_est(ig, l+1) = zw2(ig, l)
1371 END IF
1372 IF (l>2) THEN
1373 IF ((w_est(ig,l+1)>w_est(ig,l)) .AND. (zlev(ig, &
1374 l+1)<zmax_sec(ig)) .AND. (zqla(ig,l-1)<1.E-10)) THEN
1375 detr_star(ig, l) = max(0., (rhobarz(ig, &
1376 l+1)*sqrt(w_est(ig,l+1))*sqrt(nu(ig,l)* &
1377 zlev(ig,l+1))-rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(nu(ig,l)* &
1378 zlev(ig,l)))/(r_aspect*zmax_sec(ig)))
1379 ELSE IF ((zlev(ig,l+1)<zmax_sec(ig)) .AND. (zqla(ig, &
1380 l-1)<1.E-10)) THEN
1381 detr_star(ig, l) = -f0(ig)*f_star(ig, lmix(ig))/(rhobarz(ig, &
1382 lmix(ig))*wmaxa(ig))*(rhobarz(ig,l+1)*sqrt(w_est(ig, &
1383 l+1))*((zmax_sec(ig)-zlev(ig,l+1))/((zmax_sec(ig)-zlev(ig, &
1384 lmix(ig)))))**2.-rhobarz(ig,l)*sqrt(w_est(ig, &
1385 l))*((zmax_sec(ig)-zlev(ig,l))/((zmax_sec(ig)-zlev(ig,lmix(ig &
1386 )))))**2.)
1387 ELSE
1388 detr_star(ig, l) = 0.002*f0(ig)*f_star(ig, l)* &
1389 (zlev(ig,l+1)-zlev(ig,l))
1390
1391 END IF
1392 ELSE
1393 detr_star(ig, l) = 0.
1394 END IF
1395
1396 detr_star(ig, l) = detr_star(ig, l)/f0(ig)
1397 IF (nuage) THEN
1398 entr_star(ig, l) = 0.4*detr_star(ig, l)
1399 ELSE
1400 entr_star(ig, l) = 0.4*detr_star(ig, l)
1401 END IF
1402
1403 IF ((detr_star(ig,l))>f_star(ig,l)) THEN
1404 detr_star(ig, l) = f_star(ig, l)
1405 ! entr_star(ig,l)=0.
1406 END IF
1407
1408 IF ((l<lentr(ig))) THEN
1409 entr_star(ig, l) = 0.
1410 ! detr_star(ig,l)=0.
1411 END IF
1412
1413 ! print*,'ok detr_star'
1414 END IF
1415 ! prise en compte du detrainement dans le calcul du flux
1416 f_star(ig, l+1) = f_star(ig, l) + alim_star(ig, l) + &
1417 entr_star(ig, l) - detr_star(ig, l)
1418 ! test
1419 ! if (f_star(ig,l+1).lt.0.) then
1420 ! f_star(ig,l+1)=0.
1421 ! entr_star(ig,l)=0.
1422 ! detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)
1423 ! endif
1424 ! test sur le signe de f_star
1425 IF (f_star(ig,l+1)>1.E-10) THEN
1426 ! then
1427 ! test
1428 ! if (((f_star(ig,l+1)+detr_star(ig,l)).gt.1.e-10)) then
1429 ! AM on melange Tl et qt du thermique
1430 ! on rajoute un exc�s de T dans la couche alim
1431 ! if (l.lt.lentr(ig)) then
1432 ! ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+
1433 ! s
1434 ! (alim_star(ig,l)+entr_star(ig,l))*(zthl(ig,l)+dtheta(ig,l)))
1435 ! s /(f_star(ig,l+1)+detr_star(ig,l))
1436 ! else
1437 ztla(ig, l) = (f_star(ig,l)*ztla(ig,l-1)+(alim_star(ig, &
1438 l)+entr_star(ig,l))*zthl(ig,l))/(f_star(ig,l+1)+detr_star(ig,l))
1439 ! s /(f_star(ig,l+1))
1440 ! endif
1441 ! on rajoute un exc�s de q dans la couche alim
1442 ! if (l.lt.lentr(ig)) then
1443 ! zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+
1444 ! s (alim_star(ig,l)+entr_star(ig,l))*(po(ig,l)+0.001))
1445 ! s /(f_star(ig,l+1)+detr_star(ig,l))
1446 ! else
1447 zqta(ig, l) = (f_star(ig,l)*zqta(ig,l-1)+(alim_star(ig, &
1448 l)+entr_star(ig,l))*po(ig,l))/(f_star(ig,l+1)+detr_star(ig,l))
1449 ! s /(f_star(ig,l+1))
1450 ! endif
1451 ! AM on en deduit thetav et ql du thermique
1452 ! CR test
1453 ! Tbef(ig)=ztla(ig,l)*zpspsk(ig,l)
1454 tbef(ig) = ztla(ig, l)*zpspsk(ig, l)
1455 zdelta = max(0., sign(1.,rtt-tbef(ig)))
1456 qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, l)
1457 qsatbef(ig) = min(0.5, qsatbef(ig))
1458 zcor = 1./(1.-retv*qsatbef(ig))
1459 qsatbef(ig) = qsatbef(ig)*zcor
1460 zsat(ig) = (max(0.,zqta(ig,l)-qsatbef(ig))>1.E-10)
1461
1462 IF (zsat(ig) .AND. (1==1)) THEN
1463 qlbef = max(0., zqta(ig,l)-qsatbef(ig))
1464 dt = 0.5*rlvcp*qlbef
1465 ! write(17,*)'DT0=',DT
1466 DO WHILE (abs(dt)>ddt0)
1467 ! print*,'aie'
1468 tbef(ig) = tbef(ig) + dt
1469 zdelta = max(0., sign(1.,rtt-tbef(ig)))
1470 qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, l)
1471 qsatbef(ig) = min(0.5, qsatbef(ig))
1472 zcor = 1./(1.-retv*qsatbef(ig))
1473 qsatbef(ig) = qsatbef(ig)*zcor
1474 qlbef = zqta(ig, l) - qsatbef(ig)
1475
1476 zdelta = max(0., sign(1.,rtt-tbef(ig)))
1477 zcvm5 = r5les*(1.-zdelta) + r5ies*zdelta
1478 zcor = 1./(1.-retv*qsatbef(ig))
1479 dqsat_dt = foede(tbef(ig), zdelta, zcvm5, qsatbef(ig), zcor)
1480 num = -tbef(ig) + ztla(ig, l)*zpspsk(ig, l) + rlvcp*qlbef
1481 denom = 1. + rlvcp*dqsat_dt
1482 IF (denom<1.E-10) THEN
1483 PRINT *, 'pb denom'
1484 END IF
1485 dt = num/denom
1486 ! write(17,*)'DT=',DT
1487 END DO
1488 zqla(ig, l) = max(0., zqta(ig,l)-qsatbef(ig))
1489 zqla(ig, l) = max(0., qlbef)
1490 ! zqla(ig,l)=0.
1491 END IF
1492 ! zqla(ig,l) = max(0.,zqta(ig,l)-qsatbef(ig))
1493
1494 ! on ecrit de maniere conservative (sat ou non)
1495 ! T = Tl +Lv/Cp ql
1496 ! CR rq utilisation de humidite specifique ou rapport de melange?
1497 ztva(ig, l) = ztla(ig, l)*zpspsk(ig, l) + rlvcp*zqla(ig, l)
1498 ztva(ig, l) = ztva(ig, l)/zpspsk(ig, l)
1499 ! on rajoute le calcul de zha pour diagnostiques (temp potentielle)
1500 zha(ig, l) = ztva(ig, l)
1501 ! if (l.lt.lentr(ig)) then
1502 ! ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)
1503 ! s -zqla(ig,l))-zqla(ig,l)) + 0.1
1504 ! else
1505 ztva(ig, l) = ztva(ig, l)*(1.+retv*(zqta(ig,l)-zqla(ig, &
1506 l))-zqla(ig,l))
1507 ! endif
1508 ! ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
1509 ! s /(1.-retv*zqla(ig,l))
1510 ! ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1511 ! ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)
1512 ! s /(1.-retv*zqta(ig,l))
1513 ! s -zqla(ig,l)/(1.-retv*zqla(ig,l)))
1514 ! s -zqla(ig,l)/(1.-retv*zqla(ig,l)))
1515 ! write(13,*)zqla(ig,l),zqla(ig,l)/(1.-retv*zqla(ig,l))
1516 ! on ecrit zqsat
1517 zqsatth(ig, l) = qsatbef(ig)
1518 ! enddo
1519 ! DO ig=1,ngrid
1520 ! if (zw2(ig,l).ge.1.e-10.and.
1521 ! s f_star(ig,l)+entr_star(ig,l).gt.1.e-10) then
1522 ! mise a jour de la vitesse ascendante (l'air entraine de la couche
1523 ! consideree commence avec une vitesse nulle).
1524
1525 ! if (f_star(ig,l+1).gt.1.e-10) then
1526 zw2(ig, l+1) = zw2(ig, l)* & ! s
1527 ! ((f_star(ig,l)-detr_star(ig,l))**2)
1528 ! s /f_star(ig,l+1)**2+
1529 ((f_star(ig,l))**2)/(f_star(ig,l+1)+detr_star(ig,l))**2 + & ! s
1530 ! /(f_star(ig,l+1))**2+
1531 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
1532 ! s *(f_star(ig,l)/f_star(ig,l+1))**2
1533
1534 END IF
1535 END IF
1536
1537 IF (zw2(ig,l+1)<0.) THEN
1538 linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
1539 ig,l))
1540 zw2(ig, l+1) = 0.
1541 ! print*,'linter=',linter(ig)
1542 ! else if ((zw2(ig,l+1).lt.1.e-10).and.(zw2(ig,l+1).ge.0.)) then
1543 ! linter(ig)=l+1
1544 ! print*,'linter=l',zw2(ig,l),zw2(ig,l+1)
1545 ELSE
1546 wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
1547 ! wa_moy(ig,l+1)=zw2(ig,l+1)
1548 END IF
1549 IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
1550 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
1551 lmix(ig) = l + 1
1552 wmaxa(ig) = wa_moy(ig, l+1)
1553 END IF
1554 END DO
1555 END DO
1556 PRINT *, 'fin calcul zw2'
1557
1558 ! Calcul de la couche correspondant a la hauteur du thermique
1559 DO ig = 1, ngrid
1560 lmax(ig) = lentr(ig)
1561 END DO
1562 DO ig = 1, ngrid
1563 DO l = nlay, lentr(ig) + 1, -1
1564 IF (zw2(ig,l)<=1.E-10) THEN
1565 lmax(ig) = l - 1
1566 END IF
1567 END DO
1568 END DO
1569 ! pas de thermique si couche 1 stable
1570 DO ig = 1, ngrid
1571 IF (lmin(ig)>1) THEN
1572 lmax(ig) = 1
1573 lmin(ig) = 1
1574 lentr(ig) = 1
1575 END IF
1576 END DO
1577
1578 ! Determination de zw2 max
1579 DO ig = 1, ngrid
1580 wmax(ig) = 0.
1581 END DO
1582
1583 DO l = 1, nlay
1584 DO ig = 1, ngrid
1585 IF (l<=lmax(ig)) THEN
1586 IF (zw2(ig,l)<0.) THEN
1587 PRINT *, 'pb2 zw2<0'
1588 END IF
1589 zw2(ig, l) = sqrt(zw2(ig,l))
1590 wmax(ig) = max(wmax(ig), zw2(ig,l))
1591 ELSE
1592 zw2(ig, l) = 0.
1593 END IF
1594 END DO
1595 END DO
1596
1597 ! Longueur caracteristique correspondant a la hauteur des thermiques.
1598 DO ig = 1, ngrid
1599 zmax(ig) = 0.
1600 zlevinter(ig) = zlev(ig, 1)
1601 END DO
1602 DO ig = 1, ngrid
1603 ! calcul de zlevinter
1604 zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
1605 zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
1606 ! pour le cas ou on prend tjs lmin=1
1607 ! zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
1608 zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,1))
1609 zmax0(ig) = zmax(ig)
1610 WRITE (11, *) 'ig,lmax,linter', ig, lmax(ig), linter(ig)
1611 WRITE (12, *) 'ig,zlevinter,zmax', ig, zmax(ig), zlevinter(ig)
1612 END DO
1613
1614 ! Calcul de zmax_sec et wmax_sec
1615 CALL fermeture_seche(ngrid, nlay, pplay, pplev, pphi, zlev, rhobarz, f0, &
1616 zpspsk, alim, zh, zo, lentr, lmin, nu_min, nu_max, r_aspect, zmax_sec2, &
1617 wmax_sec2)
1618
1619 PRINT *, 'avant fermeture'
1620 ! Fermeture,determination de f
1621 ! en lmax f=d-e
1622 DO ig = 1, ngrid
1623 ! entr_star(ig,lmax(ig))=0.
1624 ! f_star(ig,lmax(ig)+1)=0.
1625 ! detr_star(ig,lmax(ig))=f_star(ig,lmax(ig))+entr_star(ig,lmax(ig))
1626 ! s +alim_star(ig,lmax(ig))
1627 END DO
1628
1629 DO ig = 1, ngrid
1630 alim_star2(ig) = 0.
1631 END DO
1632 ! calcul de entr_star_tot
1633 DO ig = 1, ngrid
1634 DO k = 1, lmix(ig)
1635 entr_star_tot(ig) = entr_star_tot(ig) & ! s
1636 ! +entr_star(ig,k)
1637 +alim_star(ig, k)
1638 ! s -detr_star(ig,k)
1639 detr_star_tot(ig) = detr_star_tot(ig) & ! s
1640 ! +alim_star(ig,k)
1641 -detr_star(ig, k) + entr_star(ig, k)
1642 END DO
1643 END DO
1644
1645 DO ig = 1, ngrid
1646 IF (alim_star_tot(ig)<1.E-10) THEN
1647 f(ig) = 0.
1648 ELSE
1649 ! do k=lmin(ig),lentr(ig)
1650 DO k = 1, lentr(ig)
1651 alim_star2(ig) = alim_star2(ig) + alim_star(ig, k)**2/(rho(ig,k)*( &
1652 zlev(ig,k+1)-zlev(ig,k)))
1653 END DO
1654 IF ((zmax_sec(ig)>1.E-10) .AND. (1==1)) THEN
1655 f(ig) = wmax_sec(ig)/(max(500.,zmax_sec(ig))*r_aspect*alim_star2(ig))
1656 f(ig) = f(ig) + (f0(ig)-f(ig))*exp((-ptimestep/zmax_sec(ig))*wmax_sec &
1657 (ig))
1658 ELSE
1659 f(ig) = wmax(ig)/(max(500.,zmax(ig))*r_aspect*alim_star2(ig))
1660 f(ig) = f(ig) + (f0(ig)-f(ig))*exp((-ptimestep/zmax(ig))*wmax(ig))
1661 END IF
1662 END IF
1663 f0(ig) = f(ig)
1664 END DO
1665 PRINT *, 'apres fermeture'
1666 ! Calcul de l'entrainement
1667 DO ig = 1, ngrid
1668 DO k = 1, klev
1669 alim(ig, k) = f(ig)*alim_star(ig, k)
1670 END DO
1671 END DO
1672 ! CR:test pour entrainer moins que la masse
1673 ! do ig=1,ngrid
1674 ! do l=1,lentr(ig)
1675 ! if ((alim(ig,l)*ptimestep).gt.(0.9*masse(ig,l))) then
1676 ! alim(ig,l+1)=alim(ig,l+1)+alim(ig,l)
1677 ! s -0.9*masse(ig,l)/ptimestep
1678 ! alim(ig,l)=0.9*masse(ig,l)/ptimestep
1679 ! endif
1680 ! enddo
1681 ! enddo
1682 ! calcul du d�trainement
1683 DO ig = 1, klon
1684 DO k = 1, klev
1685 detr(ig, k) = f(ig)*detr_star(ig, k)
1686 IF (detr(ig,k)<0.) THEN
1687 ! print*,'detr1<0!!!'
1688 END IF
1689 END DO
1690 DO k = 1, klev
1691 entr(ig, k) = f(ig)*entr_star(ig, k)
1692 IF (entr(ig,k)<0.) THEN
1693 ! print*,'entr1<0!!!'
1694 END IF
1695 END DO
1696 END DO
1697
1698 ! do ig=1,ngrid
1699 ! do l=1,klev
1700 ! if (((detr(ig,l)+entr(ig,l)+alim(ig,l))*ptimestep).gt.
1701 ! s (masse(ig,l))) then
1702 ! print*,'d2+e2+a2>m2','ig=',ig,'l=',l,'lmax(ig)=',lmax(ig),'d+e+a='
1703 ! s,(detr(ig,l)+entr(ig,l)+alim(ig,l))*ptimestep,'m=',masse(ig,l)
1704 ! endif
1705 ! enddo
1706 ! enddo
1707 ! Calcul des flux
1708
1709 DO ig = 1, ngrid
1710 DO l = 1, lmax(ig)
1711 ! do l=1,klev
1712 ! fmc(ig,l+1)=f(ig)*f_star(ig,l+1)
1713 fmc(ig, l+1) = fmc(ig, l) + alim(ig, l) + entr(ig, l) - detr(ig, l)
1714 ! print*,'??!!','ig=',ig,'l=',l,'lmax=',lmax(ig),'lmix=',lmix(ig),
1715 ! s 'e=',entr(ig,l),'d=',detr(ig,l),'a=',alim(ig,l),'f=',fmc(ig,l),
1716 ! s 'f+1=',fmc(ig,l+1)
1717 IF (fmc(ig,l+1)<0.) THEN
1718 PRINT *, 'fmc1<0', l + 1, lmax(ig), fmc(ig, l+1)
1719 fmc(ig, l+1) = fmc(ig, l)
1720 detr(ig, l) = alim(ig, l) + entr(ig, l)
1721 ! fmc(ig,l+1)=0.
1722 ! print*,'fmc1<0',l+1,lmax(ig),fmc(ig,l+1)
1723 END IF
1724 ! if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1725 ! f_old=fmc(ig,l+1)
1726 ! fmc(ig,l+1)=fmc(ig,l)
1727 ! detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l+1)
1728 ! endif
1729
1730 ! if ((fmc(ig,l+1).gt.fmc(ig,l)).and.(l.gt.lentr(ig))) then
1731 ! f_old=fmc(ig,l+1)
1732 ! fmc(ig,l+1)=fmc(ig,l)
1733 ! detr(ig,l)=detr(ig,l)+f_old-fmc(ig,l)
1734 ! endif
1735 ! rajout du test sur alpha croissant
1736 ! if test
1737 ! if (1.eq.0) then
1738
1739 IF (l==klev) THEN
1740 PRINT *, 'THERMCELL PB ig=', ig, ' l=', l
1741 abort_message = 'THERMCELL PB'
1742 CALL abort_physic(modname, abort_message, 1)
1743 END IF
1744 ! if ((zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10).and.
1745 ! s (l.ge.lentr(ig)).and.
1746 IF ((zw2(ig,l+1)>1.E-10) .AND. (zw2(ig,l)>1.E-10) .AND. (l>=lentr(ig))) &
1747 THEN
1748 IF (((fmc(ig,l+1)/(rhobarz(ig,l+1)*zw2(ig,l+1)))>(fmc(ig,l)/ &
1749 (rhobarz(ig,l)*zw2(ig,l))))) THEN
1750 f_old = fmc(ig, l+1)
1751 fmc(ig, l+1) = fmc(ig, l)*rhobarz(ig, l+1)*zw2(ig, l+1)/ &
1752 (rhobarz(ig,l)*zw2(ig,l))
1753 detr(ig, l) = detr(ig, l) + f_old - fmc(ig, l+1)
1754 ! detr(ig,l)=(fmc(ig,l+1)-fmc(ig,l))/(0.4-1.)
1755 ! entr(ig,l)=0.4*detr(ig,l)
1756 ! entr(ig,l)=fmc(ig,l+1)-fmc(ig,l)+detr(ig,l)
1757 END IF
1758 END IF
1759 IF ((fmc(ig,l+1)>fmc(ig,l)) .AND. (l>lentr(ig))) THEN
1760 f_old = fmc(ig, l+1)
1761 fmc(ig, l+1) = fmc(ig, l)
1762 detr(ig, l) = detr(ig, l) + f_old - fmc(ig, l+1)
1763 END IF
1764 IF (detr(ig,l)>fmc(ig,l)) THEN
1765 detr(ig, l) = fmc(ig, l)
1766 entr(ig, l) = fmc(ig, l+1) - alim(ig, l)
1767 END IF
1768 IF (fmc(ig,l+1)<0.) THEN
1769 detr(ig, l) = detr(ig, l) + fmc(ig, l+1)
1770 fmc(ig, l+1) = 0.
1771 PRINT *, 'fmc2<0', l + 1, lmax(ig)
1772 END IF
1773
1774 ! test pour ne pas avoir f=0 et d=e/=0
1775 ! if (fmc(ig,l+1).lt.1.e-10) then
1776 ! detr(ig,l+1)=0.
1777 ! entr(ig,l+1)=0.
1778 ! zqla(ig,l+1)=0.
1779 ! zw2(ig,l+1)=0.
1780 ! lmax(ig)=l+1
1781 ! zmax(ig)=zlev(ig,lmax(ig))
1782 ! endif
1783 IF (zw2(ig,l+1)>1.E-10) THEN
1784 IF ((((fmc(ig,l+1))/(rhobarz(ig,l+1)*zw2(ig,l+1)))>1.)) THEN
1785 f_old = fmc(ig, l+1)
1786 fmc(ig, l+1) = rhobarz(ig, l+1)*zw2(ig, l+1)
1787 zw2(ig, l+1) = 0.
1788 zqla(ig, l+1) = 0.
1789 detr(ig, l) = detr(ig, l) + f_old - fmc(ig, l+1)
1790 lmax(ig) = l + 1
1791 zmax(ig) = zlev(ig, lmax(ig))
1792 PRINT *, 'alpha>1', l + 1, lmax(ig)
1793 END IF
1794 END IF
1795 ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
1796 ! endif test
1797 ! endif
1798 END DO
1799 END DO
1800 DO ig = 1, ngrid
1801 ! if (fmc(ig,lmax(ig)+1).ne.0.) then
1802 fmc(ig, lmax(ig)+1) = 0.
1803 entr(ig, lmax(ig)) = 0.
1804 detr(ig, lmax(ig)) = fmc(ig, lmax(ig)) + entr(ig, lmax(ig)) + &
1805 alim(ig, lmax(ig))
1806 ! endif
1807 END DO
1808 ! test sur le signe de fmc
1809 DO ig = 1, ngrid
1810 DO l = 1, klev + 1
1811 IF (fmc(ig,l)<0.) THEN
1812 PRINT *, 'fm1<0!!!', 'ig=', ig, 'l=', l, 'a=', alim(ig, l-1), 'e=', &
1813 entr(ig, l-1), 'f=', fmc(ig, l-1), 'd=', detr(ig, l-1), 'f+1=', &
1814 fmc(ig, l)
1815 END IF
1816 END DO
1817 END DO
1818 ! test de verification
1819 DO ig = 1, ngrid
1820 DO l = 1, lmax(ig)
1821 IF ((abs(fmc(ig,l+1)-fmc(ig,l)-alim(ig,l)-entr(ig,l)+ &
1822 detr(ig,l)))>1.E-4) THEN
1823 ! print*,'pbcm!!','ig=',ig,'l=',l,'lmax=',lmax(ig),'lmix=',lmix(ig),
1824 ! s 'e=',entr(ig,l),'d=',detr(ig,l),'a=',alim(ig,l),'f=',fmc(ig,l),
1825 ! s 'f+1=',fmc(ig,l+1)
1826 END IF
1827 IF (detr(ig,l)<0.) THEN
1828 PRINT *, 'detrdemi<0!!!'
1829 END IF
1830 END DO
1831 END DO
1832
1833 ! RC
1834 ! CR def de zmix continu (profil parabolique des vitesses)
1835 DO ig = 1, ngrid
1836 IF (lmix(ig)>1.) THEN
1837 ! test
1838 IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
1839 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
1840 zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- &
1841 (zlev(ig,lmix(ig)))))>1E-10) THEN
1842
1843 zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) &
1844 )**2-(zlev(ig,lmix(ig)+1))**2)-(zw2(ig,lmix(ig))-zw2(ig, &
1845 lmix(ig)+1))*((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))/ &
1846 (2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
1847 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
1848 zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
1849 ELSE
1850 zmix(ig) = zlev(ig, lmix(ig))
1851 PRINT *, 'pb zmix'
1852 END IF
1853 ELSE
1854 zmix(ig) = 0.
1855 END IF
1856 ! test
1857 IF ((zmax(ig)-zmix(ig))<=0.) THEN
1858 zmix(ig) = 0.9*zmax(ig)
1859 ! print*,'pb zmix>zmax'
1860 END IF
1861 END DO
1862 DO ig = 1, klon
1863 zmix0(ig) = zmix(ig)
1864 END DO
1865
1866 ! calcul du nouveau lmix correspondant
1867 DO ig = 1, ngrid
1868 DO l = 1, klev
1869 IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
1870 lmix(ig) = l
1871 END IF
1872 END DO
1873 END DO
1874
1875 ! ne devrait pas arriver!!!!!
1876 DO ig = 1, ngrid
1877 DO l = 1, klev
1878 IF (detr(ig,l)>(fmc(ig,l)+alim(ig,l))+entr(ig,l)) THEN
1879 PRINT *, 'detr2>fmc2!!!', 'ig=', ig, 'l=', l, 'd=', detr(ig, l), &
1880 'f=', fmc(ig, l), 'lmax=', lmax(ig)
1881 ! detr(ig,l)=fmc(ig,l)+alim(ig,l)+entr(ig,l)
1882 ! entr(ig,l)=0.
1883 ! fmc(ig,l+1)=0.
1884 ! zw2(ig,l+1)=0.
1885 ! zqla(ig,l+1)=0.
1886 PRINT *, 'pb!fm=0 et f_star>0', l, lmax(ig)
1887 ! lmax(ig)=l
1888 END IF
1889 END DO
1890 END DO
1891 DO ig = 1, ngrid
1892 DO l = lmax(ig) + 1, klev + 1
1893 ! fmc(ig,l)=0.
1894 ! detr(ig,l)=0.
1895 ! entr(ig,l)=0.
1896 ! zw2(ig,l)=0.
1897 ! zqla(ig,l)=0.
1898 END DO
1899 END DO
1900
1901 ! Calcul du detrainement lors du premier passage
1902 ! print*,'9 OK convect8'
1903 ! print*,'WA1 ',wa_moy
1904
1905 ! determination de l'indice du debut de la mixed layer ou w decroit
1906
1907 ! calcul de la largeur de chaque ascendance dans le cas conservatif.
1908 ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
1909 ! d'une couche est �gale � la hauteur de la couche alimentante.
1910 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
1911 ! de la vitesse d'entrainement horizontal dans la couche alimentante.
1912
1913 DO l = 2, nlay
1914 DO ig = 1, ngrid
1915 IF (l<=lmax(ig) .AND. (test(ig)==1)) THEN
1916 zw = max(wa_moy(ig,l), 1.E-10)
1917 larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw)
1918 END IF
1919 END DO
1920 END DO
1921
1922 DO l = 2, nlay
1923 DO ig = 1, ngrid
1924 IF (l<=lmax(ig) .AND. (test(ig)==1)) THEN
1925 ! if (idetr.eq.0) then
1926 ! cette option est finalement en dur.
1927 IF ((l_mix*zlev(ig,l))<0.) THEN
1928 PRINT *, 'pb l_mix*zlev<0'
1929 END IF
1930 ! CR: test: nouvelle def de lambda
1931 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
1932 IF (zw2(ig,l)>1.E-10) THEN
1933 larg_detr(ig, l) = sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
1934 ELSE
1935 larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
1936 END IF
1937 ! else if (idetr.eq.1) then
1938 ! larg_detr(ig,l)=larg_cons(ig,l)
1939 ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
1940 ! else if (idetr.eq.2) then
1941 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
1942 ! s *sqrt(wa_moy(ig,l))
1943 ! else if (idetr.eq.4) then
1944 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
1945 ! s *wa_moy(ig,l)
1946 ! endif
1947 END IF
1948 END DO
1949 END DO
1950
1951 ! print*,'10 OK convect8'
1952 ! print*,'WA2 ',wa_moy
1953 ! cal1cul de la fraction de la maille concern�e par l'ascendance en tenant
1954 ! compte de l'epluchage du thermique.
1955
1956
1957 DO l = 2, nlay
1958 DO ig = 1, ngrid
1959 IF (larg_cons(ig,l)>1. .AND. (test(ig)==1)) THEN
1960 ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
1961 fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
1962 ! test
1963 fraca(ig, l) = max(fraca(ig,l), 0.)
1964 fraca(ig, l) = min(fraca(ig,l), 0.5)
1965 fracd(ig, l) = 1. - fraca(ig, l)
1966 fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
1967 ELSE
1968 ! wa_moy(ig,l)=0.
1969 fraca(ig, l) = 0.
1970 fracc(ig, l) = 0.
1971 fracd(ig, l) = 1.
1972 END IF
1973 END DO
1974 END DO
1975 ! CR: calcul de fracazmix
1976 DO ig = 1, ngrid
1977 IF (test(ig)==1) THEN
1978 fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
1979 (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
1980 fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca( &
1981 ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
1982 END IF
1983 END DO
1984
1985 DO l = 2, nlay
1986 DO ig = 1, ngrid
1987 IF (larg_cons(ig,l)>1. .AND. (test(ig)==1)) THEN
1988 IF (l>lmix(ig)) THEN
1989 ! test
1990 IF (zmax(ig)-zmix(ig)<1.E-10) THEN
1991 ! print*,'pb xxx'
1992 xxx(ig, l) = (lmax(ig)+1.-l)/(lmax(ig)+1.-lmix(ig))
1993 ELSE
1994 xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
1995 END IF
1996 IF (idetr==0) THEN
1997 fraca(ig, l) = fracazmix(ig)
1998 ELSE IF (idetr==1) THEN
1999 fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
2000 ELSE IF (idetr==2) THEN
2001 fraca(ig, l) = fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
2002 ELSE
2003 fraca(ig, l) = fracazmix(ig)*xxx(ig, l)**2
2004 END IF
2005 ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
2006 fraca(ig, l) = max(fraca(ig,l), 0.)
2007 fraca(ig, l) = min(fraca(ig,l), 0.5)
2008 fracd(ig, l) = 1. - fraca(ig, l)
2009 fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
2010 END IF
2011 END IF
2012 END DO
2013 END DO
2014
2015 PRINT *, 'fin calcul fraca'
2016 ! print*,'11 OK convect8'
2017 ! print*,'Ea3 ',wa_moy
2018 ! ------------------------------------------------------------------
2019 ! Calcul de fracd, wd
2020 ! somme wa - wd = 0
2021 ! ------------------------------------------------------------------
2022
2023
2024 DO ig = 1, ngrid
2025 fm(ig, 1) = 0.
2026 fm(ig, nlay+1) = 0.
2027 END DO
2028
2029 DO l = 2, nlay
2030 DO ig = 1, ngrid
2031 IF (test(ig)==1) THEN
2032 fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
2033 ! CR:test
2034 IF (alim(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) &
2035 THEN
2036 fm(ig, l) = fm(ig, l-1)
2037 ! write(1,*)'ajustement fm, l',l
2038 END IF
2039 ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
2040 ! RC
2041 END IF
2042 END DO
2043 DO ig = 1, ngrid
2044 IF (fracd(ig,l)<0.1 .AND. (test(ig)==1)) THEN
2045 abort_message = 'fracd trop petit'
2046 CALL abort_physic(modname, abort_message, 1)
2047 ELSE
2048 ! vitesse descendante "diagnostique"
2049 wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
2050 END IF
2051 END DO
2052 END DO
2053
2054 DO l = 1, nlay + 1
2055 DO ig = 1, ngrid
2056 IF (test(ig)==0) THEN
2057 fm(ig, l) = fmc(ig, l)
2058 END IF
2059 END DO
2060 END DO
2061
2062 ! fin du first
2063 DO l = 1, nlay
2064 DO ig = 1, ngrid
2065 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
2066 masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
2067 END DO
2068 END DO
2069
2070 ! print*,'12 OK convect8'
2071 ! print*,'WA4 ',wa_moy
2072 ! c------------------------------------------------------------------
2073 ! calcul du transport vertical
2074 ! ------------------------------------------------------------------
2075
2076 GO TO 4444
2077 ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
2078 DO l = 2, nlay - 1
2079 DO ig = 1, ngrid
2080 IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( &
2081 ig,l+1)) THEN
2082 PRINT *, 'WARN!!! FM>M ig=', ig, ' l=', l, ' FM=', &
2083 fm(ig, l+1)*ptimestep, ' M=', masse(ig, l), masse(ig, l+1)
2084 END IF
2085 END DO
2086 END DO
2087
2088 DO l = 1, nlay
2089 DO ig = 1, ngrid
2090 IF ((alim(ig,l)+entr(ig,l))*ptimestep>masse(ig,l)) THEN
2091 PRINT *, 'WARN!!! E>M ig=', ig, ' l=', l, ' E==', &
2092 (entr(ig,l)+alim(ig,l))*ptimestep, ' M=', masse(ig, l)
2093 END IF
2094 END DO
2095 END DO
2096
2097 DO l = 1, nlay
2098 DO ig = 1, ngrid
2099 IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
2100 ! print*,'WARN!!! fm exagere ig=',ig,' l=',l
2101 ! s ,' FM=',fm(ig,l)
2102 END IF
2103 IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
2104 ! print*,'WARN!!! masse exagere ig=',ig,' l=',l
2105 ! s ,' M=',masse(ig,l)
2106 ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
2107 ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
2108 ! print*,'zlev(ig,l+1),zlev(ig,l)'
2109 ! s ,zlev(ig,l+1),zlev(ig,l)
2110 ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
2111 ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
2112 END IF
2113 IF (.NOT. alim(ig,l)>=0. .OR. .NOT. alim(ig,l)<=10.) THEN
2114 ! print*,'WARN!!! entr exagere ig=',ig,' l=',l
2115 ! s ,' E=',entr(ig,l)
2116 END IF
2117 END DO
2118 END DO
2119
2120 4444 CONTINUE
2121
2122 ! CR:redefinition du entr
2123 ! CR:test:on ne change pas la def du entr mais la def du fm
2124 DO l = 1, nlay
2125 DO ig = 1, ngrid
2126 IF (test(ig)==1) THEN
2127 detr(ig, l) = fm(ig, l) + alim(ig, l) - fm(ig, l+1)
2128 IF (detr(ig,l)<0.) THEN
2129 ! entr(ig,l)=entr(ig,l)-detr(ig,l)
2130 fm(ig, l+1) = fm(ig, l) + alim(ig, l)
2131 detr(ig, l) = 0.
2132 ! write(11,*)'l,ig,entr',l,ig,entr(ig,l)
2133 ! print*,'WARNING !!! detrainement negatif ',ig,l
2134 END IF
2135 END IF
2136 END DO
2137 END DO
2138 ! RC
2139
2140 IF (w2di==1) THEN
2141 fm0 = fm0 + ptimestep*(fm-fm0)/tho
2142 entr0 = entr0 + ptimestep*(alim+entr-entr0)/tho
2143 ELSE
2144 fm0 = fm
2145 entr0 = alim + entr
2146 detr0 = detr
2147 alim0 = alim
2148 ! zoa=zqta
2149 ! entr0=alim
2150 END IF
2151
2152 IF (1==1) THEN
2153 ! call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2154 ! . ,zh,zdhadj,zha)
2155 ! call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
2156 ! . ,zo,pdoadj,zoa)
2157 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zthl, &
2158 zdthladj, zta)
2159 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, po, pdoadj, &
2160 zoa)
2161 ELSE
2162 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
2163 zdhadj, zha)
2164 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
2165 pdoadj, zoa)
2166 END IF
2167
2168 IF (1==0) THEN
2169 CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
2170 zu, zv, pduadj, pdvadj, zua, zva)
2171 ELSE
2172 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
2173 zua)
2174 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
2175 zva)
2176 END IF
2177
2178 ! Calcul des moments
2179 ! do l=1,nlay
2180 ! do ig=1,ngrid
2181 ! zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
2182 ! zf2=zf/(1.-zf)
2183 ! thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
2184 ! wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
2185 ! enddo
2186 ! enddo
2187
2188
2189
2190
2191
2192
2193 ! print*,'13 OK convect8'
2194 ! print*,'WA5 ',wa_moy
2195 DO l = 1, nlay
2196 DO ig = 1, ngrid
2197 ! pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
2198 pdtadj(ig, l) = zdthladj(ig, l)*zpspsk(ig, l)
2199 END DO
2200 END DO
2201
2202
2203 ! do l=1,nlay
2204 ! do ig=1,ngrid
2205 ! if(abs(pdtadj(ig,l))*86400..gt.500.) then
2206 ! print*,'WARN!!! ig=',ig,' l=',l
2207 ! s ,' pdtadj=',pdtadj(ig,l)
2208 ! endif
2209 ! if(abs(pdoadj(ig,l))*86400..gt.1.) then
2210 ! print*,'WARN!!! ig=',ig,' l=',l
2211 ! s ,' pdoadj=',pdoadj(ig,l)
2212 ! endif
2213 ! enddo
2214 ! enddo
2215
2216 ! print*,'14 OK convect8'
2217 ! ------------------------------------------------------------------
2218 ! Calculs pour les sorties
2219 ! ------------------------------------------------------------------
2220 ! calcul de fraca pour les sorties
2221 DO l = 2, klev
2222 DO ig = 1, klon
2223 IF (zw2(ig,l)>1.E-10) THEN
2224 fraca(ig, l) = fm(ig, l)/(rhobarz(ig,l)*zw2(ig,l))
2225 ELSE
2226 fraca(ig, l) = 0.
2227 END IF
2228 END DO
2229 END DO
2230 IF (sorties) THEN
2231 DO l = 1, nlay
2232 DO ig = 1, ngrid
2233 zla(ig, l) = (1.-fracd(ig,l))*zmax(ig)
2234 zld(ig, l) = fracd(ig, l)*zmax(ig)
2235 IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ &
2236 (1.-fracd(ig,l))
2237 END DO
2238 END DO
2239 ! CR calcul du niveau de condensation
2240 ! initialisation
2241 DO ig = 1, ngrid
2242 nivcon(ig) = 0.
2243 zcon(ig) = 0.
2244 END DO
2245 DO k = nlay, 1, -1
2246 DO ig = 1, ngrid
2247 IF (zqla(ig,k)>1E-10) THEN
2248 nivcon(ig) = k
2249 zcon(ig) = zlev(ig, k)
2250 END IF
2251 ! if (zcon(ig).gt.1.e-10) then
2252 ! nuage=.true.
2253 ! else
2254 ! nuage=.false.
2255 ! endif
2256 END DO
2257 END DO
2258
2259 DO l = 1, nlay
2260 DO ig = 1, ngrid
2261 zf = fraca(ig, l)
2262 zf2 = zf/(1.-zf)
2263 thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
2264 wth2(ig, l) = zf2*(zw2(ig,l))**2
2265 ! print*,'wth2=',wth2(ig,l)
2266 wth3(ig, l) = zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))*zw2(ig, l)* &
2267 zw2(ig, l)*zw2(ig, l)
2268 q2(ig, l) = zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
2269 ! test: on calcul q2/po=ratqsc
2270 ! if (nuage) then
2271 ratqscth(ig, l) = sqrt(q2(ig,l))/(po(ig,l)*1000.)
2272 ! else
2273 ! ratqscth(ig,l)=0.
2274 ! endif
2275 END DO
2276 END DO
2277 ! calcul du ratqscdiff
2278 sum = 0.
2279 sumdiff = 0.
2280 ratqsdiff(:, :) = 0.
2281 DO ig = 1, ngrid
2282 DO l = 1, lentr(ig)
2283 sum = sum + alim_star(ig, l)*zqta(ig, l)*1000.
2284 END DO
2285 END DO
2286 DO ig = 1, ngrid
2287 DO l = 1, lentr(ig)
2288 zf = fraca(ig, l)
2289 zf2 = zf/(1.-zf)
2290 sumdiff = sumdiff + alim_star(ig, l)*(zqta(ig,l)*1000.-sum)**2
2291 ! ratqsdiff=ratqsdiff+alim_star(ig,l)*
2292 ! s (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
2293 END DO
2294 END DO
2295 DO l = 1, klev
2296 DO ig = 1, ngrid
2297 ratqsdiff(ig, l) = sqrt(sumdiff)/(po(ig,l)*1000.)
2298 ! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
2299 END DO
2300 END DO
2301
2302 END IF
2303
2304 ! print*,'19 OK convect8'
2305 RETURN
2306 END SUBROUTINE thermcell_cld
2307
2308 SUBROUTINE thermcell_eau(ngrid, nlay, ptimestep, pplay, pplev, pphi, pu, pv, &
2309 pt, po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0 & ! s
2310 ! ,pu_therm,pv_therm
2311 , r_aspect, l_mix, w2di, tho)
2312
2313 USE dimphy
2314 IMPLICIT NONE
2315
2316 ! =======================================================================
2317
2318 ! Calcul du transport verticale dans la couche limite en presence
2319 ! de "thermiques" explicitement representes
2320
2321 ! R��criture � partir d'un listing papier � Habas, le 14/02/00
2322
2323 ! le thermique est suppos� homog�ne et dissip� par m�lange avec
2324 ! son environnement. la longueur l_mix contr�le l'efficacit� du
2325 ! m�lange
2326
2327 ! Le calcul du transport des diff�rentes esp�ces se fait en prenant
2328 ! en compte:
2329 ! 1. un flux de masse montant
2330 ! 2. un flux de masse descendant
2331 ! 3. un entrainement
2332 ! 4. un detrainement
2333
2334 ! =======================================================================
2335
2336 ! -----------------------------------------------------------------------
2337 ! declarations:
2338 ! -------------
2339
2340 include "YOMCST.h"
2341 include "YOETHF.h"
2342 include "FCTTRE.h"
2343
2344 ! arguments:
2345 ! ----------
2346
2347 INTEGER ngrid, nlay, w2di
2348 REAL tho
2349 REAL ptimestep, l_mix, r_aspect
2350 REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
2351 REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
2352 REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
2353 REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
2354 REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
2355 REAL pphi(ngrid, nlay)
2356
2357 INTEGER idetr
2358 SAVE idetr
2359 DATA idetr/3/
2360 !$OMP THREADPRIVATE(idetr)
2361
2362 ! local:
2363 ! ------
2364
2365 INTEGER ig, k, l, lmaxa(klon), lmix(klon)
2366 REAL zsortie1d(klon)
2367 ! CR: on remplace lmax(klon,klev+1)
2368 INTEGER lmax(klon), lmin(klon), lentr(klon)
2369 REAL linter(klon)
2370 REAL zmix(klon), fracazmix(klon)
2371 ! RC
2372 REAL zmax(klon), zw, zz, zw2(klon, klev+1), ztva(klon, klev), zzz
2373
2374 REAL zlev(klon, klev+1), zlay(klon, klev)
2375 REAL zh(klon, klev), zdhadj(klon, klev)
2376 REAL zthl(klon, klev), zdthladj(klon, klev)
2377 REAL ztv(klon, klev)
2378 REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
2379 REAL zl(klon, klev)
2380 REAL wh(klon, klev+1)
2381 REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
2382 REAL zla(klon, klev+1)
2383 REAL zwa(klon, klev+1)
2384 REAL zld(klon, klev+1)
2385 REAL zwd(klon, klev+1)
2386 REAL zsortie(klon, klev)
2387 REAL zva(klon, klev)
2388 REAL zua(klon, klev)
2389 REAL zoa(klon, klev)
2390
2391 REAL zta(klon, klev)
2392 REAL zha(klon, klev)
2393 REAL wa_moy(klon, klev+1)
2394 REAL fraca(klon, klev+1)
2395 REAL fracc(klon, klev+1)
2396 REAL zf, zf2
2397 REAL thetath2(klon, klev), wth2(klon, klev)
2398 ! common/comtherm/thetath2,wth2
2399
2400 REAL count_time
2401 INTEGER ialt
2402
2403 LOGICAL sorties
2404 REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
2405 REAL zpspsk(klon, klev)
2406
2407 ! real wmax(klon,klev),wmaxa(klon)
2408 REAL wmax(klon), wmaxa(klon)
2409 REAL wa(klon, klev, klev+1)
2410 REAL wd(klon, klev+1)
2411 REAL larg_part(klon, klev, klev+1)
2412 REAL fracd(klon, klev+1)
2413 REAL xxx(klon, klev+1)
2414 REAL larg_cons(klon, klev+1)
2415 REAL larg_detr(klon, klev+1)
2416 REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
2417 REAL pu_therm(klon, klev), pv_therm(klon, klev)
2418 REAL fm(klon, klev+1), entr(klon, klev)
2419 REAL fmc(klon, klev+1)
2420
2421 REAL zcor, zdelta, zcvm5, qlbef
2422 REAL tbef(klon), qsatbef(klon)
2423 REAL dqsat_dt, dt, num, denom
2424 REAL reps, rlvcp, ddt0
2425 REAL ztla(klon, klev), zqla(klon, klev), zqta(klon, klev)
2426
2427 PARAMETER (ddt0=.01)
2428
2429 ! CR:nouvelles variables
2430 REAL f_star(klon, klev+1), entr_star(klon, klev)
2431 REAL entr_star_tot(klon), entr_star2(klon)
2432 REAL f(klon), f0(klon)
2433 REAL zlevinter(klon)
2434 LOGICAL first
2435 DATA first/.FALSE./
2436 SAVE first
2437 !$OMP THREADPRIVATE(first)
2438
2439 ! RC
2440
2441 CHARACTER *2 str2
2442 CHARACTER *10 str10
2443
2444 CHARACTER (LEN=20) :: modname = 'thermcell_eau'
2445 CHARACTER (LEN=80) :: abort_message
2446
2447 LOGICAL vtest(klon), down
2448 LOGICAL zsat(klon)
2449
2450 EXTERNAL scopy
2451
2452 INTEGER ncorrec, ll
2453 SAVE ncorrec
2454 DATA ncorrec/0/
2455 !$OMP THREADPRIVATE(ncorrec)
2456
2457
2458
2459 ! -----------------------------------------------------------------------
2460 ! initialisation:
2461 ! ---------------
2462
2463 sorties = .TRUE.
2464 IF (ngrid/=klon) THEN
2465 PRINT *
2466 PRINT *, 'STOP dans convadj'
2467 PRINT *, 'ngrid =', ngrid
2468 PRINT *, 'klon =', klon
2469 END IF
2470
2471 ! Initialisation
2472 rlvcp = rlvtt/rcpd
2473 reps = rd/rv
2474
2475 ! -----------------------------------------------------------------------
2476 ! AM Calcul de T,q,ql a partir de Tl et qT
2477 ! ---------------------------------------------------
2478
2479 ! Pr Tprec=Tl calcul de qsat
2480 ! Si qsat>qT T=Tl, q=qT
2481 ! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt)
2482 ! On cherche DDT < DDT0
2483
2484 ! defaut
2485 DO ll = 1, nlay
2486 DO ig = 1, ngrid
2487 zo(ig, ll) = po(ig, ll)
2488 zl(ig, ll) = 0.
2489 zh(ig, ll) = pt(ig, ll)
2490 END DO
2491 END DO
2492 DO ig = 1, ngrid
2493 zsat(ig) = .FALSE.
2494 END DO
2495
2496
2497 DO ll = 1, nlay
2498 ! les points insatures sont definitifs
2499 DO ig = 1, ngrid
2500 tbef(ig) = pt(ig, ll)
2501 zdelta = max(0., sign(1.,rtt-tbef(ig)))
2502 qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, ll)
2503 qsatbef(ig) = min(0.5, qsatbef(ig))
2504 zcor = 1./(1.-retv*qsatbef(ig))
2505 qsatbef(ig) = qsatbef(ig)*zcor
2506 zsat(ig) = (max(0.,po(ig,ll)-qsatbef(ig))>0.00001)
2507 END DO
2508
2509 DO ig = 1, ngrid
2510 IF (zsat(ig)) THEN
2511 qlbef = max(0., po(ig,ll)-qsatbef(ig))
2512 ! si sature: ql est surestime, d'ou la sous-relax
2513 dt = 0.5*rlvcp*qlbef
2514 ! on pourra enchainer 2 ou 3 calculs sans Do while
2515 DO WHILE (dt>ddt0)
2516 ! il faut verifier si c,a conserve quand on repasse en insature ...
2517 tbef(ig) = tbef(ig) + dt
2518 zdelta = max(0., sign(1.,rtt-tbef(ig)))
2519 qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, ll)
2520 qsatbef(ig) = min(0.5, qsatbef(ig))
2521 zcor = 1./(1.-retv*qsatbef(ig))
2522 qsatbef(ig) = qsatbef(ig)*zcor
2523 ! on veut le signe de qlbef
2524 qlbef = po(ig, ll) - qsatbef(ig)
2525 ! dqsat_dT
2526 zdelta = max(0., sign(1.,rtt-tbef(ig)))
2527 zcvm5 = r5les*(1.-zdelta) + r5ies*zdelta
2528 zcor = 1./(1.-retv*qsatbef(ig))
2529 dqsat_dt = foede(tbef(ig), zdelta, zcvm5, qsatbef(ig), zcor)
2530 num = -tbef(ig) + pt(ig, ll) + rlvcp*qlbef
2531 denom = 1. + rlvcp*dqsat_dt
2532 dt = num/denom
2533 END DO
2534 ! on ecrit de maniere conservative (sat ou non)
2535 zl(ig, ll) = max(0., qlbef)
2536 ! T = Tl +Lv/Cp ql
2537 zh(ig, ll) = pt(ig, ll) + rlvcp*zl(ig, ll)
2538 zo(ig, ll) = po(ig, ll) - zl(ig, ll)
2539 END IF
2540 END DO
2541 END DO
2542 ! AM fin
2543
2544 ! -----------------------------------------------------------------------
2545 ! incrementation eventuelle de tendances precedentes:
2546 ! ---------------------------------------------------
2547
2548 ! print*,'0 OK convect8'
2549
2550 DO l = 1, nlay
2551 DO ig = 1, ngrid
2552 zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa
2553 ! zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
2554 zu(ig, l) = pu(ig, l)
2555 zv(ig, l) = pv(ig, l)
2556 ! zo(ig,l)=po(ig,l)
2557 ! ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
2558 ! AM attention zh est maintenant le profil de T et plus le profil de
2559 ! theta !
2560
2561 ! T-> Theta
2562 ztv(ig, l) = zh(ig, l)/zpspsk(ig, l)
2563 ! AM Theta_v
2564 ztv(ig, l) = ztv(ig, l)*(1.+retv*(zo(ig,l))-zl(ig,l))
2565 ! AM Thetal
2566 zthl(ig, l) = pt(ig, l)/zpspsk(ig, l)
2567
2568 END DO
2569 END DO
2570
2571 ! print*,'1 OK convect8'
2572 ! --------------------
2573
2574
2575 ! + + + + + + + + + + +
2576
2577
2578 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
2579 ! wh,wt,wo ...
2580
2581 ! + + + + + + + + + + + zh,zu,zv,zo,rho
2582
2583
2584 ! -------------------- zlev(1)
2585 ! \\\\\\\\\\\\\\\\\\\
2586
2587
2588 ! -----------------------------------------------------------------------
2589 ! Calcul des altitudes des couches
2590 ! -----------------------------------------------------------------------
2591
2592 DO l = 2, nlay
2593 DO ig = 1, ngrid
2594 zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
2595 END DO
2596 END DO
2597 DO ig = 1, ngrid
2598 zlev(ig, 1) = 0.
2599 zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
2600 END DO
2601 DO l = 1, nlay
2602 DO ig = 1, ngrid
2603 zlay(ig, l) = pphi(ig, l)/rg
2604 END DO
2605 END DO
2606
2607 ! print*,'2 OK convect8'
2608 ! -----------------------------------------------------------------------
2609 ! Calcul des densites
2610 ! -----------------------------------------------------------------------
2611
2612 DO l = 1, nlay
2613 DO ig = 1, ngrid
2614 ! rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
2615 rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*ztv(ig,l))
2616 END DO
2617 END DO
2618
2619 DO l = 2, nlay
2620 DO ig = 1, ngrid
2621 rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
2622 END DO
2623 END DO
2624
2625 DO k = 1, nlay
2626 DO l = 1, nlay + 1
2627 DO ig = 1, ngrid
2628 wa(ig, k, l) = 0.
2629 END DO
2630 END DO
2631 END DO
2632
2633 ! print*,'3 OK convect8'
2634 ! ------------------------------------------------------------------
2635 ! Calcul de w2, quarre de w a partir de la cape
2636 ! a partir de w2, on calcule wa, vitesse de l'ascendance
2637
2638 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
2639 ! w2 est stoke dans wa
2640
2641 ! ATTENTION: dans convect8, on n'utilise le calcule des wa
2642 ! independants par couches que pour calculer l'entrainement
2643 ! a la base et la hauteur max de l'ascendance.
2644
2645 ! Indicages:
2646 ! l'ascendance provenant du niveau k traverse l'interface l avec
2647 ! une vitesse wa(k,l).
2648
2649 ! --------------------
2650
2651 ! + + + + + + + + + +
2652
2653 ! wa(k,l) ---- -------------------- l
2654 ! / ! /||\ + + + + + + + + + +
2655 ! ||
2656 ! || --------------------
2657 ! ||
2658 ! || + + + + + + + + + +
2659 ! ||
2660 ! || --------------------
2661 ! ||__
2662 ! |___ + + + + + + + + + + k
2663
2664 ! --------------------
2665
2666
2667
2668 ! ------------------------------------------------------------------
2669
2670 ! CR: ponderation entrainement des couches instables
2671 ! def des entr_star tels que entr=f*entr_star
2672 DO l = 1, klev
2673 DO ig = 1, ngrid
2674 entr_star(ig, l) = 0.
2675 END DO
2676 END DO
2677 ! determination de la longueur de la couche d entrainement
2678 DO ig = 1, ngrid
2679 lentr(ig) = 1
2680 END DO
2681
2682 ! on ne considere que les premieres couches instables
2683 DO k = nlay - 1, 1, -1
2684 DO ig = 1, ngrid
2685 IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<ztv(ig,k+2)) THEN
2686 lentr(ig) = k
2687 END IF
2688 END DO
2689 END DO
2690
2691 ! determination du lmin: couche d ou provient le thermique
2692 DO ig = 1, ngrid
2693 lmin(ig) = 1
2694 END DO
2695 DO ig = 1, ngrid
2696 DO l = nlay, 2, -1
2697 IF (ztv(ig,l-1)>ztv(ig,l)) THEN
2698 lmin(ig) = l - 1
2699 END IF
2700 END DO
2701 END DO
2702
2703 ! definition de l'entrainement des couches
2704 DO l = 1, klev - 1
2705 DO ig = 1, ngrid
2706 IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<=lentr(ig)) THEN
2707 entr_star(ig, l) = (ztv(ig,l)-ztv(ig,l+1))*(zlev(ig,l+1)