Context Navigation

thermcell_main_mars.F90 @ 501

Last change on this file since 501 was 496, checked in by acolaitis, 14 years ago

M 489 libf/phymars/thermcell_main_mars.F90
^{--------------------> New parameters for thermals, following latest version of the relevant article}

M 489 libf/phymars/physiq.F
^{--------------------> Minor modifications}

D 489 libf/phymars/surflayer_interpol.F
A 0 libf/phymars/pbl_parameters.F
^{--------------------> Replaced surflayer_interpol.F by a cleaner and richer version called pbl_parameters.F}

This routine is the base of what will later be implemented in the MCD and as a tool to
compute surface properties from output fields (so that it will ultimately disappear from libf
and might become an utility)

M 489 libf/phymars/vdif_cd.F
^{--------------------> Minor modification to the Richardson computation}

M 489 libf/phymars/vdifc.F
^{--------------------> Minor modification to the aerodynamic conductances computation, replacing wmax by hfmax, which}

is a better proxy for convective activity. Might be replaced by w_star later.

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

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

Download in other formats:

Original Format