Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_cu_gd.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
Property svn:executable set to ``*
File size: 150.1 KB

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

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

Download in other formats:

Original Format