Context Navigation

thermcell_plume.F90 @ 2186

Last change on this file since 2186 was 2159, checked in by jyg, 10 years ago
1/ Splitting of the boundary layer : the climbing down and up of Pbl_surface is split between the off-wake and wake regions ; the thermal scheme is applied only to the off-wake region. 2/ Elimination of wake_scal and calwake_scal.
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: 42.1 KB

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

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

Download in other formats:

Original Format