Context Navigation

thermcellV0_main.F90 @ 5451

Last change on this file since 5451 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

Rev	Line
[1294]	1	!
[1299]	2	! $Id$
[1294]	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
[2311]	15	USE print_control_mod, ONLY: prt_level,lunout
[1294]	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
[1861]	43	include "YOMCST.h"
	44	include "YOETHF.h"
	45	include "FCTTRE.h"
[1294]	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
[1371]	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
[1294]	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
[1299]	174	character (len=20) :: modname='thermcellV0_main'
	175	character (len=80) :: abort_message
	176
[1294]	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
[1299]	464	abort_message = 'Dans thermcell_main f0(1).lt.0 '
[2311]	465	CALL abort_physic (modname,abort_message,1)
[1294]	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
[1738]	503	call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse, &
[1294]	504	& zthl,zdthladj,zta,lev_out)
[1738]	505	call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse, &
[1294]	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
[1738]	545	call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse &
[1294]	546	& ,zu,pduadj,zua,lev_out)
[1738]	547	call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse &
[1294]	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
[1372]	577	!IM do k=1,nlay
	578	do k=1,nlay-1
[1294]	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
[1372]	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 '
[2311]	591	CALL abort_physic (modname,abort_message,1)
[1372]	592	endif
	593	enddo
[1294]	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
[1330]	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
[1294]	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)
[2311]	790	USE print_control_mod, ONLY: prt_level
[1294]	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	!-------------------------------------------------------------------------
[2311]	833	USE print_control_mod, ONLY: prt_level,lunout
[1294]	834	IMPLICIT NONE
	835
[1861]	836	include "thermcell.h"
[1294]	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
[1299]	855	character (len=20) :: modname='thermcellV0_main'
	856	character (len=80) :: abort_message
	857
[1294]	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)
[1299]	879	abort_message = 'zdenom<1.e-14'
[2311]	880	CALL abort_physic (modname,abort_message,1)
[1294]	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
[2311]	913	USE print_control_mod, ONLY: prt_level
[1294]	914	IMPLICIT NONE
	915
[1861]	916	include "YOMCST.h"
	917	include "YOETHF.h"
	918	include "FCTTRE.h"
	919	include "thermcell.h"
[1294]	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), &
[1299]	1333	! & REAL(1)/f0(ig)
[1294]	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	!--------------------------------------------------------------------------
[2311]	1711	USE print_control_mod, ONLY: prt_level
[1294]	1712	IMPLICIT NONE
[1861]	1713	include "YOMCST.h"
[1294]	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	!----------------------------------------------------------------------
[2311]	1923	USE print_control_mod, ONLY: prt_level
[1294]	1924	IMPLICIT NONE
[1861]	1925	include "thermcell.h"
[1294]	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