Context Navigation

thermcell_main.F90 @ 3543

Last change on this file since 3543 was 3451, checked in by fhourdin, 6 years ago

Saving old version of thermcell_plume in thermcell_plume_6A.F90
It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
to the 5B and 6A versions used for CMIP5 and CMIP6.
The latest was previously named thermcellV1_plume.
The new thermcell_plume is a clean version (removing obsolete
options) of thermcell_plume_6A.
The 3 versions are controled by
flag_thermals_ed <= 9 thermcell_plume_6A

<= 19 thermcell_plume_5B
else thermcell_plume (default 20 for convergence with 6A)

Fredho

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: 32.7 KB

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

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

Download in other formats:

Original Format