Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

cv_routines.F @ 5306

Last change on this file since 5306 was 416, checked in by lmdzadmin, 22 years ago
Inclusion initiale
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 55.0 KB

Rev	Line
[416]	1	SUBROUTINE cv_param(nd)
	2	implicit none
	3
	4	c------------------------------------------------------------
	5	c Set parameters for convectL
	6	c (includes microphysical parameters and parameters that
	7	c control the rate of approach to quasi-equilibrium)
	8	c------------------------------------------------------------
	9
	10	C * ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) *
	11	C * TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- *
	12	C * CONVERSION THRESHOLD IS ASSUMED TO BE ZERO *
	13	C * (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY *
	14	C * BETWEEN 0 C AND TLCRIT) *
	15	C * ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT *
	16	C * FORMULATION *
	17	C * SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *
	18	C * SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *
	19	C * OF CLOUD *
	20	C * OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN *
	21	C * OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW *
	22	C * COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *
	23	C * OF RAIN *
	24	C * COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *
	25	C * OF SNOW *
	26	C * CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM *
	27	C * TRANSPORT *
	28	C * DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION *
	29	C * A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC *
	30	C * ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF *
	31	C * APPROACH TO QUASI-EQUILIBRIUM *
	32	C * (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) *
	33	C * (DAMP MUST BE LESS THAN 1) *
	34
	35	#include "cvparam.h"
	36	integer nd
	37
	38	c noff: integer limit for convection (nd-noff)
	39	c minorig: First level of convection
	40
	41	noff = 2
	42	minorig = 2
	43
	44	nl=nd-noff
	45	nlp=nl+1
	46	nlm=nl-1
	47
	48	elcrit=0.0011
	49	tlcrit=-55.0
	50	entp=1.5
	51	sigs=0.12
	52	sigd=0.05
	53	omtrain=50.0
	54	omtsnow=5.5
	55	coeffr=1.0
	56	coeffs=0.8
	57	dtmax=0.9
	58	c
	59	cu=0.70
	60	c
	61	betad=10.0
	62	c
	63	damp=0.1
	64	alpha=0.2
	65	c
	66	delta=0.01 ! cld
	67	c
	68	return
	69	end
	70
	71	SUBROUTINE cv_prelim(len,nd,ndp1,t,q,p,ph
	72	: ,lv,cpn,tv,gz,h,hm)
	73	implicit none
	74
	75	!=====================================================================
	76	! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
	77	!=====================================================================
	78
	79	c inputs:
	80	integer len, nd, ndp1
	81	real t(len,nd), q(len,nd), p(len,nd), ph(len,ndp1)
	82
	83	c outputs:
	84	real lv(len,nd), cpn(len,nd), tv(len,nd)
	85	real gz(len,nd), h(len,nd), hm(len,nd)
	86
	87	c local variables:
	88	integer k, i
	89	real cpx(len,nd)
	90
	91	#include "cvthermo.h"
	92	#include "cvparam.h"
	93
	94
	95	do 110 k=1,nlp
	96	do 100 i=1,len
	97	lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
	98	cpn(i,k)=cpd(1.0-q(i,k))+cpvq(i,k)
	99	cpx(i,k)=cpd(1.0-q(i,k))+clq(i,k)
	100	tv(i,k)=t(i,k)(1.0+q(i,k)epsim1)
	101	100 continue
	102	110 continue
	103	c
	104	c gz = phi at the full levels (same as p).
	105	c
	106	do 120 i=1,len
	107	gz(i,1)=0.0
	108	120 continue
	109	do 140 k=2,nlp
	110	do 130 i=1,len
	111	gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
	112	& *(p(i,k-1)-p(i,k))/ph(i,k)
	113	130 continue
	114	140 continue
	115	c
	116	c h = phi + cpT (dry static energy).
	117	c hm = phi + cp(T-Tbase)+Lq
	118	c
	119	do 170 k=1,nlp
	120	do 160 i=1,len
	121	h(i,k)=gz(i,k)+cpn(i,k)*t(i,k)
	122	hm(i,k)=gz(i,k)+cpx(i,k)(t(i,k)-t(i,1))+lv(i,k)q(i,k)
	123	160 continue
	124	170 continue
	125
	126	return
	127	end
	128
	129	SUBROUTINE cv_feed(len,nd,t,q,qs,p,hm,gz
	130	: ,nk,icb,icbmax,iflag,tnk,qnk,gznk,plcl)
	131	implicit none
	132
	133	C================================================================
	134	C Purpose: CONVECTIVE FEED
	135	C================================================================
	136
	137	#include "cvparam.h"
	138
	139	c inputs:
	140	integer len, nd
	141	real t(len,nd), q(len,nd), qs(len,nd), p(len,nd)
	142	real hm(len,nd), gz(len,nd)
	143
	144	c outputs:
	145	integer iflag(len), nk(len), icb(len), icbmax
	146	real tnk(len), qnk(len), gznk(len), plcl(len)
	147
	148	c local variables:
	149	integer i, k
	150	integer ihmin(len)
	151	real work(len)
	152	real pnk(len), qsnk(len), rh(len), chi(len)
	153
	154	!-------------------------------------------------------------------
	155	! --- Find level of minimum moist static energy
	156	! --- If level of minimum moist static energy coincides with
	157	! --- or is lower than minimum allowable parcel origin level,
	158	! --- set iflag to 6.
	159	!-------------------------------------------------------------------
	160
	161	do 180 i=1,len
	162	work(i)=1.0e12
	163	ihmin(i)=nl
	164	180 continue
	165	do 200 k=2,nlp
	166	do 190 i=1,len
	167	if((hm(i,k).lt.work(i)).and.
	168	& (hm(i,k).lt.hm(i,k-1)))then
	169	work(i)=hm(i,k)
	170	ihmin(i)=k
	171	endif
	172	190 continue
	173	200 continue
	174	do 210 i=1,len
	175	ihmin(i)=min(ihmin(i),nlm)
	176	if(ihmin(i).le.minorig)then
	177	iflag(i)=6
	178	endif
	179	210 continue
	180	c
	181	!-------------------------------------------------------------------
	182	! --- Find that model level below the level of minimum moist static
	183	! --- energy that has the maximum value of moist static energy
	184	!-------------------------------------------------------------------
	185
	186	do 220 i=1,len
	187	work(i)=hm(i,minorig)
	188	nk(i)=minorig
	189	220 continue
	190	do 240 k=minorig+1,nl
	191	do 230 i=1,len
	192	if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then
	193	work(i)=hm(i,k)
	194	nk(i)=k
	195	endif
	196	230 continue
	197	240 continue
	198	!-------------------------------------------------------------------
	199	! --- Check whether parcel level temperature and specific humidity
	200	! --- are reasonable
	201	!-------------------------------------------------------------------
	202	do 250 i=1,len
	203	if(((t(i,nk(i)).lt.250.0).or.
	204	& (q(i,nk(i)).le.0.0).or.
	205	& (p(i,ihmin(i)).lt.400.0)).and.
	206	& (iflag(i).eq.0))iflag(i)=7
	207	250 continue
	208	!-------------------------------------------------------------------
	209	! --- Calculate lifted condensation level of air at parcel origin level
	210	! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
	211	!-------------------------------------------------------------------
	212	do 260 i=1,len
	213	tnk(i)=t(i,nk(i))
	214	qnk(i)=q(i,nk(i))
	215	gznk(i)=gz(i,nk(i))
	216	pnk(i)=p(i,nk(i))
	217	qsnk(i)=qs(i,nk(i))
	218	c
	219	rh(i)=qnk(i)/qsnk(i)
	220	rh(i)=min(1.0,rh(i))
	221	chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
	222	plcl(i)=pnk(i)(rh(i)*chi(i))
	223	if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0))
	224	& .and.(iflag(i).eq.0))iflag(i)=8
	225	260 continue
	226	!-------------------------------------------------------------------
	227	! --- Calculate first level above lcl (=icb)
	228	!-------------------------------------------------------------------
	229	do 270 i=1,len
	230	icb(i)=nlm
	231	270 continue
	232	c
	233	do 290 k=minorig,nl
	234	do 280 i=1,len
	235	if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
	236	& icb(i)=min(icb(i),k)
	237	280 continue
	238	290 continue
	239	c
	240	do 300 i=1,len
	241	if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
	242	300 continue
	243	c
	244	c Compute icbmax.
	245	c
	246	icbmax=2
	247	do 310 i=1,len
	248	icbmax=max(icbmax,icb(i))
	249	310 continue
	250
	251	return
	252	end
	253
	254	SUBROUTINE cv_undilute1(len,nd,t,q,qs,gz,p,nk,icb,icbmax
	255	: ,tp,tvp,clw)
	256	implicit none
	257
	258	#include "cvthermo.h"
	259	#include "cvparam.h"
	260
	261	c inputs:
	262	integer len, nd
	263	integer nk(len), icb(len), icbmax
	264	real t(len,nd), q(len,nd), qs(len,nd), gz(len,nd)
	265	real p(len,nd)
	266
	267	c outputs:
	268	real tp(len,nd), tvp(len,nd), clw(len,nd)
	269
	270	c local variables:
	271	integer i, k
	272	real tg, qg, alv, s, ahg, tc, denom, es, rg
	273	real ah0(len), cpp(len)
	274	real tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
	275
	276	!-------------------------------------------------------------------
	277	! --- Calculates the lifted parcel virtual temperature at nk,
	278	! --- the actual temperature, and the adiabatic
	279	! --- liquid water content. The procedure is to solve the equation.
	280	! cptp+Lqp+phi=cptnk+Lqnk+gznk.
	281	!-------------------------------------------------------------------
	282
	283	do 320 i=1,len
	284	tnk(i)=t(i,nk(i))
	285	qnk(i)=q(i,nk(i))
	286	gznk(i)=gz(i,nk(i))
	287	ticb(i)=t(i,icb(i))
	288	gzicb(i)=gz(i,icb(i))
	289	320 continue
	290	c
	291	c * Calculate certain parcel quantities, including static energy *
	292	c
	293	do 330 i=1,len
	294	ah0(i)=(cpd(1.-qnk(i))+clqnk(i))*tnk(i)
	295	& +qnk(i)(lv0-clmcpv(tnk(i)-273.15))+gznk(i)
	296	cpp(i)=cpd(1.-qnk(i))+qnk(i)cpv
	297	330 continue
	298	c
	299	c * Calculate lifted parcel quantities below cloud base *
	300	c
	301	do 350 k=minorig,icbmax-1
	302	do 340 i=1,len
	303	tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))/cpp(i)
	304	tvp(i,k)=tp(i,k)(1.+qnk(i)epsi)
	305	340 continue
	306	350 continue
	307	c
	308	c * Find lifted parcel quantities above cloud base *
	309	c
	310	do 360 i=1,len
	311	tg=ticb(i)
	312	qg=qs(i,icb(i))
	313	alv=lv0-clmcpv*(ticb(i)-t0)
	314	c
	315	c First iteration.
	316	c
	317	s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	318	s=1./s
	319	ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	320	tg=tg+s*(ah0(i)-ahg)
	321	tg=max(tg,35.0)
	322	tc=tg-t0
	323	denom=243.5+tc
	324	if(tc.ge.0.0)then
	325	es=6.112exp(17.67tc/denom)
	326	else
	327	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	328	endif
	329	qg=epses/(p(i,icb(i))-es(1.-eps))
	330	c
	331	c Second iteration.
	332	c
	333	s=cpd+alvalvqg/(rrvticb(i)ticb(i))
	334	s=1./s
	335	ahg=cpdtg+(cl-cpd)qnk(i)ticb(i)+alvqg+gzicb(i)
	336	tg=tg+s*(ah0(i)-ahg)
	337	tg=max(tg,35.0)
	338	tc=tg-t0
	339	denom=243.5+tc
	340	if(tc.ge.0.0)then
	341	es=6.112exp(17.67tc/denom)
	342	else
	343	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	344	end if
	345	qg=epses/(p(i,icb(i))-es(1.-eps))
	346	c
	347	alv=lv0-clmcpv*(ticb(i)-273.15)
	348	tp(i,icb(i))=(ah0(i)-(cl-cpd)qnk(i)ticb(i)
	349	& -gz(i,icb(i))-alv*qg)/cpd
	350	clw(i,icb(i))=qnk(i)-qg
	351	clw(i,icb(i))=max(0.0,clw(i,icb(i)))
	352	rg=qg/(1.-qnk(i))
	353	tvp(i,icb(i))=tp(i,icb(i))(1.+rgepsi)
	354	360 continue
	355	c
	356	do 380 k=minorig,icbmax
	357	do 370 i=1,len
	358	tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
	359	370 continue
	360	380 continue
	361	c
	362	return
	363	end
	364
	365	SUBROUTINE cv_trigger(len,nd,icb,cbmf,tv,tvp,iflag)
	366	implicit none
	367
	368	!-------------------------------------------------------------------
	369	! --- Test for instability.
	370	! --- If there was no convection at last time step and parcel
	371	! --- is stable at icb, then set iflag to 4.
	372	!-------------------------------------------------------------------
	373
	374	#include "cvparam.h"
	375
	376	c inputs:
	377	integer len, nd, icb(len)
	378	real cbmf(len), tv(len,nd), tvp(len,nd)
	379
	380	c outputs:
	381	integer iflag(len) ! also an input
	382
	383	c local variables:
	384	integer i
	385
	386
	387	do 390 i=1,len
	388	if((cbmf(i).eq.0.0) .and.(iflag(i).eq.0).and.
	389	& (tvp(i,icb(i)).le.(tv(i,icb(i))-dtmax)))iflag(i)=4
	390	390 continue
	391
	392	return
	393	end
	394
	395	SUBROUTINE cv_compress( len,nloc,ncum,nd
	396	: ,iflag1,nk1,icb1
	397	: ,cbmf1,plcl1,tnk1,qnk1,gznk1
	398	: ,t1,q1,qs1,u1,v1,gz1
	399	: ,h1,lv1,cpn1,p1,ph1,tv1,tp1,tvp1,clw1
	400	o ,iflag,nk,icb
	401	o ,cbmf,plcl,tnk,qnk,gznk
	402	o ,t,q,qs,u,v,gz,h,lv,cpn,p,ph,tv,tp,tvp,clw
	403	o ,dph )
	404	implicit none
	405
	406	#include "cvparam.h"
	407
	408	c inputs:
	409	integer len,ncum,nd,nloc
	410	integer iflag1(len),nk1(len),icb1(len)
	411	real cbmf1(len),plcl1(len),tnk1(len),qnk1(len),gznk1(len)
	412	real t1(len,nd),q1(len,nd),qs1(len,nd),u1(len,nd),v1(len,nd)
	413	real gz1(len,nd),h1(len,nd),lv1(len,nd),cpn1(len,nd)
	414	real p1(len,nd),ph1(len,nd+1),tv1(len,nd),tp1(len,nd)
	415	real tvp1(len,nd),clw1(len,nd)
	416
	417	c outputs:
	418	integer iflag(nloc),nk(nloc),icb(nloc)
	419	real cbmf(nloc),plcl(nloc),tnk(nloc),qnk(nloc),gznk(nloc)
	420	real t(nloc,nd),q(nloc,nd),qs(nloc,nd),u(nloc,nd),v(nloc,nd)
	421	real gz(nloc,nd),h(nloc,nd),lv(nloc,nd),cpn(nloc,nd)
	422	real p(nloc,nd),ph(nloc,nd+1),tv(nloc,nd),tp(nloc,nd)
	423	real tvp(nloc,nd),clw(nloc,nd)
	424	real dph(nloc,nd)
	425
	426	c local variables:
	427	integer i,k,nn
	428
	429
	430	do 110 k=1,nl+1
	431	nn=0
	432	do 100 i=1,len
	433	if(iflag1(i).eq.0)then
	434	nn=nn+1
	435	t(nn,k)=t1(i,k)
	436	q(nn,k)=q1(i,k)
	437	qs(nn,k)=qs1(i,k)
	438	u(nn,k)=u1(i,k)
	439	v(nn,k)=v1(i,k)
	440	gz(nn,k)=gz1(i,k)
	441	h(nn,k)=h1(i,k)
	442	lv(nn,k)=lv1(i,k)
	443	cpn(nn,k)=cpn1(i,k)
	444	p(nn,k)=p1(i,k)
	445	ph(nn,k)=ph1(i,k)
	446	tv(nn,k)=tv1(i,k)
	447	tp(nn,k)=tp1(i,k)
	448	tvp(nn,k)=tvp1(i,k)
	449	clw(nn,k)=clw1(i,k)
	450	endif
	451	100 continue
	452	110 continue
	453
	454	if (nn.ne.ncum) then
	455	print*,'strange! nn not equal to ncum: ',nn,ncum
	456	stop
	457	endif
	458
	459	nn=0
	460	do 150 i=1,len
	461	if(iflag1(i).eq.0)then
	462	nn=nn+1
	463	cbmf(nn)=cbmf1(i)
	464	plcl(nn)=plcl1(i)
	465	tnk(nn)=tnk1(i)
	466	qnk(nn)=qnk1(i)
	467	gznk(nn)=gznk1(i)
	468	nk(nn)=nk1(i)
	469	icb(nn)=icb1(i)
	470	iflag(nn)=iflag1(i)
	471	endif
	472	150 continue
	473
	474	do 170 k=1,nl
	475	do 160 i=1,ncum
	476	dph(i,k)=ph(i,k)-ph(i,k+1)
	477	160 continue
	478	170 continue
	479
	480	return
	481	end
	482
	483	SUBROUTINE cv_undilute2(nloc,ncum,nd,icb,nk
	484	: ,tnk,qnk,gznk,t,q,qs,gz
	485	: ,p,dph,h,tv,lv
	486	o ,inb,inb1,tp,tvp,clw,hp,ep,sigp,frac)
	487	implicit none
	488
	489	C---------------------------------------------------------------------
	490	C Purpose:
	491	C FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	492	C &
	493	C COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
	494	C FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
	495	C &
	496	C FIND THE LEVEL OF NEUTRAL BUOYANCY
	497	C---------------------------------------------------------------------
	498
	499	#include "cvthermo.h"
	500	#include "cvparam.h"
	501
	502	c inputs:
	503	integer ncum, nd, nloc
	504	integer icb(nloc), nk(nloc)
	505	real t(nloc,nd), q(nloc,nd), qs(nloc,nd), gz(nloc,nd)
	506	real p(nloc,nd), dph(nloc,nd)
	507	real tnk(nloc), qnk(nloc), gznk(nloc)
	508	real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)
	509
	510	c outputs:
	511	integer inb(nloc), inb1(nloc)
	512	real tp(nloc,nd), tvp(nloc,nd), clw(nloc,nd)
	513	real ep(nloc,nd), sigp(nloc,nd), hp(nloc,nd)
	514	real frac(nloc)
	515
	516	c local variables:
	517	integer i, k
	518	real tg,qg,ahg,alv,s,tc,es,denom,rg,tca,elacrit
	519	real by, defrac
	520	real ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
	521	logical lcape(nloc)
	522
	523	!=====================================================================
	524	! --- SOME INITIALIZATIONS
	525	!=====================================================================
	526
	527	do 170 k=1,nl
	528	do 160 i=1,ncum
	529	ep(i,k)=0.0
	530	sigp(i,k)=sigs
	531	160 continue
	532	170 continue
	533
	534	!=====================================================================
	535	! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
	536	!=====================================================================
	537	c
	538	c --- The procedure is to solve the equation.
	539	c cptp+Lqp+phi=cptnk+Lqnk+gznk.
	540	c
	541	c * Calculate certain parcel quantities, including static energy *
	542	c
	543	c
	544	do 240 i=1,ncum
	545	ah0(i)=(cpd(1.-qnk(i))+clqnk(i))*tnk(i)
	546	& +qnk(i)(lv0-clmcpv(tnk(i)-t0))+gznk(i)
	547	240 continue
	548	c
	549	c
	550	c * Find lifted parcel quantities above cloud base *
	551	c
	552	c
	553	do 300 k=minorig+1,nl
	554	do 290 i=1,ncum
	555	if(k.ge.(icb(i)+1))then
	556	tg=t(i,k)
	557	qg=qs(i,k)
	558	alv=lv0-clmcpv*(t(i,k)-t0)
	559	c
	560	c First iteration.
	561	c
	562	s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
	563	s=1./s
	564	ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
	565	tg=tg+s*(ah0(i)-ahg)
	566	tg=max(tg,35.0)
	567	tc=tg-t0
	568	denom=243.5+tc
	569	if(tc.ge.0.0)then
	570	es=6.112exp(17.67tc/denom)
	571	else
	572	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	573	endif
	574	qg=epses/(p(i,k)-es(1.-eps))
	575	c
	576	c Second iteration.
	577	c
	578	s=cpd+alvalvqg/(rrvt(i,k)t(i,k))
	579	s=1./s
	580	ahg=cpdtg+(cl-cpd)qnk(i)t(i,k)+alvqg+gz(i,k)
	581	tg=tg+s*(ah0(i)-ahg)
	582	tg=max(tg,35.0)
	583	tc=tg-t0
	584	denom=243.5+tc
	585	if(tc.ge.0.0)then
	586	es=6.112exp(17.67tc/denom)
	587	else
	588	es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
	589	endif
	590	qg=epses/(p(i,k)-es(1.-eps))
	591	c
	592	alv=lv0-clmcpv*(t(i,k)-t0)
	593	c print*,'cpd dans convect2 ',cpd
	594	c print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
	595	c print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
	596	tp(i,k)=(ah0(i)-(cl-cpd)qnk(i)t(i,k)-gz(i,k)-alv*qg)/cpd
	597	c if (.not.cpd.gt.1000.) then
	598	c print*,'CPD=',cpd
	599	c stop
	600	c endif
	601	clw(i,k)=qnk(i)-qg
	602	clw(i,k)=max(0.0,clw(i,k))
	603	rg=qg/(1.-qnk(i))
	604	tvp(i,k)=tp(i,k)(1.+rgepsi)
	605	endif
	606	290 continue
	607	300 continue
	608	c
	609	!=====================================================================
	610	! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
	611	! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
	612	! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
	613	!=====================================================================
	614	c
	615	do 320 k=minorig+1,nl
	616	do 310 i=1,ncum
	617	if(k.ge.(nk(i)+1))then
	618	tca=tp(i,k)-t0
	619	if(tca.ge.0.0)then
	620	elacrit=elcrit
	621	else
	622	elacrit=elcrit*(1.0-tca/tlcrit)
	623	endif
	624	elacrit=max(elacrit,0.0)
	625	ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8)
	626	ep(i,k)=max(ep(i,k),0.0 )
	627	ep(i,k)=min(ep(i,k),1.0 )
	628	sigp(i,k)=sigs
	629	endif
	630	310 continue
	631	320 continue
	632	c
	633	!=====================================================================
	634	! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
	635	! --- VIRTUAL TEMPERATURE
	636	!=====================================================================
	637	c
	638	do 340 k=minorig+1,nl
	639	do 330 i=1,ncum
	640	if(k.ge.(icb(i)+1))then
	641	tvp(i,k)=tvp(i,k)(1.0-qnk(i)+ep(i,k)clw(i,k))
	642	c print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
	643	c print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
	644	endif
	645	330 continue
	646	340 continue
	647	do 350 i=1,ncum
	648	tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
	649	350 continue
	650	c
	651	c=====================================================================
	652	c --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
	653	c --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
	654	c --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
	655	c=====================================================================
	656	c
	657	do 510 i=1,ncum
	658	cape(i)=0.0
	659	capem(i)=0.0
	660	inb(i)=icb(i)+1
	661	inb1(i)=inb(i)
	662	510 continue
	663	c
	664	c Originial Code
	665	c
	666	c do 530 k=minorig+1,nl-1
	667	c do 520 i=1,ncum
	668	c if(k.ge.(icb(i)+1))then
	669	c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	670	c byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	671	c cape(i)=cape(i)+by
	672	c if(by.ge.0.0)inb1(i)=k+1
	673	c if(cape(i).gt.0.0)then
	674	c inb(i)=k+1
	675	c capem(i)=cape(i)
	676	c endif
	677	c endif
	678	c520 continue
	679	c530 continue
	680	c do 540 i=1,ncum
	681	c byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
	682	c cape(i)=capem(i)+byp
	683	c defrac=capem(i)-cape(i)
	684	c defrac=max(defrac,0.001)
	685	c frac(i)=-cape(i)/defrac
	686	c frac(i)=min(frac(i),1.0)
	687	c frac(i)=max(frac(i),0.0)
	688	c540 continue
	689	c
	690	c K Emanuel fix
	691	c
	692	c call zilch(byp,ncum)
	693	c do 530 k=minorig+1,nl-1
	694	c do 520 i=1,ncum
	695	c if(k.ge.(icb(i)+1))then
	696	c by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	697	c cape(i)=cape(i)+by
	698	c if(by.ge.0.0)inb1(i)=k+1
	699	c if(cape(i).gt.0.0)then
	700	c inb(i)=k+1
	701	c capem(i)=cape(i)
	702	c byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	703	c endif
	704	c endif
	705	c520 continue
	706	c530 continue
	707	c do 540 i=1,ncum
	708	c inb(i)=max(inb(i),inb1(i))
	709	c cape(i)=capem(i)+byp(i)
	710	c defrac=capem(i)-cape(i)
	711	c defrac=max(defrac,0.001)
	712	c frac(i)=-cape(i)/defrac
	713	c frac(i)=min(frac(i),1.0)
	714	c frac(i)=max(frac(i),0.0)
	715	c540 continue
	716	c
	717	c J Teixeira fix
	718	c
	719	call zilch(byp,ncum)
	720	do 515 i=1,ncum
	721	lcape(i)=.true.
	722	515 continue
	723	do 530 k=minorig+1,nl-1
	724	do 520 i=1,ncum
	725	if(cape(i).lt.0.0)lcape(i)=.false.
	726	if((k.ge.(icb(i)+1)).and.lcape(i))then
	727	by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
	728	byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
	729	cape(i)=cape(i)+by
	730	if(by.ge.0.0)inb1(i)=k+1
	731	if(cape(i).gt.0.0)then
	732	inb(i)=k+1
	733	capem(i)=cape(i)
	734	endif
	735	endif
	736	520 continue
	737	530 continue
	738	do 540 i=1,ncum
	739	cape(i)=capem(i)+byp(i)
	740	defrac=capem(i)-cape(i)
	741	defrac=max(defrac,0.001)
	742	frac(i)=-cape(i)/defrac
	743	frac(i)=min(frac(i),1.0)
	744	frac(i)=max(frac(i),0.0)
	745	540 continue
	746	c
	747	c=====================================================================
	748	c --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
	749	c=====================================================================
	750	c
	751	c initialization:
	752	do i=1,ncum*nlp
	753	hp(i,1)=h(i,1)
	754	enddo
	755
	756	do 600 k=minorig+1,nl
	757	do 590 i=1,ncum
	758	if((k.ge.icb(i)).and.(k.le.inb(i)))then
	759	hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)t(i,k))ep(i,k)*clw(i,k)
	760	endif
	761	590 continue
	762	600 continue
	763	c
	764	return
	765	end
	766	c
	767	SUBROUTINE cv_closure(nloc,ncum,nd,nk,icb
	768	: ,tv,tvp,p,ph,dph,plcl,cpn
	769	: ,iflag,cbmf)
	770	implicit none
	771
	772	c inputs:
	773	integer ncum, nd, nloc
	774	integer nk(nloc), icb(nloc)
	775	real tv(nloc,nd), tvp(nloc,nd), p(nloc,nd), dph(nloc,nd)
	776	real ph(nloc,nd+1) ! caution nd instead ndp1 to be consistent...
	777	real plcl(nloc), cpn(nloc,nd)
	778
	779	c outputs:
	780	integer iflag(nloc)
	781	real cbmf(nloc) ! also an input
	782
	783	c local variables:
	784	integer i, k, icbmax
	785	real dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
	786	real work(nloc)
	787
	788	#include "cvthermo.h"
	789	#include "cvparam.h"
	790
	791	c-------------------------------------------------------------------
	792	c Compute icbmax.
	793	c-------------------------------------------------------------------
	794
	795	icbmax=2
	796	do 230 i=1,ncum
	797	icbmax=max(icbmax,icb(i))
	798	230 continue
	799
	800	c=====================================================================
	801	c --- CALCULATE CLOUD BASE MASS FLUX
	802	c=====================================================================
	803	c
	804	c tvpplcl = parcel temperature lifted adiabatically from level
	805	c icb-1 to the LCL.
	806	c tvaplcl = virtual temperature at the LCL.
	807	c
	808	do 610 i=1,ncum
	809	dtpbl(i)=0.0
	810	tvpplcl(i)=tvp(i,icb(i)-1)
	811	& -rrdtvp(i,icb(i)-1)(p(i,icb(i)-1)-plcl(i))
	812	& /(cpn(i,icb(i)-1)*p(i,icb(i)-1))
	813	tvaplcl(i)=tv(i,icb(i))
	814	& +(tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i,icb(i)))
	815	& /(p(i,icb(i))-p(i,icb(i)+1))
	816	610 continue
	817
	818	c-------------------------------------------------------------------
	819	c --- Interpolate difference between lifted parcel and
	820	c --- environmental temperatures to lifted condensation level
	821	c-------------------------------------------------------------------
	822	c
	823	c dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
	824	c
	825	do 630 k=minorig,icbmax
	826	do 620 i=1,ncum
	827	if((k.ge.nk(i)).and.(k.le.(icb(i)-1)))then
	828	dtpbl(i)=dtpbl(i)+(tvp(i,k)-tv(i,k))*dph(i,k)
	829	endif
	830	620 continue
	831	630 continue
	832	do 640 i=1,ncum
	833	dtpbl(i)=dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
	834	dtmin(i)=tvpplcl(i)-tvaplcl(i)+dtmax+dtpbl(i)
	835	640 continue
	836	c
	837	c-------------------------------------------------------------------
	838	c --- Adjust cloud base mass flux
	839	c-------------------------------------------------------------------
	840	c
	841	do 650 i=1,ncum
	842	work(i)=cbmf(i)
	843	cbmf(i)=max(0.0,(1.0-damp)cbmf(i)+0.1alpha*dtmin(i))
	844	if((work(i).eq.0.0).and.(cbmf(i).eq.0.0))then
	845	iflag(i)=3
	846	endif
	847	650 continue
	848
	849	return
	850	end
	851
	852	SUBROUTINE cv_mixing(nloc,ncum,nd,icb,nk,inb,inb1
	853	: ,ph,t,q,qs,u,v,h,lv,qnk
	854	: ,hp,tv,tvp,ep,clw,cbmf
	855	: ,m,ment,qent,uent,vent,nent,sij,elij)
	856	implicit none
	857
	858	#include "cvthermo.h"
	859	#include "cvparam.h"
	860
	861	c inputs:
	862	integer ncum, nd, nloc
	863	integer icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
	864	real cbmf(nloc), qnk(nloc)
	865	real ph(nloc,nd+1)
	866	real t(nloc,nd), q(nloc,nd), qs(nloc,nd), lv(nloc,nd)
	867	real u(nloc,nd), v(nloc,nd), h(nloc,nd), hp(nloc,nd)
	868	real tv(nloc,nd), tvp(nloc,nd), ep(nloc,nd), clw(nloc,nd)
	869
	870	c outputs:
	871	integer nent(nloc,nd)
	872	real m(nloc,nd), ment(nloc,nd,nd), qent(nloc,nd,nd)
	873	real uent(nloc,nd,nd), vent(nloc,nd,nd)
	874	real sij(nloc,nd,nd), elij(nloc,nd,nd)
	875
	876	c local variables:
	877	integer i, j, k, ij
	878	integer num1, num2
	879	real dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
	880	real alt, qp1, smid, sjmin, sjmax, delp, delm
	881	real work(nloc), asij(nloc), smin(nloc), scrit(nloc)
	882	real bsum(nloc,nd)
	883	logical lwork(nloc)
	884
	885	c=====================================================================
	886	c --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
	887	c=====================================================================
	888	c
	889	do 360 i=1,ncum*nlp
	890	nent(i,1)=0
	891	m(i,1)=0.0
	892	360 continue
	893	c
	894	do 400 k=1,nlp
	895	do 390 j=1,nlp
	896	do 385 i=1,ncum
	897	qent(i,k,j)=q(i,j)
	898	uent(i,k,j)=u(i,j)
	899	vent(i,k,j)=v(i,j)
	900	elij(i,k,j)=0.0
	901	ment(i,k,j)=0.0
	902	sij(i,k,j)=0.0
	903	385 continue
	904	390 continue
	905	400 continue
	906	c
	907	c-------------------------------------------------------------------
	908	c --- Calculate rates of mixing, m(i)
	909	c-------------------------------------------------------------------
	910	c
	911	call zilch(work,ncum)
	912	c
	913	do 670 j=minorig+1,nl
	914	do 660 i=1,ncum
	915	if((j.ge.(icb(i)+1)).and.(j.le.inb(i)))then
	916	k=min(j,inb1(i))
	917	dbo=abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1))
	918	& +entp0.04(ph(i,k)-ph(i,k+1))
	919	work(i)=work(i)+dbo
	920	m(i,j)=cbmf(i)*dbo
	921	endif
	922	660 continue
	923	670 continue
	924	do 690 k=minorig+1,nl
	925	do 680 i=1,ncum
	926	if((k.ge.(icb(i)+1)).and.(k.le.inb(i)))then
	927	m(i,k)=m(i,k)/work(i)
	928	endif
	929	680 continue
	930	690 continue
	931	c
	932	c
	933	c=====================================================================
	934	c --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
	935	c --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
	936	c --- FRACTION (sij)
	937	c=====================================================================
	938	c
	939	c
	940	do 750 i=minorig+1,nl
	941	do 710 j=minorig+1,nl
	942	do 700 ij=1,ncum
	943	if((i.ge.(icb(ij)+1)).and.(j.ge.icb(ij))
	944	& .and.(i.le.inb(ij)).and.(j.le.inb(ij)))then
	945	qti=qnk(ij)-ep(ij,i)*clw(ij,i)
	946	bf2=1.+lv(ij,j)lv(ij,j)qs(ij,j)
	947	& /(rrvt(ij,j)t(ij,j)*cpd)
	948	anum=h(ij,j)-hp(ij,i)+(cpv-cpd)t(ij,j)(qti-q(ij,j))
	949	denom=h(ij,i)-hp(ij,i)+(cpd-cpv)(q(ij,i)-qti)t(ij,j)
	950	dei=denom
	951	if(abs(dei).lt.0.01)dei=0.01
	952	sij(ij,i,j)=anum/dei
	953	sij(ij,i,i)=1.0
	954	altem=sij(ij,i,j)q(ij,i)+(1.-sij(ij,i,j))qti-qs(ij,j)
	955	altem=altem/bf2
	956	cwat=clw(ij,j)*(1.-ep(ij,j))
	957	stemp=sij(ij,i,j)
	958	if((stemp.lt.0.0.or.stemp.gt.1.0.or.
	959	1 altem.gt.cwat).and.j.gt.i)then
	960	anum=anum-lv(ij,j)(qti-qs(ij,j)-cwatbf2)
	961	denom=denom+lv(ij,j)*(q(ij,i)-qti)
	962	if(abs(denom).lt.0.01)denom=0.01
	963	sij(ij,i,j)=anum/denom
	964	altem=sij(ij,i,j)q(ij,i)+(1.-sij(ij,i,j))qti-qs(ij,j)
	965	altem=altem-(bf2-1.)*cwat
	966	endif
	967	if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
	968	qent(ij,i,j)=sij(ij,i,j)*q(ij,i)
	969	& +(1.-sij(ij,i,j))*qti
	970	uent(ij,i,j)=sij(ij,i,j)*u(ij,i)
	971	& +(1.-sij(ij,i,j))*u(ij,nk(ij))
	972	vent(ij,i,j)=sij(ij,i,j)*v(ij,i)
	973	& +(1.-sij(ij,i,j))*v(ij,nk(ij))
	974	elij(ij,i,j)=altem
	975	elij(ij,i,j)=max(0.0,elij(ij,i,j))
	976	ment(ij,i,j)=m(ij,i)/(1.-sij(ij,i,j))
	977	nent(ij,i)=nent(ij,i)+1
	978	endif
	979	sij(ij,i,j)=max(0.0,sij(ij,i,j))
	980	sij(ij,i,j)=min(1.0,sij(ij,i,j))
	981	endif
	982	700 continue
	983	710 continue
	984	c
	985	c * If no air can entrain at level i assume that updraft detrains *
	986	c * at that level and calculate detrained air flux and properties *
	987	c
	988	do 740 ij=1,ncum
	989	if((i.ge.(icb(ij)+1)).and.(i.le.inb(ij))
	990	& .and.(nent(ij,i).eq.0))then
	991	ment(ij,i,i)=m(ij,i)
	992	qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
	993	uent(ij,i,i)=u(ij,nk(ij))
	994	vent(ij,i,i)=v(ij,nk(ij))
	995	elij(ij,i,i)=clw(ij,i)
	996	sij(ij,i,i)=1.0
	997	endif
	998	740 continue
	999	750 continue
	1000	c
	1001	do 770 i=1,ncum
	1002	sij(i,inb(i),inb(i))=1.0
	1003	770 continue
	1004	c
	1005	c=====================================================================
	1006	c --- NORMALIZE ENTRAINED AIR MASS FLUXES
	1007	c --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
	1008	c=====================================================================
	1009	c
	1010	call zilch(bsum,ncum*nlp)
	1011	do 780 ij=1,ncum
	1012	lwork(ij)=.false.
	1013	780 continue
	1014	do 789 i=minorig+1,nl
	1015	c
	1016	num1=0
	1017	do 779 ij=1,ncum
	1018	if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))num1=num1+1
	1019	779 continue
	1020	if(num1.le.0)go to 789
	1021	c
	1022	do 781 ij=1,ncum
	1023	if((i.ge.icb(ij)+1).and.(i.le.inb(ij)))then
	1024	lwork(ij)=(nent(ij,i).ne.0)
	1025	qp1=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
	1026	anum=h(ij,i)-hp(ij,i)-lv(ij,i)*(qp1-qs(ij,i))
	1027	denom=h(ij,i)-hp(ij,i)+lv(ij,i)*(q(ij,i)-qp1)
	1028	if(abs(denom).lt.0.01)denom=0.01
	1029	scrit(ij)=anum/denom
	1030	alt=qp1-qs(ij,i)+scrit(ij)*(q(ij,i)-qp1)
	1031	if(scrit(ij).lt.0.0.or.alt.lt.0.0)scrit(ij)=1.0
	1032	asij(ij)=0.0
	1033	smin(ij)=1.0
	1034	endif
	1035	781 continue
	1036	do 783 j=minorig,nl
	1037	c
	1038	num2=0
	1039	do 778 ij=1,ncum
	1040	if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
	1041	& .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
	1042	& .and.lwork(ij))num2=num2+1
	1043	778 continue
	1044	if(num2.le.0)go to 783
	1045	c
	1046	do 782 ij=1,ncum
	1047	if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
	1048	& .and.(j.ge.icb(ij)).and.(j.le.inb(ij)).and.lwork(ij))then
	1049	if(sij(ij,i,j).gt.0.0.and.sij(ij,i,j).lt.0.9)then
	1050	if(j.gt.i)then
	1051	smid=min(sij(ij,i,j),scrit(ij))
	1052	sjmax=smid
	1053	sjmin=smid
	1054	if(smid.lt.smin(ij)
	1055	& .and.sij(ij,i,j+1).lt.smid)then
	1056	smin(ij)=smid
	1057	sjmax=min(sij(ij,i,j+1),sij(ij,i,j),scrit(ij))
	1058	sjmin=max(sij(ij,i,j-1),sij(ij,i,j))
	1059	sjmin=min(sjmin,scrit(ij))
	1060	endif
	1061	else
	1062	sjmax=max(sij(ij,i,j+1),scrit(ij))
	1063	smid=max(sij(ij,i,j),scrit(ij))
	1064	sjmin=0.0
	1065	if(j.gt.1)sjmin=sij(ij,i,j-1)
	1066	sjmin=max(sjmin,scrit(ij))
	1067	endif
	1068	delp=abs(sjmax-smid)
	1069	delm=abs(sjmin-smid)
	1070	asij(ij)=asij(ij)+(delp+delm)
	1071	& *(ph(ij,j)-ph(ij,j+1))
	1072	ment(ij,i,j)=ment(ij,i,j)*(delp+delm)
	1073	& *(ph(ij,j)-ph(ij,j+1))
	1074	endif
	1075	endif
	1076	782 continue
	1077	783 continue
	1078	do 784 ij=1,ncum
	1079	if((i.ge.icb(ij)+1).and.(i.le.inb(ij)).and.lwork(ij))then
	1080	asij(ij)=max(1.0e-21,asij(ij))
	1081	asij(ij)=1.0/asij(ij)
	1082	bsum(ij,i)=0.0
	1083	endif
	1084	784 continue
	1085	do 786 j=minorig,nl+1
	1086	do 785 ij=1,ncum
	1087	if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
	1088	& .and.(j.ge.icb(ij)).and.(j.le.inb(ij))
	1089	& .and.lwork(ij))then
	1090	ment(ij,i,j)=ment(ij,i,j)*asij(ij)
	1091	bsum(ij,i)=bsum(ij,i)+ment(ij,i,j)
	1092	endif
	1093	785 continue
	1094	786 continue
	1095	do 787 ij=1,ncum
	1096	if((i.ge.icb(ij)+1).and.(i.le.inb(ij))
	1097	& .and.(bsum(ij,i).lt.1.0e-18).and.lwork(ij))then
	1098	nent(ij,i)=0
	1099	ment(ij,i,i)=m(ij,i)
	1100	qent(ij,i,i)=q(ij,nk(ij))-ep(ij,i)*clw(ij,i)
	1101	uent(ij,i,i)=u(ij,nk(ij))
	1102	vent(ij,i,i)=v(ij,nk(ij))
	1103	elij(ij,i,i)=clw(ij,i)
	1104	sij(ij,i,i)=1.0
	1105	endif
	1106	787 continue
	1107	789 continue
	1108	c
	1109	return
	1110	end
	1111
	1112	SUBROUTINE cv_unsat(nloc,ncum,nd,inb,t,q,qs,gz,u,v,p,ph
	1113	: ,h,lv,ep,sigp,clw,m,ment,elij
	1114	: ,iflag,mp,qp,up,vp,wt,water,evap)
	1115	implicit none
	1116
	1117
	1118	#include "cvthermo.h"
	1119	#include "cvparam.h"
	1120
	1121	c inputs:
	1122	integer ncum, nd, nloc
	1123	integer inb(nloc)
	1124	real t(nloc,nd), q(nloc,nd), qs(nloc,nd)
	1125	real gz(nloc,nd), u(nloc,nd), v(nloc,nd)
	1126	real p(nloc,nd), ph(nloc,nd+1), h(nloc,nd)
	1127	real lv(nloc,nd), ep(nloc,nd), sigp(nloc,nd), clw(nloc,nd)
	1128	real m(nloc,nd), ment(nloc,nd,nd), elij(nloc,nd,nd)
	1129
	1130	c outputs:
	1131	integer iflag(nloc) ! also an input
	1132	real mp(nloc,nd), qp(nloc,nd), up(nloc,nd), vp(nloc,nd)
	1133	real water(nloc,nd), evap(nloc,nd), wt(nloc,nd)
	1134
	1135	c local variables:
	1136	integer i,j,k,ij,num1
	1137	integer jtt(nloc)
	1138	real awat, coeff, qsm, afac, sigt, b6, c6, revap
	1139	real dhdp, fac, qstm, rat
	1140	real wdtrain(nloc)
	1141	logical lwork(nloc)
	1142
	1143	c=====================================================================
	1144	c --- PRECIPITATING DOWNDRAFT CALCULATION
	1145	c=====================================================================
	1146	c
	1147	c Initializations:
	1148	c
	1149	do i = 1, ncum
	1150	do k = 1, nl+1
	1151	wt(i,k) = omtsnow
	1152	mp(i,k) = 0.0
	1153	evap(i,k) = 0.0
	1154	water(i,k) = 0.0
	1155	enddo
	1156	enddo
	1157
	1158	do 420 i=1,ncum
	1159	qp(i,1)=q(i,1)
	1160	up(i,1)=u(i,1)
	1161	vp(i,1)=v(i,1)
	1162	420 continue
	1163
	1164	do 440 k=2,nl+1
	1165	do 430 i=1,ncum
	1166	qp(i,k)=q(i,k-1)
	1167	up(i,k)=u(i,k-1)
	1168	vp(i,k)=v(i,k-1)
	1169	430 continue
	1170	440 continue
	1171
	1172
	1173	c * Check whether ep(inb)=0, if so, skip precipitating *
	1174	c * downdraft calculation *
	1175	c
	1176	c
	1177	c * Integrate liquid water equation to find condensed water *
	1178	c * and condensed water flux *
	1179	c
	1180	c
	1181	do 890 i=1,ncum
	1182	jtt(i)=2
	1183	if(ep(i,inb(i)).le.0.0001)iflag(i)=2
	1184	if(iflag(i).eq.0)then
	1185	lwork(i)=.true.
	1186	else
	1187	lwork(i)=.false.
	1188	endif
	1189	890 continue
	1190	c
	1191	c * Begin downdraft loop *
	1192	c
	1193	c
	1194	call zilch(wdtrain,ncum)
	1195	do 899 i=nl+1,1,-1
	1196	c
	1197	num1=0
	1198	do 879 ij=1,ncum
	1199	if((i.le.inb(ij)).and.lwork(ij))num1=num1+1
	1200	879 continue
	1201	if(num1.le.0)go to 899
	1202	c
	1203	c
	1204	c * Calculate detrained precipitation *
	1205	c
	1206	do 891 ij=1,ncum
	1207	if((i.le.inb(ij)).and.(lwork(ij)))then
	1208	wdtrain(ij)=gep(ij,i)m(ij,i)*clw(ij,i)
	1209	endif
	1210	891 continue
	1211	c
	1212	if(i.gt.1)then
	1213	do 893 j=1,i-1
	1214	do 892 ij=1,ncum
	1215	if((i.le.inb(ij)).and.(lwork(ij)))then
	1216	awat=elij(ij,j,i)-(1.-ep(ij,i))*clw(ij,i)
	1217	awat=max(0.0,awat)
	1218	wdtrain(ij)=wdtrain(ij)+gawatment(ij,j,i)
	1219	endif
	1220	892 continue
	1221	893 continue
	1222	endif
	1223	c
	1224	c * Find rain water and evaporation using provisional *
	1225	c * estimates of qp(i)and qp(i-1) *
	1226	c
	1227	c
	1228	c * Value of terminal velocity and coeffecient of evaporation for snow *
	1229	c
	1230	do 894 ij=1,ncum
	1231	if((i.le.inb(ij)).and.(lwork(ij)))then
	1232	coeff=coeffs
	1233	wt(ij,i)=omtsnow
	1234	c
	1235	c * Value of terminal velocity and coeffecient of evaporation for rain *
	1236	c
	1237	if(t(ij,i).gt.273.0)then
	1238	coeff=coeffr
	1239	wt(ij,i)=omtrain
	1240	endif
	1241	qsm=0.5*(q(ij,i)+qp(ij,i+1))
	1242	afac=coeffph(ij,i)(qs(ij,i)-qsm)
	1243	& /(1.0e4+2.0e3ph(ij,i)qs(ij,i))
	1244	afac=max(afac,0.0)
	1245	sigt=sigp(ij,i)
	1246	sigt=max(0.0,sigt)
	1247	sigt=min(1.0,sigt)
	1248	b6=100.(ph(ij,i)-ph(ij,i+1))sigt*afac/wt(ij,i)
	1249	c6=(water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij,i)
	1250	revap=0.5(-b6+sqrt(b6b6+4.*c6))
	1251	evap(ij,i)=sigtafacrevap
	1252	water(ij,i)=revap*revap
	1253	c
	1254	c * Calculate precipitating downdraft mass flux under *
	1255	c * hydrostatic approximation *
	1256	c
	1257	if(i.gt.1)then
	1258	dhdp=(h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
	1259	dhdp=max(dhdp,10.0)
	1260	mp(ij,i)=100.ginvlv(ij,i)sigdevap(ij,i)/dhdp
	1261	mp(ij,i)=max(mp(ij,i),0.0)
	1262	c
	1263	c * Add small amount of inertia to downdraft *
	1264	c
	1265	fac=20.0/(ph(ij,i-1)-ph(ij,i))
	1266	mp(ij,i)=(fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
	1267	c
	1268	c * Force mp to decrease linearly to zero *
	1269	c * between about 950 mb and the surface *
	1270	c
	1271	if(p(ij,i).gt.(0.949*p(ij,1)))then
	1272	jtt(ij)=max(jtt(ij),i)
	1273	mp(ij,i)=mp(ij,jtt(ij))*(p(ij,1)-p(ij,i))
	1274	& /(p(ij,1)-p(ij,jtt(ij)))
	1275	endif
	1276	endif
	1277	c
	1278	c * Find mixing ratio of precipitating downdraft *
	1279	c
	1280	if(i.ne.inb(ij))then
	1281	if(i.eq.1)then
	1282	qstm=qs(ij,1)
	1283	else
	1284	qstm=qs(ij,i-1)
	1285	endif
	1286	if(mp(ij,i).gt.mp(ij,i+1))then
	1287	rat=mp(ij,i+1)/mp(ij,i)
	1288	qp(ij,i)=qp(ij,i+1)rat+q(ij,i)(1.0-rat)+100.ginv
	1289	& sigd(ph(ij,i)-ph(ij,i+1))(evap(ij,i)/mp(ij,i))
	1290	up(ij,i)=up(ij,i+1)rat+u(ij,i)(1.-rat)
	1291	vp(ij,i)=vp(ij,i+1)rat+v(ij,i)(1.-rat)
	1292	else
	1293	if(mp(ij,i+1).gt.0.0)then
	1294	qp(ij,i)=(gz(ij,i+1)-gz(ij,i)
	1295	& +qp(ij,i+1)*(lv(ij,i+1)+t(ij,i+1)
	1296	& (cl-cpd))+cpd(t(ij,i+1)-t(ij,i)))
	1297	& /(lv(ij,i)+t(ij,i)*(cl-cpd))
	1298	up(ij,i)=up(ij,i+1)
	1299	vp(ij,i)=vp(ij,i+1)
	1300	endif
	1301	endif
	1302	qp(ij,i)=min(qp(ij,i),qstm)
	1303	qp(ij,i)=max(qp(ij,i),0.0)
	1304	endif
	1305	endif
	1306	894 continue
	1307	899 continue
	1308	c
	1309	return
	1310	end
	1311
	1312	SUBROUTINE cv_yield(nloc,ncum,nd,nk,icb,inb,delt
	1313	: ,t,q,u,v,gz,p,ph,h,hp,lv,cpn
	1314	: ,ep,clw,frac,m,mp,qp,up,vp
	1315	: ,wt,water,evap
	1316	: ,ment,qent,uent,vent,nent,elij
	1317	: ,tv,tvp
	1318	o ,iflag,wd,qprime,tprime
	1319	o ,precip,cbmf,ft,fq,fu,fv,Ma,qcondc)
	1320	implicit none
	1321
	1322	#include "cvthermo.h"
	1323	#include "cvparam.h"
	1324
	1325	c inputs
	1326	integer ncum, nd, nloc
	1327	integer nk(nloc), icb(nloc), inb(nloc)
	1328	integer nent(nloc,nd)
	1329	real delt
	1330	real t(nloc,nd), q(nloc,nd), u(nloc,nd), v(nloc,nd)
	1331	real gz(nloc,nd)
	1332	real p(nloc,nd), ph(nloc,nd+1), h(nloc,nd)
	1333	real hp(nloc,nd), lv(nloc,nd)
	1334	real cpn(nloc,nd), ep(nloc,nd), clw(nloc,nd), frac(nloc)
	1335	real m(nloc,nd), mp(nloc,nd), qp(nloc,nd)
	1336	real up(nloc,nd), vp(nloc,nd)
	1337	real wt(nloc,nd), water(nloc,nd), evap(nloc,nd)
	1338	real ment(nloc,nd,nd), qent(nloc,nd,nd), elij(nloc,nd,nd)
	1339	real uent(nloc,nd,nd), vent(nloc,nd,nd)
	1340	real tv(nloc,nd), tvp(nloc,nd)
	1341
	1342	c outputs
	1343	integer iflag(nloc) ! also an input
	1344	real cbmf(nloc) ! also an input
	1345	real wd(nloc), tprime(nloc), qprime(nloc)
	1346	real precip(nloc)
	1347	real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
	1348	real Ma(nloc,nd)
	1349	real qcondc(nloc,nd)
	1350
	1351	c local variables
	1352	integer i,j,ij,k,num1
	1353	real dpinv,cpinv,awat,fqold,ftold,fuold,fvold,delti
	1354	real work(nloc), am(nloc),amp1(nloc),ad(nloc)
	1355	real ents(nloc), uav(nloc),vav(nloc),lvcp(nloc,nd)
	1356	real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd) ! cld
	1357	real siga(nloc,nd), ax(nloc,nd), mac(nloc,nd) ! cld
	1358
	1359
	1360	c -- initializations:
	1361
	1362	delti = 1.0/delt
	1363
	1364	do 160 i=1,ncum
	1365	precip(i)=0.0
	1366	wd(i)=0.0
	1367	tprime(i)=0.0
	1368	qprime(i)=0.0
	1369	do 170 k=1,nl+1
	1370	ft(i,k)=0.0
	1371	fu(i,k)=0.0
	1372	fv(i,k)=0.0
	1373	fq(i,k)=0.0
	1374	lvcp(i,k)=lv(i,k)/cpn(i,k)
	1375	qcondc(i,k)=0.0 ! cld
	1376	qcond(i,k)=0.0 ! cld
	1377	nqcond(i,k)=0.0 ! cld
	1378	170 continue
	1379	160 continue
	1380
	1381	c
	1382	c * Calculate surface precipitation in mm/day *
	1383	c
	1384	do 1190 i=1,ncum
	1385	if(iflag(i).le.1)then
	1386	cc precip(i)=precip(i)+wt(i,1)sigdwater(i,1)3600.24000.
	1387	cc & /(rowl*g)
	1388	cc precip(i)=precip(i)*delt/86400.
	1389	precip(i) = wt(i,1)sigdwater(i,1)*86400/g
	1390	endif
	1391	1190 continue
	1392	c
	1393	c
	1394	c * Calculate downdraft velocity scale and surface temperature and *
	1395	c * water vapor fluctuations *
	1396	c
	1397	do i=1,ncum
	1398	wd(i)=betadabs(mp(i,icb(i)))0.01rrdt(i,icb(i))
	1399	: /(sigd*p(i,icb(i)))
	1400	qprime(i)=0.5*(qp(i,1)-q(i,1))
	1401	tprime(i)=lv0*qprime(i)/cpd
	1402	enddo
	1403	c
	1404	c * Calculate tendencies of lowest level potential temperature *
	1405	c * and mixing ratio *
	1406	c
	1407	do 1200 i=1,ncum
	1408	work(i)=0.01/(ph(i,1)-ph(i,2))
	1409	am(i)=0.0
	1410	1200 continue
	1411	do 1220 k=2,nl
	1412	do 1210 i=1,ncum
	1413	if((nk(i).eq.1).and.(k.le.inb(i)).and.(nk(i).eq.1))then
	1414	am(i)=am(i)+m(i,k)
	1415	endif
	1416	1210 continue
	1417	1220 continue
	1418	do 1240 i=1,ncum
	1419	if((gwork(i)am(i)).ge.delti)iflag(i)=1
	1420	ft(i,1)=ft(i,1)+gwork(i)am(i)*(t(i,2)-t(i,1)
	1421	& +(gz(i,2)-gz(i,1))/cpn(i,1))
	1422	ft(i,1)=ft(i,1)-lvcp(i,1)sigdevap(i,1)
	1423	ft(i,1)=ft(i,1)+sigdwt(i,2)(cl-cpd)water(i,2)(t(i,2)
	1424	& -t(i,1))*work(i)/cpn(i,1)
	1425	fq(i,1)=fq(i,1)+gmp(i,2)(qp(i,2)-q(i,1))*
	1426	& work(i)+sigd*evap(i,1)
	1427	fq(i,1)=fq(i,1)+gam(i)(q(i,2)-q(i,1))*work(i)
	1428	fu(i,1)=fu(i,1)+gwork(i)(mp(i,2)*(up(i,2)-u(i,1))
	1429	& +am(i)*(u(i,2)-u(i,1)))
	1430	fv(i,1)=fv(i,1)+gwork(i)(mp(i,2)*(vp(i,2)-v(i,1))
	1431	& +am(i)*(v(i,2)-v(i,1)))
	1432	1240 continue
	1433	do 1260 j=2,nl
	1434	do 1250 i=1,ncum
	1435	if(j.le.inb(i))then
	1436	fq(i,1)=fq(i,1)
	1437	& +gwork(i)ment(i,j,1)*(qent(i,j,1)-q(i,1))
	1438	fu(i,1)=fu(i,1)
	1439	& +gwork(i)ment(i,j,1)*(uent(i,j,1)-u(i,1))
	1440	fv(i,1)=fv(i,1)
	1441	& +gwork(i)ment(i,j,1)*(vent(i,j,1)-v(i,1))
	1442	endif
	1443	1250 continue
	1444	1260 continue
	1445	c
	1446	c * Calculate tendencies of potential temperature and mixing ratio *
	1447	c * at levels above the lowest level *
	1448	c
	1449	c * First find the net saturated updraft and downdraft mass fluxes *
	1450	c * through each level *
	1451	c
	1452	do 1500 i=2,nl+1
	1453	c
	1454	num1=0
	1455	do 1265 ij=1,ncum
	1456	if(i.le.inb(ij))num1=num1+1
	1457	1265 continue
	1458	if(num1.le.0)go to 1500
	1459	c
	1460	call zilch(amp1,ncum)
	1461	call zilch(ad,ncum)
	1462	c
	1463	do 1280 k=i+1,nl+1
	1464	do 1270 ij=1,ncum
	1465	if((i.ge.nk(ij)).and.(i.le.inb(ij))
	1466	& .and.(k.le.(inb(ij)+1)))then
	1467	amp1(ij)=amp1(ij)+m(ij,k)
	1468	endif
	1469	1270 continue
	1470	1280 continue
	1471	c
	1472	do 1310 k=1,i
	1473	do 1300 j=i+1,nl+1
	1474	do 1290 ij=1,ncum
	1475	if((j.le.(inb(ij)+1)).and.(i.le.inb(ij)))then
	1476	amp1(ij)=amp1(ij)+ment(ij,k,j)
	1477	endif
	1478	1290 continue
	1479	1300 continue
	1480	1310 continue
	1481	do 1340 k=1,i-1
	1482	do 1330 j=i,nl+1
	1483	do 1320 ij=1,ncum
	1484	if((i.le.inb(ij)).and.(j.le.inb(ij)))then
	1485	ad(ij)=ad(ij)+ment(ij,j,k)
	1486	endif
	1487	1320 continue
	1488	1330 continue
	1489	1340 continue
	1490	c
	1491	do 1350 ij=1,ncum
	1492	if(i.le.inb(ij))then
	1493	dpinv=0.01/(ph(ij,i)-ph(ij,i+1))
	1494	cpinv=1.0/cpn(ij,i)
	1495	c
	1496	ft(ij,i)=ft(ij,i)
	1497	& +gdpinv(amp1(ij)*(t(ij,i+1)-t(ij,i)
	1498	& +(gz(ij,i+1)-gz(ij,i))*cpinv)
	1499	& -ad(ij)(t(ij,i)-t(ij,i-1)+(gz(ij,i)-gz(ij,i-1))cpinv))
	1500	& -sigdlvcp(ij,i)evap(ij,i)
	1501	ft(ij,i)=ft(ij,i)+gdpinvment(ij,i,i)*(hp(ij,i)-h(ij,i)+
	1502	& t(ij,i)(cpv-cpd)(q(ij,i)-qent(ij,i,i)))*cpinv
	1503	ft(ij,i)=ft(ij,i)+sigdwt(ij,i+1)(cl-cpd)water(ij,i+1)
	1504	& (t(ij,i+1)-t(ij,i))dpinvcpinv
	1505	fq(ij,i)=fq(ij,i)+gdpinv(amp1(ij)*(q(ij,i+1)-q(ij,i))-
	1506	& ad(ij)*(q(ij,i)-q(ij,i-1)))
	1507	fu(ij,i)=fu(ij,i)+gdpinv(amp1(ij)*(u(ij,i+1)-u(ij,i))-
	1508	& ad(ij)*(u(ij,i)-u(ij,i-1)))
	1509	fv(ij,i)=fv(ij,i)+gdpinv(amp1(ij)*(v(ij,i+1)-v(ij,i))-
	1510	& ad(ij)*(v(ij,i)-v(ij,i-1)))
	1511	endif
	1512	1350 continue
	1513	do 1370 k=1,i-1
	1514	do 1360 ij=1,ncum
	1515	if(i.le.inb(ij))then
	1516	awat=elij(ij,k,i)-(1.-ep(ij,i))*clw(ij,i)
	1517	awat=max(awat,0.0)
	1518	fq(ij,i)=fq(ij,i)
	1519	& +gdpinvment(ij,k,i)*(qent(ij,k,i)-awat-q(ij,i))
	1520	fu(ij,i)=fu(ij,i)
	1521	& +gdpinvment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
	1522	fv(ij,i)=fv(ij,i)
	1523	& +gdpinvment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
	1524	c (saturated updrafts resulting from mixing) ! cld
	1525	qcond(ij,i)=qcond(ij,i)+(elij(ij,k,i)-awat) ! cld
	1526	nqcond(ij,i)=nqcond(ij,i)+1. ! cld
	1527	endif
	1528	1360 continue
	1529	1370 continue
	1530	do 1390 k=i,nl+1
	1531	do 1380 ij=1,ncum
	1532	if((i.le.inb(ij)).and.(k.le.inb(ij)))then
	1533	fq(ij,i)=fq(ij,i)
	1534	& +gdpinvment(ij,k,i)*(qent(ij,k,i)-q(ij,i))
	1535	fu(ij,i)=fu(ij,i)
	1536	& +gdpinvment(ij,k,i)*(uent(ij,k,i)-u(ij,i))
	1537	fv(ij,i)=fv(ij,i)
	1538	& +gdpinvment(ij,k,i)*(vent(ij,k,i)-v(ij,i))
	1539	endif
	1540	1380 continue
	1541	1390 continue
	1542	do 1400 ij=1,ncum
	1543	if(i.le.inb(ij))then
	1544	fq(ij,i)=fq(ij,i)
	1545	& +sigdevap(ij,i)+g(mp(ij,i+1)*
	1546	& (qp(ij,i+1)-q(ij,i))
	1547	& -mp(ij,i)(qp(ij,i)-q(ij,i-1)))dpinv
	1548	fu(ij,i)=fu(ij,i)
	1549	& +g(mp(ij,i+1)(up(ij,i+1)-u(ij,i))-mp(ij,i)*
	1550	& (up(ij,i)-u(ij,i-1)))*dpinv
	1551	fv(ij,i)=fv(ij,i)
	1552	& +g(mp(ij,i+1)(vp(ij,i+1)-v(ij,i))-mp(ij,i)*
	1553	& (vp(ij,i)-v(ij,i-1)))*dpinv
	1554	C (saturated downdrafts resulting from mixing) ! cld
	1555	do k=i+1,inb(ij) ! cld
	1556	qcond(ij,i)=qcond(ij,i)+elij(ij,k,i) ! cld
	1557	nqcond(ij,i)=nqcond(ij,i)+1. ! cld
	1558	enddo ! cld
	1559	C (particular case: no detraining level is found) ! cld
	1560	if (nent(ij,i).eq.0) then ! cld
	1561	qcond(ij,i)=qcond(ij,i)+(1.-ep(ij,i))*clw(ij,i) ! cld
	1562	nqcond(ij,i)=nqcond(ij,i)+1. ! cld
	1563	endif ! cld
	1564	if (nqcond(ij,i).ne.0.) then ! cld
	1565	qcond(ij,i)=qcond(ij,i)/nqcond(ij,i) ! cld
	1566	endif ! cld
	1567	endif
	1568	1400 continue
	1569	1500 continue
	1570	c
	1571	c * Adjust tendencies at top of convection layer to reflect *
	1572	c * actual position of the level zero cape *
	1573	c
	1574	do 503 ij=1,ncum
	1575	fqold=fq(ij,inb(ij))
	1576	fq(ij,inb(ij))=fq(ij,inb(ij))*(1.-frac(ij))
	1577	fq(ij,inb(ij)-1)=fq(ij,inb(ij)-1)
	1578	& +frac(ij)fqold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
	1579	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*lv(ij,inb(ij))
	1580	& /lv(ij,inb(ij)-1)
	1581	ftold=ft(ij,inb(ij))
	1582	ft(ij,inb(ij))=ft(ij,inb(ij))*(1.-frac(ij))
	1583	ft(ij,inb(ij)-1)=ft(ij,inb(ij)-1)
	1584	& +frac(ij)ftold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
	1585	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))*cpn(ij,inb(ij))
	1586	& /cpn(ij,inb(ij)-1)
	1587	fuold=fu(ij,inb(ij))
	1588	fu(ij,inb(ij))=fu(ij,inb(ij))*(1.-frac(ij))
	1589	fu(ij,inb(ij)-1)=fu(ij,inb(ij)-1)
	1590	& +frac(ij)fuold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
	1591	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
	1592	fvold=fv(ij,inb(ij))
	1593	fv(ij,inb(ij))=fv(ij,inb(ij))*(1.-frac(ij))
	1594	fv(ij,inb(ij)-1)=fv(ij,inb(ij)-1)
	1595	& +frac(ij)fvold((ph(ij,inb(ij))-ph(ij,inb(ij)+1))/
	1596	1 (ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
	1597	503 continue
	1598	c
	1599	c * Very slightly adjust tendencies to force exact *
	1600	c * enthalpy, momentum and tracer conservation *
	1601	c
	1602	do 682 ij=1,ncum
	1603	ents(ij)=0.0
	1604	uav(ij)=0.0
	1605	vav(ij)=0.0
	1606	do 681 i=1,inb(ij)
	1607	ents(ij)=ents(ij)
	1608	& +(cpn(ij,i)ft(ij,i)+lv(ij,i)fq(ij,i))*(ph(ij,i)-ph(ij,i+1))
	1609	uav(ij)=uav(ij)+fu(ij,i)*(ph(ij,i)-ph(ij,i+1))
	1610	vav(ij)=vav(ij)+fv(ij,i)*(ph(ij,i)-ph(ij,i+1))
	1611	681 continue
	1612	682 continue
	1613	do 683 ij=1,ncum
	1614	ents(ij)=ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
	1615	uav(ij)=uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
	1616	vav(ij)=vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
	1617	683 continue
	1618	do 642 ij=1,ncum
	1619	do 641 i=1,inb(ij)
	1620	ft(ij,i)=ft(ij,i)-ents(ij)/cpn(ij,i)
	1621	fu(ij,i)=(1.-cu)*(fu(ij,i)-uav(ij))
	1622	fv(ij,i)=(1.-cu)*(fv(ij,i)-vav(ij))
	1623	641 continue
	1624	642 continue
	1625	c
	1626	do 1810 k=1,nl+1
	1627	do 1800 i=1,ncum
	1628	if((q(i,k)+delt*fq(i,k)).lt.0.0)iflag(i)=10
	1629	1800 continue
	1630	1810 continue
	1631	c
	1632	c
	1633	do 1900 i=1,ncum
	1634	if(iflag(i).gt.2)then
	1635	precip(i)=0.0
	1636	cbmf(i)=0.0
	1637	endif
	1638	1900 continue
	1639	do 1920 k=1,nl
	1640	do 1910 i=1,ncum
	1641	if(iflag(i).gt.2)then
	1642	ft(i,k)=0.0
	1643	fq(i,k)=0.0
	1644	fu(i,k)=0.0
	1645	fv(i,k)=0.0
	1646	qcondc(i,k)=0.0 ! cld
	1647	endif
	1648	1910 continue
	1649	1920 continue
	1650
	1651	do k=1,nl+1
	1652	do i=1,ncum
	1653	Ma(i,k) = 0.
	1654	enddo
	1655	enddo
	1656	do k=nl,1,-1
	1657	do i=1,ncum
	1658	Ma(i,k) = Ma(i,k+1)+m(i,k)
	1659	enddo
	1660	enddo
	1661
	1662	c
	1663	c * diagnose the in-cloud mixing ratio * ! cld
	1664	c * of condensed water * ! cld
	1665	c ! cld
	1666	DO ij=1,ncum ! cld
	1667	do i=1,nd ! cld
	1668	mac(ij,i)=0.0 ! cld
	1669	wa(ij,i)=0.0 ! cld
	1670	siga(ij,i)=0.0 ! cld
	1671	enddo ! cld
	1672	do i=nk(ij),inb(ij) ! cld
	1673	do k=i+1,inb(ij)+1 ! cld
	1674	mac(ij,i)=mac(ij,i)+m(ij,k) ! cld
	1675	enddo ! cld
	1676	enddo ! cld
	1677	do i=icb(ij),inb(ij)-1 ! cld
	1678	ax(ij,i)=0. ! cld
	1679	do j=icb(ij),i ! cld
	1680	ax(ij,i)=ax(ij,i)+rrd*(tvp(ij,j)-tv(ij,j)) ! cld
	1681	: *(ph(ij,j)-ph(ij,j+1))/p(ij,j) ! cld
	1682	enddo ! cld
	1683	if (ax(ij,i).gt.0.0) then ! cld
	1684	wa(ij,i)=sqrt(2.*ax(ij,i)) ! cld
	1685	endif ! cld
	1686	enddo ! cld
	1687	do i=1,nl ! cld
	1688	if (wa(ij,i).gt.0.0) ! cld
	1689	: siga(ij,i)=mac(ij,i)/wa(ij,i) ! cld
	1690	: rrdtvp(ij,i)/p(ij,i)/100./delta ! cld
	1691	siga(ij,i) = min(siga(ij,i),1.0) ! cld
	1692	qcondc(ij,i)=siga(ij,i)clw(ij,i)(1.-ep(ij,i)) ! cld
	1693	: + (1.-siga(ij,i))*qcond(ij,i) ! cld
	1694	enddo ! cld
	1695	ENDDO ! cld
	1696
	1697	return
	1698	end
	1699
	1700	SUBROUTINE cv_uncompress(nloc,len,ncum,nd,idcum
	1701	: ,iflag
	1702	: ,precip,cbmf
	1703	: ,ft,fq,fu,fv
	1704	: ,Ma,qcondc
	1705	: ,iflag1
	1706	: ,precip1,cbmf1
	1707	: ,ft1,fq1,fu1,fv1
	1708	: ,Ma1,qcondc1
	1709	: )
	1710	implicit none
	1711
	1712	#include "cvparam.h"
	1713
	1714	c inputs:
	1715	integer len, ncum, nd, nloc
	1716	integer idcum(nloc)
	1717	integer iflag(nloc)
	1718	real precip(nloc), cbmf(nloc)
	1719	real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
	1720	real Ma(nloc,nd)
	1721	real qcondc(nloc,nd) !cld
	1722
	1723	c outputs:
	1724	integer iflag1(len)
	1725	real precip1(len), cbmf1(len)
	1726	real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
	1727	real Ma1(len,nd)
	1728	real qcondc1(len,nd) !cld
	1729
	1730	c local variables:
	1731	integer i,k
	1732
	1733	do 2000 i=1,ncum
	1734	precip1(idcum(i))=precip(i)
	1735	cbmf1(idcum(i))=cbmf(i)
	1736	iflag1(idcum(i))=iflag(i)
	1737	2000 continue
	1738
	1739	do 2020 k=1,nl
	1740	do 2010 i=1,ncum
	1741	ft1(idcum(i),k)=ft(i,k)
	1742	fq1(idcum(i),k)=fq(i,k)
	1743	fu1(idcum(i),k)=fu(i,k)
	1744	fv1(idcum(i),k)=fv(i,k)
	1745	Ma1(idcum(i),k)=Ma(i,k)
	1746	qcondc1(idcum(i),k)=qcondc(i,k)
	1747	2010 continue
	1748	2020 continue
	1749
	1750	return
	1751	end
	1752

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

Download in other formats:

Original Format