Context Navigation

thermcellV0_main.F90 @ 2179

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

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

Download in other formats:

Original Format