Context Navigation

lmdz_thermcell_main.F90 @ 4686

Last change on this file since 4686 was 4684, checked in by fhourdin, 14 months ago
Noms de variables plus explicites
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.4 KB

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

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

Download in other formats:

Original Format