Context Navigation

lmdz_thermcell_plume_6A.F90 @ 5408

Last change on this file since 5408 was 5390, checked in by yann meurdesoif, 10 days ago

Remove UTF8 character that inihibit fortran parsing with GPU morphosis
Add missing END SUBROUTINE instead of simple END, that inhibit correct parsing with regulat expression parser (quick and dirty parsing)

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.8 KB

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

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

Download in other formats:

Original Format