Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_cu_g3.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 165.4 KB

Rev	Line
[2759]	1	!WRF:MODEL_LAYER:PHYSICS
	2	!
	3
	4	MODULE module_cu_g3
	5
	6	CONTAINS
	7
	8	!-------------------------------------------------------------
	9	SUBROUTINE G3DRV( &
	10	DT,itimestep,DX &
	11	,rho,RAINCV,PRATEC &
	12	,U,V,t,W,q,p,pi &
	13	,dz8w,p8w,XLV,CP,G,r_v &
	14	,STEPCU,htop,hbot &
	15	,CU_ACT_FLAG,warm_rain &
	16	,APR_GR,APR_W,APR_MC,APR_ST,APR_AS &
	17	,APR_CAPMA,APR_CAPME,APR_CAPMI &
	18	,MASS_FLUX,XF_ENS,PR_ENS,HT,XLAND,gsw,edt_out &
	19	,GDC,GDC2 &
	20	,cugd_tten,cugd_qvten ,cugd_qcten &
	21	,cugd_ttens,cugd_qvtens,cugd_avedx,imomentum &
	22	,ensdim,maxiens,maxens,maxens2,maxens3,ichoice &
	23	,ids,ide, jds,jde, kds,kde &
	24	,ims,ime, jms,jme, kms,kme &
	25	,ips,ipe, jps,jpe, kps,kpe &
	26	,its,ite, jts,jte, kts,kte &
	27	,periodic_x,periodic_y &
	28	,RQVCUTEN,RQCCUTEN,RQICUTEN &
	29	,RQVFTEN,RTHFTEN,RTHCUTEN &
	30	,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
	31	)
	32	!-------------------------------------------------------------
	33	IMPLICIT NONE
	34	!-------------------------------------------------------------
	35	INTEGER, INTENT(IN ) :: &
	36	ids,ide, jds,jde, kds,kde, &
	37	ims,ime, jms,jme, kms,kme, &
	38	ips,ipe, jps,jpe, kps,kpe, &
	39	its,ite, jts,jte, kts,kte
	40	LOGICAL periodic_x,periodic_y
	41	integer, parameter :: ens4_spread = 3 ! max(3,cugd_avedx)
	42	integer, parameter :: ens4=ens4_spread*ens4_spread
	43
	44	integer, intent (in ) :: &
	45	ensdim,maxiens,maxens,maxens2,maxens3,ichoice
	46
	47	INTEGER, INTENT(IN ) :: STEPCU, ITIMESTEP,cugd_avedx,imomentum
	48	LOGICAL, INTENT(IN ) :: warm_rain
	49
	50	REAL, INTENT(IN ) :: XLV, R_v
	51	REAL, INTENT(IN ) :: CP,G
	52
	53	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
	54	INTENT(IN ) :: &
	55	U, &
	56	V, &
	57	W, &
	58	pi, &
	59	t, &
	60	q, &
	61	p, &
	62	dz8w, &
	63	p8w, &
	64	rho
	65	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
	66	OPTIONAL , &
	67	INTENT(INOUT ) :: &
	68	GDC,GDC2
	69
	70	REAL, DIMENSION( ims:ime , jms:jme ),INTENT(IN) :: GSW,HT,XLAND
	71	!
	72	REAL, INTENT(IN ) :: DT, DX
	73	!
	74
	75	REAL, DIMENSION( ims:ime , jms:jme ), &
	76	INTENT(INOUT) :: pratec,RAINCV, MASS_FLUX, &
	77	APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
	78	edt_out,APR_CAPMA,APR_CAPME,APR_CAPMI,htop,hbot
	79	!+lxz
	80	! REAL, DIMENSION( ims:ime , jms:jme ) :: & !, INTENT(INOUT) :: &
	81	! HTOP, &! highest model layer penetrated by cumulus since last reset in radiation_driver
	82	! HBOT ! lowest model layer penetrated by cumulus since last reset in radiation_driver
	83	! ! HBOT>HTOP follow physics leveling convention
	84
	85	LOGICAL, DIMENSION( ims:ime , jms:jme ), &
	86	INTENT(INOUT) :: CU_ACT_FLAG
	87
	88	!
	89	! Optionals
	90	!
	91	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
	92	OPTIONAL, &
	93	INTENT(INOUT) :: RTHFTEN, &
	94	cugd_tten,cugd_qvten,cugd_qcten, &
	95	cugd_ttens,cugd_qvtens, &
	96	RQVFTEN
	97
	98	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
	99	OPTIONAL, &
	100	INTENT(INOUT) :: &
	101	RTHCUTEN, &
	102	RQVCUTEN, &
	103	RQCCUTEN, &
	104	RQICUTEN
	105	!
	106	! Flags relating to the optional tendency arrays declared above
	107	! Models that carry the optional tendencies will provdide the
	108	! optional arguments at compile time; these flags all the model
	109	! to determine at run-time whether a particular tracer is in
	110	! use or not.
	111	!
	112	LOGICAL, OPTIONAL :: &
	113	F_QV &
	114	,F_QC &
	115	,F_QR &
	116	,F_QI &
	117	,F_QS
	118
	119
	120
	121	! LOCAL VARS
	122	real, dimension(ims:ime,jms:jme,1:ensdim),intent(inout) :: &
	123	xf_ens,pr_ens
	124	real, dimension ( its:ite , jts:jte , 1:ensdim) :: &
	125	massflni,xfi_ens,pri_ens
	126	REAL, DIMENSION( its:ite , jts:jte ) :: MASSI_FLX, &
	127	APRi_GR,APRi_W,APRi_MC,APRi_ST,APRi_AS, &
	128	edti_out,APRi_CAPMA,APRi_CAPME,APRi_CAPMI,gswi
	129	real, dimension (its:ite,kts:kte) :: &
	130	SUBT,SUBQ,OUTT,OUTQ,OUTQC,phh,subm,cupclw
	131	real, dimension (its:ite,kts:kte+1) :: phf
	132	real, dimension (its:ite) :: &
	133	pret, ter11, aa0, fp,xlandi
	134	!+lxz
	135	integer, dimension (its:ite) :: &
	136	kbcon, ktop
	137	!.lxz
	138	integer, dimension (its:ite,jts:jte) :: &
	139	iact_old_gr
	140	integer :: iens,ibeg,iend,jbeg,jend,n,nn,ens4n
	141	integer :: ibegh,iendh,jbegh,jendh
	142	integer :: ibegc,iendc,jbegc,jendc
	143
	144	!
	145	! basic environmental input includes moisture convergence (mconv)
	146	! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
	147	! convection for this call only and at that particular gridpoint
	148	!
	149	real, dimension (its:ite,kts:kte) :: &
	150	T2d,q2d,PO,P2d,US,VS,tn,qo
	151	real, dimension (ips-2:ipe+2,kps:kpe,jps-2:jpe+2) :: &
	152	ave_f_t,ave_f_q
	153	real, dimension (its:ite,kts:kte,1:ens4) :: &
	154	omeg,tx,qx
	155	real, dimension (its:ite) :: &
	156	Z1,PSUR,AAEQ,direction,cuten,umean,vmean,pmean
	157	real, dimension (its:ite,1:ens4) :: &
	158	mconv
	159
	160	INTEGER :: i,j,k,ICLDCK,ipr,jpr
	161	REAL :: tcrit,dp,dq,sub_spread,subcenter
	162	INTEGER :: itf,jtf,ktf,iss,jss,nbegin,nend
	163	INTEGER :: high_resolution
	164	REAL :: rkbcon,rktop !-lxz
	165	! ruc variable
	166	real, dimension (its:ite) :: tkm
	167
	168	high_resolution=0
	169	if(cugd_avedx.gt.1) high_resolution=1
	170	subcenter=0.
	171	! subcenter=1./float(cugd_avedx)
	172	sub_spread=max(1.,float(cugd_avedx*cugd_avedx-1))
	173	sub_spread=(1.-subcenter)/sub_spread
	174	iens=1
	175	ipr=37
	176	jpr=1
	177	ipr=0
	178	jpr=0
	179	! if(itimestep.eq.8)then
	180	! ipr=37
	181	! jpr=16
	182	! endif
	183	IF ( periodic_x ) THEN
	184	ibeg=max(its,ids)
	185	iend=min(ite,ide-1)
	186	ibegc=max(its,ids)
	187	iendc=min(ite,ide-1)
	188	ELSE
	189	ibeg=max(its,ids)
	190	iend=min(ite,ide-1)
	191	ibegc=max(its,ids+4)
	192	iendc=min(ite,ide-5)
	193	END IF
	194	IF ( periodic_y ) THEN
	195	jbeg=max(jts,jds)
	196	jend=min(jte,jde-1)
	197	jbegc=max(jts,jds)
	198	jendc=min(jte,jde-1)
	199	ELSE
	200	jbeg=max(jts,jds)
	201	jend=min(jte,jde-1)
	202	jbegc=max(jts,jds+4)
	203	jendc=min(jte,jde-5)
	204	END IF
	205	tcrit=258.
	206	ave_f_t=0.
	207	ave_f_q=0.
	208
	209	itf=MIN(ite,ide-1)
	210	ktf=MIN(kte,kde-1)
	211	jtf=MIN(jte,jde-1)
	212	!
	213	#if ( EM_CORE == 1 )
	214	if(high_resolution.eq.1)then
	215	!
	216	! calculate these on the halo...the incominh tendencies have been exchanged on a 24pt halo
	217	! only neede for high resolution run
	218	!
	219	ibegh=its
	220	jbegh=jts
	221	iendh=ite
	222	jendh=jte
	223	if(its.eq.ips)ibegh=max(its-1,ids)
	224	if(jts.eq.jps)jbegh=max(jts-1,jds)
	225	if(jte.eq.jpe)jendh=min(jte+1,jde-1)
	226	if(ite.eq.ipe)iendh=min(ite+1,ide-1)
	227	DO J = jbegh,jendh
	228	DO k= kts,ktf
	229	DO I= ibegh,iendh
	230	ave_f_t(i,k,j)=(rthften(i-1,k,j-1)+rthften(i-1,k,j) + rthften(i-1,k,j+1)+ &
	231	rthften(i,k,j-1) +rthften(i,k,j) +rthften(i,k,j+1)+ &
	232	rthften(i+1,k,j-1) +rthften(i+1,k,j) +rthften(i+1,k,j+1))/9.
	233	ave_f_q(i,k,j)=(rqvften(i-1,k,j-1)+rqvften(i-1,k,j) + rqvften(i-1,k,j+1)+ &
	234	rqvften(i,k,j-1) +rqvften(i,k,j) +rqvften(i,k,j+1)+ &
	235	rqvften(i+1,k,j-1) +rqvften(i+1,k,j) +rqvften(i+1,k,j+1))/9.
	236	! ave_f_t(i,k,j)=rthften(i,k,j)
	237	! ave_f_q(i,k,j)=rqvften(i,k,j)
	238	ENDDO
	239	ENDDO
	240	ENDDO
	241	endif
	242	#endif
	243	DO 100 J = jts,jtf
	244	DO n= 1,ensdim
	245	DO I= its,itf
	246	xfi_ens(i,j,n)=0.
	247	pri_ens(i,j,n)=0.
	248	! xfi_ens(i,j,n)=xf_ens(i,j,n)
	249	! pri_ens(i,j,n)=pr_ens(i,j,n)
	250	ENDDO
	251	ENDDO
	252	DO I= its,itf
	253	kbcon(i)=0
	254	ktop(i)=0
	255	tkm(i)=0.
	256	iact_old_gr(i,j)=0
	257	mass_flux(i,j)=0.
	258	massi_flx(i,j)=0.
	259	raincv(i,j)=0.
	260	pratec (i,j)=0.
	261	edt_out(i,j)=0.
	262	edti_out(i,j)=0.
	263	gswi(i,j)=gsw(i,j)
	264	xlandi(i)=xland(i,j)
	265	APRi_GR(i,j)=apr_gr(i,j)
	266	APRi_w(i,j)=apr_w(i,j)
	267	APRi_mc(i,j)=apr_mc(i,j)
	268	APRi_st(i,j)=apr_st(i,j)
	269	APRi_as(i,j)=apr_as(i,j)
	270	APRi_capma(i,j)=apr_capma(i,j)
	271	APRi_capme(i,j)=apr_capme(i,j)
	272	APRi_capmi(i,j)=apr_capmi(i,j)
	273	CU_ACT_FLAG(i,j) = .true.
	274	ENDDO
	275	do k=kts,kte
	276	DO I= its,itf
	277	cugd_tten(i,k,j)=0.
	278	cugd_ttens(i,k,j)=0.
	279	cugd_qvten(i,k,j)=0.
	280	cugd_qvtens(i,k,j)=0.
	281	cugd_qcten(i,k,j)=0.
	282	ENDDO
	283	ENDDO
	284	DO n=1,ens4
	285	DO I= its,itf
	286	mconv(i,n)=0.
	287	ENDDO
	288	do k=kts,kte
	289	DO I= its,itf
	290	omeg(i,k,n)=0.
	291	tx(i,k,n)=0.
	292	qx(i,k,n)=0.
	293	ENDDO
	294	ENDDO
	295	ENDDO
	296	DO k=1,ensdim
	297	DO I= its,itf
	298	massflni(i,j,k)=0.
	299	ENDDO
	300	ENDDO
	301	#if ( EM_CORE == 1 )
	302	! hydrostatic pressure, first on full levels
	303	DO I=ITS,ITF
	304	phf(i,1) = p8w(i,1,j)
	305	ENDDO
	306	! integrate up, dp = -rho * g * dz
	307	DO K=kts+1,ktf+1
	308	DO I=ITS,ITF
	309	phf(i,k) = phf(i,k-1) - rho(i,k-1,j) * g * dz8w(i,k-1,j)
	310	ENDDO
	311	ENDDO
	312	! scale factor so that pressure is not zero after integration
	313	DO I=ITS,ITF
	314	fp(i) = (p8w(i,kts,j)-p8w(i,kte,j))/(phf(i,kts)-phf(i,kte))
	315	ENDDO
	316	! re-integrate up, dp = -rho * g * dz * scale_factor
	317	DO K=kts+1,ktf+1
	318	DO I=ITS,ITF
	319	phf(i,k) = phf(i,k-1) - rho(i,k-1,j) * g * dz8w(i,k-1,j) * fp(i)
	320	ENDDO
	321	ENDDO
	322	! put hydrostatic pressure on half levels
	323	DO K=kts,ktf
	324	DO I=ITS,ITF
	325	phh(i,k) = (phf(i,k) + phf(i,k+1))*0.5
	326	ENDDO
	327	ENDDO
	328
	329	#endif
	330	DO I=ITS,ITF
	331	#if ( EM_CORE == 1 )
	332	PSUR(I)=p8w(I,1,J)*.01
	333	#endif
	334	#if ( NMM_CORE == 1 )
	335	PSUR(I)=p(I,1,J)*.01
	336	#endif
	337	! PSUR(I)=p(I,1,J)*.01
	338	TER11(I)=HT(i,j)
	339	aaeq(i)=0.
	340	direction(i)=0.
	341	pret(i)=0.
	342	umean(i)=0.
	343	vmean(i)=0.
	344	pmean(i)=0.
	345	ENDDO
	346	DO K=kts,ktf
	347	DO I=ITS,ITF
	348	#if ( EM_CORE == 1 )
	349	po(i,k)=phh(i,k)*.01
	350	#endif
	351
	352	#if ( NMM_CORE == 1 )
	353	po(i,k)=p(i,k,j)*.01
	354	#endif
	355	subm(i,k)=0.
	356	P2d(I,K)=PO(i,k)
	357	US(I,K) =u(i,k,j)
	358	VS(I,K) =v(i,k,j)
	359	T2d(I,K)=t(i,k,j)
	360	q2d(I,K)=q(i,k,j)
	361	IF(Q2d(I,K).LT.1.E-08)Q2d(I,K)=1.E-08
	362	SUBT(I,K)=0.
	363	SUBQ(I,K)=0.
	364	OUTT(I,K)=0.
	365	OUTQ(I,K)=0.
	366	OUTQC(I,K)=0.
	367	TN(I,K)=t2d(i,k)+RTHFTEN(i,k,j)*dt
	368	QO(I,K)=q2d(i,k)+RQVFTEN(i,k,j)*dt
	369	if(high_resolution.eq.1)then
	370	TN(I,K)=t2d(i,k)+ave_f_t(i,k,j)*dt
	371	QO(I,K)=q2d(i,k)+ave_f_q(i,k,j)*dt
	372	endif
	373	IF(TN(I,K).LT.200.)TN(I,K)=T2d(I,K)
	374	IF(QO(I,K).LT.1.E-08)QO(I,K)=1.E-08
	375	ENDDO
	376	ENDDO
	377	ens4n=0
	378	nbegin=0
	379	nend=0
	380	if(ens4_spread.gt.1)then
	381	nbegin=-ens4_spread/2
	382	nend=ens4_spread/2
	383	endif
	384	do nn=nbegin,nend,1
	385	jss=max(j+nn,jds+0)
	386	jss=min(jss,jde-1)
	387	do n=nbegin,nend,1
	388	ens4n=ens4n+1
	389	DO K=kts,ktf
	390	DO I=ITS,ITF
	391	iss=max(i+n,ids+0)
	392	iss=min(iss,ide-1)
	393	omeg(I,K,ens4n)= -grho(i,k,j)w(iss,k,jss)
	394	! omeg(I,K,ens4n)= -grho(i,k,j)w(i,k,j)
	395	Tx(I,K,ens4n)=t2d(i,k)+RTHFTEN(iss,k,jss)*dt
	396	! Tx(I,K,ens4n)=t2d(i,k)+RTHFTEN(i,k,j)*dt
	397	if(high_resolution.eq.1)Tx(I,K,ens4n)=t2d(i,k)+ave_f_t(iss,k,jss)*dt
	398	IF(Tx(I,K,ens4n).LT.200.)Tx(I,K,ens4n)=T2d(I,K)
	399	Qx(I,K,ens4n)=q2d(i,k)+RQVFTEN(iss,k,jss)*dt
	400	Qx(I,K,ens4n)=q2d(i,k)+RQVFTEN(i,k,j)*dt
	401	if(high_resolution.eq.1)qx(I,K,ens4n)=q2d(i,k)+ave_f_q(iss,k,jss)*dt
	402	IF(Qx(I,K,ens4n).LT.1.E-08)Qx(I,K,ens4n)=1.E-08
	403	enddo
	404	enddo
	405	enddo !n
	406	enddo !nn
	407	do k= kts+1,ktf-1
	408	DO I = its,itf
	409	if((p2d(i,1)-p2d(i,k)).gt.150.and.p2d(i,k).gt.300)then
	410	dp=-.5*(p2d(i,k+1)-p2d(i,k-1))
	411	umean(i)=umean(i)+us(i,k)*dp
	412	vmean(i)=vmean(i)+vs(i,k)*dp
	413	pmean(i)=pmean(i)+dp
	414	endif
	415	enddo
	416	enddo
	417	DO I = its,itf
	418	umean(i)=umean(i)/pmean(i)
	419	vmean(i)=vmean(i)/pmean(i)
	420	direction(i)=(atan2(umean(i),vmean(i))+3.1415926)*57.29578
	421	if(direction(i).gt.360.)direction(i)=direction(i)-360.
	422	ENDDO
	423	do n=1,ens4
	424	DO K=kts,ktf-1
	425	DO I = its,itf
	426	dq=(q2d(i,k+1)-q2d(i,k))
	427	mconv(i,n)=mconv(i,n)+omeg(i,k,n)*dq/g
	428	enddo
	429	ENDDO
	430	ENDDO
	431	do n=1,ens4
	432	DO I = its,itf
	433	if(mconv(i,n).lt.0.)mconv(i,n)=0.
	434	ENDDO
	435	ENDDO
	436	!
	437	!---- CALL CUMULUS PARAMETERIZATION
	438	!
	439	CALL CUP_enss_3d(outqc,j,AAEQ,T2d,Q2d,TER11,subm,TN,QO,PO,PRET,&
	440	P2d,OUTT,OUTQ,DT,itimestep,tkm,PSUR,US,VS,tcrit,iens,tx,qx, &
	441	mconv,massflni,iact_old_gr,omeg,direction,MASSi_FLX, &
	442	maxiens,maxens,maxens2,maxens3,ensdim, &
	443	APRi_GR,APRi_W,APRi_MC,APRi_ST,APRi_AS, &
	444	APRi_CAPMA,APRi_CAPME,APRi_CAPMI,kbcon,ktop,cupclw, &
	445	xfi_ens,pri_ens,XLANDi,gswi,edti_out,subt,subq, &
	446	! ruc lv_p,rv_p,cpd_p,g0_p,ichoice,ipr,jpr, &
	447	xlv,r_v,cp,g,ichoice,ipr,jpr,ens4,high_resolution, &
	448	itf,jtf,ktf, &
	449	its,ite, jts,jte, kts,kte )
	450
	451
	452	if(j.lt.jbegc.or.j.gt.jendc)go to 100
	453	DO I=ibegc,iendc
	454	cuten(i)=0.
	455	if(pret(i).gt.0.)then
	456	cuten(i)=1.
	457	! raincv(i,j)=pret(i)*dt
	458	endif
	459	ENDDO
	460	DO I=ibegc,iendc
	461	DO K=kts,ktf
	462	cugd_ttens(I,K,J)=subt(i,k)cuten(i)sub_spread
	463	cugd_qvtens(I,K,J)=subq(i,k)cuten(i)sub_spread
	464	cugd_tten(I,K,J)=outt(i,k)*cuten(i)
	465	cugd_qvten(I,K,J)=outq(i,k)*cuten(i)
	466	cugd_qcten(I,K,J)=outqc(i,k)*cuten(i)
	467	ENDDO
	468	ENDDO
	469	DO I=ibegc,iendc
	470	if(pret(i).gt.0.)then
	471	raincv(i,j)=pret(i)*dt
	472	pratec(i,j)=pret(i)
	473	rkbcon = kte+kts - kbcon(i)
	474	rktop = kte+kts - ktop(i)
	475	if (ktop(i) > HTOP(i,j)) HTOP(i,j) = ktop(i)+.001
	476	if (kbcon(i) < HBOT(i,j)) HBOT(i,j) = kbcon(i)+.001
	477	endif
	478	ENDDO
	479	DO n= 1,ensdim
	480	DO I= ibegc,iendc
	481	xf_ens(i,j,n)=xfi_ens(i,j,n)
	482	pr_ens(i,j,n)=pri_ens(i,j,n)
	483	ENDDO
	484	ENDDO
	485	DO I= ibegc,iendc
	486	APR_GR(i,j)=apri_gr(i,j)
	487	APR_w(i,j)=apri_w(i,j)
	488	APR_mc(i,j)=apri_mc(i,j)
	489	APR_st(i,j)=apri_st(i,j)
	490	APR_as(i,j)=apri_as(i,j)
	491	APR_capma(i,j)=apri_capma(i,j)
	492	APR_capme(i,j)=apri_capme(i,j)
	493	APR_capmi(i,j)=apri_capmi(i,j)
	494	mass_flux(i,j)=massi_flx(i,j)
	495	edt_out(i,j)=edti_out(i,j)
	496	ENDDO
	497	IF(PRESENT(RQCCUTEN)) THEN
	498	IF ( F_QC ) THEN
	499	DO K=kts,ktf
	500	DO I=ibegc,iendc
	501	RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
	502	IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
	503	IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=0.
	504	ENDDO
	505	ENDDO
	506	ENDIF
	507	ENDIF
	508
	509	!...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
	510
	511	IF(PRESENT(RQICUTEN).AND.PRESENT(RQCCUTEN))THEN
	512	IF (F_QI) THEN
	513	DO K=kts,ktf
	514	DO I=ibegc,iendc
	515	if(t2d(i,k).lt.258.)then
	516	RQICUTEN(I,K,J)=outqc(I,K)*cuten(i)
	517	cugd_qcten(i,k,j)=0.
	518	RQCCUTEN(I,K,J)=0.
	519	IF ( PRESENT( GDC2 ) ) GDC2(I,K,J)=CUPCLW(I,K)*cuten(i)
	520	else
	521	RQICUTEN(I,K,J)=0.
	522	RQCCUTEN(I,K,J)=outqc(I,K)*cuten(i)
	523	IF ( PRESENT( GDC ) ) GDC(I,K,J)=CUPCLW(I,K)*cuten(i)
	524	endif
	525	ENDDO
	526	ENDDO
	527	ENDIF
	528	ENDIF
	529
	530	100 continue
	531
	532	END SUBROUTINE G3DRV
	533
	534	SUBROUTINE CUP_enss_3d(OUTQC,J,AAEQ,T,Q,Z1,sub_mas, &
	535	TN,QO,PO,PRE,P,OUTT,OUTQ,DTIME,ktau,tkmax,PSUR,US,VS, &
	536	TCRIT,iens,tx,qx,mconv,massfln,iact, &
	537	omeg,direction,massflx,maxiens, &
	538	maxens,maxens2,maxens3,ensdim, &
	539	APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
	540	APR_CAPMA,APR_CAPME,APR_CAPMI,kbcon,ktop,cupclw, & !-lxz
	541	xf_ens,pr_ens,xland,gsw,edt_out,subt,subq, &
	542	xl,rv,cp,g,ichoice,ipr,jpr,ens4,high_resolution, &
	543	itf,jtf,ktf, &
	544	its,ite, jts,jte, kts,kte )
	545
	546	IMPLICIT NONE
	547
	548	integer &
	549	,intent (in ) :: &
	550	itf,jtf,ktf,ktau, &
	551	its,ite, jts,jte, kts,kte,ipr,jpr,ens4,high_resolution
	552	integer, intent (in ) :: &
	553	j,ensdim,maxiens,maxens,maxens2,maxens3,ichoice,iens
	554	!
	555	!
	556	!
	557	real, dimension (its:ite,jts:jte,1:ensdim) &
	558	,intent (inout) :: &
	559	massfln,xf_ens,pr_ens
	560	real, dimension (its:ite,jts:jte) &
	561	,intent (inout ) :: &
	562	APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
	563	APR_CAPME,APR_CAPMI,massflx,edt_out
	564	real, dimension (its:ite,jts:jte) &
	565	,intent (in ) :: &
	566	gsw
	567	integer, dimension (its:ite,jts:jte) &
	568	,intent (in ) :: &
	569	iact
	570	! outtem = output temp tendency (per s)
	571	! outq = output q tendency (per s)
	572	! outqc = output qc tendency (per s)
	573	! pre = output precip
	574	real, dimension (its:ite,kts:kte) &
	575	,intent (inout ) :: &
	576	OUTT,OUTQ,OUTQC,subt,subq,sub_mas,cupclw
	577	real, dimension (its:ite) &
	578	,intent (out ) :: &
	579	pre
	580	!+lxz
	581	integer, dimension (its:ite) &
	582	,intent (out ) :: &
	583	kbcon,ktop
	584	!.lxz
	585	!
	586	! basic environmental input includes moisture convergence (mconv)
	587	! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
	588	! convection for this call only and at that particular gridpoint
	589	!
	590	real, dimension (its:ite,kts:kte) &
	591	,intent (in ) :: &
	592	T,PO,P,US,VS,tn
	593	real, dimension (its:ite,kts:kte,1:ens4) &
	594	,intent (inout ) :: &
	595	omeg,tx,qx
	596	real, dimension (its:ite,kts:kte) &
	597	,intent (inout) :: &
	598	Q,QO
	599	real, dimension (its:ite) &
	600	,intent (in ) :: &
	601	Z1,PSUR,AAEQ,direction,tkmax,xland
	602	real, dimension (its:ite,1:ens4) &
	603	,intent (in ) :: &
	604	mconv
	605
	606
	607	real &
	608	,intent (in ) :: &
	609	dtime,tcrit,xl,cp,rv,g
	610
	611
	612	!
	613	! local ensemble dependent variables in this routine
	614	!
	615	real, dimension (its:ite,1:maxens) :: &
	616	xaa0_ens
	617	real, dimension (1:maxens) :: &
	618	mbdt_ens
	619	real, dimension (1:maxens2) :: &
	620	edt_ens
	621	real, dimension (its:ite,1:maxens2) :: &
	622	edtc
	623	real, dimension (its:ite,kts:kte,1:maxens2) :: &
	624	dellat_ens,dellaqc_ens,dellaq_ens,pwo_ens,subt_ens,subq_ens
	625	!
	626	!
	627	!
	628	!***************** the following are your basic environmental
	629	! variables. They carry a "_cup" if they are
	630	! on model cloud levels (staggered). They carry
	631	! an "o"-ending (z becomes zo), if they are the forced
	632	! variables. They are preceded by x (z becomes xz)
	633	! to indicate modification by some typ of cloud
	634	!
	635	! z = heights of model levels
	636	! q = environmental mixing ratio
	637	! qes = environmental saturation mixing ratio
	638	! t = environmental temp
	639	! p = environmental pressure
	640	! he = environmental moist static energy
	641	! hes = environmental saturation moist static energy
	642	! z_cup = heights of model cloud levels
	643	! q_cup = environmental q on model cloud levels
	644	! qes_cup = saturation q on model cloud levels
	645	! t_cup = temperature (Kelvin) on model cloud levels
	646	! p_cup = environmental pressure
	647	! he_cup = moist static energy on model cloud levels
	648	! hes_cup = saturation moist static energy on model cloud levels
	649	! gamma_cup = gamma on model cloud levels
	650	!
	651	!
	652	! hcd = moist static energy in downdraft
	653	! zd normalized downdraft mass flux
	654	! dby = buoancy term
	655	! entr = entrainment rate
	656	! zd = downdraft normalized mass flux
	657	! entr= entrainment rate
	658	! hcd = h in model cloud
	659	! bu = buoancy term
	660	! zd = normalized downdraft mass flux
	661	! gamma_cup = gamma on model cloud levels
	662	! mentr_rate = entrainment rate
	663	! qcd = cloud q (including liquid water) after entrainment
	664	! qrch = saturation q in cloud
	665	! pwd = evaporate at that level
	666	! pwev = total normalized integrated evaoprate (I2)
	667	! entr= entrainment rate
	668	! z1 = terrain elevation
	669	! entr = downdraft entrainment rate
	670	! jmin = downdraft originating level
	671	! kdet = level above ground where downdraft start detraining
	672	! psur = surface pressure
	673	! z1 = terrain elevation
	674	! pr_ens = precipitation ensemble
	675	! xf_ens = mass flux ensembles
	676	! massfln = downdraft mass flux ensembles used in next timestep
	677	! omeg = omega from large scale model
	678	! mconv = moisture convergence from large scale model
	679	! zd = downdraft normalized mass flux
	680	! zu = updraft normalized mass flux
	681	! dir = "storm motion"
	682	! mbdt = arbitrary numerical parameter
	683	! dtime = dt over which forcing is applied
	684	! iact_gr_old = flag to tell where convection was active
	685	! kbcon = LFC of parcel from k22
	686	! k22 = updraft originating level
	687	! icoic = flag if only want one closure (usually set to zero!)
	688	! dby = buoancy term
	689	! ktop = cloud top (output)
	690	! xmb = total base mass flux
	691	! hc = cloud moist static energy
	692	! hkb = moist static energy at originating level
	693	! mentr_rate = entrainment rate
	694
	695	real, dimension (its:ite,kts:kte) :: &
	696	he,hes,qes,z, &
	697	heo,heso,qeso,zo, &
	698	xhe,xhes,xqes,xz,xt,xq, &
	699
	700	qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
	701	qeso_cup,qo_cup,heo_cup,heso_cup,zo_cup,po_cup,gammao_cup, &
	702	tn_cup, &
	703	xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup,xp_cup,xgamma_cup, &
	704	xt_cup, &
	705
	706	dby,qc,qrcd,pwd,pw,hcd,qcd,dbyd,hc,qrc,zu,zd,clw_all, &
	707	dbyo,qco,qrcdo,pwdo,pwo,hcdo,qcdo,dbydo,hco,qrco,zuo,zdo, &
	708	xdby,xqc,xqrcd,xpwd,xpw,xhcd,xqcd,xhc,xqrc,xzu,xzd, &
	709
	710	! cd = detrainment function for updraft
	711	! cdd = detrainment function for downdraft
	712	! dellat = change of temperature per unit mass flux of cloud ensemble
	713	! dellaq = change of q per unit mass flux of cloud ensemble
	714	! dellaqc = change of qc per unit mass flux of cloud ensemble
	715
	716	cd,cdd,scr1,DELLAH,DELLAQ,DELLAT,DELLAQC,dsubt,dsubq
	717
	718	! aa0 cloud work function for downdraft
	719	! edt = epsilon
	720	! aa0 = cloud work function without forcing effects
	721	! aa1 = cloud work function with forcing effects
	722	! xaa0 = cloud work function with cloud effects (ensemble dependent)
	723	! edt = epsilon
	724	real, dimension (its:ite) :: &
	725	edt,edto,edtx,AA1,AA0,XAA0,HKB,HKBO,aad,XHKB,QKB,QKBO, &
	726	XMB,XPWAV,XPWEV,PWAV,PWEV,PWAVO,PWEVO,BU,BUO,cap_max,xland1, &
	727	cap_max_increment,closure_n
	728	real, dimension (its:ite,1:ens4) :: &
	729	axx
	730	integer, dimension (its:ite) :: &
	731	kzdown,KDET,K22,KB,JMIN,kstabi,kstabm,K22x, & !-lxz
	732	KBCONx,KBx,KTOPx,ierr,ierr2,ierr3,KBMAX
	733
	734	integer :: &
	735	nall,iedt,nens,nens3,ki,I,K,KK,iresult
	736	real :: &
	737	day,dz,mbdt,entr_rate,radius,entrd_rate,mentr_rate,mentrd_rate, &
	738	zcutdown,edtmax,edtmin,depth_min,zkbmax,z_detr,zktop, &
	739	massfld,dh,cap_maxs,trash
	740
	741	integer :: jmini
	742	logical :: keep_going
	743
	744
	745
	746	day=86400.
	747	do i=its,itf
	748	closure_n(i)=16.
	749	xland1(i)=1.
	750	if(xland(i).gt.1.5)xland1(i)=0.
	751	! cap_max_increment(i)=50.
	752	cap_max_increment(i)=25.
	753	enddo
	754	!
	755	!--- specify entrainmentrate and detrainmentrate
	756	!
	757	if(iens.le.4)then
	758	radius=14000.-float(iens)*2000.
	759	else
	760	radius=12000.
	761	endif
	762	!
	763	!--- gross entrainment rate (these may be changed later on in the
	764	!--- program, depending what your detrainment is!!)
	765	!
	766	entr_rate=.2/radius
	767	!
	768	!--- entrainment of mass
	769	!
	770	mentrd_rate=0.
	771	mentr_rate=entr_rate
	772	!
	773	!--- initial detrainmentrates
	774	!
	775	do k=kts,ktf
	776	do i=its,itf
	777	cupclw(i,k)=0.
	778	cd(i,k)=0.01*entr_rate
	779	cdd(i,k)=0.
	780	enddo
	781	enddo
	782	!
	783	!--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
	784	! base mass flux
	785	!
	786	edtmax=1.
	787	edtmin=.2
	788	!
	789	!--- minimum depth (m), clouds must have
	790	!
	791	depth_min=500.
	792	!
	793	!--- maximum depth (mb) of capping
	794	!--- inversion (larger cap = no convection)
	795	!
	796	! cap_maxs=125.
	797	cap_maxs=75.
	798	DO i=its,itf
	799	kbmax(i)=1
	800	aa0(i)=0.
	801	aa1(i)=0.
	802	aad(i)=0.
	803	edt(i)=0.
	804	kstabm(i)=ktf-1
	805	IERR(i)=0
	806	IERR2(i)=0
	807	IERR3(i)=0
	808	if(aaeq(i).lt.-0.1)then
	809	ierr(i)=20
	810	endif
	811	enddo
	812	!
	813	!--- first check for upstream convection
	814	!
	815	do i=its,itf
	816	cap_max(i)=cap_maxs
	817	if(gsw(i,j).lt.1.or.high_resolution.eq.1)cap_max(i)=25.
	818	iresult=0
	819
	820	enddo
	821	!
	822	!--- max height(m) above ground where updraft air can originate
	823	!
	824	zkbmax=4000.
	825	!
	826	!--- height(m) above which no downdrafts are allowed to originate
	827	!
	828	zcutdown=3000.
	829	!
	830	!--- depth(m) over which downdraft detrains all its mass
	831	!
	832	z_detr=1250.
	833	!
	834	do nens=1,maxens
	835	mbdt_ens(nens)=(float(nens)-3.)dtime1.e-3+dtime*5.E-03
	836	enddo
	837	do nens=1,maxens2
	838	edt_ens(nens)=.95-float(nens)*.01
	839	enddo
	840	!
	841	!--- environmental conditions, FIRST HEIGHTS
	842	!
	843	do i=its,itf
	844	if(ierr(i).ne.20)then
	845	do k=1,maxensmaxens2maxens3
	846	xf_ens(i,j,(iens-1)maxensmaxens2*maxens3+k)=0.
	847	pr_ens(i,j,(iens-1)maxensmaxens2*maxens3+k)=0.
	848	enddo
	849	endif
	850	enddo
	851	!
	852	!--- calculate moist static energy, heights, qes
	853	!
	854	call cup_env(z,qes,he,hes,t,q,p,z1, &
	855	psur,ierr,tcrit,0,xl,cp, &
	856	itf,jtf,ktf, &
	857	its,ite, jts,jte, kts,kte)
	858	call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
	859	psur,ierr,tcrit,0,xl,cp, &
	860	itf,jtf,ktf, &
	861	its,ite, jts,jte, kts,kte)
	862	!
	863	!--- environmental values on cloud levels
	864	!
	865	call cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup,he_cup, &
	866	hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
	867	ierr,z1,xl,rv,cp, &
	868	itf,jtf,ktf, &
	869	its,ite, jts,jte, kts,kte)
	870	call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
	871	heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur, &
	872	ierr,z1,xl,rv,cp, &
	873	itf,jtf,ktf, &
	874	its,ite, jts,jte, kts,kte)
	875	do i=its,itf
	876	if(ierr(i).eq.0)then
	877	!
	878	do k=kts,ktf
	879	if(zo_cup(i,k).gt.zkbmax+z1(i))then
	880	kbmax(i)=k
	881	go to 25
	882	endif
	883	enddo
	884	25 continue
	885	!
	886	!--- level where detrainment for downdraft starts
	887	!
	888	do k=kts,ktf
	889	if(zo_cup(i,k).gt.z_detr+z1(i))then
	890	kdet(i)=k
	891	go to 26
	892	endif
	893	enddo
	894	26 continue
	895	!
	896	endif
	897	enddo
	898	!
	899	!
	900	!
	901	!------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
	902	!
	903	CALL cup_MAXIMI(HEO_CUP,3,KBMAX,K22,ierr, &
	904	itf,jtf,ktf, &
	905	its,ite, jts,jte, kts,kte)
	906	DO 36 i=its,itf
	907	IF(ierr(I).eq.0.)THEN
	908	IF(K22(I).GE.KBMAX(i))ierr(i)=2
	909	endif
	910	36 CONTINUE
	911	!
	912	!--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
	913	!
	914	call cup_kbcon(cap_max_increment,1,k22,kbcon,heo_cup,heso_cup, &
	915	ierr,kbmax,po_cup,cap_max, &
	916	itf,jtf,ktf, &
	917	its,ite, jts,jte, kts,kte)
	918	!
	919	!--- increase detrainment in stable layers
	920	!
	921	CALL cup_minimi(HEso_cup,Kbcon,kstabm,kstabi,ierr, &
	922	itf,jtf,ktf, &
	923	its,ite, jts,jte, kts,kte)
	924	do i=its,itf
	925	IF(ierr(I).eq.0.)THEN
	926	if(kstabm(i)-1.gt.kstabi(i))then
	927	do k=kstabi(i),kstabm(i)-1
	928	cd(i,k)=cd(i,k-1)+.15*entr_rate
	929	if(cd(i,k).gt.1.0entr_rate)cd(i,k)=1.0entr_rate
	930	enddo
	931	ENDIF
	932	ENDIF
	933	ENDDO
	934	!
	935	!--- calculate incloud moist static energy
	936	!
	937	call cup_up_he(k22,hkb,z_cup,cd,mentr_rate,he_cup,hc, &
	938	kbcon,ierr,dby,he,hes_cup, &
	939	itf,jtf,ktf, &
	940	its,ite, jts,jte, kts,kte)
	941	call cup_up_he(k22,hkbo,zo_cup,cd,mentr_rate,heo_cup,hco, &
	942	kbcon,ierr,dbyo,heo,heso_cup, &
	943	itf,jtf,ktf, &
	944	its,ite, jts,jte, kts,kte)
	945
	946	!--- DETERMINE CLOUD TOP - KTOP
	947	!
	948	call cup_ktop(1,dbyo,kbcon,ktop,ierr, &
	949	itf,jtf,ktf, &
	950	its,ite, jts,jte, kts,kte)
	951	DO 37 i=its,itf
	952	kzdown(i)=0
	953	if(ierr(i).eq.0)then
	954	zktop=(zo_cup(i,ktop(i))-z1(i))*.6
	955	zktop=min(zktop+z1(i),zcutdown+z1(i))
	956	do k=kts,kte
	957	if(zo_cup(i,k).gt.zktop)then
	958	kzdown(i)=k
	959	go to 37
	960	endif
	961	enddo
	962	endif
	963	37 CONTINUE
	964	!
	965	!--- DOWNDRAFT ORIGINATING LEVEL - JMIN
	966	!
	967	call cup_minimi(HEso_cup,K22,kzdown,JMIN,ierr, &
	968	itf,jtf,ktf, &
	969	its,ite, jts,jte, kts,kte)
	970	DO 100 i=its,ite
	971	IF(ierr(I).eq.0.)THEN
	972	!
	973	!--- check whether it would have buoyancy, if there where
	974	!--- no entrainment/detrainment
	975	!
	976	jmini = jmin(i)
	977	keep_going = .TRUE.
	978	do while ( keep_going )
	979	keep_going = .FALSE.
	980	if ( jmini - 1 .lt. kdet(i) ) kdet(i) = jmini-1
	981	if ( jmini .ge. ktop(i)-1 ) jmini = ktop(i) - 2
	982	ki = jmini
	983	hcdo(i,ki)=heso_cup(i,ki)
	984	DZ=Zo_cup(i,Ki+1)-Zo_cup(i,Ki)
	985	dh=0.
	986	do k=ki-1,1,-1
	987	hcdo(i,k)=heso_cup(i,jmini)
	988	DZ=Zo_cup(i,K+1)-Zo_cup(i,K)
	989	dh=dh+dz*(HCDo(i,K)-heso_cup(i,k))
	990	if(dh.gt.0.)then
	991	jmini=jmini-1
	992	if ( jmini .gt. 3 ) then
	993	keep_going = .TRUE.
	994	else
	995	ierr(i) = 9
	996	exit
	997	endif
	998	endif
	999	enddo
	1000	enddo
	1001	jmin(i) = jmini
	1002	if ( jmini .le. 3 ) then
	1003	ierr(i)=4
	1004	endif
	1005	ENDIF
	1006	100 continue
	1007	!
	1008	! - Must have at least depth_min m between cloud convective base
	1009	! and cloud top.
	1010	!
	1011	do i=its,itf
	1012	IF(ierr(I).eq.0.)THEN
	1013	IF(-zo_cup(I,KBCON(I))+zo_cup(I,KTOP(I)).LT.depth_min)then
	1014	ierr(i)=6
	1015	endif
	1016	endif
	1017	enddo
	1018
	1019	!
	1020	!c--- normalized updraft mass flux profile
	1021	!
	1022	call cup_up_nms(zu,z_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
	1023	itf,jtf,ktf, &
	1024	its,ite, jts,jte, kts,kte)
	1025	call cup_up_nms(zuo,zo_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
	1026	itf,jtf,ktf, &
	1027	its,ite, jts,jte, kts,kte)
	1028	!
	1029	!c--- normalized downdraft mass flux profile,also work on bottom detrainment
	1030	!--- in this routine
	1031	!
	1032	call cup_dd_nms(zd,z_cup,cdd,mentrd_rate,jmin,ierr, &
	1033	0,kdet,z1, &
	1034	itf,jtf,ktf, &
	1035	its,ite, jts,jte, kts,kte)
	1036	call cup_dd_nms(zdo,zo_cup,cdd,mentrd_rate,jmin,ierr, &
	1037	1,kdet,z1, &
	1038	itf,jtf,ktf, &
	1039	its,ite, jts,jte, kts,kte)
	1040	!
	1041	!--- downdraft moist static energy
	1042	!
	1043	call cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,mentrd_rate, &
	1044	jmin,ierr,he,dbyd,he_cup, &
	1045	itf,jtf,ktf, &
	1046	its,ite, jts,jte, kts,kte)
	1047	call cup_dd_he(heso_cup,zdo,hcdo,zo_cup,cdd,mentrd_rate, &
	1048	jmin,ierr,heo,dbydo,he_cup,&
	1049	itf,jtf,ktf, &
	1050	its,ite, jts,jte, kts,kte)
	1051	!
	1052	!--- calculate moisture properties of downdraft
	1053	!
	1054	call cup_dd_moisture_3d(zd,hcd,hes_cup,qcd,qes_cup, &
	1055	pwd,q_cup,z_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
	1056	pwev,bu,qrcd,q,he,t_cup,2,xl,high_resolution, &
	1057	itf,jtf,ktf, &
	1058	its,ite, jts,jte, kts,kte)
	1059	call cup_dd_moisture_3d(zdo,hcdo,heso_cup,qcdo,qeso_cup, &
	1060	pwdo,qo_cup,zo_cup,cdd,mentrd_rate,jmin,ierr,gammao_cup, &
	1061	pwevo,bu,qrcdo,qo,heo,tn_cup,1,xl,high_resolution, &
	1062	itf,jtf,ktf, &
	1063	its,ite, jts,jte, kts,kte)
	1064	!
	1065	!--- calculate moisture properties of updraft
	1066	!
	1067	call cup_up_moisture(ierr,z_cup,qc,qrc,pw,pwav, &
	1068	kbcon,ktop,cd,dby,mentr_rate,clw_all, &
	1069	q,GAMMA_cup,zu,qes_cup,k22,q_cup,xl, &
	1070	itf,jtf,ktf, &
	1071	its,ite, jts,jte, kts,kte)
	1072	do k=kts,ktf
	1073	do i=its,itf
	1074	cupclw(i,k)=qrc(i,k)
	1075	enddo
	1076	enddo
	1077	call cup_up_moisture(ierr,zo_cup,qco,qrco,pwo,pwavo, &
	1078	kbcon,ktop,cd,dbyo,mentr_rate,clw_all, &
	1079	qo,GAMMAo_cup,zuo,qeso_cup,k22,qo_cup,xl,&
	1080	itf,jtf,ktf, &
	1081	its,ite, jts,jte, kts,kte)
	1082	!
	1083	!--- calculate workfunctions for updrafts
	1084	!
	1085	call cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
	1086	kbcon,ktop,ierr, &
	1087	itf,jtf,ktf, &
	1088	its,ite, jts,jte, kts,kte)
	1089	call cup_up_aa0(aa1,zo,zuo,dbyo,GAMMAo_CUP,tn_cup, &
	1090	kbcon,ktop,ierr, &
	1091	itf,jtf,ktf, &
	1092	its,ite, jts,jte, kts,kte)
	1093	do i=its,itf
	1094	if(ierr(i).eq.0)then
	1095	if(aa1(i).eq.0.)then
	1096	ierr(i)=17
	1097	endif
	1098	endif
	1099	enddo
	1100	call cup_axx(tcrit,kbmax,z1,p,psur,xl,rv,cp,tx,qx,axx,ierr, &
	1101	cap_max,cap_max_increment,entr_rate,mentr_rate,&
	1102	j,itf,jtf,ktf, &
	1103	its,ite, jts,jte, kts,kte,ens4)
	1104
	1105	!
	1106	!--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
	1107	!
	1108	call cup_dd_edt(ierr,us,vs,zo,ktop,kbcon,edt,po,pwavo, &
	1109	pwevo,edtmax,edtmin,maxens2,edtc, &
	1110	itf,jtf,ktf, &
	1111	its,ite, jts,jte, kts,kte)
	1112	do 250 iedt=1,maxens2
	1113	do i=its,itf
	1114	if(ierr(i).eq.0)then
	1115	edt(i)=edtc(i,iedt)
	1116	edto(i)=edtc(i,iedt)
	1117	edtx(i)=edtc(i,iedt)
	1118	edt_out(i,j)=edtc(i,2)
	1119	if(high_resolution.eq.1)then
	1120	edt(i)=edtc(i,3)
	1121	edto(i)=edtc(i,3)
	1122	edtx(i)=edtc(i,3)
	1123	edt_out(i,j)=edtc(i,3)
	1124	endif
	1125	endif
	1126	enddo
	1127	do k=kts,ktf
	1128	do i=its,itf
	1129	subt_ens(i,k,iedt)=0.
	1130	subq_ens(i,k,iedt)=0.
	1131	dellat_ens(i,k,iedt)=0.
	1132	dellaq_ens(i,k,iedt)=0.
	1133	dellaqc_ens(i,k,iedt)=0.
	1134	pwo_ens(i,k,iedt)=0.
	1135	enddo
	1136	enddo
	1137	!
	1138	if(j.eq.jpr.and.iedt.eq.1.and.ipr.gt.its.and.ipr.lt.ite)then
	1139	! if(j.eq.jpr)then
	1140	i=ipr
	1141	! write(0,*)'in 250 loop ',iedt,edt(ipr),ierr(ipr)
	1142	! if(ierr(i).eq.0.or.ierr(i).eq.3)then
	1143	write(0,*)'250',k22(I),kbcon(i),ktop(i),jmin(i)
	1144	write(0,*)edt(i),aa0(i),aa1(i)
	1145	do k=kts,ktf
	1146	write(0,*)k,z(i,k),he(i,k),hes(i,k)
	1147	enddo
	1148	write(0,*)'end 250 loop ',iedt,edt(ipr),ierr(ipr)
	1149	do k=1,ktop(i)+1
	1150	write(0,*)zu(i,k),zd(i,k),pw(i,k),pwd(i,k)
	1151	enddo
	1152	! endif
	1153	endif
	1154	do i=its,itf
	1155	aad(i)=0.
	1156	enddo
	1157	!
	1158	!--- change per unit mass that a model cloud would modify the environment
	1159	!
	1160	!--- 1. in bottom layer
	1161	!
	1162	call cup_dellabot(ipr,jpr,heo_cup,ierr,zo_cup,po,hcdo,edto, &
	1163	zdo,cdd,heo,dellah,dsubt,j,mentrd_rate,zo,g, &
	1164	itf,jtf,ktf, &
	1165	its,ite, jts,jte, kts,kte)
	1166	call cup_dellabot(ipr,jpr,qo_cup,ierr,zo_cup,po,qrcdo,edto, &
	1167	zdo,cdd,qo,dellaq,dsubq,j,mentrd_rate,zo,g,&
	1168	itf,jtf,ktf, &
	1169	its,ite, jts,jte, kts,kte)
	1170	!
	1171	!--- 2. everywhere else
	1172	!
	1173	call cup_dellas_3d(ierr,zo_cup,po_cup,hcdo,edto,zdo,cdd, &
	1174	heo,dellah,dsubt,j,mentrd_rate,zuo,g, &
	1175	cd,hco,ktop,k22,kbcon,mentr_rate,jmin,heo_cup,kdet, &
	1176	k22,ipr,jpr,'deep',high_resolution, &
	1177	itf,jtf,ktf, &
	1178	its,ite, jts,jte, kts,kte)
	1179	!
	1180	!-- take out cloud liquid water for detrainment
	1181	!
	1182	!?? do k=kts,ktf
	1183	do k=kts,ktf-1
	1184	do i=its,itf
	1185	scr1(i,k)=0.
	1186	dellaqc(i,k)=0.
	1187	if(ierr(i).eq.0)then
	1188	scr1(i,k)=qco(i,k)-qrco(i,k)
	1189	if(k.eq.ktop(i)-0)dellaqc(i,k)= &
	1190	.01zuo(i,ktop(i))qrco(i,ktop(i))* &
	1191	9.81/(po_cup(i,k)-po_cup(i,k+1))
	1192	if(k.lt.ktop(i).and.k.gt.kbcon(i))then
	1193	dz=zo_cup(i,k+1)-zo_cup(i,k)
	1194	dellaqc(i,k)=.019.81cd(i,k)dzzuo(i,k) &
	1195	.5(qrco(i,k)+qrco(i,k+1))/ &
	1196	(po_cup(i,k)-po_cup(i,k+1))
	1197	endif
	1198	endif
	1199	enddo
	1200	enddo
	1201	call cup_dellas_3d(ierr,zo_cup,po_cup,qrcdo,edto,zdo,cdd, &
	1202	qo,dellaq,dsubq,j,mentrd_rate,zuo,g, &
	1203	cd,qco,ktop,k22,kbcon,mentr_rate,jmin,qo_cup,kdet, &
	1204	k22,ipr,jpr,'deep',high_resolution, &
	1205	itf,jtf,ktf, &
	1206	its,ite, jts,jte, kts,kte )
	1207	!
	1208	!--- using dellas, calculate changed environmental profiles
	1209	!
	1210	! do 200 nens=1,maxens
	1211	mbdt=mbdt_ens(2)
	1212	do i=its,itf
	1213	xaa0_ens(i,1)=0.
	1214	xaa0_ens(i,2)=0.
	1215	xaa0_ens(i,3)=0.
	1216	enddo
	1217
	1218	if(j.eq.jpr)then
	1219	write(0,*)'xt',xl,'DELLAH(I,K),DELLAQ(I,K),dsubq(I,K),dsubt(i,k)'
	1220	endif
	1221	do k=kts,ktf
	1222	do i=its,itf
	1223	dellat(i,k)=0.
	1224	if(ierr(i).eq.0)then
	1225	trash=dsubt(i,k)
	1226	XHE(I,K)=(dsubt(i,k)+DELLAH(I,K))*MBDT+HEO(I,K)
	1227	XQ(I,K)=(dsubq(i,k)+DELLAQ(I,K))*MBDT+QO(I,K)
	1228	DELLAT(I,K)=(1./cp)(DELLAH(I,K)-xlDELLAQ(I,K))
	1229	dSUBT(I,K)=(1./cp)(dsubt(i,k)-xldsubq(i,k))
	1230	XT(I,K)= (DELLAT(I,K)+dsubt(i,k))*MBDT+TN(I,K)
	1231	IF(XQ(I,K).LE.0.)XQ(I,K)=1.E-08
	1232	if(i.eq.ipr.and.j.eq.jpr)then
	1233	write(0,*)k,trash,DELLAQ(I,K),dsubq(I,K),dsubt(i,k)
	1234	endif
	1235	ENDIF
	1236	enddo
	1237	enddo
	1238	do i=its,itf
	1239	if(ierr(i).eq.0)then
	1240	XHE(I,ktf)=HEO(I,ktf)
	1241	XQ(I,ktf)=QO(I,ktf)
	1242	XT(I,ktf)=TN(I,ktf)
	1243	IF(XQ(I,ktf).LE.0.)XQ(I,ktf)=1.E-08
	1244	endif
	1245	enddo
	1246	!
	1247	!--- calculate moist static energy, heights, qes
	1248	!
	1249	call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
	1250	psur,ierr,tcrit,2,xl,cp, &
	1251	itf,jtf,ktf, &
	1252	its,ite, jts,jte, kts,kte)
	1253	!
	1254	!--- environmental values on cloud levels
	1255	!
	1256	call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
	1257	xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur, &
	1258	ierr,z1,xl,rv,cp, &
	1259	itf,jtf,ktf, &
	1260	its,ite, jts,jte, kts,kte)
	1261	!
	1262	!
	1263	!**************************** static control
	1264	!
	1265	!--- moist static energy inside cloud
	1266	!
	1267	do i=its,itf
	1268	if(ierr(i).eq.0)then
	1269	xhkb(i)=xhe(i,k22(i))
	1270	endif
	1271	enddo
	1272	call cup_up_he(k22,xhkb,xz_cup,cd,mentr_rate,xhe_cup,xhc, &
	1273	kbcon,ierr,xdby,xhe,xhes_cup, &
	1274	itf,jtf,ktf, &
	1275	its,ite, jts,jte, kts,kte)
	1276	!
	1277	!c--- normalized mass flux profile
	1278	!
	1279	call cup_up_nms(xzu,xz_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
	1280	itf,jtf,ktf, &
	1281	its,ite, jts,jte, kts,kte)
	1282	!
	1283	!--- moisture downdraft
	1284	!
	1285	call cup_dd_nms(xzd,xz_cup,cdd,mentrd_rate,jmin,ierr, &
	1286	1,kdet,z1, &
	1287	itf,jtf,ktf, &
	1288	its,ite, jts,jte, kts,kte)
	1289	call cup_dd_he(xhes_cup,xzd,xhcd,xz_cup,cdd,mentrd_rate, &
	1290	jmin,ierr,xhe,dbyd,xhe_cup,&
	1291	itf,jtf,ktf, &
	1292	its,ite, jts,jte, kts,kte)
	1293	call cup_dd_moisture_3d(xzd,xhcd,xhes_cup,xqcd,xqes_cup, &
	1294	xpwd,xq_cup,xz_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
	1295	xpwev,bu,xqrcd,xq,xhe,xt_cup,3,xl,high_resolution, &
	1296	itf,jtf,ktf, &
	1297	its,ite, jts,jte, kts,kte)
	1298
	1299	!
	1300	!------- MOISTURE updraft
	1301	!
	1302	call cup_up_moisture(ierr,xz_cup,xqc,xqrc,xpw,xpwav, &
	1303	kbcon,ktop,cd,xdby,mentr_rate,clw_all, &
	1304	xq,GAMMA_cup,xzu,xqes_cup,k22,xq_cup,xl, &
	1305	itf,jtf,ktf, &
	1306	its,ite, jts,jte, kts,kte)
	1307	!
	1308	!--- workfunctions for updraft
	1309	!
	1310	call cup_up_aa0(xaa0,xz,xzu,xdby,GAMMA_CUP,xt_cup, &
	1311	kbcon,ktop,ierr, &
	1312	itf,jtf,ktf, &
	1313	its,ite, jts,jte, kts,kte)
	1314	do 200 nens=1,maxens
	1315	do i=its,itf
	1316	if(ierr(i).eq.0)then
	1317	xaa0_ens(i,nens)=xaa0(i)
	1318	nall=(iens-1)maxens3maxens*maxens2 &
	1319	+(iedt-1)maxensmaxens3 &
	1320	+(nens-1)*maxens3
	1321	do k=kts,ktf
	1322	if(k.le.ktop(i))then
	1323	do nens3=1,maxens3
	1324	if(nens3.eq.7)then
	1325	!--- b=0
	1326	pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3) &
	1327	+edto(i)*pwdo(i,k) &
	1328	+pwo(i,k)
	1329	!--- b=beta
	1330	else if(nens3.eq.8)then
	1331	pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
	1332	pwo(i,k)
	1333	!--- b=beta/2
	1334	else if(nens3.eq.9)then
	1335	pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3) &
	1336	+.5edto(i)pwdo(i,k) &
	1337	+ pwo(i,k)
	1338	else
	1339	pr_ens(i,j,nall+nens3)=pr_ens(i,j,nall+nens3)+ &
	1340	pwo(i,k)+edto(i)*pwdo(i,k)
	1341	endif
	1342	enddo
	1343	endif
	1344	enddo
	1345	if(pr_ens(i,j,nall+7).lt.1.e-6)then
	1346	ierr(i)=18
	1347	do nens3=1,maxens3
	1348	pr_ens(i,j,nall+nens3)=0.
	1349	enddo
	1350	endif
	1351	do nens3=1,maxens3
	1352	if(pr_ens(i,j,nall+nens3).lt.1.e-4)then
	1353	pr_ens(i,j,nall+nens3)=0.
	1354	endif
	1355	enddo
	1356	endif
	1357	enddo
	1358	200 continue
	1359	!
	1360	!--- LARGE SCALE FORCING
	1361	!
	1362	!
	1363	!------- CHECK wether aa0 should have been zero
	1364	!
	1365	!
	1366	CALL cup_MAXIMI(HEO_CUP,3,KBMAX,K22x,ierr, &
	1367	itf,jtf,ktf, &
	1368	its,ite, jts,jte, kts,kte)
	1369	do i=its,itf
	1370	ierr2(i)=ierr(i)
	1371	ierr3(i)=ierr(i)
	1372	enddo
	1373	call cup_kbcon(cap_max_increment,2,k22x,kbconx,heo_cup, &
	1374	heso_cup,ierr2,kbmax,po_cup,cap_max, &
	1375	itf,jtf,ktf, &
	1376	its,ite, jts,jte, kts,kte)
	1377	call cup_kbcon(cap_max_increment,3,k22x,kbconx,heo_cup, &
	1378	heso_cup,ierr3,kbmax,po_cup,cap_max, &
	1379	itf,jtf,ktf, &
	1380	its,ite, jts,jte, kts,kte)
	1381	!
	1382	!--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
	1383	!
	1384
	1385	call cup_forcing_ens_3d(closure_n,xland1,aa0,aa1,xaa0_ens,mbdt_ens,dtime, &
	1386	ierr,ierr2,ierr3,xf_ens,j,'deeps',axx, &
	1387	maxens,iens,iedt,maxens2,maxens3,mconv, &
	1388	po_cup,ktop,omeg,zdo,k22,zuo,pr_ens,edto,kbcon, &
	1389	massflx,iact,direction,ensdim,massfln,ichoice,edt_out, &
	1390	high_resolution,itf,jtf,ktf, &
	1391	its,ite, jts,jte, kts,kte,ens4,ktau)
	1392	!
	1393	do k=kts,ktf
	1394	do i=its,itf
	1395	if(ierr(i).eq.0)then
	1396	subt_ens(i,k,iedt)=dsubt(i,k)
	1397	subq_ens(i,k,iedt)=dsubq(i,k)
	1398	dellat_ens(i,k,iedt)=dellat(i,k)
	1399	dellaq_ens(i,k,iedt)=dellaq(i,k)
	1400	dellaqc_ens(i,k,iedt)=dellaqc(i,k)
	1401	pwo_ens(i,k,iedt)=pwo(i,k)+edt(i)*pwdo(i,k)
	1402	else
	1403	subt_ens(i,k,iedt)=0.
	1404	subq_ens(i,k,iedt)=0.
	1405	dellat_ens(i,k,iedt)=0.
	1406	dellaq_ens(i,k,iedt)=0.
	1407	dellaqc_ens(i,k,iedt)=0.
	1408	pwo_ens(i,k,iedt)=0.
	1409	endif
	1410	if(i.eq.ipr.and.j.eq.jpr)then
	1411	write(0,*)'1',iens,iedt,dellat(i,k),dellat_ens(i,k,iedt), &
	1412	dellaq(i,k), dellaqc(i,k)
	1413	write(0,*)'2',k,subt_ens(i,k,iedt),subq_ens(i,k,iedt)
	1414	endif
	1415	enddo
	1416	enddo
	1417	250 continue
	1418	!
	1419	!--- FEEDBACK
	1420	!
	1421	call cup_output_ens_3d(xf_ens,ierr,dellat_ens,dellaq_ens, &
	1422	dellaqc_ens,subt_ens,subq_ens,subt,subq,outt, &
	1423	outq,outqc,zuo,sub_mas,pre,pwo_ens,xmb,ktop, &
	1424	j,'deep',maxens2,maxens,iens,ierr2,ierr3, &
	1425	pr_ens,maxens3,ensdim,massfln, &
	1426	APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
	1427	APR_CAPMA,APR_CAPME,APR_CAPMI,closure_n,xland1, &
	1428	itf,jtf,ktf, &
	1429	its,ite, jts,jte, kts,kte)
	1430	k=1
	1431	do i=its,itf
	1432	PRE(I)=MAX(PRE(I),0.)
	1433	if(i.eq.ipr.and.j.eq.jpr)then
	1434	write(0,*)'i,j,pre(i),aa0(i),aa1(i)'
	1435	write(0,*)i,j,pre(i),aa0(i)
	1436	endif
	1437	enddo
	1438	!
	1439	!---------------------------done------------------------------
	1440	!
	1441	do i=its,itf
	1442	if(ierr(i).eq.0)then
	1443	if(i.eq.ipr.and.j.eq.jpr)then
	1444	write(0,*)'on output, pre =',pre(i),its,itf,kts,ktf
	1445	do k=kts,ktf
	1446	write(0,)z(i,k),outt(i,k)86400.,subt(i,k)*86400.
	1447	enddo
	1448	write(0,*)i,j,(axx(i,k),k=1,ens4)
	1449	endif
	1450	endif
	1451	enddo
	1452	! print *,'ierr(i) = ',ierr(i),pre(i)
	1453
	1454	END SUBROUTINE CUP_enss_3d
	1455
	1456
	1457	SUBROUTINE cup_dd_aa0(edt,ierr,aa0,jmin,gamma_cup,t_cup, &
	1458	hcd,hes_cup,z,zd, &
	1459	itf,jtf,ktf, &
	1460	its,ite, jts,jte, kts,kte )
	1461
	1462	IMPLICIT NONE
	1463	!
	1464	! on input
	1465	!
	1466
	1467	! only local wrf dimensions are need as of now in this routine
	1468
	1469	integer &
	1470	,intent (in ) :: &
	1471	itf,jtf,ktf, &
	1472	its,ite, jts,jte, kts,kte
	1473	! aa0 cloud work function for downdraft
	1474	! gamma_cup = gamma on model cloud levels
	1475	! t_cup = temperature (Kelvin) on model cloud levels
	1476	! hes_cup = saturation moist static energy on model cloud levels
	1477	! hcd = moist static energy in downdraft
	1478	! edt = epsilon
	1479	! zd normalized downdraft mass flux
	1480	! z = heights of model levels
	1481	! ierr error value, maybe modified in this routine
	1482	!
	1483	real, dimension (its:ite,kts:kte) &
	1484	,intent (in ) :: &
	1485	z,zd,gamma_cup,t_cup,hes_cup,hcd
	1486	real, dimension (its:ite) &
	1487	,intent (in ) :: &
	1488	edt
	1489	integer, dimension (its:ite) &
	1490	,intent (in ) :: &
	1491	jmin
	1492	!
	1493	! input and output
	1494	!
	1495
	1496
	1497	integer, dimension (its:ite) &
	1498	,intent (inout) :: &
	1499	ierr
	1500	real, dimension (its:ite) &
	1501	,intent (out ) :: &
	1502	aa0
	1503	!
	1504	! local variables in this routine
	1505	!
	1506
	1507	integer :: &
	1508	i,k,kk
	1509	real :: &
	1510	dz
	1511	!
	1512	do i=its,itf
	1513	aa0(i)=0.
	1514	enddo
	1515	!
	1516	!?? DO k=kts,kte-1
	1517	DO k=kts,ktf-1
	1518	do i=its,itf
	1519	IF(ierr(I).eq.0.and.k.lt.jmin(i))then
	1520	KK=JMIN(I)-K
	1521	!
	1522	!--- ORIGINAL
	1523	!
	1524	DZ=(Z(I,KK)-Z(I,KK+1))
	1525	AA0(I)=AA0(I)+zd(i,kk)EDT(I)DZ(9.81/(1004.T_cup(I,KK))) &
	1526	*((hcd(i,kk)-hes_cup(i,kk))/(1.+GAMMA_cup(i,kk)))
	1527	endif
	1528	enddo
	1529	enddo
	1530
	1531	END SUBROUTINE CUP_dd_aa0
	1532
	1533
	1534	SUBROUTINE cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
	1535	pwev,edtmax,edtmin,maxens2,edtc, &
	1536	itf,jtf,ktf, &
	1537	its,ite, jts,jte, kts,kte )
	1538
	1539	IMPLICIT NONE
	1540
	1541	integer &
	1542	,intent (in ) :: &
	1543	itf,jtf,ktf, &
	1544	its,ite, jts,jte, kts,kte
	1545	integer, intent (in ) :: &
	1546	maxens2
	1547	!
	1548	! ierr error value, maybe modified in this routine
	1549	!
	1550	real, dimension (its:ite,kts:kte) &
	1551	,intent (in ) :: &
	1552	us,vs,z,p
	1553	real, dimension (its:ite,1:maxens2) &
	1554	,intent (out ) :: &
	1555	edtc
	1556	real, dimension (its:ite) &
	1557	,intent (out ) :: &
	1558	edt
	1559	real, dimension (its:ite) &
	1560	,intent (in ) :: &
	1561	pwav,pwev
	1562	real &
	1563	,intent (in ) :: &
	1564	edtmax,edtmin
	1565	integer, dimension (its:ite) &
	1566	,intent (in ) :: &
	1567	ktop,kbcon
	1568	integer, dimension (its:ite) &
	1569	,intent (inout) :: &
	1570	ierr
	1571	!
	1572	! local variables in this routine
	1573	!
	1574
	1575	integer i,k,kk
	1576	real einc,pef,pefb,prezk,zkbc
	1577	real, dimension (its:ite) :: &
	1578	vshear,sdp,vws
	1579
	1580	!
	1581	!--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
	1582	!
	1583	! */ calculate an average wind shear over the depth of the cloud
	1584	!
	1585	do i=its,itf
	1586	edt(i)=0.
	1587	vws(i)=0.
	1588	sdp(i)=0.
	1589	vshear(i)=0.
	1590	enddo
	1591	do k=1,maxens2
	1592	do i=its,itf
	1593	edtc(i,k)=0.
	1594	enddo
	1595	enddo
	1596	do kk = kts,ktf-1
	1597	do 62 i=its,itf
	1598	IF(ierr(i).ne.0)GO TO 62
	1599	if (kk .le. min0(ktop(i),ktf) .and. kk .ge. kbcon(i)) then
	1600	vws(i) = vws(i)+ &
	1601	(abs((us(i,kk+1)-us(i,kk))/(z(i,kk+1)-z(i,kk))) &
	1602	+ abs((vs(i,kk+1)-vs(i,kk))/(z(i,kk+1)-z(i,kk)))) * &
	1603	(p(i,kk) - p(i,kk+1))
	1604	sdp(i) = sdp(i) + p(i,kk) - p(i,kk+1)
	1605	endif
	1606	if (kk .eq. ktf)vshear(i) = 1.e3 * vws(i) / sdp(i)
	1607	62 continue
	1608	end do
	1609	do i=its,itf
	1610	IF(ierr(i).eq.0)then
	1611	pef=(1.591-.639VSHEAR(I)+.0953(VSHEAR(I)**2) &
	1612	-.00496(VSHEAR(I)*3))
	1613	if(pef.gt.1.)pef=1.
	1614	if(pef.lt.0.)pef=0.
	1615	!
	1616	!--- cloud base precip efficiency
	1617	!
	1618	zkbc=z(i,kbcon(i))*3.281e-3
	1619	prezk=.02
	1620	if(zkbc.gt.3.)then
	1621	prezk=.96729352+zkbc(-.70034167+zkbc(.162179896+zkbc &
	1622	(- 1.2569798E-2+zkbc(4.2772E-4-zkbc*5.44E-6))))
	1623	endif
	1624	if(zkbc.gt.25)then
	1625	prezk=2.4
	1626	endif
	1627	pefb=1./(1.+prezk)
	1628	if(pefb.gt.1.)pefb=1.
	1629	if(pefb.lt.0.)pefb=0.
	1630	EDT(I)=1.-.5*(pefb+pef)
	1631	!--- edt here is 1-precipeff!
	1632	einc=.2*edt(i)
	1633	do k=1,maxens2
	1634	edtc(i,k)=edt(i)+float(k-2)*einc
	1635	enddo
	1636	endif
	1637	enddo
	1638	do i=its,itf
	1639	IF(ierr(i).eq.0)then
	1640	do k=1,maxens2
	1641	EDTC(I,K)=-EDTC(I,K)*PWAV(I)/PWEV(I)
	1642	IF(EDTC(I,K).GT.edtmax)EDTC(I,K)=edtmax
	1643	IF(EDTC(I,K).LT.edtmin)EDTC(I,K)=edtmin
	1644	enddo
	1645	endif
	1646	enddo
	1647
	1648	END SUBROUTINE cup_dd_edt
	1649
	1650
	1651	SUBROUTINE cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,entr, &
	1652	jmin,ierr,he,dby,he_cup, &
	1653	itf,jtf,ktf, &
	1654	its,ite, jts,jte, kts,kte )
	1655
	1656	IMPLICIT NONE
	1657	!
	1658	! on input
	1659	!
	1660
	1661	! only local wrf dimensions are need as of now in this routine
	1662
	1663	integer &
	1664	,intent (in ) :: &
	1665	itf,jtf,ktf, &
	1666	its,ite, jts,jte, kts,kte
	1667	! hcd = downdraft moist static energy
	1668	! he = moist static energy on model levels
	1669	! he_cup = moist static energy on model cloud levels
	1670	! hes_cup = saturation moist static energy on model cloud levels
	1671	! dby = buoancy term
	1672	! cdd= detrainment function
	1673	! z_cup = heights of model cloud levels
	1674	! entr = entrainment rate
	1675	! zd = downdraft normalized mass flux
	1676	!
	1677	real, dimension (its:ite,kts:kte) &
	1678	,intent (in ) :: &
	1679	he,he_cup,hes_cup,z_cup,cdd,zd
	1680	! entr= entrainment rate
	1681	real &
	1682	,intent (in ) :: &
	1683	entr
	1684	integer, dimension (its:ite) &
	1685	,intent (in ) :: &
	1686	jmin
	1687	!
	1688	! input and output
	1689	!
	1690
	1691	! ierr error value, maybe modified in this routine
	1692
	1693	integer, dimension (its:ite) &
	1694	,intent (inout) :: &
	1695	ierr
	1696
	1697	real, dimension (its:ite,kts:kte) &
	1698	,intent (out ) :: &
	1699	hcd,dby
	1700	!
	1701	! local variables in this routine
	1702	!
	1703
	1704	integer :: &
	1705	i,k,ki
	1706	real :: &
	1707	dz
	1708
	1709
	1710	do k=kts+1,ktf
	1711	do i=its,itf
	1712	dby(i,k)=0.
	1713	IF(ierr(I).eq.0)then
	1714	hcd(i,k)=hes_cup(i,k)
	1715	endif
	1716	enddo
	1717	enddo
	1718	!
	1719	do 100 i=its,itf
	1720	IF(ierr(I).eq.0)then
	1721	k=jmin(i)
	1722	hcd(i,k)=hes_cup(i,k)
	1723	dby(i,k)=hcd(i,jmin(i))-hes_cup(i,k)
	1724	!
	1725	do ki=jmin(i)-1,1,-1
	1726	DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
	1727	HCD(i,Ki)=(HCD(i,Ki+1)(1.-.5CDD(i,Ki)*DZ) &
	1728	+entrDZHE(i,Ki) &
	1729	)/(1.+entrDZ-.5CDD(i,Ki)*DZ)
	1730	dby(i,ki)=HCD(i,Ki)-hes_cup(i,ki)
	1731	enddo
	1732	!
	1733	endif
	1734	!--- end loop over i
	1735	100 continue
	1736
	1737
	1738	END SUBROUTINE cup_dd_he
	1739
	1740
	1741	SUBROUTINE cup_dd_moisture_3d(zd,hcd,hes_cup,qcd,qes_cup, &
	1742	pwd,q_cup,z_cup,cdd,entr,jmin,ierr, &
	1743	gamma_cup,pwev,bu,qrcd, &
	1744	q,he,t_cup,iloop,xl,high_resolution, &
	1745	itf,jtf,ktf, &
	1746	its,ite, jts,jte, kts,kte )
	1747
	1748	IMPLICIT NONE
	1749
	1750	integer &
	1751	,intent (in ) :: &
	1752	itf,jtf,ktf, &
	1753	its,ite, jts,jte, kts,kte,high_resolution
	1754	! cdd= detrainment function
	1755	! q = environmental q on model levels
	1756	! q_cup = environmental q on model cloud levels
	1757	! qes_cup = saturation q on model cloud levels
	1758	! hes_cup = saturation h on model cloud levels
	1759	! hcd = h in model cloud
	1760	! bu = buoancy term
	1761	! zd = normalized downdraft mass flux
	1762	! gamma_cup = gamma on model cloud levels
	1763	! mentr_rate = entrainment rate
	1764	! qcd = cloud q (including liquid water) after entrainment
	1765	! qrch = saturation q in cloud
	1766	! pwd = evaporate at that level
	1767	! pwev = total normalized integrated evaoprate (I2)
	1768	! entr= entrainment rate
	1769	!
	1770	real, dimension (its:ite,kts:kte) &
	1771	,intent (in ) :: &
	1772	zd,t_cup,hes_cup,hcd,qes_cup,q_cup,z_cup,cdd,gamma_cup,q,he
	1773	real &
	1774	,intent (in ) :: &
	1775	entr,xl
	1776	integer &
	1777	,intent (in ) :: &
	1778	iloop
	1779	integer, dimension (its:ite) &
	1780	,intent (in ) :: &
	1781	jmin
	1782	integer, dimension (its:ite) &
	1783	,intent (inout) :: &
	1784	ierr
	1785	real, dimension (its:ite,kts:kte) &
	1786	,intent (out ) :: &
	1787	qcd,qrcd,pwd
	1788	real, dimension (its:ite) &
	1789	,intent (out ) :: &
	1790	pwev,bu
	1791	!
	1792	! local variables in this routine
	1793	!
	1794
	1795	integer :: &
	1796	i,k,ki
	1797	real :: &
	1798	dh,dz,dqeva
	1799
	1800	do i=its,itf
	1801	bu(i)=0.
	1802	pwev(i)=0.
	1803	enddo
	1804	do k=kts,ktf
	1805	do i=its,itf
	1806	qcd(i,k)=0.
	1807	qrcd(i,k)=0.
	1808	pwd(i,k)=0.
	1809	enddo
	1810	enddo
	1811	!
	1812	!
	1813	!
	1814	do 100 i=its,itf
	1815	IF(ierr(I).eq.0)then
	1816	k=jmin(i)
	1817	DZ=Z_cup(i,K+1)-Z_cup(i,K)
	1818	qcd(i,k)=q_cup(i,k)
	1819	if(high_resolution.eq.1)qcd(i,k)=.5*(qes_cup(i,k)+q_cup(i,k))
	1820	qrcd(i,k)=qes_cup(i,k)
	1821	pwd(i,jmin(i))=min(0.,qcd(i,k)-qrcd(i,k))
	1822	pwev(i)=pwev(i)+pwd(i,jmin(i))
	1823	qcd(i,k)=qes_cup(i,k)
	1824	!
	1825	DH=HCD(I,k)-HES_cup(I,K)
	1826	bu(i)=dz*dh
	1827	do ki=jmin(i)-1,1,-1
	1828	DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
	1829	QCD(i,Ki)=(qCD(i,Ki+1)(1.-.5CDD(i,Ki)*DZ) &
	1830	+entrDZq(i,Ki) &
	1831	)/(1.+entrDZ-.5CDD(i,Ki)*DZ)
	1832	!
	1833	!--- to be negatively buoyant, hcd should be smaller than hes!
	1834	!
	1835	DH=HCD(I,ki)-HES_cup(I,Ki)
	1836	bu(i)=bu(i)+dz*dh
	1837	QRCD(I,Ki)=qes_cup(i,ki)+(1./XL)*(GAMMA_cup(i,ki) &
	1838	/(1.+GAMMA_cup(i,ki)))*DH
	1839	dqeva=qcd(i,ki)-qrcd(i,ki)
	1840	if(dqeva.gt.0.)dqeva=0.
	1841	pwd(i,ki)=zd(i,ki)*dqeva
	1842	qcd(i,ki)=qrcd(i,ki)
	1843	pwev(i)=pwev(i)+pwd(i,ki)
	1844	! if(iloop.eq.1.and.i.eq.102.and.j.eq.62)then
	1845	! print *,'in cup_dd_moi ', hcd(i,ki),HES_cup(I,Ki),dh,dqeva
	1846	! endif
	1847	enddo
	1848	!
	1849	!--- end loop over i
	1850	if(pwev(I).eq.0.and.iloop.eq.1)then
	1851	! print *,'problem with buoy in cup_dd_moisture',i
	1852	ierr(i)=7
	1853	endif
	1854	if(BU(I).GE.0.and.iloop.eq.1)then
	1855	! print *,'problem with buoy in cup_dd_moisture',i
	1856	ierr(i)=7
	1857	endif
	1858	endif
	1859	100 continue
	1860
	1861	END SUBROUTINE cup_dd_moisture_3d
	1862
	1863
	1864	SUBROUTINE cup_dd_nms(zd,z_cup,cdd,entr,jmin,ierr, &
	1865	itest,kdet,z1, &
	1866	itf,jtf,ktf, &
	1867	its,ite, jts,jte, kts,kte )
	1868
	1869	IMPLICIT NONE
	1870	!
	1871	! on input
	1872	!
	1873
	1874	! only local wrf dimensions are need as of now in this routine
	1875
	1876	integer &
	1877	,intent (in ) :: &
	1878	itf,jtf,ktf, &
	1879	its,ite, jts,jte, kts,kte
	1880	! z_cup = height of cloud model level
	1881	! z1 = terrain elevation
	1882	! entr = downdraft entrainment rate
	1883	! jmin = downdraft originating level
	1884	! kdet = level above ground where downdraft start detraining
	1885	! itest = flag to whether to calculate cdd
	1886
	1887	real, dimension (its:ite,kts:kte) &
	1888	,intent (in ) :: &
	1889	z_cup
	1890	real, dimension (its:ite) &
	1891	,intent (in ) :: &
	1892	z1
	1893	real &
	1894	,intent (in ) :: &
	1895	entr
	1896	integer, dimension (its:ite) &
	1897	,intent (in ) :: &
	1898	jmin,kdet
	1899	integer &
	1900	,intent (in ) :: &
	1901	itest
	1902	!
	1903	! input and output
	1904	!
	1905
	1906	! ierr error value, maybe modified in this routine
	1907
	1908	integer, dimension (its:ite) &
	1909	,intent (inout) :: &
	1910	ierr
	1911	! zd is the normalized downdraft mass flux
	1912	! cdd is the downdraft detrainmen function
	1913
	1914	real, dimension (its:ite,kts:kte) &
	1915	,intent (out ) :: &
	1916	zd
	1917	real, dimension (its:ite,kts:kte) &
	1918	,intent (inout) :: &
	1919	cdd
	1920	!
	1921	! local variables in this routine
	1922	!
	1923
	1924	integer :: &
	1925	i,k,ki
	1926	real :: &
	1927	a,perc,dz
	1928
	1929	!
	1930	!--- perc is the percentage of mass left when hitting the ground
	1931	!
	1932	perc=.03
	1933
	1934	do k=kts,ktf
	1935	do i=its,itf
	1936	zd(i,k)=0.
	1937	if(itest.eq.0)cdd(i,k)=0.
	1938	enddo
	1939	enddo
	1940	a=1.-perc
	1941	!
	1942	!
	1943	!
	1944	do 100 i=its,itf
	1945	IF(ierr(I).eq.0)then
	1946	zd(i,jmin(i))=1.
	1947	!
	1948	!--- integrate downward, specify detrainment(cdd)!
	1949	!
	1950	do ki=jmin(i)-1,1,-1
	1951	DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
	1952	if(ki.le.kdet(i).and.itest.eq.0)then
	1953	cdd(i,ki)=entr+(1.- (a*(z_cup(i,ki)-z1(i)) &
	1954	+perc*(z_cup(i,kdet(i))-z1(i)) ) &
	1955	/(a*(z_cup(i,ki+1)-z1(i)) &
	1956	+perc*(z_cup(i,kdet(i))-z1(i))))/dz
	1957	endif
	1958	zd(i,ki)=zd(i,ki+1)(1.+(entr-cdd(i,ki))dz)
	1959	enddo
	1960	!
	1961	endif
	1962	!--- end loop over i
	1963	100 continue
	1964
	1965	END SUBROUTINE cup_dd_nms
	1966
	1967
	1968	SUBROUTINE cup_dellabot(ipr,jpr,he_cup,ierr,z_cup,p_cup, &
	1969	hcd,edt,zd,cdd,he,della,subs,j,mentrd_rate,z,g, &
	1970	itf,jtf,ktf, &
	1971	its,ite, jts,jte, kts,kte )
	1972
	1973	IMPLICIT NONE
	1974
	1975	integer &
	1976	,intent (in ) :: &
	1977	itf,jtf,ktf, &
	1978	its,ite, jts,jte, kts,kte
	1979	integer, intent (in ) :: &
	1980	j,ipr,jpr
	1981	!
	1982	! ierr error value, maybe modified in this routine
	1983	!
	1984	real, dimension (its:ite,kts:kte) &
	1985	,intent (out ) :: &
	1986	della,subs
	1987	real, dimension (its:ite,kts:kte) &
	1988	,intent (in ) :: &
	1989	z_cup,p_cup,hcd,zd,cdd,he,z,he_cup
	1990	real, dimension (its:ite) &
	1991	,intent (in ) :: &
	1992	edt
	1993	real &
	1994	,intent (in ) :: &
	1995	g,mentrd_rate
	1996	integer, dimension (its:ite) &
	1997	,intent (inout) :: &
	1998	ierr
	1999	!
	2000	! local variables in this routine
	2001	!
	2002
	2003	integer i
	2004	real detdo,detdo1,detdo2,entdo,dp,dz,subin, &
	2005	totmas
	2006	!
	2007	!
	2008	do 100 i=its,itf
	2009	della(i,1)=0.
	2010	subs(i,1)=0.
	2011	if(ierr(i).ne.0)go to 100
	2012	dz=z_cup(i,2)-z_cup(i,1)
	2013	DP=100.*(p_cup(i,1)-P_cup(i,2))
	2014	detdo1=edt(i)zd(i,2)CDD(i,1)*DZ
	2015	detdo2=edt(i)*zd(i,1)
	2016	entdo=edt(i)zd(i,2)mentrd_rate*dz
	2017	subin=-EDT(I)*zd(i,2)
	2018	detdo=detdo1+detdo2-entdo+subin
	2019	DELLA(I,1)=(detdo1.5(HCD(i,1)+HCD(i,2)) &
	2020	+detdo2*hcd(i,1) &
	2021	+subin*he_cup(i,2) &
	2022	-entdohe(i,1))g/dp
	2023	SUBS(I,1)=0.
	2024	if(i.eq.ipr.and.j.eq.jpr)then
	2025	write(0,*)'db1',della(i,1),subs(i,1),subin,entdo
	2026	write(0,*)'db2',detdo1,detdo2,detdo1+detdo2-entdo+subin
	2027	endif
	2028	100 CONTINUE
	2029
	2030	END SUBROUTINE cup_dellabot
	2031
	2032
	2033	SUBROUTINE cup_dellas_3d(ierr,z_cup,p_cup,hcd,edt,zd,cdd, &
	2034	he,della,subs,j,mentrd_rate,zu,g, &
	2035	cd,hc,ktop,k22,kbcon,mentr_rate,jmin,he_cup,kdet,kpbl, &
	2036	ipr,jpr,name,high_res, &
	2037	itf,jtf,ktf, &
	2038	its,ite, jts,jte, kts,kte )
	2039
	2040	IMPLICIT NONE
	2041
	2042	integer &
	2043	,intent (in ) :: &
	2044	itf,jtf,ktf, &
	2045	its,ite, jts,jte, kts,kte
	2046	integer, intent (in ) :: &
	2047	j,ipr,jpr,high_res
	2048	!
	2049	! ierr error value, maybe modified in this routine
	2050	!
	2051	real, dimension (its:ite,kts:kte) &
	2052	,intent (out ) :: &
	2053	della,subs
	2054	real, dimension (its:ite,kts:kte) &
	2055	,intent (in ) :: &
	2056	z_cup,p_cup,hcd,zd,cdd,he,hc,cd,zu,he_cup
	2057	real, dimension (its:ite) &
	2058	,intent (in ) :: &
	2059	edt
	2060	real &
	2061	,intent (in ) :: &
	2062	g,mentrd_rate,mentr_rate
	2063	integer, dimension (its:ite) &
	2064	,intent (in ) :: &
	2065	kbcon,ktop,k22,jmin,kdet,kpbl
	2066	integer, dimension (its:ite) &
	2067	,intent (inout) :: &
	2068	ierr
	2069	character (), intent (in) :: &
	2070	name
	2071	!
	2072	! local variables in this routine
	2073	!
	2074
	2075	integer i,k
	2076	real detdo1,detdo2,entdo,dp,dz,subin,detdo,entup, &
	2077	detup,subdown,entdoj,entupk,detupk,totmas
	2078	!
	2079	i=ipr
	2080	DO K=kts+1,ktf
	2081	do i=its,itf
	2082	della(i,k)=0.
	2083	subs(i,k)=0.
	2084	enddo
	2085	enddo
	2086	!
	2087	DO 100 k=kts+1,ktf-1
	2088	DO 100 i=its,ite
	2089	IF(ierr(i).ne.0)GO TO 100
	2090	IF(K.Gt.KTOP(I))GO TO 100
	2091	!
	2092	!--- SPECIFY DETRAINMENT OF DOWNDRAFT, HAS TO BE CONSISTENT
	2093	!--- WITH ZD CALCULATIONS IN SOUNDD.
	2094	!
	2095	DZ=Z_cup(I,K+1)-Z_cup(I,K)
	2096	detdo=edt(i)CDD(i,K)DZ*ZD(i,k+1)
	2097	entdo=edt(i)mentrd_ratedz*zd(i,k+1)
	2098	!3d subin=zu(i,k+1)-zd(i,k+1)*edt(i)
	2099	subin=-zd(i,k+1)*edt(i)
	2100	entup=0.
	2101	detup=0.
	2102	if(k.ge.kbcon(i).and.k.lt.ktop(i))then
	2103	entup=mentr_ratedzzu(i,k)
	2104	detup=CD(i,K+1)DZZU(i,k)
	2105	endif
	2106	!3d subdown=(zu(i,k)-zd(i,k)*edt(i))
	2107	subdown=-zd(i,k)*edt(i)
	2108	entdoj=0.
	2109	entupk=0.
	2110	detupk=0.
	2111	!
	2112	if(k.eq.jmin(i))then
	2113	entdoj=edt(i)*zd(i,k)
	2114	endif
	2115
	2116	if(k.eq.k22(i)-1)then
	2117	entupk=zu(i,kpbl(i))
	2118	subin=zu(i,k+1)-zd(i,k+1)*edt(i)
	2119	if(high_res.eq.1)subin=-zd(i,k+1)*edt(i)
	2120	! subin=-zd(i,k+1)*edt(i)
	2121	endif
	2122
	2123	if(k.gt.kdet(i))then
	2124	detdo=0.
	2125	endif
	2126
	2127	if(k.eq.ktop(i)-0)then
	2128	detupk=zu(i,ktop(i))
	2129	subin=0.
	2130	!
	2131	! this subsidene for ktop now in subs term!
	2132	! subdown=zu(i,k)
	2133	subdown=0.
	2134	endif
	2135	if(k.lt.kbcon(i))then
	2136	detup=0.
	2137	endif
	2138	!C
	2139	!C--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
	2140	!C
	2141	totmas=subin-subdown+detup-entup-entdo+ &
	2142	detdo-entupk-entdoj+detupk
	2143	! if(j.eq.jpr.and.i.eq.ipr)print *,'k,totmas,sui,sud = ',k,
	2144	! 1 totmas,subin,subdown
	2145	! if(j.eq.jpr.and.i.eq.ipr)print *,'updr stuff = ',detup,
	2146	! 1 entup,entupk,detupk
	2147	! if(j.eq.jpr.and.i.eq.ipr)print *,'dddr stuff = ',entdo,
	2148	! 1 detdo,entdoj
	2149	if(abs(totmas).gt.1.e-6)then
	2150	! print ,'********************',i,j,k,totmas,name
	2151	! print *,kpbl(i),k22(i),kbcon(i),ktop(i)
	2152	!c print *,'updr stuff = ',subin,
	2153	!c 1 subdown,detup,entup,entupk,detupk
	2154	!c print *,'dddr stuff = ',entdo,
	2155	!c 1 detdo,entdoj
	2156	! call wrf_error_fatal ( 'totmas .gt.1.e-6' )
	2157	endif
	2158	dp=100.*(p_cup(i,k-1)-p_cup(i,k))
	2159	della(i,k)=(detup.5(HC(i,K+1)+HC(i,K)) &
	2160	+detdo.5(HCD(i,K+1)+HCD(i,K)) &
	2161	-entup*he(i,k) &
	2162	-entdo*he(i,k) &
	2163	+subin*he_cup(i,k+1) &
	2164	-subdown*he_cup(i,k) &
	2165	+detupk*(hc(i,ktop(i))-he_cup(i,ktop(i))) &
	2166	-entupk*he_cup(i,k22(i)) &
	2167	-entdoj*he_cup(i,jmin(i)) &
	2168	)*g/dp
	2169	if(high_res.eq.1)then
	2170	! the first term includes entr and detr into/from updraft as well as (entup-detup)*he(i,k) from
	2171	! neighbouring point, to make things mass consistent....
	2172	! if(k.ge.k22(i))then
	2173	della(i,k)=( &
	2174	detup.5(HC(i,K+1)+HC(i,K))-entuphe(i,k)+(entup-detup)he(i,k) &
	2175	+detdo.5(HCD(i,K+1)+HCD(i,K)) &
	2176	-entdo*he(i,k) &
	2177	+subin*he_cup(i,k+1) &
	2178	-subdown*he_cup(i,k) &
	2179	+detupk*(hc(i,ktop(i))-he(i,ktop(i))) &
	2180	-entdoj*he_cup(i,jmin(i)) &
	2181	-entupkhe_cup(i,k22(i))+entupkhe(i,k) &
	2182	)*g/dp
	2183	! else if(k.eq.k22(i)-1)then
	2184	! della(i,k)=(-entupkhe_cup(i,k22(i))+entupkhe(i,k))*g/dp
	2185	endif
	2186	!3d subin=zu(i,k+1)-zd(i,k+1)*edt(i)
	2187	!
	2188	! updraft subsidence only
	2189	!
	2190	if(k.ge.k22(i).and.k.lt.ktop(i))then
	2191	subs(i,k)=(zu(i,k+1)*he_cup(i,k+1) &
	2192	-zu(i,k)he_cup(i,k))g/dp
	2193	! else if(k.eq.ktop(i))then
	2194	! subs(i,k)=-detupkhe_cup(i,k)g/dp
	2195	endif
	2196	!
	2197	! in igh res case, subsidence terms are for meighbouring points only. This has to be
	2198	! done mass consistent with the della term
	2199	if(high_res.eq.1)then
	2200	if(k.ge.k22(i).and.k.lt.ktop(i))then
	2201	subs(i,k)=(zu(i,k+1)he_cup(i,k+1)-zu(i,k)he_cup(i,k)-(entup-detup)he(i,k))g/dp
	2202	else if(k.eq.ktop(i))then
	2203	subs(i,k)=detupk(he(i,ktop(i))-he_cup(i,ktop(i)))g/dp
	2204	else if(k.eq.k22(i)-1)then
	2205	subs(i,k)=(entupkhe(i,k)-entupkhe_cup(i,k))*g/dp
	2206	endif
	2207	endif
	2208	if(i.eq.ipr.and.j.eq.jpr)then
	2209	write(0,*)'d',k,della(i,k),subs(i,k),subin,subdown
	2210	! write(0,*)'d',detup,entup,entdo,entupk,entdoj
	2211	! print ,k,della(i,k),subinhe_cup(i,k+1),subdown*he_cup(i,k),
	2212	! 1 detdo.5(HCD(i,K+1)+HCD(i,K))
	2213	! print ,k,detup.5(HC(i,K+1)+HC(i,K)),detupkhc(i,ktop(i)),
	2214	! 1 entuphe(i,k),entdohe(i,k)
	2215	! print ,k,he_cup(i,k+1),he_cup(i,k),entupkhe_cup(i,k)
	2216	endif
	2217
	2218	100 CONTINUE
	2219
	2220	END SUBROUTINE cup_dellas_3d
	2221
	2222
	2223	SUBROUTINE cup_direction2(i,j,dir,id,massflx, &
	2224	iresult,imass,massfld, &
	2225	itf,jtf,ktf, &
	2226	its,ite, jts,jte, kts,kte )
	2227
	2228	IMPLICIT NONE
	2229
	2230	integer &
	2231	,intent (in ) :: &
	2232	itf,jtf,ktf, &
	2233	its,ite, jts,jte, kts,kte
	2234	integer, intent (in ) :: &
	2235	i,j,imass
	2236	integer, intent (out ) :: &
	2237	iresult
	2238	!
	2239	! ierr error value, maybe modified in this routine
	2240	!
	2241	integer, dimension (its:ite,jts:jte) &
	2242	,intent (in ) :: &
	2243	id
	2244	real, dimension (its:ite,jts:jte) &
	2245	,intent (in ) :: &
	2246	massflx
	2247	real, dimension (its:ite) &
	2248	,intent (inout) :: &
	2249	dir
	2250	real &
	2251	,intent (out ) :: &
	2252	massfld
	2253	!
	2254	! local variables in this routine
	2255	!
	2256
	2257	integer k,ia,ja,ib,jb
	2258	real diff
	2259	!
	2260	!
	2261	!
	2262	if(imass.eq.1)then
	2263	massfld=massflx(i,j)
	2264	endif
	2265	iresult=0
	2266	! return
	2267	diff=22.5
	2268	if(dir(i).lt.22.5)dir(i)=360.+dir(i)
	2269	if(id(i,j).eq.1)iresult=1
	2270	! ja=max(2,j-1)
	2271	! ia=max(2,i-1)
	2272	! jb=min(mjx-1,j+1)
	2273	! ib=min(mix-1,i+1)
	2274	ja=j-1
	2275	ia=i-1
	2276	jb=j+1
	2277	ib=i+1
	2278	if(dir(i).gt.90.-diff.and.dir(i).le.90.+diff)then
	2279	!--- steering flow from the east
	2280	if(id(ib,j).eq.1)then
	2281	iresult=1
	2282	if(imass.eq.1)then
	2283	massfld=max(massflx(ib,j),massflx(i,j))
	2284	endif
	2285	return
	2286	endif
	2287	else if(dir(i).gt.135.-diff.and.dir(i).le.135.+diff)then
	2288	!--- steering flow from the south-east
	2289	if(id(ib,ja).eq.1)then
	2290	iresult=1
	2291	if(imass.eq.1)then
	2292	massfld=max(massflx(ib,ja),massflx(i,j))
	2293	endif
	2294	return
	2295	endif
	2296	!--- steering flow from the south
	2297	else if(dir(i).gt.180.-diff.and.dir(i).le.180.+diff)then
	2298	if(id(i,ja).eq.1)then
	2299	iresult=1
	2300	if(imass.eq.1)then
	2301	massfld=max(massflx(i,ja),massflx(i,j))
	2302	endif
	2303	return
	2304	endif
	2305	!--- steering flow from the south west
	2306	else if(dir(i).gt.225.-diff.and.dir(i).le.225.+diff)then
	2307	if(id(ia,ja).eq.1)then
	2308	iresult=1
	2309	if(imass.eq.1)then
	2310	massfld=max(massflx(ia,ja),massflx(i,j))
	2311	endif
	2312	return
	2313	endif
	2314	!--- steering flow from the west
	2315	else if(dir(i).gt.270.-diff.and.dir(i).le.270.+diff)then
	2316	if(id(ia,j).eq.1)then
	2317	iresult=1
	2318	if(imass.eq.1)then
	2319	massfld=max(massflx(ia,j),massflx(i,j))
	2320	endif
	2321	return
	2322	endif
	2323	!--- steering flow from the north-west
	2324	else if(dir(i).gt.305.-diff.and.dir(i).le.305.+diff)then
	2325	if(id(ia,jb).eq.1)then
	2326	iresult=1
	2327	if(imass.eq.1)then
	2328	massfld=max(massflx(ia,jb),massflx(i,j))
	2329	endif
	2330	return
	2331	endif
	2332	!--- steering flow from the north
	2333	else if(dir(i).gt.360.-diff.and.dir(i).le.360.+diff)then
	2334	if(id(i,jb).eq.1)then
	2335	iresult=1
	2336	if(imass.eq.1)then
	2337	massfld=max(massflx(i,jb),massflx(i,j))
	2338	endif
	2339	return
	2340	endif
	2341	!--- steering flow from the north-east
	2342	else if(dir(i).gt.45.-diff.and.dir(i).le.45.+diff)then
	2343	if(id(ib,jb).eq.1)then
	2344	iresult=1
	2345	if(imass.eq.1)then
	2346	massfld=max(massflx(ib,jb),massflx(i,j))
	2347	endif
	2348	return
	2349	endif
	2350	endif
	2351
	2352	END SUBROUTINE cup_direction2
	2353
	2354
	2355	SUBROUTINE cup_env(z,qes,he,hes,t,q,p,z1, &
	2356	psur,ierr,tcrit,itest,xl,cp, &
	2357	itf,jtf,ktf, &
	2358	its,ite, jts,jte, kts,kte )
	2359
	2360	IMPLICIT NONE
	2361
	2362	integer &
	2363	,intent (in ) :: &
	2364	itf,jtf,ktf, &
	2365	its,ite, jts,jte, kts,kte
	2366	!
	2367	! ierr error value, maybe modified in this routine
	2368	! q = environmental mixing ratio
	2369	! qes = environmental saturation mixing ratio
	2370	! t = environmental temp
	2371	! tv = environmental virtual temp
	2372	! p = environmental pressure
	2373	! z = environmental heights
	2374	! he = environmental moist static energy
	2375	! hes = environmental saturation moist static energy
	2376	! psur = surface pressure
	2377	! z1 = terrain elevation
	2378	!
	2379	!
	2380	real, dimension (its:ite,kts:kte) &
	2381	,intent (in ) :: &
	2382	p,t
	2383	real, dimension (its:ite,kts:kte) &
	2384	,intent (out ) :: &
	2385	he,hes,qes
	2386	real, dimension (its:ite,kts:kte) &
	2387	,intent (inout) :: &
	2388	z,q
	2389	real, dimension (its:ite) &
	2390	,intent (in ) :: &
	2391	psur,z1
	2392	real &
	2393	,intent (in ) :: &
	2394	xl,cp
	2395	integer, dimension (its:ite) &
	2396	,intent (inout) :: &
	2397	ierr
	2398	integer &
	2399	,intent (in ) :: &
	2400	itest
	2401	!
	2402	! local variables in this routine
	2403	!
	2404
	2405	integer :: &
	2406	i,k,iph
	2407	real, dimension (1:2) :: AE,BE,HT
	2408	real, dimension (its:ite,kts:kte) :: tv
	2409	real :: tcrit,e,tvbar
	2410
	2411
	2412	HT(1)=XL/CP
	2413	HT(2)=2.834E6/CP
	2414	BE(1)=.622*HT(1)/.286
	2415	AE(1)=BE(1)/273.+ALOG(610.71)
	2416	BE(2)=.622*HT(2)/.286
	2417	AE(2)=BE(2)/273.+ALOG(610.71)
	2418	! print *, 'TCRIT = ', tcrit,its,ite
	2419	DO k=kts,ktf
	2420	do i=its,itf
	2421	if(ierr(i).eq.0)then
	2422	!Csgb - IPH is for phase, dependent on TCRIT (water or ice)
	2423	IPH=1
	2424	IF(T(I,K).LE.TCRIT)IPH=2
	2425	! print *, 'AE(IPH),BE(IPH) = ',AE(IPH),BE(IPH),AE(IPH)-BE(IPH),T(i,k),i,k
	2426	E=EXP(AE(IPH)-BE(IPH)/T(I,K))
	2427	! print *, 'P, E = ', P(I,K), E
	2428	QES(I,K)=.622E/(100.P(I,K)-E)
	2429	IF(QES(I,K).LE.1.E-08)QES(I,K)=1.E-08
	2430	IF(Q(I,K).GT.QES(I,K))Q(I,K)=QES(I,K)
	2431	TV(I,K)=T(I,K)+.608Q(I,K)T(I,K)
	2432	endif
	2433	enddo
	2434	enddo
	2435	!
	2436	!--- z's are calculated with changed h's and q's and t's
	2437	!--- if itest=2
	2438	!
	2439	if(itest.ne.2)then
	2440	do i=its,itf
	2441	if(ierr(i).eq.0)then
	2442	Z(I,1)=max(0.,Z1(I))-(ALOG(P(I,1))- &
	2443	ALOG(PSUR(I)))287.TV(I,1)/9.81
	2444	endif
	2445	enddo
	2446
	2447	! --- calculate heights
	2448	DO K=kts+1,ktf
	2449	do i=its,itf
	2450	if(ierr(i).eq.0)then
	2451	TVBAR=.5TV(I,K)+.5TV(I,K-1)
	2452	Z(I,K)=Z(I,K-1)-(ALOG(P(I,K))- &
	2453	ALOG(P(I,K-1)))287.TVBAR/9.81
	2454	endif
	2455	enddo
	2456	enddo
	2457	else
	2458	do k=kts,ktf
	2459	do i=its,itf
	2460	if(ierr(i).eq.0)then
	2461	z(i,k)=(he(i,k)-1004.t(i,k)-2.5e6q(i,k))/9.81
	2462	z(i,k)=max(1.e-3,z(i,k))
	2463	endif
	2464	enddo
	2465	enddo
	2466	endif
	2467	!
	2468	!--- calculate moist static energy - HE
	2469	! saturated moist static energy - HES
	2470	!
	2471	DO k=kts,ktf
	2472	do i=its,itf
	2473	if(ierr(i).eq.0)then
	2474	if(itest.eq.0)HE(I,K)=9.81Z(I,K)+1004.T(I,K)+2.5E06*Q(I,K)
	2475	HES(I,K)=9.81Z(I,K)+1004.T(I,K)+2.5E06*QES(I,K)
	2476	IF(HE(I,K).GE.HES(I,K))HE(I,K)=HES(I,K)
	2477	endif
	2478	enddo
	2479	enddo
	2480
	2481	END SUBROUTINE cup_env
	2482
	2483
	2484	SUBROUTINE cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup, &
	2485	he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
	2486	ierr,z1,xl,rv,cp, &
	2487	itf,jtf,ktf, &
	2488	its,ite, jts,jte, kts,kte )
	2489
	2490	IMPLICIT NONE
	2491
	2492	integer &
	2493	,intent (in ) :: &
	2494	itf,jtf,ktf, &
	2495	its,ite, jts,jte, kts,kte
	2496	!
	2497	! ierr error value, maybe modified in this routine
	2498	! q = environmental mixing ratio
	2499	! q_cup = environmental mixing ratio on cloud levels
	2500	! qes = environmental saturation mixing ratio
	2501	! qes_cup = environmental saturation mixing ratio on cloud levels
	2502	! t = environmental temp
	2503	! t_cup = environmental temp on cloud levels
	2504	! p = environmental pressure
	2505	! p_cup = environmental pressure on cloud levels
	2506	! z = environmental heights
	2507	! z_cup = environmental heights on cloud levels
	2508	! he = environmental moist static energy
	2509	! he_cup = environmental moist static energy on cloud levels
	2510	! hes = environmental saturation moist static energy
	2511	! hes_cup = environmental saturation moist static energy on cloud levels
	2512	! gamma_cup = gamma on cloud levels
	2513	! psur = surface pressure
	2514	! z1 = terrain elevation
	2515	!
	2516	!
	2517	real, dimension (its:ite,kts:kte) &
	2518	,intent (in ) :: &
	2519	qes,q,he,hes,z,p,t
	2520	real, dimension (its:ite,kts:kte) &
	2521	,intent (out ) :: &
	2522	qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup
	2523	real, dimension (its:ite) &
	2524	,intent (in ) :: &
	2525	psur,z1
	2526	real &
	2527	,intent (in ) :: &
	2528	xl,rv,cp
	2529	integer, dimension (its:ite) &
	2530	,intent (inout) :: &
	2531	ierr
	2532	!
	2533	! local variables in this routine
	2534	!
	2535
	2536	integer :: &
	2537	i,k
	2538
	2539
	2540	do k=kts+1,ktf
	2541	do i=its,itf
	2542	if(ierr(i).eq.0)then
	2543	qes_cup(i,k)=.5*(qes(i,k-1)+qes(i,k))
	2544	q_cup(i,k)=.5*(q(i,k-1)+q(i,k))
	2545	hes_cup(i,k)=.5*(hes(i,k-1)+hes(i,k))
	2546	he_cup(i,k)=.5*(he(i,k-1)+he(i,k))
	2547	if(he_cup(i,k).gt.hes_cup(i,k))he_cup(i,k)=hes_cup(i,k)
	2548	z_cup(i,k)=.5*(z(i,k-1)+z(i,k))
	2549	p_cup(i,k)=.5*(p(i,k-1)+p(i,k))
	2550	t_cup(i,k)=.5*(t(i,k-1)+t(i,k))
	2551	gamma_cup(i,k)=(xl/cp)(xl/(rvt_cup(i,k) &
	2552	t_cup(i,k)))qes_cup(i,k)
	2553	endif
	2554	enddo
	2555	enddo
	2556	do i=its,itf
	2557	if(ierr(i).eq.0)then
	2558	qes_cup(i,1)=qes(i,1)
	2559	q_cup(i,1)=q(i,1)
	2560	hes_cup(i,1)=hes(i,1)
	2561	he_cup(i,1)=he(i,1)
	2562	z_cup(i,1)=.5*(z(i,1)+z1(i))
	2563	p_cup(i,1)=.5*(p(i,1)+psur(i))
	2564	t_cup(i,1)=t(i,1)
	2565	gamma_cup(i,1)=xl/cp(xl/(rvt_cup(i,1) &
	2566	t_cup(i,1)))qes_cup(i,1)
	2567	endif
	2568	enddo
	2569
	2570	END SUBROUTINE cup_env_clev
	2571
	2572
	2573	SUBROUTINE cup_forcing_ens_3d(closure_n,xland,aa0,aa1,xaa0,mbdt,dtime,ierr,ierr2,ierr3,&
	2574	xf_ens,j,name,axx,maxens,iens,iedt,maxens2,maxens3,mconv, &
	2575	p_cup,ktop,omeg,zd,k22,zu,pr_ens,edt,kbcon,massflx, &
	2576	iact_old_gr,dir,ensdim,massfln,icoic,edt_out, &
	2577	high_resolution,itf,jtf,ktf, &
	2578	its,ite, jts,jte, kts,kte,ens4,ktau )
	2579
	2580	IMPLICIT NONE
	2581
	2582	integer &
	2583	,intent (in ) :: &
	2584	itf,jtf,ktf, &
	2585	its,ite, jts,jte, kts,kte,ens4,high_resolution,ktau
	2586	integer, intent (in ) :: &
	2587	j,ensdim,maxens,iens,iedt,maxens2,maxens3
	2588	!
	2589	! ierr error value, maybe modified in this routine
	2590	! pr_ens = precipitation ensemble
	2591	! xf_ens = mass flux ensembles
	2592	! massfln = downdraft mass flux ensembles used in next timestep
	2593	! omeg = omega from large scale model
	2594	! mconv = moisture convergence from large scale model
	2595	! zd = downdraft normalized mass flux
	2596	! zu = updraft normalized mass flux
	2597	! aa0 = cloud work function without forcing effects
	2598	! aa1 = cloud work function with forcing effects
	2599	! xaa0 = cloud work function with cloud effects (ensemble dependent)
	2600	! edt = epsilon
	2601	! dir = "storm motion"
	2602	! mbdt = arbitrary numerical parameter
	2603	! dtime = dt over which forcing is applied
	2604	! iact_gr_old = flag to tell where convection was active
	2605	! kbcon = LFC of parcel from k22
	2606	! k22 = updraft originating level
	2607	! icoic = flag if only want one closure (usually set to zero!)
	2608	! name = deep or shallow convection flag
	2609	!
	2610	real, dimension (its:ite,jts:jte,1:ensdim) &
	2611	,intent (inout) :: &
	2612	pr_ens
	2613	real, dimension (its:ite,jts:jte,1:ensdim) &
	2614	,intent (out ) :: &
	2615	xf_ens,massfln
	2616	real, dimension (its:ite,jts:jte) &
	2617	,intent (inout ) :: &
	2618	edt_out
	2619	real, dimension (its:ite,jts:jte) &
	2620	,intent (in ) :: &
	2621	massflx
	2622	real, dimension (its:ite,kts:kte) &
	2623	,intent (in ) :: &
	2624	zd,zu,p_cup
	2625	real, dimension (its:ite,kts:kte,1:ens4) &
	2626	,intent (in ) :: &
	2627	omeg
	2628	real, dimension (its:ite,1:maxens) &
	2629	,intent (in ) :: &
	2630	xaa0
	2631	real, dimension (its:ite) &
	2632	,intent (in ) :: &
	2633	aa1,edt,dir,xland
	2634	real, dimension (its:ite,1:ens4) &
	2635	,intent (in ) :: &
	2636	mconv,axx
	2637	real, dimension (its:ite) &
	2638	,intent (inout) :: &
	2639	aa0,closure_n
	2640	real, dimension (1:maxens) &
	2641	,intent (in ) :: &
	2642	mbdt
	2643	real &
	2644	,intent (in ) :: &
	2645	dtime
	2646	integer, dimension (its:ite,jts:jte) &
	2647	,intent (in ) :: &
	2648	iact_old_gr
	2649	integer, dimension (its:ite) &
	2650	,intent (in ) :: &
	2651	k22,kbcon,ktop
	2652	integer, dimension (its:ite) &
	2653	,intent (inout) :: &
	2654	ierr,ierr2,ierr3
	2655	integer &
	2656	,intent (in ) :: &
	2657	icoic
	2658	character (), intent (in) :: &
	2659	name
	2660	!
	2661	! local variables in this routine
	2662	!
	2663
	2664	real, dimension (1:maxens3) :: &
	2665	xff_ens3
	2666	real, dimension (1:maxens) :: &
	2667	xk
	2668	integer :: &
	2669	i,k,nall,n,ne,nens,nens3,iresult,iresultd,iresulte,mkxcrt,kclim
	2670	parameter (mkxcrt=15)
	2671	real :: &
	2672	fens4,a1,massfld,a_ave,xff0,xff00,xxx,xomg,aclim1,aclim2,aclim3,aclim4
	2673	real, dimension(1:mkxcrt) :: &
	2674	pcrit,acrit,acritt
	2675
	2676	integer :: nall2,ixxx,irandom
	2677	integer, dimension (8) :: seed
	2678
	2679
	2680	DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
	2681	350.,300.,250.,200.,150./
	2682	DATA ACRIT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
	2683	.3151,.3677,.41,.5255,.7663,1.1686,1.6851/
	2684	! GDAS DERIVED ACRIT
	2685	DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688, &
	2686	.743,.813,.886,.947,1.138,1.377,1.896/
	2687	!
	2688	seed=0
	2689	seed(2)=j
	2690	seed(3)=ktau
	2691	nens=0
	2692	irandom=1
	2693	if(high_resolution.eq.1)irandom=0
	2694	irandom=0
	2695	fens4=float(ens4)
	2696
	2697	!--- LARGE SCALE FORCING
	2698	!
	2699	DO 100 i=its,itf
	2700	if(name.eq.'deeps'.and.ierr(i).gt.995)then
	2701	aa0(i)=0.
	2702	ierr(i)=0
	2703	endif
	2704	IF(ierr(i).eq.0)then
	2705	!
	2706	!---
	2707	!
	2708	if(name.eq.'deeps')then
	2709	!
	2710	a_ave=0.
	2711	do ne=1,ens4
	2712	a_ave=a_ave+axx(i,ne)
	2713	enddo
	2714	a_ave=max(0.,a_ave/fens4)
	2715	a_ave=min(a_ave,aa1(i))
	2716	a_ave=max(0.,a_ave)
	2717	do ne=1,16
	2718	xff_ens3(ne)=0.
	2719	enddo
	2720	xff0= (AA1(I)-AA0(I))/DTIME
	2721	if(high_resolution.eq.1)xff0= (a_ave-AA0(I))/DTIME
	2722	xff_ens3(1)=(AA1(I)-AA0(I))/dtime
	2723	xff_ens3(2)=(a_ave-AA0(I))/dtime
	2724	if(irandom.eq.1)then
	2725	seed(1)=i
	2726	call random_seed (PUT=seed)
	2727	call random_number (xxx)
	2728	ixxx=min(ens4,max(1,int(fens4*xxx+1.e-8)))
	2729	xff_ens3(3)=(axx(i,ixxx)-AA0(I))/dtime
	2730	call random_number (xxx)
	2731	ixxx=min(ens4,max(1,int(fens4*xxx+1.e-8)))
	2732	xff_ens3(13)=(axx(i,ixxx)-AA0(I))/dtime
	2733	else
	2734	xff_ens3(3)=(AA1(I)-AA0(I))/dtime
	2735	xff_ens3(13)=(AA1(I)-AA0(I))/dtime
	2736	endif
	2737	if(high_resolution.eq.1)then
	2738	xff_ens3(1)=(a_ave-AA0(I))/dtime
	2739	xff_ens3(2)=(a_ave-AA0(I))/dtime
	2740	xff_ens3(3)=(a_ave-AA0(I))/dtime
	2741	xff_ens3(13)=(a_ave-AA0(I))/dtime
	2742	endif
	2743	!
	2744	!--- more original Arakawa-Schubert (climatologic value of aa0)
	2745	!
	2746	!
	2747	!--- omeg is in bar/s, mconv done with omeg in Pa/s
	2748	! more like Brown (1979), or Frank-Cohen (199?)
	2749	!
	2750	xff_ens3(14)=0.
	2751	do ne=1,ens4
	2752	xff_ens3(14)=xff_ens3(14)-omeg(i,k22(i),ne)/(fens4*9.81)
	2753	enddo
	2754	if(xff_ens3(14).lt.0.)xff_ens3(14)=0.
	2755	xff_ens3(5)=0.
	2756	do ne=1,ens4
	2757	xff_ens3(5)=xff_ens3(5)-omeg(i,kbcon(i),ne)/(fens4*9.81)
	2758	enddo
	2759	if(xff_ens3(5).lt.0.)xff_ens3(5)=0.
	2760	!
	2761	! minimum below kbcon
	2762	!
	2763	if(high_resolution.eq.0)then
	2764	xff_ens3(4)=-omeg(i,2,1)/9.81
	2765	do k=2,kbcon(i)-1
	2766	do ne=1,ens4
	2767	xomg=-omeg(i,k,ne)/9.81
	2768	if(xomg.lt.xff_ens3(4))xff_ens3(4)=xomg
	2769	enddo
	2770	enddo
	2771	if(xff_ens3(4).lt.0.)xff_ens3(4)=0.
	2772	!
	2773	! max below kbcon
	2774	xff_ens3(6)=-omeg(i,2,1)/9.81
	2775	do k=2,kbcon(i)-1
	2776	do ne=1,ens4
	2777	xomg=-omeg(i,k,ne)/9.81
	2778	if(xomg.gt.xff_ens3(6))xff_ens3(6)=xomg
	2779	enddo
	2780	enddo
	2781	if(xff_ens3(6).lt.0.)xff_ens3(6)=0.
	2782	endif
	2783	if(high_resolution.eq.1)then
	2784	xff_ens3(5)=min(xff_ens3(5),xff_ens3(14))
	2785	xff_ens3(4)=xff_ens3(5)
	2786	xff_ens3(6)=xff_ens3(5)
	2787	endif
	2788	!
	2789	!--- more like Krishnamurti et al.; pick max and average values
	2790	!
	2791	xff_ens3(7)=mconv(i,1)
	2792	xff_ens3(8)=mconv(i,1)
	2793	xff_ens3(9)=mconv(i,1)
	2794	if(ens4.gt.1)then
	2795	do ne=2,ens4
	2796	if (mconv(i,ne).gt.xff_ens3(7))xff_ens3(7)=mconv(i,ne)
	2797	enddo
	2798	do ne=2,ens4
	2799	if (mconv(i,ne).lt.xff_ens3(8))xff_ens3(8)=mconv(i,ne)
	2800	enddo
	2801	do ne=2,ens4
	2802	xff_ens3(9)=xff_ens3(9)+mconv(i,ne)
	2803	enddo
	2804	xff_ens3(9)=xff_ens3(9)/fens4
	2805	endif
	2806	if(high_resolution.eq.1)then
	2807	xff_ens3(7)=xff_ens3(9)
	2808	xff_ens3(8)=xff_ens3(9)
	2809	xff_ens3(15)=xff_ens3(9)
	2810	endif
	2811	!
	2812	if(high_resolution.eq.0)then
	2813	if(irandom.eq.1)then
	2814	seed(1)=i
	2815	call random_seed (PUT=seed)
	2816	call random_number (xxx)
	2817	ixxx=min(ens4,max(1,int(fens4*xxx+1.e-8)))
	2818	xff_ens3(15)=mconv(i,ixxx)
	2819	else
	2820	xff_ens3(15)=mconv(i,1)
	2821	endif
	2822	endif
	2823	!
	2824	!--- more like Fritsch Chappel or Kain Fritsch (plus triggers)
	2825	!
	2826	xff_ens3(10)=A_AVE/(60.*40.)
	2827	xff_ens3(11)=AA1(I)/(60.*40.)
	2828	if(irandom.eq.1)then
	2829	seed(1)=i
	2830	call random_seed (PUT=seed)
	2831	call random_number (xxx)
	2832	ixxx=min(ens4,max(1,int(fens4*xxx+1.e-8)))
	2833	xff_ens3(12)=AXX(I,ixxx)/(60.*40.)
	2834	else
	2835	xff_ens3(12)=AA1(I)/(60.*40.)
	2836	endif
	2837	if(high_resolution.eq.1)then
	2838	xff_ens3(11)=xff_ens3(10)
	2839	xff_ens3(12)=xff_ens3(10)
	2840	endif
	2841	!
	2842	!--- more original Arakawa-Schubert (climatologic value of aa0)
	2843	!
	2844	! edt_out(i,j)=xff0
	2845	if(icoic.eq.0)then
	2846	if(xff0.lt.0.)then
	2847	xff_ens3(1)=0.
	2848	xff_ens3(2)=0.
	2849	xff_ens3(3)=0.
	2850	xff_ens3(13)=0.
	2851	xff_ens3(10)=0.
	2852	xff_ens3(11)=0.
	2853	xff_ens3(12)=0.
	2854	endif
	2855	endif
	2856
	2857
	2858
	2859	do nens=1,maxens
	2860	XK(nens)=(XAA0(I,nens)-AA1(I))/MBDT(2)
	2861	if(xk(nens).le.0.and.xk(nens).gt.-1.e-6) &
	2862	xk(nens)=-1.e-6
	2863	if(xk(nens).gt.0.and.xk(nens).lt.1.e-6) &
	2864	xk(nens)=1.e-6
	2865	enddo
	2866	!
	2867	!--- add up all ensembles
	2868	!
	2869	do 350 ne=1,maxens
	2870	!
	2871	!--- for every xk, we have maxens3 xffs
	2872	!--- iens is from outermost ensemble (most expensive!
	2873	!
	2874	!--- iedt (maxens2 belongs to it)
	2875	!--- is from second, next outermost, not so expensive
	2876	!
	2877	!--- so, for every outermost loop, we have maxensmaxens23
	2878	!--- ensembles!!! nall would be 0, if everything is on first
	2879	!--- loop index, then ne would start counting, then iedt, then iens....
	2880	!
	2881	iresult=0
	2882	iresultd=0
	2883	iresulte=0
	2884	nall=(iens-1)maxens3maxens*maxens2 &
	2885	+(iedt-1)maxensmaxens3 &
	2886	+(ne-1)*maxens3
	2887	!
	2888	! over water, enfor!e small cap for some of the closures
	2889	!
	2890	if(xland(i).lt.0.1)then
	2891	if(ierr2(i).gt.0.or.ierr3(i).gt.0)then
	2892	xff_ens3(1) =0.
	2893	massfln(i,j,nall+1)=0.
	2894	xff_ens3(2) =0.
	2895	massfln(i,j,nall+2)=0.
	2896	xff_ens3(3) =0.
	2897	massfln(i,j,nall+3)=0.
	2898	xff_ens3(10) =0.
	2899	massfln(i,j,nall+10)=0.
	2900	xff_ens3(11) =0.
	2901	massfln(i,j,nall+11)=0.
	2902	xff_ens3(12) =0.
	2903	massfln(i,j,nall+12)=0.
	2904	xff_ens3(7) =0.
	2905	massfln(i,j,nall+7)=0.
	2906	xff_ens3(8) =0.
	2907	massfln(i,j,nall+8)=0.
	2908	xff_ens3(9) =0.
	2909	massfln(i,j,nall+9)=0.
	2910	closure_n(i)=closure_n(i)-1.
	2911	xff_ens3(13) =0.
	2912	massfln(i,j,nall+13)=0.
	2913	xff_ens3(15) =0.
	2914	massfln(i,j,nall+15)=0.
	2915	endif
	2916	endif
	2917	!
	2918	! end water treatment
	2919	!
	2920	!
	2921	!--- check for upwind convection
	2922	! iresult=0
	2923	massfld=0.
	2924
	2925	! call cup_direction2(i,j,dir,iact_old_gr, &
	2926	! massflx,iresult,1, &
	2927	! massfld, &
	2928	! itf,jtf,ktf, &
	2929	! ims,ime, jms,jme, kms,kme, &
	2930	! its,ite, jts,jte, kts,kte )
	2931	! if(i.eq.ipr.and.j.eq.jpr.and.iedt.eq.1.and.ne.eq.1)then
	2932	! if(iedt.eq.1.and.ne.eq.1)then
	2933	! print *,massfld,ne,iedt,iens
	2934	! print *,xk(ne),xff_ens3(1),xff_ens3(2),xff_ens3(3)
	2935	! endif
	2936	! print *,i,j,massfld,aa0(i),aa1(i)
	2937	IF(XK(ne).lt.0.and.xff0.gt.0.)iresultd=1
	2938	iresulte=max(iresult,iresultd)
	2939	iresulte=1
	2940	if(iresulte.eq.1)then
	2941	!
	2942	!--- special treatment for stability closures
	2943	!
	2944
	2945	if(xff0.ge.0.)then
	2946	xf_ens(i,j,nall+1)=massfld
	2947	xf_ens(i,j,nall+2)=massfld
	2948	xf_ens(i,j,nall+3)=massfld
	2949	xf_ens(i,j,nall+13)=massfld
	2950	if(xff_ens3(1).gt.0)xf_ens(i,j,nall+1)=max(0.,-xff_ens3(1)/xk(ne)) &
	2951	+massfld
	2952	if(xff_ens3(2).gt.0)xf_ens(i,j,nall+2)=max(0.,-xff_ens3(2)/xk(ne)) &
	2953	+massfld
	2954	if(xff_ens3(3).gt.0)xf_ens(i,j,nall+3)=max(0.,-xff_ens3(3)/xk(ne)) &
	2955	+massfld
	2956	if(xff_ens3(13).gt.0)xf_ens(i,j,nall+13)=max(0.,-xff_ens3(13)/xk(ne)) &
	2957	+massfld
	2958	! endif
	2959	else
	2960	xf_ens(i,j,nall+1)=massfld
	2961	xf_ens(i,j,nall+2)=massfld
	2962	xf_ens(i,j,nall+3)=massfld
	2963	xf_ens(i,j,nall+13)=massfld
	2964	endif
	2965	!
	2966	!--- if iresult.eq.1, following independent of xff0
	2967	!
	2968	xf_ens(i,j,nall+4)=max(0.,xff_ens3(4) &
	2969	+massfld)
	2970	xf_ens(i,j,nall+5)=max(0.,xff_ens3(5) &
	2971	+massfld)
	2972	xf_ens(i,j,nall+6)=max(0.,xff_ens3(6) &
	2973	+massfld)
	2974	xf_ens(i,j,nall+14)=max(0.,xff_ens3(14) &
	2975	+massfld)
	2976	a1=max(1.e-3,pr_ens(i,j,nall+7))
	2977	xf_ens(i,j,nall+7)=max(0.,xff_ens3(7) &
	2978	/a1)
	2979	a1=max(1.e-3,pr_ens(i,j,nall+8))
	2980	xf_ens(i,j,nall+8)=max(0.,xff_ens3(8) &
	2981	/a1)
	2982	a1=max(1.e-3,pr_ens(i,j,nall+9))
	2983	xf_ens(i,j,nall+9)=max(0.,xff_ens3(9) &
	2984	/a1)
	2985	a1=max(1.e-3,pr_ens(i,j,nall+15))
	2986	xf_ens(i,j,nall+15)=max(0.,xff_ens3(15) &
	2987	/a1)
	2988	if(XK(ne).lt.0.)then
	2989	xf_ens(i,j,nall+10)=max(0., &
	2990	-xff_ens3(10)/xk(ne)) &
	2991	+massfld
	2992	xf_ens(i,j,nall+11)=max(0., &
	2993	-xff_ens3(11)/xk(ne)) &
	2994	+massfld
	2995	xf_ens(i,j,nall+12)=max(0., &
	2996	-xff_ens3(12)/xk(ne)) &
	2997	+massfld
	2998	else
	2999	xf_ens(i,j,nall+10)=massfld
	3000	xf_ens(i,j,nall+11)=massfld
	3001	xf_ens(i,j,nall+12)=massfld
	3002	endif
	3003	if(icoic.ge.1)then
	3004	closure_n(i)=0.
	3005	xf_ens(i,j,nall+1)=xf_ens(i,j,nall+icoic)
	3006	xf_ens(i,j,nall+2)=xf_ens(i,j,nall+icoic)
	3007	xf_ens(i,j,nall+3)=xf_ens(i,j,nall+icoic)
	3008	xf_ens(i,j,nall+4)=xf_ens(i,j,nall+icoic)
	3009	xf_ens(i,j,nall+5)=xf_ens(i,j,nall+icoic)
	3010	xf_ens(i,j,nall+6)=xf_ens(i,j,nall+icoic)
	3011	xf_ens(i,j,nall+7)=xf_ens(i,j,nall+icoic)
	3012	xf_ens(i,j,nall+8)=xf_ens(i,j,nall+icoic)
	3013	xf_ens(i,j,nall+9)=xf_ens(i,j,nall+icoic)
	3014	xf_ens(i,j,nall+10)=xf_ens(i,j,nall+icoic)
	3015	xf_ens(i,j,nall+11)=xf_ens(i,j,nall+icoic)
	3016	xf_ens(i,j,nall+12)=xf_ens(i,j,nall+icoic)
	3017	xf_ens(i,j,nall+13)=xf_ens(i,j,nall+icoic)
	3018	xf_ens(i,j,nall+14)=xf_ens(i,j,nall+icoic)
	3019	xf_ens(i,j,nall+15)=xf_ens(i,j,nall+icoic)
	3020	xf_ens(i,j,nall+16)=xf_ens(i,j,nall+icoic)
	3021	endif
	3022	!
	3023	! 16 is a randon pick from the oher 15
	3024	!
	3025	if(irandom.eq.1)then
	3026	call random_number (xxx)
	3027	ixxx=min(15,max(1,int(15.*xxx+1.e-8)))
	3028	xf_ens(i,j,nall+16)=xf_ens(i,j,nall+ixxx)
	3029	else
	3030	xf_ens(i,j,nall+16)=xf_ens(i,j,nall+1)
	3031	endif
	3032	!
	3033	!
	3034	!--- store new for next time step
	3035	!
	3036	do nens3=1,maxens3
	3037	massfln(i,j,nall+nens3)=edt(i) &
	3038	*xf_ens(i,j,nall+nens3)
	3039	massfln(i,j,nall+nens3)=max(0., &
	3040	massfln(i,j,nall+nens3))
	3041	enddo
	3042	!
	3043	!
	3044	!--- do some more on the caps!!! ne=1 for 175, ne=2 for 100,....
	3045	!
	3046	! do not care for caps here for closure groups 1 and 5,
	3047	! they are fine, do not turn them off here
	3048	!
	3049	!
	3050	if(ne.eq.2.and.ierr2(i).gt.0)then
	3051	xf_ens(i,j,nall+1) =0.
	3052	xf_ens(i,j,nall+2) =0.
	3053	xf_ens(i,j,nall+3) =0.
	3054	xf_ens(i,j,nall+4) =0.
	3055	xf_ens(i,j,nall+5) =0.
	3056	xf_ens(i,j,nall+6) =0.
	3057	xf_ens(i,j,nall+7) =0.
	3058	xf_ens(i,j,nall+8) =0.
	3059	xf_ens(i,j,nall+9) =0.
	3060	xf_ens(i,j,nall+10)=0.
	3061	xf_ens(i,j,nall+11)=0.
	3062	xf_ens(i,j,nall+12)=0.
	3063	xf_ens(i,j,nall+13)=0.
	3064	xf_ens(i,j,nall+14)=0.
	3065	xf_ens(i,j,nall+15)=0.
	3066	xf_ens(i,j,nall+16)=0.
	3067	massfln(i,j,nall+1)=0.
	3068	massfln(i,j,nall+2)=0.
	3069	massfln(i,j,nall+3)=0.
	3070	massfln(i,j,nall+4)=0.
	3071	massfln(i,j,nall+5)=0.
	3072	massfln(i,j,nall+6)=0.
	3073	massfln(i,j,nall+7)=0.
	3074	massfln(i,j,nall+8)=0.
	3075	massfln(i,j,nall+9)=0.
	3076	massfln(i,j,nall+10)=0.
	3077	massfln(i,j,nall+11)=0.
	3078	massfln(i,j,nall+12)=0.
	3079	massfln(i,j,nall+13)=0.
	3080	massfln(i,j,nall+14)=0.
	3081	massfln(i,j,nall+15)=0.
	3082	massfln(i,j,nall+16)=0.
	3083	endif
	3084	if(ne.eq.3.and.ierr3(i).gt.0)then
	3085	xf_ens(i,j,nall+1) =0.
	3086	xf_ens(i,j,nall+2) =0.
	3087	xf_ens(i,j,nall+3) =0.
	3088	xf_ens(i,j,nall+4) =0.
	3089	xf_ens(i,j,nall+5) =0.
	3090	xf_ens(i,j,nall+6) =0.
	3091	xf_ens(i,j,nall+7) =0.
	3092	xf_ens(i,j,nall+8) =0.
	3093	xf_ens(i,j,nall+9) =0.
	3094	xf_ens(i,j,nall+10)=0.
	3095	xf_ens(i,j,nall+11)=0.
	3096	xf_ens(i,j,nall+12)=0.
	3097	xf_ens(i,j,nall+13)=0.
	3098	xf_ens(i,j,nall+14)=0.
	3099	xf_ens(i,j,nall+15)=0.
	3100	xf_ens(i,j,nall+16)=0.
	3101	massfln(i,j,nall+1)=0.
	3102	massfln(i,j,nall+2)=0.
	3103	massfln(i,j,nall+3)=0.
	3104	massfln(i,j,nall+4)=0.
	3105	massfln(i,j,nall+5)=0.
	3106	massfln(i,j,nall+6)=0.
	3107	massfln(i,j,nall+7)=0.
	3108	massfln(i,j,nall+8)=0.
	3109	massfln(i,j,nall+9)=0.
	3110	massfln(i,j,nall+10)=0.
	3111	massfln(i,j,nall+11)=0.
	3112	massfln(i,j,nall+12)=0.
	3113	massfln(i,j,nall+13)=0.
	3114	massfln(i,j,nall+14)=0.
	3115	massfln(i,j,nall+15)=0.
	3116	massfln(i,j,nall+16)=0.
	3117	endif
	3118
	3119	endif
	3120	350 continue
	3121	! ne=1, cap=175
	3122	!
	3123	nall=(iens-1)maxens3maxens*maxens2 &
	3124	+(iedt-1)maxensmaxens3
	3125	! ne=2, cap=100
	3126	!
	3127	nall2=(iens-1)maxens3maxens*maxens2 &
	3128	+(iedt-1)maxensmaxens3 &
	3129	+(2-1)*maxens3
	3130	xf_ens(i,j,nall+4) = xf_ens(i,j,nall2+4)
	3131	xf_ens(i,j,nall+5) =xf_ens(i,j,nall2+5)
	3132	xf_ens(i,j,nall+6) =xf_ens(i,j,nall2+6)
	3133	xf_ens(i,j,nall+14) =xf_ens(i,j,nall2+14)
	3134	xf_ens(i,j,nall+7) =xf_ens(i,j,nall2+7)
	3135	xf_ens(i,j,nall+8) =xf_ens(i,j,nall2+8)
	3136	xf_ens(i,j,nall+9) =xf_ens(i,j,nall2+9)
	3137	xf_ens(i,j,nall+15) =xf_ens(i,j,nall2+15)
	3138	xf_ens(i,j,nall+10)=xf_ens(i,j,nall2+10)
	3139	xf_ens(i,j,nall+11)=xf_ens(i,j,nall2+11)
	3140	xf_ens(i,j,nall+12)=xf_ens(i,j,nall2+12)
	3141	go to 100
	3142	endif
	3143	elseif(ierr(i).ne.20.and.ierr(i).ne.0)then
	3144	do n=1,ensdim
	3145	xf_ens(i,j,n)=0.
	3146	massfln(i,j,n)=0.
	3147	enddo
	3148	endif
	3149	100 continue
	3150
	3151	END SUBROUTINE cup_forcing_ens_3d
	3152
	3153
	3154	SUBROUTINE cup_kbcon(cap_inc,iloop,k22,kbcon,he_cup,hes_cup, &
	3155	ierr,kbmax,p_cup,cap_max, &
	3156	itf,jtf,ktf, &
	3157	its,ite, jts,jte, kts,kte )
	3158
	3159	IMPLICIT NONE
	3160	!
	3161
	3162	! only local wrf dimensions are need as of now in this routine
	3163
	3164	integer &
	3165	,intent (in ) :: &
	3166	itf,jtf,ktf, &
	3167	its,ite, jts,jte, kts,kte
	3168	!
	3169	!
	3170	!
	3171	! ierr error value, maybe modified in this routine
	3172	!
	3173	real, dimension (its:ite,kts:kte) &
	3174	,intent (in ) :: &
	3175	he_cup,hes_cup,p_cup
	3176	real, dimension (its:ite) &
	3177	,intent (in ) :: &
	3178	cap_max,cap_inc
	3179	integer, dimension (its:ite) &
	3180	,intent (in ) :: &
	3181	kbmax
	3182	integer, dimension (its:ite) &
	3183	,intent (inout) :: &
	3184	kbcon,k22,ierr
	3185	integer &
	3186	,intent (in ) :: &
	3187	iloop
	3188	!
	3189	! local variables in this routine
	3190	!
	3191
	3192	integer :: &
	3193	i
	3194	real :: &
	3195	pbcdif,plus
	3196	!
	3197	!--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
	3198	!
	3199	DO 27 i=its,itf
	3200	kbcon(i)=1
	3201	IF(ierr(I).ne.0)GO TO 27
	3202	KBCON(I)=K22(I)
	3203	GO TO 32
	3204	31 CONTINUE
	3205	KBCON(I)=KBCON(I)+1
	3206	IF(KBCON(I).GT.KBMAX(i)+2)THEN
	3207	if(iloop.lt.4)ierr(i)=3
	3208	! if(iloop.lt.4)ierr(i)=997
	3209	GO TO 27
	3210	ENDIF
	3211	32 CONTINUE
	3212	IF(HE_cup(I,K22(I)).LT.HES_cup(I,KBCON(I)))GO TO 31
	3213
	3214	! cloud base pressure and max moist static energy pressure
	3215	! i.e., the depth (in mb) of the layer of negative buoyancy
	3216	if(KBCON(I)-K22(I).eq.1)go to 27
	3217	PBCDIF=-P_cup(I,KBCON(I))+P_cup(I,K22(I))
	3218	plus=max(25.,cap_max(i)-float(iloop-1)*cap_inc(i))
	3219	if(iloop.eq.4)plus=cap_max(i)
	3220	IF(PBCDIF.GT.plus)THEN
	3221	K22(I)=K22(I)+1
	3222	KBCON(I)=K22(I)
	3223	GO TO 32
	3224	ENDIF
	3225	27 CONTINUE
	3226
	3227	END SUBROUTINE cup_kbcon
	3228
	3229
	3230	SUBROUTINE cup_ktop(ilo,dby,kbcon,ktop,ierr, &
	3231	itf,jtf,ktf, &
	3232	its,ite, jts,jte, kts,kte )
	3233
	3234	IMPLICIT NONE
	3235	!
	3236	! on input
	3237	!
	3238
	3239	! only local wrf dimensions are need as of now in this routine
	3240
	3241	integer &
	3242	,intent (in ) :: &
	3243	itf,jtf,ktf, &
	3244	its,ite, jts,jte, kts,kte
	3245	! dby = buoancy term
	3246	! ktop = cloud top (output)
	3247	! ilo = flag
	3248	! ierr error value, maybe modified in this routine
	3249	!
	3250	real, dimension (its:ite,kts:kte) &
	3251	,intent (inout) :: &
	3252	dby
	3253	integer, dimension (its:ite) &
	3254	,intent (in ) :: &
	3255	kbcon
	3256	integer &
	3257	,intent (in ) :: &
	3258	ilo
	3259	integer, dimension (its:ite) &
	3260	,intent (out ) :: &
	3261	ktop
	3262	integer, dimension (its:ite) &
	3263	,intent (inout) :: &
	3264	ierr
	3265	!
	3266	! local variables in this routine
	3267	!
	3268
	3269	integer :: &
	3270	i,k
	3271	!
	3272	DO 42 i=its,itf
	3273	ktop(i)=1
	3274	IF(ierr(I).EQ.0)then
	3275	DO 40 K=KBCON(I)+1,ktf-1
	3276	IF(DBY(I,K).LE.0.)THEN
	3277	KTOP(I)=K-1
	3278	GO TO 41
	3279	ENDIF
	3280	40 CONTINUE
	3281	if(ilo.eq.1)ierr(i)=5
	3282	! if(ilo.eq.2)ierr(i)=998
	3283	GO TO 42
	3284	41 CONTINUE
	3285	do k=ktop(i)+1,ktf
	3286	dby(i,k)=0.
	3287	enddo
	3288	endif
	3289	42 CONTINUE
	3290
	3291	END SUBROUTINE cup_ktop
	3292
	3293
	3294	SUBROUTINE cup_MAXIMI(ARRAY,KS,KE,MAXX,ierr, &
	3295	itf,jtf,ktf, &
	3296	its,ite, jts,jte, kts,kte )
	3297
	3298	IMPLICIT NONE
	3299	!
	3300	! on input
	3301	!
	3302
	3303	! only local wrf dimensions are need as of now in this routine
	3304
	3305	integer &
	3306	,intent (in ) :: &
	3307	itf,jtf,ktf, &
	3308	its,ite, jts,jte, kts,kte
	3309	! array input array
	3310	! x output array with return values
	3311	! kt output array of levels
	3312	! ks,kend check-range
	3313	real, dimension (its:ite,kts:kte) &
	3314	,intent (in ) :: &
	3315	array
	3316	integer, dimension (its:ite) &
	3317	,intent (in ) :: &
	3318	ierr,ke
	3319	integer &
	3320	,intent (in ) :: &
	3321	ks
	3322	integer, dimension (its:ite) &
	3323	,intent (out ) :: &
	3324	maxx
	3325	real, dimension (its:ite) :: &
	3326	x
	3327	real :: &
	3328	xar
	3329	integer :: &
	3330	i,k
	3331
	3332	DO 200 i=its,itf
	3333	MAXX(I)=KS
	3334	if(ierr(i).eq.0)then
	3335	X(I)=ARRAY(I,KS)
	3336	!
	3337	DO 100 K=KS,KE(i)
	3338	XAR=ARRAY(I,K)
	3339	IF(XAR.GE.X(I)) THEN
	3340	X(I)=XAR
	3341	MAXX(I)=K
	3342	ENDIF
	3343	100 CONTINUE
	3344	endif
	3345	200 CONTINUE
	3346
	3347	END SUBROUTINE cup_MAXIMI
	3348
	3349
	3350	SUBROUTINE cup_minimi(ARRAY,KS,KEND,KT,ierr, &
	3351	itf,jtf,ktf, &
	3352	its,ite, jts,jte, kts,kte )
	3353
	3354	IMPLICIT NONE
	3355	!
	3356	! on input
	3357	!
	3358
	3359	! only local wrf dimensions are need as of now in this routine
	3360
	3361	integer &
	3362	,intent (in ) :: &
	3363	itf,jtf,ktf, &
	3364	its,ite, jts,jte, kts,kte
	3365	! array input array
	3366	! x output array with return values
	3367	! kt output array of levels
	3368	! ks,kend check-range
	3369	real, dimension (its:ite,kts:kte) &
	3370	,intent (in ) :: &
	3371	array
	3372	integer, dimension (its:ite) &
	3373	,intent (in ) :: &
	3374	ierr,ks,kend
	3375	integer, dimension (its:ite) &
	3376	,intent (out ) :: &
	3377	kt
	3378	real, dimension (its:ite) :: &
	3379	x
	3380	integer :: &
	3381	i,k,kstop
	3382
	3383	DO 200 i=its,itf
	3384	KT(I)=KS(I)
	3385	if(ierr(i).eq.0)then
	3386	X(I)=ARRAY(I,KS(I))
	3387	KSTOP=MAX(KS(I)+1,KEND(I))
	3388	!
	3389	DO 100 K=KS(I)+1,KSTOP
	3390	IF(ARRAY(I,K).LT.X(I)) THEN
	3391	X(I)=ARRAY(I,K)
	3392	KT(I)=K
	3393	ENDIF
	3394	100 CONTINUE
	3395	endif
	3396	200 CONTINUE
	3397
	3398	END SUBROUTINE cup_MINIMI
	3399
	3400
	3401	SUBROUTINE cup_output_ens_3d(xf_ens,ierr,dellat,dellaq,dellaqc, &
	3402	subt_ens,subq_ens,subt,subq,outtem,outq,outqc, &
	3403	zu,sub_mas,pre,pw,xmb,ktop, &
	3404	j,name,nx,nx2,iens,ierr2,ierr3,pr_ens, &
	3405	maxens3,ensdim,massfln, &
	3406	APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
	3407	APR_CAPMA,APR_CAPME,APR_CAPMI,closure_n,xland1, &
	3408	itf,jtf,ktf, &
	3409	its,ite, jts,jte, kts,kte)
	3410
	3411	IMPLICIT NONE
	3412	!
	3413	! on input
	3414	!
	3415
	3416	! only local wrf dimensions are need as of now in this routine
	3417
	3418	integer &
	3419	,intent (in ) :: &
	3420	itf,jtf,ktf, &
	3421	its,ite, jts,jte, kts,kte
	3422	integer, intent (in ) :: &
	3423	j,ensdim,nx,nx2,iens,maxens3
	3424	! xf_ens = ensemble mass fluxes
	3425	! pr_ens = precipitation ensembles
	3426	! dellat = change of temperature per unit mass flux of cloud ensemble
	3427	! dellaq = change of q per unit mass flux of cloud ensemble
	3428	! dellaqc = change of qc per unit mass flux of cloud ensemble
	3429	! outtem = output temp tendency (per s)
	3430	! outq = output q tendency (per s)
	3431	! outqc = output qc tendency (per s)
	3432	! pre = output precip
	3433	! xmb = total base mass flux
	3434	! xfac1 = correction factor
	3435	! pw = pw -epsilon*pd (ensemble dependent)
	3436	! ierr error value, maybe modified in this routine
	3437	!
	3438	real, dimension (its:ite,jts:jte,1:ensdim) &
	3439	,intent (inout) :: &
	3440	xf_ens,pr_ens,massfln
	3441	real, dimension (its:ite,jts:jte) &
	3442	,intent (inout) :: &
	3443	APR_GR,APR_W,APR_MC,APR_ST,APR_AS,APR_CAPMA, &
	3444	APR_CAPME,APR_CAPMI
	3445
	3446	real, dimension (its:ite,kts:kte) &
	3447	,intent (out ) :: &
	3448	outtem,outq,outqc,subt,subq,sub_mas
	3449	real, dimension (its:ite,kts:kte) &
	3450	,intent (in ) :: &
	3451	zu
	3452	real, dimension (its:ite) &
	3453	,intent (out ) :: &
	3454	pre,xmb
	3455	real, dimension (its:ite) &
	3456	,intent (inout ) :: &
	3457	closure_n,xland1
	3458	real, dimension (its:ite,kts:kte,1:nx) &
	3459	,intent (in ) :: &
	3460	subt_ens,subq_ens,dellat,dellaqc,dellaq,pw
	3461	integer, dimension (its:ite) &
	3462	,intent (in ) :: &
	3463	ktop
	3464	integer, dimension (its:ite) &
	3465	,intent (inout) :: &
	3466	ierr,ierr2,ierr3
	3467	!
	3468	! local variables in this routine
	3469	!
	3470
	3471	integer :: &
	3472	i,k,n,ncount
	3473	real :: &
	3474	outtes,ddtes,dtt,dtq,dtqc,dtpw,tuning,prerate,clos_wei,xmbhelp
	3475	real :: &
	3476	dtts,dtqs
	3477	real, dimension (its:ite) :: &
	3478	xfac1,xfac2
	3479	real, dimension (its:ite):: &
	3480	xmb_ske,xmb_ave,xmb_std,xmb_cur,xmbweight
	3481	real, dimension (its:ite):: &
	3482	pr_ske,pr_ave,pr_std,pr_cur
	3483	real, dimension (its:ite,jts:jte):: &
	3484	pr_gr,pr_w,pr_mc,pr_st,pr_as,pr_capma, &
	3485	pr_capme,pr_capmi
	3486	real, dimension (5) :: weight,wm,wm1,wm2,wm3
	3487	real, dimension (its:ite,5) :: xmb_w
	3488
	3489	!
	3490	character (), intent (in) :: &
	3491	name
	3492	!
	3493	weight(1) = -999. !this will turn off weights
	3494	wm(1)=-999.
	3495
	3496	tuning=0.
	3497	!
	3498	!
	3499	DO k=kts,ktf
	3500	do i=its,itf
	3501	outtem(i,k)=0.
	3502	outq(i,k)=0.
	3503	outqc(i,k)=0.
	3504	subt(i,k)=0.
	3505	subq(i,k)=0.
	3506	sub_mas(i,k)=0.
	3507	enddo
	3508	enddo
	3509	do i=its,itf
	3510	pre(i)=0.
	3511	xmb(i)=0.
	3512	xfac1(i)=0.
	3513	xfac2(i)=0.
	3514	xmbweight(i)=1.
	3515	enddo
	3516	do i=its,itf
	3517	IF(ierr(i).eq.0)then
	3518	do n=(iens-1)nxnx2maxens3+1,iensnxnx2maxens3
	3519	if(pr_ens(i,j,n).le.0.)then
	3520	xf_ens(i,j,n)=0.
	3521	endif
	3522	enddo
	3523	endif
	3524	enddo
	3525	!
	3526	!--- calculate ensemble average mass fluxes
	3527	!
	3528	call massflx_stats(xf_ens,ensdim,nx2,nx,maxens3, &
	3529	xmb_ave,xmb_std,xmb_cur,xmb_ske,j,ierr,1, &
	3530	APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
	3531	APR_CAPMA,APR_CAPME,APR_CAPMI, &
	3532	pr_gr,pr_w,pr_mc,pr_st,pr_as, &
	3533	pr_capma,pr_capme,pr_capmi, &
	3534	itf,jtf,ktf, &
	3535	its,ite, jts,jte, kts,kte )
	3536	xmb_w=0.
	3537	call massflx_stats(pr_ens,ensdim,nx2,nx,maxens3, &
	3538	pr_ave,pr_std,pr_cur,pr_ske,j,ierr,2, &
	3539	APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
	3540	APR_CAPMA,APR_CAPME,APR_CAPMI, &
	3541	pr_gr,pr_w,pr_mc,pr_st,pr_as, &
	3542	pr_capma,pr_capme,pr_capmi, &
	3543	itf,jtf,ktf, &
	3544	its,ite, jts,jte, kts,kte )
	3545	!
	3546	!-- now do feedback
	3547	!
	3548	ddtes=100.
	3549	do i=its,itf
	3550	if(ierr(i).eq.0)then
	3551	if(xmb_ave(i).le.0.)then
	3552	ierr(i)=13
	3553	xmb_ave(i)=0.
	3554	endif
	3555	xmb(i)=max(.1xmb_ave(i),xmb_ave(i)-tuningxmb_std(i))
	3556	! --- Now use proper count of how many closures were actually
	3557	! used in cup_forcing_ens (including screening of some
	3558	! closures over water) to properly normalize xmb
	3559	clos_wei=16./max(1.,closure_n(i))
	3560	if (xland1(i).lt.0.5)xmb(i)=xmb(i)*clos_wei
	3561	if(xmb(i).eq.0.)then
	3562	ierr(i)=19
	3563	endif
	3564	if(xmb(i).gt.100.)then
	3565	ierr(i)=19
	3566	endif
	3567	xfac1(i)=xmb(i)
	3568	xfac2(i)=xmb(i)
	3569
	3570	endif
	3571	! if(weight(1).lt.-100.)xfac1(i)=xmb_ave(i)
	3572	! if(weight(1).lt.-100.)xfac2(i)=xmb_ave(i)
	3573	ENDDO
	3574	DO k=kts,ktf
	3575	do i=its,itf
	3576	dtt=0.
	3577	dtts=0.
	3578	dtq=0.
	3579	dtqs=0.
	3580	dtqc=0.
	3581	dtpw=0.
	3582	IF(ierr(i).eq.0.and.k.le.ktop(i))then
	3583	do n=1,nx
	3584	dtt=dtt+dellat(i,k,n)
	3585	dtts=dtts+subt_ens(i,k,n)
	3586	dtq=dtq+dellaq(i,k,n)
	3587	dtqs=dtqs+subq_ens(i,k,n)
	3588	dtqc=dtqc+dellaqc(i,k,n)
	3589	dtpw=dtpw+pw(i,k,n)
	3590	enddo
	3591	OUTTEM(I,K)=XMB(I)*dtt/float(nx)
	3592	SUBT(I,K)=XMB(I)*dtts/float(nx)
	3593	OUTQ(I,K)=XMB(I)*dtq/float(nx)
	3594	SUBQ(I,K)=XMB(I)*dtqs/float(nx)
	3595	OUTQC(I,K)=XMB(I)*dtqc/float(nx)
	3596	PRE(I)=PRE(I)+XMB(I)*dtpw/float(nx)
	3597	sub_mas(i,k)=zu(i,k)*xmb(i)
	3598	endif
	3599	enddo
	3600	enddo
	3601
	3602	do i=its,itf
	3603	if(ierr(i).eq.0)then
	3604	do k=(iens-1)nxnx2maxens3+1,iensnxnx2maxens3
	3605	massfln(i,j,k)=massfln(i,j,k)*xfac1(i)
	3606	xf_ens(i,j,k)=xf_ens(i,j,k)*xfac1(i)
	3607	enddo
	3608	endif
	3609	ENDDO
	3610
	3611	END SUBROUTINE cup_output_ens_3d
	3612
	3613
	3614	SUBROUTINE cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
	3615	kbcon,ktop,ierr, &
	3616	itf,jtf,ktf, &
	3617	its,ite, jts,jte, kts,kte )
	3618
	3619	IMPLICIT NONE
	3620	!
	3621	! on input
	3622	!
	3623
	3624	! only local wrf dimensions are need as of now in this routine
	3625
	3626	integer &
	3627	,intent (in ) :: &
	3628	itf,jtf,ktf, &
	3629	its,ite, jts,jte, kts,kte
	3630	! aa0 cloud work function
	3631	! gamma_cup = gamma on model cloud levels
	3632	! t_cup = temperature (Kelvin) on model cloud levels
	3633	! dby = buoancy term
	3634	! zu= normalized updraft mass flux
	3635	! z = heights of model levels
	3636	! ierr error value, maybe modified in this routine
	3637	!
	3638	real, dimension (its:ite,kts:kte) &
	3639	,intent (in ) :: &
	3640	z,zu,gamma_cup,t_cup,dby
	3641	integer, dimension (its:ite) &
	3642	,intent (in ) :: &
	3643	kbcon,ktop
	3644	!
	3645	! input and output
	3646	!
	3647
	3648
	3649	integer, dimension (its:ite) &
	3650	,intent (inout) :: &
	3651	ierr
	3652	real, dimension (its:ite) &
	3653	,intent (out ) :: &
	3654	aa0
	3655	!
	3656	! local variables in this routine
	3657	!
	3658
	3659	integer :: &
	3660	i,k
	3661	real :: &
	3662	dz,da
	3663	!
	3664	do i=its,itf
	3665	aa0(i)=0.
	3666	enddo
	3667	DO 100 k=kts+1,ktf
	3668	DO 100 i=its,itf
	3669	IF(ierr(i).ne.0)GO TO 100
	3670	IF(K.LE.KBCON(I))GO TO 100
	3671	IF(K.Gt.KTOP(I))GO TO 100
	3672	DZ=Z(I,K)-Z(I,K-1)
	3673	da=zu(i,k)DZ(9.81/(1004.*( &
	3674	(T_cup(I,K)))))*DBY(I,K-1)/ &
	3675	(1.+GAMMA_CUP(I,K))
	3676	IF(K.eq.KTOP(I).and.da.le.0.)go to 100
	3677	AA0(I)=AA0(I)+da
	3678	if(aa0(i).lt.0.)aa0(i)=0.
	3679	100 continue
	3680
	3681	END SUBROUTINE cup_up_aa0
	3682
	3683
	3684	SUBROUTINE cup_up_he(k22,hkb,z_cup,cd,entr,he_cup,hc, &
	3685	kbcon,ierr,dby,he,hes_cup, &
	3686	itf,jtf,ktf, &
	3687	its,ite, jts,jte, kts,kte )
	3688
	3689	IMPLICIT NONE
	3690	!
	3691	! on input
	3692	!
	3693
	3694	! only local wrf dimensions are need as of now in this routine
	3695
	3696	integer &
	3697	,intent (in ) :: &
	3698	itf,jtf,ktf, &
	3699	its,ite, jts,jte, kts,kte
	3700	! hc = cloud moist static energy
	3701	! hkb = moist static energy at originating level
	3702	! he = moist static energy on model levels
	3703	! he_cup = moist static energy on model cloud levels
	3704	! hes_cup = saturation moist static energy on model cloud levels
	3705	! dby = buoancy term
	3706	! cd= detrainment function
	3707	! z_cup = heights of model cloud levels
	3708	! entr = entrainment rate
	3709	!
	3710	real, dimension (its:ite,kts:kte) &
	3711	,intent (in ) :: &
	3712	he,he_cup,hes_cup,z_cup,cd
	3713	! entr= entrainment rate
	3714	real &
	3715	,intent (in ) :: &
	3716	entr
	3717	integer, dimension (its:ite) &
	3718	,intent (in ) :: &
	3719	kbcon,k22
	3720	!
	3721	! input and output
	3722	!
	3723
	3724	! ierr error value, maybe modified in this routine
	3725
	3726	integer, dimension (its:ite) &
	3727	,intent (inout) :: &
	3728	ierr
	3729
	3730	real, dimension (its:ite,kts:kte) &
	3731	,intent (out ) :: &
	3732	hc,dby
	3733	real, dimension (its:ite) &
	3734	,intent (out ) :: &
	3735	hkb
	3736	!
	3737	! local variables in this routine
	3738	!
	3739
	3740	integer :: &
	3741	i,k
	3742	real :: &
	3743	dz
	3744	!
	3745	!--- moist static energy inside cloud
	3746	!
	3747	do k=kts,ktf
	3748	do i=its,itf
	3749	hc(i,k)=0.
	3750	DBY(I,K)=0.
	3751	enddo
	3752	enddo
	3753	do i=its,itf
	3754	hkb(i)=0.
	3755	enddo
	3756	do i=its,itf
	3757	if(ierr(i).eq.0.)then
	3758	hkb(i)=he_cup(i,k22(i))
	3759	do k=1,k22(i)
	3760	hc(i,k)=he_cup(i,k)
	3761	! DBY(I,K)=0.
	3762	enddo
	3763	do k=k22(i),kbcon(i)-1
	3764	hc(i,k)=hkb(i)
	3765	! DBY(I,K)=0.
	3766	enddo
	3767	k=kbcon(i)
	3768	hc(i,k)=hkb(i)
	3769	DBY(I,Kbcon(i))=Hkb(I)-HES_cup(I,K)
	3770	endif
	3771	enddo
	3772	do k=kts+1,ktf
	3773	do i=its,itf
	3774	if(k.gt.kbcon(i).and.ierr(i).eq.0.)then
	3775	DZ=Z_cup(i,K)-Z_cup(i,K-1)
	3776	HC(i,K)=(HC(i,K-1)(1.-.5CD(i,K)DZ)+entr &
	3777	DZHE(i,K-1))/(1.+entrDZ-.5cd(i,k)dz)
	3778	DBY(I,K)=HC(I,K)-HES_cup(I,K)
	3779	endif
	3780	enddo
	3781
	3782	enddo
	3783
	3784	END SUBROUTINE cup_up_he
	3785
	3786
	3787	SUBROUTINE cup_up_moisture(ierr,z_cup,qc,qrc,pw,pwav, &
	3788	kbcon,ktop,cd,dby,mentr_rate,clw_all, &
	3789	q,GAMMA_cup,zu,qes_cup,k22,qe_cup,xl, &
	3790	itf,jtf,ktf, &
	3791	its,ite, jts,jte, kts,kte )
	3792
	3793	IMPLICIT NONE
	3794	!
	3795	! on input
	3796	!
	3797
	3798	! only local wrf dimensions are need as of now in this routine
	3799
	3800	integer &
	3801	,intent (in ) :: &
	3802	itf,jtf,ktf, &
	3803	its,ite, jts,jte, kts,kte
	3804	! cd= detrainment function
	3805	! q = environmental q on model levels
	3806	! qe_cup = environmental q on model cloud levels
	3807	! qes_cup = saturation q on model cloud levels
	3808	! dby = buoancy term
	3809	! cd= detrainment function
	3810	! zu = normalized updraft mass flux
	3811	! gamma_cup = gamma on model cloud levels
	3812	! mentr_rate = entrainment rate
	3813	!
	3814	real, dimension (its:ite,kts:kte) &
	3815	,intent (in ) :: &
	3816	q,zu,gamma_cup,qe_cup,dby,qes_cup,z_cup,cd
	3817	! entr= entrainment rate
	3818	real &
	3819	,intent (in ) :: &
	3820	mentr_rate,xl
	3821	integer, dimension (its:ite) &
	3822	,intent (in ) :: &
	3823	kbcon,ktop,k22
	3824	!
	3825	! input and output
	3826	!
	3827
	3828	! ierr error value, maybe modified in this routine
	3829
	3830	integer, dimension (its:ite) &
	3831	,intent (inout) :: &
	3832	ierr
	3833	! qc = cloud q (including liquid water) after entrainment
	3834	! qrch = saturation q in cloud
	3835	! qrc = liquid water content in cloud after rainout
	3836	! pw = condensate that will fall out at that level
	3837	! pwav = totan normalized integrated condensate (I1)
	3838	! c0 = conversion rate (cloud to rain)
	3839
	3840	real, dimension (its:ite,kts:kte) &
	3841	,intent (out ) :: &
	3842	qc,qrc,pw,clw_all
	3843	real, dimension (its:ite) &
	3844	,intent (out ) :: &
	3845	pwav
	3846	!
	3847	! local variables in this routine
	3848	!
	3849
	3850	integer :: &
	3851	iall,i,k
	3852	real :: &
	3853	dh,qrch,c0,dz,radius
	3854	!
	3855	iall=0
	3856	c0=.002
	3857	!
	3858	!--- no precip for small clouds
	3859	!
	3860	if(mentr_rate.gt.0.)then
	3861	radius=.2/mentr_rate
	3862	if(radius.lt.900.)c0=0.
	3863	! if(radius.lt.900.)iall=0
	3864	endif
	3865	do i=its,itf
	3866	pwav(i)=0.
	3867	enddo
	3868	do k=kts,ktf
	3869	do i=its,itf
	3870	pw(i,k)=0.
	3871	qc(i,k)=0.
	3872	if(ierr(i).eq.0)qc(i,k)=qes_cup(i,k)
	3873	clw_all(i,k)=0.
	3874	qrc(i,k)=0.
	3875	enddo
	3876	enddo
	3877	do i=its,itf
	3878	if(ierr(i).eq.0.)then
	3879	do k=k22(i),kbcon(i)-1
	3880	qc(i,k)=qe_cup(i,k22(i))
	3881	enddo
	3882	endif
	3883	enddo
	3884
	3885	DO 100 k=kts+1,ktf
	3886	DO 100 i=its,itf
	3887	IF(ierr(i).ne.0)GO TO 100
	3888	IF(K.Lt.KBCON(I))GO TO 100
	3889	IF(K.Gt.KTOP(I))GO TO 100
	3890	DZ=Z_cup(i,K)-Z_cup(i,K-1)
	3891	!
	3892	!------ 1. steady state plume equation, for what could
	3893	!------ be in cloud without condensation
	3894	!
	3895	!
	3896	QC(i,K)=(QC(i,K-1)(1.-.5CD(i,K)DZ)+mentr_rate &
	3897	DZQ(i,K-1))/(1.+mentr_rateDZ-.5cd(i,k)dz)
	3898	!
	3899	!--- saturation in cloud, this is what is allowed to be in it
	3900	!
	3901	QRCH=QES_cup(I,K)+(1./XL)*(GAMMA_cup(i,k) &
	3902	/(1.+GAMMA_cup(i,k)))*DBY(I,K)
	3903	!
	3904	!------- LIQUID WATER CONTENT IN cloud after rainout
	3905	!
	3906	clw_all(i,k)=QC(I,K)-QRCH
	3907	QRC(I,K)=(QC(I,K)-QRCH)/(1.+C0*DZ)
	3908	if(qrc(i,k).lt.0.)then
	3909	qrc(i,k)=0.
	3910	endif
	3911	!
	3912	!------- 3.Condensation
	3913	!
	3914	PW(i,k)=c0dzQRC(I,K)*zu(i,k)
	3915	if(iall.eq.1)then
	3916	qrc(i,k)=0.
	3917	pw(i,k)=(QC(I,K)-QRCH)*zu(i,k)
	3918	if(pw(i,k).lt.0.)pw(i,k)=0.
	3919	endif
	3920	!
	3921	!----- set next level
	3922	!
	3923	QC(I,K)=QRC(I,K)+qrch
	3924	!
	3925	!--- integrated normalized ondensate
	3926	!
	3927	PWAV(I)=PWAV(I)+PW(I,K)
	3928	100 CONTINUE
	3929
	3930	END SUBROUTINE cup_up_moisture
	3931
	3932
	3933	SUBROUTINE cup_up_nms(zu,z_cup,entr,cd,kbcon,ktop,ierr,k22, &
	3934	itf,jtf,ktf, &
	3935	its,ite, jts,jte, kts,kte )
	3936
	3937	IMPLICIT NONE
	3938
	3939	!
	3940	! on input
	3941	!
	3942
	3943	! only local wrf dimensions are need as of now in this routine
	3944
	3945	integer &
	3946	,intent (in ) :: &
	3947	itf,jtf,ktf, &
	3948	its,ite, jts,jte, kts,kte
	3949	! cd= detrainment function
	3950	real, dimension (its:ite,kts:kte) &
	3951	,intent (in ) :: &
	3952	z_cup,cd
	3953	! entr= entrainment rate
	3954	real &
	3955	,intent (in ) :: &
	3956	entr
	3957	integer, dimension (its:ite) &
	3958	,intent (in ) :: &
	3959	kbcon,ktop,k22
	3960	!
	3961	! input and output
	3962	!
	3963
	3964	! ierr error value, maybe modified in this routine
	3965
	3966	integer, dimension (its:ite) &
	3967	,intent (inout) :: &
	3968	ierr
	3969	! zu is the normalized mass flux
	3970
	3971	real, dimension (its:ite,kts:kte) &
	3972	,intent (out ) :: &
	3973	zu
	3974	!
	3975	! local variables in this routine
	3976	!
	3977
	3978	integer :: &
	3979	i,k
	3980	real :: &
	3981	dz
	3982	!
	3983	! initialize for this go around
	3984	!
	3985	do k=kts,ktf
	3986	do i=its,itf
	3987	zu(i,k)=0.
	3988	enddo
	3989	enddo
	3990	!
	3991	! do normalized mass budget
	3992	!
	3993	do i=its,itf
	3994	IF(ierr(I).eq.0)then
	3995	do k=k22(i),kbcon(i)
	3996	zu(i,k)=1.
	3997	enddo
	3998	DO K=KBcon(i)+1,KTOP(i)
	3999	DZ=Z_cup(i,K)-Z_cup(i,K-1)
	4000	ZU(i,K)=ZU(i,K-1)(1.+(entr-cd(i,k))DZ)
	4001	enddo
	4002	endif
	4003	enddo
	4004
	4005	END SUBROUTINE cup_up_nms
	4006
	4007	!====================================================================
	4008	SUBROUTINE g3init(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
	4009	MASS_FLUX,cp,restart, &
	4010	P_QC,P_QI,P_FIRST_SCALAR, &
	4011	RTHFTEN, RQVFTEN, &
	4012	APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
	4013	APR_CAPMA,APR_CAPME,APR_CAPMI, &
	4014	cugd_tten,cugd_ttens,cugd_qvten, &
	4015	cugd_qvtens,cugd_qcten, &
	4016	allowed_to_read, &
	4017	ids, ide, jds, jde, kds, kde, &
	4018	ims, ime, jms, jme, kms, kme, &
	4019	its, ite, jts, jte, kts, kte )
	4020	!--------------------------------------------------------------------
	4021	IMPLICIT NONE
	4022	!--------------------------------------------------------------------
	4023	LOGICAL , INTENT(IN) :: restart,allowed_to_read
	4024	INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
	4025	ims, ime, jms, jme, kms, kme, &
	4026	its, ite, jts, jte, kts, kte
	4027	INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
	4028	REAL, INTENT(IN) :: cp
	4029
	4030	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
	4031	CUGD_TTEN, &
	4032	CUGD_TTENS, &
	4033	CUGD_QVTEN, &
	4034	CUGD_QVTENS, &
	4035	CUGD_QCTEN
	4036	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
	4037	RTHCUTEN, &
	4038	RQVCUTEN, &
	4039	RQCCUTEN, &
	4040	RQICUTEN
	4041
	4042	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
	4043	RTHFTEN, &
	4044	RQVFTEN
	4045
	4046	REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(OUT) :: &
	4047	APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
	4048	APR_CAPMA,APR_CAPME,APR_CAPMI, &
	4049	MASS_FLUX
	4050
	4051	INTEGER :: i, j, k, itf, jtf, ktf
	4052
	4053	jtf=min0(jte,jde-1)
	4054	ktf=min0(kte,kde-1)
	4055	itf=min0(ite,ide-1)
	4056
	4057	IF(.not.restart)THEN
	4058	DO j=jts,jte
	4059	DO k=kts,kte
	4060	DO i=its,ite
	4061	RTHCUTEN(i,k,j)=0.
	4062	RQVCUTEN(i,k,j)=0.
	4063	ENDDO
	4064	ENDDO
	4065	ENDDO
	4066	DO j=jts,jte
	4067	DO k=kts,kte
	4068	DO i=its,ite
	4069	cugd_tten(i,k,j)=0.
	4070	cugd_ttens(i,k,j)=0.
	4071	cugd_qvten(i,k,j)=0.
	4072	cugd_qvtens(i,k,j)=0.
	4073	ENDDO
	4074	ENDDO
	4075	ENDDO
	4076
	4077	DO j=jts,jtf
	4078	DO k=kts,ktf
	4079	DO i=its,itf
	4080	RTHFTEN(i,k,j)=0.
	4081	RQVFTEN(i,k,j)=0.
	4082	ENDDO
	4083	ENDDO
	4084	ENDDO
	4085
	4086	IF (P_QC .ge. P_FIRST_SCALAR) THEN
	4087	DO j=jts,jtf
	4088	DO k=kts,ktf
	4089	DO i=its,itf
	4090	RQCCUTEN(i,k,j)=0.
	4091	cugd_qcten(i,k,j)=0.
	4092	ENDDO
	4093	ENDDO
	4094	ENDDO
	4095	ENDIF
	4096
	4097	IF (P_QI .ge. P_FIRST_SCALAR) THEN
	4098	DO j=jts,jtf
	4099	DO k=kts,ktf
	4100	DO i=its,itf
	4101	RQICUTEN(i,k,j)=0.
	4102	ENDDO
	4103	ENDDO
	4104	ENDDO
	4105	ENDIF
	4106
	4107	DO j=jts,jtf
	4108	DO i=its,itf
	4109	mass_flux(i,j)=0.
	4110	ENDDO
	4111	ENDDO
	4112
	4113	ENDIF
	4114	DO j=jts,jtf
	4115	DO i=its,itf
	4116	APR_GR(i,j)=0.
	4117	APR_ST(i,j)=0.
	4118	APR_W(i,j)=0.
	4119	APR_MC(i,j)=0.
	4120	APR_AS(i,j)=0.
	4121	APR_CAPMA(i,j)=0.
	4122	APR_CAPME(i,j)=0.
	4123	APR_CAPMI(i,j)=0.
	4124	ENDDO
	4125	ENDDO
	4126
	4127	END SUBROUTINE g3init
	4128
	4129
	4130	SUBROUTINE massflx_stats(xf_ens,ensdim,maxens,maxens2,maxens3, &
	4131	xt_ave,xt_std,xt_cur,xt_ske,j,ierr,itest, &
	4132	APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
	4133	APR_CAPMA,APR_CAPME,APR_CAPMI, &
	4134	pr_gr,pr_w,pr_mc,pr_st,pr_as, &
	4135	pr_capma,pr_capme,pr_capmi, &
	4136	itf,jtf,ktf, &
	4137	its,ite, jts,jte, kts,kte)
	4138
	4139	IMPLICIT NONE
	4140
	4141	integer, intent (in ) :: &
	4142	j,ensdim,maxens3,maxens,maxens2,itest
	4143	INTEGER, INTENT(IN ) :: &
	4144	itf,jtf,ktf, &
	4145	its,ite, jts,jte, kts,kte
	4146
	4147
	4148	real, dimension (its:ite) &
	4149	, intent(inout) :: &
	4150	xt_ave,xt_cur,xt_std,xt_ske
	4151	integer, dimension (its:ite), intent (in) :: &
	4152	ierr
	4153	real, dimension (its:ite,jts:jte,1:ensdim) &
	4154	, intent(in ) :: &
	4155	xf_ens
	4156	real, dimension (its:ite,jts:jte) &
	4157	, intent(inout) :: &
	4158	APR_GR,APR_W,APR_MC,APR_ST,APR_AS, &
	4159	APR_CAPMA,APR_CAPME,APR_CAPMI
	4160	real, dimension (its:ite,jts:jte) &
	4161	, intent(inout) :: &
	4162	pr_gr,pr_w,pr_mc,pr_st,pr_as, &
	4163	pr_capma,pr_capme,pr_capmi
	4164
	4165	!
	4166	! local stuff
	4167	!
	4168	real, dimension (its:ite , 1:maxens3 ) :: &
	4169	x_ave,x_cur,x_std,x_ske
	4170	real, dimension (its:ite , 1:maxens ) :: &
	4171	x_ave_cap
	4172
	4173
	4174	integer, dimension (1:maxens3) :: nc1
	4175	integer :: i,k
	4176	integer :: num,kk,num2,iedt
	4177	real :: a3,a4
	4178
	4179	num=ensdim/maxens3
	4180	num2=ensdim/maxens
	4181	if(itest.eq.1)then
	4182	do i=its,ite
	4183	pr_gr(i,j) = 0.
	4184	pr_w(i,j) = 0.
	4185	pr_mc(i,j) = 0.
	4186	pr_st(i,j) = 0.
	4187	pr_as(i,j) = 0.
	4188	pr_capma(i,j) = 0.
	4189	pr_capme(i,j) = 0.
	4190	pr_capmi(i,j) = 0.
	4191	enddo
	4192	endif
	4193
	4194	do k=1,maxens
	4195	do i=its,ite
	4196	x_ave_cap(i,k)=0.
	4197	enddo
	4198	enddo
	4199	do k=1,maxens3
	4200	do i=its,ite
	4201	x_ave(i,k)=0.
	4202	x_std(i,k)=0.
	4203	x_ske(i,k)=0.
	4204	x_cur(i,k)=0.
	4205	enddo
	4206	enddo
	4207	do i=its,ite
	4208	xt_ave(i)=0.
	4209	xt_std(i)=0.
	4210	xt_ske(i)=0.
	4211	xt_cur(i)=0.
	4212	enddo
	4213	do kk=1,num
	4214	do k=1,maxens3
	4215	do i=its,ite
	4216	if(ierr(i).eq.0)then
	4217	x_ave(i,k)=x_ave(i,k)+xf_ens(i,j,maxens3*(kk-1)+k)
	4218	endif
	4219	enddo
	4220	enddo
	4221	enddo
	4222	do iedt=1,maxens2
	4223	do k=1,maxens
	4224	do kk=1,maxens3
	4225	do i=its,ite
	4226	if(ierr(i).eq.0)then
	4227	x_ave_cap(i,k)=x_ave_cap(i,k) &
	4228	+xf_ens(i,j,maxens3(k-1)+(iedt-1)maxens*maxens3+kk)
	4229	endif
	4230	enddo
	4231	enddo
	4232	enddo
	4233	enddo
	4234	do k=1,maxens
	4235	do i=its,ite
	4236	if(ierr(i).eq.0)then
	4237	x_ave_cap(i,k)=x_ave_cap(i,k)/float(num2)
	4238	endif
	4239	enddo
	4240	enddo
	4241
	4242	do k=1,maxens3
	4243	do i=its,ite
	4244	if(ierr(i).eq.0)then
	4245	x_ave(i,k)=x_ave(i,k)/float(num)
	4246	endif
	4247	enddo
	4248	enddo
	4249	do k=1,maxens3
	4250	do i=its,ite
	4251	if(ierr(i).eq.0)then
	4252	xt_ave(i)=xt_ave(i)+x_ave(i,k)
	4253	endif
	4254	enddo
	4255	enddo
	4256	do i=its,ite
	4257	if(ierr(i).eq.0)then
	4258	xt_ave(i)=xt_ave(i)/float(maxens3)
	4259	endif
	4260	enddo
	4261	!
	4262	!--- now do std, skewness,curtosis
	4263	!
	4264	do kk=1,num
	4265	do k=1,maxens3
	4266	do i=its,ite
	4267	if(ierr(i).eq.0.and.x_ave(i,k).gt.0.)then
	4268	! print ,i,j,k,kk,x_std(i,k),xf_ens(i,j,maxens3(kk-1)+k),x_ave(i,k)
	4269	x_std(i,k)=x_std(i,k)+(xf_ens(i,j,maxens3(kk-1)+k)-x_ave(i,k))*2
	4270	x_ske(i,k)=x_ske(i,k)+(xf_ens(i,j,maxens3(kk-1)+k)-x_ave(i,k))*3
	4271	x_cur(i,k)=x_cur(i,k)+(xf_ens(i,j,maxens3(kk-1)+k)-x_ave(i,k))*4
	4272	endif
	4273	enddo
	4274	enddo
	4275	enddo
	4276	do k=1,maxens3
	4277	do i=its,ite
	4278	if(ierr(i).eq.0.and.xt_ave(i).gt.0.)then
	4279	xt_std(i)=xt_std(i)+(x_ave(i,k)-xt_ave(i))**2
	4280	xt_ske(i)=xt_ske(i)+(x_ave(i,k)-xt_ave(i))**3
	4281	xt_cur(i)=xt_cur(i)+(x_ave(i,k)-xt_ave(i))**4
	4282	endif
	4283	enddo
	4284	enddo
	4285	do k=1,maxens3
	4286	do i=its,ite
	4287	if(ierr(i).eq.0.and.x_std(i,k).gt.0.)then
	4288	x_std(i,k)=x_std(i,k)/float(num)
	4289	a3=max(1.e-6,x_std(i,k))
	4290	x_std(i,k)=sqrt(a3)
	4291	a3=max(1.e-6,x_std(i,k)**3)
	4292	a4=max(1.e-6,x_std(i,k)**4)
	4293	x_ske(i,k)=x_ske(i,k)/float(num)/a3
	4294	x_cur(i,k)=x_cur(i,k)/float(num)/a4
	4295	endif
	4296	! print*,' '
	4297	! print*,'Some statistics at gridpoint i,j, ierr',i,j,ierr(i)
	4298	! print*,'statistics for closure number ',k
	4299	! print*,'Average= ',x_ave(i,k),' Std= ',x_std(i,k)
	4300	! print*,'Skewness= ',x_ske(i,k),' Curtosis= ',x_cur(i,k)
	4301	! print*,' '
	4302
	4303	enddo
	4304	enddo
	4305	do i=its,ite
	4306	if(ierr(i).eq.0.and.xt_std(i).gt.0.)then
	4307	xt_std(i)=xt_std(i)/float(maxens3)
	4308	a3=max(1.e-6,xt_std(i))
	4309	xt_std(i)=sqrt(a3)
	4310	a3=max(1.e-6,xt_std(i)**3)
	4311	a4=max(1.e-6,xt_std(i)**4)
	4312	xt_ske(i)=xt_ske(i)/float(maxens3)/a3
	4313	xt_cur(i)=xt_cur(i)/float(maxens3)/a4
	4314	! print*,' '
	4315	! print*,'Total ensemble independent statistics at i =',i
	4316	! print*,'Average= ',xt_ave(i),' Std= ',xt_std(i)
	4317	! print*,'Skewness= ',xt_ske(i),' Curtosis= ',xt_cur(i)
	4318	! print*,' '
	4319	!
	4320	! first go around: store massflx for different closures/caps
	4321	!
	4322	if(itest.eq.1)then
	4323	pr_gr(i,j) = .25*(x_ave(i,1)+x_ave(i,2)+x_ave(i,3)+x_ave(i,13))
	4324	pr_w(i,j) = .25*(x_ave(i,4)+x_ave(i,5)+x_ave(i,6)+x_ave(i,14))
	4325	pr_mc(i,j) = .25*(x_ave(i,7)+x_ave(i,8)+x_ave(i,9)+x_ave(i,15))
	4326	pr_st(i,j) = .333*(x_ave(i,10)+x_ave(i,11)+x_ave(i,12))
	4327	pr_as(i,j) = x_ave(i,16)
	4328	pr_capma(i,j) = x_ave_cap(i,1)
	4329	pr_capme(i,j) = x_ave_cap(i,2)
	4330	pr_capmi(i,j) = x_ave_cap(i,3)
	4331	!
	4332	! second go around: store preciprates (mm/hour) for different closures/caps
	4333	!
	4334	else if (itest.eq.2)then
	4335	APR_GR(i,j)=.25(x_ave(i,1)+x_ave(i,2)+x_ave(i,3)+x_ave(i,13)) &
	4336	3600.*pr_gr(i,j) +APR_GR(i,j)
	4337	APR_W(i,j)=.25(x_ave(i,4)+x_ave(i,5)+x_ave(i,6)+x_ave(i,14)) &
	4338	3600.*pr_w(i,j) +APR_W(i,j)
	4339	APR_MC(i,j)=.25(x_ave(i,7)+x_ave(i,8)+x_ave(i,9)+x_ave(i,15)) &
	4340	3600.*pr_mc(i,j) +APR_MC(i,j)
	4341	APR_ST(i,j)=.333(x_ave(i,10)+x_ave(i,11)+x_ave(i,12)) &
	4342	3600.*pr_st(i,j) +APR_ST(i,j)
	4343	APR_AS(i,j)=x_ave(i,16)* &
	4344	3600.*pr_as(i,j) +APR_AS(i,j)
	4345	APR_CAPMA(i,j) = x_ave_cap(i,1)* &
	4346	3600.*pr_capma(i,j) +APR_CAPMA(i,j)
	4347	APR_CAPME(i,j) = x_ave_cap(i,2)* &
	4348	3600.*pr_capme(i,j) +APR_CAPME(i,j)
	4349	APR_CAPMI(i,j) = x_ave_cap(i,3)* &
	4350	3600.*pr_capmi(i,j) +APR_CAPMI(i,j)
	4351	endif
	4352	endif
	4353	enddo
	4354
	4355	END SUBROUTINE massflx_stats
	4356
	4357	SUBROUTINE cup_axx(tcrit,kbmax,z1,p,psur,xl,rv,cp,tx,qx,axx,ierr, &
	4358	cap_max,cap_max_increment,entr_rate,mentr_rate,&
	4359	j,itf,jtf,ktf, &
	4360	its,ite, jts,jte, kts,kte,ens4)
	4361	IMPLICIT NONE
	4362	INTEGER, INTENT(IN ) :: &
	4363	j,itf,jtf,ktf, &
	4364	its,ite, jts,jte, kts,kte,ens4
	4365	real, dimension (its:ite,kts:kte,1:ens4) &
	4366	, intent(inout) :: &
	4367	tx,qx
	4368	real, dimension (its:ite,kts:kte) &
	4369	, intent(in) :: &
	4370	p
	4371	real, dimension (its:ite) &
	4372	, intent(in) :: &
	4373	z1,psur,cap_max,cap_max_increment
	4374	real, intent(in) :: &
	4375	tcrit,xl,rv,cp,mentr_rate,entr_rate
	4376	real, dimension (its:ite,1:ens4) &
	4377	, intent(out) :: &
	4378	axx
	4379	integer, dimension (its:ite), intent (in) :: &
	4380	ierr,kbmax
	4381	integer, dimension (its:ite) :: &
	4382	ierrxx,k22xx,kbconxx,ktopxx,kstabm,kstabi
	4383	real, dimension (1:2) :: AE,BE,HT
	4384	real, dimension (its:ite,kts:kte) :: tv
	4385	real :: e,tvbar
	4386	integer n,i,k,iph
	4387	real, dimension (its:ite,kts:kte) :: &
	4388	he,hes,qes,z, &
	4389	qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
	4390	tn_cup, &
	4391	dby,qc,qrcd,pwd,pw,hcd,qcd,dbyd,hc,qrc,zu,zd,cd
	4392
	4393	real, dimension (its:ite) :: &
	4394	AA0,HKB,QKB, &
	4395	PWAV,BU
	4396	do n=1,ens4
	4397	do i=its,ite
	4398	axx(i,n)=0.
	4399	enddo
	4400	enddo
	4401	HT(1)=XL/CP
	4402	HT(2)=2.834E6/CP
	4403	BE(1)=.622*HT(1)/.286
	4404	AE(1)=BE(1)/273.+ALOG(610.71)
	4405	BE(2)=.622*HT(2)/.286
	4406	AE(2)=BE(2)/273.+ALOG(610.71)
	4407	!
	4408	!
	4409	do 100 n=1,ens4
	4410
	4411	do k=kts,ktf
	4412	do i=its,itf
	4413	cd(i,k)=0.1*entr_rate
	4414	enddo
	4415	enddo
	4416
	4417
	4418	do i=its,itf
	4419	ierrxx(i)=ierr(i)
	4420	k22xx(i)=1
	4421	kbconxx(i)=1
	4422	ktopxx(i)=1
	4423	kstabm(i)=ktf-1
	4424	enddo
	4425	DO k=kts,ktf
	4426	do i=its,itf
	4427	if(ierrxx(i).eq.0)then
	4428	IPH=1
	4429	IF(Tx(I,K,n).LE.TCRIT)IPH=2
	4430	E=EXP(AE(IPH)-BE(IPH)/TX(I,K,N))
	4431	QES(I,K)=.622E/(100.P(I,K)-E)
	4432	IF(QES(I,K).LE.1.E-08)QES(I,K)=1.E-08
	4433	IF(Qx(I,K,N).GT.QES(I,K))Qx(I,K,N)=QES(I,K)
	4434	TV(I,K)=Tx(I,K,N)+.608Qx(I,K,N)Tx(I,K,N)
	4435	endif
	4436	enddo
	4437	enddo
	4438	!
	4439	do i=its,itf
	4440	if(ierrxx(i).eq.0)then
	4441	Z(I,KTS)=max(0.,Z1(I))-(ALOG(P(I,KTS))- &
	4442	ALOG(PSUR(I)))287.TV(I,KTS)/9.81
	4443	endif
	4444	enddo
	4445
	4446	! --- calculate heights
	4447	DO K=kts+1,ktf
	4448	do i=its,itf
	4449	if(ierrxx(i).eq.0)then
	4450	TVBAR=.5TV(I,K)+.5TV(I,K-1)
	4451	Z(I,K)=Z(I,K-1)-(ALOG(P(I,K))- &
	4452	ALOG(P(I,K-1)))287.TVBAR/9.81
	4453	endif
	4454	enddo
	4455	enddo
	4456	!
	4457	!--- calculate moist static energy - HE
	4458	! saturated moist static energy - HES
	4459	!
	4460	DO k=kts,ktf
	4461	do i=its,itf
	4462	if(ierrxx(i).eq.0)then
	4463	HE(I,K)=9.81Z(I,K)+1004.Tx(I,K,n)+2.5E06*Qx(I,K,n)
	4464	HES(I,K)=9.81Z(I,K)+1004.Tx(I,K,n)+2.5E06*QES(I,K)
	4465	IF(HE(I,K).GE.HES(I,K))HE(I,K)=HES(I,K)
	4466	endif
	4467	enddo
	4468	enddo
	4469
	4470	! cup levels
	4471	!
	4472	do k=kts+1,ktf
	4473	do i=its,itf
	4474	if(ierrxx(i).eq.0)then
	4475	qes_cup(i,k)=.5*(qes(i,k-1)+qes(i,k))
	4476	q_cup(i,k)=.5*(qx(i,k-1,n)+qx(i,k,n))
	4477	hes_cup(i,k)=.5*(hes(i,k-1)+hes(i,k))
	4478	he_cup(i,k)=.5*(he(i,k-1)+he(i,k))
	4479	if(he_cup(i,k).gt.hes_cup(i,k))he_cup(i,k)=hes_cup(i,k)
	4480	z_cup(i,k)=.5*(z(i,k-1)+z(i,k))
	4481	p_cup(i,k)=.5*(p(i,k-1)+p(i,k))
	4482	t_cup(i,k)=.5*(tx(i,k-1,n)+tx(i,k,n))
	4483	gamma_cup(i,k)=(xl/cp)(xl/(rvt_cup(i,k) &
	4484	t_cup(i,k)))qes_cup(i,k)
	4485	endif
	4486	enddo
	4487	enddo
	4488	do i=its,itf
	4489	if(ierrxx(i).eq.0)then
	4490	qes_cup(i,1)=qes(i,1)
	4491	q_cup(i,1)=qx(i,1,n)
	4492	hes_cup(i,1)=hes(i,1)
	4493	he_cup(i,1)=he(i,1)
	4494	z_cup(i,1)=.5*(z(i,1)+z1(i))
	4495	p_cup(i,1)=.5*(p(i,1)+psur(i))
	4496	t_cup(i,1)=tx(i,1,n)
	4497	gamma_cup(i,1)=xl/cp(xl/(rvt_cup(i,1) &
	4498	t_cup(i,1)))qes_cup(i,1)
	4499	endif
	4500	enddo
	4501	!
	4502	!
	4503	!------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
	4504	!
	4505	CALL cup_MAXIMI(HE_CUP,3,KBMAX,K22XX,ierrxx, &
	4506	itf,jtf,ktf, &
	4507	its,ite, jts,jte, kts,kte)
	4508	DO 36 i=its,itf
	4509	IF(ierrxx(I).eq.0.)THEN
	4510	IF(K22xx(I).GE.KBMAX(i))ierrxx(i)=2
	4511	endif
	4512	36 CONTINUE
	4513	!
	4514	!--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
	4515	!
	4516	call cup_kbcon(cap_max_increment,1,k22xx,kbconxx,he_cup,hes_cup, &
	4517	ierrxx,kbmax,p_cup,cap_max, &
	4518	itf,jtf,ktf, &
	4519	its,ite, jts,jte, kts,kte)
	4520	!
	4521	!--- increase detrainment in stable layers
	4522	!
	4523	CALL cup_minimi(HEs_cup,Kbconxx,kstabm,kstabi,ierrxx, &
	4524	itf,jtf,ktf, &
	4525	its,ite, jts,jte, kts,kte)
	4526	do i=its,itf
	4527	IF(ierrxx(I).eq.0.)THEN
	4528	if(kstabm(i)-1.gt.kstabi(i))then
	4529	do k=kstabi(i),kstabm(i)-1
	4530	cd(i,k)=cd(i,k-1)+1.5*entr_rate
	4531	if(cd(i,k).gt.10.0entr_rate)cd(i,k)=10.0entr_rate
	4532	enddo
	4533	ENDIF
	4534	ENDIF
	4535	ENDDO
	4536	!
	4537	!--- calculate incloud moist static energy
	4538	!
	4539	call cup_up_he(k22xx,hkb,z_cup,cd,mentr_rate,he_cup,hc, &
	4540	kbconxx,ierrxx,dby,he,hes_cup, &
	4541	itf,jtf,ktf, &
	4542	its,ite, jts,jte, kts,kte)
	4543
	4544	!--- DETERMINE CLOUD TOP - KTOP
	4545	!
	4546	call cup_ktop(1,dby,kbconxx,ktopxx,ierrxx, &
	4547	itf,jtf,ktf, &
	4548	its,ite, jts,jte, kts,kte)
	4549	!
	4550	!c--- normalized updraft mass flux profile
	4551	!
	4552	call cup_up_nms(zu,z_cup,mentr_rate,cd,kbconxx,ktopxx,ierrxx,k22xx, &
	4553	itf,jtf,ktf, &
	4554	its,ite, jts,jte, kts,kte)
	4555	!
	4556	!--- calculate workfunctions for updrafts
	4557	!
	4558	call cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
	4559	kbconxx,ktopxx,ierrxx, &
	4560	itf,jtf,ktf, &
	4561	its,ite, jts,jte, kts,kte)
	4562	do i=its,itf
	4563	if(ierrxx(i).eq.0)axx(i,n)=aa0(i)
	4564	enddo
	4565	100 continue
	4566	END SUBROUTINE cup_axx
	4567
	4568	SUBROUTINE conv_grell_spread3d(rthcuten,rqvcuten,rqccuten,raincv, &
	4569	& cugd_avedx,cugd_tten,cugd_qvten,rqicuten,cugd_ttens, &
	4570	& cugd_qvtens,cugd_qcten,pi_phy,moist_qv,pratec,dt,num_tiles,&
	4571	& imomentum,F_QV ,F_QC ,F_QR ,F_QI ,F_QS, &
	4572	& ids, ide, jds, jde, kds, kde, &
	4573	& ips, ipe, jps, jpe, kps, kpe, &
	4574	& ims, ime, jms, jme, kms, kme, &
	4575	& i_start,i_end,j_start,j_end,kts,kte )
	4576
	4577	!
	4578
	4579	INTEGER, INTENT(IN ) :: num_tiles,imomentum
	4580	INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
	4581	& i_start,i_end,j_start,j_end
	4582	INTEGER, INTENT(IN ) :: ids, ide, jds, jde, kds, kde,&
	4583	ims,ime, jms,jme, kms,kme, &
	4584	ips,ipe, jps,jpe, kps,kpe, &
	4585	kts,kte,cugd_avedx
	4586	REAL, DIMENSION (ims:ime,kms:kme,jms:jme), optional,INTENT (INOUT) :: &
	4587	& rthcuten,rqvcuten,rqccuten,rqicuten,cugd_tten, &
	4588	& cugd_qvten,cugd_ttens,cugd_qvtens,cugd_qcten
	4589	REAL, DIMENSION (ims:ime,kms:kme,jms:jme),INTENT (IN) :: &
	4590	moist_qv
	4591	REAL, DIMENSION (ims:ime,kms:kme,jms:jme), INTENT (IN) :: &
	4592	PI_PHY
	4593	REAL, DIMENSION (ims:ime,jms:jme), INTENT (INOUT) :: &
	4594	RAINCV,PRATEC
	4595	REAL, INTENT(IN) :: dt
	4596	INTEGER :: ikk1,ikk2,ikk11,i,j,k,kk,nn,smoothh,smoothv
	4597	INTEGER :: ifs,ife,jfs,jfe,its,ite,jts,jte,ido,jdo,cugd_spread
	4598	LOGICAL :: new
	4599	!
	4600	! Flags relating to the optional tendency arrays declared above
	4601	! Models that carry the optional tendencies will provdide the
	4602	! optional arguments at compile time; these flags all the model
	4603	! to determine at run-time whether a particular tracer is in
	4604	! use or not.
	4605	!
	4606	LOGICAL, OPTIONAL :: &
	4607	F_QV &
	4608	,F_QC &
	4609	,F_QR &
	4610	,F_QI &
	4611	,F_QS
	4612	REAL, DIMENSION (ips-2:ipe+2,kps:kpe,jps-2:jpe+2) :: &
	4613	rthcutent,rqvcutent
	4614	real, dimension (ips-2:ipe+2,jps-2:jpe+2) :: qmem
	4615	real, dimension (ips-1:ipe+1,jps-1:jpe+1) :: smtt,smtq
	4616	real, dimension (kps:kpe) :: conv_trasht,conv_trashq
	4617	REAL :: qmem1,qmem2,qmemf,thresh
	4618	smoothh=1
	4619	smoothv=1
	4620	cugd_spread=cugd_avedx/2
	4621	! SET START AND END POINTS FOR TILES
	4622	!$OMP PARALLEL DO &
	4623	!$OMP PRIVATE ( ij ,ifs,ife,jfs,jfe,its,ite,jts,jte, i,j,k,kk,nn,ikk1,ikk2,ikk11) &
	4624	!$OMP PRIVATE ( ido,jdo,qmemf,qmem1,qmem2,qmem,thresh,conv_trasht,conv_trashq,smtt,smtq)
	4625
	4626	DO ij = 1 , num_tiles
	4627	its = i_start(ij)
	4628	ite = min(i_end(ij),ide-1)
	4629	jts = j_start(ij)
	4630	jte = min(j_end(ij),jde-1)
	4631
	4632	do j=jts-2,jte+2
	4633	do i=its-2,ite+2
	4634	qmem(i,j)=1.
	4635	enddo
	4636	enddo
	4637	do j=jts-1,jte+1
	4638	do i=its-1,ite+1
	4639	smtt(i,j)=0.
	4640	smtq(i,j)=0.
	4641	enddo
	4642	enddo
	4643	do j=jts,jte
	4644	do k=kts,kte
	4645	do i=its,ite
	4646	rthcuten(i,k,j)=0.
	4647	rqvcuten(i,k,j)=0.
	4648	enddo
	4649	enddo
	4650	enddo
	4651	do j=jts-2,jte+2
	4652	do k=kts,kte
	4653	do i=its-2,ite+2
	4654	rthcutent(i,k,j)=0.
	4655	rqvcutent(i,k,j)=0.
	4656	enddo
	4657	enddo
	4658	enddo
	4659	!
	4660	ifs=max(its,ids)
	4661	jfs=max(jts,jds)
	4662	ife=min(ite,ide-1)
	4663	jfe=min(jte,jde-1)
	4664	!
	4665	!
	4666	!
	4667	! prelims finished, now go real for every grid point
	4668	!
	4669	ifs=max(its,ids)
	4670	ife=min(ite,ide-1)
	4671	jfs=max(jts,jds)
	4672	jfe=min(jte,jde-1)
	4673	if(cugd_spread.gt.0.or.smoothh.eq.1)then
	4674	if(its.eq.ips)ifs=max(its-1,ids)
	4675	if(ite.eq.ipe)ife=min(ite+1,ide-1)
	4676	if(jts.eq.jps)jfs=max(jts-1,jds)
	4677	if(jte.eq.jpe)jfe=min(jte+1,jde-1)
	4678	endif
	4679	do j=jfs,jfe
	4680	do i=ifs,ife
	4681	!
	4682	do k=kts,kte
	4683	rthcutent(i,k,j)=cugd_tten(i,k,j)
	4684	rqvcutent(i,k,j)=cugd_qvten(i,k,j)
	4685	enddo
	4686	!
	4687	! for high res run, spread the subsidence
	4688	! this is tricky......only consider grid points where there was no rain,
	4689	! so cugd_tten and such are zero!
	4690	!
	4691	if(cugd_spread.gt.0)then
	4692	do k=kts,kte
	4693	do nn=-1,1,1
	4694	jdo=max(j+nn,jds)
	4695	jdo=min(jdo,jde-1)
	4696	do kk=-1,1,1
	4697	ido=max(i+kk,ids)
	4698	ido=min(ido,ide-1)
	4699	rthcutent(i,k,j)=rthcutent(i,k,j) &
	4700	+qmem(ido,jdo)*cugd_ttens(ido,k,jdo)
	4701	rqvcutent(i,k,j)=rqvcutent(i,k,j) &
	4702	+qmem(ido,jdo)*cugd_qvtens(ido,k,jdo)
	4703	enddo
	4704	enddo
	4705	enddo
	4706	endif
	4707	!
	4708	! end spreading
	4709
	4710	if(cugd_spread.eq.0)then
	4711	do k=kts,kte
	4712	rthcutent(i,k,j)=rthcutent(i,k,j)+cugd_ttens(i,k,j)
	4713	rqvcutent(i,k,j)=rqvcutent(i,k,j)+cugd_qvtens(i,k,j)
	4714	enddo
	4715	endif
	4716	enddo ! end j
	4717	enddo ! end i
	4718	! smooth
	4719	do k=kts,kte
	4720	if(smoothh.eq.0)then
	4721	ifs=max(its,ids+4)
	4722	ife=min(ite,ide-5)
	4723	jfs=max(jts,jds+4)
	4724	jfe=min(jte,jde-5)
	4725	do i=ifs,ife
	4726	do j=jfs,jfe
	4727	rthcuten(i,k,j)=rthcutent(i,k,j)
	4728	rqvcuten(i,k,j)=rqvcutent(i,k,j)
	4729	enddo ! end j
	4730	enddo ! end j
	4731	else if(smoothh.eq.1)then ! smooth
	4732	ifs=max(its,ids)
	4733	ife=min(ite,ide-1)
	4734	jfs=max(jts,jds)
	4735	jfe=min(jte,jde-1)
	4736	! we need an extra row for j (halo comp)
	4737	if(jts.eq.jps)jfs=max(jts-1,jds)
	4738	if(jte.eq.jpe)jfe=min(jte+1,jde-1)
	4739	do i=ifs,ife
	4740	do j=jfs,jfe
	4741	smtt(i,j)=.25(rthcutent(i-1,k,j)+2.rthcutent(i,k,j)+rthcutent(i+1,k,j))
	4742	smtq(i,j)=.25(rqvcutent(i-1,k,j)+2.rqvcutent(i,k,j)+rqvcutent(i+1,k,j))
	4743	enddo ! end j
	4744	enddo ! end j
	4745	ifs=max(its,ids+4)
	4746	ife=min(ite,ide-5)
	4747	jfs=max(jts,jds+4)
	4748	jfe=min(jte,jde-5)
	4749	do i=ifs,ife
	4750	do j=jfs,jfe
	4751	rthcuten(i,k,j)=.25(smtt(i,j-1)+2.smtt(i,j)+smtt(i,j+1))
	4752	rqvcuten(i,k,j)=.25(smtq(i,j-1)+2.smtq(i,j)+smtq(i,j+1))
	4753	enddo ! end j
	4754	enddo ! end i
	4755	endif ! smoothh
	4756
	4757	enddo ! end k
	4758	!
	4759	! check moistening rates
	4760	!
	4761	ifs=max(its,ids+4)
	4762	ife=min(ite,ide-5)
	4763	jfs=max(jts,jds+4)
	4764	jfe=min(jte,jde-5)
	4765	do j=jfs,jfe
	4766	do i=ifs,ife
	4767	qmemf=1.
	4768	thresh=1.e-20
	4769	do k=kts,kte
	4770	if(rqvcuten(i,k,j).lt.0.)then
	4771
	4772	qmem1=moist_qv(i,k,j)+rqvcuten(i,k,j)*dt
	4773	if(qmem1.lt.thresh)then
	4774	qmem1=rqvcuten(i,k,j)
	4775	qmem2=(thresh-moist_qv(i,k,j))/dt
	4776	qmemf=min(qmemf,qmem2/qmem1)
	4777	qmemf=max(0.,qmemf)
	4778	qmemf=min(1.,qmemf)
	4779	endif
	4780
	4781	endif
	4782	enddo
	4783	do k=kts,kte
	4784	rqvcuten(i,k,j)=rqvcuten(i,k,j)*qmemf
	4785	rthcuten(i,k,j)=rthcuten(i,k,j)*qmemf
	4786	enddo
	4787	if(present(rqccuten))then
	4788	if(f_qc) then
	4789	do k=kts,kte
	4790	rqccuten(i,k,j)=rqccuten(i,k,j)*qmemf
	4791	enddo
	4792	endif
	4793	endif
	4794	if(present(rqicuten))then
	4795	if(f_qi) then
	4796	do k=kts,kte
	4797	rqicuten(i,k,j)=rqicuten(i,k,j)*qmemf
	4798	enddo
	4799	endif
	4800	endif
	4801	RAINCV(I,J)=RAINCV(I,J)*qmemf
	4802	PRATEC(I,J)=PRATEC(I,J)*qmemf
	4803	!
	4804	! check heating rates
	4805
	4806	!
	4807	thresh=200.
	4808	qmemf=1.
	4809	qmem1=0.
	4810	do k=kts,kte
	4811	qmem1=abs(rthcuten(i,k,j))*86400.
	4812
	4813	if(qmem1.gt.thresh)then
	4814	qmem2=thresh/qmem1
	4815	qmemf=min(qmemf,qmem2)
	4816	qmemf=max(0.,qmemf)
	4817	endif
	4818
	4819	enddo
	4820	RAINCV(I,J)=RAINCV(I,J)*qmemf
	4821	PRATEC(I,J)=PRATEC(I,J)*qmemf
	4822	do k=kts,kte
	4823	rqvcuten(i,k,j)=rqvcuten(i,k,j)*qmemf
	4824	rthcuten(i,k,j)=rthcuten(i,k,j)*qmemf
	4825	enddo
	4826	if(present(rqccuten))then
	4827	if(f_qc) then
	4828	do k=kts,kte
	4829	rqccuten(i,k,j)=rqccuten(i,k,j)*qmemf
	4830	enddo
	4831	endif
	4832	endif
	4833	if(present(rqicuten))then
	4834	if(f_qi) then
	4835	do k=kts,kte
	4836	rqicuten(i,k,j)=rqicuten(i,k,j)*qmemf
	4837	enddo
	4838	endif
	4839	endif
	4840	if(smoothv.eq.1)then
	4841	!
	4842	! smooth for now
	4843	!
	4844	do k=kts+2,kte-2
	4845	conv_trasht(k)= .25(rthcuten(i,k-1,j)+2.rthcuten(i,k,j)+rthcuten(i,k+1,j))
	4846	conv_trashq(k)= .25(rqvcuten(i,k-1,j)+2.rqvcuten(i,k,j)+rqvcuten(i,k+1,j))
	4847	enddo
	4848	do k=kts+2,kte-2
	4849	rthcuten(i,k,j)=conv_trasht(k)
	4850	rqvcuten(i,k,j)=conv_trashq(k)
	4851	enddo
	4852	endif
	4853	do k=kts,kte
	4854	rthcuten(i,k,j)=rthcuten(i,k,j)/pi_phy(i,k,j)
	4855	enddo
	4856	enddo ! end j
	4857	enddo ! end i
	4858	ENDDO
	4859	!$OMP END PARALLEL DO
	4860
	4861
	4862	END SUBROUTINE CONV_GRELL_SPREAD3D
	4863	!-------------------------------------------------------
	4864	END MODULE module_cu_g3

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

Download in other formats:

Original Format