Context Navigation

thermcell_main.F90 @ 1235

Last change on this file since 1235 was 1142, checked in by yann meurdesoif, 16 years ago

Portage vers Vargas

Y.M + E.M

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 27.0 KB

Rev	Line
[878]	1	!
	2	! $Header$
	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 &
[972]	10	& ,r_aspect,l_mix,tau_thermals &
[927]	11	& ,Ale_bl,Alp_bl,lalim_conv,wght_th &
[1026]	12	& ,zmax0, f0,zw2,fraca)
[878]	13
[972]	14	USE dimphy
[1026]	15	USE comgeomphy , ONLY:rlond,rlatd
[878]	16	IMPLICIT NONE
	17
	18	!=======================================================================
	19	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
	20	! Version du 09.02.07
	21	! Calcul du transport vertical dans la couche limite en presence
	22	! de "thermiques" explicitement representes avec processus nuageux
	23	!
	24	! Réécriture à partir d'un listing papier à Habas, le 14/02/00
	25	!
	26	! le thermique est supposé homogène et dissipé par mélange avec
	27	! son environnement. la longueur l_mix contrôle l'efficacité du
	28	! mélange
	29	!
	30	! Le calcul du transport des différentes espèces se fait en prenant
	31	! en compte:
	32	! 1. un flux de masse montant
	33	! 2. un flux de masse descendant
	34	! 3. un entrainement
	35	! 4. un detrainement
	36	!
	37	!=======================================================================
	38
	39	!-----------------------------------------------------------------------
	40	! declarations:
	41	! -------------
	42
	43	#include "dimensions.h"
	44	#include "YOMCST.h"
	45	#include "YOETHF.h"
	46	#include "FCTTRE.h"
[938]	47	#include "iniprint.h"
[878]	48
	49	! arguments:
	50	! ----------
	51
[972]	52	!IM 140508
	53	INTEGER itap
	54
	55	INTEGER ngrid,nlay,w2di
	56	real tau_thermals
[878]	57	real ptimestep,l_mix,r_aspect
	58	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
	59	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
	60	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
	61	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
	62	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	63	real pphi(ngrid,nlay)
	64
	65	! local:
	66	! ------
	67
[972]	68	integer icount
	69	data icount/0/
	70	save icount
[987]	71	!$OMP THREADPRIVATE(icount)
[972]	72
[883]	73	integer,save :: igout=1
[987]	74	!$OMP THREADPRIVATE(igout)
[938]	75	integer,save :: lunout1=6
[987]	76	!$OMP THREADPRIVATE(lunout1)
[883]	77	integer,save :: lev_out=10
[987]	78	!$OMP THREADPRIVATE(lev_out)
[878]	79
	80	INTEGER ig,k,l,ll
	81	real zsortie1d(klon)
	82	INTEGER lmax(klon),lmin(klon),lalim(klon)
	83	INTEGER lmix(klon)
[1026]	84	INTEGER lmix_bis(klon)
[878]	85	real linter(klon)
	86	real zmix(klon)
[1026]	87	real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1)
	88	! real fraca(klon,klev)
	89
[878]	90	real zmax_sec(klon)
	91	!on garde le zmax du pas de temps precedent
	92	real zmax0(klon)
[927]	93	!FH/IM save zmax0
[878]	94
[972]	95	real lambda
	96
[878]	97	real zlev(klon,klev+1),zlay(klon,klev)
	98	real deltaz(klon,klev)
[972]	99	REAL zh(klon,klev)
[878]	100	real zthl(klon,klev),zdthladj(klon,klev)
	101	REAL ztv(klon,klev)
	102	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
	103	real zl(klon,klev)
	104	real zsortie(klon,klev)
	105	real zva(klon,klev)
	106	real zua(klon,klev)
	107	real zoa(klon,klev)
	108
	109	real zta(klon,klev)
	110	real zha(klon,klev)
	111	real fraca(klon,klev+1)
	112	real zf,zf2
	113	real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
	114	real q2(klon,klev)
[972]	115	! FH probleme de dimensionnement avec l'allocation dynamique
	116	! common/comtherm/thetath2,wth2
[878]	117
	118	real ratqscth(klon,klev)
	119	real var
	120	real vardiff
	121	real ratqsdiff(klon,klev)
	122
	123	logical sorties
[972]	124	real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
[878]	125	real zpspsk(klon,klev)
	126
	127	real wmax(klon)
	128	real wmax_sec(klon)
[972]	129	real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
	130	real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
[878]	131
	132	real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
	133	!niveau de condensation
[879]	134	integer nivcon(klon)
[878]	135	real zcon(klon)
	136	REAL CHI
	137	real zcon2(klon)
	138	real pcon(klon)
	139	real zqsat(klon,klev)
	140	real zqsatth(klon,klev)
	141
	142	real f_star(klon,klev+1),entr_star(klon,klev)
	143	real detr_star(klon,klev)
	144	real alim_star_tot(klon),alim_star2(klon)
	145	real alim_star(klon,klev)
	146	real f(klon), f0(klon)
[927]	147	!FH/IM save f0
[878]	148	real zlevinter(klon)
	149	logical debut
	150	real seuil
	151
	152	!
[879]	153	!nouvelles variables pour la convection
	154	real Ale_bl(klon)
	155	real Alp_bl(klon)
	156	real alp_int(klon)
	157	real ale_int(klon)
	158	integer n_int(klon)
	159	real fm_tot(klon)
	160	real wght_th(klon,klev)
	161	integer lalim_conv(klon)
[926]	162	!v1d logical therm
	163	!v1d save therm
[878]	164
	165	character*2 str2
	166	character*10 str10
	167
	168	EXTERNAL SCOPY
	169	!
	170
	171	!-----------------------------------------------------------------------
	172	! initialisation:
	173	! ---------------
	174	!
	175
	176	seuil=0.25
	177
[972]	178	if (debut) then
	179	fm0=0.
	180	entr0=0.
	181	detr0=0.
[1026]	182
	183
	184	! #define wrgrads_thermcell
	185	#undef wrgrads_thermcell
	186	#ifdef wrgrads_thermcell
	187	! Initialisation des sorties grads pour les thermiques.
	188	! Pour l'instant en 1D sur le point igout.
	189	! Utilise par thermcell_out3d.h
	190	str10='therm'
	191	call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
	192	& rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1., &
	193	& ptimestep,str10,'therm ')
	194	#endif
	195
	196
	197
[972]	198	endif
	199
	200	fm=0. ; entr=0. ; detr=0.
	201
	202	icount=icount+1
	203
	204	!IM 090508 beg
	205	!print*,'====================================================================='
	206	!print*,'====================================================================='
	207	!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
	208	!print*,'====================================================================='
	209	!print*,'====================================================================='
	210	!IM 090508 end
	211
[938]	212	if (prt_level.ge.1) print*,'thermcell_main V4'
[878]	213
	214	sorties=.true.
	215	IF(ngrid.NE.klon) THEN
	216	PRINT*
	217	PRINT*,'STOP dans convadj'
	218	PRINT*,'ngrid =',ngrid
	219	PRINT*,'klon =',klon
	220	ENDIF
	221	!
	222	!Initialisation
	223	!
[972]	224	! IF (1.eq.0) THEN
	225	! do ig=1,klon
[927]	226	!FH/IM 130308 if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
[972]	227	! if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
	228	! f0(ig)=1.e-5
	229	! zmax0(ig)=40.
[926]	230	!v1d therm=.false.
[972]	231	! endif
	232	! enddo
	233	! ENDIF !(1.eq.0) THEN
[1134]	234	if (prt_level.ge.10)write(lunout,*) &
	235	& 'WARNING thermcell_main f0=max(f0,1.e-2)'
[972]	236	do ig=1,klon
	237	if (prt_level.ge.20) then
	238	print*,'th_main ig f0',ig,f0(ig)
[878]	239	endif
[972]	240	f0(ig)=max(f0(ig),1.e-2)
	241	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
	242	enddo
[878]	243
	244	!-----------------------------------------------------------------------
	245	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
	246	! --------------------------------------------------------------------
	247	!
	248	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
	249	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
	250
[938]	251	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
[878]	252
	253	!------------------------------------------------------------------------
	254	! --------------------
	255	!
	256	!
	257	! + + + + + + + + + + +
	258	!
	259	!
	260	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
	261	! wh,wt,wo ...
	262	!
	263	! + + + + + + + + + + + zh,zu,zv,zo,rho
	264	!
	265	!
	266	! -------------------- zlev(1)
	267	! \\\\\\\\\\\\\\\\\\\\
	268	!
	269	!
	270
	271	!-----------------------------------------------------------------------
	272	! Calcul des altitudes des couches
	273	!-----------------------------------------------------------------------
	274
	275	do l=2,nlay
	276	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
	277	enddo
	278	zlev(:,1)=0.
	279	zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
	280	do l=1,nlay
	281	zlay(:,l)=pphi(:,l)/RG
	282	enddo
	283	!calcul de l epaisseur des couches
	284	do l=1,nlay
	285	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
	286	enddo
	287
	288	! print*,'2 OK convect8'
	289	!-----------------------------------------------------------------------
	290	! Calcul des densites
	291	!-----------------------------------------------------------------------
	292
	293	do l=1,nlay
	294	rho(:,l)=pplay(:,l)/(zpspsk(:,l)RDztv(:,l))
	295	enddo
	296
[972]	297	!IM
[1134]	298	if (prt_level.ge.10)write(lunout,*) &
	299	& 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
[972]	300	rhobarz(:,1)=rho(:,1)
	301
[878]	302	do l=2,nlay
	303	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
	304	enddo
	305
	306	!calcul de la masse
	307	do l=1,nlay
	308	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
	309	enddo
	310
[938]	311	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
[878]	312
	313	!------------------------------------------------------------------
	314	!
	315	! /\|\
	316	! -------- \| F_k+1 -------
	317	! ----> D_k
	318	! /\|\ <---- E_k , A_k
	319	! -------- \| F_k ---------
	320	! ----> D_k-1
	321	! <---- E_k-1 , A_k-1
	322	!
	323	!
	324	!
	325	!
	326	!
	327	! ---------------------------
	328	!
	329	! ----- F_lmax+1=0 ---------- \
	330	! lmax (zmax) \|
	331	! --------------------------- \|
	332	! \|
	333	! --------------------------- \|
	334	! \|
	335	! --------------------------- \|
	336	! \|
	337	! --------------------------- \|
	338	! \|
	339	! --------------------------- \|
	340	! \| E
	341	! --------------------------- \| D
	342	! \|
	343	! --------------------------- \|
	344	! \|
	345	! --------------------------- \ \|
	346	! lalim \| \|
	347	! --------------------------- \| \|
	348	! \| \|
	349	! --------------------------- \| \|
	350	! \| A \|
	351	! --------------------------- \| \|
	352	! \| \|
	353	! --------------------------- \| \|
	354	! lmin (=1 pour le moment) \| \|
	355	! ----- F_lmin=0 ------------ / /
	356	!
	357	! ---------------------------
	358	! //////////////////////////
	359	!
	360	!
	361	!=============================================================================
	362	! Calculs initiaux ne faisant pas intervenir les changements de phase
	363	!=============================================================================
	364
	365	!------------------------------------------------------------------
	366	! 1. alim_star est le profil vertical de l'alimentation à la base du
	367	! panache thermique, calculé à partir de la flotabilité de l'air sec
	368	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
	369	!------------------------------------------------------------------
	370	!
	371	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
[972]	372	CALL thermcell_init(ngrid,nlay,ztv,zlay,zlev, &
[878]	373	& lalim,lmin,alim_star,alim_star_tot,lev_out)
	374
	375	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin ')
	376	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
	377
	378
[938]	379	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
	380	if (prt_level.ge.10) then
[972]	381	write(lunout1,*) 'Dans thermcell_main 1'
	382	write(lunout1,*) 'lmin ',lmin(igout)
	383	write(lunout1,*) 'lalim ',lalim(igout)
	384	write(lunout1,*) ' ig l alim_star thetav'
	385	write(lunout1,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
[878]	386	& ,ztv(igout,l),l=1,lalim(igout)+4)
	387	endif
	388
[938]	389	!v1d do ig=1,klon
[926]	390	!v1d if (alim_star(ig,1).gt.1.e-10) then
	391	!v1d therm=.true.
	392	!v1d endif
[938]	393	!v1d enddo
[878]	394	!-----------------------------------------------------------------------------
	395	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
	396	! panache sec conservatif (e=d=0) alimente selon alim_star
	397	! Il s'agit d'un calcul de type CAPE
	398	! zmax_sec est utilisé pour déterminer la géométrie du thermique.
	399	!------------------------------------------------------------------------------
	400	!
	401	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
	402	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
	403
	404	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
	405	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
	406
[938]	407	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
	408	if (prt_level.ge.10) then
[972]	409	write(lunout1,*) 'Dans thermcell_main 1b'
	410	write(lunout1,*) 'lmin ',lmin(igout)
	411	write(lunout1,*) 'lalim ',lalim(igout)
	412	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
	413	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
[878]	414	& ,l=1,lalim(igout)+4)
	415	endif
	416
	417
	418
	419	!---------------------------------------------------------------------------------
	420	!calcul du melange et des variables dans le thermique
	421	!--------------------------------------------------------------------------------
	422	!
[972]	423	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
	424	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
	425	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
[1026]	426	& zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot, &
[878]	427	& lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, &
[1026]	428	& ztla,zqla,zqta,zha,zw2,zw_est,zqsatth,lmix,lmix_bis,linter &
	429	& ,lev_out,lunout1,igout)
[972]	430	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
	431
[878]	432	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
	433	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
	434
[938]	435	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
	436	if (prt_level.ge.10) then
[972]	437	write(lunout1,*) 'Dans thermcell_main 2'
	438	write(lunout1,*) 'lmin ',lmin(igout)
	439	write(lunout1,*) 'lalim ',lalim(igout)
	440	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
	441	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
[878]	442	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
	443	endif
	444
	445	!-------------------------------------------------------------------------------
	446	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
	447	!-------------------------------------------------------------------------------
	448	!
	449	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
	450	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
	451
	452
	453	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
	454	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
	455	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
	456	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
	457
[938]	458	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
[878]	459
	460	!-------------------------------------------------------------------------------
	461	! Fermeture,determination de f
	462	!-------------------------------------------------------------------------------
[1026]	463	!
	464	!avant closure: on redéfinit lalim, alim_star_tot et alim_star
	465	! do ig=1,klon
	466	! do l=2,lalim(ig)
	467	! alim_star(ig,l)=entr_star(ig,l)
	468	! entr_star(ig,l)=0.
	469	! enddo
	470	! enddo
[878]	471
	472	CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
[1026]	473	& zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
[878]	474
[938]	475	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
[878]	476
[972]	477	if (tau_thermals>1.) then
	478	lambda=exp(-ptimestep/tau_thermals)
	479	f0=(1.-lambda)f+lambdaf0
	480	else
	481	f0=f
	482	endif
	483
	484	! Test valable seulement en 1D mais pas genant
	485	if (.not. (f0(1).ge.0.) ) then
[1142]	486	stop 'Dans thermcell_main'
[972]	487	endif
	488
[878]	489	!-------------------------------------------------------------------------------
	490	!deduction des flux
	491	!-------------------------------------------------------------------------------
	492
[972]	493	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
[878]	494	& lalim,lmax,alim_star, &
	495	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
[972]	496	& detr,zqla,lev_out,lunout1,igout)
	497	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
[878]	498
[938]	499	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
[878]	500	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
	501	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
	502
	503	!------------------------------------------------------------------
[972]	504	! On ne prend pas directement les profils issus des calculs precedents
	505	! mais on s'autorise genereusement une relaxation vers ceci avec
	506	! une constante de temps tau_thermals (typiquement 1800s).
	507	!------------------------------------------------------------------
[878]	508
[972]	509	if (tau_thermals>1.) then
	510	lambda=exp(-ptimestep/tau_thermals)
	511	fm0=(1.-lambda)fm+lambdafm0
	512	entr0=(1.-lambda)entr+lambdaentr0
	513	! detr0=(1.-lambda)detr+lambdadetr0
[878]	514	else
	515	fm0=fm
	516	entr0=entr
	517	detr0=detr
	518	endif
	519
[972]	520	!c------------------------------------------------------------------
	521	! calcul du transport vertical
	522	!------------------------------------------------------------------
	523
[878]	524	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
	525	& zthl,zdthladj,zta,lev_out)
	526	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
	527	& po,pdoadj,zoa,lev_out)
	528
[883]	529	!------------------------------------------------------------------
	530	! Calcul de la fraction de l'ascendance
	531	!------------------------------------------------------------------
	532	do ig=1,klon
	533	fraca(ig,1)=0.
	534	fraca(ig,nlay+1)=0.
	535	enddo
	536	do l=2,nlay
	537	do ig=1,klon
	538	if (zw2(ig,l).gt.1.e-10) then
	539	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
	540	else
	541	fraca(ig,l)=0.
	542	endif
	543	enddo
	544	enddo
	545
	546	!------------------------------------------------------------------
	547	! calcul du transport vertical du moment horizontal
	548	!------------------------------------------------------------------
[878]	549
[972]	550	!IM 090508
[883]	551	if (1.eq.1) then
[972]	552	!IM 070508 vers. _dq
	553	! if (1.eq.0) then
[883]	554
	555
[878]	556	! Calcul du transport de V tenant compte d'echange par gradient
	557	! de pression horizontal avec l'environnement
	558
	559	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
	560	& ,fraca,zmax &
[972]	561	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
	562	!IM 050508 & ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
[878]	563	else
	564
	565	! calcul purement conservatif pour le transport de V
	566	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
	567	& ,zu,pduadj,zua,lev_out)
	568	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
	569	& ,zv,pdvadj,zva,lev_out)
	570	endif
	571
	572	! print*,'13 OK convect8'
	573	do l=1,nlay
	574	do ig=1,ngrid
	575	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
	576	enddo
	577	enddo
	578
[972]	579	if (prt_level.ge.1) print*,'14 OK convect8'
[878]	580	!------------------------------------------------------------------
	581	! Calculs de diagnostiques pour les sorties
	582	!------------------------------------------------------------------
	583	!calcul de fraca pour les sorties
	584
	585	if (sorties) then
[972]	586	if (prt_level.ge.1) print*,'14a OK convect8'
[878]	587	! calcul du niveau de condensation
	588	! initialisation
	589	do ig=1,ngrid
[879]	590	nivcon(ig)=0
[878]	591	zcon(ig)=0.
	592	enddo
	593	!nouveau calcul
	594	do ig=1,ngrid
	595	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
	596	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
	597	enddo
	598	do k=1,nlay
	599	do ig=1,ngrid
	600	if ((pcon(ig).le.pplay(ig,k)) &
	601	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
	602	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
	603	endif
	604	enddo
	605	enddo
[972]	606	if (prt_level.ge.1) print*,'14b OK convect8'
[878]	607	do k=nlay,1,-1
	608	do ig=1,ngrid
	609	if (zqla(ig,k).gt.1e-10) then
	610	nivcon(ig)=k
	611	zcon(ig)=zlev(ig,k)
	612	endif
	613	enddo
	614	enddo
[972]	615	if (prt_level.ge.1) print*,'14c OK convect8'
[878]	616	!calcul des moments
	617	!initialisation
	618	do l=1,nlay
	619	do ig=1,ngrid
	620	q2(ig,l)=0.
	621	wth2(ig,l)=0.
	622	wth3(ig,l)=0.
	623	ratqscth(ig,l)=0.
	624	ratqsdiff(ig,l)=0.
	625	enddo
	626	enddo
[972]	627	if (prt_level.ge.1) print*,'14d OK convect8'
[1134]	628	if (prt_level.ge.10)write(lunout,*) &
	629	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
[878]	630	do l=1,nlay
	631	do ig=1,ngrid
	632	zf=fraca(ig,l)
	633	zf2=zf/(1.-zf)
[972]	634	!
	635	if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
	636	!
	637	if (prt_level.ge.10) print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
[878]	638	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))*2
[972]	639	if(zw2(ig,l).gt.1.e-10) then
	640	wth2(ig,l)=zf2(zw2(ig,l))*2
	641	else
	642	wth2(ig,l)=0.
	643	endif
[878]	644	! print*,'wth2=',wth2(ig,l)
	645	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
	646	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
[972]	647	if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
[878]	648	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
	649	!test: on calcul q2/po=ratqsc
	650	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
	651	enddo
	652	enddo
[879]	653	!calcul de ale_bl et alp_bl
	654	!pour le calcul d'une valeur intégrée entre la surface et lmax
	655	do ig=1,ngrid
	656	alp_int(ig)=0.
	657	ale_int(ig)=0.
	658	n_int(ig)=0
	659	enddo
[972]	660	!
	661	do l=1,nlay
[879]	662	do ig=1,ngrid
[972]	663	if(l.LE.lmax(ig)) THEN
	664	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)
	665	ale_int(ig)=ale_int(ig)+0.5zw2(ig,l)*2
	666	n_int(ig)=n_int(ig)+1
	667	endif
[879]	668	enddo
	669	enddo
	670	! print*,'avant calcul ale et alp'
	671	!calcul de ALE et ALP pour la convection
	672	do ig=1,ngrid
	673	! Alp_bl(ig)=0.5rhobarz(ig,lmix_bis(ig))wth3(ig,lmix(ig))
	674	! Alp_bl(ig)=0.5rhobarz(ig,nivcon(ig))wth3(ig,nivcon(ig))
	675	! Alp_bl(ig)=0.5rhobarz(ig,lmix(ig))wth3(ig,lmix(ig))
	676	! & *0.1
	677	!valeur integree de alp_bl * 0.5:
	678	if (n_int(ig).gt.0) then
	679	Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
	680	! if (Alp_bl(ig).lt.0.) then
	681	! Alp_bl(ig)=0.
	682	endif
	683	! endif
	684	! write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
	685	! s wth3(ig,nivcon(ig)),Alp_bl(ig)
	686	! write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
	687	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
	688	! if (nivcon(ig).eq.1) then
	689	! Ale_bl(ig)=0.
	690	! else
	691	!valeur max de ale_bl:
	692	Ale_bl(ig)=0.5zw2(ig,lmix(ig))*2
	693	! & /2.
	694	! & *0.1
	695	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
	696	! if (n_int(ig).gt.0) then
	697	! Ale_bl(ig)=ale_int(ig)/n_int(ig)
	698	! Ale_bl(ig)=4.
	699	! endif
	700	! endif
	701	! Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
	702	! Ale_bl(ig)=wth2(ig,nivcon(ig))
	703	! write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
	704	enddo
	705	!test:calcul de la ponderation des couches pour KE
	706	!initialisations
	707	! print*,'ponderation'
	708	do ig=1,ngrid
	709	fm_tot(ig)=0.
	710	enddo
	711	do ig=1,ngrid
	712	do k=1,klev
	713	wght_th(ig,k)=1.
	714	enddo
	715	enddo
	716	do ig=1,ngrid
	717	! lalim_conv(ig)=lmix_bis(ig)
	718	!la hauteur de la couche alim_conv = hauteur couche alim_therm
	719	lalim_conv(ig)=lalim(ig)
	720	! zentr(ig)=zlev(ig,lalim(ig))
	721	enddo
	722	do ig=1,ngrid
	723	do k=1,lalim_conv(ig)
	724	fm_tot(ig)=fm_tot(ig)+fm(ig,k)
	725	enddo
	726	enddo
	727	do ig=1,ngrid
	728	do k=1,lalim_conv(ig)
	729	if (fm_tot(ig).gt.1.e-10) then
	730	! wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
	731	endif
	732	!on pondere chaque couche par a*
	733	if (alim_star(ig,k).gt.1.e-10) then
	734	wght_th(ig,k)=alim_star(ig,k)
	735	else
	736	wght_th(ig,k)=1.
	737	endif
	738	enddo
	739	enddo
	740	! print*,'apres wght_th'
	741	!test pour prolonger la convection
	742	do ig=1,ngrid
[926]	743	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
	744	if ((alim_star(ig,1).lt.1.e-10)) then
[879]	745	lalim_conv(ig)=1
	746	wght_th(ig,1)=1.
	747	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
	748	endif
	749	enddo
	750
[878]	751	!calcul du ratqscdiff
[972]	752	if (prt_level.ge.1) print*,'14e OK convect8'
[878]	753	var=0.
	754	vardiff=0.
	755	ratqsdiff(:,:)=0.
	756	do ig=1,ngrid
	757	do l=1,lalim(ig)
	758	var=var+alim_star(ig,l)zqta(ig,l)1000.
	759	enddo
	760	enddo
[972]	761	if (prt_level.ge.1) print*,'14f OK convect8'
[878]	762	do ig=1,ngrid
	763	do l=1,lalim(ig)
	764	zf=fraca(ig,l)
	765	zf2=zf/(1.-zf)
	766	vardiff=vardiff+alim_star(ig,l) &
	767	& (zqta(ig,l)1000.-var)**2
	768	! ratqsdiff=ratqsdiff+alim_star(ig,l)*
	769	! s (zqta(ig,l)1000.-po(ig,l)1000.)**2
	770	enddo
	771	enddo
[972]	772	if (prt_level.ge.1) print*,'14g OK convect8'
[878]	773	do l=1,nlay
	774	do ig=1,ngrid
	775	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
	776	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
	777	enddo
	778	enddo
	779	!--------------------------------------------------------------------
	780	!
	781	!ecriture des fichiers sortie
	782	! print*,'15 OK convect8'
	783
[938]	784	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
[1026]	785	#ifdef wrgrads_thermcell
[878]	786	#include "thermcell_out3d.h"
	787	#endif
	788
	789	endif
	790
[938]	791	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
[878]	792
[972]	793	! if(icount.eq.501) stop'au pas 301 dans thermcell_main'
[878]	794	return
	795	end
	796
	797	!-----------------------------------------------------------------------------
	798
	799	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
[938]	800	IMPLICIT NONE
	801	#include "iniprint.h"
[878]	802
[938]	803	integer i, k, klon,klev
[878]	804	real pplev(klon,klev+1),pplay(klon,klev)
	805	real ztv(klon,klev)
	806	real po(klon,klev)
	807	real ztva(klon,klev)
	808	real zqla(klon,klev)
	809	real f_star(klon,klev)
	810	real zw2(klon,klev)
	811	integer long(klon)
	812	real seuil
	813	character*21 comment
	814
[938]	815	if (prt_level.ge.1) THEN
	816	print*,'WARNING !!! TEST ',comment
	817	endif
[879]	818	return
	819
[878]	820	! test sur la hauteur des thermiques ...
	821	do i=1,klon
[972]	822	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
	823	if (prt_level.ge.10) then
[878]	824	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
	825	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
	826	do k=1,klev
	827	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)
	828	enddo
	829	! stop
[972]	830	endif
[878]	831	enddo
	832
	833
	834	return
	835	end
	836

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

Download in other formats:

Original Format