Context Navigation

thermcell_main.F90 @ 1678

Last change on this file since 1678 was 1638, checked in by idelkadi, 12 years ago
Introduction du declenchement stochastique de la convection
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 40.1 KB

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

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

Download in other formats:

Original Format