Context Navigation

thermcell_plume.F90 @ 409

Last change on this file since 409 was 179, checked in by lfita, 10 years ago
No Ale, Alp, ... values because pphi is at the intermediate level!!! Removing raw prints of checkig for Ale, Alp
File size: 27.7 KB

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

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

Download in other formats:

Original Format