Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_kfeta.F @ 3431

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

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

Download in other formats:

Original Format