Context Navigation

thermcell_plume.F90 @ 1997

Last change on this file since 1997 was 1982, checked in by fhourdin, 11 years ago
Suppression d'impressions
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: 28.1 KB

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

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

Download in other formats:

Original Format