Context Navigation

thermcell_main.F90 @ 4390

Last change on this file since 4390 was 4381, checked in by evignon, 21 months ago
pour les downdrafts, suppression de la cle de compilation pour protection sous flag iflag_thermals_down + stop si le flux compensatoire est vers le 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: 29.7 KB

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

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

Download in other formats:

Original Format