Context Navigation

thermcell_plume.F90 @ 2093

Last change on this file since 2093 was 2046, checked in by fhourdin, 11 years ago

Nouvelles options des thermiques en route vers les stratocumulus

+ sorties de la discretisation verticale dans le vieux disvert_old

New options in thermals for stratocumulus

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 36.0 KB

Rev	Line
[1403]	1	!
	2	! $Id: thermcell_plume.F90 2046 2014-05-18 13:18:19Z fhourdin $
	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 &
[2046]	8	& ,lev_out,lunout1,igout)
[878]	9	!--------------------------------------------------------------------------
[1998]	10	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
[878]	11	!--------------------------------------------------------------------------
	12
[1968]	13	IMPLICIT NONE
[878]	14
	15	#include "YOMCST.h"
	16	#include "YOETHF.h"
	17	#include "FCTTRE.h"
[938]	18	#include "iniprint.h"
[1026]	19	#include "thermcell.h"
[878]	20
[972]	21	INTEGER itap
[1026]	22	INTEGER lunout1,igout
[878]	23	INTEGER ngrid,klev
[972]	24	REAL ptimestep
[878]	25	REAL ztv(ngrid,klev)
	26	REAL zthl(ngrid,klev)
	27	REAL po(ngrid,klev)
	28	REAL zl(ngrid,klev)
	29	REAL rhobarz(ngrid,klev)
	30	REAL zlev(ngrid,klev+1)
	31	REAL pplev(ngrid,klev+1)
	32	REAL pphi(ngrid,klev)
	33	REAL zpspsk(ngrid,klev)
	34	REAL alim_star(ngrid,klev)
	35	REAL f0(ngrid)
	36	INTEGER lalim(ngrid)
	37	integer lev_out ! niveau pour les print
[1503]	38	integer nbpb
[1026]	39
	40	real alim_star_tot(ngrid)
[878]	41
	42	REAL ztva(ngrid,klev)
	43	REAL ztla(ngrid,klev)
	44	REAL zqla(ngrid,klev)
	45	REAL zqta(ngrid,klev)
	46	REAL zha(ngrid,klev)
	47
	48	REAL detr_star(ngrid,klev)
[972]	49	REAL coefc
[878]	50	REAL entr_star(ngrid,klev)
	51	REAL detr(ngrid,klev)
	52	REAL entr(ngrid,klev)
	53
[1403]	54	REAL csc(ngrid,klev)
	55
[878]	56	REAL zw2(ngrid,klev+1)
	57	REAL w_est(ngrid,klev+1)
	58	REAL f_star(ngrid,klev+1)
	59	REAL wa_moy(ngrid,klev+1)
	60
	61	REAL ztva_est(ngrid,klev)
[1968]	62	REAL ztv_est(ngrid,klev)
[878]	63	REAL zqla_est(ngrid,klev)
	64	REAL zqsatth(ngrid,klev)
[1026]	65	REAL zta_est(ngrid,klev)
[1968]	66	REAL ztemp(ngrid),zqsat(ngrid)
	67	REAL zdw2,zdw2bis
[1403]	68	REAL zw2modif
[1968]	69	REAL zw2fact,zw2factbis
	70	REAL zeps(ngrid,klev)
[878]	71
	72	REAL linter(ngrid)
	73	INTEGER lmix(ngrid)
[972]	74	INTEGER lmix_bis(ngrid)
[878]	75	REAL wmaxa(ngrid)
	76
[1968]	77	INTEGER ig,l,k,lt,it
[878]	78
[1968]	79	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
	80	real zbuoyjam(ngrid,klev)
	81	real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
[2046]	82	real fact_shell
	83	real ztv1,ztv2,factinv,zinv,zlmel
	84	real ztv_est1,ztv_est2
[878]	85	real zcor,zdelta,zcvm5,qlbef
[1968]	86	real betalpha,zbetalpha
	87	real eps, afact
[878]	88	REAL REPS,RLvCp,DDT0
	89	PARAMETER (DDT0=.01)
	90	logical Zsat
[1403]	91	LOGICAL active(ngrid),activetmp(ngrid)
[1968]	92	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
[1026]	93	REAL c2(ngrid,klev)
[878]	94	Zsat=.false.
	95	! Initialisation
[1968]	96
[878]	97	RLvCp = RLVTT/RCPD
[1968]	98	fact_epsilon=0.002
	99	betalpha=0.9
[2046]	100	afact=2./3.
	101	fact_shell=0.85
[1026]	102
[1968]	103	zbetalpha=betalpha/(1.+betalpha)
	104
	105
	106	! Initialisations des variables r?elles
[1403]	107	if (1==1) then
	108	ztva(:,:)=ztv(:,:)
	109	ztva_est(:,:)=ztva(:,:)
[1968]	110	ztv_est(:,:)=ztv(:,:)
[1403]	111	ztla(:,:)=zthl(:,:)
	112	zqta(:,:)=po(:,:)
[1968]	113	zqla(:,:)=0.
[1403]	114	zha(:,:) = ztva(:,:)
	115	else
	116	ztva(:,:)=0.
[1968]	117	ztv_est(:,:)=0.
[1403]	118	ztva_est(:,:)=0.
	119	ztla(:,:)=0.
	120	zqta(:,:)=0.
	121	zha(:,:) =0.
	122	endif
[878]	123
[1403]	124	zqla_est(:,:)=0.
	125	zqsatth(:,:)=0.
	126	zqla(:,:)=0.
	127	detr_star(:,:)=0.
	128	entr_star(:,:)=0.
	129	alim_star(:,:)=0.
	130	alim_star_tot(:)=0.
	131	csc(:,:)=0.
	132	detr(:,:)=0.
	133	entr(:,:)=0.
	134	zw2(:,:)=0.
[1968]	135	zbuoy(:,:)=0.
	136	zbuoyjam(:,:)=0.
	137	gamma(:,:)=0.
	138	zeps(:,:)=0.
[1403]	139	w_est(:,:)=0.
	140	f_star(:,:)=0.
	141	wa_moy(:,:)=0.
	142	linter(:)=1.
[1968]	143	! linter(:)=1.
[1403]	144	! Initialisation des variables entieres
	145	lmix(:)=1
	146	lmix_bis(:)=2
	147	wmaxa(:)=0.
	148	lalim(:)=1
[972]	149
[1968]	150
[1403]	151	!-------------------------------------------------------------------------
	152	! On ne considere comme actif que les colonnes dont les deux premieres
	153	! couches sont instables.
	154	!-------------------------------------------------------------------------
	155	active(:)=ztv(:,1)>ztv(:,2)
	156
	157	!-------------------------------------------------------------------------
	158	! Definition de l'alimentation a l'origine dans thermcell_init
	159	!-------------------------------------------------------------------------
	160	do l=1,klev-1
[878]	161	do ig=1,ngrid
[1403]	162	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
	163	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
	164	& *sqrt(zlev(ig,l+1))
[1503]	165	lalim(ig)=l+1
[1403]	166	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
[2046]	167	! print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l)
[1403]	168	endif
[878]	169	enddo
	170	enddo
[1403]	171	do l=1,klev
	172	do ig=1,ngrid
	173	if (alim_star_tot(ig) > 1.e-10 ) then
	174	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
	175	endif
	176	enddo
[878]	177	enddo
[1403]	178	alim_star_tot(:)=1.
[878]	179
[972]	180
[1968]	181
	182
[1403]	183	!------------------------------------------------------------------------------
	184	! Calcul dans la premiere couche
	185	! On decide dans cette version que le thermique n'est actif que si la premiere
	186	! couche est instable.
[1968]	187	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
[1403]	188	! dans une couche l>1
	189	!------------------------------------------------------------------------------
	190	do ig=1,ngrid
[972]	191	! Le panache va prendre au debut les caracteristiques de l'air contenu
	192	! dans cette couche.
[1403]	193	if (active(ig)) then
	194	ztla(ig,1)=zthl(ig,1)
	195	zqta(ig,1)=po(ig,1)
	196	zqla(ig,1)=zl(ig,1)
	197	!cr: attention, prise en compte de f*(1)=1
	198	f_star(ig,2)=alim_star(ig,1)
	199	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
	200	& *(zlev(ig,2)-zlev(ig,1)) &
	201	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
	202	w_est(ig,2)=zw2(ig,2)
	203	endif
	204	enddo
[878]	205	!
[972]	206
[1403]	207	!==============================================================================
	208	!boucle de calcul de la vitesse verticale dans le thermique
	209	!==============================================================================
	210	do l=2,klev-1
	211	!==============================================================================
[972]	212
[878]	213
[1403]	214	! On decide si le thermique est encore actif ou non
	215	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
	216	do ig=1,ngrid
	217	active(ig)=active(ig) &
	218	& .and. zw2(ig,l)>1.e-10 &
	219	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
	220	enddo
[878]	221
	222
	223
	224	!---------------------------------------------------------------------------
[1403]	225	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
	226	! sans tenir compte du detrainement et de l'entrainement dans cette
	227	! couche
[1968]	228	! C'est a dire qu'on suppose
	229	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
[1403]	230	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
	231	! avant) a l'alimentation pour avoir un calcul plus propre
[878]	232	!---------------------------------------------------------------------------
[972]	233
[1968]	234	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
	235	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
	236	do ig=1,ngrid
[2046]	237	! print*,'active',active(ig),ig,l
[1968]	238	if(active(ig)) then
	239	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
[878]	240	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
[1403]	241	zta_est(ig,l)=ztva_est(ig,l)
[878]	242	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
	243	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
	244	& -zqla_est(ig,l))-zqla_est(ig,l))
[1968]	245
[878]	246
[1968]	247	!------------------------------------------------
	248	!AJAM:nouveau calcul de w?
	249	!------------------------------------------------
	250	zdz=zlev(ig,l+1)-zlev(ig,l)
	251	zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
	252	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
[1403]	253
[1968]	254	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	255	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
	256	zdw2=afact*zbuoy(ig,l)/fact_epsilon
	257	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
	258	! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
	259	! w_est(ig,l+1)=(zdz/zdzbis)Max(0.0001,exp(-zw2fact) &
	260	! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
	261	! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
	262	! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))zdw2+w_est(ig,l)exp(-zw2fact))
[2046]	263
	264	!--------------------------------------------------
	265	!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
	266	!--------------------------------------------------
	267	if (iflag_thermals_ed==8) then
	268	! Ancienne version
[1968]	269	w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
	270	& (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
	271	& (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
[2046]	272
	273	! Nouvelle version Arnaud
	274	else
	275	w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
	276	& (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
	277	& (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
	278	endif
	279
	280
	281	if (iflag_thermals_ed<6) then
	282	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
	283	! fact_epsilon=0.0005/(zalpha+0.025)**0.5
	284	! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
	285	fact_epsilon=0.0002/(zalpha+0.1)
	286	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	287	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
	288	zdw2=afact*zbuoy(ig,l)/fact_epsilon
	289	zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
	290	w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
	291	& (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
	292	& (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
	293
	294	endif
	295	!--------------------------------------------------
	296	!AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
	297	!on fait max(0.0001,.....)
	298	!--------------------------------------------------
	299
	300	! if (w_est(ig,l+1).lt.0.) then
[1968]	301	! w_est(ig,l+1)=zw2(ig,l)
[2046]	302	! w_est(ig,l+1)=0.0001
	303	! endif
	304
[1403]	305	endif
	306	enddo
[972]	307
[1968]	308
[1403]	309	!-------------------------------------------------
	310	!calcul des taux d'entrainement et de detrainement
	311	!-------------------------------------------------
[972]	312
[1403]	313	do ig=1,ngrid
	314	if (active(ig)) then
[1968]	315
	316	! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
	317	zw2m=w_est(ig,l+1)
	318	! zw2m=zw2(ig,l)
[1403]	319	zdz=zlev(ig,l+1)-zlev(ig,l)
[1968]	320	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	321	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
	322	zbuoybis=zbuoy(ig,l)
	323	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
	324	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
[972]	325
[1968]	326
	327	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
	328	! & afact*zbuoybis/zw2m - fact_epsilon )
[972]	329
[1968]	330	! entr_star(ig,l)=MAX(0.,f_star(ig,l)zdzzbetalpha* &
	331	! & afact*zbuoybis/zw2m - fact_epsilon )
[972]	332
[1968]	333	!Modif AJAM
	334
[2000]	335	lmel=fact_thermals_ed_dz*zlev(ig,l)
[2046]	336	zlmel=zlev(ig,l)+lmel
[1968]	337	! lmel=2.5*(zlev(ig,l)-zlev(ig,l-1))
	338	lt=l+1
[2046]	339	zdz2=zlev(ig,lt)-zlev(ig,l)
	340	!--------------------------------------------------
	341	!AJ052014: J'ai remplac? la boucle do par un do while
	342	! afin de faire moins de calcul dans la boucle
	343	!--------------------------------------------------
	344	do while (zdz2.lt.lmel)
	345	lt=lt+1
	346	zdz2=zlev(ig,lt)-zlev(ig,l)
	347	end do
	348
	349	zdz3=zlev(ig,lt)-zlev(ig,lt-1)
	350
	351	!--------------------------------------------------
	352	!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
	353	! on cherche o? se trouve l'altitude d'inversion
	354	! en calculant ztv1 (interpolation de la valeur de
	355	! theta au niveau lt en utilisant les niveaux lt-1 et
	356	! lt-2) et ztv2 (interpolation avec les niveaux lt+1
	357	! et lt+2). Si theta r?ellement calcul?e au niveau lt
	358	! comprise entre ztv1 et ztv2, alors il y a inversion
	359	! et on calcule son altitude zinv en supposant que ztv(lt)
	360	! est une combinaison lin?aire de ztv1 et ztv2.
	361	! Ensuite, on calcule la flottabilit? en comparant
	362	! la temp?rature de la couche l ? celle de l'air situ?
	363	! l+lmel plus haut, ce qui n?cessite de savoir quel fraction
	364	! de cet air est au-dessus ou en-dessous de l'inversion
	365	!--------------------------------------------------
	366
	367
	368	if (iflag_thermals_ed.lt.8) then
	369
	370	ztv1=(ztv(ig,lt-1)-ztv(ig,lt-2))*zlev(ig,lt)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
	371	& +(ztv(ig,lt-2)zlev(ig,lt-1)-ztv(ig,lt-1)zlev(ig,lt-2)) &
	372	& /(zlev(ig,lt-1)-zlev(ig,lt-2))
	373
	374	ztv2=(ztv(ig,lt+2)-ztv(ig,lt+1))*zlev(ig,lt)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
	375	& +(ztv(ig,lt+1)zlev(ig,lt+2)-ztv(ig,lt+2)zlev(ig,lt+1)) &
	376	& /(zlev(ig,lt+2)-zlev(ig,lt+1))
	377
	378	if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then
	379
	380	factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
	381	zinv=zlev(ig,lt-1)+factinv*(zlev(ig,lt)-zlev(ig,lt-1))
	382
	383	if (zlmel+0.5*zdz.ge.zinv) then
	384	if (zlmel-0.5*zdz.ge.zinv) then
	385
	386	ztv_est(ig,l)=(ztv(ig,lt+2)-ztv(ig,lt+1))(zlmel-0.zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
	387	& +(ztv(ig,lt+1)zlev(ig,lt+2)-ztv(ig,lt+2)zlev(ig,lt+1)) &
	388	& /(zlev(ig,lt+2)-zlev(ig,lt+1))
	389
	390	zbuoyjam(ig,l)=fact_shellRG(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l)
	391
	392	else
	393
	394	ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))0.5(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
	395	& +(ztv(ig,lt+1)zlev(ig,lt+2)-ztv(ig,lt+2)zlev(ig,lt+1)) &
	396	& /(zlev(ig,lt+2)-zlev(ig,lt+1))
	397	ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))0.5(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
	398	& +(ztv(ig,lt-2)zlev(ig,lt-1)-ztv(ig,lt-1)zlev(ig,lt-2)) &
	399	& /(zlev(ig,lt-1)-zlev(ig,lt-2))
	400	zbuoyjam(ig,l)=fact_shellRG(((zlmel+0.5zdz-zinv)/zdz)(ztva_est(ig,l)- &
	401	& ztv_est1)/ztv_est1+((zinv-zlmel+0.5zdz)/zdz)(ztva_est(ig,l)- &
	402	& ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l)
	403
	404	endif
	405
	406	else
	407
	408	ztv_est(ig,l)=(ztv(ig,lt-1)-ztv(ig,lt-2))(zlmel-0.zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
	409	& +(ztv(ig,lt-2)zlev(ig,lt-1)-ztv(ig,lt-1)zlev(ig,lt-2)) &
	410	& /(zlev(ig,lt-1)-zlev(ig,lt-2))
	411
	412	zbuoyjam(ig,l)=fact_shellRG(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l)
	413	! ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))0.5(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
	414	! & +(ztv(ig,lt+1)zlev(ig,lt+2)-ztv(ig,lt+2)zlev(ig,lt+1)) &
	415	! & /(zlev(ig,lt+2)-zlev(ig,lt+1))
	416	! ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))0.5(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
	417	! & +(ztv(ig,lt-2)zlev(ig,lt-1)-ztv(ig,lt-1)zlev(ig,lt-2)) &
	418	! & /(zlev(ig,lt-1)-zlev(ig,lt-2))
	419	! zbuoyjam(ig,l)=fact_shellRG(((zlmel+0.5zdz-zinv)/zdz)(ztva_est(ig,l)- &
	420	! & ztv_est1)/ztv_est1+((zinv-zlmel+0.5zdz)/zdz)(ztva_est(ig,l)- &
	421	! & ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l)
	422
	423
	424
	425
	426	! print*,'on est pass? par l?',l,lt,zbuoyjam(ig,l),zbuoy(ig,l)
	427	endif
	428
	429
	430	else
	431
[1968]	432	! ztv_est(ig,l)=(lmel/zdz2)*(ztv(ig,lt)-ztv(ig,l))+ztv(ig,l)
	433	! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)
[972]	434
[2046]	435	zbuoyjam(ig,l)=fact_shellRG(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
[1968]	436	& ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
[2046]	437	& ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
[972]	438
[1968]	439	! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
	440	! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
[2046]	441	! & po(ig,lt-1))/po(ig,lt-1))
	442
	443	endif
	444
[1968]	445	else
[972]	446
[2046]	447	zbuoyjam(ig,l)=1.RG(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- &
	448	& ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
	449	& ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
	450
	451	endif
	452
	453
[1968]	454	! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	455
[1998]	456	! entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
	457	! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
[1968]	458
[1998]	459	! entrbis=entr_star(ig,l)
[1968]	460
[2046]	461	if (iflag_thermals_ed.lt.6) then
	462	fact_epsilon=0.0002/(zalpha+0.1)
	463	endif
[1968]	464
[2046]	465	detr_star(ig,l)=f_star(ig,l)*zdz &
[1968]	466	& MAX(1.e-4, -afactzbetalpha*zbuoyjam(ig,l)/zw2m &
[2046]	467	& + 0.012(zdqt(ig,l)/zw2m)*0.5)
[1968]	468
[1998]	469	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
[1968]	470
[1998]	471	entr_star(ig,l)=f_star(ig,l)zdzzbetalpha*MAX(0., &
[2046]	472	& afact*zbuoy(ig,l)/zw2m - fact_epsilon)
[1968]	473
[1998]	474	! entr_star(ig,l)=Max(0.,f_star(ig,l)zdzzbetalpha* &
	475	! & afact*zbuoy(ig,l)/zw2m &
	476	! & - 1.*fact_epsilon)
	477
[1968]	478
[1403]	479	! En dessous de lalim, on prend le max de alim_star et entr_star pour
	480	! alim_star et 0 sinon
	481	if (l.lt.lalim(ig)) then
	482	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
	483	entr_star(ig,l)=0.
	484	endif
[1968]	485	! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
	486	! alim_star(ig,l)=entrbis
	487	! endif
[972]	488
[2046]	489	! print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
[1403]	490	! Calcul du flux montant normalise
	491	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	492	& -detr_star(ig,l)
[972]	493
[1403]	494	endif
	495	enddo
[972]	496
[1968]	497
[1403]	498	!----------------------------------------------------------------------------
	499	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	500	!---------------------------------------------------------------------------
	501	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
	502	do ig=1,ngrid
	503	if (activetmp(ig)) then
	504	Zsat=.false.
	505	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
	506	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	507	& /(f_star(ig,l+1)+detr_star(ig,l))
	508	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
	509	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	510	& /(f_star(ig,l+1)+detr_star(ig,l))
[972]	511
[1403]	512	endif
	513	enddo
[972]	514
[1968]	515	ztemp(:)=zpspsk(:,l)*ztla(:,l)
	516	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
[1403]	517	do ig=1,ngrid
	518	if (activetmp(ig)) then
	519	! on ecrit de maniere conservative (sat ou non)
	520	! T = Tl +Lv/Cp ql
[1968]	521	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
[1403]	522	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	523	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	524	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	525	zha(ig,l) = ztva(ig,l)
	526	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	527	& -zqla(ig,l))-zqla(ig,l))
[1968]	528	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
[1403]	529	zdz=zlev(ig,l+1)-zlev(ig,l)
[1968]	530	zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
	531	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
[2046]	532	fact_epsilon=0.002
[1968]	533	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	534	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
	535	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
	536	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
	537	! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
	538	zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
	539	& (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
	540	& (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
[2046]	541
	542	if (iflag_thermals_ed.lt.6) then
	543	zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
	544	! fact_epsilon=0.0005/(zalpha+0.025)**0.5
	545	! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
	546	fact_epsilon=0.0002/(zalpha+0.1)**1
	547	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	548	zw2factbis=fact_epsilon2.zdzbis/(1.+betalpha)
	549	zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
	550	zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
	551
	552	zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)(exp(-zw2fact) &
	553	& (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
	554	& (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
	555
	556	endif
	557
	558
[1403]	559	endif
	560	enddo
	561
	562	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
[878]	563	!
[1403]	564	!---------------------------------------------------------------------------
	565	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
	566	!---------------------------------------------------------------------------
[972]	567
[1503]	568	nbpb=0
[1403]	569	do ig=1,ngrid
	570	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	571	! stop'On tombe sur le cas particulier de thermcell_dry'
[1503]	572	! print*,'On tombe sur le cas particulier de thermcell_plume'
	573	nbpb=nbpb+1
[1403]	574	zw2(ig,l+1)=0.
	575	linter(ig)=l+1
	576	endif
[972]	577
[1403]	578	if (zw2(ig,l+1).lt.0.) then
	579	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	580	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	581	zw2(ig,l+1)=0.
[1998]	582	!+CR:04/05/12:correction calcul linter pour calcul de zmax continu
	583	elseif (f_star(ig,l+1).lt.0.) then
	584	linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) &
	585	& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
	586	zw2(ig,l+1)=0.
	587	!fin CR:04/05/12
[1403]	588	endif
[972]	589
[1403]	590	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
[972]	591
[1403]	592	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	593	! lmix est le niveau de la couche ou w (wa_moy) est maximum
	594	!on rajoute le calcul de lmix_bis
	595	if (zqla(ig,l).lt.1.e-10) then
	596	lmix_bis(ig)=l+1
	597	endif
	598	lmix(ig)=l+1
	599	wmaxa(ig)=wa_moy(ig,l+1)
	600	endif
	601	enddo
[972]	602
[1503]	603	if (nbpb>0) then
	604	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
	605	endif
	606
[1403]	607	!=========================================================================
	608	! FIN DE LA BOUCLE VERTICALE
	609	enddo
	610	!=========================================================================
[972]	611
[1403]	612	!on recalcule alim_star_tot
	613	do ig=1,ngrid
	614	alim_star_tot(ig)=0.
	615	enddo
	616	do ig=1,ngrid
	617	do l=1,lalim(ig)-1
	618	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	619	enddo
	620	enddo
	621
[972]	622
[1403]	623	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
[972]	624
[2046]	625	#undef wrgrads_thermcell
	626	#ifdef wrgrads_thermcell
	627	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
	628	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
	629	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
	630	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
	631	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
	632	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
	633	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
	634	#endif
	635
	636
[1968]	637	return
	638	end
	639
	640
[1998]	641
	642
[2046]	643
	644
	645
	646
	647
	648
	649
	650
	651
	652
	653
	654
	655
	656
	657
	658
	659
	660
	661
	662
	663
	664
	665
	666
	667
	668
	669
	670
	671
	672
	673
[1403]	674	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	675	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	676	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	677	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	678	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	679	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	680	SUBROUTINE thermcellV1_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, &
	681	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	682	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	683	& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
	684	& ,lev_out,lunout1,igout)
[972]	685
[1403]	686	!--------------------------------------------------------------------------
	687	!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
	688	! Version conforme a l'article de Rio et al. 2010.
	689	! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
	690	!--------------------------------------------------------------------------
[878]	691
[1403]	692	IMPLICIT NONE
[972]	693
[1403]	694	#include "YOMCST.h"
	695	#include "YOETHF.h"
	696	#include "FCTTRE.h"
	697	#include "iniprint.h"
	698	#include "thermcell.h"
[972]	699
[1403]	700	INTEGER itap
	701	INTEGER lunout1,igout
	702	INTEGER ngrid,klev
	703	REAL ptimestep
	704	REAL ztv(ngrid,klev)
	705	REAL zthl(ngrid,klev)
	706	REAL po(ngrid,klev)
	707	REAL zl(ngrid,klev)
	708	REAL rhobarz(ngrid,klev)
	709	REAL zlev(ngrid,klev+1)
	710	REAL pplev(ngrid,klev+1)
	711	REAL pphi(ngrid,klev)
	712	REAL zpspsk(ngrid,klev)
	713	REAL alim_star(ngrid,klev)
	714	REAL f0(ngrid)
	715	INTEGER lalim(ngrid)
	716	integer lev_out ! niveau pour les print
[1503]	717	integer nbpb
[1403]	718
	719	real alim_star_tot(ngrid)
[878]	720
[1403]	721	REAL ztva(ngrid,klev)
	722	REAL ztla(ngrid,klev)
	723	REAL zqla(ngrid,klev)
	724	REAL zqta(ngrid,klev)
	725	REAL zha(ngrid,klev)
[878]	726
[1403]	727	REAL detr_star(ngrid,klev)
	728	REAL coefc
	729	REAL entr_star(ngrid,klev)
	730	REAL detr(ngrid,klev)
	731	REAL entr(ngrid,klev)
	732
	733	REAL csc(ngrid,klev)
	734
	735	REAL zw2(ngrid,klev+1)
	736	REAL w_est(ngrid,klev+1)
	737	REAL f_star(ngrid,klev+1)
	738	REAL wa_moy(ngrid,klev+1)
	739
	740	REAL ztva_est(ngrid,klev)
	741	REAL zqla_est(ngrid,klev)
	742	REAL zqsatth(ngrid,klev)
	743	REAL zta_est(ngrid,klev)
[1998]	744	REAL zbuoyjam(ngrid,klev)
[1403]	745	REAL ztemp(ngrid),zqsat(ngrid)
	746	REAL zdw2
	747	REAL zw2modif
	748	REAL zw2fact
	749	REAL zeps(ngrid,klev)
	750
	751	REAL linter(ngrid)
	752	INTEGER lmix(ngrid)
	753	INTEGER lmix_bis(ngrid)
	754	REAL wmaxa(ngrid)
	755
	756	INTEGER ig,l,k
	757
	758	real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m
	759	real zbuoybis
	760	real zcor,zdelta,zcvm5,qlbef,zdz2
	761	real betalpha,zbetalpha
	762	real eps, afact
	763	REAL REPS,RLvCp,DDT0
	764	PARAMETER (DDT0=.01)
	765	logical Zsat
	766	LOGICAL active(ngrid),activetmp(ngrid)
	767	REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2
	768	REAL c2(ngrid,klev)
	769	Zsat=.false.
	770	! Initialisation
	771
	772	RLvCp = RLVTT/RCPD
	773	fact_epsilon=0.002
	774	betalpha=0.9
	775	afact=2./3.
	776
	777	zbetalpha=betalpha/(1.+betalpha)
	778
	779
	780	! Initialisations des variables reeles
[1998]	781	if (1==1) then
[1403]	782	ztva(:,:)=ztv(:,:)
	783	ztva_est(:,:)=ztva(:,:)
	784	ztla(:,:)=zthl(:,:)
	785	zqta(:,:)=po(:,:)
	786	zha(:,:) = ztva(:,:)
	787	else
	788	ztva(:,:)=0.
	789	ztva_est(:,:)=0.
	790	ztla(:,:)=0.
	791	zqta(:,:)=0.
	792	zha(:,:) =0.
	793	endif
	794
	795	zqla_est(:,:)=0.
	796	zqsatth(:,:)=0.
	797	zqla(:,:)=0.
	798	detr_star(:,:)=0.
	799	entr_star(:,:)=0.
	800	alim_star(:,:)=0.
	801	alim_star_tot(:)=0.
	802	csc(:,:)=0.
	803	detr(:,:)=0.
	804	entr(:,:)=0.
	805	zw2(:,:)=0.
	806	zbuoy(:,:)=0.
[1998]	807	zbuoyjam(:,:)=0.
[1403]	808	gamma(:,:)=0.
	809	zeps(:,:)=0.
	810	w_est(:,:)=0.
	811	f_star(:,:)=0.
	812	wa_moy(:,:)=0.
	813	linter(:)=1.
	814	! linter(:)=1.
	815	! Initialisation des variables entieres
	816	lmix(:)=1
	817	lmix_bis(:)=2
	818	wmaxa(:)=0.
	819	lalim(:)=1
	820
	821
	822	!-------------------------------------------------------------------------
	823	! On ne considere comme actif que les colonnes dont les deux premieres
	824	! couches sont instables.
	825	!-------------------------------------------------------------------------
	826	active(:)=ztv(:,1)>ztv(:,2)
	827
	828	!-------------------------------------------------------------------------
	829	! Definition de l'alimentation a l'origine dans thermcell_init
	830	!-------------------------------------------------------------------------
	831	do l=1,klev-1
	832	do ig=1,ngrid
	833	if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
	834	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
	835	& *sqrt(zlev(ig,l+1))
	836	lalim(ig)=l+1
	837	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
[878]	838	endif
[1403]	839	enddo
	840	enddo
	841	do l=1,klev
	842	do ig=1,ngrid
	843	if (alim_star_tot(ig) > 1.e-10 ) then
	844	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
	845	endif
	846	enddo
	847	enddo
	848	alim_star_tot(:)=1.
[878]	849
	850
[972]	851
[1403]	852	!------------------------------------------------------------------------------
	853	! Calcul dans la premiere couche
	854	! On decide dans cette version que le thermique n'est actif que si la premiere
	855	! couche est instable.
[1968]	856	! Pourrait etre change si on veut que le thermiques puisse se d??clencher
[1403]	857	! dans une couche l>1
	858	!------------------------------------------------------------------------------
	859	do ig=1,ngrid
	860	! Le panache va prendre au debut les caracteristiques de l'air contenu
	861	! dans cette couche.
	862	if (active(ig)) then
	863	ztla(ig,1)=zthl(ig,1)
	864	zqta(ig,1)=po(ig,1)
	865	zqla(ig,1)=zl(ig,1)
	866	!cr: attention, prise en compte de f*(1)=1
	867	f_star(ig,2)=alim_star(ig,1)
	868	zw2(ig,2)=2.RG(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
	869	& *(zlev(ig,2)-zlev(ig,1)) &
	870	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
	871	w_est(ig,2)=zw2(ig,2)
	872	endif
	873	enddo
	874	!
[1026]	875
[1403]	876	!==============================================================================
	877	!boucle de calcul de la vitesse verticale dans le thermique
	878	!==============================================================================
	879	do l=2,klev-1
	880	!==============================================================================
[1026]	881
	882
[1403]	883	! On decide si le thermique est encore actif ou non
	884	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
	885	do ig=1,ngrid
	886	active(ig)=active(ig) &
	887	& .and. zw2(ig,l)>1.e-10 &
	888	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
	889	enddo
[1026]	890
	891
	892
[1403]	893	!---------------------------------------------------------------------------
	894	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
	895	! sans tenir compte du detrainement et de l'entrainement dans cette
	896	! couche
	897	! C'est a dire qu'on suppose
	898	! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
	899	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
	900	! avant) a l'alimentation pour avoir un calcul plus propre
	901	!---------------------------------------------------------------------------
	902
	903	ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
	904	call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
	905
	906	do ig=1,ngrid
	907	! print*,'active',active(ig),ig,l
	908	if(active(ig)) then
	909	zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
	910	ztva_est(ig,l) = ztla(ig,l-1)zpspsk(ig,l)+RLvCpzqla_est(ig,l)
	911	zta_est(ig,l)=ztva_est(ig,l)
	912	ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
	913	ztva_est(ig,l) = ztva_est(ig,l)(1.+RETV(zqta(ig,l-1) &
	914	& -zqla_est(ig,l))-zqla_est(ig,l))
	915
	916	!------------------------------------------------
[1968]	917	!AJAM:nouveau calcul de w?
[1403]	918	!------------------------------------------------
	919	zdz=zlev(ig,l+1)-zlev(ig,l)
	920	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	921
	922	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	923	zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon)
	924	w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
	925
	926
	927	if (w_est(ig,l+1).lt.0.) then
	928	w_est(ig,l+1)=zw2(ig,l)
	929	endif
[1026]	930	endif
[1403]	931	enddo
[1026]	932
	933
[1403]	934	!-------------------------------------------------
	935	!calcul des taux d'entrainement et de detrainement
	936	!-------------------------------------------------
[1026]	937
[1403]	938	do ig=1,ngrid
	939	if (active(ig)) then
[1026]	940
[1403]	941	zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
	942	zw2m=w_est(ig,l+1)
	943	zdz=zlev(ig,l+1)-zlev(ig,l)
	944	zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	945	! zbuoybis=zbuoy(ig,l)+RG*0.1/300.
	946	zbuoybis=zbuoy(ig,l)
	947	zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
	948	zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l)
[1026]	949
[1403]	950
	951	entr_star(ig,l)=f_star(ig,l)zdz zbetalpha*MAX(0., &
	952	& afact*zbuoybis/zw2m - fact_epsilon )
	953
	954
	955	detr_star(ig,l)=f_star(ig,l)*zdz &
	956	& MAX(1.e-3, -afactzbetalpha*zbuoy(ig,l)/zw2m &
	957	& + 0.012(zdqt(ig,l)/zw2m)*0.5 )
	958
	959	! En dessous de lalim, on prend le max de alim_star et entr_star pour
	960	! alim_star et 0 sinon
	961	if (l.lt.lalim(ig)) then
	962	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
	963	entr_star(ig,l)=0.
	964	endif
	965
	966	! Calcul du flux montant normalise
[878]	967	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	968	& -detr_star(ig,l)
	969
[1403]	970	endif
	971	enddo
	972
	973
[878]	974	!----------------------------------------------------------------------------
	975	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	976	!---------------------------------------------------------------------------
[1403]	977	activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
	978	do ig=1,ngrid
	979	if (activetmp(ig)) then
	980	Zsat=.false.
	981	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
[878]	982	& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) &
	983	& /(f_star(ig,l+1)+detr_star(ig,l))
[1403]	984	zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ &
[878]	985	& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) &
	986	& /(f_star(ig,l+1)+detr_star(ig,l))
	987
[1403]	988	endif
	989	enddo
	990
	991	ztemp(:)=zpspsk(:,l)*ztla(:,l)
	992	call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
	993
	994	do ig=1,ngrid
	995	if (activetmp(ig)) then
[878]	996	! on ecrit de maniere conservative (sat ou non)
	997	! T = Tl +Lv/Cp ql
[1403]	998	zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l))
[878]	999	ztva(ig,l) = ztla(ig,l)zpspsk(ig,l)+RLvCpzqla(ig,l)
	1000	ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
	1001	!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
	1002	zha(ig,l) = ztva(ig,l)
	1003	ztva(ig,l) = ztva(ig,l)(1.+RETV(zqta(ig,l) &
	1004	& -zqla(ig,l))-zqla(ig,l))
[1403]	1005	zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	1006	zdz=zlev(ig,l+1)-zlev(ig,l)
	1007	zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
[878]	1008
[1403]	1009	zw2fact=fact_epsilon2.zdz/(1.+betalpha)
	1010	zdw2=afact*zbuoy(ig,l)/(fact_epsilon)
	1011	zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
	1012	endif
	1013	enddo
[1026]	1014
[972]	1015	if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
[878]	1016	!
[1403]	1017	!---------------------------------------------------------------------------
[878]	1018	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
[1403]	1019	!---------------------------------------------------------------------------
[878]	1020
[1503]	1021	nbpb=0
[1403]	1022	do ig=1,ngrid
[878]	1023	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	1024	! stop'On tombe sur le cas particulier de thermcell_dry'
[1503]	1025	! print*,'On tombe sur le cas particulier de thermcell_plume'
	1026	nbpb=nbpb+1
[878]	1027	zw2(ig,l+1)=0.
	1028	linter(ig)=l+1
[2046]	1029	endif
[878]	1030
	1031	if (zw2(ig,l+1).lt.0.) then
	1032	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	1033	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	1034	zw2(ig,l+1)=0.
	1035	endif
	1036
	1037	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	1038
	1039	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	1040	! lmix est le niveau de la couche ou w (wa_moy) est maximum
[1026]	1041	!on rajoute le calcul de lmix_bis
	1042	if (zqla(ig,l).lt.1.e-10) then
	1043	lmix_bis(ig)=l+1
	1044	endif
[878]	1045	lmix(ig)=l+1
	1046	wmaxa(ig)=wa_moy(ig,l+1)
	1047	endif
[1403]	1048	enddo
	1049
[1503]	1050	if (nbpb>0) then
	1051	print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume'
	1052	endif
	1053
[1403]	1054	!=========================================================================
	1055	! FIN DE LA BOUCLE VERTICALE
[878]	1056	enddo
[1403]	1057	!=========================================================================
[878]	1058
[1403]	1059	!on recalcule alim_star_tot
	1060	do ig=1,ngrid
	1061	alim_star_tot(ig)=0.
	1062	enddo
	1063	do ig=1,ngrid
	1064	do l=1,lalim(ig)-1
	1065	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	1066	enddo
	1067	enddo
	1068
	1069
[972]	1070	if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
[878]	1071
[2046]	1072	#undef wrgrads_thermcell
	1073	#ifdef wrgrads_thermcell
	1074	call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta ','esta ')
	1075	call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta ','dsta ')
	1076	call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy ','buoy ')
	1077	call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt ','dqt ')
	1078	call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est ','w_est ')
	1079	call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2 ','w_es2 ')
	1080	call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A ','zw2A ')
	1081	#endif
[1998]	1082
	1083
[1403]	1084	return
	1085	end

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

Download in other formats:

Original Format