LMDZ
thermcell_plume.F90
Go to the documentation of this file.
1 !
2 ! $Id: thermcell_plume.F90 2392 2015-11-11 18:42:11Z fhourdin $
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 ! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
10 !--------------------------------------------------------------------------
11 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
12 !--------------------------------------------------------------------------
13 USE ioipsl, ONLY : getin
14 
15  USE print_control_mod, ONLY: prt_level
16  IMPLICIT NONE
17 
18 #include "YOMCST.h"
19 #include "YOETHF.h"
20 #include "FCTTRE.h"
21 #include "thermcell.h"
22 
23  INTEGER itap
24  INTEGER lunout1,igout
25  INTEGER ngrid,klev
26  REAL ptimestep
27  REAL ztv(ngrid,klev)
28  REAL zthl(ngrid,klev)
29  REAL po(ngrid,klev)
30  REAL zl(ngrid,klev)
31  REAL rhobarz(ngrid,klev)
32  REAL zlev(ngrid,klev+1)
33  REAL pplev(ngrid,klev+1)
34  REAL pphi(ngrid,klev)
35  REAL zpspsk(ngrid,klev)
36  REAL alim_star(ngrid,klev)
37  REAL f0(ngrid)
38  INTEGER lalim(ngrid)
39  integer lev_out ! niveau pour les print
40  integer nbpb
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 ztv_est(ngrid,klev)
65  REAL zqla_est(ngrid,klev)
66  REAL zqsatth(ngrid,klev)
67  REAL zta_est(ngrid,klev)
68  REAL ztemp(ngrid),zqsat(ngrid)
69  REAL zdw2,zdw2bis
70  REAL zw2modif
71  REAL zw2fact,zw2factbis
72  REAL zeps(ngrid,klev)
73 
74  REAL linter(ngrid)
75  INTEGER lmix(ngrid)
76  INTEGER lmix_bis(ngrid)
77  REAL wmaxa(ngrid)
78 
79  INTEGER ig,l,k,lt,it,lm
80 
81  real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
82  real zbuoyjam(ngrid,klev),zdqtjam(ngrid,klev)
83  real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
84  real d_temp(ngrid)
85  real ztv1,ztv2,factinv,zinv,zlmel
86  real zlmelup,zlmeldwn,zlt,zltdwn,zltup
87  real atv1,atv2,btv1,btv2
88  real ztv_est1,ztv_est2
89  real zcor,zdelta,zcvm5,qlbef
90  real zbetalpha, coefzlmel
91  real eps
92  REAL reps,rlvcp,ddt0
93  parameter(ddt0=.01)
94  logical zsat
95  LOGICAL active(ngrid),activetmp(ngrid)
96  REAL fact_gamma,fact_gamma2,fact_epsilon2
97 
98  REAL, SAVE :: fact_epsilon, fact_epsilon_omp=0.002
99  REAL, SAVE :: betalpha, betalpha_omp=0.9
100  REAL, SAVE :: afact, afact_omp=2./3.
101  REAL, SAVE :: fact_shell, fact_shell_omp=1.
102  REAL,SAVE :: detr_min,detr_min_omp=1.e-5
103  REAL,SAVE :: entr_min,entr_min_omp=1.e-5
104  REAL,SAVE :: detr_q_coef,detr_q_coef_omp=0.012
105  REAL,SAVE :: detr_q_power,detr_q_power_omp=0.5
106  REAL,SAVE :: mix0,mix0_omp=0.
107  INTEGER,SAVE :: thermals_flag_alim,thermals_flag_alim_omp=0
108 
109  LOGICAL, SAVE :: first=.true.
110 
111  REAL c2(ngrid,klev)
112 
113  if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
114  zsat=.false.
115 ! Initialisation
116 
117  rlvcp = rlvtt/rcpd
118  IF (first) THEN
119  !$omp master
120 ! FH : if ok_sync=.true. , the time axis is written at each time step
121 ! in the output files. Only at the end in the opposite case
122  CALL getin('thermals_fact_epsilon',fact_epsilon_omp)
123  CALL getin('thermals_betalpha',betalpha_omp)
124  CALL getin('thermals_afact',afact_omp)
125  CALL getin('thermals_fact_shell',fact_shell_omp)
126  CALL getin('thermals_detr_min',detr_min_omp)
127  CALL getin('thermals_entr_min',entr_min_omp)
128  CALL getin('thermals_detr_q_coef',detr_q_coef_omp)
129  CALL getin('thermals_detr_q_power',detr_q_power_omp)
130  CALL getin('thermals_mix0',mix0_omp)
131  CALL getin('thermals_flag_alim',thermals_flag_alim_omp)
132 ! CALL getin('thermals_X',X_omp)
133 ! X=X_omp
134  !$omp end master
135  !$omp barrier
136  fact_epsilon=fact_epsilon_omp
137  betalpha=betalpha_omp
138  afact=afact_omp
139  fact_shell=fact_shell_omp
140  detr_min=detr_min_omp
141  entr_min=entr_min_omp
142  detr_q_coef=detr_q_coef_omp
143  detr_q_power=detr_q_power_omp
144  mix0=mix0_omp
145  thermals_flag_alim=thermals_flag_alim_omp
146  first=.false.
147  ENDIF
148 
149  zbetalpha=betalpha/(1.+betalpha)
150 
151 
152 ! Initialisations des variables r?elles
153 if (1==1) then
154  ztva(:,:)=ztv(:,:)
155  ztva_est(:,:)=ztva(:,:)
156  ztv_est(:,:)=ztv(:,:)
157  ztla(:,:)=zthl(:,:)
158  zqta(:,:)=po(:,:)
159  zqla(:,:)=0.
160  zha(:,:) = ztva(:,:)
161 else
162  ztva(:,:)=0.
163  ztv_est(:,:)=0.
164  ztva_est(:,:)=0.
165  ztla(:,:)=0.
166  zqta(:,:)=0.
167  zha(:,:) =0.
168 endif
169 
170  zqla_est(:,:)=0.
171  zqsatth(:,:)=0.
172  zqla(:,:)=0.
173  detr_star(:,:)=0.
174  entr_star(:,:)=0.
175  alim_star(:,:)=0.
176  alim_star_tot(:)=0.
177  csc(:,:)=0.
178  detr(:,:)=0.
179  entr(:,:)=0.
180  zw2(:,:)=0.
181  zbuoy(:,:)=0.
182  zbuoyjam(:,:)=0.
183  gamma(:,:)=0.
184  zeps(:,:)=0.
185  w_est(:,:)=0.
186  f_star(:,:)=0.
187  wa_moy(:,:)=0.
188  linter(:)=1.
189 ! linter(:)=1.
190 ! Initialisation des variables entieres
191  lmix(:)=1
192  lmix_bis(:)=2
193  wmaxa(:)=0.
194 
195 
196 !-------------------------------------------------------------------------
197 ! On ne considere comme actif que les colonnes dont les deux premieres
198 ! couches sont instables.
199 !-------------------------------------------------------------------------
200 
201  active(:)=ztv(:,1)>ztv(:,2)
202  d_temp(:)=0. ! Pour activer un contraste de temperature a la base
203  ! du panache
204 ! Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
205  CALL thermcell_alim(thermals_flag_alim,ngrid,klev,ztv,d_temp,zlev,alim_star,lalim)
206 
207 !------------------------------------------------------------------------------
208 ! Calcul dans la premiere couche
209 ! On decide dans cette version que le thermique n'est actif que si la premiere
210 ! couche est instable.
211 ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
212 ! dans une couche l>1
213 !------------------------------------------------------------------------------
214 do ig=1,ngrid
215 ! Le panache va prendre au debut les caracteristiques de l'air contenu
216 ! dans cette couche.
217  if (active(ig)) then
218  ztla(ig,1)=zthl(ig,1)
219  zqta(ig,1)=po(ig,1)
220  zqla(ig,1)=zl(ig,1)
221 !cr: attention, prise en compte de f*(1)=1
222  f_star(ig,2)=alim_star(ig,1)
223  zw2(ig,2)=2.*rg*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
224 & *(zlev(ig,2)-zlev(ig,1)) &
225 & *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
226  w_est(ig,2)=zw2(ig,2)
227  endif
228 enddo
229 !
230 
231 !==============================================================================
232 !boucle de calcul de la vitesse verticale dans le thermique
233 !==============================================================================
234 do l=2,klev-1
235 !==============================================================================
236 
237 
238 ! On decide si le thermique est encore actif ou non
239 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
240  do ig=1,ngrid
241  active(ig)=active(ig) &
242 & .and. zw2(ig,l)>1.e-10 &
243 & .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
244  enddo
245 
246 
247 
248 !---------------------------------------------------------------------------
249 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
250 ! sans tenir compte du detrainement et de l'entrainement dans cette
251 ! couche
252 ! C'est a dire qu'on suppose
253 ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
254 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
255 ! avant) a l'alimentation pour avoir un calcul plus propre
256 !---------------------------------------------------------------------------
257 
258  ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
259  call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
260  do ig=1,ngrid
261 ! print*,'active',active(ig),ig,l
262  if(active(ig)) then
263  zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
264  ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+rlvcp*zqla_est(ig,l)
265  zta_est(ig,l)=ztva_est(ig,l)
266  ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
267  ztva_est(ig,l) = ztva_est(ig,l)*(1.+retv*(zqta(ig,l-1) &
268  & -zqla_est(ig,l))-zqla_est(ig,l))
269 
270 
271 !Modif AJAM
272 
273  zbuoy(ig,l)=rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
274  zdz=zlev(ig,l+1)-zlev(ig,l)
275  lmel=fact_thermals_ed_dz*zlev(ig,l)
276 ! lmel=0.09*zlev(ig,l)
277  zlmel=zlev(ig,l)+lmel
278  zlmelup=zlmel+(zdz/2)
279  zlmeldwn=zlmel-(zdz/2)
280 
281  lt=l+1
282  zlt=zlev(ig,lt)
283  zdz3=zlev(ig,lt+1)-zlt
284  zltdwn=zlt-zdz3/2
285  zltup=zlt+zdz3/2
286 
287 !=========================================================================
288 ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
289 !=========================================================================
290 
291 !--------------------------------------------------
292  if (iflag_thermals_ed.lt.8) then
293 !--------------------------------------------------
294 !AJ052014: J'ai remplac?? la boucle do par un do while
295 ! afin de faire moins de calcul dans la boucle
296 !--------------------------------------------------
297  do while (zlmelup.gt.zltup)
298  lt=lt+1
299  zlt=zlev(ig,lt)
300  zdz3=zlev(ig,lt+1)-zlt
301  zltdwn=zlt-zdz3/2
302  zltup=zlt+zdz3/2
303  enddo
304 !--------------------------------------------------
305 !AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
306 ! on cherche o?? se trouve l'altitude d'inversion
307 ! en calculant ztv1 (interpolation de la valeur de
308 ! theta au niveau lt en utilisant les niveaux lt-1 et
309 ! lt-2) et ztv2 (interpolation avec les niveaux lt+1
310 ! et lt+2). Si theta r??ellement calcul??e au niveau lt
311 ! comprise entre ztv1 et ztv2, alors il y a inversion
312 ! et on calcule son altitude zinv en supposant que ztv(lt)
313 ! est une combinaison lineaire de ztv1 et ztv2.
314 ! Ensuite, on calcule la flottabilite en comparant
315 ! la temperature de la couche l a celle de l'air situe
316 ! l+lmel plus haut, ce qui necessite de savoir quel fraction
317 ! de cet air est au-dessus ou en-dessous de l'inversion
318 !--------------------------------------------------
319  atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2))
320  btv1=(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
321  & /(zlev(ig,lt-1)-zlev(ig,lt-2))
322  atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1))
323  btv2=(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
324  & /(zlev(ig,lt+2)-zlev(ig,lt+1))
325 
326  ztv1=atv1*zlt+btv1
327  ztv2=atv2*zlt+btv2
328 
329  if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
330 
331 !--------------------------------------------------
332 !AJ052014: D??calage de zinv qui est entre le haut
333 ! et le bas de la couche lt
334 !--------------------------------------------------
335  factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
336  zinv=zltdwn+zdz3*factinv
337 
338 
339  if (zlmeldwn.ge.zinv) then
340  ztv_est(ig,l)=atv2*zlmel+btv2
341  zbuoyjam(ig,l)=fact_shell*rg*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
342  & +(1.-fact_shell)*zbuoy(ig,l)
343  elseif (zlmelup.ge.zinv) then
344  ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2
345  ztv_est1=atv1*0.5*(zinv+zlmeldwn)+btv1
346  ztv_est(ig,l)=((zlmelup-zinv)/zdz)*ztv_est2+((zinv-zlmeldwn)/zdz)*ztv_est1
347 
348  zbuoyjam(ig,l)=fact_shell*rg*(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- &
349  & ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
350  & ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
351 
352  else
353  ztv_est(ig,l)=atv1*zlmel+btv1
354  zbuoyjam(ig,l)=fact_shell*rg*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
355  & +(1.-fact_shell)*zbuoy(ig,l)
356  endif
357 
358  else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
359 
360  if (zlmeldwn.gt.zltdwn) then
361  zbuoyjam(ig,l)=fact_shell*rg*((ztva_est(ig,l)- &
362  & ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l)
363  else
364  zbuoyjam(ig,l)=fact_shell*rg*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
365  & ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
366  & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
367 
368  endif
369 
370 ! zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
371 ! & ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
372 ! & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
373 ! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
374 ! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
375 ! & po(ig,lt-1))/po(ig,lt-1))
376  endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
377 
378  else ! if (iflag_thermals_ed.lt.8) then
379  lt=l+1
380  zlt=zlev(ig,lt)
381  zdz2=zlev(ig,lt)-zlev(ig,l)
382 
383  do while (lmel.gt.zdz2)
384  lt=lt+1
385  zlt=zlev(ig,lt)
386  zdz2=zlt-zlev(ig,l)
387  enddo
388  zdz3=zlev(ig,lt+1)-zlt
389  zltdwn=zlev(ig,lt)-zdz3/2
390  zlmelup=zlmel+(zdz/2)
391  coefzlmel=min(1.,(zlmelup-zltdwn)/zdz)
392  zbuoyjam(ig,l)=1.*rg*(coefzlmel*(ztva_est(ig,l)- &
393  & ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
394  & ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
395  endif ! if (iflag_thermals_ed.lt.8) then
396 
397 !------------------------------------------------
398 !AJAM:nouveau calcul de w?
399 !------------------------------------------------
400  zdz=zlev(ig,l+1)-zlev(ig,l)
401  zdzbis=zlev(ig,l)-zlev(ig,l-1)
402  zbuoy(ig,l)=rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
403 
404  zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
405  zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
406  zdw2=afact*zbuoy(ig,l)/fact_epsilon
407  zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
408 ! zdw2bis=0.5*(zdw2+zdw2bis)
409  lm=max(1,l-2)
410 ! zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
411 ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
412 ! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
413 ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
414 ! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
415 ! w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
416 ! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
417 ! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
418 ! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
419 
420 !--------------------------------------------------
421 !AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
422 !--------------------------------------------------
423  if (iflag_thermals_ed==8) then
424 ! Ancienne version
425 ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
426 ! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
427 ! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
428 
429  w_est(ig,l+1)=max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
430 
431 ! Nouvelle version Arnaud
432  else
433 ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
434 ! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
435 ! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
436 
437  w_est(ig,l+1)=max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
438 
439 ! w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* &
440 ! & (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
441 ! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
442 
443 
444 
445 ! w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact))
446 
447 ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
448 ! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
449 
450 ! w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
451 
452  endif
453 
454 
455  if (iflag_thermals_ed<6) then
456  zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
457 ! fact_epsilon=0.0005/(zalpha+0.025)**0.5
458 ! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
459  fact_epsilon=0.0002/(zalpha+0.1)
460  zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
461  zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
462  zdw2=afact*zbuoy(ig,l)/fact_epsilon
463  zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
464 ! w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
465 
466 ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
467 ! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
468 ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
469 
470  w_est(ig,l+1)=max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
471 
472 
473  endif
474 !--------------------------------------------------
475 !AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
476 !on fait max(0.0001,.....)
477 !--------------------------------------------------
478 
479 ! if (w_est(ig,l+1).lt.0.) then
480 ! w_est(ig,l+1)=zw2(ig,l)
481 ! w_est(ig,l+1)=0.0001
482 ! endif
483 
484  endif
485  enddo
486 
487 
488 !-------------------------------------------------
489 !calcul des taux d'entrainement et de detrainement
490 !-------------------------------------------------
491 
492  do ig=1,ngrid
493  if (active(ig)) then
494 
495 ! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
496  zw2m=w_est(ig,l+1)
497 ! zw2m=zw2(ig,l)
498  zdz=zlev(ig,l+1)-zlev(ig,l)
499  zbuoy(ig,l)=rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
500 ! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
501  zbuoybis=zbuoy(ig,l)
502  zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
503  zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
504 
505 
506 ! entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., &
507 ! & afact*zbuoybis/zw2m - fact_epsilon )
508 
509 ! entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha* &
510 ! & afact*zbuoybis/zw2m - fact_epsilon )
511 
512 
513 
514 ! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
515 
516 !=========================================================================
517 ! 4. Calcul de l'entrainement et du detrainement
518 !=========================================================================
519 
520 ! entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., &
521 ! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
522 ! entrbis=entr_star(ig,l)
523 
524  if (iflag_thermals_ed.lt.6) then
525  fact_epsilon=0.0002/(zalpha+0.1)
526  endif
527 
528 
529 
530  detr_star(ig,l)=f_star(ig,l)*zdz &
531  & *( mix0 * 0.1 / (zalpha+0.001) &
532  & + max(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m &
533  & + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
534 
535 ! detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
536 ! & ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
537 
538  zbuoy(ig,l)=rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
539 
540  entr_star(ig,l)=f_star(ig,l)*zdz* ( &
541  & mix0 * 0.1 / (zalpha+0.001) &
542  & + zbetalpha*max(entr_min, &
543  & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
544 
545 
546 ! entr_star(ig,l)=f_star(ig,l)*zdz* ( &
547 ! & mix0 * 0.1 / (zalpha+0.001) &
548 ! & + MAX(entr_min, &
549 ! & zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon + &
550 ! & detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))
551 
552 
553 ! entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
554 ! & ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
555 
556 ! entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha* &
557 ! & afact*zbuoy(ig,l)/zw2m &
558 ! & - 1.*fact_epsilon)
559 
560 
561 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
562 ! alim_star et 0 sinon
563  if (l.lt.lalim(ig)) then
564  alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
565  entr_star(ig,l)=0.
566  endif
567 ! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
568 ! alim_star(ig,l)=entrbis
569 ! endif
570 
571 ! print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
572 ! Calcul du flux montant normalise
573  f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
574  & -detr_star(ig,l)
575 
576  endif
577  enddo
578 
579 
580 !============================================================================
581 ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
582 !===========================================================================
583 
584  activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
585  do ig=1,ngrid
586  if (activetmp(ig)) then
587  zsat=.false.
588  ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
589  & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
590  & /(f_star(ig,l+1)+detr_star(ig,l))
591  zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
592  & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
593  & /(f_star(ig,l+1)+detr_star(ig,l))
594 
595  endif
596  enddo
597 
598  ztemp(:)=zpspsk(:,l)*ztla(:,l)
599  call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
600  do ig=1,ngrid
601  if (activetmp(ig)) then
602 ! on ecrit de maniere conservative (sat ou non)
603 ! T = Tl +Lv/Cp ql
604  zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
605  ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+rlvcp*zqla(ig,l)
606  ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
607 !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
608  zha(ig,l) = ztva(ig,l)
609  ztva(ig,l) = ztva(ig,l)*(1.+retv*(zqta(ig,l) &
610  & -zqla(ig,l))-zqla(ig,l))
611  zbuoy(ig,l)=rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
612  zdz=zlev(ig,l+1)-zlev(ig,l)
613  zdzbis=zlev(ig,l)-zlev(ig,l-1)
614  zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
615 !!!!!!! fact_epsilon=0.002
616  zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
617  zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
618  zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
619  zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
620 ! zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
621 ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
622 ! lm=Max(1,l-2)
623 ! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
624 ! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
625 ! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
626 ! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
627 ! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
628 ! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
629 ! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
630 ! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
631 ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
632  if (iflag_thermals_ed==8) then
633  zw2(ig,l+1)=max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
634  else
635  zw2(ig,l+1)=max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
636  endif
637 ! zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* &
638 ! & (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
639 ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
640 
641 
642  if (iflag_thermals_ed.lt.6) then
643  zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
644 ! fact_epsilon=0.0005/(zalpha+0.025)**0.5
645 ! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
646  fact_epsilon=0.0002/(zalpha+0.1)**1
647  zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
648  zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
649  zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
650  zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
651 
652 ! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
653 ! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
654 ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
655 ! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
656  zw2(ig,l+1)=max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
657 
658  endif
659 
660 
661  endif
662  enddo
663 
664  if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
665 !
666 !===========================================================================
667 ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
668 !===========================================================================
669 
670  nbpb=0
671  do ig=1,ngrid
672  if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
673 ! stop'On tombe sur le cas particulier de thermcell_dry'
674 ! print*,'On tombe sur le cas particulier de thermcell_plume'
675  nbpb=nbpb+1
676  zw2(ig,l+1)=0.
677  linter(ig)=l+1
678  endif
679 
680  if (zw2(ig,l+1).lt.0.) then
681  linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
682  & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
683  zw2(ig,l+1)=0.
684 !+CR:04/05/12:correction calcul linter pour calcul de zmax continu
685  elseif (f_star(ig,l+1).lt.0.) then
686  linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
687  & -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
688  zw2(ig,l+1)=0.
689 !fin CR:04/05/12
690  endif
691 
692  wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
693 
694  if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
695 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
696 !on rajoute le calcul de lmix_bis
697  if (zqla(ig,l).lt.1.e-10) then
698  lmix_bis(ig)=l+1
699  endif
700  lmix(ig)=l+1
701  wmaxa(ig)=wa_moy(ig,l+1)
702  endif
703  enddo
704 
705  if (nbpb>0) then
706  print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
707  endif
708 
709 !=========================================================================
710 ! FIN DE LA BOUCLE VERTICALE
711  enddo
712 !=========================================================================
713 
714 !on recalcule alim_star_tot
715  do ig=1,ngrid
716  alim_star_tot(ig)=0.
717  enddo
718  do ig=1,ngrid
719  do l=1,lalim(ig)-1
720  alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
721  enddo
722  enddo
723 
724 
725  if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
726 
727 #undef wrgrads_thermcell
728 #ifdef wrgrads_thermcell
729  call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
730  call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
731  call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
732  call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
733  call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
734  call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
735  call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
736 #endif
737 
738 
739  return
740  end
741 
742 
743 
744 
745 
746 
747 
748 
749 
750 
751 
752 
753 
754 
755 
756 
757 
758 
759 
760 
761 
762 
763 
764 
765 
766 
767 
768 
769 
770 
771 
772 
773 
774 
775 
776 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
777 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
778 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
779 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
780 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
781 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
782  SUBROUTINE thermcellv1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
783 & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
784 & lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
785 & ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
786 & ,lev_out,lunout1,igout)
787 !& ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
788 
789 !--------------------------------------------------------------------------
790 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
791 ! Version conforme a l'article de Rio et al. 2010.
792 ! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
793 !--------------------------------------------------------------------------
794 
795  USE print_control_mod, ONLY: prt_level
796  IMPLICIT NONE
797 
798 #include "YOMCST.h"
799 #include "YOETHF.h"
800 #include "FCTTRE.h"
801 #include "thermcell.h"
802 
803  INTEGER itap
804  INTEGER lunout1,igout
805  INTEGER ngrid,klev
806  REAL ptimestep
807  REAL ztv(ngrid,klev)
808  REAL zthl(ngrid,klev)
809  REAL po(ngrid,klev)
810  REAL zl(ngrid,klev)
811  REAL rhobarz(ngrid,klev)
812  REAL zlev(ngrid,klev+1)
813  REAL pplev(ngrid,klev+1)
814  REAL pphi(ngrid,klev)
815  REAL zpspsk(ngrid,klev)
816  REAL alim_star(ngrid,klev)
817  REAL f0(ngrid)
818  INTEGER lalim(ngrid)
819  integer lev_out ! niveau pour les print
820  integer nbpb
821 
822  real alim_star_tot(ngrid)
823 
824  REAL ztva(ngrid,klev)
825  REAL ztla(ngrid,klev)
826  REAL zqla(ngrid,klev)
827  REAL zqta(ngrid,klev)
828  REAL zha(ngrid,klev)
829 
830  REAL detr_star(ngrid,klev)
831  REAL coefc
832  REAL entr_star(ngrid,klev)
833  REAL detr(ngrid,klev)
834  REAL entr(ngrid,klev)
835 
836  REAL csc(ngrid,klev)
837 
838  REAL zw2(ngrid,klev+1)
839  REAL w_est(ngrid,klev+1)
840  REAL f_star(ngrid,klev+1)
841  REAL wa_moy(ngrid,klev+1)
842 
843  REAL ztva_est(ngrid,klev)
844  REAL zqla_est(ngrid,klev)
845  REAL zqsatth(ngrid,klev)
846  REAL zta_est(ngrid,klev)
847  REAL zbuoyjam(ngrid,klev)
848  REAL ztemp(ngrid),zqsat(ngrid)
849  REAL zdw2
850  REAL zw2modif
851  REAL zw2fact
852  REAL zeps(ngrid,klev)
853 
854  REAL linter(ngrid)
855  INTEGER lmix(ngrid)
856  INTEGER lmix_bis(ngrid)
857  REAL wmaxa(ngrid)
858 
859  INTEGER ig,l,k
860 
861  real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
862  real zbuoybis
863  real zcor,zdelta,zcvm5,qlbef,zdz2
864  real betalpha,zbetalpha
865  real eps, afact
866  REAL reps,rlvcp,ddt0
867  parameter(ddt0=.01)
868  logical zsat
869  LOGICAL active(ngrid),activetmp(ngrid)
870  REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
871  REAL c2(ngrid,klev)
872  zsat=.false.
873 ! Initialisation
874 
875  rlvcp = rlvtt/rcpd
876  fact_epsilon=0.002
877  betalpha=0.9
878  afact=2./3.
879 
880  zbetalpha=betalpha/(1.+betalpha)
881 
882 
883 ! Initialisations des variables reeles
884 if (1==1) then
885  ztva(:,:)=ztv(:,:)
886  ztva_est(:,:)=ztva(:,:)
887  ztla(:,:)=zthl(:,:)
888  zqta(:,:)=po(:,:)
889  zha(:,:) = ztva(:,:)
890 else
891  ztva(:,:)=0.
892  ztva_est(:,:)=0.
893  ztla(:,:)=0.
894  zqta(:,:)=0.
895  zha(:,:) =0.
896 endif
897 
898  zqla_est(:,:)=0.
899  zqsatth(:,:)=0.
900  zqla(:,:)=0.
901  detr_star(:,:)=0.
902  entr_star(:,:)=0.
903  alim_star(:,:)=0.
904  alim_star_tot(:)=0.
905  csc(:,:)=0.
906  detr(:,:)=0.
907  entr(:,:)=0.
908  zw2(:,:)=0.
909  zbuoy(:,:)=0.
910  zbuoyjam(:,:)=0.
911  gamma(:,:)=0.
912  zeps(:,:)=0.
913  w_est(:,:)=0.
914  f_star(:,:)=0.
915  wa_moy(:,:)=0.
916  linter(:)=1.
917 ! linter(:)=1.
918 ! Initialisation des variables entieres
919  lmix(:)=1
920  lmix_bis(:)=2
921  wmaxa(:)=0.
922  lalim(:)=1
923 
924 
925 !-------------------------------------------------------------------------
926 ! On ne considere comme actif que les colonnes dont les deux premieres
927 ! couches sont instables.
928 !-------------------------------------------------------------------------
929  active(:)=ztv(:,1)>ztv(:,2)
930 
931 !-------------------------------------------------------------------------
932 ! Definition de l'alimentation a l'origine dans thermcell_init
933 !-------------------------------------------------------------------------
934  do l=1,klev-1
935  do ig=1,ngrid
936  if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
937  alim_star(ig,l)=max((ztv(ig,l)-ztv(ig,l+1)),0.) &
938  & *sqrt(zlev(ig,l+1))
939  lalim(ig)=l+1
940  alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
941  endif
942  enddo
943  enddo
944  do l=1,klev
945  do ig=1,ngrid
946  if (alim_star_tot(ig) > 1.e-10 ) then
947  alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
948  endif
949  enddo
950  enddo
951  alim_star_tot(:)=1.
952 
953 
954 
955 !------------------------------------------------------------------------------
956 ! Calcul dans la premiere couche
957 ! On decide dans cette version que le thermique n'est actif que si la premiere
958 ! couche est instable.
959 ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
960 ! dans une couche l>1
961 !------------------------------------------------------------------------------
962 do ig=1,ngrid
963 ! Le panache va prendre au debut les caracteristiques de l'air contenu
964 ! dans cette couche.
965  if (active(ig)) then
966  ztla(ig,1)=zthl(ig,1)
967  zqta(ig,1)=po(ig,1)
968  zqla(ig,1)=zl(ig,1)
969 !cr: attention, prise en compte de f*(1)=1
970  f_star(ig,2)=alim_star(ig,1)
971  zw2(ig,2)=2.*rg*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
972 & *(zlev(ig,2)-zlev(ig,1)) &
973 & *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
974  w_est(ig,2)=zw2(ig,2)
975  endif
976 enddo
977 !
978 
979 !==============================================================================
980 !boucle de calcul de la vitesse verticale dans le thermique
981 !==============================================================================
982 do l=2,klev-1
983 !==============================================================================
984 
985 
986 ! On decide si le thermique est encore actif ou non
987 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
988  do ig=1,ngrid
989  active(ig)=active(ig) &
990 & .and. zw2(ig,l)>1.e-10 &
991 & .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
992  enddo
993 
994 
995 
996 !---------------------------------------------------------------------------
997 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
998 ! sans tenir compte du detrainement et de l'entrainement dans cette
999 ! couche
1000 ! C'est a dire qu'on suppose
1001 ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
1002 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
1003 ! avant) a l'alimentation pour avoir un calcul plus propre
1004 !---------------------------------------------------------------------------
1005 
1006  ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
1007  call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
1008 
1009  do ig=1,ngrid
1010 ! print*,'active',active(ig),ig,l
1011  if(active(ig)) then
1012  zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
1013  ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+rlvcp*zqla_est(ig,l)
1014  zta_est(ig,l)=ztva_est(ig,l)
1015  ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
1016  ztva_est(ig,l) = ztva_est(ig,l)*(1.+retv*(zqta(ig,l-1) &
1017  & -zqla_est(ig,l))-zqla_est(ig,l))
1018 
1019 !------------------------------------------------
1020 !AJAM:nouveau calcul de w?
1021 !------------------------------------------------
1022  zdz=zlev(ig,l+1)-zlev(ig,l)
1023  zbuoy(ig,l)=rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
1024 
1025  zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
1026  zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
1027  w_est(ig,l+1)=max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
1028 
1029 
1030  if (w_est(ig,l+1).lt.0.) then
1031  w_est(ig,l+1)=zw2(ig,l)
1032  endif
1033  endif
1034  enddo
1035 
1036 
1037 !-------------------------------------------------
1038 !calcul des taux d'entrainement et de detrainement
1039 !-------------------------------------------------
1040 
1041  do ig=1,ngrid
1042  if (active(ig)) then
1043 
1044  zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
1045  zw2m=w_est(ig,l+1)
1046  zdz=zlev(ig,l+1)-zlev(ig,l)
1047  zbuoy(ig,l)=rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
1048 ! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
1049  zbuoybis=zbuoy(ig,l)
1050  zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
1051  zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
1052 
1053 
1054  entr_star(ig,l)=f_star(ig,l)*zdz* zbetalpha*max(0., &
1055  & afact*zbuoybis/zw2m - fact_epsilon )
1056 
1057 
1058  detr_star(ig,l)=f_star(ig,l)*zdz &
1059  & *max(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m &
1060  & + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
1061 
1062 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
1063 ! alim_star et 0 sinon
1064  if (l.lt.lalim(ig)) then
1065  alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
1066  entr_star(ig,l)=0.
1067  endif
1068 
1069 ! Calcul du flux montant normalise
1070  f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
1071  & -detr_star(ig,l)
1072 
1073  endif
1074  enddo
1075 
1076 
1077 !----------------------------------------------------------------------------
1078 !calcul de la vitesse verticale en melangeant Tl et qt du thermique
1079 !---------------------------------------------------------------------------
1080  activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
1081  do ig=1,ngrid
1082  if (activetmp(ig)) then
1083  zsat=.false.
1084  ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
1085  & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
1086  & /(f_star(ig,l+1)+detr_star(ig,l))
1087  zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
1088  & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
1089  & /(f_star(ig,l+1)+detr_star(ig,l))
1090 
1091  endif
1092  enddo
1093 
1094  ztemp(:)=zpspsk(:,l)*ztla(:,l)
1095  call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
1096 
1097  do ig=1,ngrid
1098  if (activetmp(ig)) then
1099 ! on ecrit de maniere conservative (sat ou non)
1100 ! T = Tl +Lv/Cp ql
1101  zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
1102  ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+rlvcp*zqla(ig,l)
1103  ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1104 !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
1105  zha(ig,l) = ztva(ig,l)
1106  ztva(ig,l) = ztva(ig,l)*(1.+retv*(zqta(ig,l) &
1107  & -zqla(ig,l))-zqla(ig,l))
1108  zbuoy(ig,l)=rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
1109  zdz=zlev(ig,l+1)-zlev(ig,l)
1110  zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
1111 
1112  zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
1113  zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
1114  zw2(ig,l+1)=max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
1115  endif
1116  enddo
1117 
1118  if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
1119 !
1120 !---------------------------------------------------------------------------
1121 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
1122 !---------------------------------------------------------------------------
1123 
1124  nbpb=0
1125  do ig=1,ngrid
1126  if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
1127 ! stop'On tombe sur le cas particulier de thermcell_dry'
1128 ! print*,'On tombe sur le cas particulier de thermcell_plume'
1129  nbpb=nbpb+1
1130  zw2(ig,l+1)=0.
1131  linter(ig)=l+1
1132  endif
1133 
1134  if (zw2(ig,l+1).lt.0.) then
1135  linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
1136  & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1137  zw2(ig,l+1)=0.
1138  elseif (f_star(ig,l+1).lt.0.) then
1139  linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
1140  & -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
1141 ! print*,"linter plume", linter(ig)
1142  zw2(ig,l+1)=0.
1143  endif
1144 
1145  wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1146 
1147  if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1148 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
1149 !on rajoute le calcul de lmix_bis
1150  if (zqla(ig,l).lt.1.e-10) then
1151  lmix_bis(ig)=l+1
1152  endif
1153  lmix(ig)=l+1
1154  wmaxa(ig)=wa_moy(ig,l+1)
1155  endif
1156  enddo
1157 
1158  if (nbpb>0) then
1159  print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
1160  endif
1161 
1162 !=========================================================================
1163 ! FIN DE LA BOUCLE VERTICALE
1164  enddo
1165 !=========================================================================
1166 
1167 !on recalcule alim_star_tot
1168  do ig=1,ngrid
1169  alim_star_tot(ig)=0.
1170  enddo
1171  do ig=1,ngrid
1172  do l=1,lalim(ig)-1
1173  alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
1174  enddo
1175  enddo
1176 
1177 
1178  if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
1179 
1180 #undef wrgrads_thermcell
1181 #ifdef wrgrads_thermcell
1182  call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
1183  call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
1184  call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
1185  call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
1186  call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
1187  call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
1188  call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
1189 #endif
1190 
1191 
1192  return
1193  end
itap
!$Id klon initialisation mois suivants day_rain itap
Definition: calcul_divers.h:18
l
c c $Id c c calculs statistiques distribution nuage ftion du regime dynamique c c Ce calcul doit etre fait a partir de valeurs mensuelles CALL nbregdyn DO kmaxm1 DO l
Definition: calcul_REGDYN.h:13
dimphy::klev
integer, save klev
Definition: dimphy.F90:7
eps
!$Id!Thermodynamical constants for t0 real clmci real eps
Definition: cvthermo.h:6
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
prt_level
!FH On elimine toutes les clefs physiques dans la dynamique prt_level
Definition: write_paramLMDZ_dyn.h:9
parameter
!$Header!integer nvarmx parameter(nfmx=10, imx=200, jmx=150, lmx=200, nvarmx=1000) real xd(imx
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
print_control_mod
Definition: print_control_mod.F90:2
iflag_thermals_ed
nsplit_thermals!nrlmd le iflag_clos_bl tau_trig_deep real::s_trig!fin nrlmd le fact_thermals_ed_dz iflag_wake iflag_thermals_closure common ctherm1 iflag_thermals_closure common ctherm2 fact_thermals_ed_dz common ctherm4 iflag_wake common ctherm5 iflag_thermals_ed
Definition: thermcell.h:12
thermcell_qsat
subroutine thermcell_qsat(klon, active, pplev, ztemp, zqta, zqsat)
Definition: thermcell_qsat.F90:2
thermcell_alim
subroutine thermcell_alim(flag, ngrid, klev, ztv, d_temp, zlev, alim_star, lalim)
Definition: thermcell_alim.F90:5
rg
real rg
Definition: comcstphy.h:1