Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

concvl.F

Last change on this file was 1398, checked in by musat, 14 years ago

Last corrections for CMIP5:

Add O3 at standard level files histmthNMC.nc
Add positive attribute "down" for vertical axes for all output files
Replace "inst" by "ave" for hist*NMC.nc files to have time_counter and bounds for time axis (Marie-Alice's hint)
Correct units for vertical axes : mb instead of hPa
Add mass flux at the bottom of clouds
Comment non initialized variables (s_capCL, s_oliqCL, s_cteiCL, s_trmb1, s_trmb2, s_trmb3) for the output files
Geopotential field phy850, phi700, phi500, etc are modified to "geopotential height and are called z850, z700, z500, etc
Meaning of specific humidity outputs - ovapinit and ovap - were interchanged
Fields albs, albslw become alb1, alb2 in output files
Correct title for rugs_* fields
Correct units for pbase and ptop are Pa (not mb)
Correct ndayrain field

FH/JLD/JYG/MAF/IM

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

File size: 15.7 KB

Rev	Line
[879]	1	SUBROUTINE concvl (iflag_con,iflag_clos,
	2	. dtime,paprs,pplay,
[1146]	3	. t,q,t_wake,q_wake,s_wake,u,v,tra,ntra,
[879]	4	. ALE,ALP,work1,work2,
	5	. d_t,d_q,d_u,d_v,d_tra,
	6	. rain, snow, kbas, ktop, sigd,
[1398]	7	. cbmf,upwd,dnwd,dnwdbis,Ma,mip,Vprecip,
[879]	8	. cape,cin,tvp,Tconv,iflag,
[524]	9	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
[879]	10	. qcondc,wd,pmflxr,pmflxs,
	11	. da,phi,mp,dd_t,dd_q,lalim_conv,wght_th)
	12	***************************************************************
	13	* *
	14	* CONCVL *
	15	* *
	16	* *
	17	* written by : Sandrine Bony-Lena, 17/05/2003, 11.16.04 *
	18	* modified by : *
	19	***************************************************************
	20	*
[524]	21	c
[940]	22	USE dimphy
[1146]	23	USE infotrac, ONLY : nbtr
[524]	24	IMPLICIT none
	25	c======================================================================
[879]	26	c Auteur(s): S. Bony-Lena (LMD/CNRS) date: ???
[524]	27	c Objet: schema de convection de Emanuel (1991) interface
	28	c======================================================================
	29	c Arguments:
	30	c dtime--input-R-pas d'integration (s)
	31	c s-------input-R-la valeur "s" pour chaque couche
	32	c sigs----input-R-la valeur "sigma" de chaque couche
	33	c sig-----input-R-la valeur de "sigma" pour chaque niveau
	34	c psolpa--input-R-la pression au sol (en Pa)
	35	C pskapa--input-R-exponentiel kappa de psolpa
	36	c h-------input-R-enthalpie potentielle (CpT/P*kappa)
	37	c q-------input-R-vapeur d'eau (en kg/kg)
	38	c
	39	c work*: input et output: deux variables de travail,
	40	c on peut les mettre a 0 au debut
	41	c ALE-----input-R-energie disponible pour soulevement
[879]	42	c ALP-----input-R-puissance disponible pour soulevement
[524]	43	c
	44	C d_h-----output-R-increment de l'enthalpie potentielle (h)
	45	c d_q-----output-R-increment de la vapeur d'eau
	46	c rain----output-R-la pluie (mm/s)
	47	c snow----output-R-la neige (mm/s)
	48	c upwd----output-R-saturated updraft mass flux (kg/m**2/s)
	49	c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
	50	c dnwd0---output-R-unsaturated downdraft mass flux (kg/m**2/s)
[879]	51	c Ma------output-R-adiabatic ascent mass flux (kg/m2/s)
	52	c mip-----output-R-mass flux shed by adiabatic ascent (kg/m2/s)
	53	c Vprecip-output-R-vertical profile of precipitations (kg/m2/s)
	54	c Tconv---output-R-environment temperature seen by convective scheme (K)
[524]	55	c Cape----output-R-CAPE (J/kg)
[879]	56	c Cin ----output-R-CIN (J/kg)
[524]	57	c Tvp-----output-R-Temperature virtuelle d'une parcelle soulevee
	58	c adiabatiquement a partir du niveau 1 (K)
	59	c deltapb-output-R-distance entre LCL et base de la colonne (<0 ; Pa)
	60	c Ice_flag-input-L-TRUE->prise en compte de la thermodynamique de la glace
[879]	61	c dd_t-----output-R-increment de la temperature du aux descentes precipitantes
	62	c dd_q-----output-R-increment de la vapeur d'eau du aux desc precip
[524]	63	c======================================================================
	64	c
	65	#include "dimensions.h"
	66	c
[879]	67	INTEGER iflag_con,iflag_clos
[524]	68	c
	69	REAL dtime, paprs(klon,klev+1),pplay(klon,klev)
	70	REAL t(klon,klev),q(klon,klev),u(klon,klev),v(klon,klev)
[879]	71	REAL t_wake(klon,klev),q_wake(klon,klev)
[1146]	72	Real s_wake(klon)
	73	REAL tra(klon,klev,nbtr)
[524]	74	INTEGER ntra
[879]	75	REAL work1(klon,klev),work2(klon,klev),ptop2(klon)
[619]	76	REAL pmflxr(klon,klev+1),pmflxs(klon,klev+1)
[879]	77	REAL ALE(klon),ALP(klon)
[524]	78	c
	79	REAL d_t(klon,klev),d_q(klon,klev),d_u(klon,klev),d_v(klon,klev)
[879]	80	REAL dd_t(klon,klev),dd_q(klon,klev)
[1146]	81	REAL d_tra(klon,klev,nbtr)
[524]	82	REAL rain(klon),snow(klon)
	83	c
	84	INTEGER kbas(klon),ktop(klon)
	85	REAL em_ph(klon,klev+1),em_p(klon,klev)
	86	REAL upwd(klon,klev),dnwd(klon,klev),dnwdbis(klon,klev)
[1334]	87
	88	!! REAL Ma(klon,klev), mip(klon,klev),Vprecip(klon,klev) !jyg
	89	REAL Ma(klon,klev), mip(klon,klev),Vprecip(klon,klev+1) !jyg
	90
[619]	91	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
[879]	92	REAL cape(klon),cin(klon),tvp(klon,klev)
	93	REAL Tconv(klon,klev)
	94	c
	95	cCR:test: on passe lentr et alim_star des thermiques
	96	INTEGER lalim_conv(klon)
	97	REAL wght_th(klon,klev)
	98	REAL em_sig1feed ! sigma at lower bound of feeding layer
	99	REAL em_sig2feed ! sigma at upper bound of feeding layer
	100	REAL em_wght(klev) ! weight density determining the feeding mixture
	101	con enleve le save
	102	c SAVE em_sig1feed,em_sig2feed,em_wght
	103	c
[524]	104	INTEGER iflag(klon)
	105	REAL rflag(klon)
	106	REAL pbase(klon),bbase(klon)
	107	REAL dtvpdt1(klon,klev),dtvpdq1(klon,klev)
	108	REAL dplcldt(klon),dplcldr(klon)
	109	REAL qcondc(klon,klev)
	110	REAL wd(klon)
[879]	111	REAL Plim1(klon),Plim2(klon),asupmax(klon,klev)
	112	REAL supmax0(klon),asupmaxmin(klon)
[524]	113	c
[879]	114	REAL sigd(klon)
[524]	115	REAL zx_t,zdelta,zx_qs,zcor
	116	c
[963]	117	! INTEGER iflag_mix
	118	! SAVE iflag_mix
[524]	119	INTEGER noff, minorig
	120	INTEGER i,k,itra
[879]	121	REAL qs(klon,klev),qs_wake(klon,klev)
[1398]	122	REAL cbmf(klon)
[940]	123	cLF SAVE cbmf
[1398]	124	cIM/JYG REAL, SAVE, ALLOCATABLE :: cbmf(:)
	125	ccc$OMP THREADPRIVATE(cbmf)!
[940]	126	REAL cbmflast(klon)
[524]	127	INTEGER ifrst
	128	SAVE ifrst
	129	DATA ifrst /0/
[766]	130	c$OMP THREADPRIVATE(ifrst)
	131
[879]	132	c
	133	C Variables supplementaires liees au bilan d'energie
	134	c Real paire(klon)
[940]	135	cLF Real ql(klon,klev)
[879]	136	c Save paire
[940]	137	cLF Save ql
	138	cLF Real t1(klon,klev),q1(klon,klev)
	139	cLF Save t1,q1
[879]	140	c Data paire /1./
[940]	141	REAL, SAVE, ALLOCATABLE :: ql(:,:), q1(:,:), t1(:,:)
	142	c$OMP THREADPRIVATE(ql, q1, t1)
[879]	143	c
	144	C Variables liees au bilan d'energie et d'enthalpi
	145	REAL ztsol(klon)
	146	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	147	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	148	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	149	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
[987]	150	c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot)
	151	c$OMP THREADPRIVATE(h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
[879]	152	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
	153	REAL d_h_vcol_phy
	154	REAL fs_bound, fq_bound
	155	SAVE d_h_vcol_phy
[987]	156	c$OMP THREADPRIVATE(d_h_vcol_phy)
[879]	157	REAL zero_v(klon)
	158	CHARACTER*15 ztit
	159	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
	160	SAVE ip_ebil
	161	DATA ip_ebil/2/
[987]	162	c$OMP THREADPRIVATE(ip_ebil)
[879]	163	INTEGER if_ebil ! level for energy conserv. dignostics
	164	SAVE if_ebil
	165	DATA if_ebil/2/
[987]	166	c$OMP THREADPRIVATE(if_ebil)
[879]	167	c+jld ec_conser
	168	REAL d_t_ec(klon,klev) ! tendance du a la conersion Ec -> E thermique
	169	REAL ZRCPD
	170	c-jld ec_conser
[940]	171	cLF
[973]	172	INTEGER nloc
[940]	173	logical, save :: first=.true.
[987]	174	c$OMP THREADPRIVATE(first)
	175	INTEGER, SAVE :: itap, igout
	176	c$OMP THREADPRIVATE(itap, igout)
[879]	177	c
[524]	178	#include "YOMCST.h"
[879]	179	#include "YOMCST2.h"
[524]	180	#include "YOETHF.h"
	181	#include "FCTTRE.h"
[973]	182	#include "iniprint.h"
[524]	183	c
[940]	184	if (first) then
	185	c Allocate some variables LF 04/2008
	186	c
[1398]	187	cIM/JYG allocate(cbmf(klon))
[940]	188	allocate(ql(klon,klev))
	189	allocate(t1(klon,klev))
	190	allocate(q1(klon,klev))
[973]	191	itap=0
	192	igout=klon/2+1/klon
[940]	193	endif
[973]	194	c Incrementer le compteur de la physique
	195	itap = itap + 1
[879]	196
	197	c Copy T into Tconv
	198	DO k = 1,klev
	199	DO i = 1,klon
	200	Tconv(i,k) = T(i,k)
	201	ENDDO
	202	ENDDO
[524]	203	c
[879]	204	IF (if_ebil.ge.1) THEN
	205	DO i=1,klon
	206	ztsol(i) = t(i,1)
	207	zero_v(i)=0.
	208	Do k = 1,klev
	209	ql(i,k) = 0.
	210	ENDDO
	211	END DO
	212	END IF
	213	c
[559]	214	cym
	215	snow(:)=0
	216
[940]	217	c IF (ifrst .EQ. 0) THEN
	218	c ifrst = 1
	219	if (first) then
	220	first=.false.
[879]	221	c
	222	C===========================================================================
	223	C READ IN PARAMETERS FOR THE CLOSURE AND THE MIXING DISTRIBUTION
	224	C===========================================================================
	225	C
	226	if (iflag_con.eq.3) then
[963]	227	c CALL cv3_inicp()
	228	CALL cv3_inip()
[879]	229	endif
	230	c
	231	C===========================================================================
	232	C READ IN PARAMETERS FOR CONVECTIVE INHIBITION BY TROPOS. DRYNESS
	233	C===========================================================================
	234	C
[987]	235	cc$$$ open (56,file='supcrit.data')
	236	cc$$$ read (56,*) Supcrit1, Supcrit2
	237	cc$$$ close (56)
[879]	238	c
	239	print*, 'supcrit1, supcrit2' ,supcrit1, supcrit2
	240	C
	241	C===========================================================================
	242	C Initialisation pour les bilans d'eau et d'energie
	243	C===========================================================================
	244	IF (if_ebil.ge.1) d_h_vcol_phy=0.
	245	c
[524]	246	DO i = 1, klon
	247	cbmf(i) = 0.
[973]	248	sigd(i) = 0.
[524]	249	ENDDO
[879]	250	ENDIF !(ifrst .EQ. 0)
[524]	251
	252	DO k = 1, klev+1
	253	DO i=1,klon
	254	em_ph(i,k) = paprs(i,k) / 100.0
[1334]	255	pmflxr(i,k)=0.
[619]	256	pmflxs(i,k)=0.
[524]	257	ENDDO
	258	ENDDO
	259	c
	260	DO k = 1, klev
	261	DO i=1,klon
	262	em_p(i,k) = pplay(i,k) / 100.0
	263	ENDDO
	264	ENDDO
[879]	265	c
	266	!
	267	! Feeding layer
	268	!
	269	em_sig1feed = 1.
	270	em_sig2feed = 0.97
	271	c em_sig2feed = 0.8
	272	! Relative Weight densities
	273	do k=1,klev
	274	em_wght(k)=1.
	275	end do
	276	cCRtest: couche alim des tehrmiques ponderee par a*
	277	c DO i = 1, klon
	278	c do k=1,lalim_conv(i)
	279	c em_wght(k)=wght_th(i,k)
	280	c print*,'em_wght=',em_wght(k),wght_th(i,k)
	281	c end do
	282	c END DO
[524]	283
	284	if (iflag_con .eq. 4) then
	285	DO k = 1, klev
	286	DO i = 1, klon
	287	zx_t = t(i,k)
	288	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
	289	zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
	290	zcor=1./(1.-retv*zx_qs)
	291	qs(i,k)=zx_qs*zcor
	292	ENDDO
[879]	293	DO i = 1, klon
	294	zx_t = t_wake(i,k)
	295	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
	296	zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
	297	zcor=1./(1.-retv*zx_qs)
	298	qs_wake(i,k)=zx_qs*zcor
	299	ENDDO
[524]	300	ENDDO
	301	else ! iflag_con=3 (modif de puristes qui fait la diffce pour la convergence numerique)
	302	DO k = 1, klev
	303	DO i = 1, klon
	304	zx_t = t(i,k)
	305	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
	306	zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
	307	zx_qs= MIN(0.5,zx_qs)
	308	zcor=1./(1.-retv*zx_qs)
	309	zx_qs=zx_qs*zcor
	310	qs(i,k)=zx_qs
	311	ENDDO
[879]	312	DO i = 1, klon
	313	zx_t = t_wake(i,k)
	314	zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
	315	zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
	316	zx_qs= MIN(0.5,zx_qs)
	317	zcor=1./(1.-retv*zx_qs)
	318	zx_qs=zx_qs*zcor
	319	qs_wake(i,k)=zx_qs
	320	ENDDO
[524]	321	ENDDO
	322	endif ! iflag_con
	323	c
	324	C------------------------------------------------------------------
	325
	326	C Main driver for convection:
[879]	327	C iflag_con=3 -> nvlle version de KE (JYG)
	328	C iflag_con = 30 -> equivalent to convect3
[524]	329	C iflag_con = 4 -> equivalent to convect1/2
[879]	330	c
	331	c
	332	if (iflag_con.eq.30) then
[524]	333
	334	CALL cv_driver(klon,klev,klev+1,ntra,iflag_con,
	335	: t,q,qs,u,v,tra,
	336	$ em_p,em_ph,iflag,
	337	$ d_t,d_q,d_u,d_v,d_tra,rain,
[1334]	338	!! $ pmflxr,cbmf,work1,work2, !jyg
	339	$ Vprecip,cbmf,work1,work2, !jyg
[619]	340	$ kbas,ktop,
	341	$ dtime,Ma,upwd,dnwd,dnwdbis,qcondc,wd,cape,
	342	$ da,phi,mp)
[1398]	343	c
	344	DO i = 1,klon
	345	cbmf(i) = Ma(i,kbas(i))
	346	ENDDO
	347	c
[879]	348	else
[524]	349
[940]	350	cLF necessary for gathered fields
	351	nloc=klon
	352	CALL cva_driver(klon,klev,klev+1,ntra,nloc,
[879]	353	$ iflag_con,iflag_mix,iflag_clos,dtime,
[1146]	354	: t,q,qs,t_wake,q_wake,qs_wake,s_wake,u,v,tra,
[879]	355	$ em_p,em_ph,
	356	. ALE,ALP,
	357	. em_sig1feed,em_sig2feed,em_wght,
	358	. iflag,d_t,d_q,d_u,d_v,d_tra,rain,kbas,ktop,
	359	$ cbmf,work1,work2,ptop2,sigd,
	360	$ Ma,mip,Vprecip,upwd,dnwd,dnwdbis,qcondc,wd,
	361	$ cape,cin,tvp,
	362	$ dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
	363	$ asupmaxmin,lalim_conv)
	364	endif
[524]	365	C------------------------------------------------------------------
	366
	367	DO i = 1,klon
	368	rain(i) = rain(i)/86400.
	369	rflag(i)=iflag(i)
	370	ENDDO
	371
	372	DO k = 1, klev
	373	DO i = 1, klon
	374	d_t(i,k) = dtime*d_t(i,k)
	375	d_q(i,k) = dtime*d_q(i,k)
	376	d_u(i,k) = dtime*d_u(i,k)
	377	d_v(i,k) = dtime*d_v(i,k)
	378	ENDDO
	379	ENDDO
[879]	380	c
	381	if (iflag_con.eq.30) then
[619]	382	DO itra = 1,ntra
	383	DO k = 1, klev
	384	DO i = 1, klon
	385	d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
	386	ENDDO
	387	ENDDO
[879]	388	ENDDO
	389	endif
	390
	391	DO k = 1, klev
	392	DO i = 1, klon
	393	t1(i,k) = t(i,k)+ d_t(i,k)
	394	q1(i,k) = q(i,k)+ d_q(i,k)
	395	ENDDO
	396	ENDDO
[1334]	397	c !jyg
	398	c--Separation neige/pluie (pour diagnostics) !jyg
	399	DO k = 1, klev !jyg
	400	DO i = 1, klon !jyg
	401	IF (t1(i,k).LT.RTT) THEN !jyg
	402	pmflxs(i,k)=Vprecip(i,k) !jyg
	403	ELSE !jyg
	404	pmflxr(i,k)=Vprecip(i,k) !jyg
	405	ENDIF !jyg
	406	ENDDO !jyg
	407	ENDDO !jyg
[879]	408	c
	409	cc IF (if_ebil.ge.2) THEN
	410	cc ztit='after convect'
	411	cc CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
	412	cc e , t1,q1,ql,qs,u,v,paprs,pplay
	413	cc s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	414	cc call diagphy(paire,ztit,ip_ebil
	415	cc e , zero_v, zero_v, zero_v, zero_v, zero_v
	416	cc e , zero_v, rain, zero_v, ztsol
	417	cc e , d_h_vcol, d_qt, d_ec
	418	cc s , fs_bound, fq_bound )
	419	cc END IF
	420	C
	421	c
[524]	422	c les traceurs ne sont pas mis dans cette version de convect4:
	423	if (iflag_con.eq.4) then
	424	DO itra = 1,ntra
	425	DO k = 1, klev
	426	DO i = 1, klon
	427	d_tra(i,k,itra) = 0.
	428	ENDDO
	429	ENDDO
	430	ENDDO
	431	endif
[938]	432	c print*, 'concvl->: dd_t,dd_q ',dd_t(1,1),dd_q(1,1)
[879]	433
[970]	434	DO k = 1, klev
	435	DO i = 1, klon
	436	dtvpdt1(i,k) = 0.
	437	dtvpdq1(i,k) = 0.
	438	ENDDO
	439	ENDDO
	440	DO i = 1, klon
	441	dplcldt(i) = 0.
	442	dplcldr(i) = 0.
	443	ENDDO
[973]	444	c
	445	if(prt_level.GE.20) THEN
	446	DO k=1,klev
	447	! print*,'physiq apres_add_con i k it d_u d_v d_t d_q qdl0',igout
	448	! .,k,itap,d_u_con(igout,k) ,d_v_con(igout,k), d_t_con(igout,k),
	449	! .d_q_con(igout,k),dql0(igout,k)
	450	! print*,'phys apres_add_con itap Ma cin ALE ALP wak t q undi t q'
	451	! .,itap,Ma(igout,k),cin(igout),ALE(igout), ALP(igout),
	452	! . t_wake(igout,k),q_wake(igout,k),t_undi(igout,k),q_undi(igout,k)
	453	! print*,'phy apres_add_con itap CON rain snow EMA wk1 wk2 Vpp mip'
	454	! .,itap,rain_con(igout),snow_con(igout),ema_work1(igout,k),
	455	! .ema_work2(igout,k),Vprecip(igout,k), mip(igout,k)
	456	! print*,'phy apres_add_con itap upwd dnwd dnwd0 cape tvp Tconv '
	457	! .,itap,upwd(igout,k),dnwd(igout,k),dnwd0(igout,k),cape(igout),
	458	! .tvp(igout,k),Tconv(igout,k)
	459	! print*,'phy apres_add_con itap dtvpdt dtvdq dplcl dplcldr qcondc'
	460	! .,itap,dtvpdt1(igout,k),dtvpdq1(igout,k),dplcldt(igout),
	461	! .dplcldr(igout),qcondc(igout,k)
	462	! print*,'phy apres_add_con itap wd pmflxr Kpmflxr Kp1 Kpmflxs Kp1'
	463	! .,itap,wd(igout),pmflxr(igout,k),pmflxr(igout,k+1),pmflxs(igout,k)
	464	! .,pmflxs(igout,k+1)
	465	! print*,'phy apres_add_con itap da phi mp ftd fqd lalim wgth',
	466	! .itap,da(igout,k),phi(igout,k,k),mp(igout,k),ftd(igout,k),
	467	! . fqd(igout,k),lalim_conv(igout),wght_th(igout,k)
	468	ENDDO
	469	endif !(prt_level.EQ.20) THEN
	470	c
[524]	471	RETURN
	472	END
	473

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

Download in other formats:

Original Format