Context Navigation

thermcellV0_main.F90 @ 3818

Last change on this file since 3818 was 3818, checked in by millour, 10 years ago
Some partial cleanup on uses of "dimensions.h" in physics. At this point 3D gcm compiles and bench seems to run fine :-) EM
File size: 68.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format