Context Navigation

thermcell_plume.F90 @ 5223

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

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

Download in other formats:

Original Format