Context Navigation

lmdz_thermcell_main.F90 @ 4590

Last change on this file since 4590 was 4590, checked in by fhourdin, 11 months ago
Passage des thermiques a la nouvelle norme.
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: 31.0 KB

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

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

Download in other formats:

Original Format