Context Navigation

thermcell_main.F90 @ 2337

Last change on this file since 2337 was 2220, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2186:2216 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 40.8 KB

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

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

Download in other formats:

Original Format