Context Navigation

module_cu_gd.F

Last change on this file was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
Property svn:executable set to ``*
File size: 150.1 KB

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

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

Download in other formats:

Original Format