Context Navigation

thermcell_main_mars.F90 @ 627

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

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

Download in other formats:

Original Format