Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_kfeta.F @ 2295

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 110.9 KB

Line
1	MODULE module_cu_kfeta
2
3	USE module_wrf_error
4
5	!
6	! V3.3: A new trigger function is added based Ma and Tan (2009):
7	! Ma, L.-M. and Z.-M. Tan, 2009: Improving the behavior of
8	! the cumulus parameterization for tropical cyclone prediction:
9	! Convection trigger. Atmospheric Research, 92, 190 - 211.
10	!
11	!--------------------------------------------------------------------
12	! Lookup table variables:
13	INTEGER, PARAMETER :: KFNT=250,KFNP=220
14	REAL, DIMENSION(KFNT,KFNP),PRIVATE, SAVE :: TTAB,QSTAB
15	REAL, DIMENSION(KFNP),PRIVATE, SAVE :: THE0K
16	REAL, DIMENSION(200),PRIVATE, SAVE :: ALU
17	REAL, PRIVATE, SAVE :: RDPR,RDTHK,PLUTOP
18	! Note: KF Lookup table is used by subroutines KF_eta_PARA, TPMIX2,
19	! TPMIX2DD, ENVIRTHT
20	! End of Lookup table variables:
21
22	CONTAINS
23
24	SUBROUTINE KF_eta_CPS( &
25	ids,ide, jds,jde, kds,kde &
26	,ims,ime, jms,jme, kms,kme &
27	,its,ite, jts,jte, kts,kte &
28	,trigger &
29	,DT,KTAU,DX,CUDT,CURR_SECS,ADAPT_STEP_FLAG &
30	,rho,RAINCV,PRATEC,NCA &
31	,U,V,TH,T,W,dz8w,Pcps,pi &
32	,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
33	,EP2,SVP1,SVP2,SVP3,SVPT0 &
34	,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT &
35	,QV &
36	! optionals
37	,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
38	,RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN &
39	,RQICUTEN,RQSCUTEN, RQVFTEN &
40	)
41	!
42	!-------------------------------------------------------------
43	IMPLICIT NONE
44	!-------------------------------------------------------------
45	INTEGER, INTENT(IN ) :: &
46	ids,ide, jds,jde, kds,kde, &
47	ims,ime, jms,jme, kms,kme, &
48	its,ite, jts,jte, kts,kte
49
50	INTEGER, INTENT(IN ) :: trigger
51	INTEGER, INTENT(IN ) :: STEPCU
52	LOGICAL, INTENT(IN ) :: warm_rain
53
54	REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1
55	REAL, INTENT(IN ) :: CP,R,G,EP1,EP2
56	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
57
58	INTEGER, INTENT(IN ) :: KTAU
59
60	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
61	INTENT(IN ) :: &
62	U, &
63	V, &
64	W, &
65	TH, &
66	T, &
67	QV, &
68	dz8w, &
69	Pcps, &
70	rho, &
71	pi
72	!
73	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
74	INTENT(INOUT) :: &
75	W0AVG
76
77	REAL, INTENT(IN ) :: DT, DX
78	REAL, INTENT(IN ) :: CUDT
79	REAL, INTENT(IN ) :: CURR_SECS
80	LOGICAL,INTENT(IN ) :: ADAPT_STEP_FLAG
81	!
82	REAL, DIMENSION( ims:ime , jms:jme ), &
83	INTENT(INOUT) :: RAINCV
84
85	REAL, DIMENSION( ims:ime , jms:jme ), &
86	INTENT(INOUT) :: PRATEC
87
88	REAL, DIMENSION( ims:ime , jms:jme ), &
89	INTENT(INOUT) :: NCA
90
91	REAL, DIMENSION( ims:ime , jms:jme ), &
92	INTENT(OUT) :: CUBOT, &
93	CUTOP
94
95	LOGICAL, DIMENSION( ims:ime , jms:jme ), &
96	INTENT(INOUT) :: CU_ACT_FLAG
97
98	!
99	! Optional arguments
100	!
101
102	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
103	OPTIONAL, &
104	INTENT(INOUT) :: &
105	RTHCUTEN, &
106	RQVCUTEN, &
107	RQCCUTEN, &
108	RQRCUTEN, &
109	RQICUTEN, &
110	RQSCUTEN, &
111	RQVFTEN
112	!
113	! Flags relating to the optional tendency arrays declared above
114	! Models that carry the optional tendencies will provdide the
115	! optional arguments at compile time; these flags all the model
116	! to determine at run-time whether a particular tracer is in
117	! use or not.
118	!
119	LOGICAL, OPTIONAL :: &
120	F_QV &
121	,F_QC &
122	,F_QR &
123	,F_QI &
124	,F_QS
125
126
127	! LOCAL VARS
128
129	LOGICAL :: flag_qr, flag_qi, flag_qs
130
131	REAL, DIMENSION( kts:kte ) :: &
132	U1D, &
133	V1D, &
134	T1D, &
135	DZ1D, &
136	QV1D, &
137	P1D, &
138	RHO1D, &
139	tpart_v1D, &
140	tpart_h1D, &
141	W0AVG1D
142
143	REAL, DIMENSION( kts:kte ):: &
144	DQDT, &
145	DQIDT, &
146	DQCDT, &
147	DQRDT, &
148	DQSDT, &
149	DTDT
150
151	REAL, DIMENSION (its-1:ite+1,kts:kte,jts-1:jte+1) :: aveh_t, aveh_q
152	REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: aveh_qmax, aveh_qmin
153	REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: avev_t, avev_q
154	REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: avev_qmax, avev_qmin
155	REAL, DIMENSION (its:ite,kts:kte,jts:jte) :: coef_v, coef_h, tpart_h, tpart_v
156	INTEGER :: ii,jj,kk
157
158	REAL :: ttop
159	REAL, DIMENSION (kts:kte) :: z0
160
161	REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp
162	integer :: ibegh,iendh,jbegh,jendh
163	integer :: istart,iend,jstart,jend
164	INTEGER :: i,j,k,NTST
165	REAL :: lastdt = -1.0
166	REAL :: W0AVGfctr, W0fctr, W0den
167	LOGICAL :: run_param
168
169	!
170	DXSQ=DX*DX
171
172	!----------------------
173	NTST=STEPCU
174	TST=float(NTST*2)
175	flag_qr = .FALSE.
176	flag_qi = .FALSE.
177	flag_qs = .FALSE.
178	IF ( PRESENT(F_QR) ) flag_qr = F_QR
179	IF ( PRESENT(F_QI) ) flag_qi = F_QI
180	IF ( PRESENT(F_QS) ) flag_qs = F_QS
181	!
182	if (lastdt < 0) then
183	lastdt = dt
184	endif
185
186	if (ADAPT_STEP_FLAG) then
187	W0AVGfctr = 2 * MAX(CUDT*60,dt) - dt
188	W0fctr = dt
189	W0den = 2 * MAX(CUDT*60,dt)
190	else
191	W0AVGfctr = (TST-1.)
192	W0fctr = 1.
193	W0den = TST
194	endif
195
196	DO J = jts,jte
197	DO K=kts,kte
198	DO I= its,ite
199	! SCR1=-5.0E-4Grho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
200	! TV=T(I,K,J)(1.+EP1QV(I,K,J))
201	! RHOE=Pcps(I,K,J)/(R*TV)
202	! W0=-101.9368*SCR1/RHOE
203	W0=0.5*(w(I,K,J)+w(I,K+1,J))
204
205	! Old:
206	!
207	! W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
208	!
209	! New, to support adaptive time step:
210	!
211	W0AVG(I,K,J) = ( W0AVG(I,K,J) * W0AVGfctr + W0 * W0fctr ) / W0den
212	ENDDO
213	ENDDO
214	ENDDO
215	lastdt = dt
216	!
217	!...CHECK FOR CONVECTIVE INITIATION EVERY 5 MINUTES (OR NTST/2)...
218	!
219	! Modified for adaptive time step
220	!
221	if (ADAPT_STEP_FLAG) then
222	if ( (KTAU .eq. 1) .or. (cudt .eq. 0) .or. &
223	( CURR_SECS + dt >= &
224	( int( CURR_SECS / ( cudt * 60 ) ) + 1 ) * cudt * 60 ) ) then
225	run_param = .TRUE.
226	else
227	run_param = .FALSE.
228	endif
229
230	else
231	if (MOD(KTAU,NTST) .EQ. 0 .or. KTAU .eq. 1) then
232	run_param = .TRUE.
233	else
234	run_param = .FALSE.
235	endif
236	endif
237
238	if (run_param) then
239
240	! New trigger function
241	IF (trigger.eq.2) THEN
242	!
243	! calculate 9-point average of moisture advection and temperature using halo (Horizontal)
244	!
245	aveh_t=-999 ! horizontal 9-point ave
246	aveh_q=-999
247	avev_t=0 ! vertical 3-level ave
248	avev_q=0
249	avev_qmax=0
250	avev_qmin=0
251	aveh_qmax=0
252	aveh_qmin=0
253	tpart_h=0
254	tpart_v=0
255	coef_h=0
256	coef_v=0
257	ibegh=max(its-1, ids+1) ! start from 2
258	jbegh=max(jts-1, jds+1)
259	iendh=min(ite+1, ide-2) ! end at ide-2
260	jendh=min(jte+1, jde-2)
261	DO J = jbegh,jendh
262	DO K = kts,kte
263	DO I = ibegh,iendh
264	aveh_t(i,k,j)=(T(i-1,k,j-1)+T(i-1,k,j) +T(i-1,k,j+1)+ &
265	T(i,k,j-1) +T(i,k,j) +T(i,k,j+1)+ &
266	T(i+1,k,j-1) +T(i+1,k,j) +T(i+1,k,j+1))/9.
267	aveh_q(i,k,j)=(rqvften(i-1,k,j-1)+rqvften(i-1,k,j) +rqvften(i-1,k,j+1)+ &
268	rqvften(i,k,j-1) +rqvften(i,k,j) +rqvften(i,k,j+1)+ &
269	rqvften(i+1,k,j-1) +rqvften(i+1,k,j) +rqvften(i+1,k,j+1))/9.
270	ENDDO
271	ENDDO
272	ENDDO
273	! boundary value ( all processors will do the following? Or just those processsors handling sub-area including boundary)
274	DO K = kts,kte
275	DO J = jts-1,jte+1
276	DO I = its-1,ite+1
277
278	if(i.eq.ids) then
279	aveh_t(i,k,j)=aveh_t(i+1,k,j)
280	aveh_q(i,k,j)=aveh_q(i+1,k,j)
281	elseif(i.eq.ide-1) then
282	aveh_t(i,k,j)=aveh_t(i-1,k,j)
283	aveh_q(i,k,j)=aveh_q(i-1,k,j)
284	endif
285
286	if(j.eq.jds) then
287	aveh_t(i,k,j)=aveh_t(i,k,j+1)
288	aveh_q(i,k,j)=aveh_q(i,k,j+1)
289	elseif(j.eq.jde-1) then
290	aveh_t(i,k,j)=aveh_t(i,k,j-1)
291	aveh_q(i,k,j)=aveh_q(i,k,j-1)
292	endif
293
294	if(j.eq.jds.and.i.eq.ids) then
295	aveh_q(i,k,j)=aveh_q(i+1,k,j+1)
296	aveh_t(i,k,j)=aveh_t(i+1,k,j+1)
297	endif
298
299	if(j.eq.jde-1.and.i.eq.ids) then
300	aveh_q(i,k,j)=aveh_q(i+1,k,j-1)
301	aveh_t(i,k,j)=aveh_t(i+1,k,j-1)
302	endif
303
304	if(j.eq.jde-1.and.i.eq.ide-1) then
305	aveh_q(i,k,j)=aveh_q(i-1,k,j-1)
306	aveh_t(i,k,j)=aveh_t(i-1,k,j-1)
307	endif
308
309	if(j.eq.jds.and.i.eq.ide-1) then
310	aveh_q(i,k,j)=aveh_q(i-1,k,j+1)
311	aveh_t(i,k,j)=aveh_t(i-1,k,j+1)
312	endif
313
314	ENDDO
315	ENDDO
316	ENDDO
317	! search for max/min moisture advection in 9-point range, calculate horizontal T-perturbation (tpart_h)
318	istart=max(its, ids+1) ! start from 2
319	jstart=max(jts, jds+1)
320	iend=min(ite, ide-2) ! end at ide-2
321	jend=min(jte, jde-2)
322	DO K = kts,kte
323	DO J = jstart,jend
324	DO I = istart,iend
325	aveh_qmax(i,k,j)=aveh_q(i,k,j)
326	aveh_qmin(i,k,j)=aveh_q(i,k,j)
327	DO ii=-1, 1
328	DO jj=-1,1
329	if(aveh_q(i+II,k,j+JJ).gt.aveh_qmax(i,k,j)) aveh_qmax(i,k,j)=aveh_q(i+II,k,j+JJ)
330	if(aveh_q(i+II,k,j+JJ).lt.aveh_qmin(i,k,j)) aveh_qmin(i,k,j)=aveh_q(i+II,k,j+JJ)
331	ENDDO
332	ENDDO
333	if(aveh_qmax(i,k,j).gt.aveh_qmin(i,k,j))then
334	coef_h(i,k,j)=(aveh_q(i,k,j)-aveh_qmin(i,k,j))/(aveh_qmax(i,k,j)-aveh_qmin(i,k,j))
335	else
336	coef_h(i,k,j)=0.
337	endif
338	coef_h(i,k,j)=amin1(coef_h(i,k,j),1.0)
339	coef_h(i,k,j)=amax1(coef_h(i,k,j),0.0)
340	tpart_h(i,k,j)=coef_h(i,k,j)*(T(i,k,j)-aveh_t(i,k,j))
341	ENDDO
342	ENDDO
343	ENDDO
344	89 continue
345	! vertical 3-layer calculation
346	DO J = jts, jte
347	DO I = its, ite
348	z0(1) = 0.5 * dz8w(i,1,j)
349	DO K = 2, kte
350	Z0(K) = Z0(K-1) + .5 * (DZ8W(i,K,j) + DZ8W(i,K-1,j))
351	ENDDO
352	DO K = kts+1,kte-1
353	ttop = t(i,k,j) + ((t(i,k,j) - t(i,k+1,j)) / (z0(k) - z0(k+1))) * (z0(k)-z0(k-1))
354	avev_t(i,k,j)=(T(i,k-1,j) + T(i,k,j) + ttop)/3.
355	! avev_t(i,k,j)=(T(i,k-1,j)+T(i,k,j) + T(i,k+1,j))/3.
356	avev_q(i,k,j)=(rqvften(i,k-1,j)+rqvften(i,k,j) + rqvften(i,k+1,j))/3.
357	ENDDO
358	avev_t(i,kts,j)=avev_t(i,kts+1,j) ! lowest level value, is it the same as avev_t(i,kds,j)=avev_t(i,kds+1,j)?
359	avev_q(i,kts,j)=avev_q(i,kts+1,j)
360	avev_t(i,kte,j)=avev_t(i,kte-1,j) ! highest level value
361	avev_q(i,kte,j)=avev_q(i,kte-1,j)
362	ENDDO
363	ENDDO
364	! max /min value
365	DO J = jts, jte
366	DO I = its, ite
367	DO K = kts+1,kte-1
368	avev_qmax(i,k,j)=avev_q(i,k,j)
369	avev_qmin(i,k,j)=avev_q(i,k,j)
370	DO kk=-1,1
371	if(avev_q(i,k+kk,j).gt.avev_qmax(i,k,j)) avev_qmax(i,k,j)=avev_q(i,k+kk,j)
372	if(avev_q(i,k+kk,j).lt.avev_qmin(i,k,j)) avev_qmin(i,k,j)=avev_q(i,k+kk,j)
373	ENDDO
374	if(avev_qmax(i,k,j).gt.avev_qmin(i,k,j)) then
375	coef_v(i,k,j)=(avev_q(i,k,j)-avev_qmin(i,k,j))/(avev_qmax(i,k,j)-avev_qmin(i,k,j))
376	else
377	coef_v(i,k,j)=0
378	endif
379	tpart_v(i,k,j)=coef_v(i,k,j)*(T(i,k,j)-avev_t(i,k,j))
380	ENDDO
381	tpart_v(i,kts,j)= tpart_v(i,kts+1,j) ! lowest level
382	tpart_v(i,kte,j)= tpart_v(i,kte-1,j) ! highest level
383	ENDDO
384	ENDDO
385	ENDIF ! endif (trigger.eq.2)
386	!
387	DO J = jts,jte
388	DO I= its,ite
389	CU_ACT_FLAG(i,j) = .true.
390	ENDDO
391	ENDDO
392
393	DO J = jts,jte
394	DO I=its,ite
395
396
397	IF ( NCA(I,J) .ge. 0.5*DT ) then
398	CU_ACT_FLAG(i,j) = .false.
399	ELSE
400
401	DO k=kts,kte
402	DQDT(k)=0.
403	DQIDT(k)=0.
404	DQCDT(k)=0.
405	DQRDT(k)=0.
406	DQSDT(k)=0.
407	DTDT(k)=0.
408	ENDDO
409	RAINCV(I,J)=0.
410	CUTOP(I,J)=KTS
411	CUBOT(I,J)=KTE+1
412	PRATEC(I,J)=0.
413	!
414	! assign vars from 3D to 1D
415
416	DO K=kts,kte
417	U1D(K) =U(I,K,J)
418	V1D(K) =V(I,K,J)
419	T1D(K) =T(I,K,J)
420	RHO1D(K) =rho(I,K,J)
421	QV1D(K)=QV(I,K,J)
422	P1D(K) =Pcps(I,K,J)
423	W0AVG1D(K) =W0AVG(I,K,J)
424	DZ1D(k)=dz8w(I,K,J)
425	IF (trigger.eq.2) THEN
426	tpart_h1D(K) =tpart_h(I,K,J)
427	tpart_v1D(K) =tpart_v(I,K,J)
428	ELSE
429	tpart_h1D(K) = 0.
430	tpart_v1D(K) = 0.
431	ENDIF
432	ENDDO
433	CALL KF_eta_PARA(I, J, &
434	U1D,V1D,T1D,QV1D,P1D,DZ1D,W0AVG1D, &
435	tpart_h1D,tpart_v1D, &
436	trigger, &
437	DT,DX,DXSQ,RHO1D, &
438	XLV0,XLV1,XLS0,XLS1,CP,R,G, &
439	EP2,SVP1,SVP2,SVP3,SVPT0, &
440	DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
441	RAINCV,PRATEC,NCA, &
442	flag_QI,flag_QS,warm_rain, &
443	CUTOP,CUBOT,CUDT, &
444	ids,ide, jds,jde, kds,kde, &
445	ims,ime, jms,jme, kms,kme, &
446	its,ite, jts,jte, kts,kte)
447	IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
448	DO K=kts,kte
449	RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
450	RQVCUTEN(I,K,J)=DQDT(K)
451	ENDDO
452	ENDIF
453
454	IF(PRESENT(rqrcuten).AND.PRESENT(rqccuten)) THEN
455	IF( F_QR )THEN
456	DO K=kts,kte
457	RQRCUTEN(I,K,J)=DQRDT(K)
458	RQCCUTEN(I,K,J)=DQCDT(K)
459	ENDDO
460	ELSE
461	! This is the case for Eta microphysics without 3d rain field
462	DO K=kts,kte
463	RQRCUTEN(I,K,J)=0.
464	RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
465	ENDDO
466	ENDIF
467	ENDIF
468
469	!...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
470
471
472	IF(PRESENT( rqicuten )) THEN
473	IF ( F_QI ) THEN
474	DO K=kts,kte
475	RQICUTEN(I,K,J)=DQIDT(K)
476	ENDDO
477	ENDIF
478	ENDIF
479
480	IF(PRESENT( rqscuten )) THEN
481	IF ( F_QS ) THEN
482	DO K=kts,kte
483	RQSCUTEN(I,K,J)=DQSDT(K)
484	ENDDO
485	ENDIF
486	ENDIF
487	!
488	ENDIF
489	ENDDO ! i-loop
490	ENDDO ! j-loop
491	ENDIF ! run_param
492	!
493	END SUBROUTINE KF_eta_CPS
494	! ****************************************************************************
495	!-----------------------------------------------------------
496	SUBROUTINE KF_eta_PARA (I, J, &
497	U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
498	TPART_H0,TPART_V0, &
499	trigger, &
500	DT,DX,DXSQ,rhoe, &
501	XLV0,XLV1,XLS0,XLS1,CP,R,G, &
502	EP2,SVP1,SVP2,SVP3,SVPT0, &
503	DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
504	RAINCV,PRATEC,NCA, &
505	F_QI,F_QS,warm_rain, &
506	CUTOP,CUBOT,CUDT, &
507	ids,ide, jds,jde, kds,kde, &
508	ims,ime, jms,jme, kms,kme, &
509	its,ite, jts,jte, kts,kte)
510	!-----------------------------------------------------------
511	!***** The KF scheme that is currently used in experimental runs of EMCs
512	!***** Eta model....jsk 8/00
513	!
514	IMPLICIT NONE
515	!-----------------------------------------------------------
516	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
517	ims,ime, jms,jme, kms,kme, &
518	its,ite, jts,jte, kts,kte, &
519	I,J
520	! ,P_QI,P_QS,P_FIRST_SCALAR
521	INTEGER, INTENT(IN ) :: trigger
522
523	LOGICAL, INTENT(IN ) :: F_QI, F_QS
524
525	LOGICAL, INTENT(IN ) :: warm_rain
526	!
527	REAL, DIMENSION( kts:kte ), &
528	INTENT(IN ) :: U0, &
529	V0, &
530	TPART_H0, &
531	TPART_V0, &
532	T0, &
533	QV0, &
534	P0, &
535	rhoe, &
536	DZQ, &
537	W0AVG1D
538	!
539	REAL, INTENT(IN ) :: DT,DX,DXSQ
540	!
541
542	REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
543	REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
544
545	!
546	REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
547	DQDT, &
548	DQIDT, &
549	DQCDT, &
550	DQRDT, &
551	DQSDT, &
552	DTDT
553
554	REAL, DIMENSION( ims:ime , jms:jme ), &
555	INTENT(INOUT) :: NCA
556
557	REAL, DIMENSION( ims:ime , jms:jme ), &
558	INTENT(INOUT) :: RAINCV
559
560	REAL, DIMENSION( ims:ime , jms:jme ), &
561	INTENT(INOUT) :: PRATEC
562
563	REAL, DIMENSION( ims:ime , jms:jme ), &
564	INTENT(OUT) :: CUBOT, &
565	CUTOP
566	REAL, INTENT(IN ) :: CUDT
567	!
568	!...DEFINE LOCAL VARIABLES...
569	!
570	REAL, DIMENSION( kts:kte ) :: &
571	Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
572	QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
573	UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
574	UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
575	THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
576	QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
577	DETLQ2,DETIC2,RATIO,RATIO2
578
579
580	REAL, DIMENSION( kts:kte ) :: &
581	DOMGDP,EXN,TVQU,DP,RH,EQFRC,WSPD, &
582	QDT,FXM,THTAG,THPA,THFXOUT, &
583	THFXIN,QPA,QFXOUT,QFXIN,QLPA,QLFXIN, &
584	QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
585	QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
586	QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
587
588
589	REAL, DIMENSION( kts:kte+1 ) :: OMG
590	REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
591	REAL, DIMENSION( kts:kte ) :: &
592	CLDHGT,QSD,DILFRC,DDILFRC,TKE,TGU,QGU,THTEEG
593
594	! LOCAL VARS
595
596	REAL :: P00,T00,RLF,RHIC,RHBC,PIE, &
597	TTFRZ,TBFRZ,C5,RATE
598	REAL :: GDRY,ROCP,ALIQ,BLIQ, &
599	CLIQ,DLIQ
600	REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
601	ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
602	CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
603	ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
604	TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
605	UPNEW,ABE,WKLCL,TTEMP,FRC1, &
606	QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
607	DZZ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
608	THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
609	UD1,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
610	THTA,VCONV,TIMEC,SHSIGN,VWS,PEF, &
611	CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
612	DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
613	DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
614	UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
615	DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
616	AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
617	DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
618	TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
619	UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
620	RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
621	DDFRC,TDC,DEFRC,RHBAR,DMFFRC,DPMIN,DILBE
622	REAL :: ASTRT,TP,VALUE,AINTRP,TKEMAX,QFRZ,&
623	QSS,PPTMLT,DTMELT,RHH,EVAC,BINC
624	!
625	INTEGER :: INDLU,NU,NUCHM,NNN,KLFS
626	REAL :: CHMIN,PM15,CHMAX,DTRH,RAD,DPPP
627	REAL :: TVDIFF,DTTOT,ABSOMG,ABSOMGTC,FRDP
628
629	INTEGER :: KX,K,KL
630	!
631	INTEGER :: NCHECK
632	INTEGER, DIMENSION (kts:kte) :: KCHECK
633
634	INTEGER :: ISTOP,ML,L5,KMIX,LOW, &
635	LC,MXLAYR,LLFC,NLAYRS,NK, &
636	KPBL,KLCL,LCL,LET,IFLAG, &
637	NK1,LTOP,NJ,LTOP1, &
638	LTOPM1,LVF,KSTART,KMIN,LFS, &
639	ND,NIC,LDB,LDT,ND1,NDK, &
640	NM,LMAX,NCOUNT,NOITR, &
641	NSTEP,NTC,NCHM,ISHALL,NSHALL
642	LOGICAL :: IPRNT
643	REAL :: u00,qslcl,rhlcl,dqssdt !jfb
644	CHARACTER*1024 message
645	!
646	DATA P00,T00/1.E5,273.16/
647	DATA RLF/3.339E5/
648	DATA RHIC,RHBC/1.,0.90/
649	DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
650	DATA RATE/0.03/ ! wrf default
651	! DATA RATE/0.01/ ! value used in NRCM
652	! DATA RATE/0.001/ ! effectively turn off autoconversion
653	!-----------------------------------------------------------
654	IPRNT=.FALSE.
655	GDRY=-G/CP
656	ROCP=R/CP
657	NSHALL = 0
658	KL=kte
659	KX=kte
660	!
661	! ALIQ = 613.3
662	! BLIQ = 17.502
663	! CLIQ = 4780.8
664	! DLIQ = 32.19
665	ALIQ = SVP1*1000.
666	BLIQ = SVP2
667	CLIQ = SVP2*SVPT0
668	DLIQ = SVP3
669	!
670	!
671	!****************************************************************************
672	! ! PPT FB MODS
673	!...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER ! PPT FB MODS
674	!...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL) ! PPT FB MODS
675	!...FIELD. "FBFRC" IS THE FRACTION OF AVAILABLE ! PPT FB MODS
676	!...PRECIPITATION TO BE FED BACK (0.0 - 1.0)... ! PPT FB MODS
677	FBFRC=0.0 ! PPT FB MODS
678	!...mods to allow shallow convection...
679	NCHM = 0
680	ISHALL = 0
681	DPMIN = 5.E3
682	!...
683	P300=P0(1)-30000.
684	!
685	!...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
686	!...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
687	!...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
688	!
689	!...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
690	!...FROM BOTTOM-UP IN THE KF SCHEME...
691	!
692	ML=0
693	!SUE tmprpsb=1./PSB(I,J)
694	!SUE CELL=PTOP*tmprpsb
695	!
696	DO K=1,KX
697	!
698	! Saturation vapor pressure (ES) is calculated following Buck (1981)
699	!...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
700	!
701	ES=ALIQEXP((BLIQT0(K)-CLIQ)/(T0(K)-DLIQ))
702	QES(K)=0.622*ES/(P0(K)-ES)
703	Q0(K)=AMIN1(QES(K),QV0(K))
704	Q0(K)=AMAX1(0.000001,Q0(K))
705	QL0(K)=0.
706	QI0(K)=0.
707	QR0(K)=0.
708	QS0(K)=0.
709	RH(K) = Q0(K)/QES(K)
710	DILFRC(K) = 1.
711	TV0(K)=T0(K)(1.+0.608Q0(K))
712	! RHOE(K)=P0(K)/(R*TV0(K))
713	! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
714	DP(K)=rhoe(k)gDZQ(k)
715	! IF Turbulent Kinetic Energy (TKE) is available from turbulent mixing scheme
716	! use it for shallow convection...For now, assume it is not available....
717	! TKE(K) = Q2(I,J,NK)
718	TKE(K) = 0.
719	CLDHGT(K) = 0.
720	! IF(P0(K).GE.500E2)L5=K
721	IF(P0(K).GE.0.5*P0(1))L5=K
722	IF(P0(K).GE.P300)LLFC=K
723	ENDDO
724	!
725	!...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
726	Z0(1)=.5*DZQ(1)
727	!cdir novector
728	DO K=2,KL
729	Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
730	DZA(K-1)=Z0(K)-Z0(K-1)
731	ENDDO
732	DZA(KL)=0.
733	!
734	!
735	! To save time, specify a pressure interval to move up in sequential
736	! check of different ~50 mb deep groups of adjacent model layers in
737	! the process of identifying updraft source layer (USL). Note that
738	! this search is terminated as soon as a buoyant parcel is found and
739	! this parcel can produce a cloud greater than specifed minimum depth
740	! (CHMIN)...For now, set interval at 15 mb...
741	!
742	NCHECK = 1
743	KCHECK(NCHECK)=1
744	PM15 = P0(1)-15.E2
745	DO K=2,LLFC
746	IF(P0(K).LT.PM15)THEN
747	NCHECK = NCHECK+1
748	KCHECK(NCHECK) = K
749	PM15 = PM15-15.E2
750	ENDIF
751	ENDDO
752	!
753	NU=0
754	NUCHM=0
755	usl: DO
756	NU = NU+1
757	IF(NU.GT.NCHECK)THEN
758	IF(ISHALL.EQ.1)THEN
759	CHMAX = 0.
760	NCHM = 0
761	DO NK = 1,NCHECK
762	NNN=KCHECK(NK)
763	IF(CLDHGT(NNN).GT.CHMAX)THEN
764	NCHM = NNN
765	NUCHM = NK
766	CHMAX = CLDHGT(NNN)
767	ENDIF
768	ENDDO
769	NU = NUCHM-1
770	FBFRC=1.
771	CYCLE usl
772	ELSE
773	RETURN
774	ENDIF
775	ENDIF
776	KMIX = KCHECK(NU)
777	LOW=KMIX
778	!...
779	LC = LOW
780	!
781	!...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
782	!...UNSTABLE AIR AT LEAST 50 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
783	!...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
784	!...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 50 mb..
785	!
786	NLAYRS=0
787	DPTHMX=0.
788	NK=LC-1
789	IF ( NK+1 .LT. KTS ) THEN
790	WRITE(message,*)'WOULD GO OFF BOTTOM: KF_ETA_PARA I,J,NK',I,J,NK
791	CALL wrf_message (TRIM(message))
792	ELSE
793	DO
794	NK=NK+1
795	IF ( NK .GT. KTE ) THEN
796	WRITE(message,*)'WOULD GO OFF TOP: KF_ETA_PARA I,J,DPTHMX,DPMIN',I,J,DPTHMX,DPMIN
797	CALL wrf_message (TRIM(message))
798	EXIT
799	ENDIF
800	DPTHMX=DPTHMX+DP(NK)
801	NLAYRS=NLAYRS+1
802	IF(DPTHMX.GT.DPMIN)THEN
803	EXIT
804	ENDIF
805	END DO
806	ENDIF
807	IF(DPTHMX.LT.DPMIN)THEN
808	RETURN
809	ENDIF
810	KPBL=LC+NLAYRS-1
811	!
812	!...********************************************************
813	!...for computational simplicity without much loss in accuracy,
814	!...mix temperature instead of theta for evaluating convective
815	!...initiation (triggering) potential...
816	! THMIX=0.
817	TMIX=0.
818	QMIX=0.
819	ZMIX=0.
820	PMIX=0.
821	!
822	!...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
823	!...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
824	!...LAYERS...
825	!
826	!cdir novector
827	DO NK=LC,KPBL
828	TMIX=TMIX+DP(NK)*T0(NK)
829	QMIX=QMIX+DP(NK)*Q0(NK)
830	ZMIX=ZMIX+DP(NK)*Z0(NK)
831	PMIX=PMIX+DP(NK)*P0(NK)
832	ENDDO
833	! THMIX=THMIX/DPTHMX
834	TMIX=TMIX/DPTHMX
835	QMIX=QMIX/DPTHMX
836	ZMIX=ZMIX/DPTHMX
837	PMIX=PMIX/DPTHMX
838	EMIX=QMIX*PMIX/(0.622+QMIX)
839	!
840	!...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL...
841	!
842	! TLOG=ALOG(EMIX/ALIQ)
843	! ...calculate dewpoint using lookup table...
844	!
845	astrt=1.e-3
846	ainc=0.075
847	a1=emix/aliq
848	tp=(a1-astrt)/ainc
849	indlu=int(tp)+1
850	value=(indlu-1)*ainc+astrt
851	aintrp=(a1-value)/ainc
852	tlog=aintrpalu(indlu+1)+(1-aintrp)alu(indlu)
853	TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
854	TLCL=TDPT-(.212+1.571E-3(TDPT-T00)-4.36E-4(TMIX-T00))*(TMIX-TDPT)
855	TLCL=AMIN1(TLCL,TMIX)
856	TVLCL=TLCL(1.+0.608QMIX)
857	ZLCL = ZMIX+(TLCL-TMIX)/GDRY
858	NK = LC-1
859	DO
860	NK = NK+1
861	KLCL=NK
862	IF(ZLCL.LE.Z0(NK) .or. NK.GT.KL)THEN
863	EXIT
864	ENDIF
865	ENDDO
866	IF(NK.GT.KL)THEN
867	RETURN
868	ENDIF
869	K=KLCL-1
870	! calculate DLP using Z instead of log(P)
871	DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
872	!
873	!...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
874	!
875	TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
876	QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
877	TVEN=TENV(1.+0.608QENV)
878	!
879	!...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
880	!...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0 IS AN
881	!...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
882	!...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
883	!...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
884	!...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
885	!...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
886	!...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
887	!...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
888	!
889	IF(ZLCL.LT.2.E3)THEN ! Kain (2004) Eq. 2
890	WKLCL=0.02*ZLCL/2.E3
891	ELSE
892	WKLCL=0.02 ! units of m/s
893	ENDIF
894	WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))DLP)DX/25.E3-WKLCL
895	IF(WKL.LT.0.0001)THEN
896	DTLCL=0.
897	ELSE
898	DTLCL=4.64WKL*0.33 ! Kain (2004) Eq. 1
899	ENDIF
900
901	! New trigger function
902	IF(trigger.eq.2) then
903	DTLCL=amax1(TPART_H0(KLCL)+TPART_V0(KLCL),0.0)
904	ENDIF
905	! end for trigger function
906	!
907	dtrh = 0.
908	if (trigger .eq. 3) then
909	!...for ETA model, give parcel an extra temperature perturbation based
910	!...the threshold RH for condensation (U00)...
911	! as described in Narita and Ohmori (2007, 12th Mesoscale Conf.
912	!
913	!...for now, just assume U00=0.75...
914	!...!!!!!! for MM5, SET DTRH = 0. !!!!!!!!
915	U00 = 0.75
916	IF(U00.lt.1.)THEN
917	QSLCL=QES(K)+(QES(KLCL)-QES(K))*DLP
918	RHLCL = QENV/QSLCL
919	DQSSDT = QMIX(CLIQ-BLIQDLIQ)/((TLCL-DLIQ)*(TLCL-DLIQ))
920	IF(RHLCL.ge.0.75 .and. RHLCL.le.0.95)then
921	DTRH = 0.25(RHLCL-0.75)QMIX/DQSSDT
922	ELSEIF(RHLCL.GT.0.95)THEN
923	DTRH = (1./RHLCL-1.)*QMIX/DQSSDT
924	ELSE
925	DTRH = 0.
926	ENDIF
927	ENDIF
928	endif ! trigger 3
929	! IF(ISHALL.EQ.1)IPRNT=.TRUE.
930	! IPRNT=.TRUE.
931	! IF(TLCL+DTLCL.GT.TENV)GOTO 45
932
933	trigger2: IF(TLCL+DTLCL+DTRH.LT.TENV)THEN
934	!
935	! Parcel not buoyant, CYCLE back to start of trigger and evaluate next potential USL...
936	!
937	CYCLE usl
938	!
939	ELSE ! Parcel is buoyant, determine updraft
940	!
941	!...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
942	!...EQUIVALENT POTENTIAL TEMPERATURE
943	!...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
944	!
945	CALL ENVIRTHT(PMIX,TMIX,QMIX,THETEU(K),ALIQ,BLIQ,CLIQ,DLIQ)
946	!
947	!...modify calculation of initial parcel vertical velocity...jsk 11/26/97
948	!
949	DTTOT = DTLCL+DTRH
950	IF(DTTOT.GT.1.E-4)THEN
951	GDT=2.GDTTOT*500./TVEN ! Kain (2004) Eq. 3 (sort of)
952	WLCL=1.+0.5*SQRT(GDT)
953	WLCL = AMIN1(WLCL,3.)
954	ELSE
955	WLCL=1.
956	ENDIF
957	PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
958	WTW=WLCL*WLCL
959	!
960	TVLCL=TLCL(1.+0.608QMIX)
961	RHOLCL=PLCL/(R*TVLCL)
962	!
963	LCL=KLCL
964	LET=LCL
965	! make RAD a function of background vertical velocity... (Kain (2004) Eq. 6)
966	IF(WKL.LT.0.)THEN
967	RAD = 1000.
968	ELSEIF(WKL.GT.0.1)THEN
969	RAD = 2000.
970	ELSE
971	RAD = 1000.+1000*WKL/0.1
972	ENDIF
973	!
974	!*******************************************************************
975	! *
976	! COMPUTE UPDRAFT PROPERTIES *
977	! *
978	!*******************************************************************
979	!
980	!
981	!...
982	!...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
983	!
984	WU(K)=WLCL
985	AU0=0.01*DXSQ
986	UMF(K)=RHOLCL*AU0
987	VMFLCL=UMF(K)
988	UPOLD=VMFLCL
989	UPNEW=UPOLD
990	!
991	!...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
992	!...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE
993	!...BUOYANT ENERGY, TRPPT IS THE TOTAL RATE OF PRECIPITATION
994	!...PRODUCTION...
995	!
996	RATIO2(K)=0.
997	UER(K)=0.
998	ABE=0.
999	TRPPT=0.
1000	TU(K)=TLCL
1001	TVU(K)=TVLCL
1002	QU(K)=QMIX
1003	EQFRC(K)=1.
1004	QLIQ(K)=0.
1005	QICE(K)=0.
1006	QLQOUT(K)=0.
1007	QICOUT(K)=0.
1008	DETLQ(K)=0.
1009	DETIC(K)=0.
1010	PPTLIQ(K)=0.
1011	PPTICE(K)=0.
1012	IFLAG=0
1013	!
1014	!...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
1015	!...PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
1016	!...FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
1017	!...INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
1018	!...PREVIOUS MODEL LEVEL...
1019	!
1020	TTEMP=TTFRZ
1021	!
1022	!...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
1023	!...MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
1024	!...MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
1025	!
1026	! **1 variables indicate the bottom of a model layer
1027	! **2 variables indicate the top of a model layer
1028	!
1029	EE1=1.
1030	UD1=0.
1031	REI = 0.
1032	DILBE = 0.
1033	updraft: DO NK=K,KL-1
1034	NK1=NK+1
1035	RATIO2(NK1)=RATIO2(NK)
1036	FRC1=0.
1037	TU(NK1)=T0(NK1)
1038	THETEU(NK1)=THETEU(NK)
1039	QU(NK1)=QU(NK)
1040	QLIQ(NK1)=QLIQ(NK)
1041	QICE(NK1)=QICE(NK)
1042	call tpmix2(p0(nk1),theteu(nk1),tu(nk1),qu(nk1),qliq(nk1), &
1043	qice(nk1),qnewlq,qnewic,XLV1,XLV0)
1044	!
1045	!
1046	!...CHECK TO SEE IF UPDRAFT TEMP IS ABOVE THE TEMPERATURE AT WHICH
1047	!...GLACIATION IS ASSUMED TO INITIATE; IF IT IS, CALCULATE THE
1048	!...FRACTION OF REMAINING LIQUID WATER TO FREEZE...TTFRZ IS THE
1049	!...TEMP AT WHICH FREEZING BEGINS, TBFRZ THE TEMP BELOW WHICH ALL
1050	!...LIQUID WATER IS FROZEN AT EACH LEVEL...
1051	!
1052	IF(TU(NK1).LE.TTFRZ)THEN
1053	IF(TU(NK1).GT.TBFRZ)THEN
1054	IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
1055	FRC1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
1056	ELSE
1057	FRC1=1.
1058	IFLAG=1
1059	ENDIF
1060	TTEMP=TU(NK1)
1061	!
1062	! DETERMINE THE EFFECTS OF LIQUID WATER FREEZING WHEN TEMPERATURE
1063	!...IS BELOW TTFRZ...
1064	!
1065	QFRZ = (QLIQ(NK1)+QNEWLQ)*FRC1
1066	QNEWIC=QNEWIC+QNEWLQ*FRC1
1067	QNEWLQ=QNEWLQ-QNEWLQ*FRC1
1068	QICE(NK1) = QICE(NK1)+QLIQ(NK1)*FRC1
1069	QLIQ(NK1) = QLIQ(NK1)-QLIQ(NK1)*FRC1
1070	CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QFRZ, &
1071	QICE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
1072	ENDIF
1073	TVU(NK1)=TU(NK1)(1.+0.608QU(NK1))
1074	!
1075	! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
1076	!
1077	IF(NK.EQ.K)THEN
1078	BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
1079	BOTERM=2.(Z0(NK1)-ZLCL)G*BE/1.5
1080	DZZ=Z0(NK1)-ZLCL
1081	ELSE
1082	BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
1083	BOTERM=2.DZA(NK)G*BE/1.5
1084	DZZ=DZA(NK)
1085	ENDIF
1086	ENTERM=2.REIWTW/UPOLD
1087
1088	CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM, &
1089	RATE,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1),G)
1090	!
1091	!...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
1092	!...IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
1093	!
1094	IF(WTW.LT.1.E-3)THEN
1095	EXIT
1096	ELSE
1097	WU(NK1)=SQRT(WTW)
1098	ENDIF
1099	!...Calculate value of THETA-E in environment to entrain into updraft...
1100	!
1101	CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
1102	!
1103	!...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
1104	!
1105	REI=VMFLCLDP(NK1)0.03/RAD ! KF (1990) Eq. 1; Kain (2004) Eq. 5
1106	TVQU(NK1)=TU(NK1)(1.+0.608QU(NK1)-QLIQ(NK1)-QICE(NK1))
1107	IF(NK.EQ.K)THEN
1108	DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
1109	ELSE
1110	DILBE=((TVQU(NK)+TVQU(NK1))/(TV0(NK)+TV0(NK1))-1.)*DZZ
1111	ENDIF
1112	IF(DILBE.GT.0.)ABE=ABE+DILBE*G
1113	!
1114	!...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, MINIMAL
1115	!...ENTRAINMENT (0.5*REI) IS IMPOSED...
1116	!
1117	IF(TVQU(NK1).LE.TV0(NK1))THEN ! Entrain/Detrain IF BLOCK
1118	EE2=0.5 ! Kain (2004) Eq. 4
1119	UD2=1.
1120	EQFRC(NK1)=0.
1121	ELSE
1122	LET=NK1
1123	TTMP=TVQU(NK1)
1124	!
1125	!...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL AIR...
1126	!
1127	F1=0.95
1128	F2=1.-F1
1129	THTTMP=F1THETEE(NK1)+F2THETEU(NK1)
1130	QTMP=F1Q0(NK1)+F2QU(NK1)
1131	TMPLIQ=F2*QLIQ(NK1)
1132	TMPICE=F2*QICE(NK1)
1133	call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
1134	qnewlq,qnewic,XLV1,XLV0)
1135	TU95=TTMP(1.+0.608QTMP-TMPLIQ-TMPICE)
1136	IF(TU95.GT.TV0(NK1))THEN
1137	EE2=1.
1138	UD2=0.
1139	EQFRC(NK1)=1.0
1140	ELSE
1141	F1=0.10
1142	F2=1.-F1
1143	THTTMP=F1THETEE(NK1)+F2THETEU(NK1)
1144	QTMP=F1Q0(NK1)+F2QU(NK1)
1145	TMPLIQ=F2*QLIQ(NK1)
1146	TMPICE=F2*QICE(NK1)
1147	call tpmix2(p0(nk1),thttmp,ttmp,qtmp,tmpliq,tmpice, &
1148	qnewlq,qnewic,XLV1,XLV0)
1149	TU10=TTMP(1.+0.608QTMP-TMPLIQ-TMPICE)
1150	TVDIFF = ABS(TU10-TVQU(NK1))
1151	IF(TVDIFF.LT.1.e-3)THEN
1152	EE2=1.
1153	UD2=0.
1154	EQFRC(NK1)=1.0
1155	ELSE
1156	EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
1157	EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
1158	EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
1159	IF(EQFRC(NK1).EQ.1)THEN
1160	EE2=1.
1161	UD2=0.
1162	ELSEIF(EQFRC(NK1).EQ.0.)THEN
1163	EE2=0.
1164	UD2=1.
1165	ELSE
1166	!
1167	!...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
1168	! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
1169	!
1170	CALL PROF5(EQFRC(NK1),EE2,UD2)
1171	ENDIF
1172	ENDIF
1173	ENDIF
1174	ENDIF ! End of Entrain/Detrain IF BLOCK
1175	!
1176	!
1177	!...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE FRACTIONAL
1178	! VALUES IN THE LAYER...
1179	!
1180	EE2 = AMAX1(EE2,0.5)
1181	UD2 = 1.5*UD2
1182	UER(NK1)=0.5REI(EE1+EE2)
1183	UDR(NK1)=0.5REI(UD1+UD2)
1184	!
1185	!...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
1186	! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS...
1187	!
1188	IF(UMF(NK)-UDR(NK1).LT.10.)THEN
1189	!
1190	!...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL UPD MASS
1191	! FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE PREVIOUS MODEL LVL..
1192	! First, correct ABE calculation if needed...
1193	!
1194	IF(DILBE.GT.0.)THEN
1195	ABE=ABE-DILBE*G
1196	ENDIF
1197	LET=NK
1198	! WRITE(98,1015)P0(NK1)/100.
1199	EXIT
1200	ELSE
1201	EE1=EE2
1202	UD1=UD2
1203	UPOLD=UMF(NK)-UDR(NK1)
1204	UPNEW=UPOLD+UER(NK1)
1205	UMF(NK1)=UPNEW
1206	DILFRC(NK1) = UPNEW/UPOLD
1207	!
1208	!...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND
1209	!...ICE IN THE DETRAINING UPDRAFT MASS...
1210	!
1211	DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
1212	DETIC(NK1)=QICE(NK1)*UDR(NK1)
1213	QDT(NK1)=QU(NK1)
1214	QU(NK1)=(UPOLDQU(NK1)+UER(NK1)Q0(NK1))/UPNEW
1215	THETEU(NK1)=(THETEU(NK1)UPOLD+THETEE(NK1)UER(NK1))/UPNEW
1216	QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
1217	QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
1218	!
1219	!...PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF
1220	!...LIQUID PRECIP AT A GIVEN MODEL LVL, PPTICE THE SAME FOR ICE,
1221	!...TRPPT IS THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE
1222	!...CURRENT MODEL LEVEL...
1223	!
1224	PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
1225	PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
1226	!
1227	TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
1228	IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
1229	ENDIF
1230	!
1231	END DO updraft
1232	!
1233	!...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIUM
1234	! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
1235	! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE BETWEEN
1236	! THE LET AND CLOUD TOP...
1237	!
1238	!...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL VELOCITY
1239	! FIRST BECOMES NEGATIVE...
1240	!
1241	LTOP=NK
1242	CLDHGT(LC)=Z0(LTOP)-ZLCL
1243	!
1244	!...Instead of using the same minimum cloud height (for deep convection)
1245	!...everywhere, try specifying minimum cloud depth as a function of TLCL...
1246	!
1247	! Kain (2004) Eq. 7
1248	!
1249	IF(TLCL.GT.293.)THEN
1250	CHMIN = 4.E3
1251	ELSEIF(TLCL.LE.293. .and. TLCL.GE.273)THEN
1252	CHMIN = 2.E3 + 100.*(TLCL-273.)
1253	ELSEIF(TLCL.LT.273.)THEN
1254	CHMIN = 2.E3
1255	ENDIF
1256
1257	!
1258	!...If cloud top height is less than the specified minimum for deep
1259	!...convection, save value to consider this level as source for
1260	!...shallow convection, go back up to check next level...
1261	!
1262	!...Try specifying minimum cloud depth as a function of TLCL...
1263	!
1264	!
1265	!...DO NOT ALLOW ANY CLOUD FROM THIS LAYER IF:
1266	!
1267	!... 1.) if there is no CAPE, or
1268	!... 2.) cloud top is at model level just above LCL, or
1269	!... 3.) cloud top is within updraft source layer, or
1270	!... 4.) cloud-top detrainment layer begins within
1271	!... updraft source layer.
1272	!
1273	IF(LTOP.LE.KLCL .or. LTOP.LE.KPBL .or. LET+1.LE.KPBL)THEN ! No Convection Allowed
1274	CLDHGT(LC)=0.
1275	DO NK=K,LTOP
1276	UMF(NK)=0.
1277	UDR(NK)=0.
1278	UER(NK)=0.
1279	DETLQ(NK)=0.
1280	DETIC(NK)=0.
1281	PPTLIQ(NK)=0.
1282	PPTICE(NK)=0.
1283	ENDDO
1284	!
1285	ELSEIF(CLDHGT(LC).GT.CHMIN .and. ABE.GT.1)THEN ! Deep Convection allowed
1286	ISHALL=0
1287	EXIT usl
1288	ELSE
1289	!
1290	!...TO DISALLOW SHALLOW CONVECTION, COMMENT OUT NEXT LINE !!!!!!!!
1291	ISHALL = 1
1292	IF(NU.EQ.NUCHM)THEN
1293	EXIT usl ! Shallow Convection from this layer
1294	ELSE
1295	! Remember this layer (by virtue of non-zero CLDHGT) as potential shallow-cloud layer
1296	DO NK=K,LTOP
1297	UMF(NK)=0.
1298	UDR(NK)=0.
1299	UER(NK)=0.
1300	DETLQ(NK)=0.
1301	DETIC(NK)=0.
1302	PPTLIQ(NK)=0.
1303	PPTICE(NK)=0.
1304	ENDDO
1305	ENDIF
1306	ENDIF
1307	ENDIF trigger2
1308	END DO usl
1309	IF(ISHALL.EQ.1)THEN
1310	KSTART=MAX0(KPBL,KLCL)
1311	LET=KSTART
1312	endif
1313	!
1314	!...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS FL
1315	! THIS LEVEL...
1316	!
1317	IF(LET.EQ.LTOP)THEN
1318	UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
1319	DETLQ(LTOP)=QLIQ(LTOP)UDR(LTOP)UPNEW/UPOLD
1320	DETIC(LTOP)=QICE(LTOP)UDR(LTOP)UPNEW/UPOLD
1321	UER(LTOP)=0.
1322	UMF(LTOP)=0.
1323	ELSE
1324	!
1325	! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
1326	!
1327	DPTT=0.
1328	DO NJ=LET+1,LTOP
1329	DPTT=DPTT+DP(NJ)
1330	ENDDO
1331	DUMFDP=UMF(LET)/DPTT
1332	!
1333	!...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
1334	! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
1335	!
1336	DO NK=LET+1,LTOP
1337	!
1338	!...entrainment is allowed at every level except for LTOP, so disallow
1339	!...entrainment at LTOP and adjust entrainment rates between LET and LTOP
1340	!...so the the dilution factor due to entyrianment is not changed but
1341	!...the actual entrainment rate will change due due forced total
1342	!...detrainment in this layer...
1343	!
1344	IF(NK.EQ.LTOP)THEN
1345	UDR(NK) = UMF(NK-1)
1346	UER(NK) = 0.
1347	DETLQ(NK) = UDR(NK)QLIQ(NK)DILFRC(NK)
1348	DETIC(NK) = UDR(NK)QICE(NK)DILFRC(NK)
1349	ELSE
1350	UMF(NK)=UMF(NK-1)-DP(NK)*DUMFDP
1351	UER(NK)=UMF(NK)*(1.-1./DILFRC(NK))
1352	UDR(NK)=UMF(NK-1)-UMF(NK)+UER(NK)
1353	DETLQ(NK)=UDR(NK)QLIQ(NK)DILFRC(NK)
1354	DETIC(NK)=UDR(NK)QICE(NK)DILFRC(NK)
1355	ENDIF
1356	IF(NK.GE.LET+2)THEN
1357	TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
1358	PPTLIQ(NK)=UMF(NK-1)*QLQOUT(NK)
1359	PPTICE(NK)=UMF(NK-1)*QICOUT(NK)
1360	TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
1361	ENDIF
1362	ENDDO
1363	ENDIF
1364	!
1365	! Initialize some arrays below cloud base and above cloud top...
1366	!
1367	DO NK=1,LTOP
1368	IF(T0(NK).GT.T00)ML=NK
1369	ENDDO
1370	DO NK=1,K
1371	IF(NK.GE.LC)THEN
1372	IF(NK.EQ.LC)THEN
1373	UMF(NK)=VMFLCL*DP(NK)/DPTHMX
1374	UER(NK)=VMFLCL*DP(NK)/DPTHMX
1375	ELSEIF(NK.LE.KPBL)THEN
1376	UER(NK)=VMFLCL*DP(NK)/DPTHMX
1377	UMF(NK)=UMF(NK-1)+UER(NK)
1378	ELSE
1379	UMF(NK)=VMFLCL
1380	UER(NK)=0.
1381	ENDIF
1382	TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
1383	QU(NK)=QMIX
1384	WU(NK)=WLCL
1385	ELSE
1386	TU(NK)=0.
1387	QU(NK)=0.
1388	UMF(NK)=0.
1389	WU(NK)=0.
1390	UER(NK)=0.
1391	ENDIF
1392	UDR(NK)=0.
1393	QDT(NK)=0.
1394	QLIQ(NK)=0.
1395	QICE(NK)=0.
1396	QLQOUT(NK)=0.
1397	QICOUT(NK)=0.
1398	PPTLIQ(NK)=0.
1399	PPTICE(NK)=0.
1400	DETLQ(NK)=0.
1401	DETIC(NK)=0.
1402	RATIO2(NK)=0.
1403	CALL ENVIRTHT(P0(NK),T0(NK),Q0(NK),THETEE(NK),ALIQ,BLIQ,CLIQ,DLIQ)
1404	EQFRC(NK)=1.0
1405	ENDDO
1406	!
1407	LTOP1=LTOP+1
1408	LTOPM1=LTOP-1
1409	!
1410	!...DEFINE VARIABLES ABOVE CLOUD TOP...
1411	!
1412	DO NK=LTOP1,KX
1413	UMF(NK)=0.
1414	UDR(NK)=0.
1415	UER(NK)=0.
1416	QDT(NK)=0.
1417	QLIQ(NK)=0.
1418	QICE(NK)=0.
1419	QLQOUT(NK)=0.
1420	QICOUT(NK)=0.
1421	DETLQ(NK)=0.
1422	DETIC(NK)=0.
1423	PPTLIQ(NK)=0.
1424	PPTICE(NK)=0.
1425	IF(NK.GT.LTOP1)THEN
1426	TU(NK)=0.
1427	QU(NK)=0.
1428	WU(NK)=0.
1429	ENDIF
1430	THTA0(NK)=0.
1431	THTAU(NK)=0.
1432	EMS(NK)=0.
1433	EMSD(NK)=0.
1434	TG(NK)=T0(NK)
1435	QG(NK)=Q0(NK)
1436	QLG(NK)=0.
1437	QIG(NK)=0.
1438	QRG(NK)=0.
1439	QSG(NK)=0.
1440	OMG(NK)=0.
1441	ENDDO
1442	OMG(KX+1)=0.
1443	DO NK=1,LTOP
1444	EMS(NK)=DP(NK)*DXSQ/G
1445	EMSD(NK)=1./EMS(NK)
1446	!
1447	!...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION SCHEME
1448	!
1449	EXN(NK)=(P00/P0(NK))*(0.2854(1.-0.28*QDT(NK)))
1450	THTAU(NK)=TU(NK)*EXN(NK)
1451	EXN(NK)=(P00/P0(NK))*(0.2854(1.-0.28*Q0(NK)))
1452	THTA0(NK)=T0(NK)*EXN(NK)
1453	DDILFRC(NK) = 1./DILFRC(NK)
1454	OMG(NK)=0.
1455	ENDDO
1456	! IF (XTIME.LT.10.)THEN
1457	! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
1458	! * TMIX-T00,PMIX,QMIX,ABE
1459	! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
1460	! * WLCL,CLDHGT
1461	! ENDIF
1462	!
1463	!...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
1464	!...AND MIDTROPOSPHERE IS USED.
1465	!
1466	WSPD(KLCL)=SQRT(U0(KLCL)U0(KLCL)+V0(KLCL)V0(KLCL))
1467	WSPD(L5)=SQRT(U0(L5)U0(L5)+V0(L5)V0(L5))
1468	WSPD(LTOP)=SQRT(U0(LTOP)U0(LTOP)+V0(LTOP)V0(LTOP))
1469	VCONV=.5*(WSPD(KLCL)+WSPD(L5))
1470	!...for ETA model, DX is a function of location...
1471	! TIMEC=DX(I,J)/VCONV
1472	TIMEC=DX/VCONV
1473	TADVEC=TIMEC
1474	TIMEC=AMAX1(1800.,TIMEC) ! 30 minutes >= TIMEC <= 60 minutes
1475	TIMEC=AMIN1(3600.,TIMEC)
1476	IF(ISHALL.EQ.1)TIMEC=2400. ! shallow convection TIMEC = 40 minutes
1477	NIC=NINT(TIMEC/DT)
1478	TIMEC=FLOAT(NIC)*DT
1479	!
1480	!...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
1481	!
1482	IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
1483	SHSIGN=1.
1484	ELSE
1485	SHSIGN=-1.
1486	ENDIF
1487	VWS=(U0(LTOP)-U0(KLCL))(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL)) &
1488	(V0(LTOP)-V0(KLCL))
1489	VWS=1.E3SHSIGNSQRT(VWS)/(Z0(LTOP)-Z0(LCL))
1490	PEF=1.591+VWS(-.639+VWS(9.53E-2-VWS*4.96E-3))
1491	PEF=AMAX1(PEF,.2)
1492	PEF=AMIN1(PEF,.9)
1493	!
1494	!...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
1495	!
1496	CBH=(ZLCL-Z0(1))*3.281E-3
1497	IF(CBH.LT.3.)THEN
1498	RCBH=.02
1499	ELSE
1500	RCBH=.96729352+CBH(-.70034167+CBH(.162179896+CBH*(- &
1501	1.2569798E-2+CBH(4.2772E-4-CBH5.44E-6))))
1502	ENDIF
1503	IF(CBH.GT.25)RCBH=2.4
1504	PEFCBH=1./(1.+RCBH)
1505	PEFCBH=AMIN1(PEFCBH,.9)
1506	!
1507	!... MEAN PEF. IS USED TO COMPUTE RAINFALL.
1508	!
1509	PEFF=.5*(PEF+PEFCBH)
1510	PEFF2 = PEFF ! JSK MODS
1511	IF(IPRNT)THEN
1512	! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1513	WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1514	CALL wrf_message( message )
1515	! call flush(98)
1516	endif
1517	! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1518	!*****************************************************************
1519	! *
1520	! COMPUTE DOWNDRAFT PROPERTIES *
1521	! *
1522	!*****************************************************************
1523	!
1524	!
1525	TDER=0.
1526	devap:IF(ISHALL.EQ.1)THEN
1527	LFS = 1
1528	ELSE
1529	!
1530	!...start downdraft about 150 mb above cloud base...
1531	!
1532	! KSTART=MAX0(KPBL,KLCL)
1533	! KSTART=KPBL ! Changed 7/23/99
1534	KSTART=KPBL+1 ! Changed 7/23/99
1535	KLFS = LET-1
1536	DO NK = KSTART+1,KL
1537	DPPP = P0(KSTART)-P0(NK)
1538	! IF(DPPP.GT.200.E2)THEN
1539	IF(DPPP.GT.150.E2)THEN
1540	KLFS = NK
1541	EXIT
1542	ENDIF
1543	ENDDO
1544	KLFS = MIN0(KLFS,LET-1)
1545	LFS = KLFS
1546	!
1547	!...if LFS is not at least 50 mb above cloud base (implying that the
1548	!...level of equil temp, LET, is just above cloud base) do not allow a
1549	!...downdraft...
1550	!
1551	IF((P0(KSTART)-P0(LFS)).GT.50.E2)THEN
1552	THETED(LFS) = THETEE(LFS)
1553	QD(LFS) = Q0(LFS)
1554	!
1555	!...call tpmix2dd to find wet-bulb temp, qv...
1556	!
1557	call tpmix2dd(p0(lfs),theted(lfs),tz(lfs),qss,i,j)
1558	THTAD(LFS)=TZ(LFS)(P00/P0(LFS))(0.2854(1.-0.28*QSS))
1559	!
1560	!...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX...
1561	!
1562	TVD(LFS)=TZ(LFS)(1.+0.608QSS)
1563	RDD=P0(LFS)/(R*TVD(LFS))
1564	A1=(1.-PEFF)*AU0
1565	DMF(LFS)=-A1*RDD
1566	DER(LFS)=DMF(LFS)
1567	DDR(LFS)=0.
1568	RHBAR = RH(LFS)*DP(LFS)
1569	DPTT = DP(LFS)
1570	DO ND = LFS-1,KSTART,-1
1571	ND1 = ND+1
1572	DER(ND)=DER(LFS)*EMS(ND)/EMS(LFS)
1573	DDR(ND)=0.
1574	DMF(ND)=DMF(ND1)+DER(ND)
1575	THETED(ND)=(THETED(ND1)DMF(ND1)+THETEE(ND)DER(ND))/DMF(ND)
1576	QD(ND)=(QD(ND1)DMF(ND1)+Q0(ND)DER(ND))/DMF(ND)
1577	DPTT = DPTT+DP(ND)
1578	RHBAR = RHBAR+RH(ND)*DP(ND)
1579	ENDDO
1580	RHBAR = RHBAR/DPTT
1581	DMFFRC = 2.*(1.-RHBAR) ! Kain (2004) eq. 11
1582	DPDD = 0.
1583	!...Calculate melting effect
1584	!... first, compute total frozen precipitation generated...
1585	!
1586	pptmlt = 0.
1587	DO NK = KLCL,LTOP
1588	PPTMLT = PPTMLT+PPTICE(NK)
1589	ENDDO
1590	if(lc.lt.ml)then
1591	!...For now, calculate melting effect as if DMF = -UMF at KLCL, i.e., as
1592	!...if DMFFRC=1. Otherwise, for small DMFFRC, DTMELT gets too large!
1593	!...12/14/98 jsk...
1594	DTMELT = RLFPPTMLT/(CPUMF(KLCL))
1595	else
1596	DTMELT = 0.
1597	endif
1598	LDT = MIN0(LFS-1,KSTART-1)
1599	!
1600	call tpmix2dd(p0(kstart),theted(kstart),tz(kstart),qss,i,j)
1601	!
1602	tz(kstart) = tz(kstart)-dtmelt
1603	ES=ALIQEXP((BLIQTZ(KSTART)-CLIQ)/(TZ(KSTART)-DLIQ))
1604	QSS=0.622*ES/(P0(KSTART)-ES)
1605	THETED(KSTART)=TZ(KSTART)(1.E5/P0(KSTART))(0.2854(1.-0.28QSS)) &
1606	EXP((3374.6525/TZ(KSTART)-2.5403)QSS(1.+0.81*QSS))
1607	!....
1608	LDT = MIN0(LFS-1,KSTART-1)
1609	DO ND = LDT,1,-1
1610	DPDD = DPDD+DP(ND)
1611	THETED(ND) = THETED(KSTART)
1612	QD(ND) = QD(KSTART)
1613	!
1614	!...call tpmix2dd to find wet bulb temp, saturation mixing ratio...
1615	!
1616	call tpmix2dd(p0(nd),theted(nd),tz(nd),qss,i,j)
1617	qsd(nd) = qss
1618	!
1619	!...specify RH decrease of 20%/km in downdraft...
1620	!
1621	RHH = 1.-0.2/1000.*(Z0(KSTART)-Z0(ND))
1622	!
1623	!...adjust downdraft TEMP, Q to specified RH:
1624	!
1625	IF(RHH.LT.1.)THEN
1626	DSSDT=(CLIQ-BLIQDLIQ)/((TZ(ND)-DLIQ)(TZ(ND)-DLIQ))
1627	RL=XLV0-XLV1*TZ(ND)
1628	DTMP=RLQSS(1.-RHH)/(CP+RLRHHQSS*DSSDT)
1629	T1RH=TZ(ND)+DTMP
1630	ES=RHHALIQEXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
1631	QSRH=0.622*ES/(P0(ND)-ES)
1632	!
1633	!...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
1634	!...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
1635	!
1636	IF(QSRH.LT.QD(ND))THEN
1637	QSRH=QD(ND)
1638	T1RH=TZ(ND)+(QSS-QSRH)*RL/CP
1639	ENDIF
1640	TZ(ND)=T1RH
1641	QSS=QSRH
1642	QSD(ND) = QSS
1643	ENDIF
1644	TVD(nd) = tz(nd)(1.+0.608qsd(nd))
1645	IF(TVD(ND).GT.TV0(ND).OR.ND.EQ.1)THEN
1646	LDB=ND
1647	EXIT
1648	ENDIF
1649	ENDDO
1650	IF((P0(LDB)-P0(LFS)) .gt. 50.E2)THEN ! minimum Downdraft depth!
1651	DO ND=LDT,LDB,-1
1652	ND1 = ND+1
1653	DDR(ND) = -DMF(KSTART)*DP(ND)/DPDD
1654	DER(ND) = 0.
1655	DMF(ND) = DMF(ND1)+DDR(ND)
1656	TDER=TDER+(QSD(nd)-QD(ND))*DDR(ND)
1657	QD(ND)=QSD(nd)
1658	THTAD(ND)=TZ(ND)(P00/P0(ND))(0.2854(1.-0.28*QD(ND)))
1659	ENDDO
1660	ENDIF
1661	ENDIF
1662	ENDIF devap
1663	!
1664	!...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
1665	!...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
1666	!
1667	d_mf: IF(TDER.LT.1.)THEN
1668	! WRITE(98,3004)I,J
1669	!3004 FORMAT(' ','No Downdraft!; I=',I3,2X,'J=',I3,'ISHALL =',I2)
1670	PPTFLX=TRPPT
1671	CPR=TRPPT
1672	TDER=0.
1673	CNDTNF=0.
1674	UPDINC=1.
1675	LDB=LFS
1676	DO NDK=1,LTOP
1677	DMF(NDK)=0.
1678	DER(NDK)=0.
1679	DDR(NDK)=0.
1680	THTAD(NDK)=0.
1681	WD(NDK)=0.
1682	TZ(NDK)=0.
1683	QD(NDK)=0.
1684	ENDDO
1685	AINCM2=100.
1686	ELSE
1687	DDINC = -DMFFRC*UMF(KLCL)/DMF(KSTART)
1688	UPDINC=1.
1689	IF(TDER*DDINC.GT.TRPPT)THEN
1690	DDINC = TRPPT/TDER
1691	ENDIF
1692	TDER = TDER*DDINC
1693	DO NK=LDB,LFS
1694	DMF(NK)=DMF(NK)*DDINC
1695	DER(NK)=DER(NK)*DDINC
1696	DDR(NK)=DDR(NK)*DDINC
1697	ENDDO
1698	CPR=TRPPT
1699	PPTFLX = TRPPT-TDER
1700	PEFF=PPTFLX/TRPPT
1701	IF(IPRNT)THEN
1702	! write(98,*)'PRECIP EFFICIENCY =',PEFF
1703	write(message,*)'PRECIP EFFICIENCY =',PEFF
1704	CALL wrf_message(message)
1705	! call flush(98)
1706	ENDIF
1707	!
1708	!
1709	!...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
1710	! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
1711	! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
1712	!
1713	! DO NK=LC,LFS
1714	! UMF(NK)=UMF(NK)*UPDINC
1715	! UDR(NK)=UDR(NK)*UPDINC
1716	! UER(NK)=UER(NK)*UPDINC
1717	! PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
1718	! PPTICE(NK)=PPTICE(NK)*UPDINC
1719	! DETLQ(NK)=DETLQ(NK)*UPDINC
1720	! DETIC(NK)=DETIC(NK)*UPDINC
1721	! ENDDO
1722	!
1723	!...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
1724	!...DOWNDRAFT...
1725	!
1726	IF(LDB.GT.1)THEN
1727	DO NK=1,LDB-1
1728	DMF(NK)=0.
1729	DER(NK)=0.
1730	DDR(NK)=0.
1731	WD(NK)=0.
1732	TZ(NK)=0.
1733	QD(NK)=0.
1734	THTAD(NK)=0.
1735	ENDDO
1736	ENDIF
1737	DO NK=LFS+1,KX
1738	DMF(NK)=0.
1739	DER(NK)=0.
1740	DDR(NK)=0.
1741	WD(NK)=0.
1742	TZ(NK)=0.
1743	QD(NK)=0.
1744	THTAD(NK)=0.
1745	ENDDO
1746	DO NK=LDT+1,LFS-1
1747	TZ(NK)=0.
1748	QD(NK)=0.
1749	THTAD(NK)=0.
1750	ENDDO
1751	ENDIF d_mf
1752	!
1753	!...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE INFLOW
1754	! INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN IS AVAILABLE
1755	! IN THAT LAYER INITIALLY...
1756	!
1757	AINCMX=1000.
1758	LMAX=MAX0(KLCL,LFS)
1759	DO NK=LC,LMAX
1760	IF((UER(NK)-DER(NK)).GT.1.e-3)THEN
1761	AINCM1=EMS(NK)/((UER(NK)-DER(NK))*TIMEC)
1762	AINCMX=AMIN1(AINCMX,AINCM1)
1763	ENDIF
1764	ENDDO
1765	AINC=1.
1766	IF(AINCMX.LT.AINC)AINC=AINCMX
1767	!
1768	!...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRAFT AND DOWNDRAFT...THEY WILL
1769	!...BE ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE STABILIZATION
1770	!...CLOSURE...
1771	!
1772	TDER2=TDER
1773	PPTFL2=PPTFLX
1774	DO NK=1,LTOP
1775	DETLQ2(NK)=DETLQ(NK)
1776	DETIC2(NK)=DETIC(NK)
1777	UDR2(NK)=UDR(NK)
1778	UER2(NK)=UER(NK)
1779	DDR2(NK)=DDR(NK)
1780	DER2(NK)=DER(NK)
1781	UMF2(NK)=UMF(NK)
1782	DMF2(NK)=DMF(NK)
1783	ENDDO
1784	FABE=1.
1785	STAB=0.95
1786	NOITR=0
1787	ISTOP=0
1788	!
1789	IF(ISHALL.EQ.1)THEN ! First for shallow convection
1790	!
1791	! No iteration for shallow convection; if turbulent kinetic energy (TKE) is available
1792	! from a turbulence parameterization, scale cloud-base updraft mass flux as a function
1793	! of TKE, but for now, just specify shallow-cloud mass flux using TKEMAX = 5...
1794	!
1795	!...find the maximum TKE value between LC and KLCL...
1796	! TKEMAX = 0.
1797	TKEMAX = 5.
1798	! DO 173 K = LC,KLCL
1799	! NK = KX-K+1
1800	! TKEMAX = AMAX1(TKEMAX,Q2(I,J,NK))
1801	! 173 CONTINUE
1802	! TKEMAX = AMIN1(TKEMAX,10.)
1803	! TKEMAX = AMAX1(TKEMAX,5.)
1804	!c TKEMAX = 10.
1805	!c...3_24_99...DPMIN was changed for shallow convection so that it is the
1806	!c... the same as for deep convection (5.E3). Since this doubles
1807	!c... (roughly) the value of DPTHMX, add a factor of 0.5 to calcu-
1808	!c... lation of EVAC...
1809	!c EVAC = TKEMAX*0.1
1810	EVAC = 0.5TKEMAX0.1
1811	! AINC = 0.1DPTHMXDXIJDXIJ/(VMFLCLG*TIMEC)
1812	! AINC = EVACDPTHMXDX(I,J)DX(I,J)/(VMFLCLG*TIMEC)
1813	AINC = EVACDPTHMXDXSQ/(VMFLCLGTIMEC)
1814	TDER=TDER2*AINC
1815	PPTFLX=PPTFL2*AINC
1816	DO NK=1,LTOP
1817	UMF(NK)=UMF2(NK)*AINC
1818	DMF(NK)=DMF2(NK)*AINC
1819	DETLQ(NK)=DETLQ2(NK)*AINC
1820	DETIC(NK)=DETIC2(NK)*AINC
1821	UDR(NK)=UDR2(NK)*AINC
1822	UER(NK)=UER2(NK)*AINC
1823	DER(NK)=DER2(NK)*AINC
1824	DDR(NK)=DDR2(NK)*AINC
1825	ENDDO
1826	ENDIF ! Otherwise for deep convection
1827	! use iterative procedure to find mass fluxes...
1828	iter: DO NCOUNT=1,10
1829	!
1830	!*****************************************************************
1831	! *
1832	! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
1833	! *
1834	!*****************************************************************
1835	!
1836	!...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
1837	!...SATISFY MASS CONTINUITY...
1838	!
1839	DTT=TIMEC
1840	DO NK=1,LTOP
1841	DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
1842	IF(NK.GT.1)THEN
1843	OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
1844	ABSOMG = ABS(OMG(NK))
1845	ABSOMGTC = ABSOMG*TIMEC
1846	FRDP = 0.75*DP(NK-1)
1847	IF(ABSOMGTC.GT.FRDP)THEN
1848	DTT1 = FRDP/ABSOMG
1849	DTT=AMIN1(DTT,DTT1)
1850	ENDIF
1851	ENDIF
1852	ENDDO
1853	DO NK=1,LTOP
1854	THPA(NK)=THTA0(NK)
1855	QPA(NK)=Q0(NK)
1856	NSTEP=NINT(TIMEC/DTT+1)
1857	DTIME=TIMEC/FLOAT(NSTEP)
1858	FXM(NK)=OMG(NK)*DXSQ/G
1859	ENDDO
1860	!
1861	!...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
1862	!
1863	DO NTC=1,NSTEP
1864	!
1865	!...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED ON
1866	!...SIGN OF OMEGA...
1867	!
1868	DO NK=1,LTOP
1869	THFXIN(NK)=0.
1870	THFXOUT(NK)=0.
1871	QFXIN(NK)=0.
1872	QFXOUT(NK)=0.
1873	ENDDO
1874	DO NK=2,LTOP
1875	IF(OMG(NK).LE.0.)THEN
1876	THFXIN(NK)=-FXM(NK)*THPA(NK-1)
1877	QFXIN(NK)=-FXM(NK)*QPA(NK-1)
1878	THFXOUT(NK-1)=THFXOUT(NK-1)+THFXIN(NK)
1879	QFXOUT(NK-1)=QFXOUT(NK-1)+QFXIN(NK)
1880	ELSE
1881	THFXOUT(NK)=FXM(NK)*THPA(NK)
1882	QFXOUT(NK)=FXM(NK)*QPA(NK)
1883	THFXIN(NK-1)=THFXIN(NK-1)+THFXOUT(NK)
1884	QFXIN(NK-1)=QFXIN(NK-1)+QFXOUT(NK)
1885	ENDIF
1886	ENDDO
1887	!
1888	!...UPDATE THE THETA AND QV VALUES AT EACH LEVEL...
1889	!
1890	DO NK=1,LTOP
1891	THPA(NK)=THPA(NK)+(THFXIN(NK)+UDR(NK)THTAU(NK)+DDR(NK) &
1892	THTAD(NK)-THFXOUT(NK)-(UER(NK)-DER(NK))THTA0(NK)) &
1893	DTIME*EMSD(NK)
1894	QPA(NK)=QPA(NK)+(QFXIN(NK)+UDR(NK)QDT(NK)+DDR(NK)QD(NK)- &
1895	QFXOUT(NK)-(UER(NK)-DER(NK))Q0(NK))DTIME*EMSD(NK)
1896	ENDDO
1897	ENDDO
1898	DO NK=1,LTOP
1899	THTAG(NK)=THPA(NK)
1900	QG(NK)=QPA(NK)
1901	ENDDO
1902	!
1903	!...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO, BORROW
1904	!...MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO...
1905	!
1906	DO NK=1,LTOP
1907	IF(QG(NK).LT.0.)THEN
1908	IF(NK.EQ.1)THEN ! JSK MODS
1909	! PRINT *,' PROBLEM WITH KF SCHEME: ' ! JSK MODS
1910	! PRINT *,'QG = 0 AT THE SURFACE!!!!!!!' ! JSK MODS
1911	CALL wrf_error_fatal ( 'QG, QG(NK).LT.0') ! JSK MODS
1912	ENDIF ! JSK MODS
1913	NK1=NK+1
1914	IF(NK.EQ.LTOP)THEN
1915	NK1=KLCL
1916	ENDIF
1917	TMA=QG(NK1)*EMS(NK1)
1918	TMB=QG(NK-1)*EMS(NK-1)
1919	TMM=(QG(NK)-1.E-9)*EMS(NK )
1920	BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
1921	ACOEFF=BCOEFF*TMA/TMB
1922	TMB=TMB*(1.-BCOEFF)
1923	TMA=TMA*(1.-ACOEFF)
1924	IF(NK.EQ.LTOP)THEN
1925	QVDIFF=(QG(NK1)-TMAEMSD(NK1))100./QG(NK1)
1926	! IF(ABS(QVDIFF).GT.1.)THEN
1927	! PRINT *,'!!!WARNING!!! CLOUD BASE WATER VAPOR CHANGES BY ', &
1928	! QVDIFF, &
1929	! '% WHEN MOISTURE IS BORROWED TO PREVENT NEGATIVE ', &
1930	! 'VALUES IN KAIN-FRITSCH'
1931	! ENDIF
1932	ENDIF
1933	QG(NK)=1.E-9
1934	QG(NK1)=TMA*EMSD(NK1)
1935	QG(NK-1)=TMB*EMSD(NK-1)
1936	ENDIF
1937	ENDDO
1938	TOPOMG=(UDR(LTOP)-UER(LTOP))DP(LTOP)EMSD(LTOP)
1939	IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
1940	! WRITE(99,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME; &
1941	! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1942	! TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1943	ISTOP=1
1944	IPRNT=.TRUE.
1945	EXIT iter
1946	ENDIF
1947	!
1948	!...CONVERT THETA TO T...
1949	!
1950	DO NK=1,LTOP
1951	EXN(NK)=(P00/P0(NK))*(0.2854(1.-0.28*QG(NK)))
1952	TG(NK)=THTAG(NK)/EXN(NK)
1953	TVG(NK)=TG(NK)(1.+0.608QG(NK))
1954	ENDDO
1955	IF(ISHALL.EQ.1)THEN
1956	EXIT iter
1957	ENDIF
1958	!
1959	!*******************************************************************
1960	! *
1961	! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
1962	! *
1963	!*******************************************************************
1964	!
1965	!...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
1966	!
1967	! THMIX=0.
1968	TMIX=0.
1969	QMIX=0.
1970	!
1971	!...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
1972	!...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
1973	!...LAYERS...
1974	!
1975	DO NK=LC,KPBL
1976	TMIX=TMIX+DP(NK)*TG(NK)
1977	QMIX=QMIX+DP(NK)*QG(NK)
1978	ENDDO
1979	TMIX=TMIX/DPTHMX
1980	QMIX=QMIX/DPTHMX
1981	ES=ALIQEXP((TMIXBLIQ-CLIQ)/(TMIX-DLIQ))
1982	QSS=0.622*ES/(PMIX-ES)
1983	!
1984	!...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
1985	!
1986	IF(QMIX.GT.QSS)THEN
1987	RL=XLV0-XLV1*TMIX
1988	CPM=CP(1.+0.887QMIX)
1989	DSSDT=QSS(CLIQ-BLIQDLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
1990	DQ=(QMIX-QSS)/(1.+RL*DSSDT/CPM)
1991	TMIX=TMIX+RL/CP*DQ
1992	QMIX=QMIX-DQ
1993	TLCL=TMIX
1994	ELSE
1995	QMIX=AMAX1(QMIX,0.)
1996	EMIX=QMIX*PMIX/(0.622+QMIX)
1997	astrt=1.e-3
1998	binc=0.075
1999	a1=emix/aliq
2000	tp=(a1-astrt)/binc
2001	indlu=int(tp)+1
2002	value=(indlu-1)*binc+astrt
2003	aintrp=(a1-value)/binc
2004	tlog=aintrpalu(indlu+1)+(1-aintrp)alu(indlu)
2005	TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2006	TLCL=TDPT-(.212+1.571E-3(TDPT-T00)-4.36E-4(TMIX-T00))*(TMIX-TDPT)
2007	TLCL=AMIN1(TLCL,TMIX)
2008	ENDIF
2009	TVLCL=TLCL(1.+0.608QMIX)
2010	ZLCL = ZMIX+(TLCL-TMIX)/GDRY
2011	DO NK = LC,KL
2012	KLCL=NK
2013	IF(ZLCL.LE.Z0(NK))THEN
2014	EXIT
2015	ENDIF
2016	ENDDO
2017	K=KLCL-1
2018	DLP=(ZLCL-Z0(K))/(Z0(KLCL)-Z0(K))
2019	!
2020	!...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
2021	!
2022	TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
2023	QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
2024	TVEN=TENV(1.+0.608QENV)
2025	PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
2026	THETEU(K)=TMIX(1.E5/PMIX)(0.2854(1.-0.28QMIX)) &
2027	EXP((3374.6525/TLCL-2.5403)QMIX(1.+0.81*QMIX))
2028	!
2029	!...COMPUTE ADJUSTED ABE(ABEG).
2030	!
2031	ABEG=0.
2032	DO NK=K,LTOPM1
2033	NK1=NK+1
2034	THETEU(NK1) = THETEU(NK)
2035	!
2036	call tpmix2dd(p0(nk1),theteu(nk1),tgu(nk1),qgu(nk1),i,j)
2037	!
2038	TVQU(NK1)=TGU(NK1)(1.+0.608QGU(NK1)-QLIQ(NK1)-QICE(NK1))
2039	IF(NK.EQ.K)THEN
2040	DZZ=Z0(KLCL)-ZLCL
2041	DILBE=((TVLCL+TVQU(NK1))/(TVEN+TVG(NK1))-1.)*DZZ
2042	ELSE
2043	DZZ=DZA(NK)
2044	DILBE=((TVQU(NK)+TVQU(NK1))/(TVG(NK)+TVG(NK1))-1.)*DZZ
2045	ENDIF
2046	IF(DILBE.GT.0.)ABEG=ABEG+DILBE*G
2047	!
2048	!...DILUTE BY ENTRAINMENT BY THE RATE AS ORIGINAL UPDRAFT...
2049	!
2050	CALL ENVIRTHT(P0(NK1),TG(NK1),QG(NK1),THTEEG(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
2051	THETEU(NK1)=THETEU(NK1)DDILFRC(NK1)+THTEEG(NK1)(1.-DDILFRC(NK1))
2052	ENDDO
2053	!
2054	!...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
2055	!...THE PERIOD TIMEC...
2056	!
2057	IF(NOITR.EQ.1)THEN
2058	! write(98,*)' '
2059	! write(98,*)'TAU, I, J, =',NTSD,I,J
2060	! WRITE(98,1060)FABE
2061	! GOTO 265
2062	EXIT iter
2063	ENDIF
2064	DABE=AMAX1(ABE-ABEG,0.1*ABE)
2065	FABE=ABEG/ABE
2066	IF(FABE.GT.1. .and. ISHALL.EQ.0)THEN
2067	! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS
2068	! *GRID POINT; NO CONVECTION ALLOWED!'
2069	RETURN
2070	ENDIF
2071	IF(NCOUNT.NE.1)THEN
2072	IF(ABS(AINC-AINCOLD).LT.0.0001)THEN
2073	NOITR=1
2074	AINC=AINCOLD
2075	CYCLE iter
2076	ENDIF
2077	DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
2078	IF(DFDA.GT.0.)THEN
2079	NOITR=1
2080	AINC=AINCOLD
2081	CYCLE iter
2082	ENDIF
2083	ENDIF
2084	AINCOLD=AINC
2085	FABEOLD=FABE
2086	IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
2087	! write(98,*)' '
2088	! write(98,*)'TAU, I, J, =',NTSD,I,J
2089	! WRITE(98,1055)FABE
2090	! GOTO 265
2091	EXIT
2092	ENDIF
2093	IF((FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB) .or. NCOUNT.EQ.10)THEN
2094	EXIT iter
2095	ELSE
2096	IF(NCOUNT.GT.10)THEN
2097	! write(98,*)' '
2098	! write(98,*)'TAU, I, J, =',NTSD,I,J
2099	! WRITE(98,1060)FABE
2100	! GOTO 265
2101	EXIT
2102	ENDIF
2103	!
2104	!...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE CONVECTIVE
2105	!...MASS FLUX BY THE FACTOR AINC:
2106	!
2107	IF(FABE.EQ.0.)THEN
2108	AINC=AINC*0.5
2109	ELSE
2110	IF(DABE.LT.1.e-4)THEN
2111	NOITR=1
2112	AINC=AINCOLD
2113	CYCLE iter
2114	ELSE
2115	AINC=AINCSTABABE/DABE
2116	ENDIF
2117	ENDIF
2118	! AINC=AMIN1(AINCMX,AINC)
2119	AINC=AMIN1(AINCMX,AINC)
2120	!...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION ! JSK MODS
2121	!...WILL BE MINIMAL SO JUST IGNORE IT... ! JSK MODS
2122	IF(AINC.LT.0.05)then
2123	RETURN ! JSK MODS
2124	ENDIF
2125	! AINC=AMAX1(AINC,0.05) ! JSK MODS
2126	TDER=TDER2*AINC
2127	PPTFLX=PPTFL2*AINC
2128	! IF (XTIME.LT.10.)THEN
2129	! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,
2130	! * FABEOLD,AINCOLD
2131	! ENDIF
2132	DO NK=1,LTOP
2133	UMF(NK)=UMF2(NK)*AINC
2134	DMF(NK)=DMF2(NK)*AINC
2135	DETLQ(NK)=DETLQ2(NK)*AINC
2136	DETIC(NK)=DETIC2(NK)*AINC
2137	UDR(NK)=UDR2(NK)*AINC
2138	UER(NK)=UER2(NK)*AINC
2139	DER(NK)=DER2(NK)*AINC
2140	DDR(NK)=DDR2(NK)*AINC
2141	ENDDO
2142	!
2143	!...GO BACK UP FOR ANOTHER ITERATION...
2144	!
2145	ENDIF
2146	ENDDO iter
2147	!
2148	!...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
2149	!
2150	!...FRC2 IS THE FRACTION OF TOTAL CONDENSATE ! PPT FB MODS
2151	!...GENERATED THAT GOES INTO PRECIPITIATION ! PPT FB MODS
2152	!
2153	! Redistribute hydormeteors according to the final mass-flux values:
2154	!
2155	IF(CPR.GT.0.)THEN
2156	FRC2=PPTFLX/(CPR*AINC) ! PPT FB MODS
2157	ELSE
2158	FRC2=0.
2159	ENDIF
2160	DO NK=1,LTOP
2161	QLPA(NK)=QL0(NK)
2162	QIPA(NK)=QI0(NK)
2163	QRPA(NK)=QR0(NK)
2164	QSPA(NK)=QS0(NK)
2165	RAINFB(NK)=PPTLIQ(NK)AINCFBFRC*FRC2 ! PPT FB MODS
2166	SNOWFB(NK)=PPTICE(NK)AINCFBFRC*FRC2 ! PPT FB MODS
2167	ENDDO
2168	DO NTC=1,NSTEP
2169	!
2170	!...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH LAYER
2171	!...BASED ON THE SIGN OF OMEGA...
2172	!
2173	DO NK=1,LTOP
2174	QLFXIN(NK)=0.
2175	QLFXOUT(NK)=0.
2176	QIFXIN(NK)=0.
2177	QIFXOUT(NK)=0.
2178	QRFXIN(NK)=0.
2179	QRFXOUT(NK)=0.
2180	QSFXIN(NK)=0.
2181	QSFXOUT(NK)=0.
2182	ENDDO
2183	DO NK=2,LTOP
2184	IF(OMG(NK).LE.0.)THEN
2185	QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
2186	QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
2187	QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
2188	QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
2189	QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
2190	QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
2191	QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
2192	QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
2193	ELSE
2194	QLFXOUT(NK)=FXM(NK)*QLPA(NK)
2195	QIFXOUT(NK)=FXM(NK)*QIPA(NK)
2196	QRFXOUT(NK)=FXM(NK)*QRPA(NK)
2197	QSFXOUT(NK)=FXM(NK)*QSPA(NK)
2198	QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
2199	QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
2200	QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
2201	QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
2202	ENDIF
2203	ENDDO
2204	!
2205	!...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
2206	!
2207	DO NK=1,LTOP
2208	QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))DTIMEEMSD(NK)
2209	QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))DTIMEEMSD(NK)
2210	QRPA(NK)=QRPA(NK)+(QRFXIN(NK)-QRFXOUT(NK)+RAINFB(NK))DTIMEEMSD(NK) ! PPT FB MODS
2211	QSPA(NK)=QSPA(NK)+(QSFXIN(NK)-QSFXOUT(NK)+SNOWFB(NK))DTIMEEMSD(NK) ! PPT FB MODS
2212	ENDDO
2213	ENDDO
2214	DO NK=1,LTOP
2215	QLG(NK)=QLPA(NK)
2216	QIG(NK)=QIPA(NK)
2217	QRG(NK)=QRPA(NK)
2218	QSG(NK)=QSPA(NK)
2219	ENDDO
2220	!
2221	!...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
2222	!...GRID POINT...
2223	!
2224	! IF (XTIME.LT.10.)THEN
2225	! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2226	! ENDIF
2227	IF(IPRNT)THEN
2228	! WRITE(98,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2229	WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2230	CALL wrf_message(message)
2231	! call flush(98)
2232	endif
2233	!
2234	!...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
2235	!
2236	!297 IF(IPRNT)then
2237	IF(IPRNT)then
2238	! if(I.eq.16 .and. J.eq.41)then
2239	! IF(ISTOP.EQ.1)THEN
2240	write(98,*)
2241	! write(98,)'At t(h), I, J =',float(NTSD)72./3600.,I,J
2242	write(message,*)'P(LC), DTP, WKL, WKLCL =',p0(LC)/100., &
2243	TLCL+DTLCL+dtrh-TENV,WKL,WKLCL
2244	call wrf_message(message)
2245	write(message,*)'TLCL, DTLCL, DTRH, TENV =',TLCL,DTLCL, &
2246	DTRH,TENV
2247	call wrf_message(message)
2248	WRITE(message,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG, &
2249	TMIX-T00,PMIX,QMIX,ABE
2250	call wrf_message(message)
2251	WRITE(message,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100., &
2252	WLCL,CLDHGT(LC)
2253	call wrf_message(message)
2254	WRITE(message,1035)PEF,PEFCBH,LC,LET,WKL,VWS
2255	call wrf_message(message)
2256	write(message,*)'PRECIP EFFICIENCY =',PEFF
2257	call wrf_message(message)
2258	WRITE(message,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
2259	call wrf_message(message)
2260	! ENDIF
2261	!!!!! HERE !!!!!!!
2262	WRITE(message,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
2263	' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
2264	,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
2265	call wrf_message(message)
2266	write(message,*)'just before DO 300...'
2267	call wrf_message(message)
2268	! call flush(98)
2269	DO NK=1,LTOP
2270	K=LTOP-NK+1
2271	DTT=(TG(K)-T0(K))*86400./TIMEC
2272	RL=XLV0-XLV1*TG(K)
2273	DR=-(QG(K)-Q0(K))RL86400./(TIMEC*CP)
2274	UDFRC=UDR(K)TIMECEMSD(K)
2275	UEFRC=UER(K)TIMECEMSD(K)
2276	DDFRC=DDR(K)TIMECEMSD(K)
2277	DEFRC=-DER(K)TIMECEMSD(K)
2278	WRITE(message,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)*1.E4, &
2279	UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC,DDFRC,EMS(K)/1.E11, &
2280	W0AVG1D(K)1.E2,DETLQ(K)TIMECEMSD(K)1.E3,DETIC(K)* &
2281	TIMECEMSD(K)1.E3
2282	call wrf_message(message)
2283	ENDDO
2284	WRITE(message,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
2285	'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
2286	call wrf_message(message)
2287	DO NK=1,KL
2288	K=KX-NK+1
2289	DTT=TG(K)-T0(K)
2290	TUC=TU(K)-T00
2291	IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
2292	TDC=TZ(K)-T00
2293	IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
2294	IF(T0(K).LT.T00)THEN
2295	ES=ALIQEXP((BLIQTG(K)-CLIQ)/(TG(K)-DLIQ))
2296	ELSE
2297	ES=ALIQEXP((BLIQTG(K)-CLIQ)/(TG(K)-DLIQ))
2298	ENDIF
2299	QGS=ES*0.622/(P0(K)-ES)
2300	RH0=Q0(K)/QES(K)
2301	RHG=QG(K)/QGS
2302	WRITE(message,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC, &
2303	TDC,Q0(K)1000.,QG(K)1000.,(QG(K)-Q0(K))1000.,QU(K) &
2304	1000.,QD(K)1000.,QLG(K)1000.,QIG(K)1000.,QRG(K)1000., &
2305	QSG(K)*1000.,RH0,RHG
2306	call wrf_message(message)
2307	ENDDO
2308	!
2309	!...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
2310	!...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
2311	!
2312	! IF(ISTOP.EQ.1 .or. ISHALL.EQ.1)THEN
2313
2314	! IF(ISHALL.NE.1)THEN
2315	! write(98,4421)i,j,iyr,imo,idy,ihr,imn
2316	! write(98)i,j,iyr,imo,idy,ihr,imn,kl
2317	! 4421 format(7i4)
2318	! write(98,4422)kl
2319	! 4422 format(i6)
2320	DO 310 NK = 1,KL
2321	k = kl - nk + 1
2322	write(98,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
2323	u0(k),v0(k),W0AVG1D(K),dp(k),tke(k)
2324	! write(98) p0,t0,q0,u0,v0,w0,dp,tke
2325	! WRITE(98,1115)Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000.,
2326	! * U0(K),V0(K),DP(K)/100.,W0AVG(I,J,K)
2327	310 CONTINUE
2328	IF(ISTOP.EQ.1)THEN
2329	CALL wrf_error_fatal ( 'KAIN-FRITSCH, istop=1, diags' )
2330	ENDIF
2331	! ENDIF
2332	4455 format(8f11.3)
2333	ENDIF
2334	CNDTNF=(1.-EQFRC(LFS))(QLIQ(LFS)+QICE(LFS))DMF(LFS)
2335	PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
2336	RAINCV(I,J)=DT*PRATEC(I,J) ! PPT FB MODS
2337	! RAINCV(I,J)=.1.5DT*PPTFLX/DXSQ ! PPT FB MODS
2338	! RNC=0.1TIMECPPTFLX/DXSQ
2339	RNC=RAINCV(I,J)*NIC
2340	IF(ISHALL.EQ.0.AND.IPRNT)write (98,909)I,J,RNC
2341
2342	! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
2343	!
2344	! EVALUATE MOISTURE BUDGET...
2345	!
2346
2347	QINIT=0.
2348	QFNL=0.
2349	DPT=0.
2350	DO 315 NK=1,LTOP
2351	DPT=DPT+DP(NK)
2352	QINIT=QINIT+Q0(NK)*EMS(NK)
2353	QFNL=QFNL+QG(NK)*EMS(NK)
2354	QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
2355	315 CONTINUE
2356	QFNL=QFNL+PPTFLXTIMEC(1.-FBFRC) ! PPT FB MODS
2357	! QFNL=QFNL+PPTFLX*TIMEC ! PPT FB MODS
2358	ERR2=(QFNL-QINIT)*100./QINIT
2359	IF(IPRNT)WRITE(98,1110)QINIT,QFNL,ERR2
2360	IF(ABS(ERR2).GT.0.05 .AND. ISTOP.EQ.0)THEN
2361	! write(99,*)'!!!!!!!! MOISTURE BUDGET ERROR IN KFPARA !!!'
2362	! WRITE(99,1110)QINIT,QFNL,ERR2
2363	IPRNT=.TRUE.
2364	ISTOP=1
2365	write(98,4422)kl
2366	4422 format(i6)
2367	DO 311 NK = 1,KL
2368	k = kl - nk + 1
2369	! write(99,4455) p0(k)/100.,t0(k)-273.16,q0(k)*1000., &
2370	! u0(k),v0(k),W0AVG1D(K),dp(k)
2371	! write(98) p0,t0,q0,u0,v0,w0,dp,tke
2372	! WRITE(98,1115)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
2373	! U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
2374	WRITE(98,4456)P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
2375	U0(K),V0(K),W0AVG1D(K),dp(k)/100.,tke(k)
2376	311 CONTINUE
2377	! call flush(98)
2378
2379	! GOTO 297
2380	! STOP 'QVERR'
2381	ENDIF
2382	1115 FORMAT (2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
2383	4456 format(8f12.3)
2384	IF(PPTFLX.GT.0.)THEN
2385	RELERR=ERR2QINIT/(PPTFLXTIMEC)
2386	ELSE
2387	RELERR=0.
2388	ENDIF
2389	IF(IPRNT)THEN
2390	WRITE(98,1120)RELERR
2391	WRITE(98,*)'TDER, CPR, TRPPT =', &
2392	TDER,CPRAINC,TRPPTAINC
2393	ENDIF
2394	!
2395	!...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
2396	!
2397	!...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
2398	!...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
2399	!
2400	IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
2401	NCA(I,J) = REAL(NIC)*DT
2402	IF(ISHALL.EQ.1)THEN
2403	TIMEC = 2400.
2404	NCA(I,J) = CUDT*60.
2405	NSHALL = NSHALL+1
2406	ENDIF
2407
2408	DO K=1,KX
2409	! IF(IMOIST(INEST).NE.2)THEN
2410	!
2411	!...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMATED
2412	!...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
2413	!...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
2414	!...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
2415	!...OF QG...
2416	!
2417	! RLC=XLV0-XLV1*TG(K)
2418	! RLS=XLS0-XLS1*TG(K)
2419	! CPM=CP(1.+0.887QG(K))
2420	! TG(K)=TG(K)-(RLC(QLG(K)+QRG(K))+RLS(QIG(K)+QSG(K)))/CPM
2421	! QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
2422	! DQLDT(I,J,NK)=0.
2423	! DQIDT(I,J,NK)=0.
2424	! DQRDT(I,J,NK)=0.
2425	! DQSDT(I,J,NK)=0.
2426	! ELSE
2427	!
2428	!...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
2429	!
2430	IF(.NOT. F_QI .and. warm_rain)THEN
2431
2432	CPM=CP(1.+0.887QG(K))
2433	TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2434	DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2435	DQIDT(K)=0.
2436	DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2437	DQSDT(K)=0.
2438	ELSEIF(.NOT. F_QI .and. .not. warm_rain)THEN
2439	!
2440	!...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROMETEORS
2441	!...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
2442	!
2443	CPM=CP(1.+0.887QG(K))
2444	IF(K.LE.ML)THEN
2445	TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
2446	ELSEIF(K.GT.ML)THEN
2447	TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
2448	ENDIF
2449	DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
2450	DQIDT(K)=0.
2451	DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
2452	DQSDT(K)=0.
2453	ELSEIF(F_QI) THEN
2454	!
2455	!...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE TENDENCIES
2456	!...OF HYDROMETEORS DIRECTLY...
2457	!
2458	DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
2459	DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
2460	DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
2461	IF (F_QS) THEN
2462	DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
2463	ELSE
2464	DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
2465	ENDIF
2466	ELSE
2467	! PRINT *,'THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED!'
2468	CALL wrf_error_fatal ( 'KAIN-FRITSCH, THIS COMBINATION OF IMOIST, IEXICE, IICE NOT ALLOWED' )
2469	ENDIF
2470	DTDT(K)=(TG(K)-T0(K))/TIMEC
2471	DQDT(K)=(QG(K)-Q0(K))/TIMEC
2472	ENDDO
2473	PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
2474	RAINCV(I,J)=DT*PRATEC(I,J)
2475	! RAINCV(I,J)=.1.5DT*PPTFLX/DXSQ ! PPT FB MODS
2476	! RNC=0.1TIMECPPTFLX/DXSQ
2477	RNC=RAINCV(I,J)*NIC
2478	909 FORMAT('AT I, J =',i3,1x,i3,' CONVECTIVE RAINFALL =',F8.4,' mm')
2479	! write (98,909)I,J,RNC
2480	! write (6,909)I,J,RNC
2481	! WRITE(98,*)'at NTSD =',NTSD,',No. of KF points activated =',
2482	! * NCCNT
2483	! call flush(98)
2484	1000 FORMAT(' ',10A8)
2485	1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
2486	1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
2487	1015 FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
2488	1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
2489	' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
2490	I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
2491	' CAPE=',0PF7.1)
2492	1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
2493	E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
2494	F8.1)
2495	1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL=' &
2496	,F6.3,'VWS=',F5.2)
2497	!1055 FORMAT('*** DEGREE OF STABILIZATION =',F5.3, &
2498	! ', NO MORE MASS FLUX IS ALLOWED!')
2499	!1060 FORMAT(' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED &
2500	! &DEGREE OF STABILIZATION! FABE= ',F6.4)
2501	1070 FORMAT (16A8)
2502	1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
2503	1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, TADVEC, NSTEP=', &
2504	2(1X,F5.0),I3,'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
2505	1085 FORMAT (A3,16A7,2A8)
2506	1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
2507	1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ',F10.0)
2508	1105 FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =',&
2509	E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'%')
2510	1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
2511	' TOTAL WATER CHANGE =',F8.2,'%')
2512	! 1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4)
2513	1120 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3,'%')
2514	!
2515	!-----------------------------------------------------------------------
2516	!--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------
2517	!-----------------------------------------------------------------------
2518	!
2519	CUTOP(I,J)=REAL(LTOP)
2520	CUBOT(I,J)=REAL(LCL)
2521	!
2522	!-----------------------------------------------------------------------
2523	END SUBROUTINE KF_eta_PARA
2524	!********************************************************************
2525	! ***********************************************************************
2526	SUBROUTINE TPMIX2(p,thes,tu,qu,qliq,qice,qnewlq,qnewic,XLV1,XLV0)
2527	!
2528	! Lookup table variables:
2529	! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
2530	! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2531	! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2532	! REAL, SAVE, DIMENSION(1:200) :: ALU
2533	! REAL, SAVE :: RDPR,RDTHK,PLUTOP
2534	! End of Lookup table variables:
2535	!-----------------------------------------------------------------------
2536	IMPLICIT NONE
2537	!-----------------------------------------------------------------------
2538	REAL, INTENT(IN ) :: P,THES,XLV1,XLV0
2539	REAL, INTENT(OUT ) :: QNEWLQ,QNEWIC
2540	REAL, INTENT(INOUT) :: TU,QU,QLIQ,QICE
2541	REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11, &
2542	TEMP,QS,QNEW,DQ,QTOT,RLL,CPP
2543	INTEGER :: IPTB,ITHTB
2544	!-----------------------------------------------------------------------
2545
2546	!c****** LOOKUP TABLE VARIABLES... **************************
2547	! parameter(kfnt=250,kfnp=220)
2548	!c
2549	! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp),
2550	! * alu(200),rdpr,rdthk,plutop
2551	!C***************************************************************
2552	!c
2553	!c***********************************************************************
2554	!c scaling pressure and tt table index
2555	!c***********************************************************************
2556	!c
2557	tp=(p-plutop)*rdpr
2558	qq=tp-aint(tp)
2559	iptb=int(tp)+1
2560
2561	!
2562	!***********************************************************************
2563	! base and scaling factor for the
2564	!***********************************************************************
2565	!
2566	! scaling the and tt table index
2567	bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
2568	tth=(thes-bth)*rdthk
2569	pp =tth-aint(tth)
2570	ithtb=int(tth)+1
2571	IF(IPTB.GE.220 .OR. IPTB.LE.1 .OR. ITHTB.GE.250 .OR. ITHTB.LE.1)THEN
2572	write(98,)'* OUT OF BOUNDS *******'
2573	! call flush(98)
2574	ENDIF
2575	!
2576	t00=ttab(ithtb ,iptb )
2577	t10=ttab(ithtb+1,iptb )
2578	t01=ttab(ithtb ,iptb+1)
2579	t11=ttab(ithtb+1,iptb+1)
2580	!
2581	q00=qstab(ithtb ,iptb )
2582	q10=qstab(ithtb+1,iptb )
2583	q01=qstab(ithtb ,iptb+1)
2584	q11=qstab(ithtb+1,iptb+1)
2585	!
2586	!***********************************************************************
2587	! parcel temperature
2588	!***********************************************************************
2589	!
2590	temp=(t00+(t10-t00)pp+(t01-t00)qq+(t00-t10-t01+t11)ppqq)
2591	!
2592	qs=(q00+(q10-q00)pp+(q01-q00)qq+(q00-q10-q01+q11)ppqq)
2593	!
2594	DQ=QS-QU
2595	IF(DQ.LE.0.)THEN
2596	QNEW=QU-QS
2597	QU=QS
2598	ELSE
2599	!
2600	! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
2601	! ADJUSTED...IF LIQUID WATER IS PRESENT, IT IS ALLOWED TO EVAPORATE
2602	!
2603	QNEW=0.
2604	QTOT=QLIQ+QICE
2605	!
2606	! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS AT ITS
2607	! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MIXING
2608	! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURATION
2609	! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPRIATE
2610	! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE MADE.
2611	!
2612	!...subsaturated values only occur in calculations involving various mixtures of
2613	!...updraft and environmental air for estimation of entrainment and detrainment.
2614	!...For these purposes, assume that reasonable estimates can be given using
2615	!...liquid water saturation calculations only - i.e., ignore the effect of the
2616	!...ice phase in this process only...will not affect conservative properties...
2617	!
2618	IF(QTOT.GE.DQ)THEN
2619	qliq=qliq-dq*qliq/(qtot+1.e-10)
2620	qice=qice-dq*qice/(qtot+1.e-10)
2621	QU=QS
2622	ELSE
2623	RLL=XLV0-XLV1*TEMP
2624	CPP=1004.5(1.+0.89QU)
2625	IF(QTOT.LT.1.E-10)THEN
2626	!
2627	!...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
2628	TEMP=TEMP+RLL*(DQ/(1.+DQ))/CPP
2629	ELSE
2630	!
2631	!...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURATION,
2632	! THE TEMPERATURE IS GIVEN BY:
2633	!
2634	TEMP=TEMP+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CPP
2635	QU=QU+QTOT
2636	QTOT=0.
2637	QLIQ=0.
2638	QICE=0.
2639	ENDIF
2640	ENDIF
2641	ENDIF
2642	TU=TEMP
2643	qnewlq=qnew
2644	qnewic=0.
2645	!
2646	END SUBROUTINE TPMIX2
2647	!******************************************************************************
2648	SUBROUTINE DTFRZNEW(TU,P,THTEU,QU,QFRZ,QICE,ALIQ,BLIQ,CLIQ,DLIQ)
2649	!-----------------------------------------------------------------------
2650	IMPLICIT NONE
2651	!-----------------------------------------------------------------------
2652	REAL, INTENT(IN ) :: P,QFRZ,ALIQ,BLIQ,CLIQ,DLIQ
2653	REAL, INTENT(INOUT) :: TU,THTEU,QU,QICE
2654	REAL :: RLC,RLS,RLF,CPP,A,DTFRZ,ES,QS,DQEVAP,PII
2655	!-----------------------------------------------------------------------
2656	!
2657	!...ALLOW THE FREEZING OF LIQUID WATER IN THE UPDRAFT TO PROCEED AS AN
2658	!...APPROXIMATELY LINEAR FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE
2659	!...TTFRZ TO TBFRZ...
2660	!...FOR COLDER TEMPERATURES, FREEZE ALL LIQUID WATER...
2661	!...THERMODYNAMIC PROPERTIES ARE STILL CALCULATED WITH RESPECT TO LIQUID WATER
2662	!...TO ALLOW THE USE OF LOOKUP TABLE TO EXTRACT TMP FROM THETAE...
2663	!
2664	RLC=2.5E6-2369.276*(TU-273.16)
2665	RLS=2833922.-259.532*(TU-273.16)
2666	RLF=RLS-RLC
2667	CPP=1004.5(1.+0.89QU)
2668	!
2669	! A = D(es)/DT IS THAT CALCULATED FROM BUCK (1981) EMPERICAL FORMULAS
2670	! FOR SATURATION VAPOR PRESSURE...
2671	!
2672	A=(CLIQ-BLIQDLIQ)/((TU-DLIQ)(TU-DLIQ))
2673	DTFRZ = RLFQFRZ/(CPP+RLSQU*A)
2674	TU = TU+DTFRZ
2675
2676	ES = ALIQEXP((BLIQTU-CLIQ)/(TU-DLIQ))
2677	QS = ES*0.622/(P-ES)
2678	!
2679	!...FREEZING WARMS THE AIR AND IT BECOMES UNSATURATED...ASSUME THAT SOME OF THE
2680	!...LIQUID WATER THAT IS AVAILABLE FOR FREEZING EVAPORATES TO MAINTAIN SATURA-
2681	!...TION...SINCE THIS WATER HAS ALREADY BEEN TRANSFERRED TO THE ICE CATEGORY,
2682	!...SUBTRACT IT FROM ICE CONCENTRATION, THEN SET UPDRAFT MIXING RATIO AT THE NEW
2683	!...TEMPERATURE TO THE SATURATION VALUE...
2684	!
2685	DQEVAP = QS-QU
2686	QICE = QICE-DQEVAP
2687	QU = QU+DQEVAP
2688	PII=(1.E5/P)*(0.2854(1.-0.28*QU))
2689	THTEU=TUPIIEXP((3374.6525/TU-2.5403)QU(1.+0.81*QU))
2690	!
2691	END SUBROUTINE DTFRZNEW
2692	! --------------------------------------------------------------------------------
2693
2694	SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
2695	QNEWIC,QLQOUT,QICOUT,G)
2696
2697	!-----------------------------------------------------------------------
2698	IMPLICIT NONE
2699	!-----------------------------------------------------------------------
2700	! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
2701	! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
2702	! CULATED USING THE EQUATION DQ=-RATEQDT, BUT TO SIMULATE A QUASI-
2703	! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
2704	! RESOLUTION THIS IS EXPRESSED AS Q=QEXP(-RATEDZ).
2705
2706	REAL, INTENT(IN ) :: G
2707	REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
2708	REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
2709	REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
2710	!
2711	QTOT=QLIQ+QICE
2712	QNEW=QNEWLQ+QNEWIC
2713	!
2714	! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY
2715	! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
2716	! LEVELS...
2717	!
2718	QEST=0.5*(QTOT+QNEW)
2719	G1=WTW+BOTERM-ENTERM-2.GDZ*QEST/1.5
2720	IF(G1.LT.0.0)G1=0.
2721	WAVG=0.5*(SQRT(WTW)+SQRT(G1))
2722	CONV=RATE*DZ/WAVG ! KF90 Eq. 9
2723	!
2724	! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
2725	! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
2726	! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
2727	! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
2728	!
2729	RATIO3=QNEWLQ/(QNEW+1.E-8)
2730	! OLDQ=QTOT
2731	QTOT=QTOT+0.6*QNEW
2732	OLDQ=QTOT
2733	RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-8)
2734	QTOT=QTOT*EXP(-CONV) ! KF90 Eq. 9
2735	!
2736	! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
2737	! PARCEL AT THIS LEVEL...
2738	!
2739	DQ=OLDQ-QTOT
2740	QLQOUT=RATIO4*DQ
2741	QICOUT=(1.-RATIO4)*DQ
2742	!
2743	! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
2744	! LATE VERTICAL VELOCITY
2745	!
2746	PPTDRG=0.5(OLDQ+QTOT-0.2QNEW)
2747	WTW=WTW+BOTERM-ENTERM-2.GDZ*PPTDRG/1.5
2748	IF(ABS(WTW).LT.1.E-4)WTW=1.E-4
2749	!
2750	! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
2751	! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
2752	!
2753	QLIQ=RATIO4QTOT+RATIO30.4*QNEW
2754	QICE=(1.-RATIO4)QTOT+(1.-RATIO3)0.4*QNEW
2755	QNEWLQ=0.
2756	QNEWIC=0.
2757
2758	END SUBROUTINE CONDLOAD
2759
2760	! ----------------------------------------------------------------------
2761	SUBROUTINE PROF5(EQ,EE,UD)
2762	!
2763	!***********************************************************************
2764	!*** GAUSSIAN TYPE MIXING PROFILE....****************************
2765	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2766	! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
2767	! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN FROM
2768	! "HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICS TABLES"
2769	! ED. BY ABRAMOWITZ AND STEGUN, NATL BUREAU OF STANDARDS APPLIED
2770	! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
2771	! JACK KAIN
2772	! 7/6/89
2773	! Solves for KF90 Eq. 2
2774	!
2775	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2776	!-----------------------------------------------------------------------
2777	IMPLICIT NONE
2778	!-----------------------------------------------------------------------
2779	REAL, INTENT(IN ) :: EQ
2780	REAL, INTENT(INOUT) :: EE,UD
2781	REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
2782
2783	DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
2784	0.9372980,0.33267,0.166666667,0.202765151/
2785	X=(EQ-0.5)/SIGMA
2786	Y=6.*EQ-3.
2787	EY=EXP(Y*Y/(-2))
2788	E45=EXP(-4.5)
2789	T2=1./(1.+P*ABS(Y))
2790	T1=0.500498
2791	C1=A1T1+A2T1T1+A3T1T1T1
2792	C2=A1T2+A2T2T2+A3T2T2T2
2793	IF(Y.GE.0.)THEN
2794	EE=SIGMA(0.5(SQRT2P-E45C1-EYC2)+SIGMA(E45-EY))-E45EQ*EQ/2.
2795	UD=SIGMA(0.5(EYC2-E45C1)+SIGMA(E45-EY))-E45(0.5+EQ*EQ/2.- &
2796	EQ)
2797	ELSE
2798	EE=SIGMA(0.5(EYC2-E45C1)+SIGMA(E45-EY))-E45EQ*EQ/2.
2799	UD=SIGMA(0.5(SQRT2P-E45C1-EYC2)+SIGMA(E45-EY))-E45(0.5+EQ* &
2800	EQ/2.-EQ)
2801	ENDIF
2802	EE=EE/FE
2803	UD=UD/FE
2804
2805	END SUBROUTINE PROF5
2806
2807	! ------------------------------------------------------------------------
2808	SUBROUTINE TPMIX2DD(p,thes,ts,qs,i,j)
2809	!
2810	! Lookup table variables:
2811	! INTEGER, PARAMETER :: (KFNT=250,KFNP=220)
2812	! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
2813	! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
2814	! REAL, SAVE, DIMENSION(1:200) :: ALU
2815	! REAL, SAVE :: RDPR,RDTHK,PLUTOP
2816	! End of Lookup table variables:
2817	!-----------------------------------------------------------------------
2818	IMPLICIT NONE
2819	!-----------------------------------------------------------------------
2820	REAL, INTENT(IN ) :: P,THES
2821	REAL, INTENT(INOUT) :: TS,QS
2822	INTEGER, INTENT(IN ) :: i,j ! avail for debugging
2823	REAL :: TP,QQ,BTH,TTH,PP,T00,T10,T01,T11,Q00,Q10,Q01,Q11
2824	INTEGER :: IPTB,ITHTB
2825	CHARACTER*256 :: MESS
2826	!-----------------------------------------------------------------------
2827
2828	!
2829	!****** LOOKUP TABLE VARIABLES (F77 format)... **************************
2830	! parameter(kfnt=250,kfnp=220)
2831	!
2832	! COMMON/KFLUT/ ttab(kfnt,kfnp),qstab(kfnt,kfnp),the0k(kfnp), &
2833	! alu(200),rdpr,rdthk,plutop
2834	!***************************************************************
2835	!
2836	!***********************************************************************
2837	! scaling pressure and tt table index
2838	!***********************************************************************
2839	!
2840	tp=(p-plutop)*rdpr
2841	qq=tp-aint(tp)
2842	iptb=int(tp)+1
2843	!
2844	!***********************************************************************
2845	! base and scaling factor for the
2846	!***********************************************************************
2847	!
2848	! scaling the and tt table index
2849	bth=(the0k(iptb+1)-the0k(iptb))*qq+the0k(iptb)
2850	tth=(thes-bth)*rdthk
2851	pp =tth-aint(tth)
2852	ithtb=int(tth)+1
2853	!
2854	t00=ttab(ithtb ,iptb )
2855	t10=ttab(ithtb+1,iptb )
2856	t01=ttab(ithtb ,iptb+1)
2857	t11=ttab(ithtb+1,iptb+1)
2858	!
2859	q00=qstab(ithtb ,iptb )
2860	q10=qstab(ithtb+1,iptb )
2861	q01=qstab(ithtb ,iptb+1)
2862	q11=qstab(ithtb+1,iptb+1)
2863	!
2864	!***********************************************************************
2865	! parcel temperature and saturation mixing ratio
2866	!***********************************************************************
2867	!
2868	ts=(t00+(t10-t00)pp+(t01-t00)qq+(t00-t10-t01+t11)ppqq)
2869	!
2870	qs=(q00+(q10-q00)pp+(q01-q00)qq+(q00-q10-q01+q11)ppqq)
2871	!
2872	END SUBROUTINE TPMIX2DD
2873
2874	! -----------------------------------------------------------------------
2875	SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,ALIQ,BLIQ,CLIQ,DLIQ)
2876	!
2877	!-----------------------------------------------------------------------
2878	IMPLICIT NONE
2879	!-----------------------------------------------------------------------
2880	REAL, INTENT(IN ) :: P1,T1,Q1,ALIQ,BLIQ,CLIQ,DLIQ
2881	REAL, INTENT(INOUT) :: THT1
2882	REAL :: EE,TLOG,ASTRT,AINC,A1,TP,VALUE,AINTRP,TDPT,TSAT,THT, &
2883	T00,P00,C1,C2,C3,C4,C5
2884	INTEGER :: INDLU
2885	!-----------------------------------------------------------------------
2886	DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834, &
2887	0.278296,1.0723E-3/
2888	!
2889	! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
2890	!
2891	! NOTE: Calculations for mixed/ice phase no longer used...jsk 8/00
2892	! For example, KF90 Eq. 10 no longer used
2893	!
2894	EE=Q1*P1/(0.622+Q1)
2895	! TLOG=ALOG(EE/ALIQ)
2896	! ...calculate LOG term using lookup table...
2897	!
2898	astrt=1.e-3
2899	ainc=0.075
2900	a1=ee/aliq
2901	tp=(a1-astrt)/ainc
2902	indlu=int(tp)+1
2903	value=(indlu-1)*ainc+astrt
2904	aintrp=(a1-value)/ainc
2905	tlog=aintrpalu(indlu+1)+(1-aintrp)alu(indlu)
2906	!
2907	TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2908	TSAT=TDPT-(.212+1.571E-3(TDPT-T00)-4.36E-4(T1-T00))*(T1-TDPT)
2909	THT=T1(P00/P1)(0.2854(1.-0.28*Q1))
2910	THT1=THTEXP((C1/TSAT-C2)Q1(1.+0.81Q1))
2911	!
2912	END SUBROUTINE ENVIRTHT
2913	! ***********************************************************************
2914	!====================================================================
2915	SUBROUTINE kf_eta_init(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
2916	RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, &
2917	SVP1,SVP2,SVP3,SVPT0, &
2918	P_FIRST_SCALAR,restart,allowed_to_read, &
2919	ids, ide, jds, jde, kds, kde, &
2920	ims, ime, jms, jme, kms, kme, &
2921	its, ite, jts, jte, kts, kte )
2922	!--------------------------------------------------------------------
2923	IMPLICIT NONE
2924	!--------------------------------------------------------------------
2925	LOGICAL , INTENT(IN) :: restart,allowed_to_read
2926	INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
2927	ims, ime, jms, jme, kms, kme, &
2928	its, ite, jts, jte, kts, kte
2929	INTEGER , INTENT(IN) :: P_QI,P_QS,P_FIRST_SCALAR
2930
2931	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
2932	RTHCUTEN, &
2933	RQVCUTEN, &
2934	RQCCUTEN, &
2935	RQRCUTEN, &
2936	RQICUTEN, &
2937	RQSCUTEN
2938
2939	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
2940
2941	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
2942
2943	INTEGER :: i, j, k, itf, jtf, ktf
2944	REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
2945
2946	jtf=min0(jte,jde-1)
2947	ktf=min0(kte,kde-1)
2948	itf=min0(ite,ide-1)
2949
2950	IF(.not.restart)THEN
2951
2952	DO j=jts,jtf
2953	DO k=kts,ktf
2954	DO i=its,itf
2955	RTHCUTEN(i,k,j)=0.
2956	RQVCUTEN(i,k,j)=0.
2957	RQCCUTEN(i,k,j)=0.
2958	RQRCUTEN(i,k,j)=0.
2959	ENDDO
2960	ENDDO
2961	ENDDO
2962
2963	IF (P_QI .ge. P_FIRST_SCALAR) THEN
2964	DO j=jts,jtf
2965	DO k=kts,ktf
2966	DO i=its,itf
2967	RQICUTEN(i,k,j)=0.
2968	ENDDO
2969	ENDDO
2970	ENDDO
2971	ENDIF
2972
2973	IF (P_QS .ge. P_FIRST_SCALAR) THEN
2974	DO j=jts,jtf
2975	DO k=kts,ktf
2976	DO i=its,itf
2977	RQSCUTEN(i,k,j)=0.
2978	ENDDO
2979	ENDDO
2980	ENDDO
2981	ENDIF
2982
2983	DO j=jts,jtf
2984	DO i=its,itf
2985	NCA(i,j)=-100.
2986	ENDDO
2987	ENDDO
2988
2989	DO j=jts,jtf
2990	DO k=kts,ktf
2991	DO i=its,itf
2992	W0AVG(i,k,j)=0.
2993	ENDDO
2994	ENDDO
2995	ENDDO
2996
2997	endif
2998
2999	CALL KF_LUTAB(SVP1,SVP2,SVP3,SVPT0)
3000
3001	END SUBROUTINE kf_eta_init
3002
3003	!-------------------------------------------------------
3004
3005	subroutine kf_lutab(SVP1,SVP2,SVP3,SVPT0)
3006	!
3007	! This subroutine is a lookup table.
3008	! Given a series of series of saturation equivalent potential
3009	! temperatures, the temperature is calculated.
3010	!
3011	!--------------------------------------------------------------------
3012	IMPLICIT NONE
3013	!--------------------------------------------------------------------
3014	! Lookup table variables
3015	! INTEGER, SAVE, PARAMETER :: KFNT=250,KFNP=220
3016	! REAL, SAVE, DIMENSION(1:KFNT,1:KFNP) :: TTAB,QSTAB
3017	! REAL, SAVE, DIMENSION(1:KFNP) :: THE0K
3018	! REAL, SAVE, DIMENSION(1:200) :: ALU
3019	! REAL, SAVE :: RDPR,RDTHK,PLUTOP
3020	! End of Lookup table variables
3021
3022	INTEGER :: KP,IT,ITCNT,I
3023	REAL :: DTH,TMIN,TOLER,PBOT,DPR, &
3024	TEMP,P,ES,QS,PI,THES,TGUES,THGUES,F0,T1,T0,THGS,F1,DT, &
3025	ASTRT,AINC,A1,THTGS
3026	! REAL :: ALIQ,BLIQ,CLIQ,DLIQ,SVP1,SVP2,SVP3,SVPT0
3027	REAL :: ALIQ,BLIQ,CLIQ,DLIQ
3028	REAL, INTENT(IN) :: SVP1,SVP2,SVP3,SVPT0
3029	!
3030	! equivalent potential temperature increment
3031	data dth/1./
3032	! minimum starting temp
3033	data tmin/150./
3034	! tolerance for accuracy of temperature
3035	data toler/0.001/
3036	! top pressure (pascals)
3037	plutop=5000.0
3038	! bottom pressure (pascals)
3039	pbot=110000.0
3040
3041	ALIQ = SVP1*1000.
3042	BLIQ = SVP2
3043	CLIQ = SVP2*SVPT0
3044	DLIQ = SVP3
3045
3046	!
3047	! compute parameters
3048	!
3049	! 1._over_(sat. equiv. theta increment)
3050	rdthk=1./dth
3051	! pressure increment
3052	!
3053	DPR=(PBOT-PLUTOP)/REAL(KFNP-1)
3054	! dpr=(pbot-plutop)/REAL(kfnp-1)
3055	! 1._over_(pressure increment)
3056	rdpr=1./dpr
3057	! compute the spread of thes
3058	! thespd=dth*(kfnt-1)
3059	!
3060	! calculate the starting sat. equiv. theta
3061	!
3062	temp=tmin
3063	p=plutop-dpr
3064	do kp=1,kfnp
3065	p=p+dpr
3066	es=aliqexp((bliqtemp-cliq)/(temp-dliq))
3067	qs=0.622*es/(p-es)
3068	pi=(1.e5/p)*(0.2854(1.-0.28*qs))
3069	the0k(kp)=temppiexp((3374.6525/temp-2.5403)qs &
3070	(1.+0.81*qs))
3071	enddo
3072	!
3073	! compute temperatures for each sat. equiv. potential temp.
3074	!
3075	p=plutop-dpr
3076	do kp=1,kfnp
3077	thes=the0k(kp)-dth
3078	p=p+dpr
3079	do it=1,kfnt
3080	! define sat. equiv. pot. temp.
3081	thes=thes+dth
3082	! iterate to find temperature
3083	! find initial guess
3084	if(it.eq.1) then
3085	tgues=tmin
3086	else
3087	tgues=ttab(it-1,kp)
3088	endif
3089	es=aliqexp((bliqtgues-cliq)/(tgues-dliq))
3090	qs=0.622*es/(p-es)
3091	pi=(1.e5/p)*(0.2854(1.-0.28*qs))
3092	thgues=tguespiexp((3374.6525/tgues-2.5403)qs &
3093	(1.+0.81*qs))
3094	f0=thgues-thes
3095	t1=tgues-0.5*f0
3096	t0=tgues
3097	itcnt=0
3098	! iteration loop
3099	do itcnt=1,11
3100	es=aliqexp((bliqt1-cliq)/(t1-dliq))
3101	qs=0.622*es/(p-es)
3102	pi=(1.e5/p)*(0.2854(1.-0.28*qs))
3103	thtgs=t1piexp((3374.6525/t1-2.5403)qs(1.+0.81*qs))
3104	f1=thtgs-thes
3105	if(abs(f1).lt.toler)then
3106	exit
3107	endif
3108	! itcnt=itcnt+1
3109	dt=f1*(t1-t0)/(f1-f0)
3110	t0=t1
3111	f0=f1
3112	t1=t1-dt
3113	enddo
3114	ttab(it,kp)=t1
3115	qstab(it,kp)=qs
3116	enddo
3117	enddo
3118	!
3119	! lookup table for tlog(emix/aliq)
3120	!
3121	! set up intial values for lookup tables
3122	!
3123	astrt=1.e-3
3124	ainc=0.075
3125	!
3126	a1=astrt-ainc
3127	do i=1,200
3128	a1=a1+ainc
3129	alu(i)=alog(a1)
3130	enddo
3131	!
3132	END SUBROUTINE KF_LUTAB
3133
3134	END MODULE module_cu_kfeta

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

Download in other formats:

Original Format