Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

convection_kfeta.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago
sources for flexwrf v3.1
File size: 98.7 KB

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

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

Download in other formats:

Original Format