Context Navigation

thermcell_plume_6A.F90 @ 4163

Last change on this file since 4163 was 4161, checked in by fhourdin, 3 years ago
For replay_param.sh
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.3 KB

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

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

Download in other formats:

Original Format