Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_gd.F @ 3667

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

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

Download in other formats:

Original Format