Context Navigation

thermcell_main_mars.F90 @ 222

Last change on this file since 222 was 190, checked in by acolaitis, 14 years ago

04/07/2011 == AC

Minor setting modification to thermals
Added new flux optimization in thermcell_main_mars.F90, to run properly

in 3D

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

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

Download in other formats:

Original Format