Context Navigation

lmdz_thermcell_main.F90 @ 5512

Last change on this file since 5512 was 5512, checked in by yann meurdesoif, 2 days ago

Implement GPU automatic port for :

Thermics
acama_gwd_rando
flott_gwd_rando

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

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

Download in other formats:

Original Format