Context Navigation

thermcell_plume.F90 @ 861

Last change on this file since 861 was 852, checked in by Laurent Fairhead, 17 years ago
Mise a jour de la physique avec thermiques avec la version de FH d'aout 2007 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 13.3 KB

Rev	Line
[814]	1	SUBROUTINE thermcell_plume(ngrid,klev,ztv,zthl,po,zl,rhobarz, &
	2	& zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star, &
[852]	3	& lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, &
[814]	4	& ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
	5
	6	!--------------------------------------------------------------------------
	7	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
	8	!--------------------------------------------------------------------------
	9
	10	IMPLICIT NONE
	11
	12	#include "YOMCST.h"
	13	#include "YOETHF.h"
	14	#include "FCTTRE.h"
	15
	16	INTEGER ngrid,klev
	17	REAL ztv(ngrid,klev)
	18	REAL zthl(ngrid,klev)
	19	REAL po(ngrid,klev)
	20	REAL zl(ngrid,klev)
	21	REAL rhobarz(ngrid,klev)
	22	REAL zlev(ngrid,klev+1)
	23	REAL pplev(ngrid,klev+1)
	24	REAL pphi(ngrid,klev)
	25	REAL zpspsk(ngrid,klev)
	26	REAL alim_star(ngrid,klev)
	27	REAL zmax_sec(ngrid)
	28	REAL f0(ngrid)
	29	REAL l_mix
	30	REAL r_aspect
[852]	31	INTEGER lalim(ngrid)
[814]	32	integer lev_out ! niveau pour les print
	33
	34	REAL ztva(ngrid,klev)
	35	REAL ztla(ngrid,klev)
	36	REAL zqla(ngrid,klev)
	37	REAL zqta(ngrid,klev)
	38	REAL zha(ngrid,klev)
	39
	40	REAL detr_star(ngrid,klev)
	41	REAL entr_star(ngrid,klev)
	42	REAL detr(ngrid,klev)
	43	REAL entr(ngrid,klev)
	44
	45	REAL zw2(ngrid,klev+1)
	46	REAL w_est(ngrid,klev+1)
	47	REAL f_star(ngrid,klev+1)
	48	REAL wa_moy(ngrid,klev+1)
	49
	50	REAL ztva_est(ngrid,klev)
	51	REAL zqla_est(ngrid,klev)
	52	REAL zqsatth(ngrid,klev)
	53
	54	REAL linter(ngrid)
	55	INTEGER lmix(ngrid)
	56	REAL wmaxa(ngrid)
	57
	58	INTEGER ig,l,k
	59
	60	real zcor,zdelta,zcvm5,qlbef
	61	real Tbef,qsatbef
	62	real dqsat_dT,DT,num,denom
	63	REAL REPS,RLvCp,DDT0
	64	PARAMETER (DDT0=.01)
	65	logical Zsat
	66
	67	Zsat=.false.
	68	! Initialisation
	69	RLvCp = RLVTT/RCPD
	70
	71	do l=1,klev
	72	do ig=1,ngrid
	73	zqla_est(ig,l)=0.
	74	ztva_est(ig,l)=ztva(ig,l)
	75	zqsatth(ig,l)=0.
	76	enddo
	77	enddo
	78
	79	!AM:initialisations du thermique
	80	do k=1,klev
	81	do ig=1,ngrid
	82	ztva(ig,k)=ztv(ig,k)
	83	ztla(ig,k)=zthl(ig,k)
	84	zqla(ig,k)=0.
	85	zqta(ig,k)=po(ig,k)
	86	enddo
	87	enddo
	88	do k=1,klev
	89	do ig=1,ngrid
	90	detr_star(ig,k)=0.
	91	entr_star(ig,k)=0.
	92	detr(ig,k)=0.
	93	entr(ig,k)=0.
	94	enddo
	95	enddo
	96	if (lev_out.ge.1) print*,'7 OK convect8'
	97	do k=1,klev+1
	98	do ig=1,ngrid
	99	zw2(ig,k)=0.
	100	w_est(ig,k)=0.
	101	f_star(ig,k)=0.
	102	wa_moy(ig,k)=0.
	103	enddo
	104	enddo
	105
	106	if (lev_out.ge.1) print*,'8 OK convect8'
	107	do ig=1,ngrid
	108	linter(ig)=1.
	109	lmix(ig)=1
	110	wmaxa(ig)=0.
	111	enddo
	112
	113	!-----------------------------------------------------------------------------------
	114	!boucle de calcul de la vitesse verticale dans le thermique
	115	!-----------------------------------------------------------------------------------
	116	do l=1,klev-1
	117	do ig=1,ngrid
	118	if (ztv(ig,l).gt.ztv(ig,l+1) &
	119	& .and.alim_star(ig,l).gt.1.e-10 &
	120	& .and.zw2(ig,l).lt.1e-10) then
	121	ztla(ig,l)=zthl(ig,l)
	122	zqta(ig,l)=po(ig,l)
	123	zqla(ig,l)=zl(ig,l)
	124	f_star(ig,l+1)=alim_star(ig,l)
	125	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1) &
	126	& *(zlev(ig,l+1)-zlev(ig,l)) &
	127	& 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
	128	w_est(ig,l+1)=zw2(ig,l+1)
	129	!
	130	else if ((zw2(ig,l).ge.1e-10).and. &
	131	& (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
	132	!estimation du detrainement a partir de la geometrie du pas precedent
	133	!tests sur la definition du detr
	134	!calcul de detr_star et entr_star
[852]	135
	136
	137
	138	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	139	! FH le test miraculeux de Catherine ? Le bout du tunel ?
	140	! w_est(ig,3)=zw2(ig,2)* &
	141	! & ((f_star(ig,2))**2) &
	142	! & /(f_star(ig,2)+alim_star(ig,2))**2+ &
	143	! & 2.RG(ztva(ig,1)-ztv(ig,2))/ztv(ig,2) &
	144	! & *(zlev(ig,3)-zlev(ig,2))
	145	! if (l.gt.2) then
	146	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	147
[814]	148	if (l.gt.2) then
[852]	149	w_est(ig,3)=zw2(ig,2)* &
	150	& ((f_star(ig,2))**2) &
	151	& /(f_star(ig,2)+alim_star(ig,2))**2+ &
	152	& 2.RG(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
	153	& *(zlev(ig,3)-zlev(ig,2))
	154
	155
[814]	156	!---------------------------------------------------------------------------
	157	!calcul de l entrainement et du detrainement lateral
	158	!---------------------------------------------------------------------------
	159	!
	160	!test:estimation de ztva_new_est sans entrainement
	161	Tbef=ztla(ig,l-1)*zpspsk(ig,l)
	162	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	163	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
	164	qsatbef=MIN(0.5,qsatbef)
	165	zcor=1./(1.-retv*qsatbef)
	166	qsatbef=qsatbef*zcor
	167	Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
	168	if (Zsat) then
	169	qlbef=max(0.,zqta(ig,l-1)-qsatbef)
	170	DT = 0.5RLvCpqlbef
	171	do while (abs(DT).gt.DDT0)
	172	Tbef=Tbef+DT
	173	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	174	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
	175	qsatbef=MIN(0.5,qsatbef)
	176	zcor=1./(1.-retv*qsatbef)
	177	qsatbef=qsatbef*zcor
	178	qlbef=zqta(ig,l-1)-qsatbef
	179
	180	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	181	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
	182	zcor=1./(1.-retv*qsatbef)
	183	dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
	184	num=-Tbef+ztla(ig,l-1)zpspsk(ig,l)+RLvCpqlbef
	185	denom=1.+RLvCp*dqsat_dT
	186	DT=num/denom
	187	enddo
	188	zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef)
	189	endif
	190	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
	191	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
	192	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
	193	& -zqla_est(ig,l))-zqla_est(ig,l))
	194
	195	w_est(ig,l+1)=zw2(ig,l)* &
	196	& ((f_star(ig,l))**2) &
	197	& /(f_star(ig,l)+alim_star(ig,l))**2+ &
	198	& 2.RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
	199	& *(zlev(ig,l+1)-zlev(ig,l))
	200	if (w_est(ig,l+1).lt.0.) then
	201	w_est(ig,l+1)=zw2(ig,l)
	202	endif
	203	!
	204	!calcul du detrainement
	205	if ((w_est(ig,l+1).gt.w_est(ig,l)).and. &
	206	& (zlev(ig,l+1).lt.zmax_sec(ig)).and. &
	207	& (zqla(ig,l-1).lt.1.e-10)) then
	208	detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1) &
	209	& sqrt(w_est(ig,l+1))sqrt(l_mix*zlev(ig,l+1)) &
	210	& -rhobarz(ig,l)sqrt(w_est(ig,l))sqrt(l_mix*zlev(ig,l))) &
	211	& /(r_aspect*zmax_sec(ig)))
	212	if (lev_out.ge.20) print*,'coucou calcul detr 1'
	213	else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and. &
	214	& (zqla(ig,l-1).lt.1.e-10)) then
	215	detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig)) &
	216	& /(rhobarz(ig,lmix(ig))wmaxa(ig)) &
	217	& (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1)) &
	218	& *((zmax_sec(ig)-zlev(ig,l+1))/ &
	219	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2. &
	220	& -rhobarz(ig,l)*sqrt(w_est(ig,l)) &
	221	& *((zmax_sec(ig)-zlev(ig,l))/ &
	222	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
	223	if (lev_out.ge.20) print*,'coucou calcul detr 2'
	224	else
	225	detr_star(ig,l)=0.002f0(ig)f_star(ig,l) &
	226	& *(zlev(ig,l+1)-zlev(ig,l))
	227	if (lev_out.ge.20) print*,'coucou calcul detr 3'
	228
	229	endif
	230	detr_star(ig,l)=detr_star(ig,l)/f0(ig)
	231	!
	232	!calcul de entr_star
	233	!
[852]	234	! Deplacement du calcul de entr_star pour eviter d'avoir aussi
	235	! entr_star > fstar.
	236	! Redeplacer suite a la transformation du cas detr>f
	237	! FH
	238	if(l.gt.lalim(ig)) then
	239	entr_star(ig,l)=0.4*detr_star(ig,l)
	240	else
	241
	242	! FH :
	243	! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
	244	! en haut de la couche d'alimentation.
	245	! A remettre en questoin a la premiere occasion mais ca peut aider a
	246	! ecrire un code robuste.
	247	! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
	248	! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
	249	! d'alimentation, ce qui n'est pas forcement heureux.
	250	entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
	251	detr_star(ig,l)=entr_star(ig,l)
	252	endif
	253
	254	!
[814]	255	if (detr_star(ig,l).gt.f_star(ig,l)) then
[852]	256
	257	! Ce test est là entre autres parce qu'on passe par des valeurs
	258	! delirantes de detr_star.
	259	! ca vaut sans doute le coup de verifier pourquoi.
	260
[814]	261	detr_star(ig,l)=f_star(ig,l)
[852]	262	if (l.gt.lalim(ig)+1) then
	263	entr_star(ig,l)=0.
	264	alim_star(ig,l)=0.
	265	! FH ajout pour forcer a stoper le thermique juste sous le sommet
	266	! de la couche (voir calcul de finter)
	267	zw2(ig,l+1)=-1.e-10
	268	linter(ig)=l+1
	269	else
	270	entr_star(ig,l)=detr_star(ig,l)
	271	endif
[814]	272	endif
	273
	274	else
	275	detr_star(ig,l)=0.
	276	entr_star(ig,l)=0.
	277	endif
[852]	278
	279	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	280	! FH inutile si on conserve comme on l'a fait plus haut entr=detr
	281	! dans la couche d'alimentation
[814]	282	!pas d entrainement dans la couche alim
[852]	283	! if ((l.le.lalim(ig))) then
	284	! entr_star(ig,l)=0.
	285	! endif
	286	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[814]	287	!
	288	!prise en compte du detrainement et de l entrainement dans le calcul du flux
	289	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	290	& -detr_star(ig,l)
	291
	292	!test sur le signe de f_star
	293	if (f_star(ig,l+1).gt.1.e-10) then
	294	!----------------------------------------------------------------------------
	295	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	296	!---------------------------------------------------------------------------
	297	!
	298	Zsat=.false.
	299	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
	300	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	301	& /(f_star(ig,l+1)+detr_star(ig,l))
	302	!
	303	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
	304	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	305	& /(f_star(ig,l+1)+detr_star(ig,l))
	306	!
	307	Tbef=ztla(ig,l)*zpspsk(ig,l)
	308	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	309	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
	310	qsatbef=MIN(0.5,qsatbef)
	311	zcor=1./(1.-retv*qsatbef)
	312	qsatbef=qsatbef*zcor
	313	Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
	314	if (Zsat) then
	315	qlbef=max(0.,zqta(ig,l)-qsatbef)
	316	DT = 0.5RLvCpqlbef
	317	do while (abs(DT).gt.DDT0)
	318	Tbef=Tbef+DT
	319	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	320	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
	321	qsatbef=MIN(0.5,qsatbef)
	322	zcor=1./(1.-retv*qsatbef)
	323	qsatbef=qsatbef*zcor
	324	qlbef=zqta(ig,l)-qsatbef
	325
	326	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	327	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
	328	zcor=1./(1.-retv*qsatbef)
	329	dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
	330	num=-Tbef+ztla(ig,l)zpspsk(ig,l)+RLvCpqlbef
	331	denom=1.+RLvCp*dqsat_dT
	332	DT=num/denom
	333	enddo
	334	zqla(ig,l) = max(0.,qlbef)
	335	endif
	336	!
	337	! on ecrit de maniere conservative (sat ou non)
	338	! T = Tl +Lv/Cp ql
	339	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	340	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	341	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	342	zha(ig,l) = ztva(ig,l)
	343	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	344	& -zqla(ig,l))-zqla(ig,l))
	345
	346	!on ecrit zqsat
	347	zqsatth(ig,l)=qsatbef
	348	!calcul de vitesse
	349	zw2(ig,l+1)=zw2(ig,l)* &
	350	& ((f_star(ig,l))**2) &
	351	& /(f_star(ig,l+1)+detr_star(ig,l))**2+ &
	352	& 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
	353	& *(zlev(ig,l+1)-zlev(ig,l))
	354
	355	endif
	356	endif
	357	!
	358	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
[852]	359
	360	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	361	! stop'On tombe sur le cas particulier de thermcell_dry'
	362	print*,'On tombe sur le cas particulier de thermcell_plume'
	363	zw2(ig,l+1)=0.
	364	linter(ig)=l+1
	365	endif
	366
	367
[814]	368	if (zw2(ig,l+1).lt.0.) then
	369	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	370	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	371	zw2(ig,l+1)=0.
[852]	372	endif
	373
[814]	374	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
[852]	375
[814]	376	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	377	! lmix est le niveau de la couche ou w (wa_moy) est maximum
	378	lmix(ig)=l+1
	379	wmaxa(ig)=wa_moy(ig,l+1)
	380	endif
	381	enddo
	382	enddo
	383
	384
	385	return
	386	end

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

Download in other formats:

Original Format