Context Navigation

thermcell_main.F90 @ 3603

Last change on this file since 3603 was 2387, checked in by fhourdin, 9 years ago
Poursuite du travail sur les thermiques. Extraction d'une routine calculant l'alimentation laterale à la base. + une correction de Jean-Yves sur des initialisations maquantes dans thermcell_main (mais qui sont maintenant dans thermcell_alp).
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.0 KB

Rev	Line
[878]	1	!
[1403]	2	! $Id: thermcell_main.F90 2387 2015-11-07 09:27:40Z musat $
[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
	442	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
[878]	443
[1403]	444	! Gestion temporaire de plusieurs appels à thermcell_plume au travers
	445	! de la variable iflag_thermals
[878]	446
[1403]	447	! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
	448	if (iflag_thermals_ed<=9) then
	449	! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
	450	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
	451	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	452	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	453	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
	454	& ,lev_out,lunout1,igout)
[878]	455
[1403]	456	elseif (iflag_thermals_ed>9) then
	457	! print*,'THERM RIO et al 2010, version d Arnaud'
	458	CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
	459	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
	460	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
	461	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
	462	& ,lev_out,lunout1,igout)
[878]	463
[1403]	464	endif
[878]	465
[972]	466	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
	467
[878]	468	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
	469	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
	470
[938]	471	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
	472	if (prt_level.ge.10) then
[972]	473	write(lunout1,*) 'Dans thermcell_main 2'
	474	write(lunout1,*) 'lmin ',lmin(igout)
	475	write(lunout1,*) 'lalim ',lalim(igout)
	476	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
	477	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
[878]	478	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
	479	endif
	480
	481	!-------------------------------------------------------------------------------
	482	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
	483	!-------------------------------------------------------------------------------
	484	!
	485	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
	486	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
[1403]	487	! Attention, w2 est transforme en sa racine carree dans cette routine
	488	! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
	489	wmax_tmp=0.
	490	do l=1,nlay
	491	wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
	492	enddo
	493	! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
[878]	494
	495
[1403]	496
[878]	497	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
	498	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
	499	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
	500	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
	501
[938]	502	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
[878]	503
	504	!-------------------------------------------------------------------------------
	505	! Fermeture,determination de f
	506	!-------------------------------------------------------------------------------
[1026]	507	!
[1403]	508	!
	509	!! write(lunout,*)'THERM NOUVEAU XXXXX'
	510	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
	511	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
[878]	512
[1998]	513
[1403]	514	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
	515	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
	516
	517	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
	518	if (prt_level.ge.10) then
	519	write(lunout1,*) 'Dans thermcell_main 1b'
	520	write(lunout1,*) 'lmin ',lmin(igout)
	521	write(lunout1,*) 'lalim ',lalim(igout)
	522	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
	523	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
	524	& ,l=1,lalim(igout)+4)
	525	endif
	526
	527
	528
	529
	530	! Choix de la fonction d'alimentation utilisee pour la fermeture.
	531	! Apparemment sans importance
	532	alim_star_clos(:,:)=alim_star(:,:)
	533	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
[1998]	534	!
	535	!CR Appel de la fermeture seche
	536	if (iflag_thermals_closure.eq.1) then
[1403]	537
[1496]	538	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
[1403]	539	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
[878]	540
[1403]	541	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	542	! Appel avec les zmax et wmax tenant compte de la condensation
	543	! Semble moins bien marcher
[1998]	544	else if (iflag_thermals_closure.eq.2) then
	545
	546	CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, &
	547	& zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
	548
	549	endif
	550
[1403]	551	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	552
[938]	553	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
[878]	554
[972]	555	if (tau_thermals>1.) then
	556	lambda=exp(-ptimestep/tau_thermals)
	557	f0=(1.-lambda)f+lambdaf0
	558	else
	559	f0=f
	560	endif
	561
	562	! Test valable seulement en 1D mais pas genant
	563	if (.not. (f0(1).ge.0.) ) then
[1403]	564	abort_message = '.not. (f0(1).ge.0.)'
[2311]	565	CALL abort_physic (modname,abort_message,1)
[972]	566	endif
	567
[878]	568	!-------------------------------------------------------------------------------
	569	!deduction des flux
	570	!-------------------------------------------------------------------------------
	571
[972]	572	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
[878]	573	& lalim,lmax,alim_star, &
	574	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
[972]	575	& detr,zqla,lev_out,lunout1,igout)
	576	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
[878]	577
[938]	578	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
[878]	579	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
	580	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
	581
	582	!------------------------------------------------------------------
[972]	583	! On ne prend pas directement les profils issus des calculs precedents
	584	! mais on s'autorise genereusement une relaxation vers ceci avec
	585	! une constante de temps tau_thermals (typiquement 1800s).
	586	!------------------------------------------------------------------
[878]	587
[972]	588	if (tau_thermals>1.) then
	589	lambda=exp(-ptimestep/tau_thermals)
	590	fm0=(1.-lambda)fm+lambdafm0
	591	entr0=(1.-lambda)entr+lambdaentr0
[1403]	592	detr0=(1.-lambda)detr+lambdadetr0
[878]	593	else
	594	fm0=fm
	595	entr0=entr
	596	detr0=detr
	597	endif
	598
[972]	599	!c------------------------------------------------------------------
	600	! calcul du transport vertical
	601	!------------------------------------------------------------------
	602
[1738]	603	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
[878]	604	& zthl,zdthladj,zta,lev_out)
[1738]	605	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, &
[878]	606	& po,pdoadj,zoa,lev_out)
	607
[883]	608	!------------------------------------------------------------------
	609	! Calcul de la fraction de l'ascendance
	610	!------------------------------------------------------------------
	611	do ig=1,klon
	612	fraca(ig,1)=0.
	613	fraca(ig,nlay+1)=0.
	614	enddo
	615	do l=2,nlay
	616	do ig=1,klon
	617	if (zw2(ig,l).gt.1.e-10) then
	618	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
	619	else
	620	fraca(ig,l)=0.
	621	endif
	622	enddo
	623	enddo
	624
	625	!------------------------------------------------------------------
	626	! calcul du transport vertical du moment horizontal
	627	!------------------------------------------------------------------
[878]	628
[972]	629	!IM 090508
[1738]	630	if (dvdq == 0 ) then
[883]	631
[878]	632	! Calcul du transport de V tenant compte d'echange par gradient
	633	! de pression horizontal avec l'environnement
	634
	635	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
[1738]	636	! & ,fraca*dvdq,zmax &
	637	& ,fraca,zmax &
[972]	638	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
[1403]	639
[878]	640	else
	641
	642	! calcul purement conservatif pour le transport de V
[1738]	643	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
[878]	644	& ,zu,pduadj,zua,lev_out)
[1738]	645	call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse &
[878]	646	& ,zv,pdvadj,zva,lev_out)
[1738]	647
[878]	648	endif
	649
	650	! print*,'13 OK convect8'
	651	do l=1,nlay
	652	do ig=1,ngrid
	653	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
	654	enddo
	655	enddo
	656
[972]	657	if (prt_level.ge.1) print*,'14 OK convect8'
[878]	658	!------------------------------------------------------------------
	659	! Calculs de diagnostiques pour les sorties
	660	!------------------------------------------------------------------
	661	!calcul de fraca pour les sorties
	662
	663	if (sorties) then
[972]	664	if (prt_level.ge.1) print*,'14a OK convect8'
[878]	665	! calcul du niveau de condensation
	666	! initialisation
	667	do ig=1,ngrid
[879]	668	nivcon(ig)=0
[878]	669	zcon(ig)=0.
	670	enddo
	671	!nouveau calcul
	672	do ig=1,ngrid
	673	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
	674	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
	675	enddo
[1403]	676	!IM do k=1,nlay
	677	do k=1,nlay-1
[878]	678	do ig=1,ngrid
	679	if ((pcon(ig).le.pplay(ig,k)) &
	680	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
	681	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
	682	endif
	683	enddo
	684	enddo
[1403]	685	!IM
[1494]	686	ierr=0
[1403]	687	do ig=1,ngrid
	688	if (pcon(ig).le.pplay(ig,nlay)) then
	689	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
[1494]	690	ierr=1
	691	endif
	692	enddo
	693	if (ierr==1) then
[1403]	694	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
[2311]	695	CALL abort_physic (modname,abort_message,1)
[1494]	696	endif
	697
[972]	698	if (prt_level.ge.1) print*,'14b OK convect8'
[878]	699	do k=nlay,1,-1
	700	do ig=1,ngrid
	701	if (zqla(ig,k).gt.1e-10) then
	702	nivcon(ig)=k
	703	zcon(ig)=zlev(ig,k)
	704	endif
	705	enddo
	706	enddo
[972]	707	if (prt_level.ge.1) print*,'14c OK convect8'
[878]	708	!calcul des moments
	709	!initialisation
	710	do l=1,nlay
	711	do ig=1,ngrid
	712	q2(ig,l)=0.
	713	wth2(ig,l)=0.
	714	wth3(ig,l)=0.
	715	ratqscth(ig,l)=0.
	716	ratqsdiff(ig,l)=0.
	717	enddo
	718	enddo
[972]	719	if (prt_level.ge.1) print*,'14d OK convect8'
[1146]	720	if (prt_level.ge.10)write(lunout,*) &
	721	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
[878]	722	do l=1,nlay
	723	do ig=1,ngrid
	724	zf=fraca(ig,l)
	725	zf2=zf/(1.-zf)
[972]	726	!
[1403]	727	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
[972]	728	if(zw2(ig,l).gt.1.e-10) then
	729	wth2(ig,l)=zf2(zw2(ig,l))*2
	730	else
	731	wth2(ig,l)=0.
	732	endif
[878]	733	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
	734	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
	735	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
	736	!test: on calcul q2/po=ratqsc
	737	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
	738	enddo
	739	enddo
[1403]	740	!calcul des flux: q, thetal et thetav
	741	do l=1,nlay
	742	do ig=1,ngrid
	743	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
	744	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
	745	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
	746	enddo
[879]	747	enddo
[972]	748	!
[2384]	749	! $Id: thermcell_main.F90 2387 2015-11-07 09:27:40Z musat $
	750	!
	751	CALL thermcell_alp(ngrid,nlay,ptimestep &
	752	& ,pplay,pplev &
	753	& ,fm0,entr0,lmax &
	754	& ,Ale_bl,Alp_bl,lalim_conv,wght_th &
	755	& ,zw2,fraca &
	756	!!! necessire en plus
	757	& ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
[1638]	758	!!! nrlmd le 10/04/2012
[2384]	759	& ,pbl_tke,pctsrf,omega,airephy &
	760	& ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
	761	& ,n2,s2,ale_bl_stat &
	762	& ,therm_tke_max,env_tke_max &
	763	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
	764	& ,alp_bl_conv,alp_bl_stat &
[1638]	765	!!! fin nrlmd le 10/04/2012
[2384]	766	& )
[1638]	767
[1494]	768
[1403]	769
[878]	770	!calcul du ratqscdiff
[972]	771	if (prt_level.ge.1) print*,'14e OK convect8'
[878]	772	var=0.
	773	vardiff=0.
	774	ratqsdiff(:,:)=0.
[1494]	775
	776	do l=1,klev
	777	do ig=1,ngrid
	778	if (l<=lalim(ig)) then
[878]	779	var=var+alim_star(ig,l)zqta(ig,l)1000.
[1494]	780	endif
[878]	781	enddo
	782	enddo
[1494]	783
[972]	784	if (prt_level.ge.1) print*,'14f OK convect8'
[1494]	785
	786	do l=1,klev
	787	do ig=1,ngrid
	788	if (l<=lalim(ig)) then
	789	zf=fraca(ig,l)
	790	zf2=zf/(1.-zf)
	791	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
	792	endif
	793	enddo
[878]	794	enddo
[1494]	795
[972]	796	if (prt_level.ge.1) print*,'14g OK convect8'
[878]	797	do l=1,nlay
	798	do ig=1,ngrid
	799	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
	800	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
	801	enddo
	802	enddo
	803	!--------------------------------------------------------------------
	804	!
	805	!ecriture des fichiers sortie
[1494]	806	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
[878]	807
	808	endif
	809
[938]	810	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
[878]	811
	812	return
	813	end
	814
	815	!-----------------------------------------------------------------------------
	816
	817	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
[2311]	818	USE print_control_mod, ONLY: prt_level
[938]	819	IMPLICIT NONE
[878]	820
[938]	821	integer i, k, klon,klev
[878]	822	real pplev(klon,klev+1),pplay(klon,klev)
	823	real ztv(klon,klev)
	824	real po(klon,klev)
	825	real ztva(klon,klev)
	826	real zqla(klon,klev)
	827	real f_star(klon,klev)
	828	real zw2(klon,klev)
	829	integer long(klon)
	830	real seuil
	831	character*21 comment
	832
[938]	833	if (prt_level.ge.1) THEN
	834	print*,'WARNING !!! TEST ',comment
	835	endif
[879]	836	return
	837
[878]	838	! test sur la hauteur des thermiques ...
	839	do i=1,klon
[972]	840	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
	841	if (prt_level.ge.10) then
[878]	842	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
	843	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
	844	do k=1,klev
	845	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)
	846	enddo
[972]	847	endif
[878]	848	enddo
	849
	850
	851	return
	852	end
	853
[1638]	854	!!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP
	855	! On transporte pbl_tke pour donner therm_tke
	856	! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
	857	subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, &
	858	& rg,pplev,therm_tke_max)
[2311]	859	USE print_control_mod, ONLY: prt_level
[1638]	860	implicit none
	861
	862	!=======================================================================
	863	!
	864	! Calcul du transport verticale dans la couche limite en presence
	865	! de "thermiques" explicitement representes
	866	! calcul du dq/dt une fois qu'on connait les ascendances
	867	!
	868	!=======================================================================
	869
	870	integer ngrid,nlay,nsrf
	871
	872	real ptimestep
	873	real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
	874	real entr0(ngrid,nlay),rg
	875	real therm_tke_max(ngrid,nlay)
	876	real detr0(ngrid,nlay)
	877
	878
	879	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	880	real entr(ngrid,nlay)
	881	real q(ngrid,nlay)
	882	integer lev_out ! niveau pour les print
	883
	884	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
	885
	886	real zzm
	887
	888	integer ig,k
	889	integer isrf
	890
	891
	892	lev_out=0
	893
	894
	895	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
	896
	897	! calcul du detrainement
	898	do k=1,nlay
	899	detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
	900	masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
	901	enddo
	902
	903
	904	! Decalage vertical des entrainements et detrainements.
	905	masse(:,1)=0.5*masse0(:,1)
	906	entr(:,1)=0.5*entr0(:,1)
	907	detr(:,1)=0.5*detr0(:,1)
	908	fm(:,1)=0.
	909	do k=1,nlay-1
	910	masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
	911	entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
	912	detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
	913	fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
	914	enddo
	915	fm(:,nlay+1)=0.
	916
	917	!!! nrlmd le 16/09/2010
	918	! calcul de la valeur dans les ascendances
	919	! do ig=1,ngrid
	920	! qa(ig,1)=q(ig,1)
	921	! enddo
	922	!!!
	923
	924	!do isrf=1,nsrf
	925
	926	! q(:,:)=therm_tke(:,:,isrf)
	927	q(:,:)=therm_tke_max(:,:)
	928	!!! nrlmd le 16/09/2010
	929	do ig=1,ngrid
	930	qa(ig,1)=q(ig,1)
	931	enddo
	932	!!!
	933
	934	if (1==1) then
	935	do k=2,nlay
	936	do ig=1,ngrid
	937	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
	938	& 1.e-5*masse(ig,k)) then
	939	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
	940	& /(fm(ig,k+1)+detr(ig,k))
	941	else
	942	qa(ig,k)=q(ig,k)
	943	endif
	944	if (qa(ig,k).lt.0.) then
	945	! print*,'qa<0!!!'
	946	endif
	947	if (q(ig,k).lt.0.) then
	948	! print*,'q<0!!!'
	949	endif
	950	enddo
	951	enddo
	952
	953	! Calcul du flux subsident
	954
	955	do k=2,nlay
	956	do ig=1,ngrid
	957	wqd(ig,k)=fm(ig,k)*q(ig,k)
	958	if (wqd(ig,k).lt.0.) then
	959	! print*,'wqd<0!!!'
	960	endif
	961	enddo
	962	enddo
	963	do ig=1,ngrid
	964	wqd(ig,1)=0.
	965	wqd(ig,nlay+1)=0.
	966	enddo
	967
	968	! Calcul des tendances
	969	do k=1,nlay
	970	do ig=1,ngrid
	971	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
	972	& -wqd(ig,k)+wqd(ig,k+1)) &
	973	& *ptimestep/masse(ig,k)
	974	enddo
	975	enddo
	976
	977	endif
	978
	979	therm_tke_max(:,:)=q(:,:)
	980
	981	return
	982	!!! fin nrlmd le 10/04/2012
	983	end
	984

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

Download in other formats:

Original Format