Context Navigation

thermcell_main_mars.F90 @ 173

Last change on this file since 173 was 173, checked in by acolaitis, 14 years ago
Minor modifications related to thermals
File size: 50.5 KB

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

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

Download in other formats:

Original Format