My Project
 All Classes Files Functions Variables Macros
thermcell_plume.F90
Go to the documentation of this file.
1 !
2 ! $Id: thermcell_plume.F90 1503 2011-03-23 11:57:52Z idelkadi $
3 !
4  SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
5  & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
6  & lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
7  & ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
8  & ,lev_out,lunout1,igout)
9 
10 !--------------------------------------------------------------------------
11 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
12 !--------------------------------------------------------------------------
13 
14  IMPLICIT NONE
15 
16 #include "YOMCST.h"
17 #include "YOETHF.h"
18 #include "FCTTRE.h"
19 #include "iniprint.h"
20 #include "thermcell.h"
21 
22  INTEGER itap
23  INTEGER lunout1,igout
24  INTEGER ngrid,klev
25  REAL ptimestep
26  REAL ztv(ngrid,klev)
27  REAL zthl(ngrid,klev)
28  REAL po(ngrid,klev)
29  REAL zl(ngrid,klev)
30  REAL rhobarz(ngrid,klev)
31  REAL zlev(ngrid,klev+1)
32  REAL pplev(ngrid,klev+1)
33  REAL pphi(ngrid,klev)
34  REAL zpspsk(ngrid,klev)
35  REAL alim_star(ngrid,klev)
36  REAL f0(ngrid)
37  INTEGER lalim(ngrid)
38  integer lev_out ! niveau pour les print
39  integer nbpb
40  real zcon2(ngrid)
41 
42  real alim_star_tot(ngrid)
43 
44  REAL ztva(ngrid,klev)
45  REAL ztla(ngrid,klev)
46  REAL zqla(ngrid,klev)
47  REAL zqta(ngrid,klev)
48  REAL zha(ngrid,klev)
49 
50  REAL detr_star(ngrid,klev)
51  REAL coefc
52  REAL entr_star(ngrid,klev)
53  REAL detr(ngrid,klev)
54  REAL entr(ngrid,klev)
55 
56  REAL csc(ngrid,klev)
57 
58  REAL zw2(ngrid,klev+1)
59  REAL w_est(ngrid,klev+1)
60  REAL f_star(ngrid,klev+1)
61  REAL wa_moy(ngrid,klev+1)
62 
63  REAL ztva_est(ngrid,klev)
64  REAL zqla_est(ngrid,klev)
65  REAL zqsatth(ngrid,klev)
66  REAL zta_est(ngrid,klev)
67  REAL zdw2
68  REAL zw2modif
69  REAL zeps
70 
71  REAL linter(ngrid)
72  INTEGER lmix(ngrid)
73  INTEGER lmix_bis(ngrid)
74  REAL wmaxa(ngrid)
75 
76  INTEGER ig,l,k
77 
78  real zdz,zfact,zbuoy,zalpha,zdrag
79  real zcor,zdelta,zcvm5,qlbef
80  real tbef,qsatbef
81  real dqsat_dt,dt,num,denom
82  REAL reps,rlvcp,ddt0
83  parameter(ddt0=.01)
84  logical zsat
85  LOGICAL active(ngrid),activetmp(ngrid)
86  REAL fact_gamma,fact_epsilon,fact_gamma2
87  REAL c2(ngrid,klev)
88  REAL a1,m
89 
90  REAL zw2fact,expa
91  zsat=.false.
92 ! Initialisation
93  rlvcp = rlvtt/rcpd
94 
95 
96  fact_epsilon=0.002
97  a1=2./3.
98  fact_gamma=0.9
99  zfact=fact_gamma/(1+fact_gamma)
100  fact_gamma2=zfact
101  expa=0.
102 
103 
104 ! Initialisations des variables reeles
105 if (1==1) then
106  ztva(:,:)=ztv(:,:)
107  ztva_est(:,:)=ztva(:,:)
108  ztla(:,:)=zthl(:,:)
109  zqta(:,:)=po(:,:)
110  zha(:,:) = ztva(:,:)
111 else
112  ztva(:,:)=0.
113  ztva_est(:,:)=0.
114  ztla(:,:)=0.
115  zqta(:,:)=0.
116  zha(:,:) =0.
117 endif
118 
119  zqla_est(:,:)=0.
120  zqsatth(:,:)=0.
121  zqla(:,:)=0.
122  detr_star(:,:)=0.
123  entr_star(:,:)=0.
124  alim_star(:,:)=0.
125  alim_star_tot(:)=0.
126  csc(:,:)=0.
127  detr(:,:)=0.
128  entr(:,:)=0.
129  zw2(:,:)=0.
130  w_est(:,:)=0.
131  f_star(:,:)=0.
132  wa_moy(:,:)=0.
133  linter(:)=1.
134  linter(:)=1.
135 
136 ! Initialisation des variables entieres
137  lmix(:)=1
138  lmix_bis(:)=2
139  wmaxa(:)=0.
140  lalim(:)=1
141 
142 !-------------------------------------------------------------------------
143 ! On ne considere comme actif que les colonnes dont les deux premieres
144 ! couches sont instables.
145 !-------------------------------------------------------------------------
146  active(:)=ztv(:,1)>ztv(:,2)
147 
148 !-------------------------------------------------------------------------
149 ! Definition de l'alimentation a l'origine dans thermcell_init
150 !-------------------------------------------------------------------------
151  do l=1,klev-1
152  do ig=1,ngrid
153  if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
154  alim_star(ig,l)=max((ztv(ig,l)-ztv(ig,l+1)),0.) &
155  & *sqrt(zlev(ig,l+1))
156  lalim(ig)=l+1
157  alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
158  endif
159  enddo
160  enddo
161  do l=1,klev
162  do ig=1,ngrid
163  if (alim_star_tot(ig) > 1.e-10 ) then
164  alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
165  endif
166  enddo
167  enddo
168  alim_star_tot(:)=1.
169 
170 
171 !------------------------------------------------------------------------------
172 ! Calcul dans la premiere couche
173 ! On decide dans cette version que le thermique n'est actif que si la premiere
174 ! couche est instable.
175 ! Pourrait etre change si on veut que le thermiques puisse se déclencher
176 ! dans une couche l>1
177 !------------------------------------------------------------------------------
178 do ig=1,ngrid
179 ! Le panache va prendre au debut les caracteristiques de l'air contenu
180 ! dans cette couche.
181  if (active(ig)) then
182  ztla(ig,1)=zthl(ig,1)
183  zqta(ig,1)=po(ig,1)
184  zqla(ig,1)=zl(ig,1)
185 !cr: attention, prise en compte de f*(1)=1
186  f_star(ig,2)=alim_star(ig,1)
187  zw2(ig,2)=2.*rg*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
188 & *(zlev(ig,2)-zlev(ig,1)) &
189 & *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
190  w_est(ig,2)=zw2(ig,2)
191  endif
192 enddo
193 !
194 
195 !==============================================================================
196 !boucle de calcul de la vitesse verticale dans le thermique
197 !==============================================================================
198 do l=2,klev-1
199 !==============================================================================
200 
201 
202 ! On decide si le thermique est encore actif ou non
203 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
204  do ig=1,ngrid
205  active(ig)=active(ig) &
206 & .and. zw2(ig,l)>1.e-10 &
207 & .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
208  enddo
209 
210 
211 
212 ! Premier calcul de la vitesse verticale a partir de la temperature
213 ! potentielle virtuelle
214 ! if (1.eq.1) then
215 ! w_est(ig,3)=zw2(ig,2)* &
216 ! & ((f_star(ig,2))**2) &
217 ! & /(f_star(ig,2)+alim_star(ig,2))**2+ &
218 ! & 2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
219 ! & *(zlev(ig,3)-zlev(ig,2))
220 ! endif
221 
222 
223 !---------------------------------------------------------------------------
224 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
225 ! sans tenir compte du detrainement et de l'entrainement dans cette
226 ! couche
227 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
228 ! avant) a l'alimentation pour avoir un calcul plus propre
229 !---------------------------------------------------------------------------
230 
231  call thermcell_condens(ngrid,active, &
232 & zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
233 
234  do ig=1,ngrid
235  if(active(ig)) then
236  ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+rlvcp*zqla_est(ig,l)
237  zta_est(ig,l)=ztva_est(ig,l)
238  ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
239  ztva_est(ig,l) = ztva_est(ig,l)*(1.+retv*(zqta(ig,l-1) &
240  & -zqla_est(ig,l))-zqla_est(ig,l))
241 
242  if (1.eq.0) then
243 !calcul de w_est sans prendre en compte le drag
244  w_est(ig,l+1)=zw2(ig,l)* &
245  & ((f_star(ig,l))**2) &
246  & /(f_star(ig,l)+alim_star(ig,l))**2+ &
247  & 2.*rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
248  & *(zlev(ig,l+1)-zlev(ig,l))
249  else
250 
251  zdz=zlev(ig,l+1)-zlev(ig,l)
252  zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
253  zbuoy=rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
254  zdrag=fact_epsilon/(zalpha**expa)
255  zw2fact=zbuoy/zdrag*a1
256  w_est(ig,l+1)=(w_est(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) &
257  & +zw2fact
258 
259  endif
260 
261  if (w_est(ig,l+1).lt.0.) then
262  w_est(ig,l+1)=zw2(ig,l)
263  endif
264  endif
265  enddo
266 
267 !-------------------------------------------------
268 !calcul des taux d'entrainement et de detrainement
269 !-------------------------------------------------
270 
271  do ig=1,ngrid
272  if (active(ig)) then
273  zdz=zlev(ig,l+1)-zlev(ig,l)
274  zbuoy=rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
275 
276 ! estimation de la fraction couverte par les thermiques
277  zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
278 
279 !calcul de la soumission papier
280 ! Calcul du taux d'entrainement entr_star (epsilon)
281  entr_star(ig,l)=f_star(ig,l)*zdz * ( zfact * max(0., &
282  & a1*zbuoy/w_est(ig,l+1) &
283  & - fact_epsilon/zalpha**expa ) &
284  & +0. )
285 
286 !calcul du taux de detrainment (delta)
287 ! detr_star(ig,l)=f_star(ig,l)*zdz * ( &
288 ! & MAX(1.e-3, &
289 ! & -fact_gamma2*a1*zbuoy/w_est(ig,l+1) &
290 ! & +0.01*(max(zqta(ig,l-1)-po(ig,l),0.)/(po(ig,l))/(w_est(ig,l+1)))**0.5 &
291 ! & +0. ))
292 
293  m=0.5
294 
295  detr_star(ig,l)=1.*f_star(ig,l)*zdz * &
296  & max(5.e-4,-fact_gamma2*a1*(1./w_est(ig,l+1))*((1.-(1.-m)/(1.+70*zqta(ig,l-1)))*zbuoy &
297  & -40*(1.-m)*(max(zqta(ig,l-1)-po(ig,l),0.))/(1.+70*zqta(ig,l-1)) ) )
298 
299 ! detr_star(ig,l)=f_star(ig,l)*zdz * ( &
300 ! & MAX(0.0, &
301 ! & -fact_gamma2*a1*zbuoy/w_est(ig,l+1) &
302 ! & +20*(max(zqta(ig,l-1)-po(ig,l),0.))**1*(zalpha/w_est(ig,l+1))**0.5 &
303 ! & +0. ))
304 
305 
306 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
307 ! alim_star et 0 sinon
308  if (l.lt.lalim(ig)) then
309  alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
310  entr_star(ig,l)=0.
311  endif
312 
313 !attention test
314 ! if (detr_star(ig,l).gt.(f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l))) then
315 ! detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
316 ! endif
317 ! Calcul du flux montant normalise
318  f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
319  & -detr_star(ig,l)
320 
321  endif
322  enddo
323 
324 !----------------------------------------------------------------------------
325 !calcul de la vitesse verticale en melangeant Tl et qt du thermique
326 !---------------------------------------------------------------------------
327  activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
328  do ig=1,ngrid
329  if (activetmp(ig)) then
330  zsat=.false.
331  ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
332  & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
333  & /(f_star(ig,l+1)+detr_star(ig,l))
334  zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
335  & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
336  & /(f_star(ig,l+1)+detr_star(ig,l))
337 
338  endif
339  enddo
340 
341  call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
342 
343 
344  do ig=1,ngrid
345  if (activetmp(ig)) then
346 ! on ecrit de maniere conservative (sat ou non)
347 ! T = Tl +Lv/Cp ql
348  ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+rlvcp*zqla(ig,l)
349  ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
350 !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
351  zha(ig,l) = ztva(ig,l)
352  ztva(ig,l) = ztva(ig,l)*(1.+retv*(zqta(ig,l) &
353  & -zqla(ig,l))-zqla(ig,l))
354 
355 !on ecrit zqsat
356  zqsatth(ig,l)=qsatbef
357 
358 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
359 ! zw2(ig,l+1)=&
360 ! & zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*(zlev(ig,l+1)-zlev(ig,l))) &
361 ! & +2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
362 ! & *1./(1.+fact_gamma) &
363 ! & *(zlev(ig,l+1)-zlev(ig,l))
364 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
365 ! La meme en plus modulaire :
366  zbuoy=rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
367  zdz=zlev(ig,l+1)-zlev(ig,l)
368 
369 
370  zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
371 
372 !if (1==0) then
373 ! zw2modif=zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*zdz)
374 ! zdw2=2.*zbuoy/(1.+fact_gamma)*zdz
375 ! zw2(ig,l+1)=zw2modif+zdw2
376 !else
377  zdrag=fact_epsilon/(zalpha**expa)
378  zw2fact=zbuoy/zdrag*a1
379  zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) &
380  & +zw2fact
381 !endif
382 
383  endif
384  enddo
385 
386  if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
387 !
388 !---------------------------------------------------------------------------
389 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
390 !---------------------------------------------------------------------------
391 
392  nbpb=0
393  do ig=1,ngrid
394  if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
395 ! stop'On tombe sur le cas particulier de thermcell_dry'
396 ! print*,'On tombe sur le cas particulier de thermcell_plume'
397  nbpb=nbpb+1
398  zw2(ig,l+1)=0.
399  linter(ig)=l+1
400  endif
401 
402  if (zw2(ig,l+1).lt.0.) then
403  linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
404  & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
405  zw2(ig,l+1)=0.
406  endif
407 
408  wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
409 
410  if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
411 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
412 !on rajoute le calcul de lmix_bis
413  if (zqla(ig,l).lt.1.e-10) then
414  lmix_bis(ig)=l+1
415  endif
416  lmix(ig)=l+1
417  wmaxa(ig)=wa_moy(ig,l+1)
418  endif
419  enddo
420 
421  if (nbpb>0) then
422  print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
423  endif
424 
425 !=========================================================================
426 ! FIN DE LA BOUCLE VERTICALE
427  enddo
428 !=========================================================================
429 
430 !on recalcule alim_star_tot
431  do ig=1,ngrid
432  alim_star_tot(ig)=0.
433  enddo
434  do ig=1,ngrid
435  do l=1,lalim(ig)-1
436  alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
437  enddo
438  enddo
439 
440 
441  if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
442 
443 
444  return
445  end
446 
447 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
448 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
449 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
450 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
451 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
452 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
453  SUBROUTINE thermcellv1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
454 & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
455 & lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
456 & ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
457 & ,lev_out,lunout1,igout)
458 
459 !--------------------------------------------------------------------------
460 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
461 ! Version conforme a l'article de Rio et al. 2010.
462 ! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
463 !--------------------------------------------------------------------------
464 
465  IMPLICIT NONE
466 
467 #include "YOMCST.h"
468 #include "YOETHF.h"
469 #include "FCTTRE.h"
470 #include "iniprint.h"
471 #include "thermcell.h"
472 
473  INTEGER itap
474  INTEGER lunout1,igout
475  INTEGER ngrid,klev
476  REAL ptimestep
477  REAL ztv(ngrid,klev)
478  REAL zthl(ngrid,klev)
479  REAL po(ngrid,klev)
480  REAL zl(ngrid,klev)
481  REAL rhobarz(ngrid,klev)
482  REAL zlev(ngrid,klev+1)
483  REAL pplev(ngrid,klev+1)
484  REAL pphi(ngrid,klev)
485  REAL zpspsk(ngrid,klev)
486  REAL alim_star(ngrid,klev)
487  REAL f0(ngrid)
488  INTEGER lalim(ngrid)
489  integer lev_out ! niveau pour les print
490  integer nbpb
491 
492  real alim_star_tot(ngrid)
493 
494  REAL ztva(ngrid,klev)
495  REAL ztla(ngrid,klev)
496  REAL zqla(ngrid,klev)
497  REAL zqta(ngrid,klev)
498  REAL zha(ngrid,klev)
499 
500  REAL detr_star(ngrid,klev)
501  REAL coefc
502  REAL entr_star(ngrid,klev)
503  REAL detr(ngrid,klev)
504  REAL entr(ngrid,klev)
505 
506  REAL csc(ngrid,klev)
507 
508  REAL zw2(ngrid,klev+1)
509  REAL w_est(ngrid,klev+1)
510  REAL f_star(ngrid,klev+1)
511  REAL wa_moy(ngrid,klev+1)
512 
513  REAL ztva_est(ngrid,klev)
514  REAL zqla_est(ngrid,klev)
515  REAL zqsatth(ngrid,klev)
516  REAL zta_est(ngrid,klev)
517  REAL ztemp(ngrid),zqsat(ngrid)
518  REAL zdw2
519  REAL zw2modif
520  REAL zw2fact
521  REAL zeps(ngrid,klev)
522 
523  REAL linter(ngrid)
524  INTEGER lmix(ngrid)
525  INTEGER lmix_bis(ngrid)
526  REAL wmaxa(ngrid)
527 
528  INTEGER ig,l,k
529 
530  real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
531  real zbuoybis
532  real zcor,zdelta,zcvm5,qlbef,zdz2
533  real betalpha,zbetalpha
534  real eps, afact
535  REAL reps,rlvcp,ddt0
536  parameter(ddt0=.01)
537  logical zsat
538  LOGICAL active(ngrid),activetmp(ngrid)
539  REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
540  REAL c2(ngrid,klev)
541  zsat=.false.
542 ! Initialisation
543 
544  rlvcp = rlvtt/rcpd
545  fact_epsilon=0.002
546  betalpha=0.9
547  afact=2./3.
548 
549  zbetalpha=betalpha/(1.+betalpha)
550 
551 
552 ! Initialisations des variables reeles
553 if (1==0) then
554  ztva(:,:)=ztv(:,:)
555  ztva_est(:,:)=ztva(:,:)
556  ztla(:,:)=zthl(:,:)
557  zqta(:,:)=po(:,:)
558  zha(:,:) = ztva(:,:)
559 else
560  ztva(:,:)=0.
561  ztva_est(:,:)=0.
562  ztla(:,:)=0.
563  zqta(:,:)=0.
564  zha(:,:) =0.
565 endif
566 
567  zqla_est(:,:)=0.
568  zqsatth(:,:)=0.
569  zqla(:,:)=0.
570  detr_star(:,:)=0.
571  entr_star(:,:)=0.
572  alim_star(:,:)=0.
573  alim_star_tot(:)=0.
574  csc(:,:)=0.
575  detr(:,:)=0.
576  entr(:,:)=0.
577  zw2(:,:)=0.
578  zbuoy(:,:)=0.
579  gamma(:,:)=0.
580  zeps(:,:)=0.
581  w_est(:,:)=0.
582  f_star(:,:)=0.
583  wa_moy(:,:)=0.
584  linter(:)=1.
585 ! linter(:)=1.
586 ! Initialisation des variables entieres
587  lmix(:)=1
588  lmix_bis(:)=2
589  wmaxa(:)=0.
590  lalim(:)=1
591 
592 
593 !-------------------------------------------------------------------------
594 ! On ne considere comme actif que les colonnes dont les deux premieres
595 ! couches sont instables.
596 !-------------------------------------------------------------------------
597  active(:)=ztv(:,1)>ztv(:,2)
598 
599 !-------------------------------------------------------------------------
600 ! Definition de l'alimentation a l'origine dans thermcell_init
601 !-------------------------------------------------------------------------
602  do l=1,klev-1
603  do ig=1,ngrid
604  if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
605  alim_star(ig,l)=max((ztv(ig,l)-ztv(ig,l+1)),0.) &
606  & *sqrt(zlev(ig,l+1))
607  lalim(ig)=l+1
608  alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
609  endif
610  enddo
611  enddo
612  do l=1,klev
613  do ig=1,ngrid
614  if (alim_star_tot(ig) > 1.e-10 ) then
615  alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
616  endif
617  enddo
618  enddo
619  alim_star_tot(:)=1.
620 
621 
622 
623 !------------------------------------------------------------------------------
624 ! Calcul dans la premiere couche
625 ! On decide dans cette version que le thermique n'est actif que si la premiere
626 ! couche est instable.
627 ! Pourrait etre change si on veut que le thermiques puisse se déclencher
628 ! dans une couche l>1
629 !------------------------------------------------------------------------------
630 do ig=1,ngrid
631 ! Le panache va prendre au debut les caracteristiques de l'air contenu
632 ! dans cette couche.
633  if (active(ig)) then
634  ztla(ig,1)=zthl(ig,1)
635  zqta(ig,1)=po(ig,1)
636  zqla(ig,1)=zl(ig,1)
637 !cr: attention, prise en compte de f*(1)=1
638  f_star(ig,2)=alim_star(ig,1)
639  zw2(ig,2)=2.*rg*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
640 & *(zlev(ig,2)-zlev(ig,1)) &
641 & *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
642  w_est(ig,2)=zw2(ig,2)
643  endif
644 enddo
645 !
646 
647 !==============================================================================
648 !boucle de calcul de la vitesse verticale dans le thermique
649 !==============================================================================
650 do l=2,klev-1
651 !==============================================================================
652 
653 
654 ! On decide si le thermique est encore actif ou non
655 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
656  do ig=1,ngrid
657  active(ig)=active(ig) &
658 & .and. zw2(ig,l)>1.e-10 &
659 & .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
660  enddo
661 
662 
663 
664 !---------------------------------------------------------------------------
665 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
666 ! sans tenir compte du detrainement et de l'entrainement dans cette
667 ! couche
668 ! C'est a dire qu'on suppose
669 ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
670 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
671 ! avant) a l'alimentation pour avoir un calcul plus propre
672 !---------------------------------------------------------------------------
673 
674  ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
675  call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
676 
677  do ig=1,ngrid
678 ! print*,'active',active(ig),ig,l
679  if(active(ig)) then
680  zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
681  ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+rlvcp*zqla_est(ig,l)
682  zta_est(ig,l)=ztva_est(ig,l)
683  ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
684  ztva_est(ig,l) = ztva_est(ig,l)*(1.+retv*(zqta(ig,l-1) &
685  & -zqla_est(ig,l))-zqla_est(ig,l))
686 
687 !------------------------------------------------
688 !AJAM:nouveau calcul de w²
689 !------------------------------------------------
690  zdz=zlev(ig,l+1)-zlev(ig,l)
691  zbuoy(ig,l)=rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
692 
693  zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
694  zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
695  w_est(ig,l+1)=max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
696 
697 
698  if (w_est(ig,l+1).lt.0.) then
699  w_est(ig,l+1)=zw2(ig,l)
700  endif
701  endif
702  enddo
703 
704 
705 !-------------------------------------------------
706 !calcul des taux d'entrainement et de detrainement
707 !-------------------------------------------------
708 
709  do ig=1,ngrid
710  if (active(ig)) then
711 
712  zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
713  zw2m=w_est(ig,l+1)
714  zdz=zlev(ig,l+1)-zlev(ig,l)
715  zbuoy(ig,l)=rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
716 ! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
717  zbuoybis=zbuoy(ig,l)
718  zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
719  zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
720 
721 
722  entr_star(ig,l)=f_star(ig,l)*zdz* zbetalpha*max(0., &
723  & afact*zbuoybis/zw2m - fact_epsilon )
724 
725 
726  detr_star(ig,l)=f_star(ig,l)*zdz &
727  & *max(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m &
728  & + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
729 
730 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
731 ! alim_star et 0 sinon
732  if (l.lt.lalim(ig)) then
733  alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
734  entr_star(ig,l)=0.
735  endif
736 
737 ! Calcul du flux montant normalise
738  f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
739  & -detr_star(ig,l)
740 
741  endif
742  enddo
743 
744 
745 !----------------------------------------------------------------------------
746 !calcul de la vitesse verticale en melangeant Tl et qt du thermique
747 !---------------------------------------------------------------------------
748  activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
749  do ig=1,ngrid
750  if (activetmp(ig)) then
751  zsat=.false.
752  ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
753  & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
754  & /(f_star(ig,l+1)+detr_star(ig,l))
755  zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
756  & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
757  & /(f_star(ig,l+1)+detr_star(ig,l))
758 
759  endif
760  enddo
761 
762  ztemp(:)=zpspsk(:,l)*ztla(:,l)
763  call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
764 
765  do ig=1,ngrid
766  if (activetmp(ig)) then
767 ! on ecrit de maniere conservative (sat ou non)
768 ! T = Tl +Lv/Cp ql
769  zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
770  ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+rlvcp*zqla(ig,l)
771  ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
772 !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
773  zha(ig,l) = ztva(ig,l)
774  ztva(ig,l) = ztva(ig,l)*(1.+retv*(zqta(ig,l) &
775  & -zqla(ig,l))-zqla(ig,l))
776  zbuoy(ig,l)=rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
777  zdz=zlev(ig,l+1)-zlev(ig,l)
778  zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
779 
780  zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
781  zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
782  zw2(ig,l+1)=max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
783  endif
784  enddo
785 
786  if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
787 !
788 !---------------------------------------------------------------------------
789 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
790 !---------------------------------------------------------------------------
791 
792  nbpb=0
793  do ig=1,ngrid
794  if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
795 ! stop'On tombe sur le cas particulier de thermcell_dry'
796 ! print*,'On tombe sur le cas particulier de thermcell_plume'
797  nbpb=nbpb+1
798  zw2(ig,l+1)=0.
799  linter(ig)=l+1
800  endif
801 
802  if (zw2(ig,l+1).lt.0.) then
803  linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
804  & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
805  zw2(ig,l+1)=0.
806  endif
807 
808  wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
809 
810  if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
811 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
812 !on rajoute le calcul de lmix_bis
813  if (zqla(ig,l).lt.1.e-10) then
814  lmix_bis(ig)=l+1
815  endif
816  lmix(ig)=l+1
817  wmaxa(ig)=wa_moy(ig,l+1)
818  endif
819  enddo
820 
821  if (nbpb>0) then
822  print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
823  endif
824 
825 !=========================================================================
826 ! FIN DE LA BOUCLE VERTICALE
827  enddo
828 !=========================================================================
829 
830 !on recalcule alim_star_tot
831  do ig=1,ngrid
832  alim_star_tot(ig)=0.
833  enddo
834  do ig=1,ngrid
835  do l=1,lalim(ig)-1
836  alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
837  enddo
838  enddo
839 
840 
841  if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
842 
843 #undef wrgrads_thermcell
844 #ifdef wrgrads_thermcell
845  call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
846  call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
847  call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
848  call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
849  call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
850  call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
851  call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
852 #endif
853 
854 
855  return
856  end
857