Context Navigation

thermcellV0_main.F90 @ 2056

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

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

Download in other formats:

Original Format