Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_gd.F @ 1

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

Property svn:executable set to *

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

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

Download in other formats:

Original Format