My Project
 All Classes Files Functions Variables Macros
thermcellV0_main.F90
Go to the documentation of this file.
1 !
2 ! $Id$
3 !
4  SUBROUTINE thermcellv0_main(itap,ngrid,nlay,ptimestep &
5  & ,pplay,pplev,pphi,debut &
6  & ,pu,pv,pt,po &
7  & ,pduadj,pdvadj,pdtadj,pdoadj &
8  & ,fm0,entr0,detr0,zqta,zqla,lmax &
9  & ,ratqscth,ratqsdiff,zqsatth &
10  & ,r_aspect,l_mix,tau_thermals &
11  & ,ale_bl,alp_bl,lalim_conv,wght_th &
12  & ,zmax0, f0,zw2,fraca)
13 
14  USE dimphy
15  USE comgeomphy , ONLY:rlond,rlatd
16  IMPLICIT NONE
17 
18 !=======================================================================
19 ! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
20 ! Version du 09.02.07
21 ! Calcul du transport vertical dans la couche limite en presence
22 ! de "thermiques" explicitement representes avec processus nuageux
23 !
24 ! Réécriture à partir d'un listing papier à Habas, le 14/02/00
25 !
26 ! le thermique est supposé homogène et dissipé par mélange avec
27 ! son environnement. la longueur l_mix contrôle l'efficacité du
28 ! mélange
29 !
30 ! Le calcul du transport des différentes espèces se fait en prenant
31 ! en compte:
32 ! 1. un flux de masse montant
33 ! 2. un flux de masse descendant
34 ! 3. un entrainement
35 ! 4. un detrainement
36 !
37 !=======================================================================
38 
39 !-----------------------------------------------------------------------
40 ! declarations:
41 ! -------------
42 
43 #include "dimensions.h"
44 #include "YOMCST.h"
45 #include "YOETHF.h"
46 #include "FCTTRE.h"
47 #include "iniprint.h"
48 
49 ! arguments:
50 ! ----------
51 
52 !IM 140508
53  INTEGER itap
54 
55  INTEGER ngrid,nlay,w2di
56  real tau_thermals
57  real ptimestep,l_mix,r_aspect
58  REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
59  REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
60  REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
61  REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
62  REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
63  real pphi(ngrid,nlay)
64 
65 ! local:
66 ! ------
67 
68  integer icount
69  data icount/0/
70  save icount
71 !$OMP THREADPRIVATE(icount)
72 
73  integer,save :: igout=1
74 !$OMP THREADPRIVATE(igout)
75  integer,save :: lunout1=6
76 !$OMP THREADPRIVATE(lunout1)
77  integer,save :: lev_out=10
78 !$OMP THREADPRIVATE(lev_out)
79 
80  INTEGER ig,k,l,ll
81  real zsortie1d(klon)
82  INTEGER lmax(klon),lmin(klon),lalim(klon)
83  INTEGER lmix(klon)
84  INTEGER lmix_bis(klon)
85  real linter(klon)
86  real zmix(klon)
87  real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1)
88 ! real fraca(klon,klev)
89 
90  real zmax_sec(klon)
91 !on garde le zmax du pas de temps precedent
92  real zmax0(klon)
93 !FH/IM save zmax0
94 
95  real lambda
96 
97  real zlev(klon,klev+1),zlay(klon,klev)
98  real deltaz(klon,klev)
99  REAL zh(klon,klev)
100  real zthl(klon,klev),zdthladj(klon,klev)
101  REAL ztv(klon,klev)
102  real zu(klon,klev),zv(klon,klev),zo(klon,klev)
103  real zl(klon,klev)
104  real zsortie(klon,klev)
105  real zva(klon,klev)
106  real zua(klon,klev)
107  real zoa(klon,klev)
108 
109  real zta(klon,klev)
110  real zha(klon,klev)
111  real fraca(klon,klev+1)
112  real zf,zf2
113  real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
114  real q2(klon,klev)
115 ! FH probleme de dimensionnement avec l'allocation dynamique
116 ! common/comtherm/thetath2,wth2
117 
118  real ratqscth(klon,klev)
119  real var
120  real vardiff
121  real ratqsdiff(klon,klev)
122 
123  logical sorties
124  real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
125  real zpspsk(klon,klev)
126 
127  real wmax(klon)
128  real wmax_sec(klon)
129  real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
130  real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
131 
132  real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
133 !niveau de condensation
134  integer nivcon(klon)
135  real zcon(klon)
136  REAL chi
137  real zcon2(klon)
138  real pcon(klon)
139  real zqsat(klon,klev)
140  real zqsatth(klon,klev)
141 
142  real f_star(klon,klev+1),entr_star(klon,klev)
143  real detr_star(klon,klev)
144  real alim_star_tot(klon),alim_star2(klon)
145  real alim_star(klon,klev)
146  real f(klon), f0(klon)
147 !FH/IM save f0
148  real zlevinter(klon)
149  logical debut
150  real seuil
151 
152 ! Declaration uniquement pour les sorties dans thermcell_out3d.
153 ! Inutilise en 3D
154  real wthl(klon,klev)
155  real wthv(klon,klev)
156  real wq(klon,klev)
157 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
158 
159 
160 !
161  !nouvelles variables pour la convection
162  real ale_bl(klon)
163  real alp_bl(klon)
164  real alp_int(klon)
165  real ale_int(klon)
166  integer n_int(klon)
167  real fm_tot(klon)
168  real wght_th(klon,klev)
169  integer lalim_conv(klon)
170 !v1d logical therm
171 !v1d save therm
172 
173  character*2 str2
174  character*10 str10
175 
176  character (len=20) :: modname='thermcellV0_main'
177  character (len=80) :: abort_message
178 
179  EXTERNAL scopy
180 !
181 
182 !-----------------------------------------------------------------------
183 ! initialisation:
184 ! ---------------
185 !
186 
187  seuil=0.25
188 
189  if (debut) then
190  fm0=0.
191  entr0=0.
192  detr0=0.
193 
194 
195 #undef wrgrads_thermcell
196 #ifdef wrgrads_thermcell
197 ! Initialisation des sorties grads pour les thermiques.
198 ! Pour l'instant en 1D sur le point igout.
199 ! Utilise par thermcell_out3d.h
200  str10='therm'
201  call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
202  & rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1., &
203  & ptimestep,str10,'therm ')
204 #endif
205 
206 
207 
208  endif
209 
210  fm=0. ; entr=0. ; detr=0.
211 
212  icount=icount+1
213 
214 !IM 090508 beg
215 !print*,'====================================================================='
216 !print*,'====================================================================='
217 !print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
218 !print*,'====================================================================='
219 !print*,'====================================================================='
220 !IM 090508 end
221 
222  if (prt_level.ge.1) print*,'thermcell_main V4'
223 
224  sorties=.true.
225  IF(ngrid.NE.klon) THEN
226  print*
227  print*,'STOP dans convadj'
228  print*,'ngrid =',ngrid
229  print*,'klon =',klon
230  ENDIF
231 !
232 !Initialisation
233 !
234  if (prt_level.ge.10)write(lunout,*) &
235  & 'WARNING thermcell_main f0=max(f0,1.e-2)'
236  do ig=1,klon
237  f0(ig)=max(f0(ig),1.e-2)
238  enddo
239 
240 !-----------------------------------------------------------------------
241 ! Calcul de T,q,ql a partir de Tl et qT dans l environnement
242 ! --------------------------------------------------------------------
243 !
244  CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
245  & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
246 
247  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
248 
249 !------------------------------------------------------------------------
250 ! --------------------
251 !
252 !
253 ! + + + + + + + + + + +
254 !
255 !
256 ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
257 ! wh,wt,wo ...
258 !
259 ! + + + + + + + + + + + zh,zu,zv,zo,rho
260 !
261 !
262 ! -------------------- zlev(1)
263 ! \\\\\\\\\\\\\\\\\\\\
264 !
265 !
266 
267 !-----------------------------------------------------------------------
268 ! Calcul des altitudes des couches
269 !-----------------------------------------------------------------------
270 
271  do l=2,nlay
272  zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/rg
273  enddo
274  zlev(:,1)=0.
275  zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/rg
276  do l=1,nlay
277  zlay(:,l)=pphi(:,l)/rg
278  enddo
279 !calcul de l epaisseur des couches
280  do l=1,nlay
281  deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
282  enddo
283 
284 ! print*,'2 OK convect8'
285 !-----------------------------------------------------------------------
286 ! Calcul des densites
287 !-----------------------------------------------------------------------
288 
289  do l=1,nlay
290  rho(:,l)=pplay(:,l)/(zpspsk(:,l)*rd*ztv(:,l))
291  enddo
292 
293 !IM
294  if (prt_level.ge.10)write(lunout,*) &
295  & 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
296  rhobarz(:,1)=rho(:,1)
297 
298  do l=2,nlay
299  rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
300  enddo
301 
302 !calcul de la masse
303  do l=1,nlay
304  masse(:,l)=(pplev(:,l)-pplev(:,l+1))/rg
305  enddo
306 
307  if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
308 
309 !------------------------------------------------------------------
310 !
311 ! /|\
312 ! -------- | F_k+1 -------
313 ! ----> D_k
314 ! /|\ <---- E_k , A_k
315 ! -------- | F_k ---------
316 ! ----> D_k-1
317 ! <---- E_k-1 , A_k-1
318 !
319 !
320 !
321 !
322 !
323 ! ---------------------------
324 !
325 ! ----- F_lmax+1=0 ---------- \
326 ! lmax (zmax) |
327 ! --------------------------- |
328 ! |
329 ! --------------------------- |
330 ! |
331 ! --------------------------- |
332 ! |
333 ! --------------------------- |
334 ! |
335 ! --------------------------- |
336 ! | E
337 ! --------------------------- | D
338 ! |
339 ! --------------------------- |
340 ! |
341 ! --------------------------- \ |
342 ! lalim | |
343 ! --------------------------- | |
344 ! | |
345 ! --------------------------- | |
346 ! | A |
347 ! --------------------------- | |
348 ! | |
349 ! --------------------------- | |
350 ! lmin (=1 pour le moment) | |
351 ! ----- F_lmin=0 ------------ / /
352 !
353 ! ---------------------------
354 ! //////////////////////////
355 !
356 !
357 !=============================================================================
358 ! Calculs initiaux ne faisant pas intervenir les changements de phase
359 !=============================================================================
360 
361 !------------------------------------------------------------------
362 ! 1. alim_star est le profil vertical de l'alimentation à la base du
363 ! panache thermique, calculé à partir de la flotabilité de l'air sec
364 ! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
365 !------------------------------------------------------------------
366 !
367  entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
368  CALL thermcellv0_init(ngrid,nlay,ztv,zlay,zlev, &
369  & lalim,lmin,alim_star,alim_star_tot,lev_out)
370 
371 call testV0_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin ')
372 call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
373 
374 
375  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
376  if (prt_level.ge.10) then
377  write(lunout1,*) 'Dans thermcell_main 1'
378  write(lunout1,*) 'lmin ',lmin(igout)
379  write(lunout1,*) 'lalim ',lalim(igout)
380  write(lunout1,*) ' ig l alim_star thetav'
381  write(lunout1,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
382  & ,ztv(igout,l),l=1,lalim(igout)+4)
383  endif
384 
385 !v1d do ig=1,klon
386 !v1d if (alim_star(ig,1).gt.1.e-10) then
387 !v1d therm=.true.
388 !v1d endif
389 !v1d enddo
390 !-----------------------------------------------------------------------------
391 ! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
392 ! panache sec conservatif (e=d=0) alimente selon alim_star
393 ! Il s'agit d'un calcul de type CAPE
394 ! zmax_sec est utilisé pour déterminer la géométrie du thermique.
395 !------------------------------------------------------------------------------
396 !
397  CALL thermcellv0_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
398  & lalim,lmin,zmax_sec,wmax_sec,lev_out)
399 
400 call testV0_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
401 call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
402 
403  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
404  if (prt_level.ge.10) then
405  write(lunout1,*) 'Dans thermcell_main 1b'
406  write(lunout1,*) 'lmin ',lmin(igout)
407  write(lunout1,*) 'lalim ',lalim(igout)
408  write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
409  write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
410  & ,l=1,lalim(igout)+4)
411  endif
412 
413 
414 
415 !---------------------------------------------------------------------------------
416 !calcul du melange et des variables dans le thermique
417 !--------------------------------------------------------------------------------
418 !
419  if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
420 !IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
421  CALL thermcellv0_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
422  & zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot, &
423  & lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, &
424  & ztla,zqla,zqta,zha,zw2,zw_est,zqsatth,lmix,lmix_bis,linter &
425  & ,lev_out,lunout1,igout)
426  if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
427 
428  call testv0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
429  call testv0_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
430 
431  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
432  if (prt_level.ge.10) then
433  write(lunout1,*) 'Dans thermcell_main 2'
434  write(lunout1,*) 'lmin ',lmin(igout)
435  write(lunout1,*) 'lalim ',lalim(igout)
436  write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
437  write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
438  & ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
439  endif
440 
441 !-------------------------------------------------------------------------------
442 ! Calcul des caracteristiques du thermique:zmax,zmix,wmax
443 !-------------------------------------------------------------------------------
444 !
445  CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
446  & zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
447 
448 
449  call testv0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
450  call testv0_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
451  call testv0_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
452  call testv0_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
453 
454  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
455 
456 !-------------------------------------------------------------------------------
457 ! Fermeture,determination de f
458 !-------------------------------------------------------------------------------
459 !
460 !avant closure: on redéfinit lalim, alim_star_tot et alim_star
461 ! do ig=1,klon
462 ! do l=2,lalim(ig)
463 ! alim_star(ig,l)=entr_star(ig,l)
464 ! entr_star(ig,l)=0.
465 ! enddo
466 ! enddo
467 
468  CALL thermcellv0_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
469  & zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
470 
471  if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
472 
473  if (tau_thermals>1.) then
474  lambda=exp(-ptimestep/tau_thermals)
475  f0=(1.-lambda)*f+lambda*f0
476  else
477  f0=f
478  endif
479 
480 ! Test valable seulement en 1D mais pas genant
481  if (.not. (f0(1).ge.0.) ) then
482  abort_message = .lt.'Dans thermcell_main f0(1)0 '
483  CALL abort_gcm(modname,abort_message,1)
484  endif
485 
486 !-------------------------------------------------------------------------------
487 !deduction des flux
488 !-------------------------------------------------------------------------------
489 
490  CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
491  & lalim,lmax,alim_star, &
492  & entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
493  & detr,zqla,lev_out,lunout1,igout)
494 !IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
495 
496  if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
497  call testv0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
498  call testv0_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
499 
500 !------------------------------------------------------------------
501 ! On ne prend pas directement les profils issus des calculs precedents
502 ! mais on s'autorise genereusement une relaxation vers ceci avec
503 ! une constante de temps tau_thermals (typiquement 1800s).
504 !------------------------------------------------------------------
505 
506  if (tau_thermals>1.) then
507  lambda=exp(-ptimestep/tau_thermals)
508  fm0=(1.-lambda)*fm+lambda*fm0
509  entr0=(1.-lambda)*entr+lambda*entr0
510 ! detr0=(1.-lambda)*detr+lambda*detr0
511  else
512  fm0=fm
513  entr0=entr
514  detr0=detr
515  endif
516 
517 !c------------------------------------------------------------------
518 ! calcul du transport vertical
519 !------------------------------------------------------------------
520 
521  call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse, &
522  & zthl,zdthladj,zta,lev_out)
523  call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse, &
524  & po,pdoadj,zoa,lev_out)
525 
526 !------------------------------------------------------------------
527 ! Calcul de la fraction de l'ascendance
528 !------------------------------------------------------------------
529  do ig=1,klon
530  fraca(ig,1)=0.
531  fraca(ig,nlay+1)=0.
532  enddo
533  do l=2,nlay
534  do ig=1,klon
535  if (zw2(ig,l).gt.1.e-10) then
536  fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
537  else
538  fraca(ig,l)=0.
539  endif
540  enddo
541  enddo
542 
543 !------------------------------------------------------------------
544 ! calcul du transport vertical du moment horizontal
545 !------------------------------------------------------------------
546 
547 !IM 090508
548  if (1.eq.1) then
549 !IM 070508 vers. _dq
550 ! if (1.eq.0) then
551 
552 
553 ! Calcul du transport de V tenant compte d'echange par gradient
554 ! de pression horizontal avec l'environnement
555 
556  call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
557  & ,fraca,zmax &
558  & ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
559 !IM 050508 & ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
560  else
561 
562 ! calcul purement conservatif pour le transport de V
563  call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse &
564  & ,zu,pduadj,zua,lev_out)
565  call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse &
566  & ,zv,pdvadj,zva,lev_out)
567  endif
568 
569 ! print*,'13 OK convect8'
570  do l=1,nlay
571  do ig=1,ngrid
572  pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
573  enddo
574  enddo
575 
576  if (prt_level.ge.1) print*,'14 OK convect8'
577 !------------------------------------------------------------------
578 ! Calculs de diagnostiques pour les sorties
579 !------------------------------------------------------------------
580 !calcul de fraca pour les sorties
581 
582  if (sorties) then
583  if (prt_level.ge.1) print*,'14a OK convect8'
584 ! calcul du niveau de condensation
585 ! initialisation
586  do ig=1,ngrid
587  nivcon(ig)=0
588  zcon(ig)=0.
589  enddo
590 !nouveau calcul
591  do ig=1,ngrid
592  chi=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
593  pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**chi
594  enddo
595 !IM do k=1,nlay
596  do k=1,nlay-1
597  do ig=1,ngrid
598  if ((pcon(ig).le.pplay(ig,k)) &
599  & .and.(pcon(ig).gt.pplay(ig,k+1))) then
600  zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(rg*rho(ig,k))/100.
601  endif
602  enddo
603  enddo
604 !IM
605  do ig=1,ngrid
606  if (pcon(ig).le.pplay(ig,nlay)) then
607  zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(rg*rho(ig,nlay))/100.
608  abort_message = 'thermcellV0_main: les thermiques vont trop haut '
609  CALL abort_gcm(modname,abort_message,1)
610  endif
611  enddo
612  if (prt_level.ge.1) print*,'14b OK convect8'
613  do k=nlay,1,-1
614  do ig=1,ngrid
615  if (zqla(ig,k).gt.1e-10) then
616  nivcon(ig)=k
617  zcon(ig)=zlev(ig,k)
618  endif
619  enddo
620  enddo
621  if (prt_level.ge.1) print*,'14c OK convect8'
622 !calcul des moments
623 !initialisation
624  do l=1,nlay
625  do ig=1,ngrid
626  q2(ig,l)=0.
627  wth2(ig,l)=0.
628  wth3(ig,l)=0.
629  ratqscth(ig,l)=0.
630  ratqsdiff(ig,l)=0.
631  enddo
632  enddo
633  if (prt_level.ge.1) print*,'14d OK convect8'
634  if (prt_level.ge.10)write(lunout,*) &
635  & 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
636  do l=1,nlay
637  do ig=1,ngrid
638  zf=fraca(ig,l)
639  zf2=zf/(1.-zf)
640 !
641  thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
642  if(zw2(ig,l).gt.1.e-10) then
643  wth2(ig,l)=zf2*(zw2(ig,l))**2
644  else
645  wth2(ig,l)=0.
646  endif
647 ! print*,'wth2=',wth2(ig,l)
648  wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l)) &
649  & *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
650  q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
651 !test: on calcul q2/po=ratqsc
652  ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
653  enddo
654  enddo
655 
656  if (prt_level.ge.10) then
657  print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
658  ig=igout
659  do l=1,nlay
660  print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
661  enddo
662  do l=1,nlay
663  print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
664  enddo
665  endif
666 
667  do ig=1,ngrid
668  alp_int(ig)=0.
669  ale_int(ig)=0.
670  n_int(ig)=0
671  enddo
672 !
673  do l=1,nlay
674  do ig=1,ngrid
675  if(l.LE.lmax(ig)) THEN
676  alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)
677  ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2
678  n_int(ig)=n_int(ig)+1
679  endif
680  enddo
681  enddo
682 ! print*,'avant calcul ale et alp'
683 !calcul de ALE et ALP pour la convection
684  do ig=1,ngrid
685 ! Alp_bl(ig)=0.5*rhobarz(ig,lmix_bis(ig))*wth3(ig,lmix(ig))
686 ! Alp_bl(ig)=0.5*rhobarz(ig,nivcon(ig))*wth3(ig,nivcon(ig))
687 ! Alp_bl(ig)=0.5*rhobarz(ig,lmix(ig))*wth3(ig,lmix(ig))
688 ! & *0.1
689 !valeur integree de alp_bl * 0.5:
690  if (n_int(ig).gt.0) then
691  alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
692 ! if (Alp_bl(ig).lt.0.) then
693 ! Alp_bl(ig)=0.
694  endif
695 ! endif
696 ! write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
697 ! s wth3(ig,nivcon(ig)),Alp_bl(ig)
698 ! write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
699 ! Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
700 ! if (nivcon(ig).eq.1) then
701 ! Ale_bl(ig)=0.
702 ! else
703 !valeur max de ale_bl:
704  ale_bl(ig)=0.5*zw2(ig,lmix(ig))**2
705 ! & /2.
706 ! & *0.1
707 ! Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
708 ! if (n_int(ig).gt.0) then
709 ! Ale_bl(ig)=ale_int(ig)/n_int(ig)
710 ! Ale_bl(ig)=4.
711 ! endif
712 ! endif
713 ! Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
714 ! Ale_bl(ig)=wth2(ig,nivcon(ig))
715 ! write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
716  enddo
717 !test:calcul de la ponderation des couches pour KE
718 !initialisations
719 ! print*,'ponderation'
720  do ig=1,ngrid
721  fm_tot(ig)=0.
722  enddo
723  do ig=1,ngrid
724  do k=1,klev
725  wght_th(ig,k)=1.
726  enddo
727  enddo
728  do ig=1,ngrid
729 ! lalim_conv(ig)=lmix_bis(ig)
730 !la hauteur de la couche alim_conv = hauteur couche alim_therm
731  lalim_conv(ig)=lalim(ig)
732 ! zentr(ig)=zlev(ig,lalim(ig))
733  enddo
734  do ig=1,ngrid
735  do k=1,lalim_conv(ig)
736  fm_tot(ig)=fm_tot(ig)+fm(ig,k)
737  enddo
738  enddo
739  do ig=1,ngrid
740  do k=1,lalim_conv(ig)
741  if (fm_tot(ig).gt.1.e-10) then
742 ! wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
743  endif
744 !on pondere chaque couche par a*
745  if (alim_star(ig,k).gt.1.e-10) then
746  wght_th(ig,k)=alim_star(ig,k)
747  else
748  wght_th(ig,k)=1.
749  endif
750  enddo
751  enddo
752 ! print*,'apres wght_th'
753 !test pour prolonger la convection
754  do ig=1,ngrid
755 !v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
756  if ((alim_star(ig,1).lt.1.e-10)) then
757  lalim_conv(ig)=1
758  wght_th(ig,1)=1.
759 ! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
760  endif
761  enddo
762 
763 !calcul du ratqscdiff
764  if (prt_level.ge.1) print*,'14e OK convect8'
765  var=0.
766  vardiff=0.
767  ratqsdiff(:,:)=0.
768  do ig=1,ngrid
769  do l=1,lalim(ig)
770  var=var+alim_star(ig,l)*zqta(ig,l)*1000.
771  enddo
772  enddo
773  if (prt_level.ge.1) print*,'14f OK convect8'
774  do ig=1,ngrid
775  do l=1,lalim(ig)
776  zf=fraca(ig,l)
777  zf2=zf/(1.-zf)
778  vardiff=vardiff+alim_star(ig,l) &
779  & *(zqta(ig,l)*1000.-var)**2
780 ! ratqsdiff=ratqsdiff+alim_star(ig,l)*
781 ! s (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
782  enddo
783  enddo
784  if (prt_level.ge.1) print*,'14g OK convect8'
785  do l=1,nlay
786  do ig=1,ngrid
787  ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
788 ! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
789  enddo
790  enddo
791 !--------------------------------------------------------------------
792 !
793 !ecriture des fichiers sortie
794 ! print*,'15 OK convect8'
795 
796  if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
797 #ifdef wrgrads_thermcell
798 #include "thermcell_out3d.h"
799 #endif
800 
801  endif
802 
803  if (prt_level.ge.1) print*,'thermcell_main FIN OK'
804 
805 ! if(icount.eq.501) stop'au pas 301 dans thermcell_main'
806  return
807  end
808 
809 !-----------------------------------------------------------------------------
810 
811  subroutine testv0_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
812  IMPLICIT NONE
813 #include "iniprint.h"
814 
815  integer i, k, klon,klev
816  real pplev(klon,klev+1),pplay(klon,klev)
817  real ztv(klon,klev)
818  real po(klon,klev)
819  real ztva(klon,klev)
820  real zqla(klon,klev)
821  real f_star(klon,klev)
822  real zw2(klon,klev)
823  integer long(klon)
824  real seuil
825  character*21 comment
826 
827  if (prt_level.ge.1) THEN
828  print*,'WARNING !!! TEST ',comment
829  endif
830  return
831 
832 ! test sur la hauteur des thermiques ...
833  do i=1,klon
834 !IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
835  if (prt_level.ge.10) then
836  print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
837  print*,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F* W2'
838  do k=1,klev
839  write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
840  enddo
841  endif
842  enddo
843 
844 
845  return
846  end
847 
848 !==============================================================================
849  SUBROUTINE thermcellv0_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
850  & zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
851 
852 !-------------------------------------------------------------------------
853 !thermcell_closure: fermeture, determination de f
854 !-------------------------------------------------------------------------
855  IMPLICIT NONE
856 
857 #include "iniprint.h"
858 #include "thermcell.h"
859  INTEGER ngrid,nlay
860  INTEGER ig,k
861  REAL r_aspect,ptimestep
862  integer lev_out ! niveau pour les print
863 
864  INTEGER lalim(ngrid)
865  REAL alim_star(ngrid,nlay)
866  REAL alim_star_tot(ngrid)
867  REAL rho(ngrid,nlay)
868  REAL zlev(ngrid,nlay)
869  REAL zmax(ngrid),zmax_sec(ngrid)
870  REAL wmax(ngrid),wmax_sec(ngrid)
871  real zdenom
872 
873  REAL alim_star2(ngrid)
874 
875  REAL f(ngrid)
876 
877  character (len=20) :: modname='thermcellV0_main'
878  character (len=80) :: abort_message
879 
880  do ig=1,ngrid
881  alim_star2(ig)=0.
882  enddo
883  do ig=1,ngrid
884  if (alim_star(ig,1).LT.1.e-10) then
885  f(ig)=0.
886  else
887  do k=1,lalim(ig)
888  alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2 &
889  & /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
890  enddo
891  zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig)
892  if (zdenom<1.e-14) then
893  print*,'ig=',ig
894  print*,'alim_star2',alim_star2(ig)
895  print*,'zmax',zmax(ig)
896  print*,'r_aspect',r_aspect
897  print*,'zdenom',zdenom
898  print*,'alim_star',alim_star(ig,:)
899  print*,'zmax_sec',zmax_sec(ig)
900  print*,'wmax_sec',wmax_sec(ig)
901  abort_message = 'zdenom<1.e-14'
902  CALL abort_gcm(modname,abort_message,1)
903  endif
904  if ((zmax_sec(ig).gt.1.e-10).and.(iflag_thermals_ed.eq.0)) then
905  f(ig)=wmax_sec(ig)*alim_star_tot(ig)/(max(500.,zmax_sec(ig))*r_aspect &
906  & *alim_star2(ig))
907 ! f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/ &
908 ! & zmax_sec(ig))*wmax_sec(ig))
909  if(prt_level.GE.10) write(lunout,*)'closure dry',f(ig),wmax_sec(ig),alim_star_tot(ig),zmax_sec(ig)
910  else
911  f(ig)=wmax(ig)*alim_star_tot(ig)/zdenom
912 ! f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/ &
913 ! & zmax(ig))*wmax(ig))
914  if(prt_level.GE.10) print*,'closure moist',f(ig),wmax(ig),alim_star_tot(ig),zmax(ig)
915  endif
916  endif
917 ! f0(ig)=f(ig)
918  enddo
919  if (prt_level.ge.1) print*,'apres fermeture'
920 
921 !
922  return
923  end
924 !==============================================================================
925  SUBROUTINE thermcellv0_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
926  & zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot, &
927  & lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, &
928  & ztla,zqla,zqta,zha,zw2,w_est,zqsatth,lmix,lmix_bis,linter &
929  & ,lev_out,lunout1,igout)
930 
931 !--------------------------------------------------------------------------
932 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
933 !--------------------------------------------------------------------------
934 
935  IMPLICIT NONE
936 
937 #include "YOMCST.h"
938 #include "YOETHF.h"
939 #include "FCTTRE.h"
940 #include "iniprint.h"
941 #include "thermcell.h"
942 
943  INTEGER itap
944  INTEGER lunout1,igout
945  INTEGER ngrid,klev
946  REAL ptimestep
947  REAL ztv(ngrid,klev)
948  REAL zthl(ngrid,klev)
949  REAL po(ngrid,klev)
950  REAL zl(ngrid,klev)
951  REAL rhobarz(ngrid,klev)
952  REAL zlev(ngrid,klev+1)
953  REAL pplev(ngrid,klev+1)
954  REAL pphi(ngrid,klev)
955  REAL zpspsk(ngrid,klev)
956  REAL alim_star(ngrid,klev)
957  REAL zmax_sec(ngrid)
958  REAL f0(ngrid)
959  REAL l_mix
960  REAL r_aspect
961  INTEGER lalim(ngrid)
962  integer lev_out ! niveau pour les print
963  real zcon2(ngrid)
964 
965  real alim_star_tot(ngrid)
966 
967  REAL ztva(ngrid,klev)
968  REAL ztla(ngrid,klev)
969  REAL zqla(ngrid,klev)
970  REAL zqla0(ngrid,klev)
971  REAL zqta(ngrid,klev)
972  REAL zha(ngrid,klev)
973 
974  REAL detr_star(ngrid,klev)
975  REAL coefc
976  REAL detr_stara(ngrid,klev)
977  REAL detr_starb(ngrid,klev)
978  REAL detr_starc(ngrid,klev)
979  REAL detr_star0(ngrid,klev)
980  REAL detr_star1(ngrid,klev)
981  REAL detr_star2(ngrid,klev)
982 
983  REAL entr_star(ngrid,klev)
984  REAL entr_star1(ngrid,klev)
985  REAL entr_star2(ngrid,klev)
986  REAL detr(ngrid,klev)
987  REAL entr(ngrid,klev)
988 
989  REAL zw2(ngrid,klev+1)
990  REAL w_est(ngrid,klev+1)
991  REAL f_star(ngrid,klev+1)
992  REAL wa_moy(ngrid,klev+1)
993 
994  REAL ztva_est(ngrid,klev)
995  REAL zqla_est(ngrid,klev)
996  REAL zqsatth(ngrid,klev)
997  REAL zta_est(ngrid,klev)
998 
999  REAL linter(ngrid)
1000  INTEGER lmix(ngrid)
1001  INTEGER lmix_bis(ngrid)
1002  REAL wmaxa(ngrid)
1003 
1004  INTEGER ig,l,k
1005 
1006  real zcor,zdelta,zcvm5,qlbef
1007  real tbef,qsatbef
1008  real dqsat_dt,dt,num,denom
1009  REAL reps,rlvcp,ddt0
1010  parameter(ddt0=.01)
1011  logical zsat
1012  REAL fact_gamma,fact_epsilon
1013  REAL c2(ngrid,klev)
1014 
1015  zsat=.false.
1016 ! Initialisation
1017  rlvcp = rlvtt/rcpd
1018 
1019  if (iflag_thermals_ed==0) then
1020  fact_gamma=1.
1021  fact_epsilon=1.
1022  else if (iflag_thermals_ed==1) then
1023  fact_gamma=1.
1024  fact_epsilon=1.
1025  else if (iflag_thermals_ed==2) then
1026  fact_gamma=1.
1027  fact_epsilon=2.
1028  endif
1029 
1030  do l=1,klev
1031  do ig=1,ngrid
1032  zqla_est(ig,l)=0.
1033  ztva_est(ig,l)=ztva(ig,l)
1034  zqsatth(ig,l)=0.
1035  enddo
1036  enddo
1037 
1038 !CR: attention test couche alim
1039 ! do l=2,klev
1040 ! do ig=1,ngrid
1041 ! alim_star(ig,l)=0.
1042 ! enddo
1043 ! enddo
1044 !AM:initialisations du thermique
1045  do k=1,klev
1046  do ig=1,ngrid
1047  ztva(ig,k)=ztv(ig,k)
1048  ztla(ig,k)=zthl(ig,k)
1049  zqla(ig,k)=0.
1050  zqta(ig,k)=po(ig,k)
1051 !
1052  ztva(ig,k) = ztla(ig,k)*zpspsk(ig,k)+rlvcp*zqla(ig,k)
1053  ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k)
1054  zha(ig,k) = ztva(ig,k)
1055 !
1056  enddo
1057  enddo
1058  do k=1,klev
1059  do ig=1,ngrid
1060  detr_star(ig,k)=0.
1061  entr_star(ig,k)=0.
1062 
1063  detr_stara(ig,k)=0.
1064  detr_starb(ig,k)=0.
1065  detr_starc(ig,k)=0.
1066  detr_star0(ig,k)=0.
1067  zqla0(ig,k)=0.
1068  detr_star1(ig,k)=0.
1069  detr_star2(ig,k)=0.
1070  entr_star1(ig,k)=0.
1071  entr_star2(ig,k)=0.
1072 
1073  detr(ig,k)=0.
1074  entr(ig,k)=0.
1075  enddo
1076  enddo
1077  if (prt_level.ge.1) print*,'7 OK convect8'
1078  do k=1,klev+1
1079  do ig=1,ngrid
1080  zw2(ig,k)=0.
1081  w_est(ig,k)=0.
1082  f_star(ig,k)=0.
1083  wa_moy(ig,k)=0.
1084  enddo
1085  enddo
1086 
1087  if (prt_level.ge.1) print*,'8 OK convect8'
1088  do ig=1,ngrid
1089  linter(ig)=1.
1090  lmix(ig)=1
1091  lmix_bis(ig)=2
1092  wmaxa(ig)=0.
1093  enddo
1094 
1095 !-----------------------------------------------------------------------------------
1096 !boucle de calcul de la vitesse verticale dans le thermique
1097 !-----------------------------------------------------------------------------------
1098  do l=1,klev-1
1099  do ig=1,ngrid
1100 
1101 
1102 
1103 ! Calcul dans la premiere couche active du thermique (ce qu'on teste
1104 ! en disant que la couche est instable et que w2 en bas de la couche
1105 ! est nulle.
1106 
1107  if (ztv(ig,l).gt.ztv(ig,l+1) &
1108  & .and.alim_star(ig,l).gt.1.e-10 &
1109  & .and.zw2(ig,l).lt.1e-10) then
1110 
1111 
1112 ! Le panache va prendre au debut les caracteristiques de l'air contenu
1113 ! dans cette couche.
1114  ztla(ig,l)=zthl(ig,l)
1115  zqta(ig,l)=po(ig,l)
1116  zqla(ig,l)=zl(ig,l)
1117  f_star(ig,l+1)=alim_star(ig,l)
1118 
1119  zw2(ig,l+1)=2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1) &
1120  & *(zlev(ig,l+1)-zlev(ig,l)) &
1121  & *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
1122  w_est(ig,l+1)=zw2(ig,l+1)
1123 !
1124 
1125 
1126  else if ((zw2(ig,l).ge.1e-10).and. &
1127  & (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
1128 !estimation du detrainement a partir de la geometrie du pas precedent
1129 !tests sur la definition du detr
1130 !calcul de detr_star et entr_star
1131 
1132 
1133 
1134 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1135 ! FH le test miraculeux de Catherine ? Le bout du tunel ?
1136 ! w_est(ig,3)=zw2(ig,2)* &
1137 ! & ((f_star(ig,2))**2) &
1138 ! & /(f_star(ig,2)+alim_star(ig,2))**2+ &
1139 ! & 2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2) &
1140 ! & *(zlev(ig,3)-zlev(ig,2))
1141 ! if (l.gt.2) then
1142 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1143 
1144 
1145 
1146 ! Premier calcul de la vitesse verticale a partir de la temperature
1147 ! potentielle virtuelle
1148 
1149 ! FH CESTQUOI CA ????
1150 #define int1d2
1151 !#undef int1d2
1152 #ifdef int1d2
1153  if (l.ge.2) then
1154 #else
1155  if (l.gt.2) then
1156 #endif
1157 
1158  if (1.eq.1) then
1159  w_est(ig,3)=zw2(ig,2)* &
1160  & ((f_star(ig,2))**2) &
1161  & /(f_star(ig,2)+alim_star(ig,2))**2+ &
1162  & 2.*rg*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
1163 ! & *1./3. &
1164  & *(zlev(ig,3)-zlev(ig,2))
1165  endif
1166 
1167 
1168 !---------------------------------------------------------------------------
1169 !calcul de l entrainement et du detrainement lateral
1170 !---------------------------------------------------------------------------
1171 !
1172 !test:estimation de ztva_new_est sans entrainement
1173 
1174  tbef=ztla(ig,l-1)*zpspsk(ig,l)
1175  zdelta=max(0.,sign(1.,rtt-tbef))
1176  qsatbef= r2es * foeew(tbef,zdelta)/pplev(ig,l)
1177  qsatbef=min(0.5,qsatbef)
1178  zcor=1./(1.-retv*qsatbef)
1179  qsatbef=qsatbef*zcor
1180  zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
1181  if (zsat) then
1182  qlbef=max(0.,zqta(ig,l-1)-qsatbef)
1183  dt = 0.5*rlvcp*qlbef
1184  do while (abs(dt).gt.ddt0)
1185  tbef=tbef+dt
1186  zdelta=max(0.,sign(1.,rtt-tbef))
1187  qsatbef= r2es * foeew(tbef,zdelta)/pplev(ig,l)
1188  qsatbef=min(0.5,qsatbef)
1189  zcor=1./(1.-retv*qsatbef)
1190  qsatbef=qsatbef*zcor
1191  qlbef=zqta(ig,l-1)-qsatbef
1192 
1193  zdelta=max(0.,sign(1.,rtt-tbef))
1194  zcvm5=r5les*(1.-zdelta) + r5ies*zdelta
1195  zcor=1./(1.-retv*qsatbef)
1196  dqsat_dt=foede(tbef,zdelta,zcvm5,qsatbef,zcor)
1197  num=-tbef+ztla(ig,l-1)*zpspsk(ig,l)+rlvcp*qlbef
1198  denom=1.+rlvcp*dqsat_dt
1199  dt=num/denom
1200  enddo
1201  zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef)
1202  endif
1203  ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+rlvcp*zqla_est(ig,l)
1204  ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
1205  zta_est(ig,l)=ztva_est(ig,l)
1206  ztva_est(ig,l) = ztva_est(ig,l)*(1.+retv*(zqta(ig,l-1) &
1207  & -zqla_est(ig,l))-zqla_est(ig,l))
1208 
1209  w_est(ig,l+1)=zw2(ig,l)* &
1210  & ((f_star(ig,l))**2) &
1211  & /(f_star(ig,l)+alim_star(ig,l))**2+ &
1212  & 2.*rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
1213 ! & *1./3. &
1214  & *(zlev(ig,l+1)-zlev(ig,l))
1215  if (w_est(ig,l+1).lt.0.) then
1216  w_est(ig,l+1)=zw2(ig,l)
1217  endif
1218 !
1219 !calcul du detrainement
1220 !=======================
1221 
1222 !CR:on vire les modifs
1223  if (iflag_thermals_ed==0) then
1224 
1225 ! Modifications du calcul du detrainement.
1226 ! Dans la version de la these de Catherine, on passe brusquement
1227 ! de la version seche a la version nuageuse pour le detrainement
1228 ! ce qui peut occasioner des oscillations.
1229 ! dans la nouvelle version, on commence par calculer un detrainement sec.
1230 ! Puis un autre en cas de nuages.
1231 ! Puis on combine les deux lineairement en fonction de la quantite d'eau.
1232 
1233 #define int1d3
1234 !#undef int1d3
1235 #define RIO_TH
1236 #ifdef RIO_TH
1237 !1. Cas non nuageux
1238 ! 1.1 on est sous le zmax_sec et w croit
1239  if ((w_est(ig,l+1).gt.w_est(ig,l)).and. &
1240  & (zlev(ig,l+1).lt.zmax_sec(ig)).and. &
1241 #ifdef int1d3
1242  & (zqla_est(ig,l).lt.1.e-10)) then
1243 #else
1244  & (zqla(ig,l-1).lt.1.e-10)) then
1245 #endif
1246  detr_star(ig,l)=max(0.,(rhobarz(ig,l+1) &
1247  & *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1)) &
1248  & -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l))) &
1249  & /(r_aspect*zmax_sec(ig)))
1250  detr_stara(ig,l)=detr_star(ig,l)
1251 
1252  if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l
1253 
1254 ! 1.2 on est sous le zmax_sec et w decroit
1255  else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and. &
1256 #ifdef int1d3
1257  & (zqla_est(ig,l).lt.1.e-10)) then
1258 #else
1259  & (zqla(ig,l-1).lt.1.e-10)) then
1260 #endif
1261  detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig)) &
1262  & /(rhobarz(ig,lmix(ig))*wmaxa(ig))* &
1263  & (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1)) &
1264  & *((zmax_sec(ig)-zlev(ig,l+1))/ &
1265  & ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2. &
1266  & -rhobarz(ig,l)*sqrt(w_est(ig,l)) &
1267  & *((zmax_sec(ig)-zlev(ig,l))/ &
1268  & ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
1269  detr_starb(ig,l)=detr_star(ig,l)
1270 
1271  if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l
1272 
1273  else
1274 
1275 ! 1.3 dans les autres cas
1276  detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l) &
1277  & *(zlev(ig,l+1)-zlev(ig,l))
1278  detr_starc(ig,l)=detr_star(ig,l)
1279 
1280  if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
1281 
1282  endif
1283 
1284 #else
1285 
1286 ! 1.1 on est sous le zmax_sec et w croit
1287  if ((w_est(ig,l+1).gt.w_est(ig,l)).and. &
1288  & (zlev(ig,l+1).lt.zmax_sec(ig)) ) then
1289  detr_star(ig,l)=max(0.,(rhobarz(ig,l+1) &
1290  & *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1)) &
1291  & -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l))) &
1292  & /(r_aspect*zmax_sec(ig)))
1293 
1294  if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
1295 
1296 ! 1.2 on est sous le zmax_sec et w decroit
1297  else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then
1298  detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig)) &
1299  & /(rhobarz(ig,lmix(ig))*wmaxa(ig))* &
1300  & (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1)) &
1301  & *((zmax_sec(ig)-zlev(ig,l+1))/ &
1302  & ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2. &
1303  & -rhobarz(ig,l)*sqrt(w_est(ig,l)) &
1304  & *((zmax_sec(ig)-zlev(ig,l))/ &
1305  & ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
1306  if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
1307 
1308  else
1309  detr_star=0.
1310  endif
1311 
1312 ! 1.3 dans les autres cas
1313  detr_starc(ig,l)=0.002*f0(ig)*f_star(ig,l) &
1314  & *(zlev(ig,l+1)-zlev(ig,l))
1315 
1316  coefc=min(zqla(ig,l-1)/1.e-3,1.)
1317  if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1.
1318  coefc=1.
1319 ! il semble qu'il soit important de baser le calcul sur
1320 ! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l)
1321  detr_star(ig,l)=detr_starc(ig,l)*coefc+detr_star(ig,l)*(1.-coefc)
1322 
1323  if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l
1324 
1325 #endif
1326 
1327 
1328  if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l
1329 !IM 730508 beg
1330 ! if(itap.GE.7200) THEN
1331 ! print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l)
1332 ! endif
1333 !IM 730508 end
1334 
1335  zqla0(ig,l)=zqla_est(ig,l)
1336  detr_star0(ig,l)=detr_star(ig,l)
1337 !IM 060508 beg
1338 ! if(detr_star(ig,l).GT.1.) THEN
1339 ! print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) &
1340 ! & ,detr_starc(ig,l),coefc
1341 ! endif
1342 !IM 060508 end
1343 !IM 160508 beg
1344 !IM 160508 IF (f0(ig).NE.0.) THEN
1345  detr_star(ig,l)=detr_star(ig,l)/f0(ig)
1346 !IM 160508 ELSE IF(detr_star(ig,l).EQ.0.) THEN
1347 !IM 160508 print*,'WARNING1 : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap
1348 !IM 160508 ELSE
1349 !IM 160508 print*,'WARNING2 : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l)
1350 !IM 160508 ENDIF
1351 !IM 160508 end
1352 !IM 060508 beg
1353 ! if(detr_star(ig,l).GT.1.) THEN
1354 ! print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), &
1355 ! & REAL(1)/f0(ig)
1356 ! endif
1357 !IM 060508 end
1358  if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l
1359 !
1360 !calcul de entr_star
1361 
1362 ! #undef test2
1363 ! #ifdef test2
1364 ! La version test2 destabilise beaucoup le modele.
1365 ! Il semble donc que ca aide d'avoir un entrainement important sous
1366 ! le nuage.
1367 ! if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then
1368 ! entr_star(ig,l)=0.4*detr_star(ig,l)
1369 ! else
1370 ! entr_star(ig,l)=0.
1371 ! endif
1372 ! #else
1373 !
1374 ! Deplacement du calcul de entr_star pour eviter d'avoir aussi
1375 ! entr_star > fstar.
1376 ! Redeplacer suite a la transformation du cas detr>f
1377 ! FH
1378 
1379  if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l
1380 #define int1d
1381 !FH 070508 #define int1d4
1382 !#undef int1d4
1383 ! L'option int1d4 correspond au choix dans le cas ou le detrainement
1384 ! devient trop grand.
1385 
1386 #ifdef int1d
1387 
1388 #ifdef int1d4
1389 #else
1390  detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l))
1391 !FH 070508 plus
1392  detr_star(ig,l)=min(detr_star(ig,l),1.)
1393 #endif
1394 
1395  entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.)
1396 
1397  if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l
1398 #ifdef int1d4
1399 ! Si le detrainement excede le flux en bas + l'entrainement, le thermique
1400 ! doit disparaitre.
1401  if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then
1402  detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l)
1403  f_star(ig,l+1)=0.
1404  linter(ig)=l+1
1405  zw2(ig,l+1)=-1.e-10
1406  endif
1407 #endif
1408 
1409 
1410 #else
1411 
1412  if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l
1413  if(l.gt.lalim(ig)) then
1414  entr_star(ig,l)=0.4*detr_star(ig,l)
1415  else
1416 
1417 ! FH :
1418 ! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
1419 ! en haut de la couche d'alimentation.
1420 ! A remettre en questoin a la premiere occasion mais ca peut aider a
1421 ! ecrire un code robuste.
1422 ! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
1423 ! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
1424 ! d'alimentation, ce qui n'est pas forcement heureux.
1425 
1426  if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l
1427 #undef pre_int1c
1428 #ifdef pre_int1c
1429  entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
1430  detr_star(ig,l)=entr_star(ig,l)
1431 #else
1432  entr_star(ig,l)=0.
1433 #endif
1434 
1435  endif
1436 
1437 #endif
1438 
1439  if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l
1440  entr_star1(ig,l)=entr_star(ig,l)
1441  detr_star1(ig,l)=detr_star(ig,l)
1442 !
1443 
1444 #ifdef int1d
1445 #else
1446  if (detr_star(ig,l).gt.f_star(ig,l)) then
1447 
1448 ! Ce test est là entre autres parce qu'on passe par des valeurs
1449 ! delirantes de detr_star.
1450 ! ca vaut sans doute le coup de verifier pourquoi.
1451 
1452  detr_star(ig,l)=f_star(ig,l)
1453 #ifdef pre_int1c
1454  if (l.gt.lalim(ig)+1) then
1455  entr_star(ig,l)=0.
1456  alim_star(ig,l)=0.
1457 ! FH ajout pour forcer a stoper le thermique juste sous le sommet
1458 ! de la couche (voir calcul de finter)
1459  zw2(ig,l+1)=-1.e-10
1460  linter(ig)=l+1
1461  else
1462  entr_star(ig,l)=0.4*detr_star(ig,l)
1463  endif
1464 #else
1465  entr_star(ig,l)=0.4*detr_star(ig,l)
1466 #endif
1467  endif
1468 #endif
1469 
1470  else !l > 2
1471  detr_star(ig,l)=0.
1472  entr_star(ig,l)=0.
1473  endif
1474 
1475  entr_star2(ig,l)=entr_star(ig,l)
1476  detr_star2(ig,l)=detr_star(ig,l)
1477  if (prt_level.ge.20) print*,'coucou calcul detr 450: ig, l', ig, l
1478 
1479  endif ! iflag_thermals_ed==0
1480 
1481 !CR:nvlle def de entr_star et detr_star
1482  if (iflag_thermals_ed>=1) then
1483 ! if (l.lt.lalim(ig)) then
1484 ! if (l.lt.2) then
1485 ! entr_star(ig,l)=0.
1486 ! detr_star(ig,l)=0.
1487 ! else
1488 ! if (0.001.gt.(RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))/(2.*w_est(ig,l+1)))) then
1489 ! entr_star(ig,l)=0.001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1490 ! else
1491 ! entr_star(ig,l)= &
1492 ! & f_star(ig,l)/(2.*w_est(ig,l+1)) &
1493 ! & *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)) &
1494 ! & *(zlev(ig,l+1)-zlev(ig,l))
1495 
1496 
1497  entr_star(ig,l)=max(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)), &
1498  & f_star(ig,l)/(2.*w_est(ig,l+1)) &
1499  & *rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
1500  & *(zlev(ig,l+1)-zlev(ig,l))) &
1501  & +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1502 
1503  if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
1504  alim_star_tot(ig)=alim_star_tot(ig)+entr_star(ig,l)
1505  lalim(ig)=lmix_bis(ig)
1506  if(prt_level.GE.10) print*,'alim_star_tot',alim_star_tot(ig),entr_star(ig,l)
1507  endif
1508 
1509  if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
1510 ! c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
1511  c2(ig,l)=0.001
1512  detr_star(ig,l)=max(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)), &
1513  & c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
1514  & -f_star(ig,l)/(2.*w_est(ig,l+1)) &
1515  & *rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
1516  & *(zlev(ig,l+1)-zlev(ig,l))) &
1517  & +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1518 
1519  else
1520 ! c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
1521  c2(ig,l)=0.003
1522 
1523  detr_star(ig,l)=max(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)), &
1524  & c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
1525  & -f_star(ig,l)/(2.*w_est(ig,l+1)) &
1526  & *rg*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
1527  & *(zlev(ig,l+1)-zlev(ig,l))) &
1528  & +0.0002*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1529  endif
1530 
1531 
1532 ! detr_star(ig,l)=detr_star(ig,l)*3.
1533 ! if (l.lt.lalim(ig)) then
1534 ! entr_star(ig,l)=0.
1535 ! endif
1536 ! if (l.lt.2) then
1537 ! entr_star(ig,l)=0.
1538 ! detr_star(ig,l)=0.
1539 ! endif
1540 
1541 
1542 ! endif
1543 ! else if ((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10) then
1544 ! entr_star(ig,l)=MAX(0.,0.8*f_star(ig,l)/(2.*w_est(ig,l+1)) &
1545 ! & *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)) &
1546 ! & *(zlev(ig,l+1)-zlev(ig,l))
1547 ! detr_star(ig,l)=0.002*f_star(ig,l) &
1548 ! & *(zlev(ig,l+1)-zlev(ig,l))
1549 ! else
1550 ! entr_star(ig,l)=0.001*f_star(ig,l) &
1551 ! & *(zlev(ig,l+1)-zlev(ig,l))
1552 ! detr_star(ig,l)=MAX(0.,-0.2*f_star(ig,l)/(2.*w_est(ig,l+1)) &
1553 ! & *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)) &
1554 ! & *(zlev(ig,l+1)-zlev(ig,l)) &
1555 ! & +0.002*f_star(ig,l) &
1556 ! & *(zlev(ig,l+1)-zlev(ig,l))
1557 ! endif
1558 
1559  endif ! iflag_thermals_ed==1
1560 
1561 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1562 ! FH inutile si on conserve comme on l'a fait plus haut entr=detr
1563 ! dans la couche d'alimentation
1564 !pas d entrainement dans la couche alim
1565 ! if ((l.le.lalim(ig))) then
1566 ! entr_star(ig,l)=0.
1567 ! endif
1568 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1569 !
1570 !prise en compte du detrainement et de l entrainement dans le calcul du flux
1571 
1572  f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
1573  & -detr_star(ig,l)
1574 
1575 !test sur le signe de f_star
1576  if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l
1577  if (f_star(ig,l+1).gt.1.e-10) then
1578 !----------------------------------------------------------------------------
1579 !calcul de la vitesse verticale en melangeant Tl et qt du thermique
1580 !---------------------------------------------------------------------------
1581 !
1582  zsat=.false.
1583  ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
1584  & (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
1585  & /(f_star(ig,l+1)+detr_star(ig,l))
1586 !
1587  zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
1588  & (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
1589  & /(f_star(ig,l+1)+detr_star(ig,l))
1590 !
1591  tbef=ztla(ig,l)*zpspsk(ig,l)
1592  zdelta=max(0.,sign(1.,rtt-tbef))
1593  qsatbef= r2es * foeew(tbef,zdelta)/pplev(ig,l)
1594  qsatbef=min(0.5,qsatbef)
1595  zcor=1./(1.-retv*qsatbef)
1596  qsatbef=qsatbef*zcor
1597  zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
1598  if (zsat) then
1599  qlbef=max(0.,zqta(ig,l)-qsatbef)
1600  dt = 0.5*rlvcp*qlbef
1601  do while (abs(dt).gt.ddt0)
1602  tbef=tbef+dt
1603  zdelta=max(0.,sign(1.,rtt-tbef))
1604  qsatbef= r2es * foeew(tbef,zdelta)/pplev(ig,l)
1605  qsatbef=min(0.5,qsatbef)
1606  zcor=1./(1.-retv*qsatbef)
1607  qsatbef=qsatbef*zcor
1608  qlbef=zqta(ig,l)-qsatbef
1609 
1610  zdelta=max(0.,sign(1.,rtt-tbef))
1611  zcvm5=r5les*(1.-zdelta) + r5ies*zdelta
1612  zcor=1./(1.-retv*qsatbef)
1613  dqsat_dt=foede(tbef,zdelta,zcvm5,qsatbef,zcor)
1614  num=-tbef+ztla(ig,l)*zpspsk(ig,l)+rlvcp*qlbef
1615  denom=1.+rlvcp*dqsat_dt
1616  dt=num/denom
1617  enddo
1618  zqla(ig,l) = max(0.,qlbef)
1619  endif
1620 !
1621  if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
1622 ! on ecrit de maniere conservative (sat ou non)
1623 ! T = Tl +Lv/Cp ql
1624  ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+rlvcp*zqla(ig,l)
1625  ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
1626 !on rajoute le calcul de zha pour diagnostiques (temp potentielle)
1627  zha(ig,l) = ztva(ig,l)
1628  ztva(ig,l) = ztva(ig,l)*(1.+retv*(zqta(ig,l) &
1629  & -zqla(ig,l))-zqla(ig,l))
1630 
1631 !on ecrit zqsat
1632  zqsatth(ig,l)=qsatbef
1633 !calcul de vitesse
1634  zw2(ig,l+1)=zw2(ig,l)* &
1635  & ((f_star(ig,l))**2) &
1636 ! Tests de Catherine
1637 ! & /(f_star(ig,l+1)+detr_star(ig,l))**2+ &
1638  & /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-fact_epsilon))**2+ &
1639  & 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
1640  & *fact_gamma &
1641  & *(zlev(ig,l+1)-zlev(ig,l))
1642 !prise en compte des forces de pression que qd flottabilité<0
1643 ! zw2(ig,l+1)=zw2(ig,l)* &
1644 ! & 1./(1.+2.*entr_star(ig,l)/f_star(ig,l)) + &
1645 ! & (f_star(ig,l))**2 &
1646 ! & /(f_star(ig,l)+entr_star(ig,l))**2+ &
1647 ! & (f_star(ig,l)-2.*entr_star(ig,l))**2/(f_star(ig,l)+2.*entr_star(ig,l))**2+ &
1648 ! & /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-2.))**2+ &
1649 ! & /(f_star(ig,l)**2+2.*2.*detr_star(ig,l)*f_star(ig,l)+2.*entr_star(ig,l)*f_star(ig,l))+ &
1650 ! & 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
1651 ! & *1./3. &
1652 ! & *(zlev(ig,l+1)-zlev(ig,l))
1653 
1654 ! write(30,*),l+1,zw2(ig,l+1)-zw2(ig,l), &
1655 ! & -2.*entr_star(ig,l)/f_star(ig,l)*zw2(ig,l), &
1656 ! & 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
1657 
1658 
1659 ! zw2(ig,l+1)=zw2(ig,l)* &
1660 ! & (2.-2.*entr_star(ig,l)/f_star(ig,l)) &
1661 ! & -zw2(ig,l-1)+ &
1662 ! & 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
1663 ! & *1./3. &
1664 ! & *(zlev(ig,l+1)-zlev(ig,l))
1665 
1666  endif
1667  endif
1668  if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
1669 !
1670 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
1671 
1672  if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
1673  print*,'On tombe sur le cas particulier de thermcell_plume'
1674  zw2(ig,l+1)=0.
1675  linter(ig)=l+1
1676  endif
1677 
1678 ! if ((zw2(ig,l).gt.0.).and. (zw2(ig,l+1).le.0.)) then
1679  if (zw2(ig,l+1).lt.0.) then
1680  linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
1681  & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1682  zw2(ig,l+1)=0.
1683  endif
1684 
1685  wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1686 
1687  if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1688 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
1689 !on rajoute le calcul de lmix_bis
1690  if (zqla(ig,l).lt.1.e-10) then
1691  lmix_bis(ig)=l+1
1692  endif
1693  lmix(ig)=l+1
1694  wmaxa(ig)=wa_moy(ig,l+1)
1695  endif
1696  enddo
1697  enddo
1698 
1699 !on remplace a* par e* ds premiere couche
1700 ! if (iflag_thermals_ed.ge.1) then
1701 ! do ig=1,ngrid
1702 ! do l=2,klev
1703 ! if (l.lt.lalim(ig)) then
1704 ! alim_star(ig,l)=entr_star(ig,l)
1705 ! endif
1706 ! enddo
1707 ! enddo
1708 ! do ig=1,ngrid
1709 ! lalim(ig)=lmix_bis(ig)
1710 ! enddo
1711 ! endif
1712  if (iflag_thermals_ed.ge.1) then
1713  do ig=1,ngrid
1714  do l=2,lalim(ig)
1715  alim_star(ig,l)=entr_star(ig,l)
1716  entr_star(ig,l)=0.
1717  enddo
1718  enddo
1719  endif
1720  if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
1721 
1722 ! print*,'thermcell_plume OK'
1723 
1724  return
1725  end
1726 !==============================================================================
1727  SUBROUTINE thermcellv0_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
1728  & lalim,lmin,zmax,wmax,lev_out)
1729 
1730 !--------------------------------------------------------------------------
1731 !thermcell_dry: calcul de zmax et wmax du thermique sec
1732 !--------------------------------------------------------------------------
1733  IMPLICIT NONE
1734 #include "YOMCST.h"
1735 #include "iniprint.h"
1736  INTEGER l,ig
1737 
1738  INTEGER ngrid,nlay
1739  REAL zlev(ngrid,nlay+1)
1740  REAL pphi(ngrid,nlay)
1741  REAl ztv(ngrid,nlay)
1742  REAL alim_star(ngrid,nlay)
1743  INTEGER lalim(ngrid)
1744  integer lev_out ! niveau pour les print
1745 
1746  REAL zmax(ngrid)
1747  REAL wmax(ngrid)
1748 
1749 !variables locales
1750  REAL zw2(ngrid,nlay+1)
1751  REAL f_star(ngrid,nlay+1)
1752  REAL ztva(ngrid,nlay+1)
1753  REAL wmaxa(ngrid)
1754  REAL wa_moy(ngrid,nlay+1)
1755  REAL linter(ngrid),zlevinter(ngrid)
1756  INTEGER lmix(ngrid),lmax(ngrid),lmin(ngrid)
1757 
1758 !initialisations
1759  do ig=1,ngrid
1760  do l=1,nlay+1
1761  zw2(ig,l)=0.
1762  wa_moy(ig,l)=0.
1763  enddo
1764  enddo
1765  do ig=1,ngrid
1766  do l=1,nlay
1767  ztva(ig,l)=ztv(ig,l)
1768  enddo
1769  enddo
1770  do ig=1,ngrid
1771  wmax(ig)=0.
1772  wmaxa(ig)=0.
1773  enddo
1774 !calcul de la vitesse a partir de la CAPE en melangeant thetav
1775 
1776 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1777 ! A eliminer
1778 ! Ce if complique etait fait pour reperer la premiere couche instable
1779 ! Ici, c'est lmin.
1780 !
1781 ! do l=1,nlay-2
1782 ! do ig=1,ngrid
1783 ! if (ztv(ig,l).gt.ztv(ig,l+1) &
1784 ! & .and.alim_star(ig,l).gt.1.e-10 &
1785 ! & .and.zw2(ig,l).lt.1e-10) then
1786 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1787 
1788 
1789 ! Calcul des F^*, integrale verticale de E^*
1790  f_star(:,1)=0.
1791  do l=1,nlay
1792  f_star(:,l+1)=f_star(:,l)+alim_star(:,l)
1793  enddo
1794 
1795 ! niveau (reel) auquel zw2 s'annule FH :n'etait pas initialise
1796  linter(:)=0.
1797 
1798 ! couche la plus haute concernee par le thermique.
1799  lmax(:)=1
1800 
1801 ! Le niveau linter est une variable continue qui se trouve dans la couche
1802 ! lmax
1803 
1804  do l=1,nlay-2
1805  do ig=1,ngrid
1806  if (l.eq.lmin(ig).and.lalim(ig).gt.1) then
1807 
1808 !------------------------------------------------------------------------
1809 ! Calcul de la vitesse en haut de la premiere couche instable.
1810 ! Premiere couche du panache thermique
1811 !------------------------------------------------------------------------
1812  zw2(ig,l+1)=2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1) &
1813  & *(zlev(ig,l+1)-zlev(ig,l)) &
1814  & *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
1815 
1816 !------------------------------------------------------------------------
1817 ! Tant que la vitesse en bas de la couche et la somme du flux de masse
1818 ! et de l'entrainement (c'est a dire le flux de masse en haut) sont
1819 ! positifs, on calcul
1820 ! 1. le flux de masse en haut f_star(ig,l+1)
1821 ! 2. la temperature potentielle virtuelle dans la couche ztva(ig,l)
1822 ! 3. la vitesse au carré en haut zw2(ig,l+1)
1823 !------------------------------------------------------------------------
1824 
1825 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1826 ! A eliminer : dans cette version, si zw2 est > 0 on a un therique.
1827 ! et donc, au dessus, f_star(ig,l+1) est forcement suffisamment
1828 ! grand puisque on n'a pas de detrainement.
1829 ! f_star est une fonction croissante.
1830 ! c'est donc vraiment sur zw2 uniquement qu'il faut faire le test.
1831 ! else if ((zw2(ig,l).ge.1e-10).and. &
1832 ! & (f_star(ig,l)+alim_star(ig,l).gt.1.e-10)) then
1833 ! f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)
1834 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1835 
1836  else if (zw2(ig,l).ge.1e-10) then
1837 
1838  ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+alim_star(ig,l) &
1839  & *ztv(ig,l))/f_star(ig,l+1)
1840  zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+ &
1841  & 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
1842  & *(zlev(ig,l+1)-zlev(ig,l))
1843  endif
1844 ! determination de zmax continu par interpolation lineaire
1845 !------------------------------------------------------------------------
1846 
1847  if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
1848 ! print*,'On tombe sur le cas particulier de thermcell_dry'
1849  zw2(ig,l+1)=0.
1850  linter(ig)=l+1
1851  lmax(ig)=l
1852  endif
1853 
1854  if (zw2(ig,l+1).lt.0.) then
1855  linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
1856  & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
1857  zw2(ig,l+1)=0.
1858  lmax(ig)=l
1859  endif
1860 
1861  wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
1862 
1863  if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
1864 ! lmix est le niveau de la couche ou w (wa_moy) est maximum
1865  lmix(ig)=l+1
1866  wmaxa(ig)=wa_moy(ig,l+1)
1867  endif
1868  enddo
1869  enddo
1870  if (prt_level.ge.1) print*,'fin calcul zw2'
1871 !
1872 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1873 ! A eliminer :
1874 ! Ce calcul de lmax est fait en meme temps que celui de linter, plus haut
1875 ! Calcul de la couche correspondant a la hauteur du thermique
1876 ! do ig=1,ngrid
1877 ! lmax(ig)=lalim(ig)
1878 ! enddo
1879 ! do ig=1,ngrid
1880 ! do l=nlay,lalim(ig)+1,-1
1881 ! if (zw2(ig,l).le.1.e-10) then
1882 ! lmax(ig)=l-1
1883 ! endif
1884 ! enddo
1885 ! enddo
1886 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1887 
1888 !
1889 ! Determination de zw2 max
1890  do ig=1,ngrid
1891  wmax(ig)=0.
1892  enddo
1893 
1894  do l=1,nlay
1895  do ig=1,ngrid
1896  if (l.le.lmax(ig)) then
1897  zw2(ig,l)=sqrt(zw2(ig,l))
1898  wmax(ig)=max(wmax(ig),zw2(ig,l))
1899  else
1900  zw2(ig,l)=0.
1901  endif
1902  enddo
1903  enddo
1904 
1905 ! Longueur caracteristique correspondant a la hauteur des thermiques.
1906  do ig=1,ngrid
1907  zmax(ig)=0.
1908  zlevinter(ig)=zlev(ig,1)
1909  enddo
1910  do ig=1,ngrid
1911 ! calcul de zlevinter
1912 
1913 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1914 ! FH A eliminer
1915 ! Simplification
1916 ! zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))* &
1917 ! & linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1) &
1918 ! & -zlev(ig,lmax(ig)))
1919 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1920 
1921  zlevinter(ig)=zlev(ig,lmax(ig)) + &
1922  & (linter(ig)-lmax(ig))*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
1923  zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
1924  enddo
1925 
1926 ! Verification que lalim<=lmax
1927  do ig=1,ngrid
1928  if(lalim(ig)>lmax(ig)) then
1929  if ( prt_level > 1 ) THEN
1930  print*,'WARNING thermcell_dry ig=',ig,' lalim=',lalim(ig),' lmax(ig)=',lmax(ig)
1931  endif
1932  lmax(ig)=lalim(ig)
1933  endif
1934  enddo
1935 
1936  RETURN
1937  END
1938 !==============================================================================
1939  SUBROUTINE thermcellv0_init(ngrid,nlay,ztv,zlay,zlev, &
1940  & lalim,lmin,alim_star,alim_star_tot,lev_out)
1941 
1942 !----------------------------------------------------------------------
1943 !thermcell_init: calcul du profil d alimentation du thermique
1944 !----------------------------------------------------------------------
1945  IMPLICIT NONE
1946 #include "iniprint.h"
1947 #include "thermcell.h"
1948 
1949  INTEGER l,ig
1950 !arguments d entree
1951  INTEGER ngrid,nlay
1952  REAL ztv(ngrid,nlay)
1953  REAL zlay(ngrid,nlay)
1954  REAL zlev(ngrid,nlay+1)
1955 !arguments de sortie
1956  INTEGER lalim(ngrid)
1957  INTEGER lmin(ngrid)
1958  REAL alim_star(ngrid,nlay)
1959  REAL alim_star_tot(ngrid)
1960  integer lev_out ! niveau pour les print
1961 
1962  REAL zzalim(ngrid)
1963 !CR: ponderation entrainement des couches instables
1964 !def des alim_star tels que alim=f*alim_star
1965 
1966  do l=1,nlay
1967  do ig=1,ngrid
1968  alim_star(ig,l)=0.
1969  enddo
1970  enddo
1971 ! determination de la longueur de la couche d entrainement
1972  do ig=1,ngrid
1973  lalim(ig)=1
1974  enddo
1975 
1976  if (iflag_thermals_ed.ge.1) then
1977 !si la première couche est instable, on declenche un thermique
1978  do ig=1,ngrid
1979  if (ztv(ig,1).gt.ztv(ig,2)) then
1980  lmin(ig)=1
1981  lalim(ig)=2
1982  alim_star(ig,1)=1.
1983  alim_star_tot(ig)=alim_star(ig,1)
1984  if(prt_level.GE.10) print*,'init',alim_star(ig,1),alim_star_tot(ig)
1985  else
1986  lmin(ig)=1
1987  lalim(ig)=1
1988  alim_star(ig,1)=0.
1989  alim_star_tot(ig)=0.
1990  endif
1991  enddo
1992 
1993  else
1994 !else iflag_thermals_ed=0 ancienne def de l alim
1995 
1996 !on ne considere que les premieres couches instables
1997  do l=nlay-2,1,-1
1998  do ig=1,ngrid
1999  if (ztv(ig,l).gt.ztv(ig,l+1).and. &
2000  & ztv(ig,l+1).le.ztv(ig,l+2)) then
2001  lalim(ig)=l+1
2002  endif
2003  enddo
2004  enddo
2005 
2006 ! determination du lmin: couche d ou provient le thermique
2007 
2008  do ig=1,ngrid
2009 ! FH initialisation de lmin a nlay plutot que 1.
2010 ! lmin(ig)=nlay
2011  lmin(ig)=1
2012  enddo
2013  do l=nlay,2,-1
2014  do ig=1,ngrid
2015  if (ztv(ig,l-1).gt.ztv(ig,l)) then
2016  lmin(ig)=l-1
2017  endif
2018  enddo
2019  enddo
2020 !
2021  zzalim(:)=0.
2022  do l=1,nlay-1
2023  do ig=1,ngrid
2024  if (l<lalim(ig)) then
2025  zzalim(ig)=zzalim(ig)+zlay(ig,l)*(ztv(ig,l)-ztv(ig,l+1))
2026  endif
2027  enddo
2028  enddo
2029  do ig=1,ngrid
2030  if (lalim(ig)>1) then
2031  zzalim(ig)=zlay(ig,1)+zzalim(ig)/(ztv(ig,1)-ztv(ig,lalim(ig)))
2032  else
2033  zzalim(ig)=zlay(ig,1)
2034  endif
2035  enddo
2036 
2037  if(prt_level.GE.10) print*,'ZZALIM LALIM ',zzalim,lalim,zlay(1,lalim(1))
2038 
2039 ! definition de l'entrainement des couches
2040  if (1.eq.1) then
2041  do l=1,nlay-1
2042  do ig=1,ngrid
2043  if (ztv(ig,l).gt.ztv(ig,l+1).and. &
2044  & l.ge.lmin(ig).and.l.lt.lalim(ig)) then
2045 !def possibles pour alim_star: zdthetadz, dthetadz, zdtheta
2046  alim_star(ig,l)=max((ztv(ig,l)-ztv(ig,l+1)),0.) &
2047  & *sqrt(zlev(ig,l+1))
2048  endif
2049  enddo
2050  enddo
2051  else
2052  do l=1,nlay-1
2053  do ig=1,ngrid
2054  if (ztv(ig,l).gt.ztv(ig,l+1).and. &
2055  & l.ge.lmin(ig).and.l.lt.lalim(ig)) then
2056  alim_star(ig,l)=max(3.*zzalim(ig)-zlay(ig,l),0.) &
2057  & *(zlev(ig,l+1)-zlev(ig,l))
2058  endif
2059  enddo
2060  enddo
2061  endif
2062 
2063 ! pas de thermique si couche 1 stable
2064  do ig=1,ngrid
2065 !CRnouveau test
2066  if (alim_star(ig,1).lt.1.e-10) then
2067  do l=1,nlay
2068  alim_star(ig,l)=0.
2069  enddo
2070  lmin(ig)=1
2071  endif
2072  enddo
2073 ! calcul de l alimentation totale
2074  do ig=1,ngrid
2075  alim_star_tot(ig)=0.
2076  enddo
2077  do l=1,nlay
2078  do ig=1,ngrid
2079  alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
2080  enddo
2081  enddo
2082 !
2083 ! Calcul entrainement normalise
2084  do l=1,nlay
2085  do ig=1,ngrid
2086  if (alim_star_tot(ig).gt.1.e-10) then
2087  alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
2088  endif
2089  enddo
2090  enddo
2091 
2092 !on remet alim_star_tot a 1
2093  do ig=1,ngrid
2094  alim_star_tot(ig)=1.
2095  enddo
2096 
2097  endif
2098 !endif iflag_thermals_ed
2099  return
2100  end
2101 !==============================================================================