Context Navigation

thermcell_main.F90 @ 4586

Last change on this file since 4586 was 4441, checked in by fhourdin, 17 months ago
Parametre controlant les downdrafts
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: 30.2 KB

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

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

Download in other formats:

Original Format