Context Navigation

lmdz_thermcell_main.F90 @ 5342

Last change on this file since 5342 was 4843, checked in by crio, 10 months ago
Nouvelle formulation du strig et correction thermiques montent trop haut
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: 35.1 KB

Rev	Line
[4590]	1	MODULE lmdz_thermcell_main
[1403]	2	! $Id: lmdz_thermcell_main.F90 4843 2024-03-04 17:58:03Z abarral $
[878]	3	!
[4678]	4	! A REGARDER !!!!!!!!!!!!!!!!!
	5	! ATTENTION : zpspsk est inout et out mais c'est pas forcement pour de bonnes raisons (FH, 2023)
[4690]	6	! ATTENTION : dans thermcell_env, on condense potentiellement de l'eau. Mais comme on ne mélange pas l'eau liquide supposant qu'il n'y en n'a pas, c'est potentiellement un souci
[4590]	7	CONTAINS
	8
[4089]	9	subroutine thermcell_main(itap,ngrid,nlay,ptimestep &
[878]	10	& ,pplay,pplev,pphi,debut &
[4690]	11	& ,puwind,pvwind,ptemp,p_o,ptemp_env, po_env &
[878]	12	& ,pduadj,pdvadj,pdtadj,pdoadj &
[1026]	13	& ,fm0,entr0,detr0,zqta,zqla,lmax &
[878]	14	& ,ratqscth,ratqsdiff,zqsatth &
[1403]	15	& ,zmax0, f0,zw2,fraca,ztv &
[4089]	16	& ,zpspsk,ztla,zthl,ztva &
[4843]	17	& ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,zcong &
[4089]	18	#ifdef ISO
	19	& ,xtpo,xtpdoadj &
	20	#endif
	21	& )
[878]	22
[4089]	23
[4690]	24
	25	! USE necessaires pour les lignes importees de thermcell_env
[4590]	26	USE lmdz_thermcell_ini, ONLY: thermcell_ini,dqimpl,dvdq,prt_level,lunout,prt_level
	27	USE lmdz_thermcell_ini, ONLY: iflag_thermals_closure,iflag_thermals_ed,tau_thermals,r_aspect_thermals
	28	USE lmdz_thermcell_ini, ONLY: iflag_thermals_down,fact_thermals_down
[4690]	29	USE lmdz_thermcell_ini, ONLY: iflag_thermals_tenv
[4590]	30	USE lmdz_thermcell_ini, ONLY: RD,RG
[4089]	31
[4590]	32	USE lmdz_thermcell_down, ONLY: thermcell_updown_dq
	33	USE lmdz_thermcell_closure, ONLY: thermcell_closure
	34	USE lmdz_thermcell_dq, ONLY: thermcell_dq
	35	USE lmdz_thermcell_dry, ONLY: thermcell_dry
	36	USE lmdz_thermcell_dv2, ONLY: thermcell_dv2
	37	USE lmdz_thermcell_env, ONLY: thermcell_env
	38	USE lmdz_thermcell_flux2, ONLY: thermcell_flux2
	39	USE lmdz_thermcell_height, ONLY: thermcell_height
	40	USE lmdz_thermcell_plume, ONLY: thermcell_plume
	41	USE lmdz_thermcell_plume_6A, ONLY: thermcell_plume_6A,thermcell_plume_5B
	42
[4690]	43	! USE necessaires pour les lignes importees de thermcell_env
	44	USE lmdz_thermcell_ini, ONLY : RLvCp,RKAPPA,RETV
	45	USE lmdz_thermcell_qsat, ONLY : thermcell_qsat
	46
	47
[4089]	48	#ifdef ISO
[4143]	49	USE infotrac_phy, ONLY : ntiso
[4089]	50	#ifdef ISOVERIF
	51	USE isotopes_mod, ONLY : iso_eau,iso_HDO
	52	USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
	53	iso_verif_aberrant_encadre
	54	#endif
	55	#endif
	56
	57
[878]	58	IMPLICIT NONE
	59
	60	!=======================================================================
	61	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
	62	! Version du 09.02.07
	63	! Calcul du transport vertical dans la couche limite en presence
	64	! de "thermiques" explicitement representes avec processus nuageux
	65	!
[1403]	66	! Reecriture a partir d'un listing papier a Habas, le 14/02/00
[878]	67	!
[1403]	68	! le thermique est suppose homogene et dissipe par melange avec
	69	! son environnement. la longueur l_mix controle l'efficacite du
	70	! melange
[878]	71	!
[1403]	72	! Le calcul du transport des differentes especes se fait en prenant
[878]	73	! en compte:
	74	! 1. un flux de masse montant
	75	! 2. un flux de masse descendant
	76	! 3. un entrainement
	77	! 4. un detrainement
	78	!
[1738]	79	! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
	80	! Introduction of an implicit computation of vertical advection in
	81	! the environment of thermal plumes in thermcell_dq
	82	! impl = 0 : explicit, 1 : implicit, -1 : old version
	83	! controled by iflag_thermals =
	84	! 15, 16 run with impl=-1 : numerical convergence with NPv3
	85	! 17, 18 run with impl=1 : more stable
	86	! 15 and 17 correspond to the activation of the stratocumulus "bidouille"
	87	!
[4133]	88	! Using
	89	! abort_physic
	90	! iso_verif_aberrant_encadre
	91	! iso_verif_egalite
	92	! test_ltherm
	93	! thermcell_closure
	94	! thermcell_dq
	95	! thermcell_dry
	96	! thermcell_dv2
	97	! thermcell_env
	98	! thermcell_flux2
	99	! thermcell_height
	100	! thermcell_plume
	101	! thermcell_plume_5B
	102	! thermcell_plume_6A
	103	!
[878]	104	!=======================================================================
	105
[1738]	106
[878]	107	!-----------------------------------------------------------------------
	108	! declarations:
	109	! -------------
	110
	111
	112	! arguments:
	113	! ----------
[4089]	114	integer, intent(in) :: itap,ngrid,nlay
	115	real, intent(in) :: ptimestep
[4690]	116	real, intent(in), dimension(ngrid,nlay) :: ptemp,puwind,pvwind,pplay,pphi,ptemp_env,po_env
[4678]	117	! ATTENTION : zpspsk est inout et out mais c'est pas forcement pour de bonnes raisons (FH, 2023)
	118	real, intent(in), dimension(ngrid,nlay) :: p_o
[4590]	119	real, intent(out), dimension(ngrid,nlay) :: zpspsk
[4089]	120	real, intent(in), dimension(ngrid,nlay+1) :: pplev
[4133]	121	integer, intent(out), dimension(ngrid) :: lmax
[4089]	122	real, intent(out), dimension(ngrid,nlay) :: pdtadj,pduadj,pdvadj,pdoadj,entr0,detr0
	123	real, intent(out), dimension(ngrid,nlay) :: ztla,zqla,zqta,zqsatth,zthl
	124	real, intent(out), dimension(ngrid,nlay+1) :: fm0,zw2,fraca
[4133]	125	real, intent(inout), dimension(ngrid) :: zmax0,f0
[4089]	126	real, intent(out), dimension(ngrid,nlay) :: ztva,ztv
	127	logical, intent(in) :: debut
[4133]	128	real,intent(out), dimension(ngrid,nlay) :: ratqscth,ratqsdiff
[878]	129
[4089]	130	real, intent(out), dimension(ngrid) :: pcon
	131	real, intent(out), dimension(ngrid,nlay) :: rhobarz,wth3
	132	real, intent(out), dimension(ngrid) :: wmax_sec
	133	integer,intent(out), dimension(ngrid) :: lalim
	134	real, intent(out), dimension(ngrid,nlay+1) :: fm
	135	real, intent(out), dimension(ngrid,nlay) :: alim_star
[4843]	136	real, intent(out), dimension(ngrid) :: zmax,zcong
[972]	137
[878]	138	! local:
	139	! ------
	140
[1738]	141
[883]	142	integer,save :: igout=1
[987]	143	!$OMP THREADPRIVATE(igout)
[938]	144	integer,save :: lunout1=6
[987]	145	!$OMP THREADPRIVATE(lunout1)
[883]	146	integer,save :: lev_out=10
[987]	147	!$OMP THREADPRIVATE(lev_out)
[878]	148
[4089]	149	real lambda, zf,zf2,var,vardiff,CHI
	150	integer ig,k,l,ierr,ll
[878]	151	logical sorties
[4843]	152	real, dimension(ngrid) :: linter,zmix, zmax_sec,lintercong
	153	integer,dimension(ngrid) :: lmin,lmix,lmix_bis,nivcon, lcong
[4089]	154	real, dimension(ngrid,nlay) :: ztva_est
[4690]	155	real, dimension(ngrid,nlay) :: deltaz,zlay,zdthladj,zu,zv,z_o,zl,zva,zua,z_oa
	156	real, dimension(ngrid,nlay) :: ztemp_env ! temperarure liquide de l'environnement
[4089]	157	real, dimension(ngrid,nlay) :: zta,zha,q2,wq,wthl,wthv,thetath2,wth2
[4133]	158	real, dimension(ngrid,nlay) :: rho,masse
[4089]	159	real, dimension(ngrid,nlay+1) :: zw_est,zlev
	160	real, dimension(ngrid) :: wmax,wmax_tmp
	161	real, dimension(ngrid,nlay+1) :: f_star
	162	real, dimension(ngrid,nlay) :: entr,detr,entr_star,detr_star,alim_star_clos
	163	real, dimension(ngrid,nlay) :: zqsat,csc
	164	real, dimension(ngrid) :: zcon,zcon2,alim_star_tot,f
[4413]	165	real, dimension(ngrid,nlay) :: entrdn,detrdn
[4690]	166	logical, dimension(ngrid,nlay) :: mask
[878]	167
[4089]	168	character (len=20) :: modname='thermcell_main'
	169	character (len=80) :: abort_message
[878]	170
	171
[4089]	172	#ifdef ISO
[4143]	173	REAL xtpo(ntiso,ngrid,nlay),xtpdoadj(ntiso,ngrid,nlay)
	174	REAL xtzo(ntiso,ngrid,nlay)
[4089]	175	REAL xtpdoadj_tmp(ngrid,nlay)
	176	REAL xtpo_tmp(ngrid,nlay)
	177	REAL xtzo_tmp(ngrid,nlay)
	178	integer ixt
	179	#endif
[878]	180
	181	!
	182
	183	!-----------------------------------------------------------------------
	184	! initialisation:
	185	! ---------------
	186	!
[1943]	187	fm=0. ; entr=0. ; detr=0.
[972]	188
[938]	189	if (prt_level.ge.1) print*,'thermcell_main V4'
[878]	190
	191	sorties=.true.
[4089]	192	IF(ngrid.NE.ngrid) THEN
[878]	193	PRINT*
	194	PRINT*,'STOP dans convadj'
	195	PRINT*,'ngrid =',ngrid
[4089]	196	PRINT*,'ngrid =',ngrid
[878]	197	ENDIF
	198	!
[4692]	199	!print*,'thermcell_main debut'
[1403]	200	! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
[4089]	201	do ig=1,ngrid
[972]	202	f0(ig)=max(f0(ig),1.e-2)
[1403]	203	zmax0(ig)=max(zmax0(ig),40.)
[972]	204	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
	205	enddo
[878]	206
[1494]	207	if (prt_level.ge.20) then
	208	do ig=1,ngrid
	209	print*,'th_main ig f0',ig,f0(ig)
	210	enddo
	211	endif
[4690]	212
[878]	213	!-----------------------------------------------------------------------
	214	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
	215	! --------------------------------------------------------------------
	216	!
[4690]	217	! On condense l'eau liquide si besoin.
	218	! En fait on arrive ici d'habitude (jusque 6A) après réévaporation
	219	! Dans une nouvelle mouture, on passe les profiles
	220	! avant la couche limite : iflag_thermals_tenv=1
	221	! dés le début de la physique : iflag_thermals_tenv=2
	222	! Mais même pour 2) on ne veut sans doute pas réévaporer.
	223	! On veut comparer thetav dans le thermique, après condensation,
	224	! avec le theta_v effectif de l'environnement.
	225
	226	if (iflag_thermals_tenv - 10 * ( iflag_thermals_tenv / 10 ) == 0) then
	227
	228	CALL thermcell_env(ngrid,nlay,p_o,ptemp_env,puwind,pvwind,pplay, &
[4843]	229	& pplev,z_o,ztemp_env,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lcong,lintercong,lev_out)
[4690]	230
	231	else
	232
	233	! Chantier en cours : ne pas effacer (Fredho). 15 septembre 2023
	234	! Dans la version originale de thermcell_env, on condense l'eau de l'environnement
	235	! pour calculer une temperature potentielle liquide.
	236	! On en déduit un Theta v.
	237
	238	! ...
	239	! contenu de thermcell_env
	240	! SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
	241	! & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out)
	242	! contenu thermcell_env : call thermcell_qsat(ngrid*nlay,mask,pplev,pt,po,pqsat)
	243	! contenu thermcell_env : do ll=1,nlay
	244	! contenu thermcell_env : do ig=1,ngrid
	245	! contenu thermcell_env : zl(ig,ll) = max(0.,po(ig,ll)-pqsat(ig,ll))
	246	! contenu thermcell_env : zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll) ! T = Tl + Lv/Cp ql
	247	! contenu thermcell_env : zo(ig,ll) = po(ig,ll)-zl(ig,ll)
	248	! contenu thermcell_env : enddo
	249	! contenu thermcell_env : enddo
	250	! contenu thermcell_env : do ll=1,nlay
	251	! contenu thermcell_env : do ig=1,ngrid
	252	! contenu thermcell_env : zpspsk(ig,ll)=(pplay(ig,ll)/100000.)**RKAPPA
	253	! contenu thermcell_env : zu(ig,ll)=pu(ig,ll)
	254	! contenu thermcell_env : zv(ig,ll)=pv(ig,ll)
	255	! contenu thermcell_env : ztv(ig,ll)=zh(ig,ll)/zpspsk(ig,ll)
	256	! contenu thermcell_env : ztv(ig,ll)=ztv(ig,ll)(1.+RETV(zo(ig,ll))-zl(ig,ll))
	257	! contenu thermcell_env : zthl(ig,ll)=pt(ig,ll)/zpspsk(ig,ll)
	258	! contenu thermcell_env : enddo
	259	! contenu thermcell_env : enddo
	260
	261	do l=1,nlay
	262	do ig=1,ngrid
	263	zl(ig,l)=0.
	264	zu(ig,l)=puwind(ig,l)
	265	zv(ig,l)=pvwind(ig,l)
	266	ztemp_env(ig,l)=ptemp_env(ig,l)
	267	zpspsk(ig,l)=(pplay(ig,l)/100000.)**RKAPPA
	268	ztv(ig,l)=ztemp_env(ig,l)/zpspsk(ig,l)
	269	ztv(ig,l)=ztv(ig,l)(1.+RETVpo_env(ig,l))
	270	zthl(ig,l)=ptemp(ig,l)/zpspsk(ig,l)
	271	mask(ig,l)=.true.
	272	enddo
	273	enddo
	274	call thermcell_qsat(ngrid*nlay,mask,pplev,ptemp_env,p_o,zqsat)
	275
	276	endif
[878]	277
[938]	278	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
[878]	279
	280	!------------------------------------------------------------------------
	281	! --------------------
	282	!
	283	!
	284	! + + + + + + + + + + +
	285	!
	286	!
	287	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
	288	! wh,wt,wo ...
	289	!
[4678]	290	! + + + + + + + + + + + zh,zu,zv,z_o,rho
[878]	291	!
	292	!
	293	! -------------------- zlev(1)
	294	! \\\\\\\\\\\\\\\\\\\\
	295	!
	296	!
	297
	298	!-----------------------------------------------------------------------
	299	! Calcul des altitudes des couches
	300	!-----------------------------------------------------------------------
	301
	302	do l=2,nlay
	303	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
	304	enddo
[4089]	305	zlev(:,1)=0.
	306	zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG
[878]	307	do l=1,nlay
	308	zlay(:,l)=pphi(:,l)/RG
	309	enddo
	310	do l=1,nlay
	311	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
	312	enddo
	313
	314	!-----------------------------------------------------------------------
[4089]	315	! Calcul des densites et masses
[878]	316	!-----------------------------------------------------------------------
	317
[4089]	318	rho(:,:)=pplay(:,:)/(zpspsk(:,:)RDztv(:,:))
	319	if (prt_level.ge.10) write(lunout,*) 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
[972]	320	rhobarz(:,1)=rho(:,1)
[878]	321	do l=2,nlay
	322	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
	323	enddo
	324	do l=1,nlay
	325	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
	326	enddo
[938]	327	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
[878]	328
	329	!------------------------------------------------------------------
	330	!
	331	! /\|\
	332	! -------- \| F_k+1 -------
	333	! ----> D_k
	334	! /\|\ <---- E_k , A_k
	335	! -------- \| F_k ---------
	336	! ----> D_k-1
	337	! <---- E_k-1 , A_k-1
	338	!
	339	!
	340	!
	341	!
	342	!
	343	! ---------------------------
	344	!
	345	! ----- F_lmax+1=0 ---------- \
	346	! lmax (zmax) \|
	347	! --------------------------- \|
	348	! \|
	349	! --------------------------- \|
	350	! \|
	351	! --------------------------- \|
	352	! \|
	353	! --------------------------- \|
	354	! \|
	355	! --------------------------- \|
	356	! \| E
	357	! --------------------------- \| D
	358	! \|
	359	! --------------------------- \|
	360	! \|
	361	! --------------------------- \ \|
	362	! lalim \| \|
	363	! --------------------------- \| \|
	364	! \| \|
	365	! --------------------------- \| \|
	366	! \| A \|
	367	! --------------------------- \| \|
	368	! \| \|
	369	! --------------------------- \| \|
	370	! lmin (=1 pour le moment) \| \|
	371	! ----- F_lmin=0 ------------ / /
	372	!
	373	! ---------------------------
	374	! //////////////////////////
	375	!
	376	!
	377	!=============================================================================
	378	! Calculs initiaux ne faisant pas intervenir les changements de phase
	379	!=============================================================================
	380
	381	!------------------------------------------------------------------
[1403]	382	! 1. alim_star est le profil vertical de l'alimentation a la base du
	383	! panache thermique, calcule a partir de la flotabilite de l'air sec
[878]	384	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
	385	!------------------------------------------------------------------
	386	!
	387	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
[1403]	388	lmin=1
[878]	389
	390	!-----------------------------------------------------------------------------
	391	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
	392	! panache sec conservatif (e=d=0) alimente selon alim_star
	393	! Il s'agit d'un calcul de type CAPE
[1403]	394	! zmax_sec est utilise pour determiner la geometrie du thermique.
[878]	395	!------------------------------------------------------------------------------
[1403]	396	!---------------------------------------------------------------------------------
	397	!calcul du melange et des variables dans le thermique
	398	!--------------------------------------------------------------------------------
[878]	399	!
[1403]	400	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
[878]	401
[3451]	402	!=====================================================================
	403	! Old version of thermcell_plume in thermcell_plume_6A.F90
	404	! It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
	405	! to the 5B and 6A versions used for CMIP5 and CMIP6.
	406	! The latest was previously named thermcellV1_plume.
	407	! The new thermcell_plume is a clean version (removing obsolete
	408	! options) of thermcell_plume_6A.
	409	! The 3 versions are controled by
	410	! flag_thermals_ed <= 9 thermcell_plume_6A
	411	! <= 19 thermcell_plume_5B
	412	! else thermcell_plume (default 20 for convergence with 6A)
	413	! Fredho
	414	!=====================================================================
[878]	415
[1403]	416	if (iflag_thermals_ed<=9) then
	417	! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
[4678]	418	CALL thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,p_o,zl,rhobarz,&
[1403]	419	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	420	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	421	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
	422	& ,lev_out,lunout1,igout)
[878]	423
[3451]	424	elseif (iflag_thermals_ed<=19) then
[1403]	425	! print*,'THERM RIO et al 2010, version d Arnaud'
[4678]	426	CALL thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,p_o,zl,rhobarz,&
[1403]	427	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	428	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	429	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
	430	& ,lev_out,lunout1,igout)
[3451]	431	else
[4678]	432	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,p_o,zl,rhobarz,&
[3451]	433	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	434	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	435	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
	436	& ,lev_out,lunout1,igout)
[1403]	437	endif
[878]	438
[972]	439	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
	440
[4678]	441	call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
	442	call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
[878]	443
[938]	444	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
	445	if (prt_level.ge.10) then
[972]	446	write(lunout1,*) 'Dans thermcell_main 2'
	447	write(lunout1,*) 'lmin ',lmin(igout)
	448	write(lunout1,*) 'lalim ',lalim(igout)
	449	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
	450	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
[878]	451	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
	452	endif
	453
	454	!-------------------------------------------------------------------------------
	455	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
	456	!-------------------------------------------------------------------------------
	457	!
[4843]	458	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lcong,lintercong,lmix,zw2, &
	459	& zlev,lmax,zmax,zmax0,zmix,wmax,zcong)
[1403]	460	! Attention, w2 est transforme en sa racine carree dans cette routine
[4143]	461	! Le probleme vient du fait que linter et lmix sont souvent egaux a 1.
[1403]	462	wmax_tmp=0.
	463	do l=1,nlay
	464	wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
	465	enddo
	466	! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
[878]	467
	468
[1403]	469
[4678]	470	call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
	471	call test_ltherm(ngrid,nlay,pplay,lmin ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
	472	call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
	473	call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
[878]	474
[938]	475	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
[878]	476
	477	!-------------------------------------------------------------------------------
	478	! Fermeture,determination de f
	479	!-------------------------------------------------------------------------------
[1026]	480	!
[1403]	481	!
	482	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
[4094]	483	& lalim,lmin,zmax_sec,wmax_sec)
[878]	484
[1998]	485
[4678]	486	call test_ltherm(ngrid,nlay,pplay,lmin,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
	487	call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
[1403]	488
	489	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
	490	if (prt_level.ge.10) then
	491	write(lunout1,*) 'Dans thermcell_main 1b'
	492	write(lunout1,*) 'lmin ',lmin(igout)
	493	write(lunout1,*) 'lalim ',lalim(igout)
	494	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
	495	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
	496	& ,l=1,lalim(igout)+4)
	497	endif
	498
	499
	500
	501
	502	! Choix de la fonction d'alimentation utilisee pour la fermeture.
	503	! Apparemment sans importance
	504	alim_star_clos(:,:)=alim_star(:,:)
	505	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
[1998]	506	!
	507	!CR Appel de la fermeture seche
	508	if (iflag_thermals_closure.eq.1) then
[1403]	509
[4094]	510	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
	511	& zlev,lalim,alim_star_clos,zmax_sec,wmax_sec,f)
[878]	512
[1403]	513	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	514	! Appel avec les zmax et wmax tenant compte de la condensation
	515	! Semble moins bien marcher
[1998]	516	else if (iflag_thermals_closure.eq.2) then
	517
	518	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
[4094]	519	& zlev,lalim,alim_star,zmax,wmax,f)
[1998]	520
[4094]	521
[1998]	522	endif
	523
[1403]	524	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	525
[938]	526	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
[878]	527
[972]	528	if (tau_thermals>1.) then
	529	lambda=exp(-ptimestep/tau_thermals)
	530	f0=(1.-lambda)f+lambdaf0
	531	else
	532	f0=f
	533	endif
	534
	535	! Test valable seulement en 1D mais pas genant
	536	if (.not. (f0(1).ge.0.) ) then
[1403]	537	abort_message = '.not. (f0(1).ge.0.)'
[2311]	538	CALL abort_physic (modname,abort_message,1)
[972]	539	endif
	540
[878]	541	!-------------------------------------------------------------------------------
	542	!deduction des flux
	543
[972]	544	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
[878]	545	& lalim,lmax,alim_star, &
	546	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
[972]	547	& detr,zqla,lev_out,lunout1,igout)
[4318]	548
[972]	549	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
[878]	550
[938]	551	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
[4678]	552	call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
	553	call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
[878]	554
	555	!------------------------------------------------------------------
[972]	556	! On ne prend pas directement les profils issus des calculs precedents
	557	! mais on s'autorise genereusement une relaxation vers ceci avec
	558	! une constante de temps tau_thermals (typiquement 1800s).
	559	!------------------------------------------------------------------
[878]	560
[972]	561	if (tau_thermals>1.) then
	562	lambda=exp(-ptimestep/tau_thermals)
	563	fm0=(1.-lambda)fm+lambdafm0
	564	entr0=(1.-lambda)entr+lambdaentr0
[1403]	565	detr0=(1.-lambda)detr+lambdadetr0
[878]	566	else
	567	fm0=fm
	568	entr0=entr
	569	detr0=detr
	570	endif
	571
[4396]	572	!------------------------------------------------------------------
	573	! Calcul de la fraction de l'ascendance
	574	!------------------------------------------------------------------
	575	do ig=1,ngrid
	576	fraca(ig,1)=0.
	577	fraca(ig,nlay+1)=0.
	578	enddo
	579	do l=2,nlay
	580	do ig=1,ngrid
	581	if (zw2(ig,l).gt.1.e-10) then
	582	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
	583	else
	584	fraca(ig,l)=0.
	585	endif
	586	enddo
	587	enddo
	588
[972]	589	!c------------------------------------------------------------------
	590	! calcul du transport vertical
	591	!------------------------------------------------------------------
[4381]	592	IF (iflag_thermals_down .GT. 0) THEN
[4438]	593	if (debut) print*,'WARNING !!! routine thermcell_down en cours de developpement'
[4441]	594	entrdn=fact_thermals_down*detr0
	595	detrdn=fact_thermals_down*entr0
[4413]	596	! we want to transport potential temperature, total water and momentum
	597	CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zthl,zdthladj)
[4678]	598	CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,p_o,pdoadj)
[4413]	599	CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zu,pduadj)
	600	CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zv,pdvadj)
[4396]	601	ELSE
[4381]	602	!--------------------------------------------------------------
[972]	603
[4690]	604	! Temperature potentielle liquide effectivement mélangée par les thermiques
	605	do ll=1,nlay
	606	do ig=1,ngrid
	607	zthl(ig,ll)=ptemp(ig,ll)/zpspsk(ig,ll)
	608	enddo
	609	enddo
[4396]	610	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
	611	& zthl,zdthladj,zta,lev_out)
[4678]	612
	613	do ll=1,nlay
	614	do ig=1,ngrid
	615	z_o(ig,ll)=p_o(ig,ll)
	616	enddo
	617	enddo
[4396]	618	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
[4678]	619	& z_o,pdoadj,z_oa,lev_out)
[878]	620
[4089]	621	#ifdef ISO
[4143]	622	! C Risi: on utilise directement la meme routine
	623	do ixt=1,ntiso
[4089]	624	do ll=1,nlay
	625	DO ig=1,ngrid
	626	xtpo_tmp(ig,ll)=xtpo(ixt,ig,ll)
	627	xtzo_tmp(ig,ll)=xtzo(ixt,ig,ll)
	628	enddo
	629	enddo
	630	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
	631	& xtpo_tmp,xtpdoadj_tmp,xtzo_tmp,lev_out)
	632	do ll=1,nlay
	633	DO ig=1,ngrid
	634	xtpdoadj(ixt,ig,ll)=xtpdoadj_tmp(ig,ll)
	635	enddo
	636	enddo
[4143]	637	enddo
[4089]	638	#endif
	639
	640	#ifdef ISO
	641	#ifdef ISOVERIF
	642	DO ll=1,nlay
	643	DO ig=1,ngrid
	644	if (iso_eau.gt.0) then
	645	call iso_verif_egalite(xtpo(iso_eau,ig,ll), &
[4680]	646	& p_o(ig,ll),'thermcell_main 594')
[4089]	647	call iso_verif_egalite(xtpdoadj(iso_eau,ig,ll), &
	648	& pdoadj(ig,ll),'thermcell_main 596')
	649	endif
	650	if (iso_HDO.gt.0) then
	651	call iso_verif_aberrant_encadre(xtpo(iso_hdo,ig,ll) &
[4678]	652	& /p_o(ig,ll),'thermcell_main 610')
[4089]	653	endif
	654	enddo
	655	enddo !DO ll=1,nlay
	656	write(,) 'thermcell_main 600 tmp: apres thermcell_dq'
	657	#endif
	658	#endif
	659
	660
[883]	661	!------------------------------------------------------------------
	662	! calcul du transport vertical du moment horizontal
	663	!------------------------------------------------------------------
[878]	664
[972]	665	!IM 090508
[1738]	666	if (dvdq == 0 ) then
[883]	667
[878]	668	! Calcul du transport de V tenant compte d'echange par gradient
	669	! de pression horizontal avec l'environnement
	670
	671	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
[1738]	672	! & ,fraca*dvdq,zmax &
	673	& ,fraca,zmax &
[972]	674	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
[1403]	675
[878]	676	else
	677
	678	! calcul purement conservatif pour le transport de V
[1738]	679	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
[878]	680	& ,zu,pduadj,zua,lev_out)
[1738]	681	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
[878]	682	& ,zv,pdvadj,zva,lev_out)
[1738]	683
[878]	684	endif
[4396]	685	ENDIF
[878]	686
	687	! print*,'13 OK convect8'
	688	do l=1,nlay
	689	do ig=1,ngrid
	690	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
	691	enddo
	692	enddo
	693
[972]	694	if (prt_level.ge.1) print*,'14 OK convect8'
[878]	695	!------------------------------------------------------------------
	696	! Calculs de diagnostiques pour les sorties
	697	!------------------------------------------------------------------
	698	!calcul de fraca pour les sorties
	699
	700	if (sorties) then
[972]	701	if (prt_level.ge.1) print*,'14a OK convect8'
[878]	702	! calcul du niveau de condensation
	703	! initialisation
	704	do ig=1,ngrid
[879]	705	nivcon(ig)=0
[878]	706	zcon(ig)=0.
	707	enddo
	708	!nouveau calcul
	709	do ig=1,ngrid
[4690]	710	! WARNING !!! verifier que c'est bien ztemp_env qu'on veut là
	711	CHI=ztemp_env(ig,1)/(1669.0-122.0*z_o(ig,1)/zqsat(ig,1)-ztemp_env(ig,1))
[4678]	712	pcon(ig)=pplay(ig,1)(z_o(ig,1)/zqsat(ig,1))*CHI
[878]	713	enddo
[1403]	714	!IM do k=1,nlay
	715	do k=1,nlay-1
[878]	716	do ig=1,ngrid
	717	if ((pcon(ig).le.pplay(ig,k)) &
	718	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
	719	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
	720	endif
	721	enddo
	722	enddo
[1403]	723	!IM
[1494]	724	ierr=0
[1403]	725	do ig=1,ngrid
	726	if (pcon(ig).le.pplay(ig,nlay)) then
	727	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
[1494]	728	ierr=1
	729	endif
	730	enddo
[4843]	731	! if (ierr==1) then
	732	! abort_message = 'thermcellV0_main: les thermiques vont trop haut '
	733	! CALL abort_physic (modname,abort_message,1)
	734	! endif
[1494]	735
[972]	736	if (prt_level.ge.1) print*,'14b OK convect8'
[878]	737	do k=nlay,1,-1
	738	do ig=1,ngrid
	739	if (zqla(ig,k).gt.1e-10) then
	740	nivcon(ig)=k
	741	zcon(ig)=zlev(ig,k)
	742	endif
	743	enddo
	744	enddo
[972]	745	if (prt_level.ge.1) print*,'14c OK convect8'
[878]	746	!calcul des moments
	747	!initialisation
	748	do l=1,nlay
	749	do ig=1,ngrid
	750	q2(ig,l)=0.
	751	wth2(ig,l)=0.
	752	wth3(ig,l)=0.
	753	ratqscth(ig,l)=0.
	754	ratqsdiff(ig,l)=0.
	755	enddo
	756	enddo
[972]	757	if (prt_level.ge.1) print*,'14d OK convect8'
[1146]	758	if (prt_level.ge.10)write(lunout,*) &
	759	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
[878]	760	do l=1,nlay
	761	do ig=1,ngrid
	762	zf=fraca(ig,l)
	763	zf2=zf/(1.-zf)
[972]	764	!
[1403]	765	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
[972]	766	if(zw2(ig,l).gt.1.e-10) then
	767	wth2(ig,l)=zf2(zw2(ig,l))*2
	768	else
	769	wth2(ig,l)=0.
	770	endif
[878]	771	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
	772	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
[4678]	773	q2(ig,l)=zf2(zqta(ig,l)1000.-p_o(ig,l)1000.)*2
	774	!test: on calcul q2/p_o=ratqsc
	775	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(p_o(ig,l)*1000.))
[878]	776	enddo
	777	enddo
[1403]	778	!calcul des flux: q, thetal et thetav
	779	do l=1,nlay
	780	do ig=1,ngrid
[4678]	781	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-p_o(ig,l)1000.)
[1403]	782	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
	783	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
	784	enddo
[879]	785	enddo
[1638]	786
[878]	787	!calcul du ratqscdiff
[972]	788	if (prt_level.ge.1) print*,'14e OK convect8'
[878]	789	var=0.
	790	vardiff=0.
	791	ratqsdiff(:,:)=0.
[1494]	792
[4089]	793	do l=1,nlay
[1494]	794	do ig=1,ngrid
	795	if (l<=lalim(ig)) then
[878]	796	var=var+alim_star(ig,l)zqta(ig,l)1000.
[1494]	797	endif
[878]	798	enddo
	799	enddo
[1494]	800
[972]	801	if (prt_level.ge.1) print*,'14f OK convect8'
[1494]	802
[4089]	803	do l=1,nlay
[1494]	804	do ig=1,ngrid
	805	if (l<=lalim(ig)) then
	806	zf=fraca(ig,l)
	807	zf2=zf/(1.-zf)
	808	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
	809	endif
	810	enddo
[878]	811	enddo
[1494]	812
[972]	813	if (prt_level.ge.1) print*,'14g OK convect8'
[4089]	814	do l=1,nlay
	815	do ig=1,ngrid
[4678]	816	ratqsdiff(ig,l)=sqrt(vardiff)/(p_o(ig,l)*1000.)
[4089]	817	enddo
	818	enddo
[878]	819	endif
	820
[938]	821	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
[878]	822
[4692]	823	!print*,'thermcell_main fin'
[4094]	824	RETURN
[4089]	825	end subroutine thermcell_main
[878]	826
[4089]	827	!=============================================================================
	828	!/////////////////////////////////////////////////////////////////////////////
	829	!=============================================================================
[4678]	830	subroutine test_ltherm(ngrid,nlay,pplay,long,ztv,p_o,ztva, & ! in
[4089]	831	& zqla,f_star,zw2,comment) ! in
	832	!=============================================================================
[4590]	833	USE lmdz_thermcell_ini, ONLY: prt_level
[938]	834	IMPLICIT NONE
[878]	835
[4089]	836	integer i, k, ngrid,nlay
[4678]	837	real, intent(in), dimension(ngrid,nlay) :: pplay,ztv,p_o,ztva,zqla
[4089]	838	real, intent(in), dimension(ngrid,nlay) :: f_star,zw2
	839	integer, intent(in), dimension(ngrid) :: long
[878]	840	real seuil
	841	character*21 comment
[4133]	842
[4089]	843	seuil=0.25
[878]	844
[938]	845	if (prt_level.ge.1) THEN
	846	print*,'WARNING !!! TEST ',comment
	847	endif
[879]	848	return
	849
[878]	850	! test sur la hauteur des thermiques ...
[4089]	851	do i=1,ngrid
[972]	852	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
	853	if (prt_level.ge.10) then
[878]	854	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
	855	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
[4089]	856	do k=1,nlay
[4678]	857	write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000p_o(i,k),ztva(i,k),1000zqla(i,k),f_star(i,k),zw2(i,k)
[878]	858	enddo
[972]	859	endif
[878]	860	enddo
	861
	862
	863	return
	864	end
	865
[4089]	866	! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
	867	! On transporte pbl_tke pour donner therm_tke
[4143]	868	! Copie conforme de la subroutine DTKE dans physiq.F ecrite par Frederic Hourdin
[4089]	869
	870	!=======================================================================
	871	!///////////////////////////////////////////////////////////////////////
	872	!=======================================================================
	873
	874	subroutine thermcell_tke_transport( &
	875	& ngrid,nlay,ptimestep,fm0,entr0,rg,pplev, & ! in
	876	& therm_tke_max) ! out
[4590]	877	USE lmdz_thermcell_ini, ONLY: prt_level
[1638]	878	implicit none
	879
	880	!=======================================================================
	881	!
	882	! Calcul du transport verticale dans la couche limite en presence
	883	! de "thermiques" explicitement representes
	884	! calcul du dq/dt une fois qu'on connait les ascendances
	885	!
	886	!=======================================================================
	887
[4089]	888	integer ngrid,nlay
[1638]	889
[4089]	890	real, intent(in) :: ptimestep
	891	real, intent(in), dimension(ngrid,nlay+1) :: fm0,pplev
	892	real, intent(in), dimension(ngrid,nlay) :: entr0
	893	real, intent(in) :: rg
	894	real, intent(out), dimension(ngrid,nlay) :: therm_tke_max
	895
[1638]	896	real detr0(ngrid,nlay)
[4089]	897	real masse0(ngrid,nlay)
[1638]	898	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	899	real entr(ngrid,nlay)
	900	real q(ngrid,nlay)
	901	integer lev_out ! niveau pour les print
	902
	903	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
	904	integer ig,k
	905
	906
	907	lev_out=0
	908
	909
	910	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
	911
	912	! calcul du detrainement
	913	do k=1,nlay
	914	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
	915	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
	916	enddo
	917
	918
	919	! Decalage vertical des entrainements et detrainements.
	920	masse(:,1)=0.5*masse0(:,1)
	921	entr(:,1)=0.5*entr0(:,1)
	922	detr(:,1)=0.5*detr0(:,1)
	923	fm(:,1)=0.
	924	do k=1,nlay-1
	925	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
	926	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
	927	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
	928	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
	929	enddo
	930	fm(:,nlay+1)=0.
	931
	932
	933	q(:,:)=therm_tke_max(:,:)
	934	!!! nrlmd le 16/09/2010
	935	do ig=1,ngrid
	936	qa(ig,1)=q(ig,1)
	937	enddo
	938	!!!
	939
	940	if (1==1) then
	941	do k=2,nlay
	942	do ig=1,ngrid
	943	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
	944	& 1.e-5*masse(ig,k)) then
	945	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
	946	& /(fm(ig,k+1)+detr(ig,k))
	947	else
	948	qa(ig,k)=q(ig,k)
	949	endif
	950	if (qa(ig,k).lt.0.) then
	951	! print*,'qa<0!!!'
	952	endif
	953	if (q(ig,k).lt.0.) then
	954	! print*,'q<0!!!'
	955	endif
	956	enddo
	957	enddo
	958
	959	! Calcul du flux subsident
	960
	961	do k=2,nlay
	962	do ig=1,ngrid
	963	wqd(ig,k)=fm(ig,k)*q(ig,k)
	964	if (wqd(ig,k).lt.0.) then
	965	! print*,'wqd<0!!!'
	966	endif
	967	enddo
	968	enddo
	969	do ig=1,ngrid
	970	wqd(ig,1)=0.
	971	wqd(ig,nlay+1)=0.
	972	enddo
	973
	974	! Calcul des tendances
	975	do k=1,nlay
	976	do ig=1,ngrid
	977	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
	978	& -wqd(ig,k)+wqd(ig,k+1)) &
	979	& *ptimestep/masse(ig,k)
	980	enddo
	981	enddo
	982
	983	endif
	984
	985	therm_tke_max(:,:)=q(:,:)
	986
	987	return
	988	!!! fin nrlmd le 10/04/2012
	989	end
	990
[4590]	991	END MODULE lmdz_thermcell_main

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

Download in other formats:

Original Format