Context Navigation

thermcell_main.F90 @ 4413

Last change on this file since 4413 was 4413, checked in by evignon, 2 years ago
pour les mini-projets: controle des coefficients pour le detrainement et l'entrainement des downdrafts dans les .def
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.3 KB

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

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

Download in other formats:

Original Format