Context Navigation

thermcellV0_main.F90 @ 5446

Last change on this file since 5446 was 2351, checked in by Ehouarn Millour, 9 years ago

More on physics/dynamics separation and cleanup:

Set things up so that all physics-related initializations are done via iniphysiq.
Created a "geometry_mod.F90" module in phy_common to store information on the loacl grid (i.e. replaces comgeomphy) and moreover give these variables more obvious names (e.g.: rlond => longitude, rlatd => latitude, airephy => cell_area).
removed obsolete comgeomphy.h and comgeomphy.F90

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

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

Download in other formats:

Original Format