Context Navigation

thermcell_main_mars.F90 @ 426

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

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

Download in other formats:

Original Format