Context Navigation

thermcell_main.F90 @ 4094

Last change on this file since 4094 was 4094, checked in by fhourdin, 2 years ago
Nettoyage thermiques (suite) et replay1d
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.0 KB

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

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

Download in other formats:

Original Format