Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_g3.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 182.1 KB

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

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

Download in other formats:

Original Format