Context Navigation

thermcell_main_mars.F90 @ 314

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

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

Download in other formats:

Original Format