Context Navigation

thermcell_main.F90 @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 27.0 KB

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