Context Navigation

thermcell_plume.F90 @ 2128

Last change on this file since 2128 was 2106, checked in by fhourdin, 10 years ago
Version moidifiÃe des thermiques pour une meilleure reprÃ©snation des stratocumulus + contrÃ´le par lecture de .def d'un ertain nombre de paramÃ¨tres internes. Modified version of thermals for stratocumulus + reading new free parameters. Modification by Arnaud Jam and Frederic Hourdin
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: 40.2 KB

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

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

Download in other formats:

Original Format