Context Navigation

thermcell_main_mars.F90 @ 1009

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

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

Download in other formats:

Original Format