Context Navigation

thermcell_plume.F90 @ 1664

Last change on this file since 1664 was 1503, checked in by idelkadi, 14 years ago
Correction de bug et vectorisation de thermcell_plume.F90
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.4 KB

Rev	Line
[1403]	1	!
	2	! $Id: thermcell_plume.F90 1503 2011-03-23 11:57:52Z fairhead $
	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	!--------------------------------------------------------------------------
	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
[1503]	39	integer nbpb
[972]	40	real zcon2(ngrid)
[1026]	41
	42	real alim_star_tot(ngrid)
[878]	43
	44	REAL ztva(ngrid,klev)
	45	REAL ztla(ngrid,klev)
	46	REAL zqla(ngrid,klev)
	47	REAL zqta(ngrid,klev)
	48	REAL zha(ngrid,klev)
	49
	50	REAL detr_star(ngrid,klev)
[972]	51	REAL coefc
[878]	52	REAL entr_star(ngrid,klev)
	53	REAL detr(ngrid,klev)
	54	REAL entr(ngrid,klev)
	55
[1403]	56	REAL csc(ngrid,klev)
	57
[878]	58	REAL zw2(ngrid,klev+1)
	59	REAL w_est(ngrid,klev+1)
	60	REAL f_star(ngrid,klev+1)
	61	REAL wa_moy(ngrid,klev+1)
	62
	63	REAL ztva_est(ngrid,klev)
	64	REAL zqla_est(ngrid,klev)
	65	REAL zqsatth(ngrid,klev)
[1026]	66	REAL zta_est(ngrid,klev)
[1403]	67	REAL zdw2
	68	REAL zw2modif
	69	REAL zeps
[878]	70
	71	REAL linter(ngrid)
	72	INTEGER lmix(ngrid)
[972]	73	INTEGER lmix_bis(ngrid)
[878]	74	REAL wmaxa(ngrid)
	75
	76	INTEGER ig,l,k
	77
[1403]	78	real zdz,zfact,zbuoy,zalpha,zdrag
[878]	79	real zcor,zdelta,zcvm5,qlbef
	80	real Tbef,qsatbef
	81	real dqsat_dT,DT,num,denom
	82	REAL REPS,RLvCp,DDT0
	83	PARAMETER (DDT0=.01)
	84	logical Zsat
[1403]	85	LOGICAL active(ngrid),activetmp(ngrid)
	86	REAL fact_gamma,fact_epsilon,fact_gamma2
[1026]	87	REAL c2(ngrid,klev)
[1403]	88	REAL a1,m
[878]	89
[1403]	90	REAL zw2fact,expa
[878]	91	Zsat=.false.
	92	! Initialisation
	93	RLvCp = RLVTT/RCPD
	94
[1403]	95
	96	fact_epsilon=0.002
	97	a1=2./3.
	98	fact_gamma=0.9
	99	zfact=fact_gamma/(1+fact_gamma)
	100	fact_gamma2=zfact
	101	expa=0.
	102
[1026]	103
[1403]	104	! Initialisations des variables reeles
	105	if (1==1) then
	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
[878]	118
[1403]	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
[1403]	136	! Initialisation des variables entieres
	137	lmix(:)=1
	138	lmix_bis(:)=2
	139	wmaxa(:)=0.
	140	lalim(:)=1
[972]	141
[1403]	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
[1403]	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))
[1503]	156	lalim(ig)=l+1
[1403]	157	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	158	endif
[878]	159	enddo
	160	enddo
[1403]	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
[1403]	168	alim_star_tot(:)=1.
[878]	169
[972]	170
[1403]	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.
	175	! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
	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.
[1403]	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
[1403]	195	!==============================================================================
	196	!boucle de calcul de la vitesse verticale dans le thermique
	197	!==============================================================================
	198	do l=2,klev-1
	199	!==============================================================================
[972]	200
[878]	201
[1403]	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
[1403]	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	!---------------------------------------------------------------------------
[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
	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
[1403]	231	call thermcell_condens(ngrid,active, &
	232	& zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
[878]	233
[1403]	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)
[1403]	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
[1403]	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))
[1403]	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
[1403]	264	endif
	265	enddo
[972]	266
[1403]	267	!-------------------------------------------------
	268	!calcul des taux d'entrainement et de detrainement
	269	!-------------------------------------------------
[972]	270
[1403]	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
[1403]	279	!calcul de la soumission papier
	280	! Calcul du taux d'entrainement entr_star (epsilon)
	281	entr_star(ig,l)=f_star(ig,l)zdz ( zfact * MAX(0., &
	282	& a1*zbuoy/w_est(ig,l+1) &
	283	& - fact_epsilon/zalpha**expa ) &
	284	& +0. )
	285
	286	!calcul du taux de detrainment (delta)
	287	! detr_star(ig,l)=f_star(ig,l)zdz ( &
	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
[1403]	293	m=0.5
[972]	294
[1403]	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
[1403]	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
[1403]	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.
	311	endif
[972]	312
[1403]	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
	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
[1403]	321	endif
	322	enddo
[972]	323
[1403]	324	!----------------------------------------------------------------------------
	325	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	326	!---------------------------------------------------------------------------
	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)+ &
	332	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	333	& /(f_star(ig,l+1)+detr_star(ig,l))
	334	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
	335	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	336	& /(f_star(ig,l+1)+detr_star(ig,l))
[972]	337
[1403]	338	endif
	339	enddo
[972]	340
[1403]	341	call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
[972]	342
	343
[1403]	344	do ig=1,ngrid
	345	if (activetmp(ig)) then
	346	! on ecrit de maniere conservative (sat ou non)
	347	! T = Tl +Lv/Cp ql
	348	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	349	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	350	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	351	zha(ig,l) = ztva(ig,l)
	352	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	353	& -zqla(ig,l))-zqla(ig,l))
[972]	354
[1403]	355	!on ecrit zqsat
	356	zqsatth(ig,l)=qsatbef
[972]	357
[1403]	358	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	359	! zw2(ig,l+1)=&
	360	! & zw2(ig,l)(1-fact_epsilon/(1.+fact_gamma)2.*(zlev(ig,l+1)-zlev(ig,l))) &
	361	! & +2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
	362	! & *1./(1.+fact_gamma) &
	363	! & *(zlev(ig,l+1)-zlev(ig,l))
	364	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	365	! La meme en plus modulaire :
	366	zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	367	zdz=zlev(ig,l+1)-zlev(ig,l)
[972]	368
	369
[1403]	370	zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
[972]	371
[1503]	372	!if (1==0) then
	373	! zw2modif=zw2(ig,l)(1-fact_epsilon/(1.+fact_gamma)2.*zdz)
	374	! zdw2=2.zbuoy/(1.+fact_gamma)zdz
	375	! zw2(ig,l+1)=zw2modif+zdw2
	376	!else
[1403]	377	zdrag=fact_epsilon/(zalpha**expa)
	378	zw2fact=zbuoy/zdrag*a1
	379	zw2(ig,l+1)=(zw2(ig,l)-zw2fact)exp(-2.zdrag/(1+fact_gamma)*zdz) &
	380	& +zw2fact
[1503]	381	!endif
[972]	382
[1403]	383	endif
	384	enddo
	385
	386	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
[878]	387	!
[1403]	388	!---------------------------------------------------------------------------
	389	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
	390	!---------------------------------------------------------------------------
[972]	391
[1503]	392	nbpb=0
[1403]	393	do ig=1,ngrid
	394	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	395	! stop'On tombe sur le cas particulier de thermcell_dry'
[1503]	396	! print*,'On tombe sur le cas particulier de thermcell_plume'
	397	nbpb=nbpb+1
[1403]	398	zw2(ig,l+1)=0.
	399	linter(ig)=l+1
	400	endif
[972]	401
[1403]	402	if (zw2(ig,l+1).lt.0.) then
	403	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	404	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	405	zw2(ig,l+1)=0.
	406	endif
[972]	407
[1403]	408	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
[972]	409
[1403]	410	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	411	! lmix est le niveau de la couche ou w (wa_moy) est maximum
	412	!on rajoute le calcul de lmix_bis
	413	if (zqla(ig,l).lt.1.e-10) then
	414	lmix_bis(ig)=l+1
	415	endif
	416	lmix(ig)=l+1
	417	wmaxa(ig)=wa_moy(ig,l+1)
	418	endif
	419	enddo
[972]	420
[1503]	421	if (nbpb>0) then
	422	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
	423	endif
	424
[1403]	425	!=========================================================================
	426	! FIN DE LA BOUCLE VERTICALE
	427	enddo
	428	!=========================================================================
[972]	429
[1403]	430	!on recalcule alim_star_tot
	431	do ig=1,ngrid
	432	alim_star_tot(ig)=0.
	433	enddo
	434	do ig=1,ngrid
	435	do l=1,lalim(ig)-1
	436	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	437	enddo
	438	enddo
	439
[972]	440
[1403]	441	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
[972]	442
[878]	443
[1403]	444	return
	445	end
[972]	446
[1403]	447	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	448	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	449	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	450	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	451	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	452	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	453	SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
	454	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	455	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	456	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
	457	& ,lev_out,lunout1,igout)
[972]	458
[1403]	459	!--------------------------------------------------------------------------
	460	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
	461	! Version conforme a l'article de Rio et al. 2010.
	462	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
	463	!--------------------------------------------------------------------------
[878]	464
[1403]	465	IMPLICIT NONE
[972]	466
[1403]	467	#include "YOMCST.h"
	468	#include "YOETHF.h"
	469	#include "FCTTRE.h"
	470	#include "iniprint.h"
	471	#include "thermcell.h"
[972]	472
[1403]	473	INTEGER itap
	474	INTEGER lunout1,igout
	475	INTEGER ngrid,klev
	476	REAL ptimestep
	477	REAL ztv(ngrid,klev)
	478	REAL zthl(ngrid,klev)
	479	REAL po(ngrid,klev)
	480	REAL zl(ngrid,klev)
	481	REAL rhobarz(ngrid,klev)
	482	REAL zlev(ngrid,klev+1)
	483	REAL pplev(ngrid,klev+1)
	484	REAL pphi(ngrid,klev)
	485	REAL zpspsk(ngrid,klev)
	486	REAL alim_star(ngrid,klev)
	487	REAL f0(ngrid)
	488	INTEGER lalim(ngrid)
	489	integer lev_out ! niveau pour les print
[1503]	490	integer nbpb
[1403]	491
	492	real alim_star_tot(ngrid)
[878]	493
[1403]	494	REAL ztva(ngrid,klev)
	495	REAL ztla(ngrid,klev)
	496	REAL zqla(ngrid,klev)
	497	REAL zqta(ngrid,klev)
	498	REAL zha(ngrid,klev)
[878]	499
[1403]	500	REAL detr_star(ngrid,klev)
	501	REAL coefc
	502	REAL entr_star(ngrid,klev)
	503	REAL detr(ngrid,klev)
	504	REAL entr(ngrid,klev)
	505
	506	REAL csc(ngrid,klev)
	507
	508	REAL zw2(ngrid,klev+1)
	509	REAL w_est(ngrid,klev+1)
	510	REAL f_star(ngrid,klev+1)
	511	REAL wa_moy(ngrid,klev+1)
	512
	513	REAL ztva_est(ngrid,klev)
	514	REAL zqla_est(ngrid,klev)
	515	REAL zqsatth(ngrid,klev)
	516	REAL zta_est(ngrid,klev)
	517	REAL ztemp(ngrid),zqsat(ngrid)
	518	REAL zdw2
	519	REAL zw2modif
	520	REAL zw2fact
	521	REAL zeps(ngrid,klev)
	522
	523	REAL linter(ngrid)
	524	INTEGER lmix(ngrid)
	525	INTEGER lmix_bis(ngrid)
	526	REAL wmaxa(ngrid)
	527
	528	INTEGER ig,l,k
	529
	530	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
	531	real zbuoybis
	532	real zcor,zdelta,zcvm5,qlbef,zdz2
	533	real betalpha,zbetalpha
	534	real eps, afact
	535	REAL REPS,RLvCp,DDT0
	536	PARAMETER (DDT0=.01)
	537	logical Zsat
	538	LOGICAL active(ngrid),activetmp(ngrid)
	539	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
	540	REAL c2(ngrid,klev)
	541	Zsat=.false.
	542	! Initialisation
	543
	544	RLvCp = RLVTT/RCPD
	545	fact_epsilon=0.002
	546	betalpha=0.9
	547	afact=2./3.
	548
	549	zbetalpha=betalpha/(1.+betalpha)
	550
	551
	552	! Initialisations des variables reeles
	553	if (1==0) then
	554	ztva(:,:)=ztv(:,:)
	555	ztva_est(:,:)=ztva(:,:)
	556	ztla(:,:)=zthl(:,:)
	557	zqta(:,:)=po(:,:)
	558	zha(:,:) = ztva(:,:)
	559	else
	560	ztva(:,:)=0.
	561	ztva_est(:,:)=0.
	562	ztla(:,:)=0.
	563	zqta(:,:)=0.
	564	zha(:,:) =0.
	565	endif
	566
	567	zqla_est(:,:)=0.
	568	zqsatth(:,:)=0.
	569	zqla(:,:)=0.
	570	detr_star(:,:)=0.
	571	entr_star(:,:)=0.
	572	alim_star(:,:)=0.
	573	alim_star_tot(:)=0.
	574	csc(:,:)=0.
	575	detr(:,:)=0.
	576	entr(:,:)=0.
	577	zw2(:,:)=0.
	578	zbuoy(:,:)=0.
	579	gamma(:,:)=0.
	580	zeps(:,:)=0.
	581	w_est(:,:)=0.
	582	f_star(:,:)=0.
	583	wa_moy(:,:)=0.
	584	linter(:)=1.
	585	! linter(:)=1.
	586	! Initialisation des variables entieres
	587	lmix(:)=1
	588	lmix_bis(:)=2
	589	wmaxa(:)=0.
	590	lalim(:)=1
	591
	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))
	607	lalim(ig)=l+1
	608	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
[878]	609	endif
[1403]	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.
[878]	620
	621
[972]	622
[1403]	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	!
[1026]	646
[1403]	647	!==============================================================================
	648	!boucle de calcul de la vitesse verticale dans le thermique
	649	!==============================================================================
	650	do l=2,klev-1
	651	!==============================================================================
[1026]	652
	653
[1403]	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
[1026]	661
	662
	663
[1403]	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
	668	! C'est a dire qu'on suppose
	669	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
	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
	674	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
	675	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
	676
	677	do ig=1,ngrid
	678	! print*,'active',active(ig),ig,l
	679	if(active(ig)) then
	680	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
	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
[1026]	701	endif
[1403]	702	enddo
[1026]	703
	704
[1403]	705	!-------------------------------------------------
	706	!calcul des taux d'entrainement et de detrainement
	707	!-------------------------------------------------
[1026]	708
[1403]	709	do ig=1,ngrid
	710	if (active(ig)) then
[1026]	711
[1403]	712	zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
	713	zw2m=w_est(ig,l+1)
	714	zdz=zlev(ig,l+1)-zlev(ig,l)
	715	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	716	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
	717	zbuoybis=zbuoy(ig,l)
	718	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
	719	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
[1026]	720
[1403]	721
	722	entr_star(ig,l)=f_star(ig,l)zdz zbetalpha*MAX(0., &
	723	& afact*zbuoybis/zw2m - fact_epsilon )
	724
	725
	726	detr_star(ig,l)=f_star(ig,l)*zdz &
	727	& MAX(1.e-3, -afactzbetalpha*zbuoy(ig,l)/zw2m &
	728	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
	729
	730	! En dessous de lalim, on prend le max de alim_star et entr_star pour
	731	! alim_star et 0 sinon
	732	if (l.lt.lalim(ig)) then
	733	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
	734	entr_star(ig,l)=0.
	735	endif
	736
	737	! Calcul du flux montant normalise
[878]	738	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	739	& -detr_star(ig,l)
	740
[1403]	741	endif
	742	enddo
	743
	744
[878]	745	!----------------------------------------------------------------------------
	746	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	747	!---------------------------------------------------------------------------
[1403]	748	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
	749	do ig=1,ngrid
	750	if (activetmp(ig)) then
	751	Zsat=.false.
	752	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
[878]	753	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	754	& /(f_star(ig,l+1)+detr_star(ig,l))
[1403]	755	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
[878]	756	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	757	& /(f_star(ig,l+1)+detr_star(ig,l))
	758
[1403]	759	endif
	760	enddo
	761
	762	ztemp(:)=zpspsk(:,l)*ztla(:,l)
	763	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
	764
	765	do ig=1,ngrid
	766	if (activetmp(ig)) then
[878]	767	! on ecrit de maniere conservative (sat ou non)
	768	! T = Tl +Lv/Cp ql
[1403]	769	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
[878]	770	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	771	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	772	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	773	zha(ig,l) = ztva(ig,l)
	774	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	775	& -zqla(ig,l))-zqla(ig,l))
[1403]	776	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	777	zdz=zlev(ig,l+1)-zlev(ig,l)
	778	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
[878]	779
[1403]	780	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	781	zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
	782	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
	783	endif
	784	enddo
[1026]	785
[972]	786	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
[878]	787	!
[1403]	788	!---------------------------------------------------------------------------
[878]	789	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
[1403]	790	!---------------------------------------------------------------------------
[878]	791
[1503]	792	nbpb=0
[1403]	793	do ig=1,ngrid
[878]	794	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	795	! stop'On tombe sur le cas particulier de thermcell_dry'
[1503]	796	! print*,'On tombe sur le cas particulier de thermcell_plume'
	797	nbpb=nbpb+1
[878]	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
[1026]	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
[878]	816	lmix(ig)=l+1
	817	wmaxa(ig)=wa_moy(ig,l+1)
	818	endif
[1403]	819	enddo
	820
[1503]	821	if (nbpb>0) then
	822	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
	823	endif
	824
[1403]	825	!=========================================================================
	826	! FIN DE LA BOUCLE VERTICALE
[878]	827	enddo
[1403]	828	!=========================================================================
[878]	829
[1403]	830	!on recalcule alim_star_tot
	831	do ig=1,ngrid
	832	alim_star_tot(ig)=0.
	833	enddo
	834	do ig=1,ngrid
	835	do l=1,lalim(ig)-1
	836	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	837	enddo
	838	enddo
	839
	840
[972]	841	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
[878]	842
[1403]	843	#undef wrgrads_thermcell
	844	#ifdef wrgrads_thermcell
	845	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
	846	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
	847	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
	848	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
	849	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
	850	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
	851	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
	852	#endif
[972]	853
[1403]	854
	855	return
	856	end
	857

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

Download in other formats:

Original Format