Context Navigation

thermcell_plume.F90 @ 2343

Last change on this file since 2343 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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.6 KB

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

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

Download in other formats:

Original Format