Context Navigation

cv_routines.F @ 965

Last change on this file since 965 was 868, checked in by Laurent Fairhead, 17 years ago
Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4 par la version de la physique avec thermique. On garde le repertoire phylmd pour un petit moment pour que les utilisateurs ne soient pas trop perdus ... phy_IPCC_AR4 = phylmd LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 55.0 KB

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

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

Download in other formats:

Original Format