Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_kfeta.F @ 3547

Last change on this file since 3547 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 101.2 KB

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

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

Download in other formats:

Original Format