Context Navigation

thermcell_main.F90 @ 1085

Last change on this file since 1085 was 1026, checked in by lmdzadmin, 16 years ago
Modifs thermiques FH/IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 26.7 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
	234	print*,'WARNING thermcell_main f0=max(f0,1.e-2)'
	235	do ig=1,klon
	236	if (prt_level.ge.20) then
	237	print*,'th_main ig f0',ig,f0(ig)
[878]	238	endif
[972]	239	f0(ig)=max(f0(ig),1.e-2)
	240	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
	241	enddo
[878]	242
	243	!-----------------------------------------------------------------------
	244	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
	245	! --------------------------------------------------------------------
	246	!
	247	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
	248	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
	249
[938]	250	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
[878]	251
	252	!------------------------------------------------------------------------
	253	! --------------------
	254	!
	255	!
	256	! + + + + + + + + + + +
	257	!
	258	!
	259	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
	260	! wh,wt,wo ...
	261	!
	262	! + + + + + + + + + + + zh,zu,zv,zo,rho
	263	!
	264	!
	265	! -------------------- zlev(1)
	266	! \\\\\\\\\\\\\\\\\\\\
	267	!
	268	!
	269
	270	!-----------------------------------------------------------------------
	271	! Calcul des altitudes des couches
	272	!-----------------------------------------------------------------------
	273
	274	do l=2,nlay
	275	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
	276	enddo
	277	zlev(:,1)=0.
	278	zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
	279	do l=1,nlay
	280	zlay(:,l)=pphi(:,l)/RG
	281	enddo
	282	!calcul de l epaisseur des couches
	283	do l=1,nlay
	284	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
	285	enddo
	286
	287	! print*,'2 OK convect8'
	288	!-----------------------------------------------------------------------
	289	! Calcul des densites
	290	!-----------------------------------------------------------------------
	291
	292	do l=1,nlay
	293	rho(:,l)=pplay(:,l)/(zpspsk(:,l)RDztv(:,l))
	294	enddo
	295
[972]	296	!IM
	297	print*,'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
	298	rhobarz(:,1)=rho(:,1)
	299
[878]	300	do l=2,nlay
	301	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
	302	enddo
	303
	304	!calcul de la masse
	305	do l=1,nlay
	306	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
	307	enddo
	308
[938]	309	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
[878]	310
	311	!------------------------------------------------------------------
	312	!
	313	! /\|\
	314	! -------- \| F_k+1 -------
	315	! ----> D_k
	316	! /\|\ <---- E_k , A_k
	317	! -------- \| F_k ---------
	318	! ----> D_k-1
	319	! <---- E_k-1 , A_k-1
	320	!
	321	!
	322	!
	323	!
	324	!
	325	! ---------------------------
	326	!
	327	! ----- F_lmax+1=0 ---------- \
	328	! lmax (zmax) \|
	329	! --------------------------- \|
	330	! \|
	331	! --------------------------- \|
	332	! \|
	333	! --------------------------- \|
	334	! \|
	335	! --------------------------- \|
	336	! \|
	337	! --------------------------- \|
	338	! \| E
	339	! --------------------------- \| D
	340	! \|
	341	! --------------------------- \|
	342	! \|
	343	! --------------------------- \ \|
	344	! lalim \| \|
	345	! --------------------------- \| \|
	346	! \| \|
	347	! --------------------------- \| \|
	348	! \| A \|
	349	! --------------------------- \| \|
	350	! \| \|
	351	! --------------------------- \| \|
	352	! lmin (=1 pour le moment) \| \|
	353	! ----- F_lmin=0 ------------ / /
	354	!
	355	! ---------------------------
	356	! //////////////////////////
	357	!
	358	!
	359	!=============================================================================
	360	! Calculs initiaux ne faisant pas intervenir les changements de phase
	361	!=============================================================================
	362
	363	!------------------------------------------------------------------
	364	! 1. alim_star est le profil vertical de l'alimentation à la base du
	365	! panache thermique, calculé à partir de la flotabilité de l'air sec
	366	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
	367	!------------------------------------------------------------------
	368	!
	369	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
[972]	370	CALL thermcell_init(ngrid,nlay,ztv,zlay,zlev, &
[878]	371	& lalim,lmin,alim_star,alim_star_tot,lev_out)
	372
	373	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin ')
	374	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
	375
	376
[938]	377	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
	378	if (prt_level.ge.10) then
[972]	379	write(lunout1,*) 'Dans thermcell_main 1'
	380	write(lunout1,*) 'lmin ',lmin(igout)
	381	write(lunout1,*) 'lalim ',lalim(igout)
	382	write(lunout1,*) ' ig l alim_star thetav'
	383	write(lunout1,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
[878]	384	& ,ztv(igout,l),l=1,lalim(igout)+4)
	385	endif
	386
[938]	387	!v1d do ig=1,klon
[926]	388	!v1d if (alim_star(ig,1).gt.1.e-10) then
	389	!v1d therm=.true.
	390	!v1d endif
[938]	391	!v1d enddo
[878]	392	!-----------------------------------------------------------------------------
	393	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
	394	! panache sec conservatif (e=d=0) alimente selon alim_star
	395	! Il s'agit d'un calcul de type CAPE
	396	! zmax_sec est utilisé pour déterminer la géométrie du thermique.
	397	!------------------------------------------------------------------------------
	398	!
	399	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
	400	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
	401
	402	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
	403	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
	404
[938]	405	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
	406	if (prt_level.ge.10) then
[972]	407	write(lunout1,*) 'Dans thermcell_main 1b'
	408	write(lunout1,*) 'lmin ',lmin(igout)
	409	write(lunout1,*) 'lalim ',lalim(igout)
	410	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
	411	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
[878]	412	& ,l=1,lalim(igout)+4)
	413	endif
	414
	415
	416
	417	!---------------------------------------------------------------------------------
	418	!calcul du melange et des variables dans le thermique
	419	!--------------------------------------------------------------------------------
	420	!
[972]	421	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
	422	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
	423	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
[1026]	424	& zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot, &
[878]	425	& lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, &
[1026]	426	& ztla,zqla,zqta,zha,zw2,zw_est,zqsatth,lmix,lmix_bis,linter &
	427	& ,lev_out,lunout1,igout)
[972]	428	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
	429
[878]	430	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
	431	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
	432
[938]	433	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
	434	if (prt_level.ge.10) then
[972]	435	write(lunout1,*) 'Dans thermcell_main 2'
	436	write(lunout1,*) 'lmin ',lmin(igout)
	437	write(lunout1,*) 'lalim ',lalim(igout)
	438	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
	439	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
[878]	440	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
	441	endif
	442
	443	!-------------------------------------------------------------------------------
	444	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
	445	!-------------------------------------------------------------------------------
	446	!
	447	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
	448	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
	449
	450
	451	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
	452	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
	453	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
	454	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
	455
[938]	456	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
[878]	457
	458	!-------------------------------------------------------------------------------
	459	! Fermeture,determination de f
	460	!-------------------------------------------------------------------------------
[1026]	461	!
	462	!avant closure: on redéfinit lalim, alim_star_tot et alim_star
	463	! do ig=1,klon
	464	! do l=2,lalim(ig)
	465	! alim_star(ig,l)=entr_star(ig,l)
	466	! entr_star(ig,l)=0.
	467	! enddo
	468	! enddo
[878]	469
	470	CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
[1026]	471	& zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
[878]	472
[938]	473	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
[878]	474
[972]	475	if (tau_thermals>1.) then
	476	lambda=exp(-ptimestep/tau_thermals)
	477	f0=(1.-lambda)f+lambdaf0
	478	else
	479	f0=f
	480	endif
	481
	482	! Test valable seulement en 1D mais pas genant
	483	if (.not. (f0(1).ge.0.) ) then
	484	stop'Dans thermcell_main'
	485	endif
	486
[878]	487	!-------------------------------------------------------------------------------
	488	!deduction des flux
	489	!-------------------------------------------------------------------------------
	490
[972]	491	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
[878]	492	& lalim,lmax,alim_star, &
	493	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
[972]	494	& detr,zqla,lev_out,lunout1,igout)
	495	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
[878]	496
[938]	497	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
[878]	498	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
	499	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
	500
	501	!------------------------------------------------------------------
[972]	502	! On ne prend pas directement les profils issus des calculs precedents
	503	! mais on s'autorise genereusement une relaxation vers ceci avec
	504	! une constante de temps tau_thermals (typiquement 1800s).
	505	!------------------------------------------------------------------
[878]	506
[972]	507	if (tau_thermals>1.) then
	508	lambda=exp(-ptimestep/tau_thermals)
	509	fm0=(1.-lambda)fm+lambdafm0
	510	entr0=(1.-lambda)entr+lambdaentr0
	511	! detr0=(1.-lambda)detr+lambdadetr0
[878]	512	else
	513	fm0=fm
	514	entr0=entr
	515	detr0=detr
	516	endif
	517
[972]	518	!c------------------------------------------------------------------
	519	! calcul du transport vertical
	520	!------------------------------------------------------------------
	521
[878]	522	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
	523	& zthl,zdthladj,zta,lev_out)
	524	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
	525	& po,pdoadj,zoa,lev_out)
	526
[883]	527	!------------------------------------------------------------------
	528	! Calcul de la fraction de l'ascendance
	529	!------------------------------------------------------------------
	530	do ig=1,klon
	531	fraca(ig,1)=0.
	532	fraca(ig,nlay+1)=0.
	533	enddo
	534	do l=2,nlay
	535	do ig=1,klon
	536	if (zw2(ig,l).gt.1.e-10) then
	537	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
	538	else
	539	fraca(ig,l)=0.
	540	endif
	541	enddo
	542	enddo
	543
	544	!------------------------------------------------------------------
	545	! calcul du transport vertical du moment horizontal
	546	!------------------------------------------------------------------
[878]	547
[972]	548	!IM 090508
[883]	549	if (1.eq.1) then
[972]	550	!IM 070508 vers. _dq
	551	! if (1.eq.0) then
[883]	552
	553
[878]	554	! Calcul du transport de V tenant compte d'echange par gradient
	555	! de pression horizontal avec l'environnement
	556
	557	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
	558	& ,fraca,zmax &
[972]	559	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
	560	!IM 050508 & ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
[878]	561	else
	562
	563	! calcul purement conservatif pour le transport de V
	564	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
	565	& ,zu,pduadj,zua,lev_out)
	566	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
	567	& ,zv,pdvadj,zva,lev_out)
	568	endif
	569
	570	! print*,'13 OK convect8'
	571	do l=1,nlay
	572	do ig=1,ngrid
	573	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
	574	enddo
	575	enddo
	576
[972]	577	if (prt_level.ge.1) print*,'14 OK convect8'
[878]	578	!------------------------------------------------------------------
	579	! Calculs de diagnostiques pour les sorties
	580	!------------------------------------------------------------------
	581	!calcul de fraca pour les sorties
	582
	583	if (sorties) then
[972]	584	if (prt_level.ge.1) print*,'14a OK convect8'
[878]	585	! calcul du niveau de condensation
	586	! initialisation
	587	do ig=1,ngrid
[879]	588	nivcon(ig)=0
[878]	589	zcon(ig)=0.
	590	enddo
	591	!nouveau calcul
	592	do ig=1,ngrid
	593	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
	594	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
	595	enddo
	596	do k=1,nlay
	597	do ig=1,ngrid
	598	if ((pcon(ig).le.pplay(ig,k)) &
	599	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
	600	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
	601	endif
	602	enddo
	603	enddo
[972]	604	if (prt_level.ge.1) print*,'14b OK convect8'
[878]	605	do k=nlay,1,-1
	606	do ig=1,ngrid
	607	if (zqla(ig,k).gt.1e-10) then
	608	nivcon(ig)=k
	609	zcon(ig)=zlev(ig,k)
	610	endif
	611	enddo
	612	enddo
[972]	613	if (prt_level.ge.1) print*,'14c OK convect8'
[878]	614	!calcul des moments
	615	!initialisation
	616	do l=1,nlay
	617	do ig=1,ngrid
	618	q2(ig,l)=0.
	619	wth2(ig,l)=0.
	620	wth3(ig,l)=0.
	621	ratqscth(ig,l)=0.
	622	ratqsdiff(ig,l)=0.
	623	enddo
	624	enddo
[972]	625	if (prt_level.ge.1) print*,'14d OK convect8'
	626	print*,'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
[878]	627	do l=1,nlay
	628	do ig=1,ngrid
	629	zf=fraca(ig,l)
	630	zf2=zf/(1.-zf)
[972]	631	!
	632	if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
	633	!
	634	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]	635	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))*2
[972]	636	if(zw2(ig,l).gt.1.e-10) then
	637	wth2(ig,l)=zf2(zw2(ig,l))*2
	638	else
	639	wth2(ig,l)=0.
	640	endif
[878]	641	! print*,'wth2=',wth2(ig,l)
	642	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
	643	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
[972]	644	if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
[878]	645	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
	646	!test: on calcul q2/po=ratqsc
	647	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
	648	enddo
	649	enddo
[879]	650	!calcul de ale_bl et alp_bl
	651	!pour le calcul d'une valeur intégrée entre la surface et lmax
	652	do ig=1,ngrid
	653	alp_int(ig)=0.
	654	ale_int(ig)=0.
	655	n_int(ig)=0
	656	enddo
[972]	657	!
	658	do l=1,nlay
[879]	659	do ig=1,ngrid
[972]	660	if(l.LE.lmax(ig)) THEN
	661	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)
	662	ale_int(ig)=ale_int(ig)+0.5zw2(ig,l)*2
	663	n_int(ig)=n_int(ig)+1
	664	endif
[879]	665	enddo
	666	enddo
	667	! print*,'avant calcul ale et alp'
	668	!calcul de ALE et ALP pour la convection
	669	do ig=1,ngrid
	670	! Alp_bl(ig)=0.5rhobarz(ig,lmix_bis(ig))wth3(ig,lmix(ig))
	671	! Alp_bl(ig)=0.5rhobarz(ig,nivcon(ig))wth3(ig,nivcon(ig))
	672	! Alp_bl(ig)=0.5rhobarz(ig,lmix(ig))wth3(ig,lmix(ig))
	673	! & *0.1
	674	!valeur integree de alp_bl * 0.5:
	675	if (n_int(ig).gt.0) then
	676	Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
	677	! if (Alp_bl(ig).lt.0.) then
	678	! Alp_bl(ig)=0.
	679	endif
	680	! endif
	681	! write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
	682	! s wth3(ig,nivcon(ig)),Alp_bl(ig)
	683	! write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
	684	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
	685	! if (nivcon(ig).eq.1) then
	686	! Ale_bl(ig)=0.
	687	! else
	688	!valeur max de ale_bl:
	689	Ale_bl(ig)=0.5zw2(ig,lmix(ig))*2
	690	! & /2.
	691	! & *0.1
	692	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
	693	! if (n_int(ig).gt.0) then
	694	! Ale_bl(ig)=ale_int(ig)/n_int(ig)
	695	! Ale_bl(ig)=4.
	696	! endif
	697	! endif
	698	! Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
	699	! Ale_bl(ig)=wth2(ig,nivcon(ig))
	700	! write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
	701	enddo
	702	!test:calcul de la ponderation des couches pour KE
	703	!initialisations
	704	! print*,'ponderation'
	705	do ig=1,ngrid
	706	fm_tot(ig)=0.
	707	enddo
	708	do ig=1,ngrid
	709	do k=1,klev
	710	wght_th(ig,k)=1.
	711	enddo
	712	enddo
	713	do ig=1,ngrid
	714	! lalim_conv(ig)=lmix_bis(ig)
	715	!la hauteur de la couche alim_conv = hauteur couche alim_therm
	716	lalim_conv(ig)=lalim(ig)
	717	! zentr(ig)=zlev(ig,lalim(ig))
	718	enddo
	719	do ig=1,ngrid
	720	do k=1,lalim_conv(ig)
	721	fm_tot(ig)=fm_tot(ig)+fm(ig,k)
	722	enddo
	723	enddo
	724	do ig=1,ngrid
	725	do k=1,lalim_conv(ig)
	726	if (fm_tot(ig).gt.1.e-10) then
	727	! wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
	728	endif
	729	!on pondere chaque couche par a*
	730	if (alim_star(ig,k).gt.1.e-10) then
	731	wght_th(ig,k)=alim_star(ig,k)
	732	else
	733	wght_th(ig,k)=1.
	734	endif
	735	enddo
	736	enddo
	737	! print*,'apres wght_th'
	738	!test pour prolonger la convection
	739	do ig=1,ngrid
[926]	740	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
	741	if ((alim_star(ig,1).lt.1.e-10)) then
[879]	742	lalim_conv(ig)=1
	743	wght_th(ig,1)=1.
	744	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
	745	endif
	746	enddo
	747
[878]	748	!calcul du ratqscdiff
[972]	749	if (prt_level.ge.1) print*,'14e OK convect8'
[878]	750	var=0.
	751	vardiff=0.
	752	ratqsdiff(:,:)=0.
	753	do ig=1,ngrid
	754	do l=1,lalim(ig)
	755	var=var+alim_star(ig,l)zqta(ig,l)1000.
	756	enddo
	757	enddo
[972]	758	if (prt_level.ge.1) print*,'14f OK convect8'
[878]	759	do ig=1,ngrid
	760	do l=1,lalim(ig)
	761	zf=fraca(ig,l)
	762	zf2=zf/(1.-zf)
	763	vardiff=vardiff+alim_star(ig,l) &
	764	& (zqta(ig,l)1000.-var)**2
	765	! ratqsdiff=ratqsdiff+alim_star(ig,l)*
	766	! s (zqta(ig,l)1000.-po(ig,l)1000.)**2
	767	enddo
	768	enddo
[972]	769	if (prt_level.ge.1) print*,'14g OK convect8'
[878]	770	do l=1,nlay
	771	do ig=1,ngrid
	772	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
	773	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
	774	enddo
	775	enddo
	776	!--------------------------------------------------------------------
	777	!
	778	!ecriture des fichiers sortie
	779	! print*,'15 OK convect8'
	780
[938]	781	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
[1026]	782	#ifdef wrgrads_thermcell
[878]	783	#include "thermcell_out3d.h"
	784	#endif
	785
	786	endif
	787
[938]	788	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
[878]	789
[972]	790	! if(icount.eq.501) stop'au pas 301 dans thermcell_main'
[878]	791	return
	792	end
	793
	794	!-----------------------------------------------------------------------------
	795
	796	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
[938]	797	IMPLICIT NONE
	798	#include "iniprint.h"
[878]	799
[938]	800	integer i, k, klon,klev
[878]	801	real pplev(klon,klev+1),pplay(klon,klev)
	802	real ztv(klon,klev)
	803	real po(klon,klev)
	804	real ztva(klon,klev)
	805	real zqla(klon,klev)
	806	real f_star(klon,klev)
	807	real zw2(klon,klev)
	808	integer long(klon)
	809	real seuil
	810	character*21 comment
	811
[938]	812	if (prt_level.ge.1) THEN
	813	print*,'WARNING !!! TEST ',comment
	814	endif
[879]	815	return
	816
[878]	817	! test sur la hauteur des thermiques ...
	818	do i=1,klon
[972]	819	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
	820	if (prt_level.ge.10) then
[878]	821	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
	822	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
	823	do k=1,klev
	824	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)
	825	enddo
	826	! stop
[972]	827	endif
[878]	828	enddo
	829
	830
	831	return
	832	end
	833

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

Download in other formats:

Original Format