Context Navigation

thermcell_plume.F90 @ 5225

Last change on this file since 5225 was 1399, checked in by idelkadi, 14 years ago

Modifications relatives a l'inclusion des nouveaux nuages d'Arnaud Jam

calltherm.F90 : remontee de variables suplementaires.

Passage de ztv, zpspsk, zthla, zthl en argument

fisrtilp.F : le programme de nuages avec appel a cloudth

Ajout des arguments : zqta, fraca, ztv, zpspsk, ztla, ztlh,
iflag_cldcon

physiq.F : appel modifie a calltherm et fisrtilp
thermcell_main.F90 : changement dans les flag_thermals_ed
thermcell_plume.F90 : version modifiee entrainement et detrainement

(on retrouve la precendente pour iflag_thermals_ed=10)

cloudth.F90 : PDF d'eau bi-gaussiennes

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 27.5 KB

Rev	Line
[1294]	1	!
[1299]	2	! $Id: thermcell_plume.F90 1399 2010-06-08 08:29:11Z abarral $
[1294]	3	!
[972]	4	SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
[1294]	5	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	6	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	7	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
[1026]	8	& ,lev_out,lunout1,igout)
[878]	9
	10	!--------------------------------------------------------------------------
	11	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
	12	!--------------------------------------------------------------------------
	13
	14	IMPLICIT NONE
	15
	16	#include "YOMCST.h"
	17	#include "YOETHF.h"
	18	#include "FCTTRE.h"
[938]	19	#include "iniprint.h"
[1026]	20	#include "thermcell.h"
[878]	21
[972]	22	INTEGER itap
[1026]	23	INTEGER lunout1,igout
[878]	24	INTEGER ngrid,klev
[972]	25	REAL ptimestep
[878]	26	REAL ztv(ngrid,klev)
	27	REAL zthl(ngrid,klev)
	28	REAL po(ngrid,klev)
	29	REAL zl(ngrid,klev)
	30	REAL rhobarz(ngrid,klev)
	31	REAL zlev(ngrid,klev+1)
	32	REAL pplev(ngrid,klev+1)
	33	REAL pphi(ngrid,klev)
	34	REAL zpspsk(ngrid,klev)
	35	REAL alim_star(ngrid,klev)
	36	REAL f0(ngrid)
	37	INTEGER lalim(ngrid)
	38	integer lev_out ! niveau pour les print
[972]	39	real zcon2(ngrid)
[1026]	40
	41	real alim_star_tot(ngrid)
[878]	42
	43	REAL ztva(ngrid,klev)
	44	REAL ztla(ngrid,klev)
	45	REAL zqla(ngrid,klev)
	46	REAL zqta(ngrid,klev)
	47	REAL zha(ngrid,klev)
	48
	49	REAL detr_star(ngrid,klev)
[972]	50	REAL coefc
[878]	51	REAL entr_star(ngrid,klev)
	52	REAL detr(ngrid,klev)
	53	REAL entr(ngrid,klev)
	54
[1294]	55	REAL csc(ngrid,klev)
	56
[878]	57	REAL zw2(ngrid,klev+1)
	58	REAL w_est(ngrid,klev+1)
	59	REAL f_star(ngrid,klev+1)
	60	REAL wa_moy(ngrid,klev+1)
	61
	62	REAL ztva_est(ngrid,klev)
	63	REAL zqla_est(ngrid,klev)
	64	REAL zqsatth(ngrid,klev)
[1026]	65	REAL zta_est(ngrid,klev)
[1294]	66	REAL zdw2
	67	REAL zw2modif
	68	REAL zeps
[878]	69
	70	REAL linter(ngrid)
	71	INTEGER lmix(ngrid)
[972]	72	INTEGER lmix_bis(ngrid)
[878]	73	REAL wmaxa(ngrid)
	74
	75	INTEGER ig,l,k
	76
[1399]	77	real zdz,zfact,zbuoy,zalpha,zdrag
[878]	78	real zcor,zdelta,zcvm5,qlbef
[1399]	79	real Tbef,qsatbef
	80	real dqsat_dT,DT,num,denom
[878]	81	REAL REPS,RLvCp,DDT0
	82	PARAMETER (DDT0=.01)
	83	logical Zsat
[1294]	84	LOGICAL active(ngrid),activetmp(ngrid)
	85	REAL fact_gamma,fact_epsilon,fact_gamma2
[1026]	86	REAL c2(ngrid,klev)
[1399]	87	REAL a1,m
[878]	88
[1294]	89	REAL zw2fact,expa
[878]	90	Zsat=.false.
	91	! Initialisation
	92	RLvCp = RLVTT/RCPD
	93
[1399]	94
[1294]	95	fact_epsilon=0.002
[1399]	96	a1=2./3.
[1294]	97	fact_gamma=0.9
[1399]	98	zfact=fact_gamma/(1+fact_gamma)
	99	fact_gamma2=zfact
	100	expa=0.
	101
	102	print*,'THERM 31H '
[1294]	103
	104	! Initialisations des variables reeles
[1399]	105	if (1==1) then
[1294]	106	ztva(:,:)=ztv(:,:)
	107	ztva_est(:,:)=ztva(:,:)
	108	ztla(:,:)=zthl(:,:)
	109	zqta(:,:)=po(:,:)
	110	zha(:,:) = ztva(:,:)
	111	else
	112	ztva(:,:)=0.
	113	ztva_est(:,:)=0.
	114	ztla(:,:)=0.
	115	zqta(:,:)=0.
	116	zha(:,:) =0.
	117	endif
[972]	118
[1294]	119	zqla_est(:,:)=0.
	120	zqsatth(:,:)=0.
	121	zqla(:,:)=0.
	122	detr_star(:,:)=0.
	123	entr_star(:,:)=0.
	124	alim_star(:,:)=0.
	125	alim_star_tot(:)=0.
	126	csc(:,:)=0.
	127	detr(:,:)=0.
	128	entr(:,:)=0.
	129	zw2(:,:)=0.
	130	w_est(:,:)=0.
	131	f_star(:,:)=0.
	132	wa_moy(:,:)=0.
	133	linter(:)=1.
	134	linter(:)=1.
[972]	135
[1294]	136	! Initialisation des variables entieres
	137	lmix(:)=1
	138	lmix_bis(:)=2
	139	wmaxa(:)=0.
	140	lalim(:)=1
	141
	142	!-------------------------------------------------------------------------
	143	! On ne considere comme actif que les colonnes dont les deux premieres
	144	! couches sont instables.
	145	!-------------------------------------------------------------------------
	146	active(:)=ztv(:,1)>ztv(:,2)
	147
	148	!-------------------------------------------------------------------------
	149	! Definition de l'alimentation a l'origine dans thermcell_init
	150	!-------------------------------------------------------------------------
	151	do l=1,klev-1
[878]	152	do ig=1,ngrid
[1294]	153	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
	154	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
	155	& *sqrt(zlev(ig,l+1))
[1399]	156	lalim(:)=l+1
[1294]	157	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	158	endif
[878]	159	enddo
	160	enddo
[1294]	161	do l=1,klev
	162	do ig=1,ngrid
	163	if (alim_star_tot(ig) > 1.e-10 ) then
	164	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
	165	endif
	166	enddo
[878]	167	enddo
[1294]	168	alim_star_tot(:)=1.
[878]	169
[972]	170
[1294]	171	!------------------------------------------------------------------------------
	172	! Calcul dans la premiere couche
	173	! On decide dans cette version que le thermique n'est actif que si la premiere
	174	! couche est instable.
[1399]	175	! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
[1294]	176	! dans une couche l>1
	177	!------------------------------------------------------------------------------
	178	do ig=1,ngrid
[972]	179	! Le panache va prendre au debut les caracteristiques de l'air contenu
	180	! dans cette couche.
[1294]	181	if (active(ig)) then
	182	ztla(ig,1)=zthl(ig,1)
	183	zqta(ig,1)=po(ig,1)
	184	zqla(ig,1)=zl(ig,1)
	185	!cr: attention, prise en compte de f*(1)=1
	186	f_star(ig,2)=alim_star(ig,1)
	187	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
	188	& *(zlev(ig,2)-zlev(ig,1)) &
	189	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
	190	w_est(ig,2)=zw2(ig,2)
	191	endif
	192	enddo
[878]	193	!
[972]	194
[1294]	195	!==============================================================================
	196	!boucle de calcul de la vitesse verticale dans le thermique
	197	!==============================================================================
	198	do l=2,klev-1
	199	!==============================================================================
[972]	200
[878]	201
[1294]	202	! On decide si le thermique est encore actif ou non
	203	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
	204	do ig=1,ngrid
	205	active(ig)=active(ig) &
	206	& .and. zw2(ig,l)>1.e-10 &
	207	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
	208	enddo
[878]	209
	210
	211
[972]	212	! Premier calcul de la vitesse verticale a partir de la temperature
	213	! potentielle virtuelle
[1294]	214	! if (1.eq.1) then
	215	! w_est(ig,3)=zw2(ig,2)* &
	216	! & ((f_star(ig,2))**2) &
	217	! & /(f_star(ig,2)+alim_star(ig,2))**2+ &
	218	! & 2.RG(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
	219	! & *(zlev(ig,3)-zlev(ig,2))
	220	! endif
[972]	221
	222
[878]	223	!---------------------------------------------------------------------------
[1294]	224	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
	225	! sans tenir compte du detrainement et de l'entrainement dans cette
	226	! couche
	227	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
	228	! avant) a l'alimentation pour avoir un calcul plus propre
[878]	229	!---------------------------------------------------------------------------
[972]	230
[1399]	231	call thermcell_condens(ngrid,active, &
	232	& zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
[878]	233
[1294]	234	do ig=1,ngrid
	235	if(active(ig)) then
[878]	236	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
[1294]	237	zta_est(ig,l)=ztva_est(ig,l)
[878]	238	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
	239	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
	240	& -zqla_est(ig,l))-zqla_est(ig,l))
	241
[1399]	242	if (1.eq.0) then
	243	!calcul de w_est sans prendre en compte le drag
[878]	244	w_est(ig,l+1)=zw2(ig,l)* &
	245	& ((f_star(ig,l))**2) &
	246	& /(f_star(ig,l)+alim_star(ig,l))**2+ &
	247	& 2.RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
	248	& *(zlev(ig,l+1)-zlev(ig,l))
[1399]	249	else
	250
	251	zdz=zlev(ig,l+1)-zlev(ig,l)
	252	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
	253	zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	254	zdrag=fact_epsilon/(zalpha**expa)
	255	zw2fact=zbuoy/zdrag*a1
	256	w_est(ig,l+1)=(w_est(ig,l)-zw2fact)exp(-2.zdrag/(1+fact_gamma)*zdz) &
	257	& +zw2fact
	258
	259	endif
	260
[878]	261	if (w_est(ig,l+1).lt.0.) then
	262	w_est(ig,l+1)=zw2(ig,l)
	263	endif
[1294]	264	endif
	265	enddo
[972]	266
[1294]	267	!-------------------------------------------------
	268	!calcul des taux d'entrainement et de detrainement
	269	!-------------------------------------------------
[972]	270
[1294]	271	do ig=1,ngrid
	272	if (active(ig)) then
	273	zdz=zlev(ig,l+1)-zlev(ig,l)
	274	zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	275
	276	! estimation de la fraction couverte par les thermiques
	277	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
[972]	278
[1399]	279	!calcul de la soumission papier
[1294]	280	! Calcul du taux d'entrainement entr_star (epsilon)
	281	entr_star(ig,l)=f_star(ig,l)zdz ( zfact * MAX(0., &
[1399]	282	& a1*zbuoy/w_est(ig,l+1) &
	283	& - fact_epsilon/zalpha**expa ) &
	284	& +0. )
	285
	286	!calcul du taux de detrainment (delta)
[1294]	287	! detr_star(ig,l)=f_star(ig,l)zdz ( &
[1399]	288	! & MAX(1.e-3, &
	289	! & -fact_gamma2a1zbuoy/w_est(ig,l+1) &
	290	! & +0.01(max(zqta(ig,l-1)-po(ig,l),0.)/(po(ig,l))/(w_est(ig,l+1)))*0.5 &
	291	! & +0. ))
[972]	292
[1399]	293	m=0.5
[972]	294
[1399]	295	detr_star(ig,l)=1.f_star(ig,l)zdz * &
	296	& MAX(5.e-4,-fact_gamma2a1(1./w_est(ig,l+1))((1.-(1.-m)/(1.+70zqta(ig,l-1)))*zbuoy &
	297	& -40(1.-m)(max(zqta(ig,l-1)-po(ig,l),0.))/(1.+70*zqta(ig,l-1)) ) )
[972]	298
[1399]	299	! detr_star(ig,l)=f_star(ig,l)zdz ( &
	300	! & MAX(0.0, &
	301	! & -fact_gamma2a1zbuoy/w_est(ig,l+1) &
	302	! & +20(max(zqta(ig,l-1)-po(ig,l),0.))1(zalpha/w_est(ig,l+1))**0.5 &
	303	! & +0. ))
[972]	304
	305
[1294]	306	! En dessous de lalim, on prend le max de alim_star et entr_star pour
	307	! alim_star et 0 sinon
	308	if (l.lt.lalim(ig)) then
	309	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
	310	entr_star(ig,l)=0.
[878]	311	endif
	312
[1399]	313	!attention test
	314	! if (detr_star(ig,l).gt.(f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l))) then
	315	! detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
	316	! endif
[1294]	317	! Calcul du flux montant normalise
	318	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	319	& -detr_star(ig,l)
[972]	320
[878]	321	endif
[1294]	322	enddo
[878]	323
	324	!----------------------------------------------------------------------------
	325	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	326	!---------------------------------------------------------------------------
[1294]	327	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
	328	do ig=1,ngrid
	329	if (activetmp(ig)) then
	330	Zsat=.false.
	331	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
[878]	332	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	333	& /(f_star(ig,l+1)+detr_star(ig,l))
[1294]	334	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
[878]	335	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	336	& /(f_star(ig,l+1)+detr_star(ig,l))
	337
[1294]	338	endif
	339	enddo
	340
[1399]	341	call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
[1294]	342
[1399]	343
[1294]	344	do ig=1,ngrid
	345	if (activetmp(ig)) then
[972]	346	if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
[878]	347	! on ecrit de maniere conservative (sat ou non)
	348	! T = Tl +Lv/Cp ql
	349	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	350	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	351	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	352	zha(ig,l) = ztva(ig,l)
	353	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	354	& -zqla(ig,l))-zqla(ig,l))
	355
	356	!on ecrit zqsat
[1399]	357	zqsatth(ig,l)=qsatbef
[1294]	358
	359	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	360	! zw2(ig,l+1)=&
	361	! & zw2(ig,l)(1-fact_epsilon/(1.+fact_gamma)2.*(zlev(ig,l+1)-zlev(ig,l))) &
	362	! & +2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
	363	! & *1./(1.+fact_gamma) &
[1026]	364	! & *(zlev(ig,l+1)-zlev(ig,l))
[1294]	365	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	366	! La meme en plus modulaire :
	367	zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	368	zdz=zlev(ig,l+1)-zlev(ig,l)
[1026]	369
	370
[1294]	371	zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
	372
	373	if (1==0) then
	374	zw2modif=zw2(ig,l)(1-fact_epsilon/(1.+fact_gamma)2.*zdz)
	375	zdw2=2.zbuoy/(1.+fact_gamma)zdz
	376	zw2(ig,l+1)=zw2modif+zdw2
	377	else
[1399]	378	zdrag=fact_epsilon/(zalpha**expa)
	379	zw2fact=zbuoy/zdrag*a1
	380	zw2(ig,l+1)=(zw2(ig,l)-zw2fact)exp(-2.zdrag/(1+fact_gamma)*zdz) &
[1294]	381	& +zw2fact
	382
	383
	384	endif
[1399]	385
[1294]	386	endif
	387	enddo
	388
[972]	389	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
[878]	390	!
[1294]	391	!---------------------------------------------------------------------------
[878]	392	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
[1294]	393	!---------------------------------------------------------------------------
[878]	394
[1294]	395	do ig=1,ngrid
[878]	396	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	397	! stop'On tombe sur le cas particulier de thermcell_dry'
[1399]	398	print*,'On tombe sur le cas particulier de thermcell_plume'
[878]	399	zw2(ig,l+1)=0.
	400	linter(ig)=l+1
	401	endif
	402
	403	if (zw2(ig,l+1).lt.0.) then
	404	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	405	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	406	zw2(ig,l+1)=0.
	407	endif
	408
	409	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	410
	411	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	412	! lmix est le niveau de la couche ou w (wa_moy) est maximum
[1026]	413	!on rajoute le calcul de lmix_bis
	414	if (zqla(ig,l).lt.1.e-10) then
	415	lmix_bis(ig)=l+1
	416	endif
[878]	417	lmix(ig)=l+1
	418	wmaxa(ig)=wa_moy(ig,l+1)
	419	endif
[1294]	420	enddo
	421
	422	!=========================================================================
	423	! FIN DE LA BOUCLE VERTICALE
[878]	424	enddo
[1294]	425	!=========================================================================
[878]	426
[1294]	427	!on recalcule alim_star_tot
	428	do ig=1,ngrid
	429	alim_star_tot(ig)=0.
	430	enddo
	431	do ig=1,ngrid
	432	do l=1,lalim(ig)-1
	433	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	434	enddo
	435	enddo
	436
	437
[972]	438	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
[878]	439
[972]	440	! print*,'thermcell_plume OK'
	441
[878]	442	return
	443	end
[1311]	444
	445	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	446	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	447	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	448	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	449	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	450	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	451	SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
	452	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	453	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	454	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
	455	& ,lev_out,lunout1,igout)
	456
	457	!--------------------------------------------------------------------------
	458	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
	459	! Version conforme a l'article de Rio et al. 2010.
	460	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
	461	!--------------------------------------------------------------------------
	462
	463	IMPLICIT NONE
	464
	465	#include "YOMCST.h"
	466	#include "YOETHF.h"
	467	#include "FCTTRE.h"
	468	#include "iniprint.h"
	469	#include "thermcell.h"
	470
	471	INTEGER itap
	472	INTEGER lunout1,igout
	473	INTEGER ngrid,klev
	474	REAL ptimestep
	475	REAL ztv(ngrid,klev)
	476	REAL zthl(ngrid,klev)
	477	REAL po(ngrid,klev)
	478	REAL zl(ngrid,klev)
	479	REAL rhobarz(ngrid,klev)
	480	REAL zlev(ngrid,klev+1)
	481	REAL pplev(ngrid,klev+1)
	482	REAL pphi(ngrid,klev)
	483	REAL zpspsk(ngrid,klev)
	484	REAL alim_star(ngrid,klev)
	485	REAL f0(ngrid)
	486	INTEGER lalim(ngrid)
	487	integer lev_out ! niveau pour les print
	488
	489	real alim_star_tot(ngrid)
	490
	491	REAL ztva(ngrid,klev)
	492	REAL ztla(ngrid,klev)
	493	REAL zqla(ngrid,klev)
	494	REAL zqta(ngrid,klev)
	495	REAL zha(ngrid,klev)
	496
	497	REAL detr_star(ngrid,klev)
	498	REAL coefc
	499	REAL entr_star(ngrid,klev)
	500	REAL detr(ngrid,klev)
	501	REAL entr(ngrid,klev)
	502
	503	REAL csc(ngrid,klev)
	504
	505	REAL zw2(ngrid,klev+1)
	506	REAL w_est(ngrid,klev+1)
	507	REAL f_star(ngrid,klev+1)
	508	REAL wa_moy(ngrid,klev+1)
	509
	510	REAL ztva_est(ngrid,klev)
	511	REAL zqla_est(ngrid,klev)
	512	REAL zqsatth(ngrid,klev)
	513	REAL zta_est(ngrid,klev)
[1338]	514	REAL ztemp(ngrid),zqsat(ngrid)
[1311]	515	REAL zdw2
	516	REAL zw2modif
	517	REAL zw2fact
	518	REAL zeps(ngrid,klev)
	519
	520	REAL linter(ngrid)
	521	INTEGER lmix(ngrid)
	522	INTEGER lmix_bis(ngrid)
	523	REAL wmaxa(ngrid)
	524
	525	INTEGER ig,l,k
	526
	527	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
	528	real zbuoybis
	529	real zcor,zdelta,zcvm5,qlbef,zdz2
	530	real betalpha,zbetalpha
[1338]	531	real eps, afact
[1311]	532	REAL REPS,RLvCp,DDT0
	533	PARAMETER (DDT0=.01)
	534	logical Zsat
	535	LOGICAL active(ngrid),activetmp(ngrid)
	536	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
	537	REAL c2(ngrid,klev)
	538	Zsat=.false.
	539	! Initialisation
	540
	541	RLvCp = RLVTT/RCPD
	542	fact_epsilon=0.002
	543	betalpha=0.9
	544	afact=2./3.
	545
	546	zbetalpha=betalpha/(1.+betalpha)
	547
	548	! print*,'THERM 31B'
[1371]	549	print*,'THERMCELL_PLUME OPTIMISE OOOO'
[1311]	550
	551	! Initialisations des variables reeles
	552	if (1==0) then
	553	ztva(:,:)=ztv(:,:)
	554	ztva_est(:,:)=ztva(:,:)
	555	ztla(:,:)=zthl(:,:)
	556	zqta(:,:)=po(:,:)
	557	zha(:,:) = ztva(:,:)
	558	else
	559	ztva(:,:)=0.
	560	ztva_est(:,:)=0.
	561	ztla(:,:)=0.
	562	zqta(:,:)=0.
	563	zha(:,:) =0.
	564	endif
	565
	566	zqla_est(:,:)=0.
	567	zqsatth(:,:)=0.
	568	zqla(:,:)=0.
	569	detr_star(:,:)=0.
	570	entr_star(:,:)=0.
	571	alim_star(:,:)=0.
	572	alim_star_tot(:)=0.
	573	csc(:,:)=0.
	574	detr(:,:)=0.
	575	entr(:,:)=0.
	576	zw2(:,:)=0.
	577	zbuoy(:,:)=0.
	578	gamma(:,:)=0.
	579	zeps(:,:)=0.
	580	w_est(:,:)=0.
	581	f_star(:,:)=0.
	582	wa_moy(:,:)=0.
	583	linter(:)=1.
	584	! linter(:)=1.
	585	! Initialisation des variables entieres
	586	lmix(:)=1
	587	lmix_bis(:)=2
	588	wmaxa(:)=0.
	589	lalim(:)=1
	590
[1338]	591	print*,'THERMCELL PLUME QSAT2 NDDDDN'
[1311]	592
	593	!-------------------------------------------------------------------------
	594	! On ne considere comme actif que les colonnes dont les deux premieres
	595	! couches sont instables.
	596	!-------------------------------------------------------------------------
	597	active(:)=ztv(:,1)>ztv(:,2)
	598
	599	!-------------------------------------------------------------------------
	600	! Definition de l'alimentation a l'origine dans thermcell_init
	601	!-------------------------------------------------------------------------
	602	do l=1,klev-1
	603	do ig=1,ngrid
	604	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
	605	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
	606	& *sqrt(zlev(ig,l+1))
[1373]	607	lalim(ig)=l+1
[1311]	608	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	609	endif
	610	enddo
	611	enddo
	612	do l=1,klev
	613	do ig=1,ngrid
	614	if (alim_star_tot(ig) > 1.e-10 ) then
	615	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
	616	endif
	617	enddo
	618	enddo
	619	alim_star_tot(:)=1.
	620
	621
	622
	623	!------------------------------------------------------------------------------
	624	! Calcul dans la premiere couche
	625	! On decide dans cette version que le thermique n'est actif que si la premiere
	626	! couche est instable.
	627	! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
	628	! dans une couche l>1
	629	!------------------------------------------------------------------------------
	630	do ig=1,ngrid
	631	! Le panache va prendre au debut les caracteristiques de l'air contenu
	632	! dans cette couche.
	633	if (active(ig)) then
	634	ztla(ig,1)=zthl(ig,1)
	635	zqta(ig,1)=po(ig,1)
	636	zqla(ig,1)=zl(ig,1)
	637	!cr: attention, prise en compte de f*(1)=1
	638	f_star(ig,2)=alim_star(ig,1)
	639	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
	640	& *(zlev(ig,2)-zlev(ig,1)) &
	641	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
	642	w_est(ig,2)=zw2(ig,2)
	643	endif
	644	enddo
	645	!
	646
	647	!==============================================================================
	648	!boucle de calcul de la vitesse verticale dans le thermique
	649	!==============================================================================
	650	do l=2,klev-1
	651	!==============================================================================
	652
	653
	654	! On decide si le thermique est encore actif ou non
	655	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
	656	do ig=1,ngrid
	657	active(ig)=active(ig) &
	658	& .and. zw2(ig,l)>1.e-10 &
	659	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
	660	enddo
	661
	662
	663
	664	!---------------------------------------------------------------------------
	665	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
	666	! sans tenir compte du detrainement et de l'entrainement dans cette
	667	! couche
[1338]	668	! C'est a dire qu'on suppose
	669	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
[1311]	670	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
	671	! avant) a l'alimentation pour avoir un calcul plus propre
	672	!---------------------------------------------------------------------------
	673
[1338]	674	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
	675	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
[1311]	676
	677	do ig=1,ngrid
	678	! print*,'active',active(ig),ig,l
	679	if(active(ig)) then
[1338]	680	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
[1311]	681	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
	682	zta_est(ig,l)=ztva_est(ig,l)
	683	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
	684	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
	685	& -zqla_est(ig,l))-zqla_est(ig,l))
	686
	687	!------------------------------------------------
	688	!AJAM:nouveau calcul de w²
	689	!------------------------------------------------
	690	zdz=zlev(ig,l+1)-zlev(ig,l)
	691	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	692
	693	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	694	zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
	695	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
	696
	697
	698	if (w_est(ig,l+1).lt.0.) then
	699	w_est(ig,l+1)=zw2(ig,l)
	700	endif
	701	endif
	702	enddo
	703
	704
	705	!-------------------------------------------------
	706	!calcul des taux d'entrainement et de detrainement
	707	!-------------------------------------------------
	708
[1338]	709	print*,'THERM V1 SANS DQ'
[1311]	710	do ig=1,ngrid
	711	if (active(ig)) then
	712
	713	zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
	714	zw2m=w_est(ig,l+1)
	715	zdz=zlev(ig,l+1)-zlev(ig,l)
	716	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	717	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
	718	zbuoybis=zbuoy(ig,l)
	719	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
	720	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
	721
	722
	723	entr_star(ig,l)=f_star(ig,l)zdz zbetalpha*MAX(0., &
	724	& afact*zbuoybis/zw2m - fact_epsilon )
	725
	726
	727	detr_star(ig,l)=f_star(ig,l)*zdz &
	728	& MAX(1.e-3, -afactzbetalpha*zbuoy(ig,l)/zw2m &
	729	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
	730
	731	! En dessous de lalim, on prend le max de alim_star et entr_star pour
	732	! alim_star et 0 sinon
	733	if (l.lt.lalim(ig)) then
	734	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
	735	entr_star(ig,l)=0.
	736	endif
	737
	738	! Calcul du flux montant normalise
	739	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	740	& -detr_star(ig,l)
	741
	742	endif
	743	enddo
	744
	745
	746	!----------------------------------------------------------------------------
	747	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	748	!---------------------------------------------------------------------------
	749	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
	750	do ig=1,ngrid
	751	if (activetmp(ig)) then
	752	Zsat=.false.
	753	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
	754	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	755	& /(f_star(ig,l+1)+detr_star(ig,l))
	756	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
	757	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	758	& /(f_star(ig,l+1)+detr_star(ig,l))
	759
	760	endif
	761	enddo
	762
[1338]	763	ztemp(:)=zpspsk(:,l)*ztla(:,l)
	764	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
[1311]	765
	766	do ig=1,ngrid
	767	if (activetmp(ig)) then
	768	if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
	769	! on ecrit de maniere conservative (sat ou non)
	770	! T = Tl +Lv/Cp ql
[1338]	771	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
[1311]	772	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	773	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	774	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	775	zha(ig,l) = ztva(ig,l)
	776	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	777	& -zqla(ig,l))-zqla(ig,l))
	778	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	779	zdz=zlev(ig,l+1)-zlev(ig,l)
	780	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
	781
	782	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	783	zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
	784	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
	785	endif
	786	enddo
	787
	788	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
	789	!
	790	!---------------------------------------------------------------------------
	791	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
	792	!---------------------------------------------------------------------------
	793
	794	do ig=1,ngrid
	795	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	796	! stop'On tombe sur le cas particulier de thermcell_dry'
	797	print*,'On tombe sur le cas particulier de thermcell_plume'
	798	zw2(ig,l+1)=0.
	799	linter(ig)=l+1
	800	endif
	801
	802	if (zw2(ig,l+1).lt.0.) then
	803	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	804	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	805	zw2(ig,l+1)=0.
	806	endif
	807
	808	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	809
	810	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	811	! lmix est le niveau de la couche ou w (wa_moy) est maximum
	812	!on rajoute le calcul de lmix_bis
	813	if (zqla(ig,l).lt.1.e-10) then
	814	lmix_bis(ig)=l+1
	815	endif
	816	lmix(ig)=l+1
	817	wmaxa(ig)=wa_moy(ig,l+1)
	818	endif
	819	enddo
	820
	821	!=========================================================================
	822	! FIN DE LA BOUCLE VERTICALE
	823	enddo
	824	!=========================================================================
	825
	826	! print*,'THERMCELL PLUME ARNAUD DEDANS 7'
	827	!on recalcule alim_star_tot
	828	do ig=1,ngrid
	829	alim_star_tot(ig)=0.
	830	enddo
	831	do ig=1,ngrid
	832	do l=1,lalim(ig)-1
	833	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	834	enddo
	835	enddo
	836
	837
	838	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
	839
	840	#undef wrgrads_thermcell
	841	#ifdef wrgrads_thermcell
	842	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
	843	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
	844	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
	845	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
	846	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
	847	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
	848	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
	849	#endif
	850
	851
	852	! print*,'THERMCELL PLUME ARNAUD DEDANS 8'
	853	return
	854	end
	855

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

Download in other formats:

Original Format