Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

lmdz_thermcell_plume_6A.f90 @ 5445

Last change on this file since 5445 was 5434, checked in by fhourdin, 3 days ago

Superessing CPP in lmdz_*

Not possible for lmdz_thermcell_main because of isotopes

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: 39.6 KB

Rev	Line
[4590]	1	MODULE lmdz_thermcell_plume_6A
[1403]	2	!
	3	! $Id: lmdz_thermcell_plume_6A.f90 5434 2024-12-20 10:48:05Z fhourdin $
	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
[4094]	686	RETURN
[5390]	687	END SUBROUTINE thermcell_plume_6A
[1968]	688
	689
[1998]	690
	691
[2046]	692
	693
	694
[1403]	695	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	696	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	697	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	698	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	699	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	700	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[4094]	701	SUBROUTINE thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
[1403]	702	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	703	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	704	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
[2149]	705	& ,lev_out,lunout1,igout)
	706	!& ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
[972]	707
[1403]	708	!--------------------------------------------------------------------------
	709	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
	710	! Version conforme a l'article de Rio et al. 2010.
	711	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
	712	!--------------------------------------------------------------------------
[878]	713
[4590]	714	USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG
	715	USE lmdz_thermcell_qsat, ONLY : thermcell_qsat
[1403]	716	IMPLICIT NONE
[972]	717
[1403]	718	INTEGER itap
	719	INTEGER lunout1,igout
[4094]	720	INTEGER ngrid,nlay
[1403]	721	REAL ptimestep
[4094]	722	REAL ztv(ngrid,nlay)
	723	REAL zthl(ngrid,nlay)
[4678]	724	REAL, INTENT(IN) :: po(ngrid,nlay)
[4094]	725	REAL zl(ngrid,nlay)
	726	REAL rhobarz(ngrid,nlay)
	727	REAL zlev(ngrid,nlay+1)
	728	REAL pplev(ngrid,nlay+1)
	729	REAL pphi(ngrid,nlay)
	730	REAL zpspsk(ngrid,nlay)
	731	REAL alim_star(ngrid,nlay)
[1403]	732	REAL f0(ngrid)
	733	INTEGER lalim(ngrid)
	734	integer lev_out ! niveau pour les print
[1503]	735	integer nbpb
[1403]	736
	737	real alim_star_tot(ngrid)
[878]	738
[4094]	739	REAL ztva(ngrid,nlay)
	740	REAL ztla(ngrid,nlay)
	741	REAL zqla(ngrid,nlay)
	742	REAL zqta(ngrid,nlay)
	743	REAL zha(ngrid,nlay)
[878]	744
[4094]	745	REAL detr_star(ngrid,nlay)
[1403]	746	REAL coefc
[4094]	747	REAL entr_star(ngrid,nlay)
	748	REAL detr(ngrid,nlay)
	749	REAL entr(ngrid,nlay)
[1403]	750
[4094]	751	REAL csc(ngrid,nlay)
[1403]	752
[4094]	753	REAL zw2(ngrid,nlay+1)
	754	REAL w_est(ngrid,nlay+1)
	755	REAL f_star(ngrid,nlay+1)
	756	REAL wa_moy(ngrid,nlay+1)
[1403]	757
[4094]	758	REAL ztva_est(ngrid,nlay)
	759	REAL zqla_est(ngrid,nlay)
	760	REAL zqsatth(ngrid,nlay)
	761	REAL zta_est(ngrid,nlay)
	762	REAL zbuoyjam(ngrid,nlay)
[1403]	763	REAL ztemp(ngrid),zqsat(ngrid)
	764	REAL zdw2
	765	REAL zw2modif
	766	REAL zw2fact
[4094]	767	REAL zeps(ngrid,nlay)
[1403]	768
	769	REAL linter(ngrid)
	770	INTEGER lmix(ngrid)
	771	INTEGER lmix_bis(ngrid)
	772	REAL wmaxa(ngrid)
	773
	774	INTEGER ig,l,k
	775
[4094]	776	real zdz,zbuoy(ngrid,nlay),zalpha,gamma(ngrid,nlay),zdqt(ngrid,nlay),zw2m
[1403]	777	real zbuoybis
	778	real zcor,zdelta,zcvm5,qlbef,zdz2
	779	real betalpha,zbetalpha
	780	real eps, afact
	781	logical Zsat
	782	LOGICAL active(ngrid),activetmp(ngrid)
	783	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
[4094]	784	REAL c2(ngrid,nlay)
[1403]	785	Zsat=.false.
	786	! Initialisation
	787
	788	fact_epsilon=0.002
	789	betalpha=0.9
	790	afact=2./3.
	791
	792	zbetalpha=betalpha/(1.+betalpha)
	793
	794
	795	! Initialisations des variables reeles
[1998]	796	if (1==1) then
[1403]	797	ztva(:,:)=ztv(:,:)
	798	ztva_est(:,:)=ztva(:,:)
	799	ztla(:,:)=zthl(:,:)
	800	zqta(:,:)=po(:,:)
	801	zha(:,:) = ztva(:,:)
	802	else
	803	ztva(:,:)=0.
	804	ztva_est(:,:)=0.
	805	ztla(:,:)=0.
	806	zqta(:,:)=0.
	807	zha(:,:) =0.
	808	endif
	809
	810	zqla_est(:,:)=0.
	811	zqsatth(:,:)=0.
	812	zqla(:,:)=0.
	813	detr_star(:,:)=0.
	814	entr_star(:,:)=0.
	815	alim_star(:,:)=0.
	816	alim_star_tot(:)=0.
	817	csc(:,:)=0.
	818	detr(:,:)=0.
	819	entr(:,:)=0.
	820	zw2(:,:)=0.
	821	zbuoy(:,:)=0.
[1998]	822	zbuoyjam(:,:)=0.
[1403]	823	gamma(:,:)=0.
	824	zeps(:,:)=0.
	825	w_est(:,:)=0.
	826	f_star(:,:)=0.
	827	wa_moy(:,:)=0.
	828	linter(:)=1.
	829	! linter(:)=1.
	830	! Initialisation des variables entieres
	831	lmix(:)=1
	832	lmix_bis(:)=2
	833	wmaxa(:)=0.
	834	lalim(:)=1
	835
	836
	837	!-------------------------------------------------------------------------
	838	! On ne considere comme actif que les colonnes dont les deux premieres
	839	! couches sont instables.
	840	!-------------------------------------------------------------------------
	841	active(:)=ztv(:,1)>ztv(:,2)
	842
	843	!-------------------------------------------------------------------------
[4089]	844	! Definition de l'alimentation
[1403]	845	!-------------------------------------------------------------------------
[4094]	846	do l=1,nlay-1
[1403]	847	do ig=1,ngrid
	848	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
	849	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
	850	& *sqrt(zlev(ig,l+1))
	851	lalim(ig)=l+1
	852	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
[878]	853	endif
[1403]	854	enddo
	855	enddo
[4094]	856	do l=1,nlay
[1403]	857	do ig=1,ngrid
	858	if (alim_star_tot(ig) > 1.e-10 ) then
	859	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
	860	endif
	861	enddo
	862	enddo
	863	alim_star_tot(:)=1.
[878]	864
	865
[972]	866
[1403]	867	!------------------------------------------------------------------------------
	868	! Calcul dans la premiere couche
	869	! On decide dans cette version que le thermique n'est actif que si la premiere
	870	! couche est instable.
[1968]	871	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
[1403]	872	! dans une couche l>1
	873	!------------------------------------------------------------------------------
	874	do ig=1,ngrid
	875	! Le panache va prendre au debut les caracteristiques de l'air contenu
	876	! dans cette couche.
	877	if (active(ig)) then
	878	ztla(ig,1)=zthl(ig,1)
	879	zqta(ig,1)=po(ig,1)
	880	zqla(ig,1)=zl(ig,1)
	881	!cr: attention, prise en compte de f*(1)=1
	882	f_star(ig,2)=alim_star(ig,1)
	883	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
	884	& *(zlev(ig,2)-zlev(ig,1)) &
	885	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
	886	w_est(ig,2)=zw2(ig,2)
	887	endif
	888	enddo
	889	!
[1026]	890
[1403]	891	!==============================================================================
	892	!boucle de calcul de la vitesse verticale dans le thermique
	893	!==============================================================================
[4094]	894	do l=2,nlay-1
[1403]	895	!==============================================================================
[1026]	896
	897
[1403]	898	! On decide si le thermique est encore actif ou non
	899	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
	900	do ig=1,ngrid
	901	active(ig)=active(ig) &
	902	& .and. zw2(ig,l)>1.e-10 &
	903	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
	904	enddo
[1026]	905
	906
	907
[1403]	908	!---------------------------------------------------------------------------
	909	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
	910	! sans tenir compte du detrainement et de l'entrainement dans cette
	911	! couche
	912	! C'est a dire qu'on suppose
	913	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
	914	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
	915	! avant) a l'alimentation pour avoir un calcul plus propre
	916	!---------------------------------------------------------------------------
	917
	918	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
	919	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
	920
	921	do ig=1,ngrid
	922	! print*,'active',active(ig),ig,l
	923	if(active(ig)) then
	924	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
	925	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
	926	zta_est(ig,l)=ztva_est(ig,l)
	927	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
	928	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
	929	& -zqla_est(ig,l))-zqla_est(ig,l))
	930
	931	!------------------------------------------------
[1968]	932	!AJAM:nouveau calcul de w?
[1403]	933	!------------------------------------------------
	934	zdz=zlev(ig,l+1)-zlev(ig,l)
	935	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	936
	937	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	938	zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
	939	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
	940
	941
	942	if (w_est(ig,l+1).lt.0.) then
	943	w_est(ig,l+1)=zw2(ig,l)
	944	endif
[1026]	945	endif
[1403]	946	enddo
[1026]	947
	948
[1403]	949	!-------------------------------------------------
	950	!calcul des taux d'entrainement et de detrainement
	951	!-------------------------------------------------
[1026]	952
[1403]	953	do ig=1,ngrid
	954	if (active(ig)) then
[1026]	955
[1403]	956	zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
	957	zw2m=w_est(ig,l+1)
	958	zdz=zlev(ig,l+1)-zlev(ig,l)
	959	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	960	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
	961	zbuoybis=zbuoy(ig,l)
	962	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
	963	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
[1026]	964
[1403]	965
	966	entr_star(ig,l)=f_star(ig,l)zdz zbetalpha*MAX(0., &
	967	& afact*zbuoybis/zw2m - fact_epsilon )
	968
	969
	970	detr_star(ig,l)=f_star(ig,l)*zdz &
	971	& MAX(1.e-3, -afactzbetalpha*zbuoy(ig,l)/zw2m &
	972	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
	973
	974	! En dessous de lalim, on prend le max de alim_star et entr_star pour
	975	! alim_star et 0 sinon
	976	if (l.lt.lalim(ig)) then
	977	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
	978	entr_star(ig,l)=0.
	979	endif
	980
	981	! Calcul du flux montant normalise
[878]	982	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	983	& -detr_star(ig,l)
	984
[1403]	985	endif
	986	enddo
	987
	988
[878]	989	!----------------------------------------------------------------------------
	990	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	991	!---------------------------------------------------------------------------
[1403]	992	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
	993	do ig=1,ngrid
	994	if (activetmp(ig)) then
	995	Zsat=.false.
	996	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
[878]	997	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	998	& /(f_star(ig,l+1)+detr_star(ig,l))
[1403]	999	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
[878]	1000	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	1001	& /(f_star(ig,l+1)+detr_star(ig,l))
	1002
[1403]	1003	endif
	1004	enddo
	1005
	1006	ztemp(:)=zpspsk(:,l)*ztla(:,l)
	1007	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
	1008
	1009	do ig=1,ngrid
	1010	if (activetmp(ig)) then
[878]	1011	! on ecrit de maniere conservative (sat ou non)
	1012	! T = Tl +Lv/Cp ql
[1403]	1013	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
[878]	1014	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	1015	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	1016	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	1017	zha(ig,l) = ztva(ig,l)
	1018	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	1019	& -zqla(ig,l))-zqla(ig,l))
[1403]	1020	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	1021	zdz=zlev(ig,l+1)-zlev(ig,l)
	1022	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
[878]	1023
[1403]	1024	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	1025	zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
	1026	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
	1027	endif
	1028	enddo
[1026]	1029
[972]	1030	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
[878]	1031	!
[1403]	1032	!---------------------------------------------------------------------------
[878]	1033	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
[1403]	1034	!---------------------------------------------------------------------------
[878]	1035
[1503]	1036	nbpb=0
[1403]	1037	do ig=1,ngrid
[878]	1038	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	1039	! stop'On tombe sur le cas particulier de thermcell_dry'
[1503]	1040	! print*,'On tombe sur le cas particulier de thermcell_plume'
	1041	nbpb=nbpb+1
[878]	1042	zw2(ig,l+1)=0.
	1043	linter(ig)=l+1
[2046]	1044	endif
[878]	1045
	1046	if (zw2(ig,l+1).lt.0.) then
	1047	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	1048	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	1049	zw2(ig,l+1)=0.
[2106]	1050	elseif (f_star(ig,l+1).lt.0.) then
	1051	linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
	1052	& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
[2159]	1053	! print*,"linter plume", linter(ig)
[2106]	1054	zw2(ig,l+1)=0.
[878]	1055	endif
	1056
	1057	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	1058
	1059	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	1060	! lmix est le niveau de la couche ou w (wa_moy) est maximum
[1026]	1061	!on rajoute le calcul de lmix_bis
	1062	if (zqla(ig,l).lt.1.e-10) then
	1063	lmix_bis(ig)=l+1
	1064	endif
[878]	1065	lmix(ig)=l+1
	1066	wmaxa(ig)=wa_moy(ig,l+1)
	1067	endif
[1403]	1068	enddo
	1069
[1503]	1070	if (nbpb>0) then
	1071	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
	1072	endif
	1073
[1403]	1074	!=========================================================================
	1075	! FIN DE LA BOUCLE VERTICALE
[878]	1076	enddo
[1403]	1077	!=========================================================================
[878]	1078
[1403]	1079	!on recalcule alim_star_tot
	1080	do ig=1,ngrid
	1081	alim_star_tot(ig)=0.
	1082	enddo
	1083	do ig=1,ngrid
	1084	do l=1,lalim(ig)-1
	1085	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	1086	enddo
	1087	enddo
	1088
	1089
[972]	1090	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
[878]	1091
[1403]	1092	return
[5390]	1093	END SUBROUTINE thermcell_plume_5B
[4590]	1094	END MODULE lmdz_thermcell_plume_6A

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

Download in other formats:

Original Format