Context Navigation

thermcell_plume.F90 @ 2237

Last change on this file since 2237 was 2189, checked in by Laurent Fairhead, 10 years ago

Suppression print superflus

Deleting some superfluous prints

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: 41.3 KB

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

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

Download in other formats:

Original Format