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