Context Navigation

thermcell_plume_6A.F90 @ 4440

Last change on this file since 4440 was 4302, checked in by fhourdin, 2 years ago
Initialisations à 0 pour la fonctionalité replay
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 Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 40.5 KB

Rev	Line
[1403]	1	!
	2	! $Id: thermcell_plume_6A.F90 4302 2022-10-21 17:16:36Z fhourdin $
	3	!
[4094]	4	SUBROUTINE thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
[1403]	5	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	6	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	7	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
[2267]	8	& ,lev_out,lunout1,igout)
	9	! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
[878]	10	!--------------------------------------------------------------------------
[1998]	11	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
[878]	12	!--------------------------------------------------------------------------
	13
[4089]	14	USE thermcell_ini_mod, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
	15	USE thermcell_ini_mod, ONLY: fact_epsilon, betalpha, afact, fact_shell
	16	USE thermcell_ini_mod, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
	17	USE thermcell_ini_mod, ONLY: mix0, thermals_flag_alim
	18
[1968]	19	IMPLICIT NONE
[878]	20
[4094]	21	integer,intent(in) :: itap,lev_out,lunout1,igout,ngrid,nlay
[4089]	22	real,intent(in) :: ptimestep
[4094]	23	real,intent(in),dimension(ngrid,nlay) :: ztv
	24	real,intent(in),dimension(ngrid,nlay) :: zthl
	25	real,intent(in),dimension(ngrid,nlay) :: po
	26	real,intent(in),dimension(ngrid,nlay) :: zl
	27	real,intent(in),dimension(ngrid,nlay) :: rhobarz
	28	real,intent(in),dimension(ngrid,nlay+1) :: zlev
	29	real,intent(in),dimension(ngrid,nlay+1) :: pplev
	30	real,intent(in),dimension(ngrid,nlay) :: pphi
	31	real,intent(in),dimension(ngrid,nlay) :: zpspsk
[4089]	32	real,intent(in),dimension(ngrid) :: f0
[878]	33
[4089]	34	integer,intent(out) :: lalim(ngrid)
[4094]	35	real,intent(out),dimension(ngrid,nlay) :: alim_star
[4089]	36	real,intent(out),dimension(ngrid) :: alim_star_tot
[4094]	37	real,intent(out),dimension(ngrid,nlay) :: detr_star
	38	real,intent(out),dimension(ngrid,nlay) :: entr_star
	39	real,intent(out),dimension(ngrid,nlay+1) :: f_star
	40	real,intent(out),dimension(ngrid,nlay) :: csc
	41	real,intent(out),dimension(ngrid,nlay) :: ztva
	42	real,intent(out),dimension(ngrid,nlay) :: ztla
	43	real,intent(out),dimension(ngrid,nlay) :: zqla
	44	real,intent(out),dimension(ngrid,nlay) :: zqta
	45	real,intent(out),dimension(ngrid,nlay) :: zha
	46	real,intent(out),dimension(ngrid,nlay+1) :: zw2
	47	real,intent(out),dimension(ngrid,nlay+1) :: w_est
	48	real,intent(out),dimension(ngrid,nlay) :: ztva_est
	49	real,intent(out),dimension(ngrid,nlay) :: zqsatth
	50	integer,intent(out),dimension(ngrid) :: lmix
	51	integer,intent(out),dimension(ngrid) :: lmix_bis
	52	real,intent(out),dimension(ngrid) :: linter
[878]	53
[1968]	54	REAL zdw2,zdw2bis
[1403]	55	REAL zw2modif
[1968]	56	REAL zw2fact,zw2factbis
[4161]	57	REAL,dimension(ngrid,nlay) :: zeps
[878]	58
[4161]	59	REAL, dimension(ngrid) :: wmaxa(ngrid)
[878]	60
[2140]	61	INTEGER ig,l,k,lt,it,lm
[4089]	62	integer nbpb
[878]	63
[4094]	64	real,dimension(ngrid,nlay) :: detr
	65	real,dimension(ngrid,nlay) :: entr
	66	real,dimension(ngrid,nlay+1) :: wa_moy
	67	real,dimension(ngrid,nlay) :: ztv_est
[4089]	68	real,dimension(ngrid) :: ztemp,zqsat
[4094]	69	real,dimension(ngrid,nlay) :: zqla_est
	70	real,dimension(ngrid,nlay) :: zta_est
[4089]	71
[4161]	72	real,dimension(ngrid,nlay) :: zbuoy,gamma,zdqt
	73	real zdz,zalpha,zw2m
	74	real,dimension(ngrid,nlay) :: zbuoyjam,zdqtjam
[1968]	75	real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
[4161]	76	real, dimension(ngrid) :: d_temp
[2046]	77	real ztv1,ztv2,factinv,zinv,zlmel
[2106]	78	real zlmelup,zlmeldwn,zlt,zltdwn,zltup
	79	real atv1,atv2,btv1,btv2
[2046]	80	real ztv_est1,ztv_est2
[878]	81	real zcor,zdelta,zcvm5,qlbef
[2267]	82	real zbetalpha, coefzlmel
[2106]	83	real eps
[878]	84	logical Zsat
[4161]	85	LOGICAL,dimension(ngrid) :: active,activetmp
[2106]	86	REAL fact_gamma,fact_gamma2,fact_epsilon2
[4089]	87	REAL coefc
[4161]	88	REAL,dimension(ngrid,nlay) :: c2
[2106]	89
	90	if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11'
[878]	91	Zsat=.false.
	92	! Initialisation
[1968]	93
[2406]	94
[1968]	95	zbetalpha=betalpha/(1.+betalpha)
	96
	97
	98	! Initialisations des variables r?elles
[1403]	99	if (1==1) then
	100	ztva(:,:)=ztv(:,:)
	101	ztva_est(:,:)=ztva(:,:)
[1968]	102	ztv_est(:,:)=ztv(:,:)
[1403]	103	ztla(:,:)=zthl(:,:)
	104	zqta(:,:)=po(:,:)
[1968]	105	zqla(:,:)=0.
[1403]	106	zha(:,:) = ztva(:,:)
	107	else
	108	ztva(:,:)=0.
[1968]	109	ztv_est(:,:)=0.
[1403]	110	ztva_est(:,:)=0.
	111	ztla(:,:)=0.
	112	zqta(:,:)=0.
	113	zha(:,:) =0.
	114	endif
[878]	115
[1403]	116	zqla_est(:,:)=0.
	117	zqsatth(:,:)=0.
	118	zqla(:,:)=0.
	119	detr_star(:,:)=0.
	120	entr_star(:,:)=0.
	121	alim_star(:,:)=0.
	122	alim_star_tot(:)=0.
	123	csc(:,:)=0.
	124	detr(:,:)=0.
	125	entr(:,:)=0.
	126	zw2(:,:)=0.
[1968]	127	zbuoy(:,:)=0.
	128	zbuoyjam(:,:)=0.
	129	gamma(:,:)=0.
	130	zeps(:,:)=0.
[1403]	131	w_est(:,:)=0.
	132	f_star(:,:)=0.
	133	wa_moy(:,:)=0.
	134	linter(:)=1.
[1968]	135	! linter(:)=1.
[1403]	136	! Initialisation des variables entieres
	137	lmix(:)=1
	138	lmix_bis(:)=2
	139	wmaxa(:)=0.
[972]	140
[4302]	141	! Initialisation a 0 en cas de sortie dans replay
	142	zqsat(:)=0.
	143	zta_est(:,:)=0.
	144	zdqt(:,:)=0.
	145	zdqtjam(:,:)=0.
	146	c2(:,:)=0.
[1968]	147
[4302]	148
[1403]	149	!-------------------------------------------------------------------------
	150	! On ne considere comme actif que les colonnes dont les deux premieres
	151	! couches sont instables.
	152	!-------------------------------------------------------------------------
[2387]	153
[1403]	154	active(:)=ztv(:,1)>ztv(:,2)
[2387]	155	d_temp(:)=0. ! Pour activer un contraste de temperature a la base
	156	! du panache
	157	! Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
[4094]	158	CALL thermcell_alim(thermals_flag_alim,ngrid,nlay,ztv,d_temp,zlev,alim_star,lalim)
[1403]	159
	160	!------------------------------------------------------------------------------
	161	! Calcul dans la premiere couche
	162	! On decide dans cette version que le thermique n'est actif que si la premiere
	163	! couche est instable.
[1968]	164	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
[1403]	165	! dans une couche l>1
	166	!------------------------------------------------------------------------------
	167	do ig=1,ngrid
[972]	168	! Le panache va prendre au debut les caracteristiques de l'air contenu
	169	! dans cette couche.
[1403]	170	if (active(ig)) then
	171	ztla(ig,1)=zthl(ig,1)
	172	zqta(ig,1)=po(ig,1)
	173	zqla(ig,1)=zl(ig,1)
	174	!cr: attention, prise en compte de f*(1)=1
	175	f_star(ig,2)=alim_star(ig,1)
	176	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
	177	& *(zlev(ig,2)-zlev(ig,1)) &
	178	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
	179	w_est(ig,2)=zw2(ig,2)
	180	endif
	181	enddo
[878]	182	!
[972]	183
[1403]	184	!==============================================================================
	185	!boucle de calcul de la vitesse verticale dans le thermique
	186	!==============================================================================
[4094]	187	do l=2,nlay-1
[1403]	188	!==============================================================================
[972]	189
[878]	190
[1403]	191	! On decide si le thermique est encore actif ou non
	192	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
	193	do ig=1,ngrid
	194	active(ig)=active(ig) &
	195	& .and. zw2(ig,l)>1.e-10 &
	196	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
	197	enddo
[878]	198
	199
	200
	201	!---------------------------------------------------------------------------
[1403]	202	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
	203	! sans tenir compte du detrainement et de l'entrainement dans cette
	204	! couche
[1968]	205	! C'est a dire qu'on suppose
	206	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
[1403]	207	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
	208	! avant) a l'alimentation pour avoir un calcul plus propre
[878]	209	!---------------------------------------------------------------------------
[972]	210
[1968]	211	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
	212	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
	213	do ig=1,ngrid
[2046]	214	! print*,'active',active(ig),ig,l
[1968]	215	if(active(ig)) then
	216	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
[878]	217	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
[1403]	218	zta_est(ig,l)=ztva_est(ig,l)
[878]	219	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
	220	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
	221	& -zqla_est(ig,l))-zqla_est(ig,l))
[1968]	222
[878]	223
[2140]	224	!Modif AJAM
[1403]	225
[2140]	226	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	227	zdz=zlev(ig,l+1)-zlev(ig,l)
[2106]	228	lmel=fact_thermals_ed_dz*zlev(ig,l)
	229	! lmel=0.09*zlev(ig,l)
[2046]	230	zlmel=zlev(ig,l)+lmel
[2106]	231	zlmelup=zlmel+(zdz/2)
	232	zlmeldwn=zlmel-(zdz/2)
	233
[1968]	234	lt=l+1
[2106]	235	zlt=zlev(ig,lt)
	236	zdz3=zlev(ig,lt+1)-zlt
	237	zltdwn=zlt-zdz3/2
	238	zltup=zlt+zdz3/2
	239
	240	!=========================================================================
	241	! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
	242	!=========================================================================
	243
[2046]	244	!--------------------------------------------------
[2106]	245	if (iflag_thermals_ed.lt.8) then
	246	!--------------------------------------------------
	247	!AJ052014: J'ai remplac?? la boucle do par un do while
[2046]	248	! afin de faire moins de calcul dans la boucle
	249	!--------------------------------------------------
[2106]	250	do while (zlmelup.gt.zltup)
	251	lt=lt+1
	252	zlt=zlev(ig,lt)
	253	zdz3=zlev(ig,lt+1)-zlt
	254	zltdwn=zlt-zdz3/2
	255	zltup=zlt+zdz3/2
	256	enddo
[2046]	257	!--------------------------------------------------
	258	!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
[2106]	259	! on cherche o?? se trouve l'altitude d'inversion
[2046]	260	! en calculant ztv1 (interpolation de la valeur de
	261	! theta au niveau lt en utilisant les niveaux lt-1 et
	262	! lt-2) et ztv2 (interpolation avec les niveaux lt+1
[2106]	263	! et lt+2). Si theta r??ellement calcul??e au niveau lt
[2046]	264	! comprise entre ztv1 et ztv2, alors il y a inversion
	265	! et on calcule son altitude zinv en supposant que ztv(lt)
[2106]	266	! est une combinaison lineaire de ztv1 et ztv2.
	267	! Ensuite, on calcule la flottabilite en comparant
	268	! la temperature de la couche l a celle de l'air situe
	269	! l+lmel plus haut, ce qui necessite de savoir quel fraction
[2046]	270	! de cet air est au-dessus ou en-dessous de l'inversion
	271	!--------------------------------------------------
[2106]	272	atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2))
	273	btv1=(ztv(ig,lt-2)zlev(ig,lt-1)-ztv(ig,lt-1)zlev(ig,lt-2)) &
[2046]	274	& /(zlev(ig,lt-1)-zlev(ig,lt-2))
[2106]	275	atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1))
	276	btv2=(ztv(ig,lt+1)zlev(ig,lt+2)-ztv(ig,lt+2)zlev(ig,lt+1)) &
[2046]	277	& /(zlev(ig,lt+2)-zlev(ig,lt+1))
	278
[2106]	279	ztv1=atv1*zlt+btv1
	280	ztv2=atv2*zlt+btv2
[2046]	281
[2106]	282	if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
[2046]	283
[2106]	284	!--------------------------------------------------
	285	!AJ052014: D??calage de zinv qui est entre le haut
	286	! et le bas de la couche lt
	287	!--------------------------------------------------
	288	factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
	289	zinv=zltdwn+zdz3*factinv
[2046]	290
[2106]	291
	292	if (zlmeldwn.ge.zinv) then
	293	ztv_est(ig,l)=atv2*zlmel+btv2
	294	zbuoyjam(ig,l)=fact_shellRG(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
	295	& +(1.-fact_shell)*zbuoy(ig,l)
	296	elseif (zlmelup.ge.zinv) then
	297	ztv_est2=atv20.5(zlmelup+zinv)+btv2
	298	ztv_est1=atv10.5(zinv+zlmeldwn)+btv1
	299	ztv_est(ig,l)=((zlmelup-zinv)/zdz)ztv_est2+((zinv-zlmeldwn)/zdz)ztv_est1
[2046]	300
[2106]	301	zbuoyjam(ig,l)=fact_shellRG(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- &
	302	& ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
	303	& ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
[2046]	304
[2106]	305	else
	306	ztv_est(ig,l)=atv1*zlmel+btv1
	307	zbuoyjam(ig,l)=fact_shellRG(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
	308	& +(1.-fact_shell)*zbuoy(ig,l)
	309	endif
[2046]	310
[2106]	311	else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
[2046]	312
[2106]	313	if (zlmeldwn.gt.zltdwn) then
	314	zbuoyjam(ig,l)=fact_shellRG((ztva_est(ig,l)- &
	315	& ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l)
	316	else
	317	zbuoyjam(ig,l)=fact_shellRG(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
	318	& ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
	319	& ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
[2046]	320
[2106]	321	endif
[2140]	322
[2106]	323	! zbuoyjam(ig,l)=fact_shellRG(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
	324	! & ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
	325	! & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
[1968]	326	! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
	327	! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
[2046]	328	! & po(ig,lt-1))/po(ig,lt-1))
[2106]	329	endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
[2046]	330
[2106]	331	else ! if (iflag_thermals_ed.lt.8) then
	332	lt=l+1
[2267]	333	zlt=zlev(ig,lt)
[2106]	334	zdz2=zlev(ig,lt)-zlev(ig,l)
[2046]	335
[2106]	336	do while (lmel.gt.zdz2)
	337	lt=lt+1
	338	zlt=zlev(ig,lt)
	339	zdz2=zlt-zlev(ig,l)
	340	enddo
	341	zdz3=zlev(ig,lt+1)-zlt
	342	zltdwn=zlev(ig,lt)-zdz3/2
[2267]	343	zlmelup=zlmel+(zdz/2)
	344	coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz)
	345	zbuoyjam(ig,l)=1.RG(coefzlmel*(ztva_est(ig,l)- &
	346	& ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- &
[2046]	347	& ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
[2106]	348	endif ! if (iflag_thermals_ed.lt.8) then
[2046]	349
[2140]	350	!------------------------------------------------
	351	!AJAM:nouveau calcul de w?
	352	!------------------------------------------------
	353	zdz=zlev(ig,l+1)-zlev(ig,l)
	354	zdzbis=zlev(ig,l)-zlev(ig,l-1)
	355	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	356
	357	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	358	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
	359	zdw2=afact*zbuoy(ig,l)/fact_epsilon
	360	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
[2267]	361	! zdw2bis=0.5*(zdw2+zdw2bis)
[2140]	362	lm=Max(1,l-2)
	363	! zdw2=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l) &
	364	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
	365	! zdw2bis=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l-1) &
	366	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
	367	! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
	368	! w_est(ig,l+1)=(zdz/zdzbis)Max(0.0001,exp(-zw2fact) &
	369	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
	370	! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
	371	! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))zdw2+w_est(ig,l)exp(-zw2fact))
	372
	373	!--------------------------------------------------
	374	!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
	375	!--------------------------------------------------
	376	if (iflag_thermals_ed==8) then
	377	! Ancienne version
	378	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
	379	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
	380	! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
	381
[4094]	382	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
[2140]	383
	384	! Nouvelle version Arnaud
	385	else
	386	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
	387	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
	388	! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
	389
[4094]	390	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
[2140]	391
	392	! w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))(exp(-zw2fact) &
	393	! & (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
	394	! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))
	395
	396
	397
	398	! w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2biszw2fact)exp(-zw2fact))
	399
	400	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(zw2(ig,l)+zdw2zw2fact)*exp(-zw2fact)+ &
	401	! & (zdzbis-zdz)/zdzbis(zw2(ig,l-1)+zdw2biszw2factbis)*exp(-zw2factbis))
	402
	403	! w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
	404
	405	endif
	406
	407
	408	if (iflag_thermals_ed<6) then
	409	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
	410	! fact_epsilon=0.0005/(zalpha+0.025)**0.5
	411	! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
	412	fact_epsilon=0.0002/(zalpha+0.1)
	413	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	414	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
	415	zdw2=afact*zbuoy(ig,l)/fact_epsilon
	416	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
	417	! w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2zw2fact)exp(-zw2fact))
	418
	419	! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
	420	! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
	421	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
	422
[4094]	423	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
[2140]	424
	425
	426	endif
	427	!--------------------------------------------------
	428	!AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
	429	!on fait max(0.0001,.....)
	430	!--------------------------------------------------
	431
	432	! if (w_est(ig,l+1).lt.0.) then
	433	! w_est(ig,l+1)=zw2(ig,l)
	434	! w_est(ig,l+1)=0.0001
	435	! endif
	436
	437	endif
	438	enddo
	439
	440
	441	!-------------------------------------------------
	442	!calcul des taux d'entrainement et de detrainement
	443	!-------------------------------------------------
	444
	445	do ig=1,ngrid
	446	if (active(ig)) then
	447
	448	! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
	449	zw2m=w_est(ig,l+1)
	450	! zw2m=zw2(ig,l)
	451	zdz=zlev(ig,l+1)-zlev(ig,l)
	452	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	453	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
	454	zbuoybis=zbuoy(ig,l)
	455	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
	456	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
	457
	458
	459	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
	460	! & afact*zbuoybis/zw2m - fact_epsilon )
	461
	462	! entr_star(ig,l)=MAX(0.,f_star(ig,l)zdzzbetalpha* &
	463	! & afact*zbuoybis/zw2m - fact_epsilon )
	464
	465
	466
[1968]	467	! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	468
[2106]	469	!=========================================================================
	470	! 4. Calcul de l'entrainement et du detrainement
	471	!=========================================================================
	472
[1998]	473	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
	474	! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
	475	! entrbis=entr_star(ig,l)
[1968]	476
[2106]	477	if (iflag_thermals_ed.lt.6) then
[2046]	478	fact_epsilon=0.0002/(zalpha+0.1)
[2106]	479	endif
	480
[1968]	481
[2106]	482
[2046]	483	detr_star(ig,l)=f_star(ig,l)*zdz &
[2106]	484	& ( mix0 0.1 / (zalpha+0.001) &
	485	& + MAX(detr_min, -afactzbetalphazbuoyjam(ig,l)/zw2m &
	486	& + detr_q_coef(zdqt(ig,l)/zw2m)*detr_q_power))
[1968]	487
[4094]	488	! detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
	489	! & ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)
[2140]	490
[1998]	491	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
[1968]	492
[2106]	493	entr_star(ig,l)=f_star(ig,l)zdz ( &
	494	& mix0 * 0.1 / (zalpha+0.001) &
	495	& + zbetalpha*MAX(entr_min, &
[2267]	496	& afact*zbuoyjam(ig,l)/zw2m - fact_epsilon))
[1968]	497
[2140]	498
	499	! entr_star(ig,l)=f_star(ig,l)zdz ( &
	500	! & mix0 * 0.1 / (zalpha+0.001) &
	501	! & + MAX(entr_min, &
	502	! & zbetalphaafactzbuoyjam(ig,l)/zw2m - fact_epsilon + &
	503	! & detr_q_coef(zdqt(ig,l)/zw2m)*detr_q_power))
	504
	505
[4094]	506	! entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
	507	! & ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)
[2140]	508
[1998]	509	! entr_star(ig,l)=Max(0.,f_star(ig,l)zdzzbetalpha* &
	510	! & afact*zbuoy(ig,l)/zw2m &
	511	! & - 1.*fact_epsilon)
	512
[1968]	513
[1403]	514	! En dessous de lalim, on prend le max de alim_star et entr_star pour
	515	! alim_star et 0 sinon
	516	if (l.lt.lalim(ig)) then
	517	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
	518	entr_star(ig,l)=0.
	519	endif
[1968]	520	! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
	521	! alim_star(ig,l)=entrbis
	522	! endif
[972]	523
[2106]	524	! print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
[1403]	525	! Calcul du flux montant normalise
	526	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	527	& -detr_star(ig,l)
[972]	528
[1403]	529	endif
	530	enddo
[972]	531
[1968]	532
[2106]	533	!============================================================================
	534	! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
	535	!===========================================================================
	536
[1403]	537	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
	538	do ig=1,ngrid
	539	if (activetmp(ig)) then
	540	Zsat=.false.
	541	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
	542	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	543	& /(f_star(ig,l+1)+detr_star(ig,l))
	544	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
	545	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	546	& /(f_star(ig,l+1)+detr_star(ig,l))
[972]	547
[1403]	548	endif
	549	enddo
[972]	550
[1968]	551	ztemp(:)=zpspsk(:,l)*ztla(:,l)
	552	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
[1403]	553	do ig=1,ngrid
	554	if (activetmp(ig)) then
	555	! on ecrit de maniere conservative (sat ou non)
	556	! T = Tl +Lv/Cp ql
[1968]	557	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
[1403]	558	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	559	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	560	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	561	zha(ig,l) = ztva(ig,l)
	562	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	563	& -zqla(ig,l))-zqla(ig,l))
[1968]	564	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
[1403]	565	zdz=zlev(ig,l+1)-zlev(ig,l)
[2140]	566	zdzbis=zlev(ig,l)-zlev(ig,l-1)
[1968]	567	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
[2106]	568	!!!!!!! fact_epsilon=0.002
[1968]	569	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	570	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
	571	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
	572	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
[2140]	573	! zdw2=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l) &
	574	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
	575	! lm=Max(1,l-2)
	576	! zdw2bis=(afact/fact_epsilon)((zdz/zdzbis)zbuoy(ig,l-1) &
	577	! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
[2106]	578	! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
	579	! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(zw2(ig,l)+zdw2zw2fact)*exp(-zw2fact)+ &
	580	! & (zdzbis-zdz)/zdzbis(zw2(ig,l-1)+zdw2biszw2factbis)*exp(-zw2factbis))
[2140]	581	! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2zw2fact)exp(-zw2fact))
[2106]	582	! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
	583	! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
[2140]	584	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
[2267]	585	if (iflag_thermals_ed==8) then
[4094]	586	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
[2267]	587	else
[4094]	588	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
[2267]	589	endif
[2140]	590	! zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))(exp(-zw2fact) &
[2267]	591	! & (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
[2106]	592	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))
[2046]	593
[2140]	594
[2106]	595	if (iflag_thermals_ed.lt.6) then
[2046]	596	zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
	597	! fact_epsilon=0.0005/(zalpha+0.025)**0.5
	598	! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
	599	fact_epsilon=0.0002/(zalpha+0.1)**1
	600	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	601	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
	602	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
	603	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
	604
[2106]	605	! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
	606	! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
	607	! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
[2140]	608	! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2zw2fact)exp(-zw2fact))
[4094]	609	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
[2046]	610
[2106]	611	endif
[2046]	612
	613
[1403]	614	endif
	615	enddo
	616
	617	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
[878]	618	!
[2106]	619	!===========================================================================
	620	! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
	621	!===========================================================================
[972]	622
[1503]	623	nbpb=0
[1403]	624	do ig=1,ngrid
	625	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	626	! stop'On tombe sur le cas particulier de thermcell_dry'
[1503]	627	! print*,'On tombe sur le cas particulier de thermcell_plume'
	628	nbpb=nbpb+1
[1403]	629	zw2(ig,l+1)=0.
	630	linter(ig)=l+1
	631	endif
[972]	632
[1403]	633	if (zw2(ig,l+1).lt.0.) then
	634	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	635	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	636	zw2(ig,l+1)=0.
[1998]	637	!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
	638	elseif (f_star(ig,l+1).lt.0.) then
	639	linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
	640	& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
	641	zw2(ig,l+1)=0.
	642	!fin CR:04/05/12
[1403]	643	endif
[972]	644
[1403]	645	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
[972]	646
[1403]	647	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	648	! lmix est le niveau de la couche ou w (wa_moy) est maximum
	649	!on rajoute le calcul de lmix_bis
	650	if (zqla(ig,l).lt.1.e-10) then
	651	lmix_bis(ig)=l+1
	652	endif
	653	lmix(ig)=l+1
	654	wmaxa(ig)=wa_moy(ig,l+1)
	655	endif
	656	enddo
[972]	657
[1503]	658	if (nbpb>0) then
	659	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
	660	endif
	661
[1403]	662	!=========================================================================
	663	! FIN DE LA BOUCLE VERTICALE
	664	enddo
	665	!=========================================================================
[972]	666
[1403]	667	!on recalcule alim_star_tot
	668	do ig=1,ngrid
	669	alim_star_tot(ig)=0.
	670	enddo
	671	do ig=1,ngrid
	672	do l=1,lalim(ig)-1
	673	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	674	enddo
	675	enddo
	676
[972]	677
[1403]	678	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
[972]	679
[2046]	680	#undef wrgrads_thermcell
	681	#ifdef wrgrads_thermcell
[4094]	682	call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta ','esta ')
	683	call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta ','dsta ')
	684	call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy ','buoy ')
	685	call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt ','dqt ')
	686	call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est ','w_est ')
	687	call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2 ','w_es2 ')
	688	call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A ','zw2A ')
[2046]	689	#endif
	690
	691
[4094]	692	RETURN
[1968]	693	end
	694
	695
[1998]	696
	697
[2046]	698
	699
	700
[1403]	701	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	702	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	703	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	704	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	705	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	706	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[4094]	707	SUBROUTINE thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
[1403]	708	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	709	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	710	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
[2149]	711	& ,lev_out,lunout1,igout)
	712	!& ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
[972]	713
[1403]	714	!--------------------------------------------------------------------------
	715	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
	716	! Version conforme a l'article de Rio et al. 2010.
	717	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
	718	!--------------------------------------------------------------------------
[878]	719
[4089]	720	USE thermcell_ini_mod, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
[1403]	721	IMPLICIT NONE
[972]	722
[1403]	723	INTEGER itap
	724	INTEGER lunout1,igout
[4094]	725	INTEGER ngrid,nlay
[1403]	726	REAL ptimestep
[4094]	727	REAL ztv(ngrid,nlay)
	728	REAL zthl(ngrid,nlay)
	729	REAL po(ngrid,nlay)
	730	REAL zl(ngrid,nlay)
	731	REAL rhobarz(ngrid,nlay)
	732	REAL zlev(ngrid,nlay+1)
	733	REAL pplev(ngrid,nlay+1)
	734	REAL pphi(ngrid,nlay)
	735	REAL zpspsk(ngrid,nlay)
	736	REAL alim_star(ngrid,nlay)
[1403]	737	REAL f0(ngrid)
	738	INTEGER lalim(ngrid)
	739	integer lev_out ! niveau pour les print
[1503]	740	integer nbpb
[1403]	741
	742	real alim_star_tot(ngrid)
[878]	743
[4094]	744	REAL ztva(ngrid,nlay)
	745	REAL ztla(ngrid,nlay)
	746	REAL zqla(ngrid,nlay)
	747	REAL zqta(ngrid,nlay)
	748	REAL zha(ngrid,nlay)
[878]	749
[4094]	750	REAL detr_star(ngrid,nlay)
[1403]	751	REAL coefc
[4094]	752	REAL entr_star(ngrid,nlay)
	753	REAL detr(ngrid,nlay)
	754	REAL entr(ngrid,nlay)
[1403]	755
[4094]	756	REAL csc(ngrid,nlay)
[1403]	757
[4094]	758	REAL zw2(ngrid,nlay+1)
	759	REAL w_est(ngrid,nlay+1)
	760	REAL f_star(ngrid,nlay+1)
	761	REAL wa_moy(ngrid,nlay+1)
[1403]	762
[4094]	763	REAL ztva_est(ngrid,nlay)
	764	REAL zqla_est(ngrid,nlay)
	765	REAL zqsatth(ngrid,nlay)
	766	REAL zta_est(ngrid,nlay)
	767	REAL zbuoyjam(ngrid,nlay)
[1403]	768	REAL ztemp(ngrid),zqsat(ngrid)
	769	REAL zdw2
	770	REAL zw2modif
	771	REAL zw2fact
[4094]	772	REAL zeps(ngrid,nlay)
[1403]	773
	774	REAL linter(ngrid)
	775	INTEGER lmix(ngrid)
	776	INTEGER lmix_bis(ngrid)
	777	REAL wmaxa(ngrid)
	778
	779	INTEGER ig,l,k
	780
[4094]	781	real zdz,zbuoy(ngrid,nlay),zalpha,gamma(ngrid,nlay),zdqt(ngrid,nlay),zw2m
[1403]	782	real zbuoybis
	783	real zcor,zdelta,zcvm5,qlbef,zdz2
	784	real betalpha,zbetalpha
	785	real eps, afact
	786	logical Zsat
	787	LOGICAL active(ngrid),activetmp(ngrid)
	788	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
[4094]	789	REAL c2(ngrid,nlay)
[1403]	790	Zsat=.false.
	791	! Initialisation
	792
	793	fact_epsilon=0.002
	794	betalpha=0.9
	795	afact=2./3.
	796
	797	zbetalpha=betalpha/(1.+betalpha)
	798
	799
	800	! Initialisations des variables reeles
[1998]	801	if (1==1) then
[1403]	802	ztva(:,:)=ztv(:,:)
	803	ztva_est(:,:)=ztva(:,:)
	804	ztla(:,:)=zthl(:,:)
	805	zqta(:,:)=po(:,:)
	806	zha(:,:) = ztva(:,:)
	807	else
	808	ztva(:,:)=0.
	809	ztva_est(:,:)=0.
	810	ztla(:,:)=0.
	811	zqta(:,:)=0.
	812	zha(:,:) =0.
	813	endif
	814
	815	zqla_est(:,:)=0.
	816	zqsatth(:,:)=0.
	817	zqla(:,:)=0.
	818	detr_star(:,:)=0.
	819	entr_star(:,:)=0.
	820	alim_star(:,:)=0.
	821	alim_star_tot(:)=0.
	822	csc(:,:)=0.
	823	detr(:,:)=0.
	824	entr(:,:)=0.
	825	zw2(:,:)=0.
	826	zbuoy(:,:)=0.
[1998]	827	zbuoyjam(:,:)=0.
[1403]	828	gamma(:,:)=0.
	829	zeps(:,:)=0.
	830	w_est(:,:)=0.
	831	f_star(:,:)=0.
	832	wa_moy(:,:)=0.
	833	linter(:)=1.
	834	! linter(:)=1.
	835	! Initialisation des variables entieres
	836	lmix(:)=1
	837	lmix_bis(:)=2
	838	wmaxa(:)=0.
	839	lalim(:)=1
	840
	841
	842	!-------------------------------------------------------------------------
	843	! On ne considere comme actif que les colonnes dont les deux premieres
	844	! couches sont instables.
	845	!-------------------------------------------------------------------------
	846	active(:)=ztv(:,1)>ztv(:,2)
	847
	848	!-------------------------------------------------------------------------
[4089]	849	! Definition de l'alimentation
[1403]	850	!-------------------------------------------------------------------------
[4094]	851	do l=1,nlay-1
[1403]	852	do ig=1,ngrid
	853	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
	854	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
	855	& *sqrt(zlev(ig,l+1))
	856	lalim(ig)=l+1
	857	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
[878]	858	endif
[1403]	859	enddo
	860	enddo
[4094]	861	do l=1,nlay
[1403]	862	do ig=1,ngrid
	863	if (alim_star_tot(ig) > 1.e-10 ) then
	864	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
	865	endif
	866	enddo
	867	enddo
	868	alim_star_tot(:)=1.
[878]	869
	870
[972]	871
[1403]	872	!------------------------------------------------------------------------------
	873	! Calcul dans la premiere couche
	874	! On decide dans cette version que le thermique n'est actif que si la premiere
	875	! couche est instable.
[1968]	876	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
[1403]	877	! dans une couche l>1
	878	!------------------------------------------------------------------------------
	879	do ig=1,ngrid
	880	! Le panache va prendre au debut les caracteristiques de l'air contenu
	881	! dans cette couche.
	882	if (active(ig)) then
	883	ztla(ig,1)=zthl(ig,1)
	884	zqta(ig,1)=po(ig,1)
	885	zqla(ig,1)=zl(ig,1)
	886	!cr: attention, prise en compte de f*(1)=1
	887	f_star(ig,2)=alim_star(ig,1)
	888	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
	889	& *(zlev(ig,2)-zlev(ig,1)) &
	890	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
	891	w_est(ig,2)=zw2(ig,2)
	892	endif
	893	enddo
	894	!
[1026]	895
[1403]	896	!==============================================================================
	897	!boucle de calcul de la vitesse verticale dans le thermique
	898	!==============================================================================
[4094]	899	do l=2,nlay-1
[1403]	900	!==============================================================================
[1026]	901
	902
[1403]	903	! On decide si le thermique est encore actif ou non
	904	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
	905	do ig=1,ngrid
	906	active(ig)=active(ig) &
	907	& .and. zw2(ig,l)>1.e-10 &
	908	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
	909	enddo
[1026]	910
	911
	912
[1403]	913	!---------------------------------------------------------------------------
	914	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
	915	! sans tenir compte du detrainement et de l'entrainement dans cette
	916	! couche
	917	! C'est a dire qu'on suppose
	918	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
	919	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
	920	! avant) a l'alimentation pour avoir un calcul plus propre
	921	!---------------------------------------------------------------------------
	922
	923	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
	924	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
	925
	926	do ig=1,ngrid
	927	! print*,'active',active(ig),ig,l
	928	if(active(ig)) then
	929	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
	930	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
	931	zta_est(ig,l)=ztva_est(ig,l)
	932	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
	933	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
	934	& -zqla_est(ig,l))-zqla_est(ig,l))
	935
	936	!------------------------------------------------
[1968]	937	!AJAM:nouveau calcul de w?
[1403]	938	!------------------------------------------------
	939	zdz=zlev(ig,l+1)-zlev(ig,l)
	940	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	941
	942	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	943	zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
	944	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
	945
	946
	947	if (w_est(ig,l+1).lt.0.) then
	948	w_est(ig,l+1)=zw2(ig,l)
	949	endif
[1026]	950	endif
[1403]	951	enddo
[1026]	952
	953
[1403]	954	!-------------------------------------------------
	955	!calcul des taux d'entrainement et de detrainement
	956	!-------------------------------------------------
[1026]	957
[1403]	958	do ig=1,ngrid
	959	if (active(ig)) then
[1026]	960
[1403]	961	zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
	962	zw2m=w_est(ig,l+1)
	963	zdz=zlev(ig,l+1)-zlev(ig,l)
	964	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	965	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
	966	zbuoybis=zbuoy(ig,l)
	967	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
	968	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
[1026]	969
[1403]	970
	971	entr_star(ig,l)=f_star(ig,l)zdz zbetalpha*MAX(0., &
	972	& afact*zbuoybis/zw2m - fact_epsilon )
	973
	974
	975	detr_star(ig,l)=f_star(ig,l)*zdz &
	976	& MAX(1.e-3, -afactzbetalpha*zbuoy(ig,l)/zw2m &
	977	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
	978
	979	! En dessous de lalim, on prend le max de alim_star et entr_star pour
	980	! alim_star et 0 sinon
	981	if (l.lt.lalim(ig)) then
	982	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
	983	entr_star(ig,l)=0.
	984	endif
	985
	986	! Calcul du flux montant normalise
[878]	987	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	988	& -detr_star(ig,l)
	989
[1403]	990	endif
	991	enddo
	992
	993
[878]	994	!----------------------------------------------------------------------------
	995	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	996	!---------------------------------------------------------------------------
[1403]	997	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
	998	do ig=1,ngrid
	999	if (activetmp(ig)) then
	1000	Zsat=.false.
	1001	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
[878]	1002	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	1003	& /(f_star(ig,l+1)+detr_star(ig,l))
[1403]	1004	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
[878]	1005	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	1006	& /(f_star(ig,l+1)+detr_star(ig,l))
	1007
[1403]	1008	endif
	1009	enddo
	1010
	1011	ztemp(:)=zpspsk(:,l)*ztla(:,l)
	1012	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
	1013
	1014	do ig=1,ngrid
	1015	if (activetmp(ig)) then
[878]	1016	! on ecrit de maniere conservative (sat ou non)
	1017	! T = Tl +Lv/Cp ql
[1403]	1018	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
[878]	1019	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	1020	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	1021	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	1022	zha(ig,l) = ztva(ig,l)
	1023	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	1024	& -zqla(ig,l))-zqla(ig,l))
[1403]	1025	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	1026	zdz=zlev(ig,l+1)-zlev(ig,l)
	1027	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
[878]	1028
[1403]	1029	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	1030	zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
	1031	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
	1032	endif
	1033	enddo
[1026]	1034
[972]	1035	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
[878]	1036	!
[1403]	1037	!---------------------------------------------------------------------------
[878]	1038	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
[1403]	1039	!---------------------------------------------------------------------------
[878]	1040
[1503]	1041	nbpb=0
[1403]	1042	do ig=1,ngrid
[878]	1043	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	1044	! stop'On tombe sur le cas particulier de thermcell_dry'
[1503]	1045	! print*,'On tombe sur le cas particulier de thermcell_plume'
	1046	nbpb=nbpb+1
[878]	1047	zw2(ig,l+1)=0.
	1048	linter(ig)=l+1
[2046]	1049	endif
[878]	1050
	1051	if (zw2(ig,l+1).lt.0.) then
	1052	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	1053	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	1054	zw2(ig,l+1)=0.
[2106]	1055	elseif (f_star(ig,l+1).lt.0.) then
	1056	linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
	1057	& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
[2159]	1058	! print*,"linter plume", linter(ig)
[2106]	1059	zw2(ig,l+1)=0.
[878]	1060	endif
	1061
	1062	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	1063
	1064	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	1065	! lmix est le niveau de la couche ou w (wa_moy) est maximum
[1026]	1066	!on rajoute le calcul de lmix_bis
	1067	if (zqla(ig,l).lt.1.e-10) then
	1068	lmix_bis(ig)=l+1
	1069	endif
[878]	1070	lmix(ig)=l+1
	1071	wmaxa(ig)=wa_moy(ig,l+1)
	1072	endif
[1403]	1073	enddo
	1074
[1503]	1075	if (nbpb>0) then
	1076	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
	1077	endif
	1078
[1403]	1079	!=========================================================================
	1080	! FIN DE LA BOUCLE VERTICALE
[878]	1081	enddo
[1403]	1082	!=========================================================================
[878]	1083
[1403]	1084	!on recalcule alim_star_tot
	1085	do ig=1,ngrid
	1086	alim_star_tot(ig)=0.
	1087	enddo
	1088	do ig=1,ngrid
	1089	do l=1,lalim(ig)-1
	1090	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	1091	enddo
	1092	enddo
	1093
	1094
[972]	1095	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
[878]	1096
[2046]	1097	#undef wrgrads_thermcell
	1098	#ifdef wrgrads_thermcell
[4094]	1099	call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta ','esta ')
	1100	call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta ','dsta ')
	1101	call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy ','buoy ')
	1102	call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt ','dqt ')
	1103	call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est ','w_est ')
	1104	call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2 ','w_es2 ')
	1105	call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A ','zw2A ')
[2046]	1106	#endif
[1998]	1107
	1108
[1403]	1109	return
	1110	end

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

Download in other formats:

Original Format