Context Navigation

thermcell_plume.F90 @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 27.3 KB

Rev	Line
[1]	1	!
	2	! $Id: thermcell_plume.F90 1403 2010-07-01 09:02:53Z fairhead $
	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	real zcon2(ngrid)
	40
	41	real alim_star_tot(ngrid)
	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)
	50	REAL coefc
	51	REAL entr_star(ngrid,klev)
	52	REAL detr(ngrid,klev)
	53	REAL entr(ngrid,klev)
	54
	55	REAL csc(ngrid,klev)
	56
	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)
	65	REAL zta_est(ngrid,klev)
	66	REAL zdw2
	67	REAL zw2modif
	68	REAL zeps
	69
	70	REAL linter(ngrid)
	71	INTEGER lmix(ngrid)
	72	INTEGER lmix_bis(ngrid)
	73	REAL wmaxa(ngrid)
	74
	75	INTEGER ig,l,k
	76
	77	real zdz,zfact,zbuoy,zalpha,zdrag
	78	real zcor,zdelta,zcvm5,qlbef
	79	real Tbef,qsatbef
	80	real dqsat_dT,DT,num,denom
	81	REAL REPS,RLvCp,DDT0
	82	PARAMETER (DDT0=.01)
	83	logical Zsat
	84	LOGICAL active(ngrid),activetmp(ngrid)
	85	REAL fact_gamma,fact_epsilon,fact_gamma2
	86	REAL c2(ngrid,klev)
	87	REAL a1,m
	88
	89	REAL zw2fact,expa
	90	Zsat=.false.
	91	! Initialisation
	92	RLvCp = RLVTT/RCPD
	93
	94
	95	fact_epsilon=0.002
	96	a1=2./3.
	97	fact_gamma=0.9
	98	zfact=fact_gamma/(1+fact_gamma)
	99	fact_gamma2=zfact
	100	expa=0.
	101
	102
	103	! Initialisations des variables reeles
	104	if (1==1) then
	105	ztva(:,:)=ztv(:,:)
	106	ztva_est(:,:)=ztva(:,:)
	107	ztla(:,:)=zthl(:,:)
	108	zqta(:,:)=po(:,:)
	109	zha(:,:) = ztva(:,:)
	110	else
	111	ztva(:,:)=0.
	112	ztva_est(:,:)=0.
	113	ztla(:,:)=0.
	114	zqta(:,:)=0.
	115	zha(:,:) =0.
	116	endif
	117
	118	zqla_est(:,:)=0.
	119	zqsatth(:,:)=0.
	120	zqla(:,:)=0.
	121	detr_star(:,:)=0.
	122	entr_star(:,:)=0.
	123	alim_star(:,:)=0.
	124	alim_star_tot(:)=0.
	125	csc(:,:)=0.
	126	detr(:,:)=0.
	127	entr(:,:)=0.
	128	zw2(:,:)=0.
	129	w_est(:,:)=0.
	130	f_star(:,:)=0.
	131	wa_moy(:,:)=0.
	132	linter(:)=1.
	133	linter(:)=1.
	134
	135	! Initialisation des variables entieres
	136	lmix(:)=1
	137	lmix_bis(:)=2
	138	wmaxa(:)=0.
	139	lalim(:)=1
	140
	141	!-------------------------------------------------------------------------
	142	! On ne considere comme actif que les colonnes dont les deux premieres
	143	! couches sont instables.
	144	!-------------------------------------------------------------------------
	145	active(:)=ztv(:,1)>ztv(:,2)
	146
	147	!-------------------------------------------------------------------------
	148	! Definition de l'alimentation a l'origine dans thermcell_init
	149	!-------------------------------------------------------------------------
	150	do l=1,klev-1
	151	do ig=1,ngrid
	152	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
	153	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
	154	& *sqrt(zlev(ig,l+1))
	155	lalim(:)=l+1
	156	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	157	endif
	158	enddo
	159	enddo
	160	do l=1,klev
	161	do ig=1,ngrid
	162	if (alim_star_tot(ig) > 1.e-10 ) then
	163	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
	164	endif
	165	enddo
	166	enddo
	167	alim_star_tot(:)=1.
	168
	169
	170	!------------------------------------------------------------------------------
	171	! Calcul dans la premiere couche
	172	! On decide dans cette version que le thermique n'est actif que si la premiere
	173	! couche est instable.
	174	! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
	175	! dans une couche l>1
	176	!------------------------------------------------------------------------------
	177	do ig=1,ngrid
	178	! Le panache va prendre au debut les caracteristiques de l'air contenu
	179	! dans cette couche.
	180	if (active(ig)) then
	181	ztla(ig,1)=zthl(ig,1)
	182	zqta(ig,1)=po(ig,1)
	183	zqla(ig,1)=zl(ig,1)
	184	!cr: attention, prise en compte de f*(1)=1
	185	f_star(ig,2)=alim_star(ig,1)
	186	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
	187	& *(zlev(ig,2)-zlev(ig,1)) &
	188	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
	189	w_est(ig,2)=zw2(ig,2)
	190	endif
	191	enddo
	192	!
	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	if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
	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))
	354
	355	!on ecrit zqsat
	356	zqsatth(ig,l)=qsatbef
	357
	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)
	368
	369
	370	zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
	371
	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
	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
	381
	382
	383	endif
	384
	385	endif
	386	enddo
	387
	388	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
	389	!
	390	!---------------------------------------------------------------------------
	391	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
	392	!---------------------------------------------------------------------------
	393
	394	do ig=1,ngrid
	395	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	396	! stop'On tombe sur le cas particulier de thermcell_dry'
	397	print*,'On tombe sur le cas particulier de thermcell_plume'
	398	zw2(ig,l+1)=0.
	399	linter(ig)=l+1
	400	endif
	401
	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
	407
	408	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	409
	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
	420
	421	!=========================================================================
	422	! FIN DE LA BOUCLE VERTICALE
	423	enddo
	424	!=========================================================================
	425
	426	!on recalcule alim_star_tot
	427	do ig=1,ngrid
	428	alim_star_tot(ig)=0.
	429	enddo
	430	do ig=1,ngrid
	431	do l=1,lalim(ig)-1
	432	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	433	enddo
	434	enddo
	435
	436
	437	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
	438
	439
	440	return
	441	end
	442
	443	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	444	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	445	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	446	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	447	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	448	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	449	SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
	450	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	451	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	452	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
	453	& ,lev_out,lunout1,igout)
	454
	455	!--------------------------------------------------------------------------
	456	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
	457	! Version conforme a l'article de Rio et al. 2010.
	458	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
	459	!--------------------------------------------------------------------------
	460
	461	IMPLICIT NONE
	462
	463	#include "YOMCST.h"
	464	#include "YOETHF.h"
	465	#include "FCTTRE.h"
	466	#include "iniprint.h"
	467	#include "thermcell.h"
	468
	469	INTEGER itap
	470	INTEGER lunout1,igout
	471	INTEGER ngrid,klev
	472	REAL ptimestep
	473	REAL ztv(ngrid,klev)
	474	REAL zthl(ngrid,klev)
	475	REAL po(ngrid,klev)
	476	REAL zl(ngrid,klev)
	477	REAL rhobarz(ngrid,klev)
	478	REAL zlev(ngrid,klev+1)
	479	REAL pplev(ngrid,klev+1)
	480	REAL pphi(ngrid,klev)
	481	REAL zpspsk(ngrid,klev)
	482	REAL alim_star(ngrid,klev)
	483	REAL f0(ngrid)
	484	INTEGER lalim(ngrid)
	485	integer lev_out ! niveau pour les print
	486
	487	real alim_star_tot(ngrid)
	488
	489	REAL ztva(ngrid,klev)
	490	REAL ztla(ngrid,klev)
	491	REAL zqla(ngrid,klev)
	492	REAL zqta(ngrid,klev)
	493	REAL zha(ngrid,klev)
	494
	495	REAL detr_star(ngrid,klev)
	496	REAL coefc
	497	REAL entr_star(ngrid,klev)
	498	REAL detr(ngrid,klev)
	499	REAL entr(ngrid,klev)
	500
	501	REAL csc(ngrid,klev)
	502
	503	REAL zw2(ngrid,klev+1)
	504	REAL w_est(ngrid,klev+1)
	505	REAL f_star(ngrid,klev+1)
	506	REAL wa_moy(ngrid,klev+1)
	507
	508	REAL ztva_est(ngrid,klev)
	509	REAL zqla_est(ngrid,klev)
	510	REAL zqsatth(ngrid,klev)
	511	REAL zta_est(ngrid,klev)
	512	REAL ztemp(ngrid),zqsat(ngrid)
	513	REAL zdw2
	514	REAL zw2modif
	515	REAL zw2fact
	516	REAL zeps(ngrid,klev)
	517
	518	REAL linter(ngrid)
	519	INTEGER lmix(ngrid)
	520	INTEGER lmix_bis(ngrid)
	521	REAL wmaxa(ngrid)
	522
	523	INTEGER ig,l,k
	524
	525	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
	526	real zbuoybis
	527	real zcor,zdelta,zcvm5,qlbef,zdz2
	528	real betalpha,zbetalpha
	529	real eps, afact
	530	REAL REPS,RLvCp,DDT0
	531	PARAMETER (DDT0=.01)
	532	logical Zsat
	533	LOGICAL active(ngrid),activetmp(ngrid)
	534	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
	535	REAL c2(ngrid,klev)
	536	Zsat=.false.
	537	! Initialisation
	538
	539	RLvCp = RLVTT/RCPD
	540	fact_epsilon=0.002
	541	betalpha=0.9
	542	afact=2./3.
	543
	544	zbetalpha=betalpha/(1.+betalpha)
	545
	546
	547	! Initialisations des variables reeles
	548	if (1==0) then
	549	ztva(:,:)=ztv(:,:)
	550	ztva_est(:,:)=ztva(:,:)
	551	ztla(:,:)=zthl(:,:)
	552	zqta(:,:)=po(:,:)
	553	zha(:,:) = ztva(:,:)
	554	else
	555	ztva(:,:)=0.
	556	ztva_est(:,:)=0.
	557	ztla(:,:)=0.
	558	zqta(:,:)=0.
	559	zha(:,:) =0.
	560	endif
	561
	562	zqla_est(:,:)=0.
	563	zqsatth(:,:)=0.
	564	zqla(:,:)=0.
	565	detr_star(:,:)=0.
	566	entr_star(:,:)=0.
	567	alim_star(:,:)=0.
	568	alim_star_tot(:)=0.
	569	csc(:,:)=0.
	570	detr(:,:)=0.
	571	entr(:,:)=0.
	572	zw2(:,:)=0.
	573	zbuoy(:,:)=0.
	574	gamma(:,:)=0.
	575	zeps(:,:)=0.
	576	w_est(:,:)=0.
	577	f_star(:,:)=0.
	578	wa_moy(:,:)=0.
	579	linter(:)=1.
	580	! linter(:)=1.
	581	! Initialisation des variables entieres
	582	lmix(:)=1
	583	lmix_bis(:)=2
	584	wmaxa(:)=0.
	585	lalim(:)=1
	586
	587
	588	!-------------------------------------------------------------------------
	589	! On ne considere comme actif que les colonnes dont les deux premieres
	590	! couches sont instables.
	591	!-------------------------------------------------------------------------
	592	active(:)=ztv(:,1)>ztv(:,2)
	593
	594	!-------------------------------------------------------------------------
	595	! Definition de l'alimentation a l'origine dans thermcell_init
	596	!-------------------------------------------------------------------------
	597	do l=1,klev-1
	598	do ig=1,ngrid
	599	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
	600	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
	601	& *sqrt(zlev(ig,l+1))
	602	lalim(ig)=l+1
	603	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	604	endif
	605	enddo
	606	enddo
	607	do l=1,klev
	608	do ig=1,ngrid
	609	if (alim_star_tot(ig) > 1.e-10 ) then
	610	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
	611	endif
	612	enddo
	613	enddo
	614	alim_star_tot(:)=1.
	615
	616
	617
	618	!------------------------------------------------------------------------------
	619	! Calcul dans la premiere couche
	620	! On decide dans cette version que le thermique n'est actif que si la premiere
	621	! couche est instable.
	622	! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
	623	! dans une couche l>1
	624	!------------------------------------------------------------------------------
	625	do ig=1,ngrid
	626	! Le panache va prendre au debut les caracteristiques de l'air contenu
	627	! dans cette couche.
	628	if (active(ig)) then
	629	ztla(ig,1)=zthl(ig,1)
	630	zqta(ig,1)=po(ig,1)
	631	zqla(ig,1)=zl(ig,1)
	632	!cr: attention, prise en compte de f*(1)=1
	633	f_star(ig,2)=alim_star(ig,1)
	634	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
	635	& *(zlev(ig,2)-zlev(ig,1)) &
	636	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
	637	w_est(ig,2)=zw2(ig,2)
	638	endif
	639	enddo
	640	!
	641
	642	!==============================================================================
	643	!boucle de calcul de la vitesse verticale dans le thermique
	644	!==============================================================================
	645	do l=2,klev-1
	646	!==============================================================================
	647
	648
	649	! On decide si le thermique est encore actif ou non
	650	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
	651	do ig=1,ngrid
	652	active(ig)=active(ig) &
	653	& .and. zw2(ig,l)>1.e-10 &
	654	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
	655	enddo
	656
	657
	658
	659	!---------------------------------------------------------------------------
	660	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
	661	! sans tenir compte du detrainement et de l'entrainement dans cette
	662	! couche
	663	! C'est a dire qu'on suppose
	664	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
	665	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
	666	! avant) a l'alimentation pour avoir un calcul plus propre
	667	!---------------------------------------------------------------------------
	668
	669	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
	670	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
	671
	672	do ig=1,ngrid
	673	! print*,'active',active(ig),ig,l
	674	if(active(ig)) then
	675	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
	676	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
	677	zta_est(ig,l)=ztva_est(ig,l)
	678	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
	679	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
	680	& -zqla_est(ig,l))-zqla_est(ig,l))
	681
	682	!------------------------------------------------
	683	!AJAM:nouveau calcul de w²
	684	!------------------------------------------------
	685	zdz=zlev(ig,l+1)-zlev(ig,l)
	686	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	687
	688	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	689	zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
	690	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
	691
	692
	693	if (w_est(ig,l+1).lt.0.) then
	694	w_est(ig,l+1)=zw2(ig,l)
	695	endif
	696	endif
	697	enddo
	698
	699
	700	!-------------------------------------------------
	701	!calcul des taux d'entrainement et de detrainement
	702	!-------------------------------------------------
	703
	704	do ig=1,ngrid
	705	if (active(ig)) then
	706
	707	zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
	708	zw2m=w_est(ig,l+1)
	709	zdz=zlev(ig,l+1)-zlev(ig,l)
	710	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	711	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
	712	zbuoybis=zbuoy(ig,l)
	713	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
	714	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
	715
	716
	717	entr_star(ig,l)=f_star(ig,l)zdz zbetalpha*MAX(0., &
	718	& afact*zbuoybis/zw2m - fact_epsilon )
	719
	720
	721	detr_star(ig,l)=f_star(ig,l)*zdz &
	722	& MAX(1.e-3, -afactzbetalpha*zbuoy(ig,l)/zw2m &
	723	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
	724
	725	! En dessous de lalim, on prend le max de alim_star et entr_star pour
	726	! alim_star et 0 sinon
	727	if (l.lt.lalim(ig)) then
	728	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
	729	entr_star(ig,l)=0.
	730	endif
	731
	732	! Calcul du flux montant normalise
	733	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	734	& -detr_star(ig,l)
	735
	736	endif
	737	enddo
	738
	739
	740	!----------------------------------------------------------------------------
	741	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	742	!---------------------------------------------------------------------------
	743	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
	744	do ig=1,ngrid
	745	if (activetmp(ig)) then
	746	Zsat=.false.
	747	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
	748	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	749	& /(f_star(ig,l+1)+detr_star(ig,l))
	750	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
	751	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	752	& /(f_star(ig,l+1)+detr_star(ig,l))
	753
	754	endif
	755	enddo
	756
	757	ztemp(:)=zpspsk(:,l)*ztla(:,l)
	758	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
	759
	760	do ig=1,ngrid
	761	if (activetmp(ig)) then
	762	if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
	763	! on ecrit de maniere conservative (sat ou non)
	764	! T = Tl +Lv/Cp ql
	765	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
	766	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	767	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	768	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	769	zha(ig,l) = ztva(ig,l)
	770	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	771	& -zqla(ig,l))-zqla(ig,l))
	772	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	773	zdz=zlev(ig,l+1)-zlev(ig,l)
	774	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
	775
	776	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	777	zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
	778	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
	779	endif
	780	enddo
	781
	782	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
	783	!
	784	!---------------------------------------------------------------------------
	785	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
	786	!---------------------------------------------------------------------------
	787
	788	do ig=1,ngrid
	789	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	790	! stop'On tombe sur le cas particulier de thermcell_dry'
	791	print*,'On tombe sur le cas particulier de thermcell_plume'
	792	zw2(ig,l+1)=0.
	793	linter(ig)=l+1
	794	endif
	795
	796	if (zw2(ig,l+1).lt.0.) then
	797	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	798	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	799	zw2(ig,l+1)=0.
	800	endif
	801
	802	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	803
	804	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	805	! lmix est le niveau de la couche ou w (wa_moy) est maximum
	806	!on rajoute le calcul de lmix_bis
	807	if (zqla(ig,l).lt.1.e-10) then
	808	lmix_bis(ig)=l+1
	809	endif
	810	lmix(ig)=l+1
	811	wmaxa(ig)=wa_moy(ig,l+1)
	812	endif
	813	enddo
	814
	815	!=========================================================================
	816	! FIN DE LA BOUCLE VERTICALE
	817	enddo
	818	!=========================================================================
	819
	820	!on recalcule alim_star_tot
	821	do ig=1,ngrid
	822	alim_star_tot(ig)=0.
	823	enddo
	824	do ig=1,ngrid
	825	do l=1,lalim(ig)-1
	826	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	827	enddo
	828	enddo
	829
	830
	831	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
	832
	833	#undef wrgrads_thermcell
	834	#ifdef wrgrads_thermcell
	835	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
	836	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
	837	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
	838	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
	839	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
	840	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
	841	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
	842	#endif
	843
	844
	845	return
	846	end
	847

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

Download in other formats:

Original Format