Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_g3.F @ 3094

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

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

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

Download in other formats:

Original Format