Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

thermcellV0_main.F90 @ 3316

Last change on this file since 3316 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

IPSLCM5A configuration (comparison between trunk and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
- numerical convergence on dynamical fields over 3 days
- start files are equivalent (except for RN and PB fields)
- daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

IPSLCM5A configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

MH07 configuration:
- start files are equivalent (except for RN and PB fields)
- monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

pour une config IPSLCM5A (comparaison tronc/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux (à part sur RN et Pb)
- fichiers histoire égaux

pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
- convergence numérique sur les champs dynamiques après 3 jours
- restart et restartphy égaux
- fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

pour une config IPSLCM5A:
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

pour une config MH07
- restart et restartphy égaux (à part sur RN et PB)
- pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

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

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

Download in other formats:

Original Format