Context Navigation

convect2.F @ 3948

Last change on this file since 3948 was 268, checked in by lmdzadmin, 23 years ago
Initial version in branch LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 44.9 KB

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

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

Download in other formats:

Original Format