Context Navigation

thermcell_main_mars.F90 @ 310

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

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

Download in other formats:

Original Format