Context Navigation

thermcell_plume.F90 @ 1227

Last change on this file since 1227 was 1057, checked in by lmdzadmin, 16 years ago
Un peu de nettoyage ACA/FH/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 26.8 KB

Rev	Line
[972]	1	SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
[1026]	2	& zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot, &
[878]	3	& lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, &
[1026]	4	& ztla,zqla,zqta,zha,zw2,w_est,zqsatth,lmix,lmix_bis,linter &
	5	& ,lev_out,lunout1,igout)
[878]	6
	7	!--------------------------------------------------------------------------
	8	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
	9	!--------------------------------------------------------------------------
	10
	11	IMPLICIT NONE
	12
	13	#include "YOMCST.h"
	14	#include "YOETHF.h"
	15	#include "FCTTRE.h"
[938]	16	#include "iniprint.h"
[1026]	17	#include "thermcell.h"
[878]	18
[972]	19	INTEGER itap
[1026]	20	INTEGER lunout1,igout
[878]	21	INTEGER ngrid,klev
[972]	22	REAL ptimestep
[878]	23	REAL ztv(ngrid,klev)
	24	REAL zthl(ngrid,klev)
	25	REAL po(ngrid,klev)
	26	REAL zl(ngrid,klev)
	27	REAL rhobarz(ngrid,klev)
	28	REAL zlev(ngrid,klev+1)
	29	REAL pplev(ngrid,klev+1)
	30	REAL pphi(ngrid,klev)
	31	REAL zpspsk(ngrid,klev)
	32	REAL alim_star(ngrid,klev)
	33	REAL zmax_sec(ngrid)
	34	REAL f0(ngrid)
	35	REAL l_mix
	36	REAL r_aspect
	37	INTEGER lalim(ngrid)
	38	integer lev_out ! niveau pour les print
[972]	39	real zcon2(ngrid)
[1026]	40
	41	real alim_star_tot(ngrid)
[878]	42
	43	REAL ztva(ngrid,klev)
	44	REAL ztla(ngrid,klev)
	45	REAL zqla(ngrid,klev)
[972]	46	REAL zqla0(ngrid,klev)
[878]	47	REAL zqta(ngrid,klev)
	48	REAL zha(ngrid,klev)
	49
	50	REAL detr_star(ngrid,klev)
[972]	51	REAL coefc
	52	REAL detr_stara(ngrid,klev)
	53	REAL detr_starb(ngrid,klev)
	54	REAL detr_starc(ngrid,klev)
	55	REAL detr_star0(ngrid,klev)
	56	REAL detr_star1(ngrid,klev)
	57	REAL detr_star2(ngrid,klev)
	58
[878]	59	REAL entr_star(ngrid,klev)
[972]	60	REAL entr_star1(ngrid,klev)
	61	REAL entr_star2(ngrid,klev)
[878]	62	REAL detr(ngrid,klev)
	63	REAL entr(ngrid,klev)
	64
	65	REAL zw2(ngrid,klev+1)
	66	REAL w_est(ngrid,klev+1)
	67	REAL f_star(ngrid,klev+1)
	68	REAL wa_moy(ngrid,klev+1)
	69
	70	REAL ztva_est(ngrid,klev)
	71	REAL zqla_est(ngrid,klev)
	72	REAL zqsatth(ngrid,klev)
[1026]	73	REAL zta_est(ngrid,klev)
[878]	74
	75	REAL linter(ngrid)
	76	INTEGER lmix(ngrid)
[972]	77	INTEGER lmix_bis(ngrid)
[878]	78	REAL wmaxa(ngrid)
	79
	80	INTEGER ig,l,k
	81
	82	real zcor,zdelta,zcvm5,qlbef
	83	real Tbef,qsatbef
	84	real dqsat_dT,DT,num,denom
	85	REAL REPS,RLvCp,DDT0
	86	PARAMETER (DDT0=.01)
	87	logical Zsat
[1026]	88	REAL fact_gamma,fact_epsilon
	89	REAL c2(ngrid,klev)
[878]	90
	91	Zsat=.false.
	92	! Initialisation
	93	RLvCp = RLVTT/RCPD
	94
[1026]	95	if (iflag_thermals_ed==0) then
	96	fact_gamma=1.
	97	fact_epsilon=1.
	98	else if (iflag_thermals_ed==1) then
	99	fact_gamma=1.
	100	fact_epsilon=1.
	101	else if (iflag_thermals_ed==2) then
	102	fact_gamma=1.
	103	fact_epsilon=2.
	104	endif
	105
[878]	106	do l=1,klev
	107	do ig=1,ngrid
	108	zqla_est(ig,l)=0.
	109	ztva_est(ig,l)=ztva(ig,l)
	110	zqsatth(ig,l)=0.
	111	enddo
	112	enddo
	113
[1026]	114	!CR: attention test couche alim
	115	! do l=2,klev
	116	! do ig=1,ngrid
	117	! alim_star(ig,l)=0.
	118	! enddo
	119	! enddo
[878]	120	!AM:initialisations du thermique
	121	do k=1,klev
	122	do ig=1,ngrid
	123	ztva(ig,k)=ztv(ig,k)
	124	ztla(ig,k)=zthl(ig,k)
	125	zqla(ig,k)=0.
	126	zqta(ig,k)=po(ig,k)
[972]	127	!
	128	ztva(ig,k) = ztla(ig,k)zpspsk(ig,k)+RLvCpzqla(ig,k)
	129	ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k)
	130	zha(ig,k) = ztva(ig,k)
	131	!
[878]	132	enddo
	133	enddo
	134	do k=1,klev
	135	do ig=1,ngrid
	136	detr_star(ig,k)=0.
	137	entr_star(ig,k)=0.
[972]	138
	139	detr_stara(ig,k)=0.
	140	detr_starb(ig,k)=0.
	141	detr_starc(ig,k)=0.
	142	detr_star0(ig,k)=0.
	143	zqla0(ig,k)=0.
	144	detr_star1(ig,k)=0.
	145	detr_star2(ig,k)=0.
	146	entr_star1(ig,k)=0.
	147	entr_star2(ig,k)=0.
	148
[878]	149	detr(ig,k)=0.
	150	entr(ig,k)=0.
	151	enddo
	152	enddo
[938]	153	if (prt_level.ge.1) print*,'7 OK convect8'
[878]	154	do k=1,klev+1
	155	do ig=1,ngrid
	156	zw2(ig,k)=0.
	157	w_est(ig,k)=0.
	158	f_star(ig,k)=0.
	159	wa_moy(ig,k)=0.
	160	enddo
	161	enddo
	162
[938]	163	if (prt_level.ge.1) print*,'8 OK convect8'
[878]	164	do ig=1,ngrid
	165	linter(ig)=1.
	166	lmix(ig)=1
[1026]	167	lmix_bis(ig)=2
[878]	168	wmaxa(ig)=0.
	169	enddo
	170
	171	!-----------------------------------------------------------------------------------
	172	!boucle de calcul de la vitesse verticale dans le thermique
	173	!-----------------------------------------------------------------------------------
	174	do l=1,klev-1
	175	do ig=1,ngrid
[972]	176
	177
	178
	179	! Calcul dans la premiere couche active du thermique (ce qu'on teste
	180	! en disant que la couche est instable et que w2 en bas de la couche
	181	! est nulle.
	182
[878]	183	if (ztv(ig,l).gt.ztv(ig,l+1) &
	184	& .and.alim_star(ig,l).gt.1.e-10 &
	185	& .and.zw2(ig,l).lt.1e-10) then
[972]	186
	187
	188	! Le panache va prendre au debut les caracteristiques de l'air contenu
	189	! dans cette couche.
[878]	190	ztla(ig,l)=zthl(ig,l)
	191	zqta(ig,l)=po(ig,l)
	192	zqla(ig,l)=zl(ig,l)
	193	f_star(ig,l+1)=alim_star(ig,l)
[972]	194
[878]	195	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1) &
	196	& *(zlev(ig,l+1)-zlev(ig,l)) &
	197	& 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
	198	w_est(ig,l+1)=zw2(ig,l+1)
	199	!
[972]	200
	201
[878]	202	else if ((zw2(ig,l).ge.1e-10).and. &
	203	& (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
	204	!estimation du detrainement a partir de la geometrie du pas precedent
	205	!tests sur la definition du detr
	206	!calcul de detr_star et entr_star
	207
	208
	209
	210	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	211	! FH le test miraculeux de Catherine ? Le bout du tunel ?
	212	! w_est(ig,3)=zw2(ig,2)* &
	213	! & ((f_star(ig,2))**2) &
	214	! & /(f_star(ig,2)+alim_star(ig,2))**2+ &
	215	! & 2.RG(ztva(ig,1)-ztv(ig,2))/ztv(ig,2) &
	216	! & *(zlev(ig,3)-zlev(ig,2))
	217	! if (l.gt.2) then
	218	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	219
[972]	220
	221
	222	! Premier calcul de la vitesse verticale a partir de la temperature
	223	! potentielle virtuelle
	224
	225	! FH CESTQUOI CA ????
	226	#define int1d2
	227	!#undef int1d2
	228	#ifdef int1d2
	229	if (l.ge.2) then
	230	#else
[878]	231	if (l.gt.2) then
[972]	232	#endif
	233
	234	if (1.eq.1) then
[878]	235	w_est(ig,3)=zw2(ig,2)* &
	236	& ((f_star(ig,2))**2) &
	237	& /(f_star(ig,2)+alim_star(ig,2))**2+ &
	238	& 2.RG(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
[1026]	239	! & *1./3. &
[878]	240	& *(zlev(ig,3)-zlev(ig,2))
[972]	241	endif
[878]	242
	243
	244	!---------------------------------------------------------------------------
	245	!calcul de l entrainement et du detrainement lateral
	246	!---------------------------------------------------------------------------
	247	!
	248	!test:estimation de ztva_new_est sans entrainement
[972]	249
[878]	250	Tbef=ztla(ig,l-1)*zpspsk(ig,l)
	251	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	252	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
	253	qsatbef=MIN(0.5,qsatbef)
	254	zcor=1./(1.-retv*qsatbef)
	255	qsatbef=qsatbef*zcor
	256	Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
	257	if (Zsat) then
	258	qlbef=max(0.,zqta(ig,l-1)-qsatbef)
	259	DT = 0.5RLvCpqlbef
	260	do while (abs(DT).gt.DDT0)
	261	Tbef=Tbef+DT
	262	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	263	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
	264	qsatbef=MIN(0.5,qsatbef)
	265	zcor=1./(1.-retv*qsatbef)
	266	qsatbef=qsatbef*zcor
	267	qlbef=zqta(ig,l-1)-qsatbef
	268
	269	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	270	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
	271	zcor=1./(1.-retv*qsatbef)
	272	dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
	273	num=-Tbef+ztla(ig,l-1)zpspsk(ig,l)+RLvCpqlbef
	274	denom=1.+RLvCp*dqsat_dT
	275	DT=num/denom
	276	enddo
	277	zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef)
	278	endif
	279	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
	280	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
[1026]	281	zta_est(ig,l)=ztva_est(ig,l)
[878]	282	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
	283	& -zqla_est(ig,l))-zqla_est(ig,l))
	284
	285	w_est(ig,l+1)=zw2(ig,l)* &
	286	& ((f_star(ig,l))**2) &
	287	& /(f_star(ig,l)+alim_star(ig,l))**2+ &
	288	& 2.RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
[1026]	289	! & *1./3. &
[878]	290	& *(zlev(ig,l+1)-zlev(ig,l))
	291	if (w_est(ig,l+1).lt.0.) then
	292	w_est(ig,l+1)=zw2(ig,l)
	293	endif
	294	!
	295	!calcul du detrainement
[972]	296	!=======================
	297
[1026]	298	!CR:on vire les modifs
	299	if (iflag_thermals_ed==0) then
[972]	300
	301	! Modifications du calcul du detrainement.
	302	! Dans la version de la these de Catherine, on passe brusquement
	303	! de la version seche a la version nuageuse pour le detrainement
	304	! ce qui peut occasioner des oscillations.
	305	! dans la nouvelle version, on commence par calculer un detrainement sec.
	306	! Puis un autre en cas de nuages.
	307	! Puis on combine les deux lineairement en fonction de la quantite d'eau.
	308
	309	#define int1d3
	310	!#undef int1d3
	311	#define RIO_TH
	312	#ifdef RIO_TH
	313	!1. Cas non nuageux
	314	! 1.1 on est sous le zmax_sec et w croit
[878]	315	if ((w_est(ig,l+1).gt.w_est(ig,l)).and. &
	316	& (zlev(ig,l+1).lt.zmax_sec(ig)).and. &
[972]	317	#ifdef int1d3
	318	& (zqla_est(ig,l).lt.1.e-10)) then
	319	#else
[878]	320	& (zqla(ig,l-1).lt.1.e-10)) then
[972]	321	#endif
[878]	322	detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1) &
	323	& sqrt(w_est(ig,l+1))sqrt(l_mix*zlev(ig,l+1)) &
	324	& -rhobarz(ig,l)sqrt(w_est(ig,l))sqrt(l_mix*zlev(ig,l))) &
	325	& /(r_aspect*zmax_sec(ig)))
[972]	326	detr_stara(ig,l)=detr_star(ig,l)
	327
	328	if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l
	329
	330	! 1.2 on est sous le zmax_sec et w decroit
[878]	331	else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and. &
[972]	332	#ifdef int1d3
	333	& (zqla_est(ig,l).lt.1.e-10)) then
	334	#else
[878]	335	& (zqla(ig,l-1).lt.1.e-10)) then
[972]	336	#endif
[878]	337	detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig)) &
	338	& /(rhobarz(ig,lmix(ig))wmaxa(ig)) &
	339	& (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1)) &
	340	& *((zmax_sec(ig)-zlev(ig,l+1))/ &
	341	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2. &
	342	& -rhobarz(ig,l)*sqrt(w_est(ig,l)) &
	343	& *((zmax_sec(ig)-zlev(ig,l))/ &
	344	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
[972]	345	detr_starb(ig,l)=detr_star(ig,l)
	346
	347	if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l
	348
[878]	349	else
[972]	350
	351	! 1.3 dans les autres cas
[878]	352	detr_star(ig,l)=0.002f0(ig)f_star(ig,l) &
	353	& *(zlev(ig,l+1)-zlev(ig,l))
[972]	354	detr_starc(ig,l)=detr_star(ig,l)
	355
	356	if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
[878]	357
	358	endif
[972]	359
	360	#else
	361
	362	! 1.1 on est sous le zmax_sec et w croit
	363	if ((w_est(ig,l+1).gt.w_est(ig,l)).and. &
	364	& (zlev(ig,l+1).lt.zmax_sec(ig)) ) then
	365	detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1) &
	366	& sqrt(w_est(ig,l+1))sqrt(l_mix*zlev(ig,l+1)) &
	367	& -rhobarz(ig,l)sqrt(w_est(ig,l))sqrt(l_mix*zlev(ig,l))) &
	368	& /(r_aspect*zmax_sec(ig)))
	369
	370	if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
	371
	372	! 1.2 on est sous le zmax_sec et w decroit
	373	else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then
	374	detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig)) &
	375	& /(rhobarz(ig,lmix(ig))wmaxa(ig)) &
	376	& (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1)) &
	377	& *((zmax_sec(ig)-zlev(ig,l+1))/ &
	378	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2. &
	379	& -rhobarz(ig,l)*sqrt(w_est(ig,l)) &
	380	& *((zmax_sec(ig)-zlev(ig,l))/ &
	381	& ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
	382	if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
	383
	384	else
	385	detr_star=0.
	386	endif
	387
	388	! 1.3 dans les autres cas
	389	detr_starc(ig,l)=0.002f0(ig)f_star(ig,l) &
	390	& *(zlev(ig,l+1)-zlev(ig,l))
	391
	392	coefc=min(zqla(ig,l-1)/1.e-3,1.)
	393	if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1.
	394	coefc=1.
	395	! il semble qu'il soit important de baser le calcul sur
	396	! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l)
	397	detr_star(ig,l)=detr_starc(ig,l)coefc+detr_star(ig,l)(1.-coefc)
	398
	399	if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l
	400
	401	#endif
	402
	403
	404	if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l
	405	!IM 730508 beg
	406	! if(itap.GE.7200) THEN
	407	! print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l)
	408	! endif
	409	!IM 730508 end
	410
	411	zqla0(ig,l)=zqla_est(ig,l)
	412	detr_star0(ig,l)=detr_star(ig,l)
	413	!IM 060508 beg
	414	! if(detr_star(ig,l).GT.1.) THEN
	415	! print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) &
	416	! & ,detr_starc(ig,l),coefc
	417	! endif
	418	!IM 060508 end
	419	!IM 160508 beg
	420	!IM 160508 IF (f0(ig).NE.0.) THEN
	421	detr_star(ig,l)=detr_star(ig,l)/f0(ig)
	422	!IM 160508 ELSE IF(detr_star(ig,l).EQ.0.) THEN
	423	!IM 160508 print*,'WARNING1 : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap
	424	!IM 160508 ELSE
	425	!IM 160508 print*,'WARNING2 : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l)
	426	!IM 160508 ENDIF
	427	!IM 160508 end
	428	!IM 060508 beg
	429	! if(detr_star(ig,l).GT.1.) THEN
	430	! print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), &
	431	! & float(1)/f0(ig)
	432	! endif
	433	!IM 060508 end
	434	if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l
[878]	435	!
	436	!calcul de entr_star
[972]	437
	438	! #undef test2
	439	! #ifdef test2
	440	! La version test2 destabilise beaucoup le modele.
	441	! Il semble donc que ca aide d'avoir un entrainement important sous
	442	! le nuage.
	443	! if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then
	444	! entr_star(ig,l)=0.4*detr_star(ig,l)
	445	! else
	446	! entr_star(ig,l)=0.
	447	! endif
	448	! #else
[878]	449	!
	450	! Deplacement du calcul de entr_star pour eviter d'avoir aussi
	451	! entr_star > fstar.
	452	! Redeplacer suite a la transformation du cas detr>f
	453	! FH
[972]	454
	455	if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l
	456	#define int1d
	457	!FH 070508 #define int1d4
	458	!#undef int1d4
	459	! L'option int1d4 correspond au choix dans le cas ou le detrainement
	460	! devient trop grand.
	461
	462	#ifdef int1d
	463
	464	#ifdef int1d4
	465	#else
	466	detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l))
	467	!FH 070508 plus
	468	detr_star(ig,l)=min(detr_star(ig,l),1.)
	469	#endif
	470
	471	entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.)
	472
	473	if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l
	474	#ifdef int1d4
	475	! Si le detrainement excede le flux en bas + l'entrainement, le thermique
	476	! doit disparaitre.
	477	if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then
	478	detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l)
	479	f_star(ig,l+1)=0.
	480	linter(ig)=l+1
	481	zw2(ig,l+1)=-1.e-10
	482	endif
	483	#endif
	484
	485
	486	#else
	487
	488	if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l
[878]	489	if(l.gt.lalim(ig)) then
	490	entr_star(ig,l)=0.4*detr_star(ig,l)
	491	else
	492
	493	! FH :
	494	! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
	495	! en haut de la couche d'alimentation.
	496	! A remettre en questoin a la premiere occasion mais ca peut aider a
	497	! ecrire un code robuste.
	498	! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
	499	! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
	500	! d'alimentation, ce qui n'est pas forcement heureux.
[972]	501
	502	if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l
	503	#undef pre_int1c
	504	#ifdef pre_int1c
[878]	505	entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
	506	detr_star(ig,l)=entr_star(ig,l)
[972]	507	#else
	508	entr_star(ig,l)=0.
	509	#endif
	510
[878]	511	endif
	512
[972]	513	#endif
	514
	515	if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l
	516	entr_star1(ig,l)=entr_star(ig,l)
	517	detr_star1(ig,l)=detr_star(ig,l)
[878]	518	!
[972]	519
	520	#ifdef int1d
	521	#else
[878]	522	if (detr_star(ig,l).gt.f_star(ig,l)) then
	523
	524	! Ce test est là entre autres parce qu'on passe par des valeurs
	525	! delirantes de detr_star.
	526	! ca vaut sans doute le coup de verifier pourquoi.
	527
	528	detr_star(ig,l)=f_star(ig,l)
[972]	529	#ifdef pre_int1c
[878]	530	if (l.gt.lalim(ig)+1) then
	531	entr_star(ig,l)=0.
	532	alim_star(ig,l)=0.
	533	! FH ajout pour forcer a stoper le thermique juste sous le sommet
	534	! de la couche (voir calcul de finter)
	535	zw2(ig,l+1)=-1.e-10
	536	linter(ig)=l+1
	537	else
[972]	538	entr_star(ig,l)=0.4*detr_star(ig,l)
[878]	539	endif
[972]	540	#else
	541	entr_star(ig,l)=0.4*detr_star(ig,l)
	542	#endif
[878]	543	endif
[972]	544	#endif
[878]	545
[972]	546	else !l > 2
[878]	547	detr_star(ig,l)=0.
	548	entr_star(ig,l)=0.
	549	endif
	550
[972]	551	entr_star2(ig,l)=entr_star(ig,l)
	552	detr_star2(ig,l)=detr_star(ig,l)
	553	if (prt_level.ge.20) print*,'coucou calcul detr 450: ig, l', ig, l
	554
[1026]	555	endif ! iflag_thermals_ed==0
	556
	557	!CR:nvlle def de entr_star et detr_star
	558	if (iflag_thermals_ed>=1) then
	559	! if (l.lt.lalim(ig)) then
	560	! if (l.lt.2) then
	561	! entr_star(ig,l)=0.
	562	! detr_star(ig,l)=0.
	563	! else
	564	! if (0.001.gt.(RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))/(2.w_est(ig,l+1)))) then
	565	! entr_star(ig,l)=0.001f_star(ig,l)(zlev(ig,l+1)-zlev(ig,l))
	566	! else
	567	! entr_star(ig,l)= &
	568	! & f_star(ig,l)/(2.*w_est(ig,l+1)) &
	569	! & RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)) &
	570	! & *(zlev(ig,l+1)-zlev(ig,l))
	571
	572
	573	entr_star(ig,l)=MAX(0.f_star(ig,l)(zlev(ig,l+1)-zlev(ig,l)), &
	574	& f_star(ig,l)/(2.*w_est(ig,l+1)) &
	575	& RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
	576	& *(zlev(ig,l+1)-zlev(ig,l))) &
	577	& +0.0001f_star(ig,l)(zlev(ig,l+1)-zlev(ig,l))
	578
	579	if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
	580	alim_star_tot(ig)=alim_star_tot(ig)+entr_star(ig,l)
	581	lalim(ig)=lmix_bis(ig)
[1057]	582	if(prt_level.GE.10) print*,'alim_star_tot',alim_star_tot(ig),entr_star(ig,l)
[1026]	583	endif
	584
	585	if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
	586	! c2(ig,l)=2500000.zqla_est(ig,l)/(1004.zta_est(ig,l))
	587	c2(ig,l)=0.001
	588	detr_star(ig,l)=MAX(0.f_star(ig,l)(zlev(ig,l+1)-zlev(ig,l)), &
	589	& c2(ig,l)f_star(ig,l)(zlev(ig,l+1)-zlev(ig,l)) &
	590	& -f_star(ig,l)/(2.*w_est(ig,l+1)) &
	591	& RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
	592	& *(zlev(ig,l+1)-zlev(ig,l))) &
	593	& +0.0001f_star(ig,l)(zlev(ig,l+1)-zlev(ig,l))
	594
	595	else
	596	! c2(ig,l)=2500000.zqla_est(ig,l)/(1004.zta_est(ig,l))
	597	c2(ig,l)=0.003
	598
	599	detr_star(ig,l)=MAX(0.f_star(ig,l)(zlev(ig,l+1)-zlev(ig,l)), &
	600	& c2(ig,l)f_star(ig,l)(zlev(ig,l+1)-zlev(ig,l)) &
	601	& -f_star(ig,l)/(2.*w_est(ig,l+1)) &
	602	& RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) &
	603	& *(zlev(ig,l+1)-zlev(ig,l))) &
	604	& +0.0002f_star(ig,l)(zlev(ig,l+1)-zlev(ig,l))
	605	endif
	606
	607
	608	! detr_star(ig,l)=detr_star(ig,l)*3.
	609	! if (l.lt.lalim(ig)) then
	610	! entr_star(ig,l)=0.
	611	! endif
	612	! if (l.lt.2) then
	613	! entr_star(ig,l)=0.
	614	! detr_star(ig,l)=0.
	615	! endif
	616
	617
	618	! endif
	619	! else if ((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10) then
	620	! entr_star(ig,l)=MAX(0.,0.8f_star(ig,l)/(2.w_est(ig,l+1)) &
	621	! & RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)) &
	622	! & *(zlev(ig,l+1)-zlev(ig,l))
	623	! detr_star(ig,l)=0.002*f_star(ig,l) &
	624	! & *(zlev(ig,l+1)-zlev(ig,l))
	625	! else
	626	! entr_star(ig,l)=0.001*f_star(ig,l) &
	627	! & *(zlev(ig,l+1)-zlev(ig,l))
	628	! detr_star(ig,l)=MAX(0.,-0.2f_star(ig,l)/(2.w_est(ig,l+1)) &
	629	! & RG(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)) &
	630	! & *(zlev(ig,l+1)-zlev(ig,l)) &
	631	! & +0.002*f_star(ig,l) &
	632	! & *(zlev(ig,l+1)-zlev(ig,l))
	633	! endif
	634
	635	endif ! iflag_thermals_ed==1
	636
[878]	637	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	638	! FH inutile si on conserve comme on l'a fait plus haut entr=detr
	639	! dans la couche d'alimentation
	640	!pas d entrainement dans la couche alim
	641	! if ((l.le.lalim(ig))) then
	642	! entr_star(ig,l)=0.
	643	! endif
	644	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	645	!
	646	!prise en compte du detrainement et de l entrainement dans le calcul du flux
[1026]	647
[878]	648	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	649	& -detr_star(ig,l)
	650
	651	!test sur le signe de f_star
[972]	652	if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l
[878]	653	if (f_star(ig,l+1).gt.1.e-10) then
	654	!----------------------------------------------------------------------------
	655	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	656	!---------------------------------------------------------------------------
	657	!
	658	Zsat=.false.
	659	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
	660	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	661	& /(f_star(ig,l+1)+detr_star(ig,l))
	662	!
	663	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
	664	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	665	& /(f_star(ig,l+1)+detr_star(ig,l))
	666	!
	667	Tbef=ztla(ig,l)*zpspsk(ig,l)
	668	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	669	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
	670	qsatbef=MIN(0.5,qsatbef)
	671	zcor=1./(1.-retv*qsatbef)
	672	qsatbef=qsatbef*zcor
	673	Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
	674	if (Zsat) then
	675	qlbef=max(0.,zqta(ig,l)-qsatbef)
	676	DT = 0.5RLvCpqlbef
	677	do while (abs(DT).gt.DDT0)
	678	Tbef=Tbef+DT
	679	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	680	qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
	681	qsatbef=MIN(0.5,qsatbef)
	682	zcor=1./(1.-retv*qsatbef)
	683	qsatbef=qsatbef*zcor
	684	qlbef=zqta(ig,l)-qsatbef
	685
	686	zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
	687	zcvm5=R5LES(1.-zdelta) + R5IESzdelta
	688	zcor=1./(1.-retv*qsatbef)
	689	dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
	690	num=-Tbef+ztla(ig,l)zpspsk(ig,l)+RLvCpqlbef
	691	denom=1.+RLvCp*dqsat_dT
	692	DT=num/denom
	693	enddo
	694	zqla(ig,l) = max(0.,qlbef)
	695	endif
	696	!
[972]	697	if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
[878]	698	! on ecrit de maniere conservative (sat ou non)
	699	! T = Tl +Lv/Cp ql
	700	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	701	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	702	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	703	zha(ig,l) = ztva(ig,l)
	704	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	705	& -zqla(ig,l))-zqla(ig,l))
	706
	707	!on ecrit zqsat
	708	zqsatth(ig,l)=qsatbef
	709	!calcul de vitesse
	710	zw2(ig,l+1)=zw2(ig,l)* &
	711	& ((f_star(ig,l))**2) &
[1026]	712	! Tests de Catherine
	713	! & /(f_star(ig,l+1)+detr_star(ig,l))**2+ &
	714	& /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)(1.-fact_epsilon))*2+ &
[878]	715	& 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
[1026]	716	& *fact_gamma &
[878]	717	& *(zlev(ig,l+1)-zlev(ig,l))
[1026]	718	!prise en compte des forces de pression que qd flottabilité<0
	719	! zw2(ig,l+1)=zw2(ig,l)* &
	720	! & 1./(1.+2.*entr_star(ig,l)/f_star(ig,l)) + &
	721	! & (f_star(ig,l))**2 &
	722	! & /(f_star(ig,l)+entr_star(ig,l))**2+ &
	723	! & (f_star(ig,l)-2.entr_star(ig,l))2/(f_star(ig,l)+2.entr_star(ig,l))**2+ &
	724	! & /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)(1.-2.))*2+ &
	725	! & /(f_star(ig,l)*2+2.2.detr_star(ig,l)f_star(ig,l)+2.entr_star(ig,l)f_star(ig,l))+ &
	726	! & 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
	727	! & *1./3. &
	728	! & *(zlev(ig,l+1)-zlev(ig,l))
	729
	730	! write(30,*),l+1,zw2(ig,l+1)-zw2(ig,l), &
	731	! & -2.entr_star(ig,l)/f_star(ig,l)zw2(ig,l), &
	732	! & 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
	733
	734
	735	! zw2(ig,l+1)=zw2(ig,l)* &
	736	! & (2.-2.*entr_star(ig,l)/f_star(ig,l)) &
	737	! & -zw2(ig,l-1)+ &
	738	! & 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) &
	739	! & *1./3. &
	740	! & *(zlev(ig,l+1)-zlev(ig,l))
	741
[878]	742	endif
	743	endif
[972]	744	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
[878]	745	!
	746	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
	747
	748	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	749	! stop'On tombe sur le cas particulier de thermcell_dry'
	750	print*,'On tombe sur le cas particulier de thermcell_plume'
	751	zw2(ig,l+1)=0.
	752	linter(ig)=l+1
	753	endif
	754
[1026]	755	! if ((zw2(ig,l).gt.0.).and. (zw2(ig,l+1).le.0.)) then
[878]	756	if (zw2(ig,l+1).lt.0.) then
	757	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	758	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	759	zw2(ig,l+1)=0.
	760	endif
	761
	762	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	763
	764	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	765	! lmix est le niveau de la couche ou w (wa_moy) est maximum
[1026]	766	!on rajoute le calcul de lmix_bis
	767	if (zqla(ig,l).lt.1.e-10) then
	768	lmix_bis(ig)=l+1
	769	endif
[878]	770	lmix(ig)=l+1
	771	wmaxa(ig)=wa_moy(ig,l+1)
	772	endif
	773	enddo
	774	enddo
	775
[1026]	776	!on remplace a* par e* ds premiere couche
	777	! if (iflag_thermals_ed.ge.1) then
	778	! do ig=1,ngrid
	779	! do l=2,klev
	780	! if (l.lt.lalim(ig)) then
	781	! alim_star(ig,l)=entr_star(ig,l)
	782	! endif
	783	! enddo
	784	! enddo
	785	! do ig=1,ngrid
	786	! lalim(ig)=lmix_bis(ig)
	787	! enddo
	788	! endif
	789	if (iflag_thermals_ed.ge.1) then
	790	do ig=1,ngrid
	791	do l=2,lalim(ig)
	792	alim_star(ig,l)=entr_star(ig,l)
	793	entr_star(ig,l)=0.
	794	enddo
	795	enddo
	796	endif
[972]	797	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
[878]	798
[972]	799	! print*,'thermcell_plume OK'
	800
[878]	801	return
	802	end

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

Download in other formats:

Original Format