Context Navigation

thermcell_main.F90 @ 354

Last change on this file since 354 was 179, checked in by lfita, 10 years ago
No Ale, Alp, ... values because pphi is at the intermediate level!!! Removing raw prints of checkig for Ale, Alp
File size: 42.0 KB

Rev	Line
[1]	1	!
	2	! $Id: thermcell_main.F90 1795 2013-07-18 08:20:28Z emillour $
	3	!
	4	SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep &
	5	& ,pplay,pplev,pphi,debut &
	6	& ,pu,pv,pt,po &
	7	& ,pduadj,pdvadj,pdtadj,pdoadj &
	8	& ,fm0,entr0,detr0,zqta,zqla,lmax &
	9	& ,ratqscth,ratqsdiff,zqsatth &
	10	& ,Ale_bl,Alp_bl,lalim_conv,wght_th &
	11	& ,zmax0, f0,zw2,fraca,ztv &
	12	& ,zpspsk,ztla,zthl &
	13	!!! nrlmd le 10/04/2012
	14	& ,pbl_tke,pctsrf,omega,airephy &
	15	& ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
	16	& ,n2,s2,ale_bl_stat &
	17	& ,therm_tke_max,env_tke_max &
	18	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
	19	& ,alp_bl_conv,alp_bl_stat &
	20	!!! fin nrlmd le 10/04/2012
	21	& ,ztva )
	22
	23	USE dimphy
	24	USE ioipsl
	25	USE comgeomphy , ONLY:rlond,rlatd
	26	USE indice_sol_mod
	27	IMPLICIT NONE
	28
	29	!=======================================================================
	30	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
	31	! Version du 09.02.07
	32	! Calcul du transport vertical dans la couche limite en presence
	33	! de "thermiques" explicitement representes avec processus nuageux
	34	!
	35	! Reecriture a partir d'un listing papier a Habas, le 14/02/00
	36	!
	37	! le thermique est suppose homogene et dissipe par melange avec
	38	! son environnement. la longueur l_mix controle l'efficacite du
	39	! melange
	40	!
	41	! Le calcul du transport des differentes especes se fait en prenant
	42	! en compte:
	43	! 1. un flux de masse montant
	44	! 2. un flux de masse descendant
	45	! 3. un entrainement
	46	! 4. un detrainement
	47	!
	48	! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
	49	! Introduction of an implicit computation of vertical advection in
	50	! the environment of thermal plumes in thermcell_dq
	51	! impl = 0 : explicit, 1 : implicit, -1 : old version
	52	! controled by iflag_thermals =
	53	! 15, 16 run with impl=-1 : numerical convergence with NPv3
	54	! 17, 18 run with impl=1 : more stable
	55	! 15 and 17 correspond to the activation of the stratocumulus "bidouille"
	56	!
	57	!=======================================================================
	58
	59
	60	!-----------------------------------------------------------------------
	61	! declarations:
	62	! -------------
	63
	64	#include "dimensions.h"
	65	#include "YOMCST.h"
	66	#include "YOETHF.h"
	67	#include "FCTTRE.h"
	68	#include "iniprint.h"
	69	#include "thermcell.h"
	70
	71	! arguments:
	72	! ----------
	73
	74	!IM 140508
	75	INTEGER itap
	76
	77	INTEGER ngrid,nlay
	78	real ptimestep
	79	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
	80	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
	81	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
	82	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
	83	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	84	real pphi(ngrid,nlay)
	85
	86	! local:
	87	! ------
	88
	89	integer icount
	90
	91	integer, save :: dvdq=1,dqimpl=-1
	92	!$OMP THREADPRIVATE(dvdq,dqimpl)
	93	data icount/0/
	94	save icount
	95	!$OMP THREADPRIVATE(icount)
	96
	97	integer,save :: igout=1
	98	!$OMP THREADPRIVATE(igout)
	99	integer,save :: lunout1=6
	100	!$OMP THREADPRIVATE(lunout1)
	101	integer,save :: lev_out=10
	102	!$OMP THREADPRIVATE(lev_out)
	103
	104	REAL susqr2pi, Reuler
	105
	106	INTEGER ig,k,l,ll,ierr
[142]	107	real zsortie1d(ngrid)
	108	INTEGER lmax(ngrid),lmin(ngrid),lalim(ngrid)
	109	INTEGER lmix(ngrid)
	110	INTEGER lmix_bis(ngrid)
	111	real linter(ngrid)
	112	real zmix(ngrid)
	113	real zmax(ngrid),zw2(ngrid,nlay+1),ztva(ngrid,nlay),zw_est(ngrid,nlay+1),ztva_est(ngrid,nlay)
	114	! real fraca(ngrid,nlay)
[1]	115
[142]	116	real zmax_sec(ngrid)
[1]	117	!on garde le zmax du pas de temps precedent
[142]	118	real zmax0(ngrid)
[1]	119	!FH/IM save zmax0
	120
	121	real lambda
	122
[142]	123	real zlev(ngrid,nlay+1),zlay(ngrid,nlay)
	124	real deltaz(ngrid,nlay)
	125	REAL zh(ngrid,nlay)
	126	real zthl(ngrid,nlay),zdthladj(ngrid,nlay)
	127	REAL ztv(ngrid,nlay)
	128	real zu(ngrid,nlay),zv(ngrid,nlay),zo(ngrid,nlay)
	129	real zl(ngrid,nlay)
	130	real zsortie(ngrid,nlay)
	131	real zva(ngrid,nlay)
	132	real zua(ngrid,nlay)
	133	real zoa(ngrid,nlay)
[1]	134
[142]	135	real zta(ngrid,nlay)
	136	real zha(ngrid,nlay)
	137	real fraca(ngrid,nlay+1)
[1]	138	real zf,zf2
[142]	139	real thetath2(ngrid,nlay),wth2(ngrid,nlay),wth3(ngrid,nlay)
	140	real q2(ngrid,nlay)
[1]	141	! FH probleme de dimensionnement avec l'allocation dynamique
	142	! common/comtherm/thetath2,wth2
[142]	143	real wq(ngrid,nlay)
	144	real wthl(ngrid,nlay)
	145	real wthv(ngrid,nlay)
[1]	146
[142]	147	real ratqscth(ngrid,nlay)
[1]	148	real var
	149	real vardiff
[142]	150	real ratqsdiff(ngrid,nlay)
[1]	151
	152	logical sorties
[142]	153	real rho(ngrid,nlay),rhobarz(ngrid,nlay),masse(ngrid,nlay)
	154	real zpspsk(ngrid,nlay)
[1]	155
[142]	156	real wmax(ngrid)
	157	real wmax_tmp(ngrid)
	158	real wmax_sec(ngrid)
	159	real fm0(ngrid,nlay+1),entr0(ngrid,nlay),detr0(ngrid,nlay)
	160	real fm(ngrid,nlay+1),entr(ngrid,nlay),detr(ngrid,nlay)
[1]	161
[142]	162	real ztla(ngrid,nlay),zqla(ngrid,nlay),zqta(ngrid,nlay)
[1]	163	!niveau de condensation
[142]	164	integer nivcon(ngrid)
	165	real zcon(ngrid)
[1]	166	REAL CHI
[142]	167	real zcon2(ngrid)
	168	real pcon(ngrid)
	169	real zqsat(ngrid,nlay)
	170	real zqsatth(ngrid,nlay)
[1]	171
[142]	172	real f_star(ngrid,nlay+1),entr_star(ngrid,nlay)
	173	real detr_star(ngrid,nlay)
	174	real alim_star_tot(ngrid)
	175	real alim_star(ngrid,nlay)
	176	real alim_star_clos(ngrid,nlay)
	177	real f(ngrid), f0(ngrid)
[1]	178	!FH/IM save f0
[142]	179	real zlevinter(ngrid)
[1]	180	logical debut
	181	real seuil
[142]	182	real csc(ngrid,nlay)
[1]	183
	184	!!! nrlmd le 10/04/2012
	185
	186	!------Entrées
[142]	187	real pbl_tke(ngrid,nlay+1,nbsrf)
	188	real pctsrf(ngrid,nbsrf)
	189	real omega(ngrid,nlay)
	190	real airephy(ngrid)
[1]	191	!------Sorties
[142]	192	real zlcl(ngrid),fraca0(ngrid),w0(ngrid),w_conv(ngrid)
	193	real therm_tke_max0(ngrid),env_tke_max0(ngrid)
	194	real n2(ngrid),s2(ngrid)
	195	real ale_bl_stat(ngrid)
	196	real therm_tke_max(ngrid,nlay),env_tke_max(ngrid,nlay)
	197	real alp_bl_det(ngrid),alp_bl_fluct_m(ngrid),alp_bl_fluct_tke(ngrid),alp_bl_conv(ngrid),alp_bl_stat(ngrid)
[1]	198	!------Local
	199	integer nsrf
[142]	200	real rhobarz0(ngrid) ! Densité au LCL
	201	logical ok_lcl(ngrid) ! Existence du LCL des thermiques
	202	integer klcl(ngrid) ! Niveau du LCL
	203	real interp(ngrid) ! Coef d'interpolation pour le LCL
[1]	204	!--Triggering
	205	real Su ! Surface unité: celle d'un updraft élémentaire
	206	parameter(Su=4e4)
	207	real hcoef ! Coefficient directeur pour le calcul de s2
	208	parameter(hcoef=1)
	209	real hmincoef ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
	210	parameter(hmincoef=0.3)
	211	real eps1 ! Fraction de surface occupée par la population 1 : eps1=n1s1/(fraca0Sd)
	212	parameter(eps1=0.3)
	213	real hmin(ngrid) ! Ordonnée à l'origine pour le calcul de s2
	214	real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
	215	real zmax_moy_coef
	216	parameter(zmax_moy_coef=0.33)
[142]	217	real depth(ngrid) ! Epaisseur moyenne du cumulus
	218	real w_max(ngrid) ! Vitesse max statistique
	219	real s_max(ngrid)
[1]	220	!--Closure
[142]	221	real pbl_tke_max(ngrid,nlay) ! Profil de TKE moyenne
	222	real pbl_tke_max0(ngrid) ! TKE moyenne au LCL
	223	real w_ls(ngrid,nlay) ! Vitesse verticale grande échelle (m/s)
[1]	224	real coef_m ! On considère un rendement pour alp_bl_fluct_m
	225	parameter(coef_m=1.)
	226	real coef_tke ! On considère un rendement pour alp_bl_fluct_tke
	227	parameter(coef_tke=1.)
	228
	229	!!! fin nrlmd le 10/04/2012
	230
	231	!
	232	!nouvelles variables pour la convection
[142]	233	real Ale_bl(ngrid)
	234	real Alp_bl(ngrid)
	235	real alp_int(ngrid),dp_int(ngrid),zdp
	236	real ale_int(ngrid)
	237	integer n_int(ngrid)
	238	real fm_tot(ngrid)
	239	real wght_th(ngrid,nlay)
	240	integer lalim_conv(ngrid)
[1]	241	!v1d logical therm
	242	!v1d save therm
	243
	244	character*2 str2
	245	character*10 str10
	246
	247	character (len=20) :: modname='thermcell_main'
	248	character (len=80) :: abort_message
	249
	250	EXTERNAL SCOPY
	251	!
	252
[179]	253	! L. Fita, LMD July 2014. Initializing variables. (Should be done?)
[142]	254	! Some not initializated according to values: iflag_trig_bl, iflag_clos_bl
[178]	255	!! zdthladj = 0.
	256	!! pbl_tke_max0 = 0.
	257	!! fraca0 = 0.
	258	!! w_conv = 0.
	259	!! w0 = 0.
	260	!! therm_tke_max0 = 0.
	261	!! env_tke_max0 = 0.
	262	!! alp_bl_det = 0.
	263	!! alp_bl_fluct_m = 0.
	264	!! alp_bl_fluct_tke = 0.
	265	!! alp_bl_conv = 0.
	266	!! alp_bl_stat = 0.
	267	!! interp = 0.
	268	!! klcl = 0
	269
[1]	270	!-----------------------------------------------------------------------
	271	! initialisation:
	272	! ---------------
	273	!
	274
	275	seuil=0.25
	276
	277	if (debut) then
	278	! call getin('dvdq',dvdq)
	279	! call getin('dqimpl',dqimpl)
	280
	281	if (iflag_thermals==15.or.iflag_thermals==16) then
	282	dvdq=0
	283	dqimpl=-1
	284	else
	285	dvdq=1
	286	dqimpl=1
	287	endif
	288
	289	fm0=0.
	290	entr0=0.
	291	detr0=0.
	292
	293
	294	#undef wrgrads_thermcell
	295	#ifdef wrgrads_thermcell
	296	! Initialisation des sorties grads pour les thermiques.
	297	! Pour l'instant en 1D sur le point igout.
	298	! Utilise par thermcell_out3d.h
	299	str10='therm'
	300	call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
	301	& rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1., &
	302	& ptimestep,str10,'therm ')
	303	#endif
	304
	305
	306
	307	endif
	308
	309	fm=0. ; entr=0. ; detr=0.
	310
	311
	312	icount=icount+1
	313
	314	!IM 090508 beg
	315	!print*,'====================================================================='
	316	!print*,'====================================================================='
	317	!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
	318	!print*,'====================================================================='
	319	!print*,'====================================================================='
	320	!IM 090508 end
	321
	322	if (prt_level.ge.1) print*,'thermcell_main V4'
	323
	324	sorties=.true.
[142]	325	IF(ngrid.NE.ngrid) THEN
[1]	326	PRINT*
	327	PRINT*,'STOP dans convadj'
	328	PRINT*,'ngrid =',ngrid
[142]	329	PRINT*,'ngrid =',ngrid
[1]	330	ENDIF
	331	!
	332	! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
[142]	333	do ig=1,ngrid
[1]	334	f0(ig)=max(f0(ig),1.e-2)
	335	zmax0(ig)=max(zmax0(ig),40.)
	336	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
	337	enddo
	338
	339	if (prt_level.ge.20) then
	340	do ig=1,ngrid
	341	print*,'th_main ig f0',ig,f0(ig)
	342	enddo
	343	endif
	344	!-----------------------------------------------------------------------
	345	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
	346	! --------------------------------------------------------------------
	347	!
	348	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
	349	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
[142]	350
[1]	351	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
	352
	353	!------------------------------------------------------------------------
	354	! --------------------
	355	!
	356	!
	357	! + + + + + + + + + + +
	358	!
	359	!
	360	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
	361	! wh,wt,wo ...
	362	!
	363	! + + + + + + + + + + + zh,zu,zv,zo,rho
	364	!
	365	!
	366	! -------------------- zlev(1)
	367	! \\\\\\\\\\\\\\\\\\\\
	368	!
	369	!
	370
	371	!-----------------------------------------------------------------------
	372	! Calcul des altitudes des couches
	373	!-----------------------------------------------------------------------
	374
	375	do l=2,nlay
	376	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
	377	enddo
	378	zlev(:,1)=0.
[142]	379	zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG
[1]	380	do l=1,nlay
	381	zlay(:,l)=pphi(:,l)/RG
	382	enddo
	383	!calcul de l epaisseur des couches
	384	do l=1,nlay
	385	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
	386	enddo
	387
	388	! print*,'2 OK convect8'
	389	!-----------------------------------------------------------------------
	390	! Calcul des densites
	391	!-----------------------------------------------------------------------
	392
	393	rho(:,:)=pplay(:,:)/(zpspsk(:,:)RDztv(:,:))
	394
	395	if (prt_level.ge.10)write(lunout,*) &
	396	& 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
	397	rhobarz(:,1)=rho(:,1)
	398
	399	do l=2,nlay
	400	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
	401	enddo
	402
	403	!calcul de la masse
	404	do l=1,nlay
	405	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
	406	enddo
	407
	408	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
	409
	410	!------------------------------------------------------------------
	411	!
	412	! /\|\
	413	! -------- \| F_k+1 -------
	414	! ----> D_k
	415	! /\|\ <---- E_k , A_k
	416	! -------- \| F_k ---------
	417	! ----> D_k-1
	418	! <---- E_k-1 , A_k-1
	419	!
	420	!
	421	!
	422	!
	423	!
	424	! ---------------------------
	425	!
	426	! ----- F_lmax+1=0 ---------- \
	427	! lmax (zmax) \|
	428	! --------------------------- \|
	429	! \|
	430	! --------------------------- \|
	431	! \|
	432	! --------------------------- \|
	433	! \|
	434	! --------------------------- \|
	435	! \|
	436	! --------------------------- \|
	437	! \| E
	438	! --------------------------- \| D
	439	! \|
	440	! --------------------------- \|
	441	! \|
	442	! --------------------------- \ \|
	443	! lalim \| \|
	444	! --------------------------- \| \|
	445	! \| \|
	446	! --------------------------- \| \|
	447	! \| A \|
	448	! --------------------------- \| \|
	449	! \| \|
	450	! --------------------------- \| \|
	451	! lmin (=1 pour le moment) \| \|
	452	! ----- F_lmin=0 ------------ / /
	453	!
	454	! ---------------------------
	455	! //////////////////////////
	456	!
	457	!
	458	!=============================================================================
	459	! Calculs initiaux ne faisant pas intervenir les changements de phase
	460	!=============================================================================
	461
	462	!------------------------------------------------------------------
	463	! 1. alim_star est le profil vertical de l'alimentation a la base du
	464	! panache thermique, calcule a partir de la flotabilite de l'air sec
	465	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
	466	!------------------------------------------------------------------
	467	!
	468	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
	469	lmin=1
	470
	471	!-----------------------------------------------------------------------------
	472	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
	473	! panache sec conservatif (e=d=0) alimente selon alim_star
	474	! Il s'agit d'un calcul de type CAPE
	475	! zmax_sec est utilise pour determiner la geometrie du thermique.
	476	!------------------------------------------------------------------------------
	477	!---------------------------------------------------------------------------------
	478	!calcul du melange et des variables dans le thermique
	479	!--------------------------------------------------------------------------------
	480	!
	481	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
	482	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
	483
	484	! Gestion temporaire de plusieurs appels à thermcell_plume au travers
	485	! de la variable iflag_thermals
	486
	487	! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
	488	if (iflag_thermals_ed<=9) then
	489	! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
	490	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
	491	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	492	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	493	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
	494	& ,lev_out,lunout1,igout)
	495
	496	elseif (iflag_thermals_ed>9) then
	497	! print*,'THERM RIO et al 2010, version d Arnaud'
	498	CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
	499	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	500	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	501	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
	502	& ,lev_out,lunout1,igout)
	503
	504	endif
	505
	506	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
	507
	508	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
	509	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
	510
	511	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
	512	if (prt_level.ge.10) then
	513	write(lunout1,*) 'Dans thermcell_main 2'
	514	write(lunout1,*) 'lmin ',lmin(igout)
	515	write(lunout1,*) 'lalim ',lalim(igout)
	516	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
	517	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
	518	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
	519	endif
	520
	521	!-------------------------------------------------------------------------------
	522	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
	523	!-------------------------------------------------------------------------------
	524	!
	525	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
	526	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
	527	! Attention, w2 est transforme en sa racine carree dans cette routine
	528	! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
	529	wmax_tmp=0.
	530	do l=1,nlay
	531	wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
	532	enddo
	533	! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
	534
	535	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
	536	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
	537	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
	538	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
	539
	540	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
	541
	542	!-------------------------------------------------------------------------------
	543	! Fermeture,determination de f
	544	!-------------------------------------------------------------------------------
	545	!
	546	!
	547	!! write(lunout,*)'THERM NOUVEAU XXXXX'
	548	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
	549	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
	550
	551	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
	552	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
	553
	554	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
	555	if (prt_level.ge.10) then
	556	write(lunout1,*) 'Dans thermcell_main 1b'
	557	write(lunout1,*) 'lmin ',lmin(igout)
	558	write(lunout1,*) 'lalim ',lalim(igout)
	559	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
	560	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
	561	& ,l=1,lalim(igout)+4)
	562	endif
	563
	564	! Choix de la fonction d'alimentation utilisee pour la fermeture.
	565	! Apparemment sans importance
	566	alim_star_clos(:,:)=alim_star(:,:)
	567	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
	568
	569	! Appel avec la version seche
	570	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
	571	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
	572
	573	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	574	! Appel avec les zmax et wmax tenant compte de la condensation
	575	! Semble moins bien marcher
	576	! CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
	577	! & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
	578	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	579
	580	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
	581
	582	if (tau_thermals>1.) then
	583	lambda=exp(-ptimestep/tau_thermals)
	584	f0=(1.-lambda)f+lambdaf0
	585	else
	586	f0=f
	587	endif
	588
	589	! Test valable seulement en 1D mais pas genant
	590	if (.not. (f0(1).ge.0.) ) then
	591	abort_message = '.not. (f0(1).ge.0.)'
	592	CALL abort_gcm (modname,abort_message,1)
	593	endif
	594
	595	!-------------------------------------------------------------------------------
	596	!deduction des flux
	597	!-------------------------------------------------------------------------------
	598	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
	599	& lalim,lmax,alim_star, &
	600	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
	601	& detr,zqla,lev_out,lunout1,igout)
	602	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
	603
	604	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
	605	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
	606	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
	607
	608	!------------------------------------------------------------------
	609	! On ne prend pas directement les profils issus des calculs precedents
	610	! mais on s'autorise genereusement une relaxation vers ceci avec
	611	! une constante de temps tau_thermals (typiquement 1800s).
	612	!------------------------------------------------------------------
	613
	614	if (tau_thermals>1.) then
	615	lambda=exp(-ptimestep/tau_thermals)
	616	fm0=(1.-lambda)fm+lambdafm0
	617	entr0=(1.-lambda)entr+lambdaentr0
	618	detr0=(1.-lambda)detr+lambdadetr0
	619	else
	620	fm0=fm
	621	entr0=entr
	622	detr0=detr
	623	endif
	624
	625	!c------------------------------------------------------------------
	626	! calcul du transport vertical
	627	!------------------------------------------------------------------
	628
	629	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
	630	& zthl,zdthladj,zta,lev_out)
	631	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
	632	& po,pdoadj,zoa,lev_out)
	633	!------------------------------------------------------------------
	634	! Calcul de la fraction de l'ascendance
	635	!------------------------------------------------------------------
[142]	636	do ig=1,ngrid
[1]	637	fraca(ig,1)=0.
	638	fraca(ig,nlay+1)=0.
	639	enddo
	640	do l=2,nlay
[142]	641	do ig=1,ngrid
[1]	642	if (zw2(ig,l).gt.1.e-10) then
	643	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
	644	else
	645	fraca(ig,l)=0.
	646	endif
	647	enddo
	648	enddo
	649
	650	!------------------------------------------------------------------
	651	! calcul du transport vertical du moment horizontal
	652	!------------------------------------------------------------------
	653
	654	!IM 090508
	655	if (dvdq == 0 ) then
	656
	657	! Calcul du transport de V tenant compte d'echange par gradient
	658	! de pression horizontal avec l'environnement
	659
	660	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
	661	! & ,fraca*dvdq,zmax &
	662	& ,fraca,zmax &
	663	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
	664
	665	else
	666
	667	! calcul purement conservatif pour le transport de V
	668	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
	669	& ,zu,pduadj,zua,lev_out)
	670	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
	671	& ,zv,pdvadj,zva,lev_out)
	672
	673	endif
	674
	675	! print*,'13 OK convect8'
	676	do l=1,nlay
	677	do ig=1,ngrid
	678	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
	679	enddo
	680	enddo
	681
	682	if (prt_level.ge.1) print*,'14 OK convect8'
	683	!------------------------------------------------------------------
	684	! Calculs de diagnostiques pour les sorties
	685	!------------------------------------------------------------------
	686	!calcul de fraca pour les sorties
	687
	688	if (sorties) then
	689	if (prt_level.ge.1) print*,'14a OK convect8'
	690	! calcul du niveau de condensation
	691	! initialisation
	692	do ig=1,ngrid
	693	nivcon(ig)=0
	694	zcon(ig)=0.
	695	enddo
	696	!nouveau calcul
	697	do ig=1,ngrid
	698	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
	699	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
	700	enddo
	701	!IM do k=1,nlay
	702	do k=1,nlay-1
	703	do ig=1,ngrid
	704	if ((pcon(ig).le.pplay(ig,k)) &
	705	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
	706	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
	707	endif
	708	enddo
	709	enddo
	710	!IM
	711	ierr=0
	712	do ig=1,ngrid
	713	if (pcon(ig).le.pplay(ig,nlay)) then
	714	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
	715	ierr=1
	716	endif
	717	enddo
	718	if (ierr==1) then
	719	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
	720	CALL abort_gcm (modname,abort_message,1)
	721	endif
	722
	723	if (prt_level.ge.1) print*,'14b OK convect8'
	724	do k=nlay,1,-1
	725	do ig=1,ngrid
	726	if (zqla(ig,k).gt.1e-10) then
	727	nivcon(ig)=k
	728	zcon(ig)=zlev(ig,k)
	729	endif
	730	enddo
	731	enddo
	732	if (prt_level.ge.1) print*,'14c OK convect8'
	733	!calcul des moments
	734	!initialisation
	735	do l=1,nlay
	736	do ig=1,ngrid
	737	q2(ig,l)=0.
	738	wth2(ig,l)=0.
	739	wth3(ig,l)=0.
	740	ratqscth(ig,l)=0.
	741	ratqsdiff(ig,l)=0.
	742	enddo
	743	enddo
	744	if (prt_level.ge.1) print*,'14d OK convect8'
	745	if (prt_level.ge.10)write(lunout,*) &
	746	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
	747	do l=1,nlay
	748	do ig=1,ngrid
	749	zf=fraca(ig,l)
	750	zf2=zf/(1.-zf)
	751	!
	752	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
	753	if(zw2(ig,l).gt.1.e-10) then
	754	wth2(ig,l)=zf2(zw2(ig,l))*2
	755	else
	756	wth2(ig,l)=0.
	757	endif
	758	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
	759	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
	760	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
	761	!test: on calcul q2/po=ratqsc
	762	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
	763	enddo
	764	enddo
[142]	765
[1]	766	!calcul des flux: q, thetal et thetav
	767	do l=1,nlay
	768	do ig=1,ngrid
	769	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
	770	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
	771	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
	772	enddo
	773	enddo
	774	!
	775
	776	!!! nrlmd le 10/04/2012
	777
	778	!------------Test sur le LCL des thermiques
	779	do ig=1,ngrid
	780	ok_lcl(ig)=.false.
[142]	781	if ( (pcon(ig) .gt. pplay(ig,nlay-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
[1]	782	enddo
	783
	784	!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
	785	do l=1,nlay-1
	786	do ig=1,ngrid
	787	if (ok_lcl(ig)) then
[142]	788	! if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
	789	! L. Fita, LMD July 2014. Adding avoiding divisons by zero...
	790	if ((pplay(ig,l).ge.pcon(ig)) .and. (pplay(ig,l+1).le.pcon(ig)) .and. &
	791	(klcl(ig).gt.0).and.(klcl(ig)+1.le.nlay-1) .and. (klcl(ig)+1.gt.0) ) then
	792	klcl(ig)=l
	793	interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
[1]	794	endif
	795	endif
	796	enddo
	797	enddo
	798
	799	!------------Hauteur des thermiques
	800	!!jyg le 27/04/2012
	801	!! do ig =1,ngrid
	802	!! rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
	803	!! & -rhobarz(ig,klcl(ig)))*interp(ig)
	804	!! zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
	805	!! zmax(ig)=pphi(ig,lmax(ig))/rg
	806	!! if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
	807	!! enddo
	808	do ig =1,ngrid
	809	zmax(ig)=pphi(ig,lmax(ig))/rg
[142]	810	! if (ok_lcl(ig)) then
	811	! L. Fita, LMD July 2014. Adding avoiding divisons by zero...
	812	if ( ok_lcl(ig) .and. (klcl(ig).gt.0) .and. (klcl(ig)+1.le.nlay-1) .and. &
	813	(klcl(ig)+1.gt.0) ) then
	814	rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
[1]	815	& -rhobarz(ig,klcl(ig)))*interp(ig)
[142]	816	zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
	817	zlcl(ig)=min(zlcl(ig),zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax
[1]	818	else
	819	rhobarz0(ig)=0.
	820	zlcl(ig)=zmax(ig)
	821	endif
	822	enddo
	823	!!jyg fin
	824
	825	!------------Calcul des propriétés du thermique au LCL
	826	IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
	827
	828	!-----Initialisation de la TKE moyenne
	829	do l=1,nlay
	830	do ig=1,ngrid
	831	pbl_tke_max(ig,l)=0.
	832	enddo
	833	enddo
	834
	835	!-----Calcul de la TKE moyenne
	836	do nsrf=1,nbsrf
	837	do l=1,nlay
	838	do ig=1,ngrid
	839	pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
	840	enddo
	841	enddo
	842	enddo
	843
	844	!-----Initialisations des TKE dans et hors des thermiques
	845	do l=1,nlay
	846	do ig=1,ngrid
	847	therm_tke_max(ig,l)=pbl_tke_max(ig,l)
	848	env_tke_max(ig,l)=pbl_tke_max(ig,l)
	849	enddo
	850	enddo
	851
	852	!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
	853	call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
	854	& rg,pplev,therm_tke_max)
	855	! print *,' thermcell_tke_transport -> ' !!jyg
	856
	857	!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
	858	do l=1,nlay
	859	do ig=1,ngrid
	860	pbl_tke_max(ig,l)=fraca(ig,l)therm_tke_max(ig,l)+(1.-fraca(ig,l))env_tke_max(ig,l) ! Recalcul de TKE moyenne aprés transport de TKE_TH
	861	env_tke_max(ig,l)=(pbl_tke_max(ig,l)-fraca(ig,l)*therm_tke_max(ig,l))/(1.-fraca(ig,l)) ! Recalcul de TKE dans l'environnement aprés transport de TKE_TH
	862	w_ls(ig,l)=-1.omega(ig,l)/(RGrhobarz(ig,l)) ! Vitesse verticale de grande échelle
	863	enddo
	864	enddo
	865	! print *,' apres w_ls = ' !!jyg
	866
	867	do ig=1,ngrid
	868	if (ok_lcl(ig)) then
	869	fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
	870	& -fraca(ig,klcl(ig)))*interp(ig)
	871	w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
	872	& -zw2(ig,klcl(ig)))*interp(ig)
	873	w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
	874	& -w_ls(ig,klcl(ig)))*interp(ig)
	875	therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
	876	& +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
	877	env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
	878	& -env_tke_max(ig,klcl(ig)))*interp(ig)
	879	pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
	880	& -pbl_tke_max(ig,klcl(ig)))*interp(ig)
	881	if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
	882	if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
	883	if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
	884	else
	885	fraca0(ig)=0.
	886	w0(ig)=0.
[142]	887	! L. Fita, LMD July 2014. Adding zero value for stability issues
	888	w_conv(ig) = 0.
[1]	889	!!jyg le 27/04/2012
	890	!! zlcl(ig)=0.
	891	!!
	892	endif
	893	enddo
	894
	895	ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
	896	! print *,'ENDIF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ' !!jyg
	897
	898	!------------Triggering------------------
	899	IF (iflag_trig_bl.ge.1) THEN
	900
	901	!-----Initialisations
	902	depth(:)=0.
	903	n2(:)=0.
	904	s2(:)=0.
	905	s_max(:)=0.
	906
	907	!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
	908	do ig=1,ngrid
	909	zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
	910	depth(ig)=zmax_moy(ig)-zlcl(ig)
	911	hmin(ig)=hmincoef*zlcl(ig)
	912	if (depth(ig).ge.10.) then
	913	s2(ig)=(hcoefdepth(ig)+hmin(ig))*2
	914	n2(ig)=(1.-eps1)fraca0(ig)airephy(ig)/s2(ig)
	915	!!
	916	!!jyg le 27/04/2012
	917	!! s_max(ig)=s2(ig)*log(n2(ig))
	918	!! if (n2(ig) .lt. 1) s_max(ig)=0.
	919	s_max(ig)=s2(ig)*log(max(n2(ig),1.))
	920	!!fin jyg
	921	else
	922	s2(ig)=0.
	923	n2(ig)=0.
	924	s_max(ig)=0.
	925	endif
	926	enddo
	927	! print *,'avant Calcul de Wmax ' !!jyg
	928
	929	!-----Calcul de Wmax et ALE_BL_STAT associée
	930	!!jyg le 30/04/2012
	931	!! do ig=1,ngrid
	932	!! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
	933	!! w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/su)-log(2.3.14)-log(2.log(s_max(ig)/su)-log(2.*3.14))))
	934	!! ale_bl_stat(ig)=0.5w_max(ig)*2
	935	!! else
	936	!! w_max(ig)=0.
	937	!! ale_bl_stat(ig)=0.
	938	!! endif
	939	!! enddo
	940	susqr2pi=susqrt(2.Rpi)
	941	Reuler=exp(1.)
	942	do ig=1,ngrid
	943	if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
	944	w_max(ig)=w0(ig)(1.+sqrt(2.log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
	945	ale_bl_stat(ig)=0.5w_max(ig)*2
	946	else
	947	w_max(ig)=0.
	948	ale_bl_stat(ig)=0.
	949	endif
	950	enddo
	951
	952	ENDIF ! iflag_trig_bl
	953	! print *,'ENDIF iflag_trig_bl' !!jyg
	954	!------------Closure------------------
	955
	956	IF (iflag_clos_bl.ge.1) THEN
	957
	958	!-----Calcul de ALP_BL_STAT
	959	do ig=1,ngrid
	960	alp_bl_det(ig)=0.5coef_mrhobarz0(ig)(w0(ig)3)fraca0(ig)(1.-2.fraca0(ig))/((1.-fraca0(ig))**2)
	961	alp_bl_fluct_m(ig)=1.5rhobarz0(ig)fraca0(ig)(w_conv(ig)+coef_mw0(ig))* &
	962	& (w0(ig)**2)
	963	alp_bl_fluct_tke(ig)=3.coef_mrhobarz0(ig)w0(ig)fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
	964	& +3.rhobarz0(ig)w_conv(ig)*pbl_tke_max0(ig)
	965	if (iflag_clos_bl.ge.2) then
	966	alp_bl_conv(ig)=1.5coef_mrhobarz0(ig)fraca0(ig)(fraca0(ig)/(1.-fraca0(ig)))w_conv(ig) &
	967	& (w0(ig)**2)
	968	else
	969	alp_bl_conv(ig)=0.
	970	endif
	971	alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
	972	enddo
	973
	974	!-----Sécurité ALP infinie
	975	do ig=1,ngrid
	976	if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
	977	enddo
	978
	979	ENDIF ! (iflag_clos_bl.ge.1)
	980
	981	!!! fin nrlmd le 10/04/2012
	982
	983	if (prt_level.ge.10) then
	984	ig=igout
	985	do l=1,nlay
	986	print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
	987	print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
	988	enddo
	989	endif
	990
	991	! print*,'avant calcul ale et alp'
	992	!calcul de ALE et ALP pour la convection
	993	Alp_bl(:)=0.
	994	Ale_bl(:)=0.
	995	! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
	996	do l=1,nlay
	997	do ig=1,ngrid
	998	Alp_bl(ig)=max(Alp_bl(ig),0.5rhobarz(ig,l)wth3(ig,l) )
	999	Ale_bl(ig)=max(Ale_bl(ig),0.5zw2(ig,l)*2)
	1000	! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
	1001	enddo
	1002	enddo
	1003
	1004	!test:calcul de la ponderation des couches pour KE
	1005	!initialisations
	1006
	1007	fm_tot(:)=0.
	1008	wght_th(:,:)=1.
	1009	lalim_conv(:)=lalim(:)
	1010
[142]	1011	do k=1,nlay
[1]	1012	do ig=1,ngrid
	1013	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
	1014	enddo
	1015	enddo
	1016
	1017	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
	1018	! plus petit que 1.e-10, on prend wght_th=1.
[142]	1019	do k=1,nlay
[1]	1020	do ig=1,ngrid
	1021	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
	1022	wght_th(ig,k)=alim_star(ig,k)
	1023	endif
	1024	enddo
	1025	enddo
	1026
	1027	! print*,'apres wght_th'
	1028	!test pour prolonger la convection
	1029	do ig=1,ngrid
	1030	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
	1031	if ((alim_star(ig,1).lt.1.e-10)) then
	1032	lalim_conv(ig)=1
	1033	wght_th(ig,1)=1.
	1034	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
	1035	endif
	1036	enddo
	1037
	1038	!------------------------------------------------------------------------
	1039	! Modif CR/FH 20110310 : Alp integree sur la verticale.
	1040	! Integrale verticale de ALP.
	1041	! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
	1042	! couches
	1043	!------------------------------------------------------------------------
	1044
	1045	alp_int(:)=0.
	1046	dp_int(:)=0.
	1047	do l=2,nlay
	1048	do ig=1,ngrid
	1049	if(l.LE.lmax(ig)) THEN
	1050	zdp=pplay(ig,l-1)-pplay(ig,l)
	1051	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)*zdp
	1052	dp_int(ig)=dp_int(ig)+zdp
	1053	endif
	1054	enddo
	1055	enddo
	1056
	1057	if (iflag_coupl>=3 .and. iflag_coupl<=5) then
	1058	do ig=1,ngrid
	1059	!valeur integree de alp_bl * 0.5:
	1060	if (dp_int(ig)>0.) then
	1061	Alp_bl(ig)=alp_int(ig)/dp_int(ig)
	1062	endif
	1063	enddo!
	1064	endif
	1065
	1066
	1067	! Facteur multiplicatif sur Alp_bl
	1068	Alp_bl(:)=alp_bl_k*Alp_bl(:)
	1069
	1070	!------------------------------------------------------------------------
	1071
	1072
	1073	!calcul du ratqscdiff
	1074	if (prt_level.ge.1) print*,'14e OK convect8'
	1075	var=0.
	1076	vardiff=0.
	1077	ratqsdiff(:,:)=0.
	1078
[142]	1079	do l=1,nlay
[1]	1080	do ig=1,ngrid
	1081	if (l<=lalim(ig)) then
	1082	var=var+alim_star(ig,l)zqta(ig,l)1000.
	1083	endif
	1084	enddo
	1085	enddo
	1086
	1087	if (prt_level.ge.1) print*,'14f OK convect8'
	1088
[142]	1089	do l=1,nlay
[1]	1090	do ig=1,ngrid
	1091	if (l<=lalim(ig)) then
	1092	zf=fraca(ig,l)
	1093	zf2=zf/(1.-zf)
	1094	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
	1095	endif
	1096	enddo
	1097	enddo
	1098
	1099	if (prt_level.ge.1) print*,'14g OK convect8'
	1100	do l=1,nlay
	1101	do ig=1,ngrid
	1102	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
	1103	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
	1104	enddo
	1105	enddo
	1106	!--------------------------------------------------------------------
	1107	!
	1108	!ecriture des fichiers sortie
	1109	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
	1110
	1111	#ifdef wrgrads_thermcell
	1112	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
	1113	#include "thermcell_out3d.h"
	1114	#endif
	1115
	1116	endif
	1117
	1118	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
	1119
	1120	return
	1121	end
	1122
	1123	!-----------------------------------------------------------------------------
	1124
[142]	1125	subroutine test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
[1]	1126	IMPLICIT NONE
	1127	#include "iniprint.h"
	1128
[142]	1129	integer i, k, ngrid,nlay
	1130	real pplev(ngrid,nlay+1),pplay(ngrid,nlay)
	1131	real ztv(ngrid,nlay)
	1132	real po(ngrid,nlay)
	1133	real ztva(ngrid,nlay)
	1134	real zqla(ngrid,nlay)
	1135	real f_star(ngrid,nlay)
	1136	real zw2(ngrid,nlay)
	1137	integer long(ngrid)
[1]	1138	real seuil
	1139	character*21 comment
	1140
	1141	if (prt_level.ge.1) THEN
	1142	print*,'WARNING !!! TEST ',comment
	1143	endif
	1144	return
	1145
	1146	! test sur la hauteur des thermiques ...
[142]	1147	do i=1,ngrid
[1]	1148	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
	1149	if (prt_level.ge.10) then
	1150	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
	1151	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
[142]	1152	do k=1,nlay
[1]	1153	write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000po(i,k),ztva(i,k),1000zqla(i,k),f_star(i,k),zw2(i,k)
	1154	enddo
	1155	endif
	1156	enddo
	1157
	1158
	1159	return
	1160	end
	1161
	1162	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
	1163	! On transporte pbl_tke pour donner therm_tke
	1164	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
	1165	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
	1166	& rg,pplev,therm_tke_max)
	1167	implicit none
	1168
	1169	#include "iniprint.h"
	1170	!=======================================================================
	1171	!
	1172	! Calcul du transport verticale dans la couche limite en presence
	1173	! de "thermiques" explicitement representes
	1174	! calcul du dq/dt une fois qu'on connait les ascendances
	1175	!
	1176	!=======================================================================
	1177
	1178	integer ngrid,nlay,nsrf
	1179
	1180	real ptimestep
	1181	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
	1182	real entr0(ngrid,nlay),rg
	1183	real therm_tke_max(ngrid,nlay)
	1184	real detr0(ngrid,nlay)
	1185
	1186
	1187	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	1188	real entr(ngrid,nlay)
	1189	real q(ngrid,nlay)
	1190	integer lev_out ! niveau pour les print
	1191
	1192	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
	1193
	1194	real zzm
	1195
	1196	integer ig,k
	1197	integer isrf
	1198
	1199
	1200	lev_out=0
	1201
	1202
	1203	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
	1204
	1205	! calcul du detrainement
	1206	do k=1,nlay
	1207	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
	1208	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
	1209	enddo
	1210
	1211
	1212	! Decalage vertical des entrainements et detrainements.
	1213	masse(:,1)=0.5*masse0(:,1)
	1214	entr(:,1)=0.5*entr0(:,1)
	1215	detr(:,1)=0.5*detr0(:,1)
	1216	fm(:,1)=0.
	1217	do k=1,nlay-1
	1218	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
	1219	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
	1220	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
	1221	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
	1222	enddo
	1223	fm(:,nlay+1)=0.
	1224
	1225	!!! nrlmd le 16/09/2010
	1226	! calcul de la valeur dans les ascendances
	1227	! do ig=1,ngrid
	1228	! qa(ig,1)=q(ig,1)
	1229	! enddo
	1230	!!!
	1231
	1232	!do isrf=1,nsrf
	1233
	1234	! q(:,:)=therm_tke(:,:,isrf)
	1235	q(:,:)=therm_tke_max(:,:)
	1236	!!! nrlmd le 16/09/2010
	1237	do ig=1,ngrid
	1238	qa(ig,1)=q(ig,1)
	1239	enddo
	1240	!!!
	1241
	1242	if (1==1) then
	1243	do k=2,nlay
	1244	do ig=1,ngrid
	1245	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
	1246	& 1.e-5*masse(ig,k)) then
	1247	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
	1248	& /(fm(ig,k+1)+detr(ig,k))
	1249	else
	1250	qa(ig,k)=q(ig,k)
	1251	endif
	1252	if (qa(ig,k).lt.0.) then
	1253	! print*,'qa<0!!!'
	1254	endif
	1255	if (q(ig,k).lt.0.) then
	1256	! print*,'q<0!!!'
	1257	endif
	1258	enddo
	1259	enddo
	1260
	1261	! Calcul du flux subsident
	1262
	1263	do k=2,nlay
	1264	do ig=1,ngrid
	1265	wqd(ig,k)=fm(ig,k)*q(ig,k)
	1266	if (wqd(ig,k).lt.0.) then
	1267	! print*,'wqd<0!!!'
	1268	endif
	1269	enddo
	1270	enddo
	1271	do ig=1,ngrid
	1272	wqd(ig,1)=0.
	1273	wqd(ig,nlay+1)=0.
	1274	enddo
	1275
	1276	! Calcul des tendances
	1277	do k=1,nlay
	1278	do ig=1,ngrid
	1279	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
	1280	& -wqd(ig,k)+wqd(ig,k+1)) &
	1281	& *ptimestep/masse(ig,k)
	1282	enddo
	1283	enddo
	1284
	1285	endif
	1286
	1287	therm_tke_max(:,:)=q(:,:)
	1288
	1289	return
	1290	!!! fin nrlmd le 10/04/2012
	1291	end
	1292

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

Download in other formats:

Original Format