Context Navigation

thermcell_plume.F90 @ 2394

Last change on this file since 2394 was 2392, checked in by fhourdin, 9 years ago
Retour a 1+1=2. Probleme dans thermcell_alim. Bug fixing
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.2 KB

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

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

Download in other formats:

Original Format