Context Navigation

convect2.F @ 5306

Last change on this file since 5306 was 302, checked in by lmdz, 23 years ago
Correction bug dans KE vectorise (FH, JYG) LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 44.9 KB

Rev	Line
[302]	1	c
	2	c $Header$
	3	c
[254]	4	subroutine convect2(ncum,idcum,len,nd,ndp1,nl,minorig,
	5	& nk1,icb1,
	6	& t1,q1,qs1,u1,v1,gz1,tv1,tp1,tvp1,clw1,h1,
	7	& lv1,cpn1,p1,ph1,ft1,fq1,fu1,fv1,
	8	& tnk1,qnk1,gznk1,plcl1,
	9	& precip1,cbmf1,iflag1,
	10	& delt,cpd,cpv,cl,rv,rd,lv0,g,
	11	& sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp,
	12	& alpha,entp,coeffs,coeffr,omtrain,cu,Ma)
	13	C.............................START PROLOGUE............................
	14	C
	15	C SCCS IDENTIFICATION: @(#)convect2.f 1.2 05/18/00
	16	C 22:06:22 /h/cm/library/nogaps4/src/sub/fcst/convect2.f_v
	17	C
	18	C CONFIGURATION IDENTIFICATION: None
	19	C
	20	C MODULE NAME: convect2
	21	C
	22	C DESCRIPTION:
	23	C
	24	C convect1 The Emanuel Cumulus Convection Scheme - compute tendencies
	25	C
	26	C CONTRACT NUMBER AND TITLE: None
	27	C
	28	C REFERENCES: Programmers K. Emanuel (MIT), Timothy F. Hogan, M. Peng (NRL)
	29	C
	30	C CLASSIFICATION: Unclassified
	31	C
	32	C RESTRICTIONS: None
	33	C
	34	C COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90
	35	C
	36	C COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress"
	37	C Fortran 90: -O vector3,scalar3,task1,aggress,overindex -ei -r 2
	38	C
	39	C LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a
	40	C
	41	C USAGE: call convect2(ncum,idcum,len,nd,nl,minorig,
	42	C & nk1,icb1,
	43	C & t1,q1,qs1,u1,v1,gz1,tv1,tp1,tvp1,clw1,h1,
	44	C & lv1,cpn1,p1,ph1,ft1,fq1,fu1,fv1,
	45	C & tnk1,qnk1,gznk1,plcl1,
	46	C & precip1,cbmf1,iflag1,
	47	C & delt,cpd,cpv,cl,rv,rd,lv0,g,
	48	C & sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp,
	49	C & alpha,entp,coeffs,coeffr,omtrain,cu)
	50	C
	51	C PARAMETERS:
	52	C Name Type Usage Description
	53	C ---------- ---------- ------- ----------------------------
	54	C
	55	C ncum Integer Input number of cumulus points
	56	C idcum Integer Input index of cumulus point
	57	C len Integer Input first dimension
	58	C nd Integer Input total vertical dimension
	59	C ndp1 Integer Input nd + 1
	60	C nl Integer Input vertical dimension for convection
	61	C minorig Integer Input First level where convection is allow to begin
	62	C nk1 Integer Input First level of convection
	63	C ncb1 Integer Input Level of free convection
	64	C t1 Real Input temperature
	65	C q1 Real Input specific hum
	66	C qs1 Real Input sat specific hum
	67	C u1 Real Input u-wind
	68	C v1 Real Input v-wind
	69	C gz1 Real Inout geop
	70	C tv1 Real Input virtual temp
	71	C tp1 Real Input
	72	C clw1 Real Inout cloud liquid water
	73	C h1 Real Inout
	74	C lv1 Real Inout
	75	C cpn1 Real Inout
	76	C p1 Real Input full level pressure
	77	C ph1 Real Input half level pressure
	78	C ft1 Real Output temp tend
	79	C fq1 Real Output spec hum tend
	80	C fu1 Real Output u-wind tend
	81	C fv1 Real Output v-wind tend
	82	C precip1 Real Output prec
	83	C cbmf1 Real In/Out cumulus mass flux
	84	C iflag1 Integer Output iflag on latitude strip
	85	C delt Real Input time step
	86	C cpd Integer Input See description below
	87	C cpv Integer Input See description below
	88	C cl Integer Input See description below
	89	C rv Integer Input See description below
	90	C rd Integer Input See description below
	91	C lv0 Integer Input See description below
	92	C g Integer Input See description below
	93	C sigs Integer Input See description below
	94	C sigd Integer Input See description below
	95	C elcrit Integer Input See description below
	96	C tlcrit Integer Input See description below
	97	C omtsnow Integer Input See description below
	98	C dtmax Integer Input See description below
	99	C damp Integer Input See description below
	100	C alpha Integer Input See description below
	101	C ent Integer Input See description below
	102	C coeffs Integer Input See description below
	103	C coeffr Integer Input See description below
	104	C omtrain Integer Input See description below
	105	C cu Integer Input See description below
	106	C
	107	C COMMON BLOCKS:
	108	C Block Name Type Usage Notes
	109	C -------- -------- ---- ------ ------------------------
	110	C
	111	C FILES: None
	112	C
	113	C DATA BASES: None
	114	C
	115	C NON-FILE INPUT/OUTPUT: None
	116	C
	117	C ERROR CONDITIONS: None
	118	C
	119	C ADDITIONAL COMMENTS: None
	120	C
	121	C.................MAINTENANCE SECTION................................
	122	C
	123	C MODULES CALLED:
	124	C Name Description
	125	C zilch Zero out an array
	126	C ------- ----------------------
	127	C LOCAL VARIABLES AND
	128	C STRUCTURES:
	129	C Name Type Description
	130	C ------- ------ -----------
	131	C See Comments Below
	132	C
	133	C i Integer loop index
	134	C k Integer loop index
	135	c
	136	C METHOD:
	137	C
	138	C See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation of a
	139	C convective scheme for use in climate models.
	140	C
	141	C FILES: None
	142	C
	143	C INCLUDE FILES: None
	144	C
	145	C MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak
	146	C
	147	C..............................END PROLOGUE.............................
	148	c
	149	c
	150	implicit none
	151	c
	152	#include "dimensions.h"
	153	#include "dimphy.h"
	154	c
	155	integer kmax2,imax2,kmin2,imin2
	156	real ftmax2,ftmin2
	157	integer kmax,imax,kmin,imin
	158	real ftmax,ftmin
	159	c
	160	integer ncum
	161	integer idcum(len)
	162	integer len
	163	integer nd
	164	integer ndp1
	165	integer nl
	166	integer minorig
	167	integer nk1(len)
	168	integer icb1(len)
	169	real t1(len,nd)
	170	real q1(len,nd)
	171	real qs1(len,nd)
	172	real u1(len,nd)
	173	real v1(len,nd)
	174	real gz1(len,nd)
	175	real tv1(len,nd)
	176	real tp1(len,nd)
	177	real tvp1(len,nd)
	178	real clw1(len,nd)
	179	real h1(len,nd)
	180	real lv1(len,nd)
	181	real cpn1(len,nd)
	182	real p1(len,nd)
	183	real ph1(len,ndp1)
	184	real ft1(len,nd)
	185	real fq1(len,nd)
	186	real fu1(len,nd)
	187	real fv1(len,nd)
	188	real tnk1(len)
	189	real qnk1(len)
	190	real gznk1(len)
	191	real precip1(len)
	192	real cbmf1(len)
	193	real plcl1(len)
	194	integer iflag1(len)
	195	real delt
	196	real cpd
	197	real cpv
	198	real cl
	199	real rv
	200	real rd
	201	real lv0
	202	real g
	203	real sigs ! SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE
	204	real sigd ! SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT
	205	real elcrit ! ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm)
	206	real tlcrit ! TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-
	207	c CONVERSION THRESHOLD IS ASSUMED TO BE ZERO
	208	real omtsnow ! OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW
	209	real dtmax ! DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION
	210	c A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC.
	211	real damp
	212	real alpha
	213	real entp ! ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT FORMULATION
	214	real coeffs ! COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION OF SNOW
	215	real coeffr ! COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION OF RAIN
	216	real omtrain ! OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN
	217	real cu ! CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM TRANSPORT
	218	c
	219	real Ma(len,nd)
	220	c
	221	C
	222	C * ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) *
	223	C * TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- *
	224	C * CONVERSION THRESHOLD IS ASSUMED TO BE ZERO *
	225	C * (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY *
	226	C * BETWEEN 0 C AND TLCRIT) *
	227	C * ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT *
	228	C * FORMULATION *
	229	C * SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *
	230	C * SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *
	231	C * OF CLOUD *
	232	C * OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN *
	233	C * OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW *
	234	C * COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *
	235	C * OF RAIN *
	236	C * COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *
	237	C * OF SNOW *
	238	C * CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM *
	239	C * TRANSPORT *
	240	C * DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION *
	241	C * A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC *
	242	C * ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF *
	243	C * APPROACH TO QUASI-EQUILIBRIUM *
	244	C * (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) *
	245	C * (DAMP MUST BE LESS THAN 1) *
	246	c
	247	c Local arrays.
	248	c
	249	real work(ncum)
	250	real t(ncum,klev)
	251	real q(ncum,klev)
	252	real qs(ncum,klev)
	253	real u(ncum,klev)
	254	real v(ncum,klev)
	255	real gz(ncum,klev)
	256	real h(ncum,klev)
	257	real lv(ncum,klev)
	258	real cpn(ncum,klev)
	259	real p(ncum,klev)
	260	real ph(ncum,klev)
	261	real ft(ncum,klev)
	262	real fq(ncum,klev)
	263	real fu(ncum,klev)
	264	real fv(ncum,klev)
	265	real precip(ncum)
	266	real cbmf(ncum)
	267	real plcl(ncum)
	268	real tnk(ncum)
	269	real qnk(ncum)
	270	real gznk(ncum)
	271	real tv(ncum,klev)
	272	real tp(ncum,klev)
	273	real tvp(ncum,klev)
	274	real clw(ncum,klev)
	275	c real det(ncum,klev)
	276	real dph(ncum,klev)
	277	c real wd(ncum)
	278	c real tprime(ncum)
	279	c real qprime(ncum)
	280	real ah0(ncum)
	281	real ep(ncum,klev)
	282	real sigp(ncum,klev)
	283	integer nent(ncum,klev)
	284	real water(ncum,klev)
	285	real evap(ncum,klev)
	286	real mp(ncum,klev)
	287	real m(ncum,klev)
	288	real qti
	289	real wt(ncum,klev)
	290	real hp(ncum,klev)
	291	real lvcp(ncum,klev)
	292	real elij(ncum,klev,klev)
	293	real ment(ncum,klev,klev)
	294	real sij(ncum,klev,klev)
	295	real qent(ncum,klev,klev)
	296	real uent(ncum,klev,klev)
	297	real vent(ncum,klev,klev)
	298	real qp(ncum,klev)
	299	real up(ncum,klev)
	300	real vp(ncum,klev)
	301	real cape(ncum)
	302	real capem(ncum)
	303	real frac(ncum)
	304	real dtpbl(ncum)
	305	real tvpplcl(ncum)
	306	real tvaplcl(ncum)
	307	real dtmin(ncum)
	308	real w3d(ncum,klev)
	309	real am(ncum)
	310	real ents(ncum)
	311	real uav(ncum)
	312	real vav(ncum)
	313	c
	314	integer iflag(ncum)
	315	integer nk(ncum)
	316	integer icb(ncum)
	317	integer inb(ncum)
	318	integer inb1(ncum)
	319	integer jtt(ncum)
	320	c
	321	integer nn,i,k,n,icbmax,nlp,j
	322	integer ij
	323	integer nn2,nn3
	324	real clmcpv
	325	real clmcpd
	326	real cpdmcp
	327	real cpvmcpd
	328	real eps
	329	real epsi
	330	real epsim1
	331	real tg,qg,s,alv,tc,ahg,denom,es,rg,ginv,rowl
	332	real delti
	333	real tca,elacrit
	334	real by,defrac
	335	c real byp
	336	real byp(ncum)
	337	logical lcape(ncum)
	338	real dbo
	339	real bf2,anum,dei,altem,cwat,stemp
	340	real alt,qp1,smid,sjmax,sjmin
	341	real delp,delm
	342	real awat,coeff,afac,revap,dhdp,fac,qstm,rat
	343	real qsm,sigt,b6,c6
	344	real dpinv,cpinv
	345	real fqold,ftold,fuold,fvold
	346	real wdtrain(ncum),xxx
	347	real bsum(ncum,klev)
	348	real asij(ncum)
	349	real smin(ncum)
	350	real scrit(ncum)
	351	c real amp1,ad
	352	real amp1(ncum),ad(ncum)
	353	logical lwork(ncum)
	354	integer num1,num2
	355	c
	356	c print*,'cpd en entree de convect2 ',cpd
	357	nlp=nl+1
	358	c
	359	rowl=1000.0
	360	ginv=1.0/g
	361	delti=1.0/delt
	362	c
	363	c Define some thermodynamic variables.
	364	c
	365	clmcpv=cl-cpv
	366	clmcpd=cl-cpd
	367	cpdmcp=cpd-cpv
	368	cpvmcpd=cpv-cpd
	369	eps=rd/rv
	370	epsi=1.0/eps
	371	epsim1=epsi-1.0
	372	c
	373	c Compress the fields.
	374	c
	375	do 110 k=1,nl+1
	376	nn=0
	377	do 100 i=1,len
	378	if(iflag1(i).eq.0)then
	379	nn=nn+1
	380	t(nn,k)=t1(i,k)
	381	q(nn,k)=q1(i,k)
	382	qs(nn,k)=qs1(i,k)
	383	u(nn,k)=u1(i,k)
	384	v(nn,k)=v1(i,k)
	385	gz(nn,k)=gz1(i,k)
	386	h(nn,k)=h1(i,k)
	387	lv(nn,k)=lv1(i,k)
	388	cpn(nn,k)=cpn1(i,k)
	389	p(nn,k)=p1(i,k)
	390	ph(nn,k)=ph1(i,k)
	391	tv(nn,k)=tv1(i,k)
	392	tp(nn,k)=tp1(i,k)
	393	tvp(nn,k)=tvp1(i,k)
	394	clw(nn,k)=clw1(i,k)
	395	endif
	396	100 continue
	397	c print*,'100 ncum,nn',ncum,nn
	398	110 continue
	399	nn=0
	400	do 150 i=1,len
	401	if(iflag1(i).eq.0)then
	402	nn=nn+1
	403	cbmf(nn)=cbmf1(i)
	404	plcl(nn)=plcl1(i)
	405	tnk(nn)=tnk1(i)
	406	qnk(nn)=qnk1(i)
	407	gznk(nn)=gznk1(i)
	408	nk(nn)=nk1(i)
	409	icb(nn)=icb1(i)
	410	iflag(nn)=iflag1(i)
	411	endif
	412	150 continue
	413	c print*,'150 ncum,nn',ncum,nn
	414	c
	415	c Initialize the tendencies, det, wd, tprime, qprime.
	416	c
	417	do 170 k=1,nl
	418	do 160 i=1,ncum
	419	c det(i,k)=0.0
	420	ft(i,k)=0.0
	421	fu(i,k)=0.0
	422	fv(i,k)=0.0
	423	fq(i,k)=0.0
	424	dph(i,k)=ph(i,k)-ph(i,k+1)
	425	ep(i,k)=0.0
	426	sigp(i,k)=sigs
	427	160 continue
	428	170 continue
	429	do 180 i=1,ncum
	430	c wd(i)=0.0
	431	c tprime(i)=0.0
	432	c qprime(i)=0.0
	433	precip(i)=0.0
	434	ft(i,nl+1)=0.0
	435	fu(i,nl+1)=0.0
	436	fv(i,nl+1)=0.0
	437	fq(i,nl+1)=0.0
	438	180 continue
	439	c
	440	c Compute icbmax.
	441	c
	442	icbmax=2
	443	do 230 i=1,ncum
	444	icbmax=max(icbmax,icb(i))
	445	230 continue
	446	c
	447	c
	448	!=====================================================================
	449	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	450	!=====================================================================
	451	c
	452	c --- The procedure is to solve the equation.
	453	c cptp+Lqp+phi=cptnk+Lqnk+gznk.
	454	c
	455	c * Calculate certain parcel quantities, including static energy *
	456	c
	457	c
	458	do 240 i=1,ncum
	459	ah0(i)=(cpd(1.-qnk(i))+clqnk(i))*tnk(i)
	460	& +qnk(i)(lv0-clmcpv(tnk(i)-273.15))+gznk(i)
	461	240 continue
	462	c
	463	c
	464	c * Find lifted parcel quantities above cloud base *
	465	c
	466	c
	467	do 300 k=minorig+1,nl
	468	do 290 i=1,ncum
	469	if(k.ge.(icb(i)+1))then
	470	tg=t(i,k)
	471	qg=qs(i,k)
	472	alv=lv0-clmcpv*(t(i,k)-273.15)
	473	c
	474	c First iteration.
	475	c
	476	s=cpd+alvalvqg/(rvt(i,k)t(i,k))
	477	s=1./s
	478	ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
	479	tg=tg+s*(ah0(i)-ahg)
	480	tg=max(tg,35.0)
	481	tc=tg-273.15
	482	denom=243.5+tc
	483	if(tc.ge.0.0)then
	484	es=6.112exp(17.67tc/denom)
	485	else
	486	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	487	endif
	488	qg=epses/(p(i,k)-es(1.-eps))
	489	c
	490	c Second iteration.
	491	c
	492	s=cpd+alvalvqg/(rvt(i,k)t(i,k))
	493	s=1./s
	494	ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
	495	tg=tg+s*(ah0(i)-ahg)
	496	tg=max(tg,35.0)
	497	tc=tg-273.15
	498	denom=243.5+tc
	499	if(tc.ge.0.0)then
	500	es=6.112exp(17.67tc/denom)
	501	else
	502	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	503	endif
	504	qg=epses/(p(i,k)-es(1.-eps))
	505	c
	506	alv=lv0-clmcpv*(t(i,k)-273.15)
	507	c print*,'cpd dans convect2 ',cpd
	508	c print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
	509	c print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
	510	tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)
	511	& /cpd
	512	c if (.not.cpd.gt.1000.) then
	513	c print*,'CPD=',cpd
	514	c stop
	515	c endif
	516	clw(i,k)=qnk(i)-qg
	517	clw(i,k)=max(0.0,clw(i,k))
	518	rg=qg/(1.-qnk(i))
	519	tvp(i,k)=tp(i,k)(1.+rgepsi)
	520	endif
	521	290 continue
	522	300 continue
	523	c
	524	!=====================================================================
	525	! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
	526	! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
	527	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
	528	!=====================================================================
	529	c
	530	do 320 k=minorig+1,nl
	531	do 310 i=1,ncum
	532	if(k.ge.(nk(i)+1))then
	533	tca=tp(i,k)-273.15
	534	if(tca.ge.0.0)then
	535	elacrit=elcrit
	536	else
	537	elacrit=elcrit*(1.0-tca/tlcrit)
	538	endif
	539	elacrit=max(elacrit,0.0)
	540	ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
	541	ep(i,k)=max(ep(i,k),0.0 )
	542	ep(i,k)=min(ep(i,k),1.0 )
	543	sigp(i,k)=sigs
	544	endif
	545	310 continue
	546	320 continue
	547	c
	548	!=====================================================================
	549	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
	550	! --- VIRTUAL TEMPERATURE
	551	!=====================================================================
	552	c
	553	do 340 k=minorig+1,nl
	554	do 330 i=1,ncum
	555	if(k.ge.(icb(i)+1))then
	556	tvp(i,k)=tvp(i,k)(1.0-qnk(i)+ep(i,k)clw(i,k))
	557	c print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
	558	c print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
	559	endif
	560	330 continue
	561	340 continue
	562	do 350 i=1,ncum
	563	tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
	564	350 continue
	565	c
	566	c
	567	c=====================================================================
	568	c --- NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	569	c=====================================================================
	570	c
	571	do 360 i=1,ncum*nlp
	572	nent(i,1)=0
	573	water(i,1)=0.0
	574	evap(i,1)=0.0
	575	mp(i,1)=0.0
	576	m(i,1)=0.0
	577	wt(i,1)=omtsnow
	578	hp(i,1)=h(i,1)
	579	c if(.not.cpn(i,1).gt.900.) then
	580	c print*,'i,lv(i,1),cpn(i,1)'
	581	c print*, i,lv(i,1),cpn(i,1)
	582	c k=(i-1)/ncum+1
	583	c print*,'i,k',mod(i,ncum),k,' cpn',cpn(mod(i,ncum),k)
	584	c stop
	585	c endif
	586	lvcp(i,1)=lv(i,1)/cpn(i,1)
	587	360 continue
	588	c
	589	do 380 i=1,ncumnlpnlp
	590	380 continue
	591	c
	592	do 400 k=1,nlp
	593	do 390 j=1,nlp
	594	do 385 i=1,ncum
	595	qent(i,k,j)=q(i,j)
	596	uent(i,k,j)=u(i,j)
	597	vent(i,k,j)=v(i,j)
[301]	598	elij(i,k,j)=0.0
	599	ment(i,k,j)=0.0
	600	sij(i,k,j)=0.0
[254]	601	385 continue
	602	390 continue
	603	400 continue
	604	c
	605	do 420 i=1,ncum
	606	qp(i,1)=q(i,1)
	607	up(i,1)=u(i,1)
	608	vp(i,1)=v(i,1)
	609	420 continue
	610	do 440 k=2,nlp
	611	do 430 i=1,ncum
	612	qp(i,k)=q(i,k-1)
	613	up(i,k)=u(i,k-1)
	614	vp(i,k)=v(i,k-1)
	615	430 continue
	616	440 continue
	617	c
	618	c=====================================================================
	619	c --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
	620	c --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
	621	c --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
	622	c=====================================================================
	623	c
	624	do 510 i=1,ncum
	625	cape(i)=0.0
	626	capem(i)=0.0
	627	inb(i)=icb(i)+1
	628	inb1(i)=inb(i)
	629	510 continue
	630	c
	631	c Originial Code
	632	c
	633	c do 530 k=minorig+1,nl-1
	634	c do 520 i=1,ncum
	635	c if(k.ge.(icb(i)+1))then
	636	c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	637	c byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	638	c cape(i)=cape(i)+by
	639	c if(by.ge.0.0)inb1(i)=k+1
	640	c if(cape(i).gt.0.0)then
	641	c inb(i)=k+1
	642	c capem(i)=cape(i)
	643	c endif
	644	c endif
	645	c520 continue
	646	c530 continue
	647	c do 540 i=1,ncum
	648	c byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
	649	c cape(i)=capem(i)+byp
	650	c defrac=capem(i)-cape(i)
	651	c defrac=max(defrac,0.001)
	652	c frac(i)=-cape(i)/defrac
	653	c frac(i)=min(frac(i),1.0)
	654	c frac(i)=max(frac(i),0.0)
	655	c540 continue
	656	c
	657	c K Emanuel fix
	658	c
	659	c call zilch(byp,ncum)
	660	c do 530 k=minorig+1,nl-1
	661	c do 520 i=1,ncum
	662	c if(k.ge.(icb(i)+1))then
	663	c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	664	c cape(i)=cape(i)+by
	665	c if(by.ge.0.0)inb1(i)=k+1
	666	c if(cape(i).gt.0.0)then
	667	c inb(i)=k+1
	668	c capem(i)=cape(i)
	669	c byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	670	c endif
	671	c endif
	672	c520 continue
	673	c530 continue
	674	c do 540 i=1,ncum
	675	c inb(i)=max(inb(i),inb1(i))
	676	c cape(i)=capem(i)+byp(i)
	677	c defrac=capem(i)-cape(i)
	678	c defrac=max(defrac,0.001)
	679	c frac(i)=-cape(i)/defrac
	680	c frac(i)=min(frac(i),1.0)
	681	c frac(i)=max(frac(i),0.0)
	682	c540 continue
	683	c
	684	c J Teixeira fix
	685	c
	686	call zilch(byp,ncum)
	687	do 515 i=1,ncum
	688	lcape(i)=.true.
	689	515 continue
	690	do 530 k=minorig+1,nl-1
	691	do 520 i=1,ncum
	692	if(cape(i).lt.0.0)lcape(i)=.false.
	693	if((k.ge.(icb(i)+1)).and.lcape(i))then
	694	by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	695	byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	696	cape(i)=cape(i)+by
	697	if(by.ge.0.0)inb1(i)=k+1
	698	if(cape(i).gt.0.0)then
	699	inb(i)=k+1
	700	capem(i)=cape(i)
	701	endif
	702	endif
	703	520 continue
	704	530 continue
	705	do 540 i=1,ncum
	706	cape(i)=capem(i)+byp(i)
	707	defrac=capem(i)-cape(i)
	708	defrac=max(defrac,0.001)
	709	frac(i)=-cape(i)/defrac
	710	frac(i)=min(frac(i),1.0)
	711	frac(i)=max(frac(i),0.0)
	712	540 continue
	713	c
	714	c=====================================================================
	715	c --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
	716	c=====================================================================
	717	c
	718	do 600 k=minorig+1,nl
	719	do 590 i=1,ncum
	720	if((k.ge.icb(i)).and.(k.le.inb(i)))then
	721	hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
	722	endif
	723	590 continue
	724	600 continue
	725	c
	726	c=====================================================================
	727	c --- CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I),
	728	c --- AT EACH MODEL LEVEL
	729	c=====================================================================
	730	c
	731	c tvpplcl = parcel temperature lifted adiabatically from level
	732	c icb-1 to the LCL.
	733	c tvaplcl = virtual temperature at the LCL.
	734	c
	735	do 610 i=1,ncum
	736	dtpbl(i)=0.0
	737	tvpplcl(i)=tvp(i,icb(i)-1)
	738	& -rdtvp(i,icb(i)-1)(p(i,icb(i)-1)-plcl(i))
	739	& /(cpn(i,icb(i)-1)*p(i,icb(i)-1))
	740	tvaplcl(i)=tv(i,icb(i))
	741	& +(tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i,icb(i)))
	742	& /(p(i,icb(i))-p(i,icb(i)+1))
	743	610 continue
	744	c
	745	c-------------------------------------------------------------------
	746	c --- Interpolate difference between lifted parcel and
	747	c --- environmental temperatures to lifted condensation level
	748	c-------------------------------------------------------------------
	749	c
	750	c dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
	751	c
	752	do 630 k=minorig,icbmax
	753	do 620 i=1,ncum
	754	if((k.ge.nk(i)).and.(k.le.(icb(i)-1)))then
	755	dtpbl(i)=dtpbl(i)+(tvp(i,k)-tv(i,k))*dph(i,k)
	756	endif
	757	620 continue
	758	630 continue
	759	do 640 i=1,ncum
	760	dtpbl(i)=dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
	761	dtmin(i)=tvpplcl(i)-tvaplcl(i)+dtmax+dtpbl(i)
	762	640 continue
	763	c
	764	c-------------------------------------------------------------------
	765	c --- Adjust cloud base mass flux
	766	c-------------------------------------------------------------------
	767	c
	768	do 650 i=1,ncum
	769	work(i)=cbmf(i)
	770	cbmf(i)=max(0.0,(1.0-damp)cbmf(i)+0.1alpha*dtmin(i))
	771	if((work(i).eq.0.0).and.(cbmf(i).eq.0.0))then
	772	iflag(i)=3
	773	endif
	774	650 continue
	775	c
	776	c-------------------------------------------------------------------
	777	c --- Calculate rates of mixing, m(i)
	778	c-------------------------------------------------------------------
	779	c
	780	call zilch(work,ncum)
	781	c
	782	do 670 j=minorig+1,nl
	783	do 660 i=1,ncum
	784	if((j.ge.(icb(i)+1)).and.(j.le.inb(i)))then
	785	k=min(j,inb1(i))
	786	dbo=abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1))
	787	& +entp0.04(ph(i,k)-ph(i,k+1))
	788	work(i)=work(i)+dbo
	789	m(i,j)=cbmf(i)*dbo
	790	endif
	791	660 continue
	792	670 continue
	793	do 690 k=minorig+1,nl
	794	do 680 i=1,ncum
	795	if((k.ge.(icb(i)+1)).and.(k.le.inb(i)))then
	796	m(i,k)=m(i,k)/work(i)
	797	endif
	798	680 continue
	799	690 continue
	800	c
	801	c
	802	c=====================================================================
	803	c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
	804	c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	805	c --- FRACTION (sij)
	806	c=====================================================================
	807	c
	808	c
	809	do 750 i=minorig+1,nl
	810	do 710 j=minorig+1,nl
	811	do 700 ij=1,ncum
	812	if((i.ge.(icb(ij)+1)).and.(j.ge.icb(ij))
	813	& .and.(i.le.inb(ij)).and.(j.le.inb(ij)))then
	814	qti=qnk(ij)-ep(ij,i)*clw(ij,i)
	815	bf2=1.+lv(ij,j)lv(ij,j)qs(ij,j)
	816	& /(rvt(ij,j)t(ij,j)*cpd)
	817	anum=h(ij,j)-hp(ij,i)+(cpv-cpd)t(ij,j)(qti-q(ij,j))
	818	denom=h(ij,i)-hp(ij,i)+(cpd-cpv)(q(ij,i)-qti)t(ij,j)
	819	dei=denom
	820	if(abs(dei).lt.0.01)dei=0.01
	821	sij(ij,i,j)=anum/dei
	822	sij(ij,i,i)=1.0
	823	altem=sij(ij,i,j)q(ij,i)+(1.-sij(ij,i,j))qti-qs(ij,j)
	824	altem=altem/bf2
	825	cwat=clw(ij,j)*(1.-ep(ij,j))
	826	stemp=sij(ij,i,j)
	827	if((stemp.lt.0.0.or.stemp.gt.1.0.or.
	828	1 altem.gt.cwat).and.j.gt.i)then
	829	anum=anum-lv(ij,j)(qti-qs(ij,j)-cwatbf2)
	830	denom=denom+lv(ij,j)*(q(ij,i)-qti)
	831	if(abs(denom).lt.0.01)denom=0.01
	832	sij(ij,i,j)=anum/denom
	833	altem=sij(ij,i,j)q(ij,i)+(1.-sij(ij,i,j))qti-qs(ij,j)
	834	altem=altem-(bf2-1.)*cwat
	835	endif
	836	if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
	837	qent(ij,i,j)=sij(ij,i,j)*q(ij,i)
	838	& +(1.-sij(ij,i,j))*qti
	839	uent(ij,i,j)=sij(ij,i,j)*u(ij,i)
	840	& +(1.-sij(ij,i,j))*u(ij,nk(ij))
	841	vent(ij,i,j)=sij(ij,i,j)*v(ij,i)
	842	& +(1.-sij(ij,i,j))*v(ij,nk(ij))
	843	elij(ij,i,j)=altem
	844	elij(ij,i,j)=max(0.0,elij(ij,i,j))
	845	ment(ij,i,j)=m(ij,i)/(1.-sij(ij,i,j))
	846	nent(ij,i)=nent(ij,i)+1
	847	endif
	848	sij(ij,i,j)=max(0.0,sij(ij,i,j))
	849	sij(ij,i,j)=min(1.0,sij(ij,i,j))
	850	endif
	851	700 continue
	852	710 continue
	853	c
	854	c * If no air can entrain at level i assume that updraft detrains *
	855	c * at that level and calculate detrained air flux and properties *
	856	c
	857	do 740 ij=1,ncum
	858	if((i.ge.(icb(ij)+1)).and.(i.le.inb(ij))
	859	& .and.(nent(ij,i).eq.0))then
	860	ment(ij,i,i)=m(ij,i)
	861	qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
	862	uent(ij,i,i)=u(ij,nk(ij))
	863	vent(ij,i,i)=v(ij,nk(ij))
	864	elij(ij,i,i)=clw(ij,i)
	865	sij(ij,i,i)=1.0
	866	endif
	867	740 continue
	868	750 continue
	869	c
	870	do 770 i=1,ncum
	871	sij(i,inb(i),inb(i))=1.0
	872	770 continue
	873	c
	874	c=====================================================================
	875	c --- NORMALIZE ENTRAINED AIR MASS FLUXES
	876	c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	877	c=====================================================================
	878	c
	879	c
	880	call zilch(bsum,ncum*nlp)
	881	do 780 ij=1,ncum
	882	lwork(ij)=.false.
	883	780 continue
	884	do 789 i=minorig+1,nl
	885	c
	886	num1=0
	887	do 779 ij=1,ncum
	888	if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))num1=num1+1
	889	779 continue
	890	if(num1.le.0)go to 789
	891	c
	892	do 781 ij=1,ncum
	893	if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))then
	894	lwork(ij)=(nent(ij,i).ne.0)
	895	qp1=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
	896	anum=h(ij,i)-hp(ij,i)-lv(ij,i)*(qp1-qs(ij,i))
	897	denom=h(ij,i)-hp(ij,i)+lv(ij,i)*(q(ij,i)-qp1)
	898	if(abs(denom).lt.0.01)denom=0.01
	899	scrit(ij)=anum/denom
	900	alt=qp1-qs(ij,i)+scrit(ij)*(q(ij,i)-qp1)
	901	if(scrit(ij).lt.0.0.or.alt.lt.0.0)scrit(ij)=1.0
	902	asij(ij)=0.0
	903	smin(ij)=1.0
	904	endif
	905	781 continue
	906	do 783 j=minorig,nl
	907	c
	908	num2=0
	909	do 778 ij=1,ncum
	910	if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
	911	& .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
	912	& .and.lwork(ij))num2=num2+1
	913	778 continue
	914	if(num2.le.0)go to 783
	915	c
	916	do 782 ij=1,ncum
	917	if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
	918	& .and.(j.ge.icb(ij)).and.(j.le.inb(ij)).and.lwork(ij))then
	919	if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
	920	if(j.gt.i)then
	921	smid=min(sij(ij,i,j),scrit(ij))
	922	sjmax=smid
	923	sjmin=smid
	924	if(smid.lt.smin(ij)
	925	& .and.sij(ij,i,j+1).lt.smid)then
	926	smin(ij)=smid
	927	sjmax=min(sij(ij,i,j+1),sij(ij,i,j),scrit(ij))
	928	sjmin=max(sij(ij,i,j-1),sij(ij,i,j))
	929	sjmin=min(sjmin,scrit(ij))
	930	endif
	931	else
	932	sjmax=max(sij(ij,i,j+1),scrit(ij))
	933	smid=max(sij(ij,i,j),scrit(ij))
	934	sjmin=0.0
	935	if(j.gt.1)sjmin=sij(ij,i,j-1)
	936	sjmin=max(sjmin,scrit(ij))
	937	endif
	938	delp=abs(sjmax-smid)
	939	delm=abs(sjmin-smid)
	940	asij(ij)=asij(ij)+(delp+delm)
	941	& *(ph(ij,j)-ph(ij,j+1))
	942	ment(ij,i,j)=ment(ij,i,j)*(delp+delm)
	943	& *(ph(ij,j)-ph(ij,j+1))
	944	endif
	945	endif
	946	782 continue
	947	783 continue
	948	do 784 ij=1,ncum
	949	if((i.ge.icb(ij)+1).and.(i.le.inb(ij)).and.lwork(ij))then
	950	asij(ij)=max(1.0e-21,asij(ij))
	951	asij(ij)=1.0/asij(ij)
	952	bsum(ij,i)=0.0
	953	endif
	954	784 continue
	955	do 786 j=minorig,nl+1
	956	do 785 ij=1,ncum
	957	if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
	958	& .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
	959	& .and.lwork(ij))then
	960	ment(ij,i,j)=ment(ij,i,j)*asij(ij)
	961	bsum(ij,i)=bsum(ij,i)+ment(ij,i,j)
	962	endif
	963	785 continue
	964	786 continue
	965	do 787 ij=1,ncum
	966	if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
	967	& .and.(bsum(ij,i).lt.1.0e-18).and.lwork(ij))then
	968	nent(ij,i)=0
	969	ment(ij,i,i)=m(ij,i)
	970	qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
	971	uent(ij,i,i)=u(ij,nk(ij))
	972	vent(ij,i,i)=v(ij,nk(ij))
	973	elij(ij,i,i)=clw(ij,i)
	974	sij(ij,i,i)=1.0
	975	endif
	976	787 continue
	977	789 continue
	978	c
	979	c=====================================================================
	980	c --- PRECIPITATING DOWNDRAFT CALCULATION
	981	c=====================================================================
	982	c
	983	c * Check whether ep(inb)=0, if so, skip precipitating *
	984	c * downdraft calculation *
	985	c
	986	c
	987	c * Integrate liquid water equation to find condensed water *
	988	c * and condensed water flux *
	989	c
	990	c
	991	do 890 i=1,ncum
	992	jtt(i)=2
	993	if(ep(i,inb(i)).le.0.0001)iflag(i)=2
	994	if(iflag(i).eq.0)then
	995	lwork(i)=.true.
	996	else
	997	lwork(i)=.false.
	998	endif
	999	890 continue
	1000	c
	1001	c * Begin downdraft loop *
	1002	c
	1003	c
	1004	call zilch(wdtrain,ncum)
	1005	do 899 i=nl+1,1,-1
	1006	c
	1007	num1=0
	1008	do 879 ij=1,ncum
	1009	if((i.le.inb(ij)).and.lwork(ij))num1=num1+1
	1010	879 continue
	1011	if(num1.le.0)go to 899
	1012	c
	1013	c
	1014	c * Calculate detrained precipitation *
	1015	c
	1016	do 891 ij=1,ncum
	1017	if((i.le.inb(ij)).and.(lwork(ij)))then
	1018	wdtrain(ij)=gep(ij,i)m(ij,i)*clw(ij,i)
	1019	endif
	1020	891 continue
	1021	c
	1022	if(i.gt.1)then
	1023	do 893 j=1,i-1
	1024	do 892 ij=1,ncum
	1025	if((i.le.inb(ij)).and.(lwork(ij)))then
	1026	awat=elij(ij,j,i)-(1.-ep(ij,i))*clw(ij,i)
	1027	awat=max(0.0,awat)
	1028	wdtrain(ij)=wdtrain(ij)+gawatment(ij,j,i)
	1029	endif
	1030	892 continue
	1031	893 continue
	1032	endif
	1033	c
	1034	c * Find rain water and evaporation using provisional *
	1035	c * estimates of qp(i)and qp(i-1) *
	1036	c
	1037	c
	1038	c * Value of terminal velocity and coeffecient of evaporation for snow *
	1039	c
	1040	do 894 ij=1,ncum
	1041	if((i.le.inb(ij)).and.(lwork(ij)))then
	1042	coeff=coeffs
	1043	wt(ij,i)=omtsnow
	1044	c
	1045	c * Value of terminal velocity and coeffecient of evaporation for rain *
	1046	c
	1047	if(t(ij,i).gt.273.0)then
	1048	coeff=coeffr
	1049	wt(ij,i)=omtrain
	1050	endif
	1051	qsm=0.5*(q(ij,i)+qp(ij,i+1))
	1052	afac=coeffph(ij,i)(qs(ij,i)-qsm)
	1053	& /(1.0e4+2.0e3ph(ij,i)qs(ij,i))
	1054	afac=max(afac,0.0)
	1055	sigt=sigp(ij,i)
	1056	sigt=max(0.0,sigt)
	1057	sigt=min(1.0,sigt)
	1058	b6=100.(ph(ij,i)-ph(ij,i+1))sigt*afac/wt(ij,i)
	1059	c6=(water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij,i)
	1060	revap=0.5(-b6+sqrt(b6b6+4.*c6))
	1061	evap(ij,i)=sigtafacrevap
	1062	water(ij,i)=revap*revap
	1063	c
	1064	c * Calculate precipitating downdraft mass flux under *
	1065	c * hydrostatic approximation *
	1066	c
	1067	if(i.gt.1)then
	1068	dhdp=(h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
	1069	dhdp=max(dhdp,10.0)
	1070	mp(ij,i)=100.ginvlv(ij,i)sigdevap(ij,i)/dhdp
	1071	mp(ij,i)=max(mp(ij,i),0.0)
	1072	c
	1073	c * Add small amount of inertia to downdraft *
	1074	c
	1075	fac=20.0/(ph(ij,i-1)-ph(ij,i))
	1076	mp(ij,i)=(fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
	1077	c
	1078	c * Force mp to decrease linearly to zero *
	1079	c * between about 950 mb and the surface *
	1080	c
	1081	if(p(ij,i).gt.(0.949*p(ij,1)))then
	1082	jtt(ij)=max(jtt(ij),i)
	1083	mp(ij,i)=mp(ij,jtt(ij))*(p(ij,1)-p(ij,i))
	1084	& /(p(ij,1)-p(ij,jtt(ij)))
	1085	endif
	1086	endif
	1087	c
	1088	c * Find mixing ratio of precipitating downdraft *
	1089	c
	1090	if(i.ne.inb(ij))then
	1091	if(i.eq.1)then
	1092	qstm=qs(ij,1)
	1093	else
	1094	qstm=qs(ij,i-1)
	1095	endif
	1096	if(mp(ij,i).gt.mp(ij,i+1))then
	1097	rat=mp(ij,i+1)/mp(ij,i)
	1098	qp(ij,i)=qp(ij,i+1)rat+q(ij,i)(1.0-rat)+100.ginv
	1099	& sigd(ph(ij,i)-ph(ij,i+1))(evap(ij,i)/mp(ij,i))
	1100	up(ij,i)=up(ij,i+1)rat+u(ij,i)(1.-rat)
	1101	vp(ij,i)=vp(ij,i+1)rat+v(ij,i)(1.-rat)
	1102	else
	1103	if(mp(ij,i+1).gt.0.0)then
	1104	qp(ij,i)=(gz(ij,i+1)-gz(ij,i)
	1105	& +qp(ij,i+1)*(lv(ij,i+1)+t(ij,i+1)
	1106	& (cl-cpd))+cpd(t(ij,i+1)-t(ij,i)))
	1107	& /(lv(ij,i)+t(ij,i)*(cl-cpd))
	1108	up(ij,i)=up(ij,i+1)
	1109	vp(ij,i)=vp(ij,i+1)
	1110	endif
	1111	endif
	1112	qp(ij,i)=min(qp(ij,i),qstm)
	1113	qp(ij,i)=max(qp(ij,i),0.0)
	1114	endif
	1115	endif
	1116	894 continue
	1117	899 continue
	1118	c
	1119	c * Calculate surface precipitation in mm/day *
	1120	c
	1121	do 1190 i=1,ncum
	1122	if(iflag(i).le.1)then
	1123	cc precip(i)=precip(i)+wt(i,1)sigdwater(i,1)3600.24000.
	1124	cc & /(rowl*g)
	1125	cc precip(i)=precip(i)*delt/86400.
	1126	precip(i) = wt(i,1)sigdwater(i,1)*86400/g
	1127	endif
	1128	1190 continue
	1129	c
	1130	c
	1131	c * Calculate downdraft velocity scale and surface temperature and *
	1132	c * water vapor fluctuations *
	1133	c
	1134	c wd=betaabs(mp(icb))0.01rdt(icb)/(sigd*p(icb))
	1135	c qprime=0.5*(qp(1)-q(1))
	1136	c tprime=lv0*qprime/cpd
	1137	c
	1138	c * Calculate tendencies of lowest level potential temperature *
	1139	c * and mixing ratio *
	1140	c
	1141	do 1200 i=1,ncum
	1142	work(i)=0.01/(ph(i,1)-ph(i,2))
	1143	am(i)=0.0
	1144	1200 continue
	1145	do 1220 k=2,nl
	1146	do 1210 i=1,ncum
	1147	if((nk(i).eq.1).and.(k.le.inb(i)).and.(nk(i).eq.1))then
	1148	am(i)=am(i)+m(i,k)
	1149	endif
	1150	1210 continue
	1151	1220 continue
	1152	do 1240 i=1,ncum
	1153	if((gwork(i)am(i)).ge.delti)iflag(i)=1
	1154	ft(i,1)=ft(i,1)+gwork(i)am(i)*(t(i,2)-t(i,1)
	1155	& +(gz(i,2)-gz(i,1))/cpn(i,1))
	1156	ft(i,1)=ft(i,1)-lvcp(i,1)sigdevap(i,1)
	1157	ft(i,1)=ft(i,1)+sigdwt(i,2)(cl-cpd)water(i,2)(t(i,2)
	1158	& -t(i,1))*work(i)/cpn(i,1)
	1159	fq(i,1)=fq(i,1)+gmp(i,2)(qp(i,2)-q(i,1))*
	1160	& work(i)+sigd*evap(i,1)
	1161	fq(i,1)=fq(i,1)+gam(i)(q(i,2)-q(i,1))*work(i)
	1162	fu(i,1)=fu(i,1)+gwork(i)(mp(i,2)*(up(i,2)-u(i,1))
	1163	& +am(i)*(u(i,2)-u(i,1)))
	1164	fv(i,1)=fv(i,1)+gwork(i)(mp(i,2)*(vp(i,2)-v(i,1))
	1165	& +am(i)*(v(i,2)-v(i,1)))
	1166	1240 continue
	1167	do 1260 j=2,nl
	1168	do 1250 i=1,ncum
	1169	if(j.le.inb(i))then
	1170	fq(i,1)=fq(i,1)
	1171	& +gwork(i)ment(i,j,1)*(qent(i,j,1)-q(i,1))
	1172	fu(i,1)=fu(i,1)
	1173	& +gwork(i)ment(i,j,1)*(uent(i,j,1)-u(i,1))
	1174	fv(i,1)=fv(i,1)
	1175	& +gwork(i)ment(i,j,1)*(vent(i,j,1)-v(i,1))
	1176	endif
	1177	1250 continue
	1178	1260 continue
	1179	c
	1180	c * Calculate tendencies of potential temperature and mixing ratio *
	1181	c * at levels above the lowest level *
	1182	c
	1183	c * First find the net saturated updraft and downdraft mass fluxes *
	1184	c * through each level *
	1185	c
	1186	do 1500 i=2,nl+1
	1187	c
	1188	num1=0
	1189	do 1265 ij=1,ncum
	1190	if(i.le.inb(ij))num1=num1+1
	1191	1265 continue
	1192	if(num1.le.0)go to 1500
	1193	c
	1194	call zilch(amp1,ncum)
	1195	call zilch(ad,ncum)
	1196	c
	1197	do 1280 k=i+1,nl+1
	1198	do 1270 ij=1,ncum
	1199	if((i.ge.nk(ij)).and.(i.le.inb(ij))
	1200	& .and.(k.le.(inb(ij)+1)))then
	1201	amp1(ij)=amp1(ij)+m(ij,k)
	1202	endif
	1203	1270 continue
	1204	1280 continue
	1205	c
	1206	do 1310 k=1,i
	1207	do 1300 j=i+1,nl+1
	1208	do 1290 ij=1,ncum
	1209	if((j.le.(inb(ij)+1)).and.(i.le.inb(ij)))then
	1210	amp1(ij)=amp1(ij)+ment(ij,k,j)
	1211	endif
	1212	1290 continue
	1213	1300 continue
	1214	1310 continue
	1215	do 1340 k=1,i-1
	1216	do 1330 j=i,nl+1
	1217	do 1320 ij=1,ncum
	1218	if((i.le.inb(ij)).and.(j.le.inb(ij)))then
	1219	ad(ij)=ad(ij)+ment(ij,j,k)
	1220	endif
	1221	1320 continue
	1222	1330 continue
	1223	1340 continue
	1224	c
	1225	do 1350 ij=1,ncum
	1226	if(i.le.inb(ij))then
	1227	dpinv=0.01/(ph(ij,i)-ph(ij,i+1))
	1228	cpinv=1.0/cpn(ij,i)
	1229	c
	1230	ft(ij,i)=ft(ij,i)
	1231	& +gdpinv(amp1(ij)*(t(ij,i+1)-t(ij,i)
	1232	& +(gz(ij,i+1)-gz(ij,i))*cpinv)
	1233	& -ad(ij)(t(ij,i)-t(ij,i-1)+(gz(ij,i)-gz(ij,i-1))cpinv))
	1234	& -sigdlvcp(ij,i)evap(ij,i)
	1235	ft(ij,i)=ft(ij,i)+gdpinvment(ij,i,i)*(hp(ij,i)-h(ij,i)+
	1236	& t(ij,i)(cpv-cpd)(q(ij,i)-qent(ij,i,i)))*cpinv
	1237	ft(ij,i)=ft(ij,i)+sigdwt(ij,i+1)(cl-cpd)water(ij,i+1)
	1238	& (t(ij,i+1)-t(ij,i))dpinvcpinv
	1239	fq(ij,i)=fq(ij,i)+gdpinv(amp1(ij)*(q(ij,i+1)-q(ij,i))-
	1240	& ad(ij)*(q(ij,i)-q(ij,i-1)))
	1241	fu(ij,i)=fu(ij,i)+gdpinv(amp1(ij)*(u(ij,i+1)-u(ij,i))-
	1242	& ad(ij)*(u(ij,i)-u(ij,i-1)))
	1243	fv(ij,i)=fv(ij,i)+gdpinv(amp1(ij)*(v(ij,i+1)-v(ij,i))-
	1244	& ad(ij)*(v(ij,i)-v(ij,i-1)))
	1245	endif
	1246	1350 continue
	1247	do 1370 k=1,i-1
	1248	do 1360 ij=1,ncum
	1249	if(i.le.inb(ij))then
	1250	awat=elij(ij,k,i)-(1.-ep(ij,i))*clw(ij,i)
	1251	awat=max(awat,0.0)
	1252	fq(ij,i)=fq(ij,i)
	1253	& +gdpinvment(ij,k,i)*(qent(ij,k,i)-awat-q(ij,i))
	1254	fu(ij,i)=fu(ij,i)
	1255	& +gdpinvment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
	1256	fv(ij,i)=fv(ij,i)
	1257	& +gdpinvment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
	1258	endif
	1259	1360 continue
	1260	1370 continue
	1261	do 1390 k=i,nl+1
	1262	do 1380 ij=1,ncum
	1263	if((i.le.inb(ij)).and.(k.le.inb(ij)))then
	1264	fq(ij,i)=fq(ij,i)
	1265	& +gdpinvment(ij,k,i)*(qent(ij,k,i)-q(ij,i))
	1266	fu(ij,i)=fu(ij,i)
	1267	& +gdpinvment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
	1268	fv(ij,i)=fv(ij,i)
	1269	& +gdpinvment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
	1270	endif
	1271	1380 continue
	1272	1390 continue
	1273	do 1400 ij=1,ncum
	1274	if(i.le.inb(ij))then
	1275	fq(ij,i)=fq(ij,i)
	1276	& +sigdevap(ij,i)+g(mp(ij,i+1)*
	1277	& (qp(ij,i+1)-q(ij,i))
	1278	& -mp(ij,i)(qp(ij,i)-q(ij,i-1)))dpinv
	1279	fu(ij,i)=fu(ij,i)
	1280	& +g(mp(ij,i+1)(up(ij,i+1)-u(ij,i))-mp(ij,i)*
	1281	& (up(ij,i)-u(ij,i-1)))*dpinv
	1282	fv(ij,i)=fv(ij,i)
	1283	& +g(mp(ij,i+1)(vp(ij,i+1)-v(ij,i))-mp(ij,i)*
	1284	& (vp(ij,i)-v(ij,i-1)))*dpinv
	1285	endif
	1286	1400 continue
	1287	1500 continue
	1288	c
	1289	c * Adjust tendencies at top of convection layer to reflect *
	1290	c * actual position of the level zero cape *
	1291	c
	1292	do 503 ij=1,ncum
	1293	fqold=fq(ij,inb(ij))
	1294	fq(ij,inb(ij))=fq(ij,inb(ij))*(1.-frac(ij))
	1295	fq(ij,inb(ij)-1)=fq(ij,inb(ij)-1)
	1296	& +frac(ij)fqold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
	1297	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*lv(ij,inb(ij))
	1298	& /lv(ij,inb(ij)-1)
	1299	ftold=ft(ij,inb(ij))
	1300	ft(ij,inb(ij))=ft(ij,inb(ij))*(1.-frac(ij))
	1301	ft(ij,inb(ij)-1)=ft(ij,inb(ij)-1)
	1302	& +frac(ij)ftold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
	1303	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*cpn(ij,inb(ij))
	1304	& /cpn(ij,inb(ij)-1)
	1305	fuold=fu(ij,inb(ij))
	1306	fu(ij,inb(ij))=fu(ij,inb(ij))*(1.-frac(ij))
	1307	fu(ij,inb(ij)-1)=fu(ij,inb(ij)-1)
	1308	& +frac(ij)fuold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
	1309	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
	1310	fvold=fv(ij,inb(ij))
	1311	fv(ij,inb(ij))=fv(ij,inb(ij))*(1.-frac(ij))
	1312	fv(ij,inb(ij)-1)=fv(ij,inb(ij)-1)
	1313	& +frac(ij)fvold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
	1314	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
	1315	503 continue
	1316	c
	1317	c * Very slightly adjust tendencies to force exact *
	1318	c * enthalpy, momentum and tracer conservation *
	1319	c
	1320	do 682 ij=1,ncum
	1321	ents(ij)=0.0
	1322	uav(ij)=0.0
	1323	vav(ij)=0.0
	1324	do 681 i=1,inb(ij)
	1325	ents(ij)=ents(ij)
	1326	& +(cpn(ij,i)ft(ij,i)+lv(ij,i)fq(ij,i))*(ph(ij,i)-ph(ij,i+1))
	1327	uav(ij)=uav(ij)+fu(ij,i)*(ph(ij,i)-ph(ij,i+1))
	1328	vav(ij)=vav(ij)+fv(ij,i)*(ph(ij,i)-ph(ij,i+1))
	1329	681 continue
	1330	682 continue
	1331	do 683 ij=1,ncum
	1332	ents(ij)=ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
	1333	uav(ij)=uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
	1334	vav(ij)=vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
	1335	683 continue
	1336	do 642 ij=1,ncum
	1337	do 641 i=1,inb(ij)
	1338	ft(ij,i)=ft(ij,i)-ents(ij)/cpn(ij,i)
	1339	fu(ij,i)=(1.-cu)*(fu(ij,i)-uav(ij))
	1340	fv(ij,i)=(1.-cu)*(fv(ij,i)-vav(ij))
	1341	641 continue
	1342	642 continue
	1343	c
	1344	do 1810 k=1,nl+1
	1345	do 1800 i=1,ncum
	1346	if((q(i,k)+delt*fq(i,k)).lt.0.0)iflag(i)=10
	1347	1800 continue
	1348	1810 continue
	1349	c
	1350	c
	1351	do 1900 i=1,ncum
	1352	if(iflag(i).gt.2)then
	1353	precip(i)=0.0
	1354	cbmf(i)=0.0
	1355	endif
	1356	1900 continue
	1357	do 1920 k=1,nl
	1358	do 1910 i=1,ncum
	1359	if(iflag(i).gt.2)then
	1360	ft(i,k)=0.0
	1361	fq(i,k)=0.0
	1362	fu(i,k)=0.0
	1363	fv(i,k)=0.0
	1364	endif
	1365	1910 continue
	1366	1920 continue
	1367	do 2000 i=1,ncum
	1368	precip1(idcum(i))=precip(i)
	1369	cbmf1(idcum(i))=cbmf(i)
	1370	iflag1(idcum(i))=iflag(i)
	1371	2000 continue
	1372	do 2020 k=1,nl
	1373	do 2010 i=1,ncum
	1374	ft1(idcum(i),k)=ft(i,k)
	1375	fq1(idcum(i),k)=fq(i,k)
	1376	fu1(idcum(i),k)=fu(i,k)
	1377	fv1(idcum(i),k)=fv(i,k)
	1378	2010 continue
	1379	2020 continue
	1380	c
	1381	DO k=1,nd
	1382	DO i=1,len
	1383	Ma(i,k) = 0.
	1384	ENDDO
	1385	ENDDO
	1386	DO k=nl,1,-1
	1387	DO i=1,ncum
	1388	Ma(i,k) = Ma(i,k+1)+m(i,k)
	1389	ENDDO
	1390	ENDDO
	1391	c
	1392	return
	1393	end
	1394

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

Download in other formats:

Original Format