Context Navigation

thermcellV0_main.F90 @ 1643

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

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

Download in other formats:

Original Format