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)-zlev(ig,l))
2708 END IF
2709 END DO
2710 END DO
2711 ! pas de thermique si couche 1 stable
2712 DO ig = 1, ngrid
2713 IF (lmin(ig)>1) THEN
2714 DO l = 1, klev
2715 entr_star(ig, l) = 0.
2716 END DO
2717 END IF
2718 END DO
2719 ! calcul de l entrainement total
2720 DO ig = 1, ngrid
2721 entr_star_tot(ig) = 0.
2722 END DO
2723 DO ig = 1, ngrid
2724 DO k = 1, klev
2725 entr_star_tot(ig) = entr_star_tot(ig) + entr_star(ig, k)
2726 END DO
2727 END DO
2728
2729 DO k = 1, klev
2730 DO ig = 1, ngrid
2731 ztva(ig, k) = ztv(ig, k)
2732 END DO
2733 END DO
2734 ! RC
2735 ! AM:initialisations
2736 DO k = 1, nlay
2737 DO ig = 1, ngrid
2738 ztva(ig, k) = ztv(ig, k)
2739 ztla(ig, k) = zthl(ig, k)
2740 zqla(ig, k) = 0.
2741 zqta(ig, k) = po(ig, k)
2742 zsat(ig) = .FALSE.
2743 END DO
2744 END DO
2745
2746 ! print*,'7 OK convect8'
2747 DO k = 1, klev + 1
2748 DO ig = 1, ngrid
2749 zw2(ig, k) = 0.
2750 fmc(ig, k) = 0.
2751 ! CR
2752 f_star(ig, k) = 0.
2753 ! RC
2754 larg_cons(ig, k) = 0.
2755 larg_detr(ig, k) = 0.
2756 wa_moy(ig, k) = 0.
2757 END DO
2758 END DO
2759
2760 ! print*,'8 OK convect8'
2761 DO ig = 1, ngrid
2762 linter(ig) = 1.
2763 lmaxa(ig) = 1
2764 lmix(ig) = 1
2765 wmaxa(ig) = 0.
2766 END DO
2767
2768 ! CR:
2769 DO l = 1, nlay - 2
2770 DO ig = 1, ngrid
2771 IF (ztv(ig,l)>ztv(ig,l+1) .AND. entr_star(ig,l)>1.E-10 .AND. &
2772 zw2(ig,l)<1E-10) THEN
2773 ! AM
2774 ztla(ig, l) = zthl(ig, l)
2775 zqta(ig, l) = po(ig, l)
2776 zqla(ig, l) = zl(ig, l)
2777 ! AM
2778 f_star(ig, l+1) = entr_star(ig, l)
2779 ! test:calcul de dteta
2780 zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
2781 (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
2782 larg_detr(ig, l) = 0.
2783 ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+entr_star(ig, &
2784 l)>1.E-10)) THEN
2785 f_star(ig, l+1) = f_star(ig, l) + entr_star(ig, l)
2786
2787 ! AM on melange Tl et qt du thermique
2788 ztla(ig, l) = (f_star(ig,l)*ztla(ig,l-1)+entr_star(ig,l)*zthl(ig,l))/ &
2789 f_star(ig, l+1)
2790 zqta(ig, l) = (f_star(ig,l)*zqta(ig,l-1)+entr_star(ig,l)*po(ig,l))/ &
2791 f_star(ig, l+1)
2792
2793 ! ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
2794 ! s *ztv(ig,l))/f_star(ig,l+1)
2795
2796 ! AM on en deduit thetav et ql du thermique
2797 tbef(ig) = ztla(ig, l)*zpspsk(ig, l)
2798 zdelta = max(0., sign(1.,rtt-tbef(ig)))
2799 qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, l)
2800 qsatbef(ig) = min(0.5, qsatbef(ig))
2801 zcor = 1./(1.-retv*qsatbef(ig))
2802 qsatbef(ig) = qsatbef(ig)*zcor
2803 zsat(ig) = (max(0.,zqta(ig,l)-qsatbef(ig))>0.00001)
2804 END IF
2805 END DO
2806 DO ig = 1, ngrid
2807 IF (zsat(ig)) THEN
2808 qlbef = max(0., zqta(ig,l)-qsatbef(ig))
2809 dt = 0.5*rlvcp*qlbef
2810 DO WHILE (dt>ddt0)
2811 tbef(ig) = tbef(ig) + dt
2812 zdelta = max(0., sign(1.,rtt-tbef(ig)))
2813 qsatbef(ig) = r2es*foeew(tbef(ig), zdelta)/pplev(ig, l)
2814 qsatbef(ig) = min(0.5, qsatbef(ig))
2815 zcor = 1./(1.-retv*qsatbef(ig))
2816 qsatbef(ig) = qsatbef(ig)*zcor
2817 qlbef = zqta(ig, l) - qsatbef(ig)
2818
2819 zdelta = max(0., sign(1.,rtt-tbef(ig)))
2820 zcvm5 = r5les*(1.-zdelta) + r5ies*zdelta
2821 zcor = 1./(1.-retv*qsatbef(ig))
2822 dqsat_dt = foede(tbef(ig), zdelta, zcvm5, qsatbef(ig), zcor)
2823 num = -tbef(ig) + ztla(ig, l)*zpspsk(ig, l) + rlvcp*qlbef
2824 denom = 1. + rlvcp*dqsat_dt
2825 dt = num/denom
2826 END DO
2827 zqla(ig, l) = max(0., zqta(ig,l)-qsatbef(ig))
2828 END IF
2829 ! on ecrit de maniere conservative (sat ou non)
2830 ! T = Tl +Lv/Cp ql
2831 ztva(ig, l) = ztla(ig, l)*zpspsk(ig, l) + rlvcp*zqla(ig, l)
2832 ztva(ig, l) = ztva(ig, l)/zpspsk(ig, l)
2833 ztva(ig, l) = ztva(ig, l)*(1.+retv*(zqta(ig,l)-zqla(ig,l))-zqla(ig,l))
2834
2835 END DO
2836 DO ig = 1, ngrid
2837 IF (zw2(ig,l)>=1.E-10 .AND. f_star(ig,l)+entr_star(ig,l)>1.E-10) THEN
2838 ! mise a jour de la vitesse ascendante (l'air entraine de la couche
2839 ! consideree commence avec une vitesse nulle).
2840
2841 zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/f_star(ig,l+1))**2 + &
2842 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
2843 END IF
2844 ! determination de zmax continu par interpolation lineaire
2845 IF (zw2(ig,l+1)<0.) THEN
2846 linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
2847 ig,l))
2848 zw2(ig, l+1) = 0.
2849 lmaxa(ig) = l
2850 ELSE
2851 wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
2852 END IF
2853 IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
2854 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
2855 lmix(ig) = l + 1
2856 wmaxa(ig) = wa_moy(ig, l+1)
2857 END IF
2858 END DO
2859 END DO
2860
2861 ! Calcul de la couche correspondant a la hauteur du thermique
2862 DO ig = 1, ngrid
2863 lmax(ig) = lentr(ig)
2864 END DO
2865 DO ig = 1, ngrid
2866 DO l = nlay, lentr(ig) + 1, -1
2867 IF (zw2(ig,l)<=1.E-10) THEN
2868 lmax(ig) = l - 1
2869 END IF
2870 END DO
2871 END DO
2872 ! pas de thermique si couche 1 stable
2873 DO ig = 1, ngrid
2874 IF (lmin(ig)>1) THEN
2875 lmax(ig) = 1
2876 lmin(ig) = 1
2877 END IF
2878 END DO
2879
2880 ! Determination de zw2 max
2881 DO ig = 1, ngrid
2882 wmax(ig) = 0.
2883 END DO
2884
2885 DO l = 1, nlay
2886 DO ig = 1, ngrid
2887 IF (l<=lmax(ig)) THEN
2888 zw2(ig, l) = sqrt(zw2(ig,l))
2889 wmax(ig) = max(wmax(ig), zw2(ig,l))
2890 ELSE
2891 zw2(ig, l) = 0.
2892 END IF
2893 END DO
2894 END DO
2895
2896 ! Longueur caracteristique correspondant a la hauteur des thermiques.
2897 DO ig = 1, ngrid
2898 zmax(ig) = 500.
2899 zlevinter(ig) = zlev(ig, 1)
2900 END DO
2901 DO ig = 1, ngrid
2902 ! calcul de zlevinter
2903 zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
2904 zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
2905 zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,lmin(ig)))
2906 END DO
2907
2908 ! Fermeture,determination de f
2909 DO ig = 1, ngrid
2910 entr_star2(ig) = 0.
2911 END DO
2912 DO ig = 1, ngrid
2913 IF (entr_star_tot(ig)<1.E-10) THEN
2914 f(ig) = 0.
2915 ELSE
2916 DO k = lmin(ig), lentr(ig)
2917 entr_star2(ig) = entr_star2(ig) + entr_star(ig, k)**2/(rho(ig,k)*( &
2918 zlev(ig,k+1)-zlev(ig,k)))
2919 END DO
2920 ! Nouvelle fermeture
2921 f(ig) = wmax(ig)/(zmax(ig)*r_aspect*entr_star2(ig))*entr_star_tot(ig)
2922 ! test
2923 IF (first) THEN
2924 f(ig) = f(ig) + (f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)*wmax(ig))
2925 END IF
2926 END IF
2927 f0(ig) = f(ig)
2928 first = .TRUE.
2929 END DO
2930
2931 ! Calcul de l'entrainement
2932 DO k = 1, klev
2933 DO ig = 1, ngrid
2934 entr(ig, k) = f(ig)*entr_star(ig, k)
2935 END DO
2936 END DO
2937 ! Calcul des flux
2938 DO ig = 1, ngrid
2939 DO l = 1, lmax(ig) - 1
2940 fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
2941 END DO
2942 END DO
2943
2944 ! RC
2945
2946
2947 ! print*,'9 OK convect8'
2948 ! print*,'WA1 ',wa_moy
2949
2950 ! determination de l'indice du debut de la mixed layer ou w decroit
2951
2952 ! calcul de la largeur de chaque ascendance dans le cas conservatif.
2953 ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
2954 ! d'une couche est �gale � la hauteur de la couche alimentante.
2955 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
2956 ! de la vitesse d'entrainement horizontal dans la couche alimentante.
2957
2958 DO l = 2, nlay
2959 DO ig = 1, ngrid
2960 IF (l<=lmaxa(ig)) THEN
2961 zw = max(wa_moy(ig,l), 1.E-10)
2962 larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw)
2963 END IF
2964 END DO
2965 END DO
2966
2967 DO l = 2, nlay
2968 DO ig = 1, ngrid
2969 IF (l<=lmaxa(ig)) THEN
2970 ! if (idetr.eq.0) then
2971 ! cette option est finalement en dur.
2972 larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
2973 ! else if (idetr.eq.1) then
2974 ! larg_detr(ig,l)=larg_cons(ig,l)
2975 ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
2976 ! else if (idetr.eq.2) then
2977 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2978 ! s *sqrt(wa_moy(ig,l))
2979 ! else if (idetr.eq.4) then
2980 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
2981 ! s *wa_moy(ig,l)
2982 ! endif
2983 END IF
2984 END DO
2985 END DO
2986
2987 ! print*,'10 OK convect8'
2988 ! print*,'WA2 ',wa_moy
2989 ! calcul de la fraction de la maille concern�e par l'ascendance en tenant
2990 ! compte de l'epluchage du thermique.
2991
2992 ! CR def de zmix continu (profil parabolique des vitesses)
2993 DO ig = 1, ngrid
2994 IF (lmix(ig)>1.) THEN
2995 zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig))) &
2996 **2-(zlev(ig,lmix(ig)+1))**2)-(zw2(ig,lmix(ig))-zw2(ig, &
2997 lmix(ig)+1))*((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))/ &
2998 (2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
2999 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))*((zlev( &
3000 ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
3001 ELSE
3002 zmix(ig) = 0.
3003 END IF
3004 END DO
3005
3006 ! calcul du nouveau lmix correspondant
3007 DO ig = 1, ngrid
3008 DO l = 1, klev
3009 IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
3010 lmix(ig) = l
3011 END IF
3012 END DO
3013 END DO
3014
3015 DO l = 2, nlay
3016 DO ig = 1, ngrid
3017 IF (larg_cons(ig,l)>1.) THEN
3018 ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
3019 fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
3020 ! test
3021 fraca(ig, l) = max(fraca(ig,l), 0.)
3022 fraca(ig, l) = min(fraca(ig,l), 0.5)
3023 fracd(ig, l) = 1. - fraca(ig, l)
3024 fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
3025 ELSE
3026 ! wa_moy(ig,l)=0.
3027 fraca(ig, l) = 0.
3028 fracc(ig, l) = 0.
3029 fracd(ig, l) = 1.
3030 END IF
3031 END DO
3032 END DO
3033 ! CR: calcul de fracazmix
3034 DO ig = 1, ngrid
3035 fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
3036 (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
3037 fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca(ig &
3038 ,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
3039 END DO
3040
3041 DO l = 2, nlay
3042 DO ig = 1, ngrid
3043 IF (larg_cons(ig,l)>1.) THEN
3044 IF (l>lmix(ig)) THEN
3045 xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
3046 IF (idetr==0) THEN
3047 fraca(ig, l) = fracazmix(ig)
3048 ELSE IF (idetr==1) THEN
3049 fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
3050 ELSE IF (idetr==2) THEN
3051 fraca(ig, l) = fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
3052 ELSE
3053 fraca(ig, l) = fracazmix(ig)*xxx(ig, l)**2
3054 END IF
3055 ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
3056 fraca(ig, l) = max(fraca(ig,l), 0.)
3057 fraca(ig, l) = min(fraca(ig,l), 0.5)
3058 fracd(ig, l) = 1. - fraca(ig, l)
3059 fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
3060 END IF
3061 END IF
3062 END DO
3063 END DO
3064
3065 ! print*,'11 OK convect8'
3066 ! print*,'Ea3 ',wa_moy
3067 ! ------------------------------------------------------------------
3068 ! Calcul de fracd, wd
3069 ! somme wa - wd = 0
3070 ! ------------------------------------------------------------------
3071
3072
3073 DO ig = 1, ngrid
3074 fm(ig, 1) = 0.
3075 fm(ig, nlay+1) = 0.
3076 END DO
3077
3078 DO l = 2, nlay
3079 DO ig = 1, ngrid
3080 fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
3081 ! CR:test
3082 IF (entr(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) THEN
3083 fm(ig, l) = fm(ig, l-1)
3084 ! write(1,*)'ajustement fm, l',l
3085 END IF
3086 ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
3087 ! RC
3088 END DO
3089 DO ig = 1, ngrid
3090 IF (fracd(ig,l)<0.1) THEN
3091 abort_message = 'fracd trop petit'
3092 CALL abort_physic(modname, abort_message, 1)
3093 ELSE
3094 ! vitesse descendante "diagnostique"
3095 wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
3096 END IF
3097 END DO
3098 END DO
3099
3100 DO l = 1, nlay
3101 DO ig = 1, ngrid
3102 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
3103 masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
3104 END DO
3105 END DO
3106
3107 ! print*,'12 OK convect8'
3108 ! print*,'WA4 ',wa_moy
3109 ! c------------------------------------------------------------------
3110 ! calcul du transport vertical
3111 ! ------------------------------------------------------------------
3112
3113 GO TO 4444
3114 ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
3115 DO l = 2, nlay - 1
3116 DO ig = 1, ngrid
3117 IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( &
3118 ig,l+1)) THEN
3119 ! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
3120 ! s ,fm(ig,l+1)*ptimestep
3121 ! s ,' M=',masse(ig,l),masse(ig,l+1)
3122 END IF
3123 END DO
3124 END DO
3125
3126 DO l = 1, nlay
3127 DO ig = 1, ngrid
3128 IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
3129 ! print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
3130 ! s ,entr(ig,l)*ptimestep
3131 ! s ,' M=',masse(ig,l)
3132 END IF
3133 END DO
3134 END DO
3135
3136 DO l = 1, nlay
3137 DO ig = 1, ngrid
3138 IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
3139 ! print*,'WARN!!! fm exagere ig=',ig,' l=',l
3140 ! s ,' FM=',fm(ig,l)
3141 END IF
3142 IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
3143 ! print*,'WARN!!! masse exagere ig=',ig,' l=',l
3144 ! s ,' M=',masse(ig,l)
3145 ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
3146 ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
3147 ! print*,'zlev(ig,l+1),zlev(ig,l)'
3148 ! s ,zlev(ig,l+1),zlev(ig,l)
3149 ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
3150 ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
3151 END IF
3152 IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
3153 ! print*,'WARN!!! entr exagere ig=',ig,' l=',l
3154 ! s ,' E=',entr(ig,l)
3155 END IF
3156 END DO
3157 END DO
3158
3159 4444 CONTINUE
3160
3161 IF (w2di==1) THEN
3162 fm0 = fm0 + ptimestep*(fm-fm0)/tho
3163 entr0 = entr0 + ptimestep*(entr-entr0)/tho
3164 ELSE
3165 fm0 = fm
3166 entr0 = entr
3167 END IF
3168
3169 IF (1==1) THEN
3170 ! call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3171 ! . ,zh,zdhadj,zha)
3172 ! call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
3173 ! . ,zo,pdoadj,zoa)
3174 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zthl, &
3175 zdthladj, zta)
3176 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, po, pdoadj, &
3177 zoa)
3178 ELSE
3179 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
3180 zdhadj, zha)
3181 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
3182 pdoadj, zoa)
3183 END IF
3184
3185 IF (1==0) THEN
3186 CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
3187 zu, zv, pduadj, pdvadj, zua, zva)
3188 ELSE
3189 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
3190 zua)
3191 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
3192 zva)
3193 END IF
3194
3195 DO l = 1, nlay
3196 DO ig = 1, ngrid
3197 zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
3198 zf2 = zf/(1.-zf)
3199 thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l))**2
3200 wth2(ig, l) = zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
3201 END DO
3202 END DO
3203
3204
3205
3206 ! print*,'13 OK convect8'
3207 ! print*,'WA5 ',wa_moy
3208 DO l = 1, nlay
3209 DO ig = 1, ngrid
3210 ! pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
3211 pdtadj(ig, l) = zdthladj(ig, l)*zpspsk(ig, l)
3212 END DO
3213 END DO
3214
3215
3216 ! do l=1,nlay
3217 ! do ig=1,ngrid
3218 ! if(abs(pdtadj(ig,l))*86400..gt.500.) then
3219 ! print*,'WARN!!! ig=',ig,' l=',l
3220 ! s ,' pdtadj=',pdtadj(ig,l)
3221 ! endif
3222 ! if(abs(pdoadj(ig,l))*86400..gt.1.) then
3223 ! print*,'WARN!!! ig=',ig,' l=',l
3224 ! s ,' pdoadj=',pdoadj(ig,l)
3225 ! endif
3226 ! enddo
3227 ! enddo
3228
3229 ! print*,'14 OK convect8'
3230 ! ------------------------------------------------------------------
3231 ! Calculs pour les sorties
3232 ! ------------------------------------------------------------------
3233
3234 RETURN
3235 END SUBROUTINE thermcell_eau
3236
3237 SUBROUTINE thermcell(ngrid, nlay, ptimestep, pplay, pplev, pphi, pu, pv, pt, &
3238 po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0 & ! s
3239 ! ,pu_therm,pv_therm
3240 , r_aspect, l_mix, w2di, tho)
3241
3242 USE dimphy
3243 IMPLICIT NONE
3244
3245 ! =======================================================================
3246
3247 ! Calcul du transport verticale dans la couche limite en presence
3248 ! de "thermiques" explicitement representes
3249
3250 ! R��criture � partir d'un listing papier � Habas, le 14/02/00
3251
3252 ! le thermique est suppos� homog�ne et dissip� par m�lange avec
3253 ! son environnement. la longueur l_mix contr�le l'efficacit� du
3254 ! m�lange
3255
3256 ! Le calcul du transport des diff�rentes esp�ces se fait en prenant
3257 ! en compte:
3258 ! 1. un flux de masse montant
3259 ! 2. un flux de masse descendant
3260 ! 3. un entrainement
3261 ! 4. un detrainement
3262
3263 ! =======================================================================
3264
3265 ! -----------------------------------------------------------------------
3266 ! declarations:
3267 ! -------------
3268
3269 include "YOMCST.h"
3270
3271 ! arguments:
3272 ! ----------
3273
3274 INTEGER ngrid, nlay, w2di
3275 REAL tho
3276 REAL ptimestep, l_mix, r_aspect
3277 REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
3278 REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
3279 REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
3280 REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
3281 REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
3282 REAL pphi(ngrid, nlay)
3283
3284 INTEGER idetr
3285 SAVE idetr
3286 DATA idetr/3/
3287 !$OMP THREADPRIVATE(idetr)
3288
3289 ! local:
3290 ! ------
3291
3292 INTEGER ig, k, l, lmaxa(klon), lmix(klon)
3293 REAL zsortie1d(klon)
3294 ! CR: on remplace lmax(klon,klev+1)
3295 INTEGER lmax(klon), lmin(klon), lentr(klon)
3296 REAL linter(klon)
3297 REAL zmix(klon), fracazmix(klon)
3298 ! RC
3299 REAL zmax(klon), zw, zz, zw2(klon, klev+1), ztva(klon, klev), zzz
3300
3301 REAL zlev(klon, klev+1), zlay(klon, klev)
3302 REAL zh(klon, klev), zdhadj(klon, klev)
3303 REAL ztv(klon, klev)
3304 REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
3305 REAL wh(klon, klev+1)
3306 REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
3307 REAL zla(klon, klev+1)
3308 REAL zwa(klon, klev+1)
3309 REAL zld(klon, klev+1)
3310 REAL zwd(klon, klev+1)
3311 REAL zsortie(klon, klev)
3312 REAL zva(klon, klev)
3313 REAL zua(klon, klev)
3314 REAL zoa(klon, klev)
3315
3316 REAL zha(klon, klev)
3317 REAL wa_moy(klon, klev+1)
3318 REAL fraca(klon, klev+1)
3319 REAL fracc(klon, klev+1)
3320 REAL zf, zf2
3321 REAL thetath2(klon, klev), wth2(klon, klev)
3322 ! common/comtherm/thetath2,wth2
3323
3324 REAL count_time
3325 INTEGER ialt
3326
3327 LOGICAL sorties
3328 REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
3329 REAL zpspsk(klon, klev)
3330
3331 ! real wmax(klon,klev),wmaxa(klon)
3332 REAL wmax(klon), wmaxa(klon)
3333 REAL wa(klon, klev, klev+1)
3334 REAL wd(klon, klev+1)
3335 REAL larg_part(klon, klev, klev+1)
3336 REAL fracd(klon, klev+1)
3337 REAL xxx(klon, klev+1)
3338 REAL larg_cons(klon, klev+1)
3339 REAL larg_detr(klon, klev+1)
3340 REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
3341 REAL pu_therm(klon, klev), pv_therm(klon, klev)
3342 REAL fm(klon, klev+1), entr(klon, klev)
3343 REAL fmc(klon, klev+1)
3344
3345 ! CR:nouvelles variables
3346 REAL f_star(klon, klev+1), entr_star(klon, klev)
3347 REAL entr_star_tot(klon), entr_star2(klon)
3348 REAL f(klon), f0(klon)
3349 REAL zlevinter(klon)
3350 LOGICAL first
3351 DATA first/.FALSE./
3352 SAVE first
3353 !$OMP THREADPRIVATE(first)
3354 ! RC
3355
3356 CHARACTER *2 str2
3357 CHARACTER *10 str10
3358
3359 CHARACTER (LEN=20) :: modname = 'thermcell'
3360 CHARACTER (LEN=80) :: abort_message
3361
3362 LOGICAL vtest(klon), down
3363
3364 EXTERNAL scopy
3365
3366 INTEGER ncorrec, ll
3367 SAVE ncorrec
3368 DATA ncorrec/0/
3369 !$OMP THREADPRIVATE(ncorrec)
3370
3371
3372 ! -----------------------------------------------------------------------
3373 ! initialisation:
3374 ! ---------------
3375
3376 sorties = .TRUE.
3377 IF (ngrid/=klon) THEN
3378 PRINT *
3379 PRINT *, 'STOP dans convadj'
3380 PRINT *, 'ngrid =', ngrid
3381 PRINT *, 'klon =', klon
3382 END IF
3383
3384 ! -----------------------------------------------------------------------
3385 ! incrementation eventuelle de tendances precedentes:
3386 ! ---------------------------------------------------
3387
3388 ! print*,'0 OK convect8'
3389
3390 DO l = 1, nlay
3391 DO ig = 1, ngrid
3392 zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa
3393 zh(ig, l) = pt(ig, l)/zpspsk(ig, l)
3394 zu(ig, l) = pu(ig, l)
3395 zv(ig, l) = pv(ig, l)
3396 zo(ig, l) = po(ig, l)
3397 ztv(ig, l) = zh(ig, l)*(1.+0.61*zo(ig,l))
3398 END DO
3399 END DO
3400
3401 ! print*,'1 OK convect8'
3402 ! --------------------
3403
3404
3405 ! + + + + + + + + + + +
3406
3407
3408 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
3409 ! wh,wt,wo ...
3410
3411 ! + + + + + + + + + + + zh,zu,zv,zo,rho
3412
3413
3414 ! -------------------- zlev(1)
3415 ! \\\\\\\\\\\\\\\\\\\
3416
3417
3418 ! -----------------------------------------------------------------------
3419 ! Calcul des altitudes des couches
3420 ! -----------------------------------------------------------------------
3421
3422 DO l = 2, nlay
3423 DO ig = 1, ngrid
3424 zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
3425 END DO
3426 END DO
3427 DO ig = 1, ngrid
3428 zlev(ig, 1) = 0.
3429 zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
3430 END DO
3431 DO l = 1, nlay
3432 DO ig = 1, ngrid
3433 zlay(ig, l) = pphi(ig, l)/rg
3434 END DO
3435 END DO
3436
3437 ! print*,'2 OK convect8'
3438 ! -----------------------------------------------------------------------
3439 ! Calcul des densites
3440 ! -----------------------------------------------------------------------
3441
3442 DO l = 1, nlay
3443 DO ig = 1, ngrid
3444 rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*zh(ig,l))
3445 END DO
3446 END DO
3447
3448 DO l = 2, nlay
3449 DO ig = 1, ngrid
3450 rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
3451 END DO
3452 END DO
3453
3454 DO k = 1, nlay
3455 DO l = 1, nlay + 1
3456 DO ig = 1, ngrid
3457 wa(ig, k, l) = 0.
3458 END DO
3459 END DO
3460 END DO
3461
3462 ! print*,'3 OK convect8'
3463 ! ------------------------------------------------------------------
3464 ! Calcul de w2, quarre de w a partir de la cape
3465 ! a partir de w2, on calcule wa, vitesse de l'ascendance
3466
3467 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
3468 ! w2 est stoke dans wa
3469
3470 ! ATTENTION: dans convect8, on n'utilise le calcule des wa
3471 ! independants par couches que pour calculer l'entrainement
3472 ! a la base et la hauteur max de l'ascendance.
3473
3474 ! Indicages:
3475 ! l'ascendance provenant du niveau k traverse l'interface l avec
3476 ! une vitesse wa(k,l).
3477
3478 ! --------------------
3479
3480 ! + + + + + + + + + +
3481
3482 ! wa(k,l) ---- -------------------- l
3483 ! / ! /||\ + + + + + + + + + +
3484 ! ||
3485 ! || --------------------
3486 ! ||
3487 ! || + + + + + + + + + +
3488 ! ||
3489 ! || --------------------
3490 ! ||__
3491 ! |___ + + + + + + + + + + k
3492
3493 ! --------------------
3494
3495
3496
3497 ! ------------------------------------------------------------------
3498
3499 ! CR: ponderation entrainement des couches instables
3500 ! def des entr_star tels que entr=f*entr_star
3501 DO l = 1, klev
3502 DO ig = 1, ngrid
3503 entr_star(ig, l) = 0.
3504 END DO
3505 END DO
3506 ! determination de la longueur de la couche d entrainement
3507 DO ig = 1, ngrid
3508 lentr(ig) = 1
3509 END DO
3510
3511 ! on ne considere que les premieres couches instables
3512 DO k = nlay - 2, 1, -1
3513 DO ig = 1, ngrid
3514 IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN
3515 lentr(ig) = k
3516 END IF
3517 END DO
3518 END DO
3519
3520 ! determination du lmin: couche d ou provient le thermique
3521 DO ig = 1, ngrid
3522 lmin(ig) = 1
3523 END DO
3524 DO ig = 1, ngrid
3525 DO l = nlay, 2, -1
3526 IF (ztv(ig,l-1)>ztv(ig,l)) THEN
3527 lmin(ig) = l - 1
3528 END IF
3529 END DO
3530 END DO
3531
3532 ! definition de l'entrainement des couches
3533 DO l = 1, klev - 1
3534 DO ig = 1, ngrid
3535 IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<=lentr(ig)) THEN
3536 entr_star(ig, l) = (ztv(ig,l)-ztv(ig,l+1))*(zlev(ig,l+1)-zlev(ig,l))
3537 END IF
3538 END DO
3539 END DO
3540 ! pas de thermique si couches 1->5 stables
3541 DO ig = 1, ngrid
3542 IF (lmin(ig)>5) THEN
3543 DO l = 1, klev
3544 entr_star(ig, l) = 0.
3545 END DO
3546 END IF
3547 END DO
3548 ! calcul de l entrainement total
3549 DO ig = 1, ngrid
3550 entr_star_tot(ig) = 0.
3551 END DO
3552 DO ig = 1, ngrid
3553 DO k = 1, klev
3554 entr_star_tot(ig) = entr_star_tot(ig) + entr_star(ig, k)
3555 END DO
3556 END DO
3557
3558 PRINT *, 'fin calcul entr_star'
3559 DO k = 1, klev
3560 DO ig = 1, ngrid
3561 ztva(ig, k) = ztv(ig, k)
3562 END DO
3563 END DO
3564 ! RC
3565 ! print*,'7 OK convect8'
3566 DO k = 1, klev + 1
3567 DO ig = 1, ngrid
3568 zw2(ig, k) = 0.
3569 fmc(ig, k) = 0.
3570 ! CR
3571 f_star(ig, k) = 0.
3572 ! RC
3573 larg_cons(ig, k) = 0.
3574 larg_detr(ig, k) = 0.
3575 wa_moy(ig, k) = 0.
3576 END DO
3577 END DO
3578
3579 ! print*,'8 OK convect8'
3580 DO ig = 1, ngrid
3581 linter(ig) = 1.
3582 lmaxa(ig) = 1
3583 lmix(ig) = 1
3584 wmaxa(ig) = 0.
3585 END DO
3586
3587 ! CR:
3588 DO l = 1, nlay - 2
3589 DO ig = 1, ngrid
3590 IF (ztv(ig,l)>ztv(ig,l+1) .AND. entr_star(ig,l)>1.E-10 .AND. &
3591 zw2(ig,l)<1E-10) THEN
3592 f_star(ig, l+1) = entr_star(ig, l)
3593 ! test:calcul de dteta
3594 zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
3595 (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
3596 larg_detr(ig, l) = 0.
3597 ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+entr_star(ig, &
3598 l)>1.E-10)) THEN
3599 f_star(ig, l+1) = f_star(ig, l) + entr_star(ig, l)
3600 ztva(ig, l) = (f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)*ztv(ig,l))/ &
3601 f_star(ig, l+1)
3602 zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/f_star(ig,l+1))**2 + &
3603 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
3604 END IF
3605 ! determination de zmax continu par interpolation lineaire
3606 IF (zw2(ig,l+1)<0.) THEN
3607 ! test
3608 IF (abs(zw2(ig,l+1)-zw2(ig,l))<1E-10) THEN
3609 PRINT *, 'pb linter'
3610 END IF
3611 linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
3612 ig,l))
3613 zw2(ig, l+1) = 0.
3614 lmaxa(ig) = l
3615 ELSE
3616 IF (zw2(ig,l+1)<0.) THEN
3617 PRINT *, 'pb1 zw2<0'
3618 END IF
3619 wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
3620 END IF
3621 IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
3622 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
3623 lmix(ig) = l + 1
3624 wmaxa(ig) = wa_moy(ig, l+1)
3625 END IF
3626 END DO
3627 END DO
3628 PRINT *, 'fin calcul zw2'
3629
3630 ! Calcul de la couche correspondant a la hauteur du thermique
3631 DO ig = 1, ngrid
3632 lmax(ig) = lentr(ig)
3633 END DO
3634 DO ig = 1, ngrid
3635 DO l = nlay, lentr(ig) + 1, -1
3636 IF (zw2(ig,l)<=1.E-10) THEN
3637 lmax(ig) = l - 1
3638 END IF
3639 END DO
3640 END DO
3641 ! pas de thermique si couches 1->5 stables
3642 DO ig = 1, ngrid
3643 IF (lmin(ig)>5) THEN
3644 lmax(ig) = 1
3645 lmin(ig) = 1
3646 END IF
3647 END DO
3648
3649 ! Determination de zw2 max
3650 DO ig = 1, ngrid
3651 wmax(ig) = 0.
3652 END DO
3653
3654 DO l = 1, nlay
3655 DO ig = 1, ngrid
3656 IF (l<=lmax(ig)) THEN
3657 IF (zw2(ig,l)<0.) THEN
3658 PRINT *, 'pb2 zw2<0'
3659 END IF
3660 zw2(ig, l) = sqrt(zw2(ig,l))
3661 wmax(ig) = max(wmax(ig), zw2(ig,l))
3662 ELSE
3663 zw2(ig, l) = 0.
3664 END IF
3665 END DO
3666 END DO
3667
3668 ! Longueur caracteristique correspondant a la hauteur des thermiques.
3669 DO ig = 1, ngrid
3670 zmax(ig) = 0.
3671 zlevinter(ig) = zlev(ig, 1)
3672 END DO
3673 DO ig = 1, ngrid
3674 ! calcul de zlevinter
3675 zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
3676 zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
3677 zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,lmin(ig)))
3678 END DO
3679
3680 PRINT *, 'avant fermeture'
3681 ! Fermeture,determination de f
3682 DO ig = 1, ngrid
3683 entr_star2(ig) = 0.
3684 END DO
3685 DO ig = 1, ngrid
3686 IF (entr_star_tot(ig)<1.E-10) THEN
3687 f(ig) = 0.
3688 ELSE
3689 DO k = lmin(ig), lentr(ig)
3690 entr_star2(ig) = entr_star2(ig) + entr_star(ig, k)**2/(rho(ig,k)*( &
3691 zlev(ig,k+1)-zlev(ig,k)))
3692 END DO
3693 ! Nouvelle fermeture
3694 f(ig) = wmax(ig)/(max(500.,zmax(ig))*r_aspect*entr_star2(ig))* &
3695 entr_star_tot(ig)
3696 ! test
3697 ! if (first) then
3698 ! f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
3699 ! s *wmax(ig))
3700 ! endif
3701 END IF
3702 ! f0(ig)=f(ig)
3703 ! first=.true.
3704 END DO
3705 PRINT *, 'apres fermeture'
3706
3707 ! Calcul de l'entrainement
3708 DO k = 1, klev
3709 DO ig = 1, ngrid
3710 entr(ig, k) = f(ig)*entr_star(ig, k)
3711 END DO
3712 END DO
3713 ! Calcul des flux
3714 DO ig = 1, ngrid
3715 DO l = 1, lmax(ig) - 1
3716 fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
3717 END DO
3718 END DO
3719
3720 ! RC
3721
3722
3723 ! print*,'9 OK convect8'
3724 ! print*,'WA1 ',wa_moy
3725
3726 ! determination de l'indice du debut de la mixed layer ou w decroit
3727
3728 ! calcul de la largeur de chaque ascendance dans le cas conservatif.
3729 ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
3730 ! d'une couche est �gale � la hauteur de la couche alimentante.
3731 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
3732 ! de la vitesse d'entrainement horizontal dans la couche alimentante.
3733
3734 DO l = 2, nlay
3735 DO ig = 1, ngrid
3736 IF (l<=lmaxa(ig)) THEN
3737 zw = max(wa_moy(ig,l), 1.E-10)
3738 larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw)
3739 END IF
3740 END DO
3741 END DO
3742
3743 DO l = 2, nlay
3744 DO ig = 1, ngrid
3745 IF (l<=lmaxa(ig)) THEN
3746 ! if (idetr.eq.0) then
3747 ! cette option est finalement en dur.
3748 IF ((l_mix*zlev(ig,l))<0.) THEN
3749 PRINT *, 'pb l_mix*zlev<0'
3750 END IF
3751 larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
3752 ! else if (idetr.eq.1) then
3753 ! larg_detr(ig,l)=larg_cons(ig,l)
3754 ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
3755 ! else if (idetr.eq.2) then
3756 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3757 ! s *sqrt(wa_moy(ig,l))
3758 ! else if (idetr.eq.4) then
3759 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
3760 ! s *wa_moy(ig,l)
3761 ! endif
3762 END IF
3763 END DO
3764 END DO
3765
3766 ! print*,'10 OK convect8'
3767 ! print*,'WA2 ',wa_moy
3768 ! calcul de la fraction de la maille concern�e par l'ascendance en tenant
3769 ! compte de l'epluchage du thermique.
3770
3771 ! CR def de zmix continu (profil parabolique des vitesses)
3772 DO ig = 1, ngrid
3773 IF (lmix(ig)>1.) THEN
3774 ! test
3775 IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
3776 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
3777 zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- &
3778 (zlev(ig,lmix(ig)))))>1E-10) THEN
3779
3780 zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) &
3781 )**2-(zlev(ig,lmix(ig)+1))**2)-(zw2(ig,lmix(ig))-zw2(ig, &
3782 lmix(ig)+1))*((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))/ &
3783 (2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
3784 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
3785 zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
3786 ELSE
3787 zmix(ig) = zlev(ig, lmix(ig))
3788 PRINT *, 'pb zmix'
3789 END IF
3790 ELSE
3791 zmix(ig) = 0.
3792 END IF
3793 ! test
3794 IF ((zmax(ig)-zmix(ig))<0.) THEN
3795 zmix(ig) = 0.99*zmax(ig)
3796 ! print*,'pb zmix>zmax'
3797 END IF
3798 END DO
3799
3800 ! calcul du nouveau lmix correspondant
3801 DO ig = 1, ngrid
3802 DO l = 1, klev
3803 IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
3804 lmix(ig) = l
3805 END IF
3806 END DO
3807 END DO
3808
3809 DO l = 2, nlay
3810 DO ig = 1, ngrid
3811 IF (larg_cons(ig,l)>1.) THEN
3812 ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
3813 fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
3814 ! test
3815 fraca(ig, l) = max(fraca(ig,l), 0.)
3816 fraca(ig, l) = min(fraca(ig,l), 0.5)
3817 fracd(ig, l) = 1. - fraca(ig, l)
3818 fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
3819 ELSE
3820 ! wa_moy(ig,l)=0.
3821 fraca(ig, l) = 0.
3822 fracc(ig, l) = 0.
3823 fracd(ig, l) = 1.
3824 END IF
3825 END DO
3826 END DO
3827 ! CR: calcul de fracazmix
3828 DO ig = 1, ngrid
3829 fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
3830 (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
3831 fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca(ig &
3832 ,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
3833 END DO
3834
3835 DO l = 2, nlay
3836 DO ig = 1, ngrid
3837 IF (larg_cons(ig,l)>1.) THEN
3838 IF (l>lmix(ig)) THEN
3839 ! test
3840 IF (zmax(ig)-zmix(ig)<1.E-10) THEN
3841 ! print*,'pb xxx'
3842 xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
3843 ELSE
3844 xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
3845 END IF
3846 IF (idetr==0) THEN
3847 fraca(ig, l) = fracazmix(ig)
3848 ELSE IF (idetr==1) THEN
3849 fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
3850 ELSE IF (idetr==2) THEN
3851 fraca(ig, l) = fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
3852 ELSE
3853 fraca(ig, l) = fracazmix(ig)*xxx(ig, l)**2
3854 END IF
3855 ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
3856 fraca(ig, l) = max(fraca(ig,l), 0.)
3857 fraca(ig, l) = min(fraca(ig,l), 0.5)
3858 fracd(ig, l) = 1. - fraca(ig, l)
3859 fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
3860 END IF
3861 END IF
3862 END DO
3863 END DO
3864
3865 PRINT *, 'fin calcul fraca'
3866 ! print*,'11 OK convect8'
3867 ! print*,'Ea3 ',wa_moy
3868 ! ------------------------------------------------------------------
3869 ! Calcul de fracd, wd
3870 ! somme wa - wd = 0
3871 ! ------------------------------------------------------------------
3872
3873
3874 DO ig = 1, ngrid
3875 fm(ig, 1) = 0.
3876 fm(ig, nlay+1) = 0.
3877 END DO
3878
3879 DO l = 2, nlay
3880 DO ig = 1, ngrid
3881 fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
3882 ! CR:test
3883 IF (entr(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) THEN
3884 fm(ig, l) = fm(ig, l-1)
3885 ! write(1,*)'ajustement fm, l',l
3886 END IF
3887 ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
3888 ! RC
3889 END DO
3890 DO ig = 1, ngrid
3891 IF (fracd(ig,l)<0.1) THEN
3892 abort_message = 'fracd trop petit'
3893 CALL abort_physic(modname, abort_message, 1)
3894 ELSE
3895 ! vitesse descendante "diagnostique"
3896 wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
3897 END IF
3898 END DO
3899 END DO
3900
3901 DO l = 1, nlay
3902 DO ig = 1, ngrid
3903 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
3904 masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
3905 END DO
3906 END DO
3907
3908 ! print*,'12 OK convect8'
3909 ! print*,'WA4 ',wa_moy
3910 ! c------------------------------------------------------------------
3911 ! calcul du transport vertical
3912 ! ------------------------------------------------------------------
3913
3914 GO TO 4444
3915 ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
3916 DO l = 2, nlay - 1
3917 DO ig = 1, ngrid
3918 IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( &
3919 ig,l+1)) THEN
3920 ! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
3921 ! s ,fm(ig,l+1)*ptimestep
3922 ! s ,' M=',masse(ig,l),masse(ig,l+1)
3923 END IF
3924 END DO
3925 END DO
3926
3927 DO l = 1, nlay
3928 DO ig = 1, ngrid
3929 IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
3930 ! print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
3931 ! s ,entr(ig,l)*ptimestep
3932 ! s ,' M=',masse(ig,l)
3933 END IF
3934 END DO
3935 END DO
3936
3937 DO l = 1, nlay
3938 DO ig = 1, ngrid
3939 IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
3940 ! print*,'WARN!!! fm exagere ig=',ig,' l=',l
3941 ! s ,' FM=',fm(ig,l)
3942 END IF
3943 IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
3944 ! print*,'WARN!!! masse exagere ig=',ig,' l=',l
3945 ! s ,' M=',masse(ig,l)
3946 ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
3947 ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
3948 ! print*,'zlev(ig,l+1),zlev(ig,l)'
3949 ! s ,zlev(ig,l+1),zlev(ig,l)
3950 ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
3951 ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
3952 END IF
3953 IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
3954 ! print*,'WARN!!! entr exagere ig=',ig,' l=',l
3955 ! s ,' E=',entr(ig,l)
3956 END IF
3957 END DO
3958 END DO
3959
3960 4444 CONTINUE
3961
3962 ! CR:redefinition du entr
3963 DO l = 1, nlay
3964 DO ig = 1, ngrid
3965 detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1)
3966 IF (detr(ig,l)<0.) THEN
3967 entr(ig, l) = entr(ig, l) - detr(ig, l)
3968 detr(ig, l) = 0.
3969 ! print*,'WARNING !!! detrainement negatif ',ig,l
3970 END IF
3971 END DO
3972 END DO
3973 ! RC
3974 IF (w2di==1) THEN
3975 fm0 = fm0 + ptimestep*(fm-fm0)/tho
3976 entr0 = entr0 + ptimestep*(entr-entr0)/tho
3977 ELSE
3978 fm0 = fm
3979 entr0 = entr
3980 END IF
3981
3982 IF (1==1) THEN
3983 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, &
3984 zha)
3985 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, &
3986 zoa)
3987 ELSE
3988 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
3989 zdhadj, zha)
3990 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
3991 pdoadj, zoa)
3992 END IF
3993
3994 IF (1==0) THEN
3995 CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
3996 zu, zv, pduadj, pdvadj, zua, zva)
3997 ELSE
3998 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
3999 zua)
4000 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
4001 zva)
4002 END IF
4003
4004 DO l = 1, nlay
4005 DO ig = 1, ngrid
4006 zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
4007 zf2 = zf/(1.-zf)
4008 thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l))**2
4009 wth2(ig, l) = zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
4010 END DO
4011 END DO
4012
4013
4014
4015 ! print*,'13 OK convect8'
4016 ! print*,'WA5 ',wa_moy
4017 DO l = 1, nlay
4018 DO ig = 1, ngrid
4019 pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l)
4020 END DO
4021 END DO
4022
4023
4024 ! do l=1,nlay
4025 ! do ig=1,ngrid
4026 ! if(abs(pdtadj(ig,l))*86400..gt.500.) then
4027 ! print*,'WARN!!! ig=',ig,' l=',l
4028 ! s ,' pdtadj=',pdtadj(ig,l)
4029 ! endif
4030 ! if(abs(pdoadj(ig,l))*86400..gt.1.) then
4031 ! print*,'WARN!!! ig=',ig,' l=',l
4032 ! s ,' pdoadj=',pdoadj(ig,l)
4033 ! endif
4034 ! enddo
4035 ! enddo
4036
4037 ! print*,'14 OK convect8'
4038 ! ------------------------------------------------------------------
4039 ! Calculs pour les sorties
4040 ! ------------------------------------------------------------------
4041
4042 IF (sorties) THEN
4043 DO l = 1, nlay
4044 DO ig = 1, ngrid
4045 zla(ig, l) = (1.-fracd(ig,l))*zmax(ig)
4046 zld(ig, l) = fracd(ig, l)*zmax(ig)
4047 IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ &
4048 (1.-fracd(ig,l))
4049 END DO
4050 END DO
4051
4052 ! deja fait
4053 ! do l=1,nlay
4054 ! do ig=1,ngrid
4055 ! detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
4056 ! if (detr(ig,l).lt.0.) then
4057 ! entr(ig,l)=entr(ig,l)-detr(ig,l)
4058 ! detr(ig,l)=0.
4059 ! print*,'WARNING !!! detrainement negatif ',ig,l
4060 ! endif
4061 ! enddo
4062 ! enddo
4063
4064 ! print*,'15 OK convect8'
4065
4066
4067 ! #define und
4068 GO TO 123
4069 123 CONTINUE
4070
4071 END IF
4072
4073 ! if(wa_moy(1,4).gt.1.e-10) stop
4074
4075 ! print*,'19 OK convect8'
4076 RETURN
4077 END SUBROUTINE thermcell
4078
4079 SUBROUTINE dqthermcell(ngrid, nlay, ptimestep, fm, entr, masse, q, dq, qa)
4080 USE dimphy
4081 IMPLICIT NONE
4082
4083 ! =======================================================================
4084
4085 ! Calcul du transport verticale dans la couche limite en presence
4086 ! de "thermiques" explicitement representes
4087 ! calcul du dq/dt une fois qu'on connait les ascendances
4088
4089 ! =======================================================================
4090
4091 INTEGER ngrid, nlay
4092
4093 REAL ptimestep
4094 REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4095 REAL entr(ngrid, nlay)
4096 REAL q(ngrid, nlay)
4097 REAL dq(ngrid, nlay)
4098
4099 REAL qa(klon, klev), detr(klon, klev), wqd(klon, klev+1)
4100
4101 INTEGER ig, k
4102
4103 ! calcul du detrainement
4104
4105 DO k = 1, nlay
4106 DO ig = 1, ngrid
4107 detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4108 ! test
4109 IF (detr(ig,k)<0.) THEN
4110 entr(ig, k) = entr(ig, k) - detr(ig, k)
4111 detr(ig, k) = 0.
4112 ! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
4113 ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
4114 END IF
4115 IF (fm(ig,k+1)<0.) THEN
4116 ! print*,'fm2<0!!!'
4117 END IF
4118 IF (entr(ig,k)<0.) THEN
4119 ! print*,'entr2<0!!!'
4120 END IF
4121 END DO
4122 END DO
4123
4124 ! calcul de la valeur dans les ascendances
4125 DO ig = 1, ngrid
4126 qa(ig, 1) = q(ig, 1)
4127 END DO
4128
4129 DO k = 2, nlay
4130 DO ig = 1, ngrid
4131 IF ((fm(ig,k+1)+detr(ig,k))*ptimestep>1.E-5*masse(ig,k)) THEN
4132 qa(ig, k) = (fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))/ &
4133 (fm(ig,k+1)+detr(ig,k))
4134 ELSE
4135 qa(ig, k) = q(ig, k)
4136 END IF
4137 IF (qa(ig,k)<0.) THEN
4138 ! print*,'qa<0!!!'
4139 END IF
4140 IF (q(ig,k)<0.) THEN
4141 ! print*,'q<0!!!'
4142 END IF
4143 END DO
4144 END DO
4145
4146 DO k = 2, nlay
4147 DO ig = 1, ngrid
4148 ! wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
4149 wqd(ig, k) = fm(ig, k)*q(ig, k)
4150 IF (wqd(ig,k)<0.) THEN
4151 ! print*,'wqd<0!!!'
4152 END IF
4153 END DO
4154 END DO
4155 DO ig = 1, ngrid
4156 wqd(ig, 1) = 0.
4157 wqd(ig, nlay+1) = 0.
4158 END DO
4159
4160 DO k = 1, nlay
4161 DO ig = 1, ngrid
4162 dq(ig, k) = (detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)-wqd(ig,k)+wqd(ig,k+ &
4163 1))/masse(ig, k)
4164 ! if (dq(ig,k).lt.0.) then
4165 ! print*,'dq<0!!!'
4166 ! endif
4167 END DO
4168 END DO
4169
4170 RETURN
4171 END SUBROUTINE dqthermcell
4172 SUBROUTINE dvthermcell(ngrid, nlay, ptimestep, fm, entr, masse, fraca, larga, &
4173 u, v, du, dv, ua, va)
4174 USE dimphy
4175 IMPLICIT NONE
4176
4177 ! =======================================================================
4178
4179 ! Calcul du transport verticale dans la couche limite en presence
4180 ! de "thermiques" explicitement representes
4181 ! calcul du dq/dt une fois qu'on connait les ascendances
4182
4183 ! =======================================================================
4184
4185 INTEGER ngrid, nlay
4186
4187 REAL ptimestep
4188 REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4189 REAL fraca(ngrid, nlay+1)
4190 REAL larga(ngrid)
4191 REAL entr(ngrid, nlay)
4192 REAL u(ngrid, nlay)
4193 REAL ua(ngrid, nlay)
4194 REAL du(ngrid, nlay)
4195 REAL v(ngrid, nlay)
4196 REAL va(ngrid, nlay)
4197 REAL dv(ngrid, nlay)
4198
4199 REAL qa(klon, klev), detr(klon, klev)
4200 REAL wvd(klon, klev+1), wud(klon, klev+1)
4201 REAL gamma0, gamma(klon, klev+1)
4202 REAL dua, dva
4203 INTEGER iter
4204
4205 INTEGER ig, k
4206
4207 ! calcul du detrainement
4208
4209 DO k = 1, nlay
4210 DO ig = 1, ngrid
4211 detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4212 END DO
4213 END DO
4214
4215 ! calcul de la valeur dans les ascendances
4216 DO ig = 1, ngrid
4217 ua(ig, 1) = u(ig, 1)
4218 va(ig, 1) = v(ig, 1)
4219 END DO
4220
4221 DO k = 2, nlay
4222 DO ig = 1, ngrid
4223 IF ((fm(ig,k+1)+detr(ig,k))*ptimestep>1.E-5*masse(ig,k)) THEN
4224 ! On it�re sur la valeur du coeff de freinage.
4225 ! gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
4226 gamma0 = masse(ig, k)*sqrt(0.5*(fraca(ig,k+1)+fraca(ig, &
4227 k)))*0.5/larga(ig)
4228 ! gamma0=0.
4229 ! la premi�re fois on multiplie le coefficient de freinage
4230 ! par le module du vent dans la couche en dessous.
4231 dua = ua(ig, k-1) - u(ig, k-1)
4232 dva = va(ig, k-1) - v(ig, k-1)
4233 DO iter = 1, 5
4234 gamma(ig, k) = gamma0*sqrt(dua**2+dva**2)
4235 ua(ig, k) = (fm(ig,k)*ua(ig,k-1)+(entr(ig,k)+gamma(ig, &
4236 k))*u(ig,k))/(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
4237 va(ig, k) = (fm(ig,k)*va(ig,k-1)+(entr(ig,k)+gamma(ig, &
4238 k))*v(ig,k))/(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
4239 ! print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
4240 dua = ua(ig, k) - u(ig, k)
4241 dva = va(ig, k) - v(ig, k)
4242 END DO
4243 ELSE
4244 ua(ig, k) = u(ig, k)
4245 va(ig, k) = v(ig, k)
4246 gamma(ig, k) = 0.
4247 END IF
4248 END DO
4249 END DO
4250
4251 DO k = 2, nlay
4252 DO ig = 1, ngrid
4253 wud(ig, k) = fm(ig, k)*u(ig, k)
4254 wvd(ig, k) = fm(ig, k)*v(ig, k)
4255 END DO
4256 END DO
4257 DO ig = 1, ngrid
4258 wud(ig, 1) = 0.
4259 wud(ig, nlay+1) = 0.
4260 wvd(ig, 1) = 0.
4261 wvd(ig, nlay+1) = 0.
4262 END DO
4263
4264 DO k = 1, nlay
4265 DO ig = 1, ngrid
4266 du(ig, k) = ((detr(ig,k)+gamma(ig,k))*ua(ig,k)-(entr(ig,k)+gamma(ig, &
4267 k))*u(ig,k)-wud(ig,k)+wud(ig,k+1))/masse(ig, k)
4268 dv(ig, k) = ((detr(ig,k)+gamma(ig,k))*va(ig,k)-(entr(ig,k)+gamma(ig, &
4269 k))*v(ig,k)-wvd(ig,k)+wvd(ig,k+1))/masse(ig, k)
4270 END DO
4271 END DO
4272
4273 RETURN
4274 END SUBROUTINE dvthermcell
4275 SUBROUTINE dqthermcell2(ngrid, nlay, ptimestep, fm, entr, masse, frac, q, dq, &
4276 qa)
4277 USE dimphy
4278 IMPLICIT NONE
4279
4280 ! =======================================================================
4281
4282 ! Calcul du transport verticale dans la couche limite en presence
4283 ! de "thermiques" explicitement representes
4284 ! calcul du dq/dt une fois qu'on connait les ascendances
4285
4286 ! =======================================================================
4287
4288 INTEGER ngrid, nlay
4289
4290 REAL ptimestep
4291 REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4292 REAL entr(ngrid, nlay), frac(ngrid, nlay)
4293 REAL q(ngrid, nlay)
4294 REAL dq(ngrid, nlay)
4295
4296 REAL qa(klon, klev), detr(klon, klev), wqd(klon, klev+1)
4297 REAL qe(klon, klev), zf, zf2
4298
4299 INTEGER ig, k
4300
4301 ! calcul du detrainement
4302
4303 DO k = 1, nlay
4304 DO ig = 1, ngrid
4305 detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4306 END DO
4307 END DO
4308
4309 ! calcul de la valeur dans les ascendances
4310 DO ig = 1, ngrid
4311 qa(ig, 1) = q(ig, 1)
4312 qe(ig, 1) = q(ig, 1)
4313 END DO
4314
4315 DO k = 2, nlay
4316 DO ig = 1, ngrid
4317 IF ((fm(ig,k+1)+detr(ig,k))*ptimestep>1.E-5*masse(ig,k)) THEN
4318 zf = 0.5*(frac(ig,k)+frac(ig,k+1))
4319 zf2 = 1./(1.-zf)
4320 qa(ig, k) = (fm(ig,k)*qa(ig,k-1)+zf2*entr(ig,k)*q(ig,k))/ &
4321 (fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2)
4322 qe(ig, k) = (q(ig,k)-zf*qa(ig,k))*zf2
4323 ELSE
4324 qa(ig, k) = q(ig, k)
4325 qe(ig, k) = q(ig, k)
4326 END IF
4327 END DO
4328 END DO
4329
4330 DO k = 2, nlay
4331 DO ig = 1, ngrid
4332 ! wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
4333 wqd(ig, k) = fm(ig, k)*qe(ig, k)
4334 END DO
4335 END DO
4336 DO ig = 1, ngrid
4337 wqd(ig, 1) = 0.
4338 wqd(ig, nlay+1) = 0.
4339 END DO
4340
4341 DO k = 1, nlay
4342 DO ig = 1, ngrid
4343 dq(ig, k) = (detr(ig,k)*qa(ig,k)-entr(ig,k)*qe(ig,k)-wqd(ig,k)+wqd(ig,k &
4344 +1))/masse(ig, k)
4345 END DO
4346 END DO
4347
4348 RETURN
4349 END SUBROUTINE dqthermcell2
4350 SUBROUTINE dvthermcell2(ngrid, nlay, ptimestep, fm, entr, masse, fraca, &
4351 larga, u, v, du, dv, ua, va)
4352 USE dimphy
4353 IMPLICIT NONE
4354
4355 ! =======================================================================
4356
4357 ! Calcul du transport verticale dans la couche limite en presence
4358 ! de "thermiques" explicitement representes
4359 ! calcul du dq/dt une fois qu'on connait les ascendances
4360
4361 ! =======================================================================
4362
4363 INTEGER ngrid, nlay
4364
4365 REAL ptimestep
4366 REAL masse(ngrid, nlay), fm(ngrid, nlay+1)
4367 REAL fraca(ngrid, nlay+1)
4368 REAL larga(ngrid)
4369 REAL entr(ngrid, nlay)
4370 REAL u(ngrid, nlay)
4371 REAL ua(ngrid, nlay)
4372 REAL du(ngrid, nlay)
4373 REAL v(ngrid, nlay)
4374 REAL va(ngrid, nlay)
4375 REAL dv(ngrid, nlay)
4376
4377 REAL qa(klon, klev), detr(klon, klev), zf, zf2
4378 REAL wvd(klon, klev+1), wud(klon, klev+1)
4379 REAL gamma0, gamma(klon, klev+1)
4380 REAL ue(klon, klev), ve(klon, klev)
4381 REAL dua, dva
4382 INTEGER iter
4383
4384 INTEGER ig, k
4385
4386 ! calcul du detrainement
4387
4388 DO k = 1, nlay
4389 DO ig = 1, ngrid
4390 detr(ig, k) = fm(ig, k) - fm(ig, k+1) + entr(ig, k)
4391 END DO
4392 END DO
4393
4394 ! calcul de la valeur dans les ascendances
4395 DO ig = 1, ngrid
4396 ua(ig, 1) = u(ig, 1)
4397 va(ig, 1) = v(ig, 1)
4398 ue(ig, 1) = u(ig, 1)
4399 ve(ig, 1) = v(ig, 1)
4400 END DO
4401
4402 DO k = 2, nlay
4403 DO ig = 1, ngrid
4404 IF ((fm(ig,k+1)+detr(ig,k))*ptimestep>1.E-5*masse(ig,k)) THEN
4405 ! On it�re sur la valeur du coeff de freinage.
4406 ! gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
4407 gamma0 = masse(ig, k)*sqrt(0.5*(fraca(ig,k+1)+fraca(ig, &
4408 k)))*0.5/larga(ig)*1.
4409 ! s *0.5
4410 ! gamma0=0.
4411 zf = 0.5*(fraca(ig,k)+fraca(ig,k+1))
4412 zf = 0.
4413 zf2 = 1./(1.-zf)
4414 ! la premi�re fois on multiplie le coefficient de freinage
4415 ! par le module du vent dans la couche en dessous.
4416 dua = ua(ig, k-1) - u(ig, k-1)
4417 dva = va(ig, k-1) - v(ig, k-1)
4418 DO iter = 1, 5
4419 ! On choisit une relaxation lineaire.
4420 gamma(ig, k) = gamma0
4421 ! On choisit une relaxation quadratique.
4422 gamma(ig, k) = gamma0*sqrt(dua**2+dva**2)
4423 ua(ig, k) = (fm(ig,k)*ua(ig,k-1)+(zf2*entr(ig,k)+gamma(ig, &
4424 k))*u(ig,k))/(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2+gamma(ig,k) &
4425 )
4426 va(ig, k) = (fm(ig,k)*va(ig,k-1)+(zf2*entr(ig,k)+gamma(ig, &
4427 k))*v(ig,k))/(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2+gamma(ig,k) &
4428 )
4429 ! print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
4430 dua = ua(ig, k) - u(ig, k)
4431 dva = va(ig, k) - v(ig, k)
4432 ue(ig, k) = (u(ig,k)-zf*ua(ig,k))*zf2
4433 ve(ig, k) = (v(ig,k)-zf*va(ig,k))*zf2
4434 END DO
4435 ELSE
4436 ua(ig, k) = u(ig, k)
4437 va(ig, k) = v(ig, k)
4438 ue(ig, k) = u(ig, k)
4439 ve(ig, k) = v(ig, k)
4440 gamma(ig, k) = 0.
4441 END IF
4442 END DO
4443 END DO
4444
4445 DO k = 2, nlay
4446 DO ig = 1, ngrid
4447 wud(ig, k) = fm(ig, k)*ue(ig, k)
4448 wvd(ig, k) = fm(ig, k)*ve(ig, k)
4449 END DO
4450 END DO
4451 DO ig = 1, ngrid
4452 wud(ig, 1) = 0.
4453 wud(ig, nlay+1) = 0.
4454 wvd(ig, 1) = 0.
4455 wvd(ig, nlay+1) = 0.
4456 END DO
4457
4458 DO k = 1, nlay
4459 DO ig = 1, ngrid
4460 du(ig, k) = ((detr(ig,k)+gamma(ig,k))*ua(ig,k)-(entr(ig,k)+gamma(ig, &
4461 k))*ue(ig,k)-wud(ig,k)+wud(ig,k+1))/masse(ig, k)
4462 dv(ig, k) = ((detr(ig,k)+gamma(ig,k))*va(ig,k)-(entr(ig,k)+gamma(ig, &
4463 k))*ve(ig,k)-wvd(ig,k)+wvd(ig,k+1))/masse(ig, k)
4464 END DO
4465 END DO
4466
4467 RETURN
4468 END SUBROUTINE dvthermcell2
4469 SUBROUTINE thermcell_sec(ngrid, nlay, ptimestep, pplay, pplev, pphi, zlev, &
4470 pu, pv, pt, po, pduadj, pdvadj, pdtadj, pdoadj, fm0, entr0 & ! s
4471 ! ,pu_therm,pv_therm
4472 , r_aspect, l_mix, w2di, tho)
4473
4474 USE dimphy
4475 IMPLICIT NONE
4476
4477 ! =======================================================================
4478
4479 ! Calcul du transport verticale dans la couche limite en presence
4480 ! de "thermiques" explicitement representes
4481
4482 ! R��criture � partir d'un listing papier � Habas, le 14/02/00
4483
4484 ! le thermique est suppos� homog�ne et dissip� par m�lange avec
4485 ! son environnement. la longueur l_mix contr�le l'efficacit� du
4486 ! m�lange
4487
4488 ! Le calcul du transport des diff�rentes esp�ces se fait en prenant
4489 ! en compte:
4490 ! 1. un flux de masse montant
4491 ! 2. un flux de masse descendant
4492 ! 3. un entrainement
4493 ! 4. un detrainement
4494
4495 ! =======================================================================
4496
4497 ! -----------------------------------------------------------------------
4498 ! declarations:
4499 ! -------------
4500
4501 include "YOMCST.h"
4502
4503 ! arguments:
4504 ! ----------
4505
4506 INTEGER ngrid, nlay, w2di
4507 REAL tho
4508 REAL ptimestep, l_mix, r_aspect
4509 REAL pt(ngrid, nlay), pdtadj(ngrid, nlay)
4510 REAL pu(ngrid, nlay), pduadj(ngrid, nlay)
4511 REAL pv(ngrid, nlay), pdvadj(ngrid, nlay)
4512 REAL po(ngrid, nlay), pdoadj(ngrid, nlay)
4513 REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
4514 REAL pphi(ngrid, nlay)
4515
4516 INTEGER idetr
4517 SAVE idetr
4518 DATA idetr/3/
4519 !$OMP THREADPRIVATE(idetr)
4520
4521 ! local:
4522 ! ------
4523
4524 INTEGER ig, k, l, lmaxa(klon), lmix(klon)
4525 REAL zsortie1d(klon)
4526 ! CR: on remplace lmax(klon,klev+1)
4527 INTEGER lmax(klon), lmin(klon), lentr(klon)
4528 REAL linter(klon)
4529 REAL zmix(klon), fracazmix(klon)
4530 ! RC
4531 REAL zmax(klon), zw, zz, zw2(klon, klev+1), ztva(klon, klev), zzz
4532
4533 REAL zlev(klon, klev+1), zlay(klon, klev)
4534 REAL zh(klon, klev), zdhadj(klon, klev)
4535 REAL ztv(klon, klev)
4536 REAL zu(klon, klev), zv(klon, klev), zo(klon, klev)
4537 REAL wh(klon, klev+1)
4538 REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1)
4539 REAL zla(klon, klev+1)
4540 REAL zwa(klon, klev+1)
4541 REAL zld(klon, klev+1)
4542 REAL zwd(klon, klev+1)
4543 REAL zsortie(klon, klev)
4544 REAL zva(klon, klev)
4545 REAL zua(klon, klev)
4546 REAL zoa(klon, klev)
4547
4548 REAL zha(klon, klev)
4549 REAL wa_moy(klon, klev+1)
4550 REAL fraca(klon, klev+1)
4551 REAL fracc(klon, klev+1)
4552 REAL zf, zf2
4553 REAL thetath2(klon, klev), wth2(klon, klev)
4554 ! common/comtherm/thetath2,wth2
4555
4556 REAL count_time
4557 INTEGER ialt
4558
4559 LOGICAL sorties
4560 REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
4561 REAL zpspsk(klon, klev)
4562
4563 ! real wmax(klon,klev),wmaxa(klon)
4564 REAL wmax(klon), wmaxa(klon)
4565 REAL wa(klon, klev, klev+1)
4566 REAL wd(klon, klev+1)
4567 REAL larg_part(klon, klev, klev+1)
4568 REAL fracd(klon, klev+1)
4569 REAL xxx(klon, klev+1)
4570 REAL larg_cons(klon, klev+1)
4571 REAL larg_detr(klon, klev+1)
4572 REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
4573 REAL pu_therm(klon, klev), pv_therm(klon, klev)
4574 REAL fm(klon, klev+1), entr(klon, klev)
4575 REAL fmc(klon, klev+1)
4576
4577 ! CR:nouvelles variables
4578 REAL f_star(klon, klev+1), entr_star(klon, klev)
4579 REAL entr_star_tot(klon), entr_star2(klon)
4580 REAL f(klon), f0(klon)
4581 REAL zlevinter(klon)
4582 LOGICAL first
4583 DATA first/.FALSE./
4584 SAVE first
4585 !$OMP THREADPRIVATE(first)
4586 ! RC
4587
4588 CHARACTER *2 str2
4589 CHARACTER *10 str10
4590
4591 CHARACTER (LEN=20) :: modname = 'thermcell_sec'
4592 CHARACTER (LEN=80) :: abort_message
4593
4594 LOGICAL vtest(klon), down
4595
4596 EXTERNAL scopy
4597
4598 INTEGER ncorrec, ll
4599 SAVE ncorrec
4600 DATA ncorrec/0/
4601 !$OMP THREADPRIVATE(ncorrec)
4602
4603
4604 ! -----------------------------------------------------------------------
4605 ! initialisation:
4606 ! ---------------
4607
4608 sorties = .TRUE.
4609 IF (ngrid/=klon) THEN
4610 PRINT *
4611 PRINT *, 'STOP dans convadj'
4612 PRINT *, 'ngrid =', ngrid
4613 PRINT *, 'klon =', klon
4614 END IF
4615
4616 ! -----------------------------------------------------------------------
4617 ! incrementation eventuelle de tendances precedentes:
4618 ! ---------------------------------------------------
4619
4620 ! print*,'0 OK convect8'
4621
4622 DO l = 1, nlay
4623 DO ig = 1, ngrid
4624 zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa
4625 zh(ig, l) = pt(ig, l)/zpspsk(ig, l)
4626 zu(ig, l) = pu(ig, l)
4627 zv(ig, l) = pv(ig, l)
4628 zo(ig, l) = po(ig, l)
4629 ztv(ig, l) = zh(ig, l)*(1.+0.61*zo(ig,l))
4630 END DO
4631 END DO
4632
4633 ! print*,'1 OK convect8'
4634 ! --------------------
4635
4636
4637 ! + + + + + + + + + + +
4638
4639
4640 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
4641 ! wh,wt,wo ...
4642
4643 ! + + + + + + + + + + + zh,zu,zv,zo,rho
4644
4645
4646 ! -------------------- zlev(1)
4647 ! \\\\\\\\\\\\\\\\\\\
4648
4649
4650 ! -----------------------------------------------------------------------
4651 ! Calcul des altitudes des couches
4652 ! -----------------------------------------------------------------------
4653
4654 DO l = 2, nlay
4655 DO ig = 1, ngrid
4656 zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
4657 END DO
4658 END DO
4659 DO ig = 1, ngrid
4660 zlev(ig, 1) = 0.
4661 zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
4662 END DO
4663 DO l = 1, nlay
4664 DO ig = 1, ngrid
4665 zlay(ig, l) = pphi(ig, l)/rg
4666 END DO
4667 END DO
4668
4669 ! print*,'2 OK convect8'
4670 ! -----------------------------------------------------------------------
4671 ! Calcul des densites
4672 ! -----------------------------------------------------------------------
4673
4674 DO l = 1, nlay
4675 DO ig = 1, ngrid
4676 rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*zh(ig,l))
4677 END DO
4678 END DO
4679
4680 DO l = 2, nlay
4681 DO ig = 1, ngrid
4682 rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
4683 END DO
4684 END DO
4685
4686 DO k = 1, nlay
4687 DO l = 1, nlay + 1
4688 DO ig = 1, ngrid
4689 wa(ig, k, l) = 0.
4690 END DO
4691 END DO
4692 END DO
4693
4694 ! print*,'3 OK convect8'
4695 ! ------------------------------------------------------------------
4696 ! Calcul de w2, quarre de w a partir de la cape
4697 ! a partir de w2, on calcule wa, vitesse de l'ascendance
4698
4699 ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
4700 ! w2 est stoke dans wa
4701
4702 ! ATTENTION: dans convect8, on n'utilise le calcule des wa
4703 ! independants par couches que pour calculer l'entrainement
4704 ! a la base et la hauteur max de l'ascendance.
4705
4706 ! Indicages:
4707 ! l'ascendance provenant du niveau k traverse l'interface l avec
4708 ! une vitesse wa(k,l).
4709
4710 ! --------------------
4711
4712 ! + + + + + + + + + +
4713
4714 ! wa(k,l) ---- -------------------- l
4715 ! / ! /||\ + + + + + + + + + +
4716 ! ||
4717 ! || --------------------
4718 ! ||
4719 ! || + + + + + + + + + +
4720 ! ||
4721 ! || --------------------
4722 ! ||__
4723 ! |___ + + + + + + + + + + k
4724
4725 ! --------------------
4726
4727
4728
4729 ! ------------------------------------------------------------------
4730
4731 ! CR: ponderation entrainement des couches instables
4732 ! def des entr_star tels que entr=f*entr_star
4733 DO l = 1, klev
4734 DO ig = 1, ngrid
4735 entr_star(ig, l) = 0.
4736 END DO
4737 END DO
4738 ! determination de la longueur de la couche d entrainement
4739 DO ig = 1, ngrid
4740 lentr(ig) = 1
4741 END DO
4742
4743 ! on ne considere que les premieres couches instables
4744 DO k = nlay - 2, 1, -1
4745 DO ig = 1, ngrid
4746 IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN
4747 lentr(ig) = k
4748 END IF
4749 END DO
4750 END DO
4751
4752 ! determination du lmin: couche d ou provient le thermique
4753 DO ig = 1, ngrid
4754 lmin(ig) = 1
4755 END DO
4756 DO ig = 1, ngrid
4757 DO l = nlay, 2, -1
4758 IF (ztv(ig,l-1)>ztv(ig,l)) THEN
4759 lmin(ig) = l - 1
4760 END IF
4761 END DO
4762 END DO
4763
4764 ! definition de l'entrainement des couches
4765 DO l = 1, klev - 1
4766 DO ig = 1, ngrid
4767 IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<=lentr(ig)) THEN
4768 entr_star(ig, l) = (ztv(ig,l)-ztv(ig,l+1))** & ! s
4769 ! (zlev(ig,l+1)-zlev(ig,l))
4770 sqrt(zlev(ig,l+1))
4771 END IF
4772 END DO
4773 END DO
4774 ! pas de thermique si couche 1 stable
4775 DO ig = 1, ngrid
4776 IF (lmin(ig)>1) THEN
4777 DO l = 1, klev
4778 entr_star(ig, l) = 0.
4779 END DO
4780 END IF
4781 END DO
4782 ! calcul de l entrainement total
4783 DO ig = 1, ngrid
4784 entr_star_tot(ig) = 0.
4785 END DO
4786 DO ig = 1, ngrid
4787 DO k = 1, klev
4788 entr_star_tot(ig) = entr_star_tot(ig) + entr_star(ig, k)
4789 END DO
4790 END DO
4791
4792 ! print*,'fin calcul entr_star'
4793 DO k = 1, klev
4794 DO ig = 1, ngrid
4795 ztva(ig, k) = ztv(ig, k)
4796 END DO
4797 END DO
4798 ! RC
4799 ! print*,'7 OK convect8'
4800 DO k = 1, klev + 1
4801 DO ig = 1, ngrid
4802 zw2(ig, k) = 0.
4803 fmc(ig, k) = 0.
4804 ! CR
4805 f_star(ig, k) = 0.
4806 ! RC
4807 larg_cons(ig, k) = 0.
4808 larg_detr(ig, k) = 0.
4809 wa_moy(ig, k) = 0.
4810 END DO
4811 END DO
4812
4813 ! print*,'8 OK convect8'
4814 DO ig = 1, ngrid
4815 linter(ig) = 1.
4816 lmaxa(ig) = 1
4817 lmix(ig) = 1
4818 wmaxa(ig) = 0.
4819 END DO
4820
4821 ! CR:
4822 DO l = 1, nlay - 2
4823 DO ig = 1, ngrid
4824 IF (ztv(ig,l)>ztv(ig,l+1) .AND. entr_star(ig,l)>1.E-10 .AND. &
4825 zw2(ig,l)<1E-10) THEN
4826 f_star(ig, l+1) = entr_star(ig, l)
4827 ! test:calcul de dteta
4828 zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
4829 (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
4830 larg_detr(ig, l) = 0.
4831 ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+entr_star(ig, &
4832 l)>1.E-10)) THEN
4833 f_star(ig, l+1) = f_star(ig, l) + entr_star(ig, l)
4834 ztva(ig, l) = (f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)*ztv(ig,l))/ &
4835 f_star(ig, l+1)
4836 zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/f_star(ig,l+1))**2 + &
4837 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
4838 END IF
4839 ! determination de zmax continu par interpolation lineaire
4840 IF (zw2(ig,l+1)<0.) THEN
4841 ! test
4842 IF (abs(zw2(ig,l+1)-zw2(ig,l))<1E-10) THEN
4843 ! print*,'pb linter'
4844 END IF
4845 linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
4846 ig,l))
4847 zw2(ig, l+1) = 0.
4848 lmaxa(ig) = l
4849 ELSE
4850 IF (zw2(ig,l+1)<0.) THEN
4851 ! print*,'pb1 zw2<0'
4852 END IF
4853 wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
4854 END IF
4855 IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
4856 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
4857 lmix(ig) = l + 1
4858 wmaxa(ig) = wa_moy(ig, l+1)
4859 END IF
4860 END DO
4861 END DO
4862 ! print*,'fin calcul zw2'
4863
4864 ! Calcul de la couche correspondant a la hauteur du thermique
4865 DO ig = 1, ngrid
4866 lmax(ig) = lentr(ig)
4867 END DO
4868 DO ig = 1, ngrid
4869 DO l = nlay, lentr(ig) + 1, -1
4870 IF (zw2(ig,l)<=1.E-10) THEN
4871 lmax(ig) = l - 1
4872 END IF
4873 END DO
4874 END DO
4875 ! pas de thermique si couche 1 stable
4876 DO ig = 1, ngrid
4877 IF (lmin(ig)>1) THEN
4878 lmax(ig) = 1
4879 lmin(ig) = 1
4880 END IF
4881 END DO
4882
4883 ! Determination de zw2 max
4884 DO ig = 1, ngrid
4885 wmax(ig) = 0.
4886 END DO
4887
4888 DO l = 1, nlay
4889 DO ig = 1, ngrid
4890 IF (l<=lmax(ig)) THEN
4891 IF (zw2(ig,l)<0.) THEN
4892 ! print*,'pb2 zw2<0'
4893 END IF
4894 zw2(ig, l) = sqrt(zw2(ig,l))
4895 wmax(ig) = max(wmax(ig), zw2(ig,l))
4896 ELSE
4897 zw2(ig, l) = 0.
4898 END IF
4899 END DO
4900 END DO
4901
4902 ! Longueur caracteristique correspondant a la hauteur des thermiques.
4903 DO ig = 1, ngrid
4904 zmax(ig) = 0.
4905 zlevinter(ig) = zlev(ig, 1)
4906 END DO
4907 DO ig = 1, ngrid
4908 ! calcul de zlevinter
4909 zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
4910 zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
4911 zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,lmin(ig)))
4912 END DO
4913
4914 ! print*,'avant fermeture'
4915 ! Fermeture,determination de f
4916 DO ig = 1, ngrid
4917 entr_star2(ig) = 0.
4918 END DO
4919 DO ig = 1, ngrid
4920 IF (entr_star_tot(ig)<1.E-10) THEN
4921 f(ig) = 0.
4922 ELSE
4923 DO k = lmin(ig), lentr(ig)
4924 entr_star2(ig) = entr_star2(ig) + entr_star(ig, k)**2/(rho(ig,k)*( &
4925 zlev(ig,k+1)-zlev(ig,k)))
4926 END DO
4927 ! Nouvelle fermeture
4928 f(ig) = wmax(ig)/(max(500.,zmax(ig))*r_aspect*entr_star2(ig))* &
4929 entr_star_tot(ig)
4930 ! test
4931 ! if (first) then
4932 ! f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
4933 ! s *wmax(ig))
4934 ! endif
4935 END IF
4936 ! f0(ig)=f(ig)
4937 ! first=.true.
4938 END DO
4939 ! print*,'apres fermeture'
4940
4941 ! Calcul de l'entrainement
4942 DO k = 1, klev
4943 DO ig = 1, ngrid
4944 entr(ig, k) = f(ig)*entr_star(ig, k)
4945 END DO
4946 END DO
4947 ! CR:test pour entrainer moins que la masse
4948 DO ig = 1, ngrid
4949 DO l = 1, lentr(ig)
4950 IF ((entr(ig,l)*ptimestep)>(0.9*masse(ig,l))) THEN
4951 entr(ig, l+1) = entr(ig, l+1) + entr(ig, l) - &
4952 0.9*masse(ig, l)/ptimestep
4953 entr(ig, l) = 0.9*masse(ig, l)/ptimestep
4954 END IF
4955 END DO
4956 END DO
4957 ! CR: fin test
4958 ! Calcul des flux
4959 DO ig = 1, ngrid
4960 DO l = 1, lmax(ig) - 1
4961 fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
4962 END DO
4963 END DO
4964
4965 ! RC
4966
4967
4968 ! print*,'9 OK convect8'
4969 ! print*,'WA1 ',wa_moy
4970
4971 ! determination de l'indice du debut de la mixed layer ou w decroit
4972
4973 ! calcul de la largeur de chaque ascendance dans le cas conservatif.
4974 ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
4975 ! d'une couche est �gale � la hauteur de la couche alimentante.
4976 ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
4977 ! de la vitesse d'entrainement horizontal dans la couche alimentante.
4978
4979 DO l = 2, nlay
4980 DO ig = 1, ngrid
4981 IF (l<=lmaxa(ig)) THEN
4982 zw = max(wa_moy(ig,l), 1.E-10)
4983 larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw)
4984 END IF
4985 END DO
4986 END DO
4987
4988 DO l = 2, nlay
4989 DO ig = 1, ngrid
4990 IF (l<=lmaxa(ig)) THEN
4991 ! if (idetr.eq.0) then
4992 ! cette option est finalement en dur.
4993 IF ((l_mix*zlev(ig,l))<0.) THEN
4994 ! print*,'pb l_mix*zlev<0'
4995 END IF
4996 ! CR: test: nouvelle def de lambda
4997 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
4998 IF (zw2(ig,l)>1.E-10) THEN
4999 larg_detr(ig, l) = sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
5000 ELSE
5001 larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
5002 END IF
5003 ! RC
5004 ! else if (idetr.eq.1) then
5005 ! larg_detr(ig,l)=larg_cons(ig,l)
5006 ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
5007 ! else if (idetr.eq.2) then
5008 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5009 ! s *sqrt(wa_moy(ig,l))
5010 ! else if (idetr.eq.4) then
5011 ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
5012 ! s *wa_moy(ig,l)
5013 ! endif
5014 END IF
5015 END DO
5016 END DO
5017
5018 ! print*,'10 OK convect8'
5019 ! print*,'WA2 ',wa_moy
5020 ! calcul de la fraction de la maille concern�e par l'ascendance en tenant
5021 ! compte de l'epluchage du thermique.
5022
5023 ! CR def de zmix continu (profil parabolique des vitesses)
5024 DO ig = 1, ngrid
5025 IF (lmix(ig)>1.) THEN
5026 ! test
5027 IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
5028 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
5029 zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- &
5030 (zlev(ig,lmix(ig)))))>1E-10) THEN
5031
5032 zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) &
5033 )**2-(zlev(ig,lmix(ig)+1))**2)-(zw2(ig,lmix(ig))-zw2(ig, &
5034 lmix(ig)+1))*((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))/ &
5035 (2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
5036 (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
5037 zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
5038 ELSE
5039 zmix(ig) = zlev(ig, lmix(ig))
5040 ! print*,'pb zmix'
5041 END IF
5042 ELSE
5043 zmix(ig) = 0.
5044 END IF
5045 ! test
5046 IF ((zmax(ig)-zmix(ig))<0.) THEN
5047 zmix(ig) = 0.99*zmax(ig)
5048 ! print*,'pb zmix>zmax'
5049 END IF
5050 END DO
5051
5052 ! calcul du nouveau lmix correspondant
5053 DO ig = 1, ngrid
5054 DO l = 1, klev
5055 IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
5056 lmix(ig) = l
5057 END IF
5058 END DO
5059 END DO
5060
5061 DO l = 2, nlay
5062 DO ig = 1, ngrid
5063 IF (larg_cons(ig,l)>1.) THEN
5064 ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
5065 fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
5066 ! test
5067 fraca(ig, l) = max(fraca(ig,l), 0.)
5068 fraca(ig, l) = min(fraca(ig,l), 0.5)
5069 fracd(ig, l) = 1. - fraca(ig, l)
5070 fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
5071 ELSE
5072 ! wa_moy(ig,l)=0.
5073 fraca(ig, l) = 0.
5074 fracc(ig, l) = 0.
5075 fracd(ig, l) = 1.
5076 END IF
5077 END DO
5078 END DO
5079 ! CR: calcul de fracazmix
5080 DO ig = 1, ngrid
5081 fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
5082 (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
5083 fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca(ig &
5084 ,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
5085 END DO
5086
5087 DO l = 2, nlay
5088 DO ig = 1, ngrid
5089 IF (larg_cons(ig,l)>1.) THEN
5090 IF (l>lmix(ig)) THEN
5091 ! test
5092 IF (zmax(ig)-zmix(ig)<1.E-10) THEN
5093 ! print*,'pb xxx'
5094 xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
5095 ELSE
5096 xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
5097 END IF
5098 IF (idetr==0) THEN
5099 fraca(ig, l) = fracazmix(ig)
5100 ELSE IF (idetr==1) THEN
5101 fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
5102 ELSE IF (idetr==2) THEN
5103 fraca(ig, l) = fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
5104 ELSE
5105 fraca(ig, l) = fracazmix(ig)*xxx(ig, l)**2
5106 END IF
5107 ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
5108 fraca(ig, l) = max(fraca(ig,l), 0.)
5109 fraca(ig, l) = min(fraca(ig,l), 0.5)
5110 fracd(ig, l) = 1. - fraca(ig, l)
5111 fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
5112 END IF
5113 END IF
5114 END DO
5115 END DO
5116
5117 ! print*,'fin calcul fraca'
5118 ! print*,'11 OK convect8'
5119 ! print*,'Ea3 ',wa_moy
5120 ! ------------------------------------------------------------------
5121 ! Calcul de fracd, wd
5122 ! somme wa - wd = 0
5123 ! ------------------------------------------------------------------
5124
5125
5126 DO ig = 1, ngrid
5127 fm(ig, 1) = 0.
5128 fm(ig, nlay+1) = 0.
5129 END DO
5130
5131 DO l = 2, nlay
5132 DO ig = 1, ngrid
5133 fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
5134 ! CR:test
5135 IF (entr(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) THEN
5136 fm(ig, l) = fm(ig, l-1)
5137 ! write(1,*)'ajustement fm, l',l
5138 END IF
5139 ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
5140 ! RC
5141 END DO
5142 DO ig = 1, ngrid
5143 IF (fracd(ig,l)<0.1) THEN
5144 abort_message = 'fracd trop petit'
5145 CALL abort_physic(modname, abort_message, 1)
5146 ELSE
5147 ! vitesse descendante "diagnostique"
5148 wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
5149 END IF
5150 END DO
5151 END DO
5152
5153 DO l = 1, nlay
5154 DO ig = 1, ngrid
5155 ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
5156 masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
5157 END DO
5158 END DO
5159
5160 ! print*,'12 OK convect8'
5161 ! print*,'WA4 ',wa_moy
5162 ! c------------------------------------------------------------------
5163 ! calcul du transport vertical
5164 ! ------------------------------------------------------------------
5165
5166 GO TO 4444
5167 ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
5168 DO l = 2, nlay - 1
5169 DO ig = 1, ngrid
5170 IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( &
5171 ig,l+1)) THEN
5172 ! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
5173 ! s ,fm(ig,l+1)*ptimestep
5174 ! s ,' M=',masse(ig,l),masse(ig,l+1)
5175 END IF
5176 END DO
5177 END DO
5178
5179 DO l = 1, nlay
5180 DO ig = 1, ngrid
5181 IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
5182 ! print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
5183 ! s ,entr(ig,l)*ptimestep
5184 ! s ,' M=',masse(ig,l)
5185 END IF
5186 END DO
5187 END DO
5188
5189 DO l = 1, nlay
5190 DO ig = 1, ngrid
5191 IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
5192 ! print*,'WARN!!! fm exagere ig=',ig,' l=',l
5193 ! s ,' FM=',fm(ig,l)
5194 END IF
5195 IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
5196 ! print*,'WARN!!! masse exagere ig=',ig,' l=',l
5197 ! s ,' M=',masse(ig,l)
5198 ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
5199 ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
5200 ! print*,'zlev(ig,l+1),zlev(ig,l)'
5201 ! s ,zlev(ig,l+1),zlev(ig,l)
5202 ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
5203 ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
5204 END IF
5205 IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
5206 ! print*,'WARN!!! entr exagere ig=',ig,' l=',l
5207 ! s ,' E=',entr(ig,l)
5208 END IF
5209 END DO
5210 END DO
5211
5212 4444 CONTINUE
5213
5214 ! CR:redefinition du entr
5215 DO l = 1, nlay
5216 DO ig = 1, ngrid
5217 detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1)
5218 IF (detr(ig,l)<0.) THEN
5219 entr(ig, l) = entr(ig, l) - detr(ig, l)
5220 detr(ig, l) = 0.
5221 ! print*,'WARNING !!! detrainement negatif ',ig,l
5222 END IF
5223 END DO
5224 END DO
5225 ! RC
5226 IF (w2di==1) THEN
5227 fm0 = fm0 + ptimestep*(fm-fm0)/tho
5228 entr0 = entr0 + ptimestep*(entr-entr0)/tho
5229 ELSE
5230 fm0 = fm
5231 entr0 = entr
5232 END IF
5233
5234 IF (1==1) THEN
5235 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, &
5236 zha)
5237 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, &
5238 zoa)
5239 ELSE
5240 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
5241 zdhadj, zha)
5242 CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
5243 pdoadj, zoa)
5244 END IF
5245
5246 IF (1==0) THEN
5247 CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
5248 zu, zv, pduadj, pdvadj, zua, zva)
5249 ELSE
5250 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
5251 zua)
5252 CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
5253 zva)
5254 END IF
5255
5256 DO l = 1, nlay
5257 DO ig = 1, ngrid
5258 zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
5259 zf2 = zf/(1.-zf)
5260 thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l))**2
5261 wth2(ig, l) = zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
5262 END DO
5263 END DO
5264
5265
5266
5267 ! print*,'13 OK convect8'
5268 ! print*,'WA5 ',wa_moy
5269 DO l = 1, nlay
5270 DO ig = 1, ngrid
5271 pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l)
5272 END DO
5273 END DO
5274
5275
5276 ! do l=1,nlay
5277 ! do ig=1,ngrid
5278 ! if(abs(pdtadj(ig,l))*86400..gt.500.) then
5279 ! print*,'WARN!!! ig=',ig,' l=',l
5280 ! s ,' pdtadj=',pdtadj(ig,l)
5281 ! endif
5282 ! if(abs(pdoadj(ig,l))*86400..gt.1.) then
5283 ! print*,'WARN!!! ig=',ig,' l=',l
5284 ! s ,' pdoadj=',pdoadj(ig,l)
5285 ! endif
5286 ! enddo
5287 ! enddo
5288
5289 ! print*,'14 OK convect8'
5290 ! ------------------------------------------------------------------
5291 ! Calculs pour les sorties
5292 ! ------------------------------------------------------------------
5293
5294 RETURN
5295 END SUBROUTINE thermcell_sec
5296
5297