Context Navigation

thermcell_main_mars.F90 @ 607

Last change on this file since 607 was 593, checked in by acolaitis, 13 years ago

Minor change to internal parameters of the thermal: smoothed detrainment profiles, reduced the contribution of upward subsidence to entrainment drag, and increased pressure torque induced drag. Consequences : more precise prediction of zmax (it was overestimated by 5 to 10 percent during the late afternoon), no consequence on the thermal structure, better predictions for internal variables (much better w profile and detrainment rate profile). Tested on both high resolution and low resolution (33 level grid from me and 36 level grid from Forget)

File size: 52.9 KB

Rev	Line
[161]	1	!
	2	!
[185]	3	SUBROUTINE thermcell_main_mars(ptimestep &
[161]	4	& ,pplay,pplev,pphi,zlev,zlay &
	5	& ,pu,pv,pt,pq,pq2 &
	6	& ,pduadj,pdvadj,pdtadj,pdqadj,pdq2adj &
[185]	7	& ,fm,entr,detr,lmax,zmax &
[161]	8	& ,r_aspect &
	9	& ,zw2,fraca &
	10	& ,zpopsk,ztla,heatFlux,heatFlux_down &
	11	& ,buoyancyOut, buoyancyEst)
	12
	13	IMPLICIT NONE
	14
	15	!=======================================================================
	16	! Mars version of thermcell_main. Author : A Colaitis
	17	!=======================================================================
[185]	18
	19	#include "dimensions.h"
	20	#include "dimphys.h"
	21	#include "comcstfi.h"
[508]	22	#include "tracer.h"
	23	#include "callkeys.h"
[185]	24
[161]	25	! ============== INPUTS ==============
	26
	27	REAL, INTENT(IN) :: ptimestep,r_aspect
[185]	28	REAL, INTENT(IN) :: pt(ngridmx,nlayermx)
	29	REAL, INTENT(IN) :: pu(ngridmx,nlayermx)
	30	REAL, INTENT(IN) :: pv(ngridmx,nlayermx)
	31	REAL, INTENT(IN) :: pq(ngridmx,nlayermx,nqmx)
	32	REAL, INTENT(IN) :: pq2(ngridmx,nlayermx)
	33	REAL, INTENT(IN) :: pplay(ngridmx,nlayermx)
	34	REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1)
	35	REAL, INTENT(IN) :: pphi(ngridmx,nlayermx)
	36	REAL, INTENT(IN) :: zlay(ngridmx,nlayermx)
	37	REAL, INTENT(IN) :: zlev(ngridmx,nlayermx+1)
[161]	38
	39	! ============== OUTPUTS ==============
	40
[185]	41	REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx)
[337]	42	REAL :: pduadj(ngridmx,nlayermx)
	43	REAL :: pdvadj(ngridmx,nlayermx)
	44	REAL :: pdqadj(ngridmx,nlayermx,nqmx)
[313]	45	! REAL, INTENT(OUT) :: pdq2adj(ngridmx,nlayermx)
	46	REAL :: pdq2adj(ngridmx,nlayermx)
[185]	47	REAL, INTENT(OUT) :: zw2(ngridmx,nlayermx+1)
[161]	48
	49	! Diagnostics
[185]	50	REAL, INTENT(OUT) :: heatFlux(ngridmx,nlayermx) ! interface heatflux
[313]	51	REAL, INTENT(OUT) :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
	52	! REAL, INTENT(OUT) :: buoyancyOut(ngridmx,nlayermx) ! interlayer buoyancy term
	53	! REAL, INTENT(OUT) :: buoyancyEst(ngridmx,nlayermx) ! interlayer estimated buoyancy term
[161]	54
	55	! dummy variables when output not needed :
	56
[185]	57	! REAL :: heatFlux(ngridmx,nlayermx) ! interface heatflux
	58	! REAL :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
[313]	59	REAL :: buoyancyOut(ngridmx,nlayermx) ! interlayer buoyancy term
	60	REAL :: buoyancyEst(ngridmx,nlayermx) ! interlayer estimated buoyancy term
[161]	61
	62
	63	! ============== LOCAL ================
	64
	65	INTEGER ig,k,l,ll,iq
[185]	66	INTEGER lmax(ngridmx),lmin(ngridmx),lalim(ngridmx)
	67	REAL linter(ngridmx)
	68	REAL zmax(ngridmx)
	69	REAL ztva(ngridmx,nlayermx),zw_est(ngridmx,nlayermx+1),ztva_est(ngridmx,nlayermx)
	70	REAL zmax_sec(ngridmx)
	71	REAL zh(ngridmx,nlayermx)
	72	REAL zdthladj(ngridmx,nlayermx)
	73	REAL zdthladj_down(ngridmx,nlayermx)
	74	REAL ztvd(ngridmx,nlayermx)
	75	REAL ztv(ngridmx,nlayermx)
	76	REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx),zo(ngridmx,nlayermx)
	77	REAL zva(ngridmx,nlayermx)
	78	REAL zua(ngridmx,nlayermx)
[161]	79
[185]	80	REAL zta(ngridmx,nlayermx)
	81	REAL fraca(ngridmx,nlayermx+1)
	82	REAL q2(ngridmx,nlayermx)
	83	REAL rho(ngridmx,nlayermx),rhobarz(ngridmx,nlayermx),masse(ngridmx,nlayermx)
	84	REAL zpopsk(ngridmx,nlayermx)
[161]	85
[185]	86	REAL wmax(ngridmx)
	87	REAL wmax_sec(ngridmx)
	88	REAL fm(ngridmx,nlayermx+1),entr(ngridmx,nlayermx),detr(ngridmx,nlayermx)
[161]	89
[185]	90	REAL fm_down(ngridmx,nlayermx+1)
[161]	91
[185]	92	REAL ztla(ngridmx,nlayermx)
[161]	93
[185]	94	REAL f_star(ngridmx,nlayermx+1),entr_star(ngridmx,nlayermx)
	95	REAL detr_star(ngridmx,nlayermx)
	96	REAL alim_star_tot(ngridmx)
	97	REAL alim_star(ngridmx,nlayermx)
	98	REAL alim_star_clos(ngridmx,nlayermx)
	99	REAL f(ngridmx)
[161]	100
[185]	101	REAL teta_th_int(ngridmx,nlayermx)
	102	REAL teta_env_int(ngridmx,nlayermx)
	103	REAL teta_down_int(ngridmx,nlayermx)
[161]	104
	105	CHARACTER (LEN=20) :: modname
	106	CHARACTER (LEN=80) :: abort_message
	107
	108	! ============= PLUME VARIABLES ============
	109
[185]	110	REAL w_est(ngridmx,nlayermx+1)
	111	REAL wa_moy(ngridmx,nlayermx+1)
	112	REAL wmaxa(ngridmx)
	113	REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
[508]	114	LOGICAL activecell(ngridmx),activetmp(ngridmx)
[544]	115	REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega,adalim
[290]	116	INTEGER tic
[161]	117
	118	! ==========================================
	119
	120	! ============= HEIGHT VARIABLES ===========
	121
[185]	122	REAL num(ngridmx)
	123	REAL denom(ngridmx)
	124	REAL zlevinter(ngridmx)
[161]	125
	126	! =========================================
	127
	128	! ============= CLOSURE VARIABLES =========
	129
[185]	130	REAL zdenom(ngridmx)
	131	REAL alim_star2(ngridmx)
	132	REAL alim_star_tot_clos(ngridmx)
[161]	133	INTEGER llmax
	134
	135	! =========================================
	136
	137	! ============= FLUX2 VARIABLES ===========
	138
	139	INTEGER ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
	140	INTEGER ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
	141	REAL zfm
	142	REAL f_old,ddd0,eee0,ddd,eee,zzz
	143	REAL fomass_max,alphamax
	144
	145	! =========================================
	146
	147	! ============= DTETA VARIABLES ===========
	148
	149	! rien : on prends la divergence du flux turbulent
	150
	151	! =========================================
	152
	153	! ============= DV2 VARIABLES =============
	154	! not used for now
	155
[185]	156	REAL qa(ngridmx,nlayermx),detr_dv2(ngridmx,nlayermx),zf,zf2
	157	REAL wvd(ngridmx,nlayermx+1),wud(ngridmx,nlayermx+1)
	158	REAL gamma0(ngridmx,nlayermx+1),gamma(ngridmx,nlayermx+1)
	159	REAL ue(ngridmx,nlayermx),ve(ngridmx,nlayermx)
	160	LOGICAL ltherm(ngridmx,nlayermx)
	161	REAL dua(ngridmx,nlayermx),dva(ngridmx,nlayermx)
[161]	162	INTEGER iter
	163	INTEGER nlarga0
	164
	165	! =========================================
	166
[508]	167	! ============== Theta_M Variables ========
	168
	169	INTEGER ico2
	170	SAVE ico2
	171	REAL m_co2, m_noco2, A , B
	172	SAVE A, B
	173	LOGICAL firstcall
	174	save firstcall
	175	data firstcall/.true./
	176	REAL zhc(ngridmx,nlayermx)
	177	REAL zdzfull(ngridmx,nlayermx)
	178	REAL ratiom(ngridmx,nlayermx)
	179
	180	! =========================================
	181
[161]	182	!-----------------------------------------------------------------------
	183	! initialisation:
	184	! ---------------
	185
[336]	186	entr(:,:)=0.
	187	detr(:,:)=0.
	188	fm(:,:)=0.
[337]	189	! zu(:,:)=pu(:,:)
	190	! zv(:,:)=pv(:,:)
[508]	191	zhc(:,:)=pt(:,:)/zpopsk(:,:)
[161]	192
[508]	193	! **********************************************************************
	194	! Taking into account vertical molar mass gradients
	195	! **********************************************************************
	196
	197	if(firstcall) then
	198	ico2=0
	199	if (tracer) then
	200	! Prepare Special treatment if one of the tracers is CO2 gas
	201	do iq=1,nqmx
	202	if (noms(iq).eq."co2") then
	203	ico2=iq
	204	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
	205	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
	206	! Compute A and B coefficient use to compute
	207	! mean molecular mass Mair defined by
	208	! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
	209	! 1/Mair = A*q(ico2) + B
	210	A =(1/m_co2 - 1/m_noco2)
	211	B=1/m_noco2
	212	end if
	213	enddo
	214	endif
	215
	216	firstcall=.false.
	217	endif !of if firstcall
	218
	219	!.......................................................................
	220	! Special treatment for co2
	221	!.......................................................................
	222
	223	if (ico2.ne.0) then
	224	! Special case if one of the tracers is CO2 gas
	225	DO l=1,nlayermx
	226	DO ig=1,ngridmx
	227	ztv(ig,l) = zhc(ig,l)(Apq(ig,l,ico2)+B)
	228	ENDDO
	229	ENDDO
	230	else
	231	ztv(:,:)=zhc(:,:)
	232	end if
	233
	234
[161]	235	!------------------------------------------------------------------------
	236	! --------------------
	237	!
	238	!
	239	! + + + + + + + + + + +
	240	!
	241	!
	242	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
	243	! wh,wt,wo ...
	244	!
	245	! + + + + + + + + + + + zh,zu,zv,zo,rho
	246	!
	247	!
	248	! -------------------- zlev(1)
	249	! \\\\\\\\\\\\\\\\\\\\
	250	!
	251	!
	252
	253	!-----------------------------------------------------------------------
	254	! Calcul des altitudes des couches
	255	!-----------------------------------------------------------------------
	256
[185]	257	! do l=2,nlayermx
	258	! zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/g
[161]	259	! enddo
	260	! zlev(:,1)=0.
[185]	261	! zlev(:,nlayermx+1)=(2.*pphi(:,nlayermx)-pphi(:,nlayermx-1))/g
[161]	262
[185]	263	! zlay(:,:)=pphi(:,:)/g
[161]	264	!-----------------------------------------------------------------------
	265	! Calcul des densites
	266	!-----------------------------------------------------------------------
	267
[185]	268	rho(:,:)=pplay(:,:)/(r*pt(:,:))
[161]	269
	270	rhobarz(:,1)=rho(:,1)
	271
[185]	272	do l=2,nlayermx
[161]	273	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
	274	enddo
	275
	276	!calcul de la masse
[185]	277	do l=1,nlayermx
	278	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
[161]	279	enddo
	280
	281
	282	!------------------------------------------------------------------
	283	!
	284	! /\|\
	285	! -------- \| F_k+1 -------
	286	! ----> D_k
	287	! /\|\ <---- E_k , A_k
	288	! -------- \| F_k ---------
	289	! ----> D_k-1
	290	! <---- E_k-1 , A_k-1
	291	!
	292	!
	293	! ---------------------------
	294	!
	295	! ----- F_lmax+1=0 ---------- \
	296	! lmax (zmax) \|
	297	! --------------------------- \|
	298	! \|
	299	! --------------------------- \|
	300	! \|
	301	! --------------------------- \|
	302	! \|
	303	! --------------------------- \|
	304	! \|
	305	! --------------------------- \|
	306	! \| E
	307	! --------------------------- \| D
	308	! \|
	309	! --------------------------- \|
	310	! \|
	311	! --------------------------- \ \|
	312	! lalim \| \|
	313	! --------------------------- \| \|
	314	! \| \|
	315	! --------------------------- \| \|
	316	! \| A \|
	317	! --------------------------- \| \|
	318	! \| \|
	319	! --------------------------- \| \|
	320	! lmin (=1 pour le moment) \| \|
	321	! ----- F_lmin=0 ------------ / /
	322	!
	323	! ---------------------------
	324	! //////////////////////////
	325	!
	326
	327	!=============================================================================
	328	! Calculs initiaux ne faisant pas intervenir les changements de phase
	329	!=============================================================================
	330
	331	!------------------------------------------------------------------
	332	! 1. alim_star est le profil vertical de l'alimentation a la base du
	333	! panache thermique, calcule a partir de la flotabilite de l'air sec
	334	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
	335	!------------------------------------------------------------------
	336	!
	337	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
	338	lmin=1
	339
	340	!-----------------------------------------------------------------------------
	341	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
	342	! panache sec conservatif (e=d=0) alimente selon alim_star
	343	! Il s'agit d'un calcul de type CAPE
	344	! zmax_sec est utilise pour determiner la geometrie du thermique.
	345	!------------------------------------------------------------------------------
	346	!---------------------------------------------------------------------------------
	347	!calcul du melange et des variables dans le thermique
	348	!--------------------------------------------------------------------------------
	349
	350	! ===========================================================================
	351	! ===================== PLUME ===============================================
	352	! ===========================================================================
	353
	354	! Initialisations des variables reeles
	355	ztva(:,:)=ztv(:,:)
	356	ztva_est(:,:)=ztva(:,:)
	357	ztla(:,:)=0.
[173]	358	zdz=0.
[161]	359	zbuoy(:,:)=0.
	360	w_est(:,:)=0.
	361	f_star(:,:)=0.
	362	wa_moy(:,:)=0.
	363	linter(:)=1.
[185]	364
[496]	365	! --------------------------------------------------------------------------
	366	! -------------- MAIN PARAMETERS FOR THERMALS MODEL ------------------------
	367	! -------------- see thermiques.pro and getfit.py -------------------------
	368
[185]	369	! a1=2.5 ; b1=0.0015 ; ae=0.045 ; be = 0.6 ! svn baseline
	370
[496]	371	! Using broad downdraft selection
	372	! a1=1.60226 ; b1=0.0006 ; ae=0.0454 ; be = 0.57
	373	! ad = 0.0005114 ; bd = -0.662
	374	! fdfu = -1.9
[185]	375
[496]	376	! Using conditional sampling downdraft selection
	377	a1=1.4716 ; b1=0.0005698 ; ae=0.03683 ; be = 0.57421
	378	ad = 0.00048088 ; bd = -0.6697
	379	fdfu = -1.3
[532]	380	a1inv=a1
	381	b1inv=b1
	382	omega=0.
[544]	383	adalim=0.
[496]	384
[532]	385	! One good config for 34/35 levels
	386	! a1inv=a1
	387	! b1inv=b1
	388	! be=1.1*be
	389
[512]	390	! Best configuration for 222 levels:
[532]	391
[544]	392	! omega=0.06
	393	! b1=0.
	394	! a1=1.
	395	! a1inv=0.25*a1
	396	! b1inv=0.0002
	397	!!
	398	!!
	399	!! ae=0.9*ae
[512]	400
[532]	401	! Best config for norad 222 levels:
	402
[593]	403	! omega=0.06
	404	omega=0.04
[544]	405	a1=1.
[593]	406	! b1=0.
	407	b1=0.0001
[544]	408	a1inv=a1
	409	be=1.1*be
	410	ad = 0.0004
	411	adalim=0.
[546]	412	b1inv=0.00025
[593]	413	! b1inv = 0.0003
[532]	414
[546]	415	! omega=0.06
[532]	416	! Trying stuff :
	417
[546]	418	! ad=0.00035
	419	! ae=0.95*ae
	420	! b1=0.00055
[544]	421	! omega=0.04
[546]	422	!
	423	! ad = 0.0003
	424	! ae=0.9*ae
	425
	426	! omega=0.04
[544]	427	!! b1=0.
	428	! a1=1.
	429	! a1inv=a1
	430	! b1inv=0.0005689
	431	!! be=1.1*be
	432	!! ae=0.96*ae
	433
	434
[532]	435	! omega=0.06
	436	! a1=1.
	437	! b1=0.
	438	! be=be
	439	! a1inv=0.25*a1
	440	! b1inv=0.0002
	441	! ad=1.1*ad
	442	! ae=1.*ae
[496]	443	! --------------------------------------------------------------------------
	444	! --------------------------------------------------------------------------
	445	! --------------------------------------------------------------------------
	446
[161]	447	! Initialisation des variables entieres
	448	wmaxa(:)=0.
	449	lalim(:)=1
	450
	451	!-------------------------------------------------------------------------
	452	! On ne considere comme actif que les colonnes dont les deux premieres
	453	! couches sont instables.
	454	!-------------------------------------------------------------------------
[508]	455	activecell(:)=ztv(:,1)>ztv(:,2)
[185]	456	do ig=1,ngridmx
[161]	457	if (ztv(ig,1)>=(ztv(ig,2))) then
	458	alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.) &
[557]	459	! & *log(1.+zlev(ig,2))
	460	& *sqrt(zlev(ig,2))
	461	! & *sqrt(sqrt(zlev(ig,2)))
[290]	462	! & /sqrt(zlev(ig,2))
[557]	463	! & *zlev(ig,2)
[546]	464	! & *exp(-zlev(ig,2)/1000.)
[161]	465	lalim(ig)=2
	466	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,1)
	467	endif
	468	enddo
	469
[185]	470	do l=2,nlayermx-1
	471	! do l=2,4
	472	do ig=1,ngridmx
[513]	473	if (ztv(ig,l)>(ztv(ig,l+1)+0.) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.)) then ! .and. (zlev(ig,l+1) .lt. 1000.)) then
[161]	474	alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) &
[557]	475	! & *log(1.+zlev(ig,l+1))
	476	& *sqrt(zlev(ig,l+1))
	477	! & *sqrt(sqrt(zlev(ig,l+1)))
[544]	478	! & /sqrt(zlev(ig,l+1))
[557]	479	! & *zlev(ig,l+1)
[546]	480	! & *exp(-zlev(ig,l+1)/1000.)
[161]	481	lalim(ig)=l+1
	482	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	483	endif
	484	enddo
	485	enddo
[185]	486	do l=1,nlayermx
	487	do ig=1,ngridmx
[161]	488	if (alim_star_tot(ig) > 1.e-10 ) then
	489	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
	490	endif
	491	enddo
	492	enddo
	493
	494	alim_star_tot(:)=1.
[185]	495	! if(alim_star(1,1) .ne. 0.) then
	496	! print*, 'lalim star' ,lalim(1)
	497	! endif
[161]	498
	499	!------------------------------------------------------------------------------
	500	! Calcul dans la premiere couche
	501	! On decide dans cette version que le thermique n'est actif que si la premiere
	502	! couche est instable.
	503	! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
	504	! dans une couche l>1
	505	!------------------------------------------------------------------------------
	506
[185]	507	do ig=1,ngridmx
[161]	508	! Le panache va prendre au debut les caracteristiques de l'air contenu
	509	! dans cette couche.
[508]	510	if (activecell(ig)) then
[161]	511	ztla(ig,1)=ztv(ig,1)
	512	!cr: attention, prise en compte de f(1)=1 => AC : what ? f(1) =0. ! (d'ou f(2)=a(1)
	513	! dans un panache conservatif
	514	f_star(ig,1)=0.
[290]	515
[161]	516	f_star(ig,2)=alim_star(ig,1)
[290]	517
[185]	518	zw2(ig,2)=2.g(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) &
[161]	519	& *(zlev(ig,2)-zlev(ig,1)) &
[290]	520	& 0.4pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
[161]	521	w_est(ig,2)=zw2(ig,2)
	522
	523	endif
	524	enddo
	525
	526
	527	!==============================================================================
	528	!boucle de calcul de la vitesse verticale dans le thermique
	529	!==============================================================================
[185]	530	do l=2,nlayermx-1
[161]	531	!==============================================================================
	532
	533
	534	! On decide si le thermique est encore actif ou non
	535	! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
[185]	536	do ig=1,ngridmx
[508]	537	activecell(ig)=activecell(ig) &
[161]	538	& .and. zw2(ig,l)>1.e-10 &
	539	& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
	540	enddo
	541
	542	!---------------------------------------------------------------------------
	543	! calcul des proprietes thermodynamiques et de la vitesse de la couche l
	544	! sans tenir compte du detrainement et de l'entrainement dans cette
	545	! couche
	546	! C'est a dire qu'on suppose
	547	! ztla(l)=ztla(l-1)
	548	! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
	549	! avant) a l'alimentation pour avoir un calcul plus propre
	550	!---------------------------------------------------------------------------
	551
[185]	552	do ig=1,ngridmx
[508]	553	if(activecell(ig)) then
[165]	554	! if(l .lt. lalim(ig)) then
[290]	555	! ztva_est(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
[165]	556	! & alim_star(ig,l)*ztv(ig,l)) &
	557	! & /(f_star(ig,l)+alim_star(ig,l))
	558	! else
[161]	559	ztva_est(ig,l)=ztla(ig,l-1)
[165]	560	! endif
[161]	561
	562	zdz=zlev(ig,l+1)-zlev(ig,l)
[185]	563	zbuoy(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
	564
	565	if (((a1*zbuoy(ig,l)/w_est(ig,l)-b1) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then
	566	w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.zdza1zbuoy(ig,l)-2.zdzw_est(ig,l)b1 &
[532]	567	& -2.(1.-omega)zdzw_est(ig,l)ae(a1zbuoy(ig,l)/w_est(ig,l)-b1)**be)
[161]	568	else
[512]	569	w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.zdza1invzbuoy(ig,l)-2.zdzw_est(ig,l)b1inv)
[161]	570	endif
	571	if (w_est(ig,l+1).lt.0.) then
	572	w_est(ig,l+1)=zw2(ig,l)
	573	endif
	574	endif
	575	enddo
	576
	577	!-------------------------------------------------
	578	!calcul des taux d'entrainement et de detrainement
	579	!-------------------------------------------------
	580
[185]	581	do ig=1,ngridmx
[508]	582	if (activecell(ig)) then
[161]	583
	584	zw2m=w_est(ig,l+1)
[185]	585
	586	if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
[161]	587	entr_star(ig,l)=f_star(ig,l)zdz &
[512]	588	& MAX(0.,ae(a1(zbuoy(ig,l)/zw2m)-b1)**be)
	589	! & MAX(0.,log(1. + 0.03sqrt(a1(zbuoy(ig,l)/zw2m)-b1)))
[161]	590	else
	591	entr_star(ig,l)=0.
	592	endif
[185]	593
[161]	594	if(zbuoy(ig,l) .gt. 0.) then
	595	if(l .lt. lalim(ig)) then
[544]	596
	597	! detr_star(ig,l)=0.
	598	detr_star(ig,l) = f_star(ig,l)zdz &
	599	& adalim
[161]	600	else
[185]	601
[161]	602	! detr_star(ig,l) = f_star(ig,l)zdz &
	603	! & 0.0105((zbuoy(ig,l)/zw2m)/0.048)*(1./1.7)
	604	! detr_star(ig,l) = f_star(ig,l)zdz &
	605	! & 0.0085((zbuoy(ig,l)/zw2m)/0.05)*(1./1.55)
	606
	607	! last baseline from direct les
	608	! detr_star(ig,l) = f_star(ig,l)zdz &
	609	! & 0.065(2.5(zbuoy(ig,l)/zw2m))**0.75
	610
	611	! new param from continuity eq with a fit on dfdz
	612	detr_star(ig,l) = f_star(ig,l)zdz &
[512]	613	& ad
	614	! & Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
[161]	615
[185]	616	! & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) !svn baseline
	617	! & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0008)
	618
[161]	619	! & 0.014((zbuoy(ig,l)/zw2m)/0.05)*(1./1.35)
	620	! detr_star(ig,l) = f_star(ig,l)zdz &
	621	! & ((zbuoy(ig,l)/zw2m)/2.222)! + 0.0002)
	622
	623	endif
	624	else
	625	detr_star(ig,l)=f_star(ig,l)zdz &
[593]	626	& MAX(ad,bd*zbuoy(ig,l)/zw2m)
[512]	627	! & Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
	628	! & Max(0., Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
[161]	629
[508]	630
[185]	631	! & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) !svn baseline
	632
[161]	633	! & 5.(-afactzbetalphazbuoy(ig,l)/zw2m)
	634	! & 5.(-afact*zbuoy(ig,l)/zw2m)
	635
	636	! last baseline from direct les
	637	! & 0.065(-2.5(zbuoy(ig,l)/zw2m))**0.75
	638
	639	! new param from continuity eq with a fit on dfdz
	640
	641
	642	endif
	643
	644	! En dessous de lalim, on prend le max de alim_star et entr_star pour
	645	! alim_star et 0 sinon
	646
[512]	647	if (l.lt.lalim(ig)) then
[161]	648	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
[512]	649	entr_star(ig,l)=0.
	650	endif
[161]	651
	652	! Calcul du flux montant normalise
	653
	654	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	655	& -detr_star(ig,l)
	656
	657	endif
	658	enddo
	659
	660
	661	!----------------------------------------------------------------------------
	662	!calcul de la vitesse verticale en melangeant Tl et qt du thermique
	663	!---------------------------------------------------------------------------
	664
[593]	665	DO tic=0,5 ! internal convergence loop
[508]	666	activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10
[185]	667	do ig=1,ngridmx
[161]	668	if (activetmp(ig)) then
	669
	670	ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ &
	671	& (alim_star(ig,l)+entr_star(ig,l))*ztv(ig,l)) &
	672	& /(f_star(ig,l+1)+detr_star(ig,l))
	673
	674	endif
	675	enddo
	676
[313]	677	activetmp(:)=activetmp(:).and.(abs(ztla(:,l)-ztva(:,l)).gt.0.01)
	678
[185]	679	do ig=1,ngridmx
[161]	680	if (activetmp(ig)) then
	681	ztva(ig,l) = ztla(ig,l)
[185]	682	zbuoy(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
[161]	683
[185]	684	if (((a1*zbuoy(ig,l)/zw2(ig,l)-b1) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then
	685	zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.zdza1*zbuoy(ig,l)- &
[532]	686	& 2.zdzzw2(ig,l)b1-2.(1.-omega)zdzzw2(ig,l)ae(a1zbuoy(ig,l)/zw2(ig,l)-b1)*be)
[161]	687	else
[512]	688	zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.zdza1invzbuoy(ig,l)-2.zdzzw2(ig,l)b1inv)
[161]	689	endif
	690	endif
	691	enddo
	692
[290]	693	! ================ RECOMPUTE ENTR, DETR, and F FROM NEW W2 ===================
	694
	695	do ig=1,ngridmx
	696	if (activetmp(ig)) then
	697
	698	zw2m=zw2(ig,l+1)
	699	if(zw2m .gt. 0) then
	700	if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
	701	entr_star(ig,l)=f_star(ig,l)zdz &
	702	& MAX(0.,ae(a1(zbuoy(ig,l)/zw2m)-b1)**be)
[512]	703	! & MAX(0.,log(1. + 0.03sqrt(a1(zbuoy(ig,l)/zw2m)-b1)))
[290]	704	else
	705	entr_star(ig,l)=0.
	706	endif
	707
	708	if(zbuoy(ig,l) .gt. 0.) then
	709	if(l .lt. lalim(ig)) then
[544]	710
	711	! detr_star(ig,l)=0.
	712	detr_star(ig,l) = f_star(ig,l)zdz &
	713	& adalim
	714
[290]	715	else
	716	detr_star(ig,l) = f_star(ig,l)zdz &
	717	& ad
[512]	718	! & Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
	719
[290]	720	endif
	721	else
	722	detr_star(ig,l)=f_star(ig,l)zdz &
[593]	723	& MAX(ad,bd*zbuoy(ig,l)/zw2m)
[512]	724	! & Max(0.,Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
	725
[290]	726	endif
	727	else
	728	entr_star(ig,l)=0.
	729	detr_star(ig,l)=0.
	730	endif
	731
	732	! En dessous de lalim, on prend le max de alim_star et entr_star pour
	733	! alim_star et 0 sinon
	734
	735	if (l.lt.lalim(ig)) then
	736	alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
	737	entr_star(ig,l)=0.
	738	endif
	739
	740	! Calcul du flux montant normalise
	741
	742	f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &
	743	& -detr_star(ig,l)
	744
	745	endif
	746	enddo
[546]	747
[313]	748	ENDDO ! of tic
	749
[161]	750	!---------------------------------------------------------------------------
	751	!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
	752	!---------------------------------------------------------------------------
	753
[185]	754	do ig=1,ngridmx
[161]	755	if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
	756	print*,'On tombe sur le cas particulier de thermcell_plume'
	757	zw2(ig,l+1)=0.
	758	linter(ig)=l+1
	759	endif
	760
	761	if (zw2(ig,l+1).lt.0.) then
	762	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) &
	763	& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	764	zw2(ig,l+1)=0.
	765	endif
	766	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	767
	768	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	769	! lmix est le niveau de la couche ou w (wa_moy) est maximum
	770	wmaxa(ig)=wa_moy(ig,l+1)
	771	endif
	772	enddo
	773
	774	!=========================================================================
	775	! FIN DE LA BOUCLE VERTICALE
	776	enddo
	777	!=========================================================================
	778
	779	!on recalcule alim_star_tot
[185]	780	do ig=1,ngridmx
[161]	781	alim_star_tot(ig)=0.
	782	enddo
[185]	783	do ig=1,ngridmx
[161]	784	do l=1,lalim(ig)-1
	785	alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
	786	enddo
	787	enddo
	788
[185]	789	do l=1,nlayermx
	790	do ig=1,ngridmx
[165]	791	if (alim_star_tot(ig) > 1.e-10 ) then
	792	alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
	793	endif
	794	enddo
	795	enddo
[161]	796
	797	! ===========================================================================
	798	! ================= FIN PLUME ===============================================
	799	! ===========================================================================
	800
	801	! ===========================================================================
	802	! ================= HEIGHT ==================================================
	803	! ===========================================================================
	804
	805	! Attention, w2 est transforme en sa racine carree dans cette routine
	806
	807	!-------------------------------------------------------------------------------
[356]	808	! Calcul des caracteristiques du thermique:zmax,wmax
[161]	809	!-------------------------------------------------------------------------------
	810
	811	!calcul de la hauteur max du thermique
[185]	812	do ig=1,ngridmx
[161]	813	lmax(ig)=lalim(ig)
	814	enddo
[185]	815	do ig=1,ngridmx
	816	do l=nlayermx,lalim(ig)+1,-1
[161]	817	if (zw2(ig,l).le.1.e-10) then
[512]	818	lmax(ig)=l-1
[161]	819	endif
	820	enddo
	821	enddo
	822
	823	! On traite le cas particulier qu'il faudrait éviter ou le thermique
	824	! atteind le haut du modele ...
[185]	825	do ig=1,ngridmx
	826	if ( zw2(ig,nlayermx) > 1.e-10 ) then
[161]	827	print*,'WARNING !!!!! W2 thermiques non nul derniere couche '
[185]	828	lmax(ig)=nlayermx
[161]	829	endif
	830	enddo
	831
	832	! pas de thermique si couche 1 stable
[185]	833	! do ig=1,ngridmx
[161]	834	! if (lmin(ig).gt.1) then
	835	! lmax(ig)=1
	836	! lmin(ig)=1
	837	! lalim(ig)=1
	838	! endif
	839	! enddo
	840	!
	841	! Determination de zw2 max
[185]	842	do ig=1,ngridmx
[161]	843	wmax(ig)=0.
	844	enddo
	845
[185]	846	do l=1,nlayermx
	847	do ig=1,ngridmx
[161]	848	if (l.le.lmax(ig)) then
	849	if (zw2(ig,l).lt.0.)then
	850	print*,'pb2 zw2<0'
	851	endif
	852	zw2(ig,l)=sqrt(zw2(ig,l))
	853	wmax(ig)=max(wmax(ig),zw2(ig,l))
	854	else
	855	zw2(ig,l)=0.
	856	endif
	857	enddo
	858	enddo
	859	! Longueur caracteristique correspondant a la hauteur des thermiques.
[185]	860	do ig=1,ngridmx
[161]	861	zmax(ig)=0.
	862	zlevinter(ig)=zlev(ig,1)
	863	enddo
	864
	865	num(:)=0.
	866	denom(:)=0.
[185]	867	do ig=1,ngridmx
	868	do l=1,nlayermx
[161]	869	num(ig)=num(ig)+zw2(ig,l)zlev(ig,l)(zlev(ig,l+1)-zlev(ig,l))
	870	denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
	871	enddo
	872	enddo
[185]	873	do ig=1,ngridmx
[161]	874	if (denom(ig).gt.1.e-10) then
	875	zmax(ig)=2.*num(ig)/denom(ig)
	876	endif
	877	enddo
	878
	879	! Attention, w2 est transforme en sa racine carree dans cette routine
	880
	881	! ===========================================================================
	882	! ================= FIN HEIGHT ==============================================
	883	! ===========================================================================
	884
	885	! Choix de la fonction d'alimentation utilisee pour la fermeture.
	886
	887	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
	888
	889	! ===========================================================================
	890	! ============= CLOSURE =====================================================
	891	! ===========================================================================
	892
	893	!-------------------------------------------------------------------------------
	894	! Fermeture,determination de f
	895	!-------------------------------------------------------------------------------
	896	! Appel avec la version seche
	897
	898	alim_star2(:)=0.
	899	alim_star_tot_clos(:)=0.
	900	f(:)=0.
	901
	902	! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
	903	llmax=1
[185]	904	do ig=1,ngridmx
[161]	905	if (lalim(ig)>llmax) llmax=lalim(ig)
	906	enddo
	907
	908
	909	! Calcul des integrales sur la verticale de alim_star et de
	910	! alim_star^2/(rho dz)
	911	do k=1,llmax-1
[185]	912	do ig=1,ngridmx
[161]	913	if (k<lalim(ig)) then
[185]	914	alim_star2(ig)=alim_star2(ig)+alim_star_clos(ig,k)*alim_star_clos(ig,k) &
[161]	915	& /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
	916	alim_star_tot_clos(ig)=alim_star_tot_clos(ig)+alim_star_clos(ig,k)
	917	endif
	918	enddo
	919	enddo
[185]	920
[161]	921	! WARNING : MARS MODIF : we have added 2. : ratio of wmax/vmoy
	922	! True ratio is 3.5 but wetake into account the vmoy is the one alimentating
	923	! the thermal, so there are vs=0 into the vmoy... the true vmoy is lower. (a la louche)
	924	! And r_aspect has been changed from 2 to 1.5 from observations
[185]	925	do ig=1,ngridmx
[161]	926	if (alim_star2(ig)>1.e-10) then
[185]	927	! f(ig)=wmax_sec(ig)*alim_star_tot_clos(ig)/ &
	928	! & (max(500.,zmax_sec(ig))r_aspectalim_star2(ig))
	929	f(ig)=wmax(ig)*alim_star_tot_clos(ig)/ &
	930	& (max(500.,zmax(ig))r_aspectalim_star2(ig))
	931
[161]	932	endif
	933	enddo
	934
	935	! ===========================================================================
	936	! ============= FIN CLOSURE =================================================
	937	! ===========================================================================
	938
	939
	940	! ===========================================================================
	941	! ============= FLUX2 =======================================================
	942	! ===========================================================================
	943
	944	!-------------------------------------------------------------------------------
	945	!deduction des flux
	946	!-------------------------------------------------------------------------------
	947
	948	fomass_max=0.8
	949	alphamax=0.5
	950
	951	ncorecfm1=0
	952	ncorecfm2=0
	953	ncorecfm3=0
	954	ncorecfm4=0
	955	ncorecfm5=0
	956	ncorecfm6=0
	957	ncorecfm7=0
	958	ncorecfm8=0
	959	ncorecalpha=0
	960
	961	!-------------------------------------------------------------------------
	962	! Multiplication par le flux de masse issu de la femreture
	963	!-------------------------------------------------------------------------
	964
[185]	965	do l=1,nlayermx
[161]	966	entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
	967	detr(:,l)=f(:)*detr_star(:,l)
	968	enddo
	969
[185]	970	do l=1,nlayermx
	971	do ig=1,ngridmx
[161]	972	if (l.lt.lmax(ig)) then
	973	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
	974	elseif(l.eq.lmax(ig)) then
	975	fm(ig,l+1)=0.
	976	detr(ig,l)=fm(ig,l)+entr(ig,l)
	977	else
	978	fm(ig,l+1)=0.
	979	endif
	980	enddo
	981	enddo
	982
	983	! Test provisoire : pour comprendre pourquoi on corrige plein de fois
	984	! le cas fm6, on commence par regarder une premiere fois avant les
	985	! autres corrections.
	986
[314]	987	! do l=1,nlayermx
	988	! do ig=1,ngridmx
	989	! if (detr(ig,l).gt.fm(ig,l)) then
	990	! ncorecfm8=ncorecfm8+1
	991	! endif
	992	! enddo
	993	! enddo
[161]	994
	995	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	996	! FH Version en cours de test;
	997	! par rapport a thermcell_flux, on fait une grande boucle sur "l"
	998	! et on modifie le flux avec tous les contrï¿½les appliques d'affilee
	999	! pour la meme couche
	1000	! Momentanement, on duplique le calcule du flux pour pouvoir comparer
	1001	! les flux avant et apres modif
	1002	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1003
[185]	1004	do l=1,nlayermx
[161]	1005
[185]	1006	do ig=1,ngridmx
[161]	1007	if (l.lt.lmax(ig)) then
	1008	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
	1009	elseif(l.eq.lmax(ig)) then
	1010	fm(ig,l+1)=0.
	1011	detr(ig,l)=fm(ig,l)+entr(ig,l)
	1012	else
	1013	fm(ig,l+1)=0.
	1014	endif
	1015	enddo
	1016
	1017
	1018	!-------------------------------------------------------------------------
	1019	! Verification de la positivite des flux de masse
	1020	!-------------------------------------------------------------------------
	1021
[185]	1022	do ig=1,ngridmx
[190]	1023
[161]	1024	if (fm(ig,l+1).lt.0.) then
[190]	1025	if((l+1) .eq. lmax(ig)) then
	1026	detr(ig,l)=detr(ig,l)+fm(ig,l+1)
	1027	fm(ig,l+1)=0.
[336]	1028	entr(ig,l+1)=0.
[190]	1029	ncorecfm2=ncorecfm2+1
	1030	else
	1031	print*,'fm(l+1)<0 : ig, l+1,lmax :',ig,l+1,lmax(ig),fm(ig,l+1)
	1032	ncorecfm1=ncorecfm1+1
[161]	1033	fm(ig,l+1)=fm(ig,l)
	1034	detr(ig,l)=entr(ig,l)
[190]	1035	endif
[161]	1036	endif
[190]	1037
[161]	1038	enddo
	1039
[508]	1040	! Les "optimisations" du flux sont desactivecelles : moins de bidouilles
[161]	1041	! je considere que celles ci ne sont pas justifiees ou trop delicates
	1042	! pour MARS, d'apres des observations LES.
	1043	!-------------------------------------------------------------------------
	1044	!Test sur fraca croissant
	1045	!-------------------------------------------------------------------------
	1046	! if (iflag_thermals_optflux==0) then
[185]	1047	! do ig=1,ngridmx
[161]	1048	! if (l.ge.lalim(ig).and.l.le.lmax(ig) &
	1049	! & .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
	1050	!! zzz est le flux en l+1 a frac constant
	1051	! zzz=fm(ig,l)rhobarz(ig,l+1)zw2(ig,l+1) &
	1052	! & /(rhobarz(ig,l)*zw2(ig,l))
	1053	! if (fm(ig,l+1).gt.zzz) then
	1054	! detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz
	1055	! fm(ig,l+1)=zzz
	1056	! ncorecfm4=ncorecfm4+1
	1057	! endif
	1058	! endif
	1059	! enddo
	1060	! endif
	1061	!
	1062	!-------------------------------------------------------------------------
	1063	!test sur flux de masse croissant
	1064	!-------------------------------------------------------------------------
	1065	! if (iflag_thermals_optflux==0) then
[185]	1066	! do ig=1,ngridmx
[161]	1067	! if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
	1068	! f_old=fm(ig,l+1)
	1069	! fm(ig,l+1)=fm(ig,l)
	1070	! detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
	1071	! ncorecfm5=ncorecfm5+1
	1072	! endif
	1073	! enddo
	1074	! endif
	1075	!
	1076	!-------------------------------------------------------------------------
	1077	!detr ne peut pas etre superieur a fm
	1078	!-------------------------------------------------------------------------
	1079
[185]	1080	do ig=1,ngridmx
[161]	1081	if (detr(ig,l).gt.fm(ig,l)) then
	1082	ncorecfm6=ncorecfm6+1
	1083	detr(ig,l)=fm(ig,l)
	1084	entr(ig,l)=fm(ig,l+1)
	1085
	1086	! Dans le cas ou on est au dessus de la couche d'alimentation et que le
	1087	! detrainement est plus fort que le flux de masse, on stope le thermique.
[314]	1088	! endif
[161]	1089
[336]	1090	if(l.gt.lmax(ig)) then
	1091	! if(l.gt.lalim(ig)) then
[161]	1092	detr(ig,l)=0.
	1093	fm(ig,l+1)=0.
	1094	entr(ig,l)=0.
	1095	endif
[314]	1096
	1097	endif
	1098
[161]	1099	enddo
	1100
	1101	!-------------------------------------------------------------------------
	1102	!fm ne peut pas etre negatif
	1103	!-------------------------------------------------------------------------
	1104
[185]	1105	do ig=1,ngridmx
[161]	1106	if (fm(ig,l+1).lt.0.) then
	1107	detr(ig,l)=detr(ig,l)+fm(ig,l+1)
	1108	fm(ig,l+1)=0.
	1109	ncorecfm2=ncorecfm2+1
	1110	endif
	1111	enddo
	1112
	1113	!-----------------------------------------------------------------------
	1114	!la fraction couverte ne peut pas etre superieure a 1
	1115	!-----------------------------------------------------------------------
	1116	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1117	! FH Partie a revisiter.
	1118	! Il semble qu'etaient codees ici deux optiques dans le cas
	1119	! F/ (rho *w) > 1
	1120	! soit limiter la hauteur du thermique en considerant que c'est
	1121	! la derniere chouche, soit limiter F a rho w.
	1122	! Dans le second cas, il faut en fait limiter a un peu moins
	1123	! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
	1124	! dans thermcell_main et qu'il semble de toutes facons deraisonable
	1125	! d'avoir des fractions de 1..
	1126	! Ci dessous, et dans l'etat actuel, le premier des deux if est
	1127	! sans doute inutile.
	1128	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1129
[185]	1130	do ig=1,ngridmx
[161]	1131	if (zw2(ig,l+1).gt.1.e-10) then
	1132	zfm=rhobarz(ig,l+1)zw2(ig,l+1)alphamax
	1133	if ( fm(ig,l+1) .gt. zfm) then
	1134	f_old=fm(ig,l+1)
	1135	fm(ig,l+1)=zfm
	1136	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
	1137	ncorecalpha=ncorecalpha+1
	1138	endif
	1139	endif
	1140
	1141	enddo
	1142
	1143	! Fin de la grande boucle sur les niveaux verticaux
	1144	enddo
	1145
	1146	!-----------------------------------------------------------------------
	1147	! On fait en sorte que la quantite totale d'air entraine dans le
	1148	! panache ne soit pas trop grande comparee a la masse de la maille
	1149	!-----------------------------------------------------------------------
	1150
[185]	1151	do l=1,nlayermx-1
	1152	do ig=1,ngridmx
[161]	1153	eee0=entr(ig,l)
	1154	ddd0=detr(ig,l)
	1155	eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
	1156	ddd=detr(ig,l)-eee
	1157	if (eee.gt.0.) then
	1158	ncorecfm3=ncorecfm3+1
	1159	entr(ig,l)=entr(ig,l)-eee
	1160	if ( ddd.gt.0.) then
	1161	! l'entrainement est trop fort mais l'exces peut etre compense par une
	1162	! diminution du detrainement)
	1163	detr(ig,l)=ddd
	1164	else
	1165	! l'entrainement est trop fort mais l'exces doit etre compense en partie
	1166	! par un entrainement plus fort dans la couche superieure
	1167	if(l.eq.lmax(ig)) then
	1168	detr(ig,l)=fm(ig,l)+entr(ig,l)
	1169	else
	1170	entr(ig,l+1)=entr(ig,l+1)-ddd
	1171	detr(ig,l)=0.
	1172	fm(ig,l+1)=fm(ig,l)+entr(ig,l)
	1173	detr(ig,l)=0.
	1174	endif
	1175	endif
	1176	endif
	1177	enddo
	1178	enddo
	1179	!
	1180	! ddd=detr(ig)-entre
	1181	!on s assure que tout s annule bien en zmax
[185]	1182	do ig=1,ngridmx
[161]	1183	fm(ig,lmax(ig)+1)=0.
	1184	entr(ig,lmax(ig))=0.
	1185	detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig))
	1186	enddo
	1187
	1188	!-----------------------------------------------------------------------
	1189	! Impression du nombre de bidouilles qui ont ete necessaires
	1190	!-----------------------------------------------------------------------
	1191
	1192	!IM 090508 beg
[185]	1193	if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > ngridmx/4. ) then
[161]	1194	print*,'thermcell warning : large number of corrections'
	1195	print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
	1196	& ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
	1197	& ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
	1198	& ncorecfm6,'x fm6', &
	1199	& ncorecfm7,'x fm7', &
	1200	& ncorecfm8,'x fm8', &
	1201	& ncorecalpha,'x alpha'
	1202	endif
	1203
	1204	! ===========================================================================
	1205	! ============= FIN FLUX2 ===================================================
	1206	! ===========================================================================
	1207
	1208
	1209	! ===========================================================================
	1210	! ============= TRANSPORT ===================================================
	1211	! ===========================================================================
	1212
	1213	!------------------------------------------------------------------
	1214	! calcul du transport vertical
	1215	!------------------------------------------------------------------
	1216
	1217	! ------------------------------------------------------------------
	1218	! Transport de teta dans l'updraft : (remplace thermcell_dq)
	1219	! ------------------------------------------------------------------
	1220
	1221	zdthladj(:,:)=0.
	1222
[185]	1223	do ig=1,ngridmx
[161]	1224	if(lmax(ig) .gt. 1) then
	1225	do k=1,lmax(ig)
	1226	zdthladj(ig,k)=(1./masse(ig,k))(fm(ig,k+1)ztv(ig,k+1)- &
	1227	& fm(ig,k)ztv(ig,k)+fm(ig,k)ztva(ig,k)-fm(ig,k+1)*ztva(ig,k+1))
	1228	if (ztv(ig,k) + ptimestep*zdthladj(ig,k) .le. 0.) then
[165]	1229	print,'Teta<0 in thermcell_dTeta up: qenv .. dq : ', ztv(ig,k),ptimestepzdthladj(ig,k)
	1230	if(ztv(ig,k) .gt. 0.) then
	1231	zdthladj(ig,k)=0.
	1232	endif
[161]	1233	endif
	1234	enddo
	1235	endif
	1236	enddo
	1237
	1238	! ------------------------------------------------------------------
	1239	! Prescription des proprietes du downdraft
	1240	! ------------------------------------------------------------------
	1241
	1242	ztvd(:,:)=ztv(:,:)
	1243	fm_down(:,:)=0.
[185]	1244	do ig=1,ngridmx
[161]	1245	if (lmax(ig) .gt. 1) then
	1246	do l=1,lmax(ig)
[512]	1247	! if(zlay(ig,l) .le. 0.8*zmax(ig)) then
	1248	if(zlay(ig,l) .le. zmax(ig)) then
[496]	1249	fm_down(ig,l) =fm(ig,l)* &
	1250	& max(fdfu,-3*max(0.,(zlay(ig,l)/zmax(ig)))-0.9)
[161]	1251	endif
	1252
[512]	1253	! if(zlay(ig,l) .le. 0.06*zmax(ig)) then
	1254	! ztvd(ig,l)=ztv(ig,l)max(0.,(1.+(sqrt((zlay(ig,l)/zmax(ig))/0.122449) - 1.)(ztva(ig,l)/ztv(ig,l) - 1.)))
	1255	! elseif(zlay(ig,l) .le. 0.4*zmax(ig)) then
	1256	! ztvd(ig,l)=ztv(ig,l)max(0.,1.-0.25(ztva(ig,l)/ztv(ig,l) - 1.))
	1257	! elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
	1258	! ztvd(ig,l)=ztv(ig,l)max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/1.)(ztva(ig,l)/ztv(ig,l) - 1.)))
	1259	! else
	1260	! ztvd(ig,l)=ztv(ig,l)
	1261	! endif
	1262
	1263	! if(zlay(ig,l) .le. 0.6*zmax(ig)) then
	1264	! ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/179.848 + 0.997832)
	1265	! elseif(zlay(ig,l) .le. 0.8*zmax(ig)) then
	1266	! ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-171.74)/170.94)
	1267	! else
	1268	! ztvd(ig,l)=ztv(ig,l)
	1269	! endif
	1270
	1271
	1272	! if(zlay(ig,l) .le. 0.8*zmax(ig)) then
	1273	! ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/224.81 + 0.997832)
	1274	! elseif(zlay(ig,l) .le. zmax(ig)) then
	1275	! ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-144.685)/143.885)
	1276	! else
	1277	! ztvd(ig,l)=ztv(ig,l)
	1278	! endif
	1279
	1280
[546]	1281	! if (zbuoy(ig,l) .gt. 0.) then
	1282	! ztvd(ig,l)=ztva(ig,l)*0.9998
	1283	!! ztvd(ig,l)=ztv(ig,l)*0.997832
	1284	!! else
	1285	!! if(zlay(ig,l) .le. zmax(ig)) then
	1286	!! ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/299.7 + 0.997832)
	1287	!! endif
	1288	! endif
	1289
	1290	if(zlay(ig,l) .le. zmax(ig)) then
	1291	ztvd(ig,l)=min(ztv(ig,l),ztv(ig,l)*((zlay(ig,l)/zmax(ig))/400. + 0.997832))
	1292	! ztvd(ig,l)=min(ztv(ig,l),ztv(ig,l)*((zlay(ig,l)/zmax(ig))/299.7 + 0.997832))
[161]	1293	else
	1294	ztvd(ig,l)=ztv(ig,l)
	1295	endif
	1296
	1297	enddo
	1298	endif
	1299	enddo
	1300
	1301	! ------------------------------------------------------------------
	1302	! Transport de teta dans le downdraft : (remplace thermcell_dq)
	1303	! ------------------------------------------------------------------
	1304
	1305	zdthladj_down(:,:)=0.
	1306
[185]	1307	do ig=1,ngridmx
[161]	1308	if(lmax(ig) .gt. 1) then
[290]	1309	! No downdraft in the very-near surface layer, we begin at k=3
[496]	1310
	1311	do k=3,lmax(ig)
[161]	1312	zdthladj_down(ig,k)=(1./masse(ig,k))(fm_down(ig,k+1)ztv(ig,k+1)- &
	1313	& fm_down(ig,k)ztv(ig,k)+fm_down(ig,k)ztvd(ig,k)-fm_down(ig,k+1)*ztvd(ig,k+1))
	1314	if (ztv(ig,k) + ptimestep*zdthladj_down(ig,k) .le. 0.) then
	1315	print,'q<0 in thermcell_dTeta down: qenv .. dq : ', ztv(ig,k),ptimestepzdthladj_down(ig,k)
[165]	1316	if(ztv(ig,k) .gt. 0.) then
	1317	zdthladj(ig,k)=0.
	1318	endif
[161]	1319	endif
	1320	enddo
	1321	endif
	1322	enddo
	1323
	1324	!------------------------------------------------------------------
	1325	! Calcul de la fraction de l'ascendance
	1326	!------------------------------------------------------------------
[185]	1327	do ig=1,ngridmx
[161]	1328	fraca(ig,1)=0.
[185]	1329	fraca(ig,nlayermx+1)=0.
[161]	1330	enddo
[185]	1331	do l=2,nlayermx
	1332	do ig=1,ngridmx
[161]	1333	if (zw2(ig,l).gt.1.e-10) then
	1334	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
	1335	else
	1336	fraca(ig,l)=0.
	1337	endif
	1338	enddo
	1339	enddo
	1340
	1341
	1342
	1343	! ===========================================================================
	1344	! ============= DV2 =========================================================
	1345	! ===========================================================================
	1346	! ROUTINE OVERIDE : ne prends pas en compte le downdraft
	1347	! de plus, le gradient de pression horizontal semble tout deregler... A VOIR
	1348
	1349	if (0 .eq. 1) then
	1350
	1351	!------------------------------------------------------------------
	1352	! calcul du transport vertical du moment horizontal
	1353	!------------------------------------------------------------------
	1354
	1355	! Calcul du transport de V tenant compte d'echange par gradient
	1356	! de pression horizontal avec l'environnement
	1357
	1358	! calcul du detrainement
	1359	!---------------------------
	1360
	1361	nlarga0=0.
	1362
[185]	1363	do k=1,nlayermx
	1364	do ig=1,ngridmx
[161]	1365	detr_dv2(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
	1366	enddo
	1367	enddo
	1368
	1369	! calcul de la valeur dans les ascendances
[185]	1370	do ig=1,ngridmx
[161]	1371	zua(ig,1)=zu(ig,1)
	1372	zva(ig,1)=zv(ig,1)
	1373	ue(ig,1)=zu(ig,1)
	1374	ve(ig,1)=zv(ig,1)
	1375	enddo
	1376
[185]	1377	gamma(1:ngridmx,1)=0.
	1378	do k=2,nlayermx
	1379	do ig=1,ngridmx
[161]	1380	ltherm(ig,k)=(fm(ig,k+1)+detr_dv2(ig,k))ptimestep > 1.e-5masse(ig,k)
	1381	if(ltherm(ig,k).and.zmax(ig)>0.) then
	1382	gamma0(ig,k)=masse(ig,k) &
	1383	& sqrt( 0.5(fraca(ig,k+1)+fraca(ig,k)) ) &
	1384	& *0.5/zmax(ig) &
	1385	& *1.
	1386	else
	1387	gamma0(ig,k)=0.
	1388	endif
	1389	if (ltherm(ig,k).and.zmax(ig)<=0.) nlarga0=nlarga0+1
	1390	enddo
	1391	enddo
	1392
	1393	gamma(:,:)=0.
	1394
[185]	1395	do k=2,nlayermx
[161]	1396
[185]	1397	do ig=1,ngridmx
[161]	1398
	1399	if (ltherm(ig,k)) then
	1400	dua(ig,k)=zua(ig,k-1)-zu(ig,k-1)
	1401	dva(ig,k)=zva(ig,k-1)-zv(ig,k-1)
	1402	else
	1403	zua(ig,k)=zu(ig,k)
	1404	zva(ig,k)=zv(ig,k)
	1405	ue(ig,k)=zu(ig,k)
	1406	ve(ig,k)=zv(ig,k)
	1407	endif
	1408	enddo
	1409
	1410
	1411	! Debut des iterations
	1412	!----------------------
	1413
	1414	! AC WARNING : SALE !
	1415
	1416	do iter=1,5
[185]	1417	do ig=1,ngridmx
[161]	1418	! Pour memoire : calcul prenant en compte la fraction reelle
	1419	! zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
	1420	! zf2=1./(1.-zf)
	1421	! Calcul avec fraction infiniement petite
	1422	zf=0.
	1423	zf2=1.
	1424
	1425	! la première fois on multiplie le coefficient de freinage
	1426	! par le module du vent dans la couche en dessous.
	1427	! Mais pourquoi donc ???
	1428	if (ltherm(ig,k)) then
	1429	! On choisit une relaxation lineaire.
	1430	! gamma(ig,k)=gamma0(ig,k)
	1431	! On choisit une relaxation quadratique.
[337]	1432	gamma(ig,k)=gamma0(ig,k)sqrt(dua(ig,k)2+dva(ig,k)*2)
[161]	1433	zua(ig,k)=(fm(ig,k)*zua(ig,k-1) &
	1434	& +(zf2entr(ig,k)+gamma(ig,k))zu(ig,k)) &
	1435	& /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)zfzf2 &
	1436	& +gamma(ig,k))
	1437	zva(ig,k)=(fm(ig,k)*zva(ig,k-1) &
	1438	& +(zf2entr(ig,k)+gamma(ig,k))zv(ig,k)) &
	1439	& /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)zfzf2 &
	1440	& +gamma(ig,k))
	1441
	1442	! print*,' OUTPUT DV2 !!!!!!!!!!!!!!!!!!!!',k,zua(ig,k),zva(ig,k),zu(ig,k),zv(ig,k),dua(ig,k),dva(ig,k)
	1443	dua(ig,k)=zua(ig,k)-zu(ig,k)
	1444	dva(ig,k)=zva(ig,k)-zv(ig,k)
	1445	ue(ig,k)=(zu(ig,k)-zfzua(ig,k))zf2
	1446	ve(ig,k)=(zv(ig,k)-zfzva(ig,k))zf2
	1447	endif
	1448	enddo
	1449	! Fin des iterations
	1450	!--------------------
	1451	enddo
	1452
[185]	1453	enddo ! k=2,nlayermx
[161]	1454
	1455	! Calcul du flux vertical de moment dans l'environnement.
	1456	!---------------------------------------------------------
[185]	1457	do k=2,nlayermx
	1458	do ig=1,ngridmx
[161]	1459	wud(ig,k)=fm(ig,k)*ue(ig,k)
	1460	wvd(ig,k)=fm(ig,k)*ve(ig,k)
	1461	enddo
	1462	enddo
[185]	1463	do ig=1,ngridmx
[161]	1464	wud(ig,1)=0.
[185]	1465	wud(ig,nlayermx+1)=0.
[161]	1466	wvd(ig,1)=0.
[185]	1467	wvd(ig,nlayermx+1)=0.
[161]	1468	enddo
	1469
	1470	! calcul des tendances.
	1471	!-----------------------
[185]	1472	do k=1,nlayermx
	1473	do ig=1,ngridmx
[161]	1474	pduadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zua(ig,k) &
	1475	& -(entr(ig,k)+gamma(ig,k))*ue(ig,k) &
	1476	& -wud(ig,k)+wud(ig,k+1)) &
	1477	& /masse(ig,k)
	1478	pdvadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zva(ig,k) &
	1479	& -(entr(ig,k)+gamma(ig,k))*ve(ig,k) &
	1480	& -wvd(ig,k)+wvd(ig,k+1)) &
	1481	& /masse(ig,k)
	1482	enddo
	1483	enddo
	1484
	1485
	1486	! Sorties eventuelles.
	1487	!----------------------
	1488
	1489	! if (nlarga0>0) then
	1490	! print*,'WARNING !!!!!! DANS THERMCELL_DV2 '
	1491	! print*,nlarga0,' points pour lesquels laraga=0. dans un thermique'
	1492	! print*,'Il faudrait decortiquer ces points'
	1493	! endif
	1494
	1495	! ===========================================================================
	1496	! ============= FIN DV2 =====================================================
	1497	! ===========================================================================
	1498
	1499	else
	1500
[337]	1501	! modname='momentum'
	1502	! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, &
	1503	! & masse,zu,pduadj,ztvd,fm_down,ztv,modname,lmax)
	1504	!
	1505	! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, &
	1506	! & masse,zv,pdvadj,ztvd,fm_down,ztv,modname,lmax)
[161]	1507
	1508	endif
	1509
	1510	!------------------------------------------------------------------
[508]	1511	! calcul du transport vertical de traceurs
[161]	1512	!------------------------------------------------------------------
	1513
[508]	1514	! We only transport co2 tracer because it is coupled to the scheme through theta_m
	1515	! The rest is transported outside the sub-timestep loop
	1516
	1517	if (ico2.ne.0) then
	1518	! if (nqmx .ne. 0.) then
[185]	1519	do l=1,nlayermx
[508]	1520	zdzfull(:,l)=zlev(:,l+1)-zlev(:,l)
[161]	1521	enddo
	1522
[508]	1523	modname='tracer'
	1524	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	1525	& ,fm,entr,detr, &
	1526	& masse,pq(:,:,ico2),pdqadj(:,:,ico2),modname,zdzfull)
	1527	! endif
	1528
	1529	! Compute the ratio between theta and theta_m
	1530
	1531	do l=1,nlayermx
	1532	do ig=1,ngridmx
	1533	ratiom(ig,l)=1./(A(pq(ig,l,ico2)+pdqadj(ig,l,ico2)ptimestep)+B)
	1534	enddo
	1535	enddo
	1536	else
	1537	ratiom(:,:)=1.
	1538	endif
	1539
	1540
[161]	1541	!------------------------------------------------------------------
[508]	1542	! incrementation dt
[161]	1543	!------------------------------------------------------------------
	1544
[508]	1545	do l=1,nlayermx
	1546	do ig=1,ngridmx
	1547	pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))zpopsk(ig,l)ratiom(ig,l)
[512]	1548	! pdtadj(ig,l)=zdthladj(ig,l)zpopsk(ig,l)ratiom(ig,l)
[508]	1549	enddo
	1550	enddo
[161]	1551
	1552	!------------------------------------------------------------------
	1553	! calcul du transport vertical de la tke
	1554	!------------------------------------------------------------------
	1555
[313]	1556	! modname='tke'
	1557	! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, &
	1558	! & masse,pq2,pdq2adj,ztvd,fm_down,ztv,modname,lmax)
[161]	1559
	1560	! ===========================================================================
	1561	! ============= FIN TRANSPORT ===============================================
	1562	! ===========================================================================
	1563
	1564
	1565	!------------------------------------------------------------------
	1566	! Calculs de diagnostiques pour les sorties
	1567	!------------------------------------------------------------------
	1568	! DIAGNOSTIQUE
	1569	! We compute interface values for teta env and th. The last interface
	1570	! value does not matter, as the mass flux is 0 there.
	1571
	1572
[185]	1573	do l=1,nlayermx-1
	1574	do ig=1,ngridmx
[508]	1575	teta_th_int(ig,l)=0.5(ztva(ig,l+1)+ztva(ig,l))ratiom(ig,l)
	1576	teta_down_int(ig,l) = 0.5(ztvd(ig,l+1)+ztvd(ig,l))ratiom(ig,l)
	1577	teta_env_int(ig,l)=0.5(ztv(ig,l+1)+ztv(ig,l))ratiom(ig,l)
[161]	1578	enddo
	1579	enddo
[185]	1580	do ig=1,ngridmx
	1581	teta_th_int(ig,nlayermx)=teta_th_int(ig,nlayermx-1)
	1582	teta_env_int(ig,nlayermx)=teta_env_int(ig,nlayermx-1)
	1583	teta_down_int(ig,nlayermx)=teta_down_int(ig,nlayermx-1)
[161]	1584	enddo
[185]	1585	do l=1,nlayermx
	1586	do ig=1,ngridmx
[161]	1587	heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))
[508]	1588	buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	1589	buoyancyEst(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
[161]	1590	heatFlux_down(ig,l)=fm_down(ig,l)*(teta_down_int(ig,l)-teta_env_int(ig,l))/rhobarz(ig,l)
	1591	enddo
	1592	enddo
	1593
	1594	return
	1595	end

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

Download in other formats:

Original Format