Context Navigation

thermcell_plume_6A.F90 @ 3592

Last change on this file since 3592 was 3451, checked in by fhourdin, 6 years ago

Saving old version of thermcell_plume in thermcell_plume_6A.F90
It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
to the 5B and 6A versions used for CMIP5 and CMIP6.
The latest was previously named thermcellV1_plume.
The new thermcell_plume is a clean version (removing obsolete
options) of thermcell_plume_6A.
The 3 versions are controled by
flag_thermals_ed <= 9 thermcell_plume_6A

<= 19 thermcell_plume_5B
else thermcell_plume (default 20 for convergence with 6A)

Fredho

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

File size: 40.7 KB

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

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

Download in other formats:

Original Format