Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_kf.F @ 1

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 155.9 KB

Line
1	!WRF:MODEL_LAYER:PHYSICS
2	!
3
4	MODULE module_cu_kf
5
6	USE module_wrf_error
7
8	REAL , PARAMETER :: RAD = 1500.
9
10	CONTAINS
11
12	!-------------------------------------------------------------
13	SUBROUTINE KFCPS( &
14	ids,ide, jds,jde, kds,kde &
15	,ims,ime, jms,jme, kms,kme &
16	,its,ite, jts,jte, kts,kte &
17	,DT,KTAU,DX,CUDT,CURR_SECS,ADAPT_STEP_FLAG &
18	,rho &
19	,RAINCV,PRATEC,NCA &
20	,U,V,TH,T,W,QV,dz8w,Pcps,pi &
21	,W0AVG,XLV0,XLV1,XLS0,XLS1,CP,R,G,EP1 &
22	,EP2,SVP1,SVP2,SVP3,SVPT0 &
23	,STEPCU,CU_ACT_FLAG,warm_rain &
24	! optional arguments
25	,F_QV ,F_QC ,F_QR ,F_QI ,F_QS &
26	,RQVCUTEN,RQCCUTEN,RQRCUTEN,RQICUTEN,RQSCUTEN &
27	,RTHCUTEN &
28	)
29	!-------------------------------------------------------------
30	IMPLICIT NONE
31	!-------------------------------------------------------------
32	INTEGER, INTENT(IN ) :: &
33	ids,ide, jds,jde, kds,kde, &
34	ims,ime, jms,jme, kms,kme, &
35	its,ite, jts,jte, kts,kte
36
37	INTEGER, INTENT(IN ) :: STEPCU
38	LOGICAL, INTENT(IN ) :: warm_rain
39
40	REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1
41	REAL, INTENT(IN ) :: CP,R,G,EP1,EP2
42	REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
43
44	INTEGER, INTENT(IN ) :: KTAU
45
46	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
47	INTENT(IN ) :: &
48	U, &
49	V, &
50	W, &
51	TH, &
52	QV, &
53	T, &
54	dz8w, &
55	Pcps, &
56	rho, &
57	pi
58	!
59	REAL, INTENT(IN ) :: DT, DX
60	REAL, INTENT(IN ) :: CUDT
61	REAL, INTENT(IN ) :: CURR_SECS
62	LOGICAL,INTENT(IN ) :: ADAPT_STEP_FLAG
63
64	REAL, DIMENSION( ims:ime , jms:jme ), &
65	INTENT(INOUT) :: &
66	RAINCV &
67	,PRATEC &
68	, NCA
69
70	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
71	INTENT(INOUT) :: &
72	W0AVG
73
74	LOGICAL, DIMENSION( ims:ime , jms:jme ), &
75	INTENT(INOUT) :: CU_ACT_FLAG
76
77	!
78	! Optional arguments
79	!
80
81	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
82	OPTIONAL, &
83	INTENT(INOUT) :: &
84	RTHCUTEN &
85	,RQVCUTEN &
86	,RQCCUTEN &
87	,RQRCUTEN &
88	,RQICUTEN &
89	,RQSCUTEN
90
91	!
92	! Flags relating to the optional tendency arrays declared above
93	! Models that carry the optional tendencies will provdide the
94	! optional arguments at compile time; these flags all the model
95	! to determine at run-time whether a particular tracer is in
96	! use or not.
97	!
98
99	LOGICAL, OPTIONAL :: &
100	F_QV &
101	,F_QC &
102	,F_QR &
103	,F_QI &
104	,F_QS
105
106
107
108	! LOCAL VARS
109
110	REAL, DIMENSION( kts:kte ) :: &
111	U1D, &
112	V1D, &
113	T1D, &
114	DZ1D, &
115	QV1D, &
116	P1D, &
117	RHO1D, &
118	W0AVG1D
119
120	REAL, DIMENSION( kts:kte ):: &
121	DQDT, &
122	DQIDT, &
123	DQCDT, &
124	DQRDT, &
125	DQSDT, &
126	DTDT
127
128	REAL :: TST,tv,PRS,RHOE,W0,SCR1,DXSQ,tmp
129
130	INTEGER :: i,j,k,NTST,ICLDCK
131
132	LOGICAL :: qi_flag , qs_flag
133	! adjustable time step changes
134	REAL :: lastdt = -1.0
135	REAL :: W0AVGfctr, W0fctr, W0den
136	LOGICAL :: run_param
137
138	!----------------------------------------------------------------------
139
140	!--- CALL CUMULUS PARAMETERIZATION
141	!
142	!...TST IS THE NUMBER OF TIME STEPS IN 10 MINUTES...W0AVG IS CLOSE TO A
143	!...RUNNING MEAN VERTICAL VELOCITY...NOTE THAT IF YOU CHANGE TST, IT WIL
144	!...CHANGE THE FREQUENCY OF THE CONVECTIVE INTITIATION CHECK (SEE BELOW)
145	!...NOTE THAT THE ORDERING OF VERTICAL LAYERS MUST BE REVERSED FOR W0AVG
146	!...BECAUSE THE ORDERING IS REVERSED IN KFPARA...
147	!
148	DXSQ=DX*DX
149	qi_flag = .FALSE.
150	qs_flag = .FALSE.
151	IF ( PRESENT( F_QI ) ) qi_flag = f_qi
152	IF ( PRESENT( F_QS ) ) qs_flag = f_qs
153
154	!----------------------
155	NTST=STEPCU
156	TST=float(NTST*2)
157	!----------------------
158	! NTST=NINT(1200./(DT*2.))
159	! TST=float(NTST)
160	! NTST=NINT(0.5*TST)
161	! NTST=MAX0(NTST,1)
162	!----------------------
163	! ICLDCK=MOD(KTAU,NTST)
164	!----------------------
165	! write(0,*) 'DT = ',DT,' KTAU = ',KTAU,' DX = ',DX
166	! write(0,*) 'CUDT = ',CUDT,' CURR_SECS = ',CURR_SECS
167	! write(0,*) 'ADAPT_STEP_FLAG = ',ADAPT_STEP_FLAG,' IDS = ',IDS
168	! write(0,*) 'STEPCU = ',STEPCU,' warm_rain = ',warm_rain
169	! write(0,*) 'F_QV = ',F_QV,' F_QC = ',F_QV
170	! write(0,*) 'F_QI = ',F_QI,' F_QS = ',F_QS
171	! write(0,*) 'F_QR = ',F_QR
172	! stop
173	if (lastdt < 0) then
174	lastdt = dt
175	endif
176
177	if (ADAPT_STEP_FLAG) then
178	W0AVGfctr = 2 * MAX(CUDT*60,dt) - dt
179	W0fctr = dt
180	W0den = 2 * MAX(CUDT*60,dt)
181	else
182	W0AVGfctr = (TST-1.)
183	W0fctr = 1.
184	W0den = TST
185	endif
186
187	DO J = jts,jte
188	DO K=kts,kte
189	DO I= its,ite
190	! SCR1=-5.0E-4Grho(I,K,J)*(w(I,K,J)+w(I,K+1,J))
191	! TV=T(I,K,J)(1.+EP1QV(I,K,J))
192	! RHOE=Pcps(I,K,J)/(R*TV)
193	! W0=-101.9368*SCR1/RHOE
194	W0=0.5*(w(I,K,J)+w(I,K+1,J))
195
196	! Old:
197	!
198	! W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
199	! New, to support adaptive time step:
200	!
201	W0AVG(I,K,J) = ( W0AVG(I,K,J) * W0AVGfctr + W0 * W0fctr ) / W0den
202
203	ENDDO
204	ENDDO
205	ENDDO
206	lastdt = dt
207	!
208	!...CHECK FOR CONVECTIVE INITIATION EVERY 5 MINUTES (OR NTST/2)...
209	!
210
211	!
212	! Modified for adaptive time step
213	!
214	if (ADAPT_STEP_FLAG) then
215	if ( (KTAU .eq. 1) .or. (cudt .eq. 0) .or. &
216	( CURR_SECS + dt >= &
217	( int( CURR_SECS / ( cudt * 60 ) ) + 1 ) * cudt * 60 ) ) then
218	run_param = .TRUE.
219	else
220	run_param = .FALSE.
221	endif
222
223	else
224	if (MOD(KTAU,NTST) .EQ. 0 .or. KTAU .eq. 1) then
225	run_param = .TRUE.
226	else
227	run_param = .FALSE.
228	endif
229	endif
230
231	IF (run_param) then
232	DO J = jts,jte
233	DO I= its,ite
234	CU_ACT_FLAG(i,j) = .true.
235	ENDDO
236	ENDDO
237
238	DO J = jts,jte
239	DO I=its,ite
240	! if (i.eq. 110 .and. j .eq. 59 ) then
241	! write(0,*) 'nca = ',nca(i,j),' CU_ACT_FLAG = ',CU_ACT_FLAG(i,j)
242	! write(0,*) 'dt = ',dt,' ADAPT_STEP_FLAG = ',ADAPT_STEP_FLAG
243	! endif
244	! IF ( NINT(NCA(I,J)) .gt. 0 ) then
245	IF ( NCA(I,J) .gt. 0.5*DT ) then
246	CU_ACT_FLAG(i,j) = .false.
247	ELSE
248
249	DO k=kts,kte
250	DQDT(k)=0.
251	DQIDT(k)=0.
252	DQCDT(k)=0.
253	DQRDT(k)=0.
254	DQSDT(k)=0.
255	DTDT(k)=0.
256	ENDDO
257	RAINCV(I,J)=0.
258	PRATEC(I,J)=0.
259	!
260	! assign vars from 3D to 1D
261
262	DO K=kts,kte
263	U1D(K) =U(I,K,J)
264	V1D(K) =V(I,K,J)
265	T1D(K) =T(I,K,J)
266	RHO1D(K) =rho(I,K,J)
267	QV1D(K)=QV(I,K,J)
268	P1D(K) =Pcps(I,K,J)
269	W0AVG1D(K) =W0AVG(I,K,J)
270	DZ1D(k)=dz8w(I,K,J)
271	ENDDO
272
273	!
274	CALL KFPARA(I, J, &
275	U1D,V1D,T1D,QV1D,P1D,DZ1D, &
276	W0AVG1D,DT,DX,DXSQ,RHO1D, &
277	XLV0,XLV1,XLS0,XLS1,CP,R,G, &
278	EP2,SVP1,SVP2,SVP3,SVPT0, &
279	DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
280	RAINCV,PRATEC,NCA, &
281	warm_rain,qi_flag,qs_flag, &
282	ids,ide, jds,jde, kds,kde, &
283	ims,ime, jms,jme, kms,kme, &
284	its,ite, jts,jte, kts,kte )
285
286	IF ( PRESENT( RTHCUTEN ) .AND. PRESENT( RQVCUTEN ) ) THEN
287	DO K=kts,kte
288	RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
289	RQVCUTEN(I,K,J)=DQDT(K)
290	ENDDO
291	ENDIF
292
293	IF( PRESENT(RQRCUTEN) .AND. PRESENT(RQCCUTEN) .AND. &
294	PRESENT(F_QR) ) THEN
295	IF ( F_QR ) THEN
296	DO K=kts,kte
297	RQRCUTEN(I,K,J)=DQRDT(K)
298	RQCCUTEN(I,K,J)=DQCDT(K)
299	ENDDO
300	ELSE
301	! This is the case for Eta microphysics without 3d rain field
302	DO K=kts,kte
303	RQRCUTEN(I,K,J)=0.
304	RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
305	ENDDO
306	ENDIF
307	ENDIF
308
309	!...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
310
311	IF( PRESENT( RQICUTEN ) .AND. qi_flag )THEN
312	DO K=kts,kte
313	RQICUTEN(I,K,J)=DQIDT(K)
314	ENDDO
315	ENDIF
316
317	IF( PRESENT ( RQSCUTEN ) .AND. qs_flag )THEN
318	DO K=kts,kte
319	RQSCUTEN(I,K,J)=DQSDT(K)
320	ENDDO
321	ENDIF
322	!
323	ENDIF
324	ENDDO
325	ENDDO
326
327	ENDIF
328
329	END SUBROUTINE KFCPS
330
331	!-----------------------------------------------------------
332	SUBROUTINE KFPARA (I, J, &
333	U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
334	DT,DX,DXSQ,rho, &
335	XLV0,XLV1,XLS0,XLS1,CP,R,G, &
336	EP2,SVP1,SVP2,SVP3,SVPT0, &
337	DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
338	RAINCV,PRATEC,NCA, &
339	warm_rain,qi_flag,qs_flag, &
340	ids,ide, jds,jde, kds,kde, &
341	ims,ime, jms,jme, kms,kme, &
342	its,ite, jts,jte, kts,kte )
343	!-----------------------------------------------------------
344	IMPLICIT NONE
345	!-----------------------------------------------------------
346	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
347	ims,ime, jms,jme, kms,kme, &
348	its,ite, jts,jte, kts,kte, &
349	I,J
350	LOGICAL, INTENT(IN ) :: warm_rain
351	LOGICAL :: qi_flag, qs_flag
352
353	!
354	REAL, DIMENSION( kts:kte ), &
355	INTENT(IN ) :: U0, &
356	V0, &
357	T0, &
358	QV0, &
359	P0, &
360	rho, &
361	DZQ, &
362	W0AVG1D
363	!
364	REAL, INTENT(IN ) :: DT,DX,DXSQ
365	!
366
367	REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
368	REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
369	!
370	REAL, DIMENSION( kts:kte ), INTENT(INOUT) :: &
371	DQDT, &
372	DQIDT, &
373	DQCDT, &
374	DQRDT, &
375	DQSDT, &
376	DTDT
377
378	REAL, DIMENSION( ims:ime , jms:jme ), &
379	INTENT(INOUT) :: RAINCV, &
380	PRATEC, &
381	NCA
382	!
383	!...DEFINE LOCAL VARIABLES...
384	!
385	REAL, DIMENSION( kts:kte ) :: &
386	Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
387	QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
388	UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
389	UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
390	THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
391	QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
392	DETLQ2,DETIC2,RATIO,RATIO2
393
394	REAL, DIMENSION( kts:kte ) :: &
395	DOMGDP,EXN,RHOE,TVQU,DP,RH,EQFRC,WSPD, &
396	QDT,FXM,THTAG,THTESG,THPA,THFXTOP, &
397	THFXBOT,QPA,QFXTOP,QFXBOT,QLPA,QLFXIN, &
398	QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
399	QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
400	QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
401
402	REAL, DIMENSION( kts:kte+1 ) :: OMG
403	REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
404
405	! LOCAL VARS
406
407	REAL :: P00,T00,CV,B61,RLF,RHIC,RHBC,PIE, &
408	TTFRZ,TBFRZ,C5,RATE
409	REAL :: GDRY,ROCP,ALIQ,BLIQ, &
410	CLIQ,DLIQ,AICE,BICE,CICE,DICE
411	REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
412	ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
413	CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
414	ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
415	TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
416	UPNEW,ABE,WKLCL,THTUDL,TUDL,TTEMP,FRC1, &
417	QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
418	DZZ,WSQ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
419	THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
420	UD1,CLDHGT,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
421	THTA,P150,USR,VCONV,TIMEC,SHSIGN,VWS,PEF, &
422	CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
423	DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
424	DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
425	UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
426	DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
427	AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
428	DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
429	TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
430	UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
431	RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
432	DDFRC,TDC,DEFRC
433
434	INTEGER :: KX,K,KL
435	!
436	INTEGER :: ISTOP,ML,L5,L4,KMIX,LOW, &
437	LC,MXLAYR,LLFC,NLAYRS,NK, &
438	KPBL,KLCL,LCL,LET,IFLAG, &
439	KFRZ,NK1,LTOP,NJ,LTOP1, &
440	LTOPM1,LVF,KSTART,KMIN,LFS, &
441	ND,NIC,LDB,LDT,ND1,NDK, &
442	NM,LMAX,NCOUNT,NOITR, &
443	NSTEP,NTC
444	!
445	DATA P00,T00/1.E5,273.16/
446	DATA CV,B61,RLF/717.,0.608,3.339E5/
447	DATA RHIC,RHBC/1.,0.90/
448	DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
449	DATA RATE/0.01/
450	!-----------------------------------------------------------
451	GDRY=-G/CP
452	ROCP=R/CP
453	KL=kte
454	KX=kte
455	!
456	! ALIQ = 613.3
457	! BLIQ = 17.502
458	! CLIQ = 4780.8
459	! DLIQ = 32.19
460	ALIQ = SVP1*1000.
461	BLIQ = SVP2
462	CLIQ = SVP2*SVPT0
463	DLIQ = SVP3
464	AICE = 613.2
465	BICE = 22.452
466	CICE = 6133.0
467	DICE = 0.61
468	!
469
470	!...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER
471	!...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL)
472	!...FIELD. 'FBFRC' IS THE FRACTION OF AVAILABLE
473	!...PRECIPITATION TO BE FED BACK (0.0 - 1.0)...
474	!
475	FBFRC=0.0
476	!
477	!...SCHEME IS CALLED ONCE ON EACH NORTH-SOUTH SLICE, THE LOOP BELOW
478	!...CHECKS FOR THE POSSIBILITY OF INITIATING PARAMETERIZED
479	!...CONVECTION AT EACH POINT WITHIN THE SLICE
480	!
481	!...SEE IF IT IS NECESSARY TO CHECK FOR CONVECTIVE TRIGGERING AT THIS
482	!...GRID POINT. IF NCA>0, CONVECTION IS ALREADY ACTIVE AT THIS POINT,
483	!...JUST FEED BACK THE TENDENCIES SAVED FROM THE TIME WHEN CONVECTION
484	!...WAS INITIATED. IF NCA<0, CONVECTION IS NOT ACTIVE
485	!...AND YOU MAY WANT TO CHECK TO SEE IF IT CAN BE ACTIVATED FOR THE
486	!...CURRENT CONDITIONS. IN PREVIOUS APLICATIONS OF THIS SCHEME,
487	!...THE VARIABLE ICLDCK WAS USED BELOW TO SAVE TIME BY ONLY CHECKING
488	!...FOR THE POSSIBILITY OF CONVECTIVE INITIATION EVERY 5 OR 10
489	!...MINUTES...
490	!
491
492	! 10 CONTINUE
493	!SUE P300=1000.(PSB(I,J)A(KL)+PTOP-30.)+PP3D(I,J,KL)
494
495	P300=P0(1)-30000.
496	!
497	!...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
498	!...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
499	!...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
500	!
501	!...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
502	!...FROM BOTTOM-UP IN THE KF SCHEME...
503	!
504	ML=0
505	!SUE tmprpsb=1./PSB(I,J)
506	!SUE CELL=PTOP*tmprpsb
507
508	DO 15 K=1,KX
509	!SUE P0(K)=1.E3(A(NK)PSB(I,J)+PTOP)+PP3D(I,J,NK)
510	!
511	!...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
512	!
513	ES=ALIQEXP((BLIQT0(K)-CLIQ)/(T0(K)-DLIQ))
514	QES(K)=EP2*ES/(P0(K)-ES)
515	Q0(K)=AMIN1(QES(K),QV0(K))
516	Q0(K)=AMAX1(0.000001,Q0(K))
517	QL0(K)=0.
518	QI0(K)=0.
519	QR0(K)=0.
520	QS0(K)=0.
521
522	TV0(K)=T0(K)(1.+B61Q0(K))
523	RHOE(K)=P0(K)/(R*TV0(K))
524
525	DP(K)=rho(k)gDZQ(k)
526	!
527	!...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
528	! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
529	!
530	IF(P0(K).GE.500E2)L5=K
531	IF(P0(K).GE.400E2)L4=K
532	IF(P0(K).GE.P300)LLFC=K
533	IF(T0(K).GT.T00)ML=K
534	15 CONTINUE
535
536	Z0(1)=.5*DZQ(1)
537	DO 20 K=2,KL
538	Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
539	DZA(K-1)=Z0(K)-Z0(K-1)
540	20 CONTINUE
541	DZA(KL)=0.
542	KMIX=1
543	25 LOW=KMIX
544
545	IF(LOW.GT.LLFC)GOTO 325
546
547	LC=LOW
548	MXLAYR=0
549	!
550	!...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
551	!...UNSTABLE AIR 50 TO 100 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
552	!...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
553	!...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 60 mb..
554	!
555	NLAYRS=0
556	DPTHMX=0.
557	DO 63 NK=LC,KX
558	DPTHMX=DPTHMX+DP(NK)
559	NLAYRS=NLAYRS+1
560	63 IF(DPTHMX.GT.6.E3)GOTO 64
561	GOTO 325
562	64 KPBL=LC+NLAYRS-1
563	KMIX=LC+1
564	18 THMIX=0.
565	QMIX=0.
566	ZMIX=0.
567	PMIX=0.
568	DPTHMX=0.
569	!
570	!...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
571	!...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
572	!...LAYERS...
573	!
574	DO 17 NK=LC,KPBL
575	DPTHMX=DPTHMX+DP(NK)
576	ROCPQ=0.2854(1.-0.28Q0(NK))
577	THMIX=THMIX+DP(NK)T0(NK)(P00/P0(NK))**ROCPQ
578	QMIX=QMIX+DP(NK)*Q0(NK)
579	ZMIX=ZMIX+DP(NK)*Z0(NK)
580	17 PMIX=PMIX+DP(NK)*P0(NK)
581	THMIX=THMIX/DPTHMX
582	QMIX=QMIX/DPTHMX
583	ZMIX=ZMIX/DPTHMX
584	PMIX=PMIX/DPTHMX
585	ROCPQ=0.2854(1.-0.28QMIX)
586	TMIX=THMIX(PMIX/P00)*ROCPQ
587	EMIX=QMIX*PMIX/(EP2+QMIX)
588	!
589	!...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL, PRESSURE
590	!...LEVEL OF LCL...
591	!
592	TLOG=ALOG(EMIX/ALIQ)
593	TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
594	TLCL=TDPT-(.212+1.571E-3(TDPT-T00)-4.36E-4(TMIX-T00))*(TMIX- &
595	TDPT)
596	TLCL=AMIN1(TLCL,TMIX)
597	TVLCL=TLCL(1.+0.608QMIX)
598	CPORQ=1./ROCPQ
599	PLCL=P00(TLCL/THMIX)*CPORQ
600	DO 29 NK=LC,KL
601	KLCL=NK
602	IF(PLCL.GE.P0(NK))GOTO 35
603	29 CONTINUE
604	GOTO 325
605	35 K=KLCL-1
606	DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
607	!
608	!...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
609	!
610	TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
611	QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
612	TVEN=TENV(1.+0.608QENV)
613	TVBAR=0.5*(TV0(K)+TVEN)
614	! ZLCL=Z0(K)+RTVBARALOG(P0(K)/PLCL)/G
615	ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
616	!
617	!...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
618	!...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0AVG IS AN
619	!...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
620	!...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
621	!...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
622	!...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
623	!...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
624	!...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
625	!...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
626	!
627	WKLCL=0.02*ZLCL/2.5E3
628	WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))DLP)DX/25.E3- &
629	WKLCL
630	WABS=ABS(WKL)+1.E-10
631	WSIGNE=WKL/WABS
632	DTLCL=4.64WSIGNEWABS**0.33
633	GDT=GDTLCL(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
634	WLCL=1.+.5WSIGNESQRT(ABS(GDT)+1.E-10)
635	IF(TLCL+DTLCL.GT.TENV)GOTO 45
636	IF(KPBL.GE.LLFC)GOTO 325
637	GOTO 25
638	!
639	!...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
640	!...EQUIVALENT POTENTIAL TEMPERATURE
641	!...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
642	!
643	45 THETEU(K)=TMIX(1.E5/PMIX)(0.2854(1.-0.28QMIX)) &
644	EXP((3374.6525/TLCL-2.5403)QMIX(1.+0.81*QMIX))
645	ES=ALIQEXP((TENVBLIQ-CLIQ)/(TENV-DLIQ))
646	TVAVG=0.5(TV0(KLCL)+TENV(1.+0.608*QENV))
647	PLCL=P0(KLCL)EXP(G/(RTVAVG)*(Z0(KLCL)-ZLCL))
648	QESE=EP2*ES/(PLCL-ES)
649	GDT=GDTLCL(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
650	WLCL=1.+.5WSIGNESQRT(ABS(GDT)+1.E-10)
651	THTES(K)=TENV(1.E5/PLCL)(0.2854(1.-0.28QESE)) &
652	EXP((3374.6525/TENV-2.5403)QESE(1.+0.81*QESE))
653	WTW=WLCL*WLCL
654	IF(WLCL.LT.0.)GOTO 25
655	TVLCL=TLCL(1.+0.608QMIX)
656	RHOLCL=PLCL/(R*TVLCL)
657	!
658	LCL=KLCL
659	LET=LCL
660	!
661	!*******************************************************************
662	! *
663	! COMPUTE UPDRAFT PROPERTIES *
664	! *
665	!*******************************************************************
666	!
667	!
668	!...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
669	!
670	WU(K)=WLCL
671	AU0=PIERADRAD
672	UMF(K)=RHOLCL*AU0
673	VMFLCL=UMF(K)
674	UPOLD=VMFLCL
675	UPNEW=UPOLD
676	!
677	!...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
678	!...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE BUOYANT ENERGY,
679	! TRPPT IS THE TOTAL RATE OF PRECIPITATION PRODUCTION...
680	!
681	RATIO2(K)=0.
682	UER(K)=0.
683	ABE=0.
684	TRPPT=0.
685	TU(K)=TLCL
686	TVU(K)=TVLCL
687	QU(K)=QMIX
688	EQFRC(K)=1.
689	QLIQ(K)=0.
690	QICE(K)=0.
691	QLQOUT(K)=0.
692	QICOUT(K)=0.
693	DETLQ(K)=0.
694	DETIC(K)=0.
695	PPTLIQ(K)=0.
696	PPTICE(K)=0.
697	IFLAG=0
698	KFRZ=LC
699	!
700	!...THE AMOUNT OF CONV AVAIL POT ENERGY (CAPE) IS CALCULATED WITH
701	! RESPECT TO UNDILUTE PARCEL ASCENT; EQ POT TEMP OF UNDILUTE
702	! PARCEL IS THTUDL, UNDILUTE TEMPERATURE IS GIVEN BY TUDL...
703	!
704	THTUDL=THETEU(K)
705	TUDL=TLCL
706	!
707	!...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
708	! PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
709	! FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
710	! INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
711	! PREVIOUS MODEL LEVEL...
712	!
713	TTEMP=TTFRZ
714	!
715	!...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
716	! MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
717	! MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
718	!
719	DO 60 NK=K,KL-1
720	NK1=NK+1
721	RATIO2(NK1)=RATIO2(NK)
722	!
723	!...UPDATE UPDRAFT PROPERTIES AT THE NEXT MODEL LVL TO REFLECT
724	! ENTRAINMENT OF ENVIRONMENTAL AIR...
725	!
726	FRC1=0.
727	TU(NK1)=T0(NK1)
728	THETEU(NK1)=THETEU(NK)
729	QU(NK1)=QU(NK)
730	QLIQ(NK1)=QLIQ(NK)
731	QICE(NK1)=QICE(NK)
732
733	CALL TPMIX(P0(NK1),THETEU(NK1),TU(NK1),QU(NK1),QLIQ(NK1), &
734	QICE(NK1),QNEWLQ,QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0, &
735	XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
736	TVU(NK1)=TU(NK1)(1.+0.608QU(NK1))
737	!
738	!...CHECK TO SEE IF UPDRAFT TEMP IS WITHIN THE FREEZING INTERVAL,
739	! IF IT IS, CALCULATE THE FRACTIONAL CONVERSION TO GLACIATION
740	! AND ADJUST QNEWLQ TO REFLECT THE GRADUAL CHANGE IN THETAU
741	! SINCE THE LAST MODEL LEVEL...THE GLACIATION EFFECTS WILL BE
742	! DETERMINED AFTER THE AMOUNT OF CONDENSATE AVAILABLE AFTER
743	! PRECIP FALLOUT IS DETERMINED...TTFRZ IS THE TEMP AT WHICH
744	! GLACIATION BEGINS, TBFRZ THE TEMP AT WHICH IT ENDS...
745	!
746	IF(TU(NK1).LE.TTFRZ.AND.IFLAG.LT.1)THEN
747	IF(TU(NK1).GT.TBFRZ)THEN
748	IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
749	FRC1=(TTEMP-TU(NK1))/(TTFRZ-TBFRZ)
750	R1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
751	ELSE
752	FRC1=(TTEMP-TBFRZ)/(TTFRZ-TBFRZ)
753	R1=1.
754	IFLAG=1
755	ENDIF
756	QNWFRZ=QNEWLQ
757	QNEWIC=QNEWIC+QNEWLQR10.5
758	QNEWLQ=QNEWLQ-QNEWLQR10.5
759	EFFQ=(TTFRZ-TBFRZ)/(TTEMP-TBFRZ)
760	TTEMP=TU(NK1)
761	ENDIF
762	!
763	! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
764	!
765	IF(NK.EQ.K)THEN
766	BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
767	BOTERM=2.(Z0(NK1)-ZLCL)G*BE/1.5
768	ENTERM=0.
769	DZZ=Z0(NK1)-ZLCL
770	ELSE
771	BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
772	BOTERM=2.DZA(NK)G*BE/1.5
773	ENTERM=2.UER(NK)WTW/UPOLD
774	DZZ=DZA(NK)
775	ENDIF
776	WSQ=WTW
777	CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM,RATE, &
778	QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1), G)
779
780	!...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
781	! IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
782	!
783	IF(WTW.LE.0.)GOTO 65
784	WABS=SQRT(ABS(WTW))
785	WU(NK1)=WTW/WABS
786	!
787	! UPDATE THE ABE FOR UNDILUTE ASCENT...
788	!
789	THTES(NK1)=T0(NK1)(1.E5/P0(NK1))(0.2854(1.-0.28*QES(NK1))) &
790	* &
791	EXP((3374.6525/T0(NK1)-2.5403)QES(NK1)(1.+0.81* &
792	QES(NK1)))
793	UDLBE=((2.THTUDL)/(THTES(NK)+THTES(NK1))-1.)DZZ
794	IF(UDLBE.GT.0.)ABE=ABE+UDLBE*G
795	!
796	! DETERMINE THE EFFECTS OF CLOUD GLACIATION IF WITHIN THE SPECIFIED
797	! TEMP INTERVAL...
798	!
799	IF(FRC1.GT.1.E-6)THEN
800	CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QLIQ(NK1), &
801	QICE(NK1),RATIO2(NK1),TTFRZ,TBFRZ,QNWFRZ,RL,FRC1,EFFQ, &
802	IFLAG,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE &
803	,CICE,DICE)
804	ENDIF
805	!
806	! CALL SUBROUTINE TO CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMP.
807	! WITHIN GLACIATION INTERVAL, THETAE MUST BE CALCULATED WITH RESPECT TO
808	! SAME DEGREE OF GLACIATION FOR ALL ENTRAINING AIR...
809	!
810	CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),RATIO2(NK1), &
811	RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
812
813	!...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
814	!
815	REI=VMFLCLDP(NK1)0.03/RAD
816	TVQU(NK1)=TU(NK1)(1.+0.608QU(NK1)-QLIQ(NK1)-QICE(NK1))
817	!
818	!...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, NO
819	! ENTRAINMENT IS ALLOWED AT THIS LEVEL...
820	!
821	IF(TVQU(NK1).LE.TV0(NK1))THEN
822	UER(NK1)=0.0
823	UDR(NK1)=REI
824	EE2=0.
825	UD2=1.
826	EQFRC(NK1)=0.
827	GOTO 55
828	ENDIF
829	LET=NK1
830	TTMP=TVQU(NK1)
831	!
832	!...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL
833	! AIR FOR ESTIMATION OF ENTRAINMENT AND DETRAINMENT RATES...
834	!
835	F1=0.95
836	F2=1.-F1
837	THTTMP=F1THETEE(NK1)+F2THETEU(NK1)
838	QTMP=F1Q0(NK1)+F2QU(NK1)
839	TMPLIQ=F2*QLIQ(NK1)
840	TMPICE=F2*QICE(NK1)
841	CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
842	QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
843	DLIQ,AICE,BICE,CICE,DICE)
844	TU95=TTMP(1.+0.608QTMP-TMPLIQ-TMPICE)
845	IF(TU95.GT.TV0(NK1))THEN
846	EE2=1.
847	UD2=0.
848	EQFRC(NK1)=1.0
849	GOTO 50
850	ENDIF
851	F1=0.10
852	F2=1.-F1
853	THTTMP=F1THETEE(NK1)+F2THETEU(NK1)
854	QTMP=F1Q0(NK1)+F2QU(NK1)
855	TMPLIQ=F2*QLIQ(NK1)
856	TMPICE=F2*QICE(NK1)
857	CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
858	QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
859	DLIQ,AICE,BICE,CICE,DICE)
860	TU10=TTMP(1.+0.608QTMP-TMPLIQ-TMPICE)
861	IF(TU10.EQ.TVQU(NK1))THEN
862	EE2=1.
863	UD2=0.
864	EQFRC(NK1)=1.0
865	GOTO 50
866	ENDIF
867	EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
868	EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
869	EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
870	IF(EQFRC(NK1).EQ.1)THEN
871	EE2=1.
872	UD2=0.
873	GOTO 50
874	ELSEIF(EQFRC(NK1).EQ.0.)THEN
875	EE2=0.
876	UD2=1.
877	GOTO 50
878	ELSE
879	!
880	!...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
881	! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
882	!
883	CALL PROF5(EQFRC(NK1),EE2,UD2)
884	ENDIF
885	!
886	50 IF(NK.EQ.K)THEN
887	EE1=1.
888	UD1=0.
889	ENDIF
890	!
891	!...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE
892	! FRACTIONAL VALUES IN THE LAYER...
893	!
894	UER(NK1)=0.5REI(EE1+EE2)
895	UDR(NK1)=0.5REI(UD1+UD2)
896	!
897	!...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
898	! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATION
899	!
900	55 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
901	!
902	!...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL
903	! UPDRAFT FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE
904	! PREVIOUS MODEL
905	!
906	IF(UDLBE.GT.0.)ABE=ABE-UDLBE*G
907	LET=NK
908	! WRITE(98,1015)P0(NK1)/100.
909	GOTO 65
910	ENDIF
911	EE1=EE2
912	UD1=UD2
913	UPOLD=UMF(NK)-UDR(NK1)
914	UPNEW=UPOLD+UER(NK1)
915	UMF(NK1)=UPNEW
916	!
917	!...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND ICE IN
918	! THE DETRAINING UPDRAFT MASS...
919	!
920	DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
921	DETIC(NK1)=QICE(NK1)*UDR(NK1)
922	QDT(NK1)=QU(NK1)
923	QU(NK1)=(UPOLDQU(NK1)+UER(NK1)Q0(NK1))/UPNEW
924	THETEU(NK1)=(THETEU(NK1)UPOLD+THETEE(NK1)UER(NK1))/UPNEW
925	QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
926	QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
927	!
928	!...KFRZ IS THE HIGHEST MODEL LEVEL AT WHICH LIQUID CONDENSATE IS
929	! GENERATING PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF LIQUID
930	! PRECIP AT A GIVING MODEL LVL, PPTICE THE SAME FOR ICE, TRPPT IS
931	! THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE CURRENT MODEL LEVEL
932	!
933	IF(ABS(RATIO2(NK1)-1.).GT.1.E-6)KFRZ=NK1
934	PPTLIQ(NK1)=QLQOUT(NK1)*(UMF(NK)-UDR(NK1))
935	PPTICE(NK1)=QICOUT(NK1)*(UMF(NK)-UDR(NK1))
936	TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
937	IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
938	60 CONTINUE
939	!
940	!...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIU
941	! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
942	! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE
943	! BETWEEN THE LET AND CLOUD TOP...
944	!
945	!...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL
946	! VELOCITY FIRST BECOMES NEGATIVE...
947	!
948	65 LTOP=NK
949	CLDHGT=Z0(LTOP)-ZLCL
950	!
951	!...IF CLOUD TOP HGT IS LESS THAN SPECIFIED MINIMUM HEIGHT, GO BACK AND
952	! THE NEXT HIGHEST 60MB LAYER TO SEE IF A BIGGER CLOUD CAN BE OBTAINED
953	! THAT SOURCE AIR...
954	!
955	! IF(CLDHGT.LT.4.E3.OR.ABE.LT.1.)THEN
956	IF(CLDHGT.LT.3.E3.OR.ABE.LT.1.)THEN
957	DO 70 NK=K,LTOP
958	UMF(NK)=0.
959	UDR(NK)=0.
960	UER(NK)=0.
961	DETLQ(NK)=0.
962	DETIC(NK)=0.
963	PPTLIQ(NK)=0.
964	70 PPTICE(NK)=0.
965	GOTO 25
966	ENDIF
967	!
968	!...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS
969	! FLUX THIS LEVEL...
970	!
971	IF(LET.EQ.LTOP)THEN
972	UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
973	DETLQ(LTOP)=QLIQ(LTOP)UDR(LTOP)UPNEW/UPOLD
974	DETIC(LTOP)=QICE(LTOP)UDR(LTOP)UPNEW/UPOLD
975	TRPPT=TRPPT-(PPTLIQ(LTOP)+PPTICE(LTOP))
976	UER(LTOP)=0.
977	UMF(LTOP)=0.
978	GOTO 85
979	ENDIF
980	!
981	! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
982	!
983	DPTT=0.
984	DO 71 NJ=LET+1,LTOP
985	71 DPTT=DPTT+DP(NJ)
986	DUMFDP=UMF(LET)/DPTT
987	!
988	!...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
989	! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
990	! PTOP
991	!
992	DO 75 NK=LET+1,LTOP
993	UDR(NK)=DP(NK)*DUMFDP
994	UMF(NK)=UMF(NK-1)-UDR(NK)
995	DETLQ(NK)=QLIQ(NK)*UDR(NK)
996	DETIC(NK)=QICE(NK)*UDR(NK)
997	TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
998	PPTLIQ(NK)=(UMF(NK-1)-UDR(NK))*QLQOUT(NK)
999	PPTICE(NK)=(UMF(NK-1)-UDR(NK))*QICOUT(NK)
1000	TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
1001	75 CONTINUE
1002	!
1003	!...SEND UPDRAFT CHARACTERISTICS TO OUTPUT FILES...
1004	!
1005	85 CONTINUE
1006	!
1007	!...EXTEND THE UPDRAFT MASS FLUX PROFILE DOWN TO THE SOURCE LAYER FOR
1008	! THE UPDRAFT AIR...ALSO, DEFINE THETAE FOR LEVELS BELOW THE LCL...
1009	!
1010	DO 90 NK=1,K
1011	IF(NK.GE.LC)THEN
1012	IF(NK.EQ.LC)THEN
1013	UMF(NK)=VMFLCL*DP(NK)/DPTHMX
1014	UER(NK)=VMFLCL*DP(NK)/DPTHMX
1015	ELSEIF(NK.LE.KPBL)THEN
1016	UER(NK)=VMFLCL*DP(NK)/DPTHMX
1017	UMF(NK)=UMF(NK-1)+UER(NK)
1018	ELSE
1019	UMF(NK)=VMFLCL
1020	UER(NK)=0.
1021	ENDIF
1022	TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
1023	QU(NK)=QMIX
1024	WU(NK)=WLCL
1025	ELSE
1026	TU(NK)=0.
1027	QU(NK)=0.
1028	UMF(NK)=0.
1029	WU(NK)=0.
1030	UER(NK)=0.
1031	ENDIF
1032	UDR(NK)=0.
1033	QDT(NK)=0.
1034	QLIQ(NK)=0.
1035	QICE(NK)=0.
1036	QLQOUT(NK)=0.
1037	QICOUT(NK)=0.
1038	PPTLIQ(NK)=0.
1039	PPTICE(NK)=0.
1040	DETLQ(NK)=0.
1041	DETIC(NK)=0.
1042	RATIO2(NK)=0.
1043	EE=Q0(NK)*P0(NK)/(EP2+Q0(NK))
1044	TLOG=ALOG(EE/ALIQ)
1045	TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1046	TSAT=TDPT-(.212+1.571E-3(TDPT-T00)-4.36E-4(T0(NK)-T00))*( &
1047	T0(NK)-TDPT)
1048	THTA=T0(NK)(1.E5/P0(NK))(0.2854(1.-0.28*Q0(NK)))
1049	THETEE(NK)=THTA* &
1050	EXP((3374.6525/TSAT-2.5403)Q0(NK)(1.+0.81*Q0(NK)) &
1051	)
1052	THTES(NK)=THTA* &
1053	EXP((3374.6525/T0(NK)-2.5403)QES(NK)(1.+0.81* &
1054	QES(NK)))
1055	EQFRC(NK)=1.0
1056	90 CONTINUE
1057	!
1058	LTOP1=LTOP+1
1059	LTOPM1=LTOP-1
1060	!
1061	!...DEFINE VARIABLES ABOVE CLOUD TOP...
1062	!
1063	DO 95 NK=LTOP1,KX
1064	UMF(NK)=0.
1065	UDR(NK)=0.
1066	UER(NK)=0.
1067	QDT(NK)=0.
1068	QLIQ(NK)=0.
1069	QICE(NK)=0.
1070	QLQOUT(NK)=0.
1071	QICOUT(NK)=0.
1072	DETLQ(NK)=0.
1073	DETIC(NK)=0.
1074	PPTLIQ(NK)=0.
1075	PPTICE(NK)=0.
1076	IF(NK.GT.LTOP1)THEN
1077	TU(NK)=0.
1078	QU(NK)=0.
1079	WU(NK)=0.
1080	ENDIF
1081	THTA0(NK)=0.
1082	THTAU(NK)=0.
1083	EMS(NK)=DP(NK)*DXSQ/G
1084	EMSD(NK)=1./EMS(NK)
1085	TG(NK)=T0(NK)
1086	QG(NK)=Q0(NK)
1087	QLG(NK)=0.
1088	QIG(NK)=0.
1089	QRG(NK)=0.
1090	QSG(NK)=0.
1091	95 OMG(NK)=0.
1092	OMG(KL+1)=0.
1093	P150=P0(KLCL)-1.50E4
1094	DO 100 NK=1,LTOP
1095	THTAD(NK)=0.
1096	EMS(NK)=DP(NK)*DXSQ/G
1097	EMSD(NK)=1./EMS(NK)
1098	!
1099	!...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION
1100	! SCHEME
1101	!
1102	EXN(NK)=(P00/P0(NK))*(0.2854(1.-0.28*QDT(NK)))
1103	THTAU(NK)=TU(NK)*EXN(NK)
1104	EXN(NK)=(P00/P0(NK))*(0.2854(1.-0.28*Q0(NK)))
1105	THTA0(NK)=T0(NK)*EXN(NK)
1106	!
1107	!...LVF IS THE LEVEL AT WHICH MOISTURE FLUX IS ESTIMATED AS THE BASIS
1108	!...FOR PRECIPITATION EFFICIENCY CALCULATIONS...
1109	!
1110	IF(P0(NK).GT.P150)LVF=NK
1111	100 OMG(NK)=0.
1112	LVF=MIN0(LVF,LET)
1113	USR=UMF(LVF+1)*(QU(LVF+1)+QLIQ(LVF+1)+QICE(LVF+1))
1114	USR=AMIN1(USR,TRPPT)
1115	IF(USR.LT.1.E-8)USR=TRPPT
1116	!
1117	! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
1118	! * TMIX-T00,PMIX,QMIX,ABE
1119	! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
1120	! * WLCL,CLDHGT
1121	!
1122	!...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
1123	!...AND MIDTROPOSPHERE IS USED.
1124	!
1125	WSPD(KLCL)=SQRT(U0(KLCL)U0(KLCL)+V0(KLCL)V0(KLCL))
1126	WSPD(L5)=SQRT(U0(L5)U0(L5)+V0(L5)V0(L5))
1127	WSPD(LTOP)=SQRT(U0(LTOP)U0(LTOP)+V0(LTOP)V0(LTOP))
1128	VCONV=.5*(WSPD(KLCL)+WSPD(L5))
1129	if (VCONV .gt. 0.) then
1130	TIMEC=DX/VCONV
1131	else
1132	TIMEC=3600.
1133	endif
1134	! TIMEC=DX/VCONV
1135	TADVEC=TIMEC
1136	TIMEC=AMAX1(1800.,TIMEC)
1137	TIMEC=AMIN1(3600.,TIMEC)
1138	NIC=NINT(TIMEC/DT)
1139	TIMEC=FLOAT(NIC)*DT
1140	!
1141	!...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
1142	!
1143	! SHSIGN = CVMGT(1.,-1.,WSPD(LTOP).GT.WSPD(KLCL))
1144	IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
1145	SHSIGN=1.
1146	ELSE
1147	SHSIGN=-1.
1148	ENDIF
1149	VWS=(U0(LTOP)-U0(KLCL))(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL)) &
1150	(V0(LTOP)-V0(KLCL))
1151	VWS=1.E3SHSIGNSQRT(VWS)/(Z0(LTOP)-Z0(LCL))
1152	PEF=1.591+VWS(-.639+VWS(9.53E-2-VWS*4.96E-3))
1153	PEF=AMAX1(PEF,.2)
1154	PEF=AMIN1(PEF,.9)
1155	!
1156	!...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
1157	!
1158	CBH=(ZLCL-Z0(1))*3.281E-3
1159	IF(CBH.LT.3.)THEN
1160	RCBH=.02
1161	ELSE
1162	RCBH=.96729352+CBH(-.70034167+CBH(.162179896+CBH*(- &
1163	1.2569798E-2+CBH(4.2772E-4-CBH5.44E-6))))
1164	ENDIF
1165	IF(CBH.GT.25)RCBH=2.4
1166	PEFCBH=1./(1.+RCBH)
1167	PEFCBH=AMIN1(PEFCBH,.9)
1168	!
1169	!... MEAN PEF. IS USED TO COMPUTE RAINFALL.
1170	!
1171	PEFF=.5*(PEF+PEFCBH)
1172	PEFF2=PEFF
1173	! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1174	!
1175	!*****************************************************************
1176	! *
1177	! COMPUTE DOWNDRAFT PROPERTIES *
1178	! *
1179	!*****************************************************************
1180	!
1181	!...LET DOWNDRAFT ORIGINATE AT THE LEVEL OF MINIMUM SATURATION EQUIVALEN
1182	!...POTENTIAL TEMPERATURE (SEQT) IN THE CLOUD LAYER, EXTEND DOWNWARD TO
1183	!...SURFACE, OR TO THE LAYER BELOW CLOUD BASE AT WHICH ENVIR SEQT IS LES
1184	!...THAN MIN SEQT IN THE CLOUD LAYER...LET DOWNDRAFT DETRAIN OVER A LAYE
1185	!...OF SPECIFIED PRESSURE-DEPTH (DPDD)...
1186	!
1187	TDER=0.
1188	KSTART=MAX0(KPBL,KLCL)
1189	THTMIN=THTES(KSTART+1)
1190	KMIN=KSTART+1
1191	DO 104 NK=KSTART+2,LTOP-1
1192	THTMIN=AMIN1(THTMIN,THTES(NK))
1193	IF(THTMIN.EQ.THTES(NK))KMIN=NK
1194	104 CONTINUE
1195	LFS=KMIN
1196	IF(RATIO2(LFS).GT.0.)CALL ENVIRTHT(P0(LFS),T0(LFS),Q0(LFS), &
1197	THETEE(LFS),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1198	EQFRC(LFS)=(THTES(LFS)-THETEU(LFS))/(THETEE(LFS)-THETEU(LFS))
1199	EQFRC(LFS)=AMAX1(EQFRC(LFS),0.)
1200	EQFRC(LFS)=AMIN1(EQFRC(LFS),1.)
1201	THETED(LFS)=THTES(LFS)
1202	!
1203	!...ESTIMATE THE EFFECT OF MELTING PRECIPITATION IN THE DOWNDRAFT...
1204	!
1205	IF(ML.GT.0)THEN
1206	DTMLTD=0.5(QU(KLCL)-QU(LTOP))RLF/CP
1207	ELSE
1208	DTMLTD=0.
1209	ENDIF
1210	TZ(LFS)=T0(LFS)-DTMLTD
1211	ES=ALIQEXP((TZ(LFS)BLIQ-CLIQ)/(TZ(LFS)-DLIQ))
1212	QS=EP2*ES/(P0(LFS)-ES)
1213	QD(LFS)=EQFRC(LFS)Q0(LFS)+(1.-EQFRC(LFS))QU(LFS)
1214	THTAD(LFS)=TZ(LFS)(P00/P0(LFS))(0.2854(1.-0.28*QD(LFS)))
1215	IF(QD(LFS).GE.QS)THEN
1216	THETED(LFS)=THTAD(LFS)* &
1217	EXP((3374.6525/TZ(LFS)-2.5403)QS(1.+0.81*QS))
1218	ELSE
1219	CALL ENVIRTHT(P0(LFS),TZ(LFS),QD(LFS),THETED(LFS),0.,RL,EP2,ALIQ, &
1220	BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1221	ENDIF
1222	DO 107 NK=1,LFS
1223	ND=LFS-NK
1224	IF(THETED(LFS).GT.THTES(ND).OR.ND.EQ.1)THEN
1225	LDB=ND
1226	!
1227	!...IF DOWNDRAFT NEVER BECOMES NEGATIVELY BUOYANT OR IF IT
1228	!...IS SHALLOWER 50 mb, DON'T ALLOW IT TO OCCUR AT ALL...
1229	!
1230	IF(NK.EQ.1.OR.(P0(LDB)-P0(LFS)).LT.50.E2)GOTO 141
1231	! testing ---- no downdraft
1232	! GOTO 141
1233	GOTO 110
1234	ENDIF
1235	107 CONTINUE
1236	!
1237	!...ALLOW DOWNDRAFT TO DETRAIN IN A SINGLE LAYER, BUT WITH DOWNDRAFT AIR
1238	!...TYPICALLY FLUSHED UP INTO HIGHER LAYERS AS ALLOWED IN THE TOTAL
1239	!...VERTICAL ADVECTION CALCULATIONS FARTHER DOWN IN THE CODE...
1240	!
1241	110 DPDD=DP(LDB)
1242	LDT=LDB
1243	FRC=1.
1244	DPT=0.
1245	! DO 115 NK=LDB,LFS
1246	! DPT=DPT+DP(NK)
1247	! IF(DPT.GT.DPDD)THEN
1248	! LDT=NK
1249	! FRC=(DPDD+DP(NK)-DPT)/DP(NK)
1250	! GOTO 120
1251	! ENDIF
1252	! IF(NK.EQ.LFS-1)THEN
1253	! LDT=NK
1254	! FRC=1.
1255	! DPDD=DPT
1256	! GOTO 120
1257	! ENDIF
1258	!115 CONTINUE
1259	120 CONTINUE
1260	!
1261	!...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX..
1262	!
1263	TVD(LFS)=T0(LFS)(1.+0.608QES(LFS))
1264	RDD=P0(LFS)/(R*TVD(LFS))
1265	A1=(1.-PEFF)*AU0
1266	DMF(LFS)=-A1*RDD
1267	DER(LFS)=EQFRC(LFS)*DMF(LFS)
1268	DDR(LFS)=0.
1269	DO 140 ND=LFS-1,LDB,-1
1270	ND1=ND+1
1271	IF(ND.LE.LDT)THEN
1272	DER(ND)=0.
1273	DDR(ND)=-DMF(LDT+1)DP(ND)FRC/DPDD
1274	DMF(ND)=DMF(ND1)+DDR(ND)
1275	FRC=1.
1276	THETED(ND)=THETED(ND1)
1277	QD(ND)=QD(ND1)
1278	ELSE
1279	DER(ND)=DMF(LFS)0.03DP(ND)/RAD
1280	DDR(ND)=0.
1281	DMF(ND)=DMF(ND1)+DER(ND)
1282	IF(RATIO2(ND).GT.0.)CALL ENVIRTHT(P0(ND),T0(ND),Q0(ND), &
1283	THETEE(ND),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1284	THETED(ND)=(THETED(ND1)DMF(ND1)+THETEE(ND)DER(ND))/DMF(ND)
1285	QD(ND)=(QD(ND1)DMF(ND1)+Q0(ND)DER(ND))/DMF(ND)
1286	ENDIF
1287	140 CONTINUE
1288	TDER=0.
1289	!
1290	!...CALCULATION AN EVAPORATION RATE FOR GIVEN MASS FLUX...
1291	!
1292	DO 135 ND=LDB,LDT
1293	TZ(ND)= &
1294	TPDD(P0(ND),THETED(LDT),T0(ND),QS,QD(ND),1.0,XLV0,XLV1, &
1295	EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1296	ES=ALIQEXP((TZ(ND)BLIQ-CLIQ)/(TZ(ND)-DLIQ))
1297	QS=EP2*ES/(P0(ND)-ES)
1298	DSSDT=(CLIQ-BLIQDLIQ)/((TZ(ND)-DLIQ)(TZ(ND)-DLIQ))
1299	RL=XLV0-XLV1*TZ(ND)
1300	DTMP=RLQS(1.-RHBC)/(CP+RLRHBCQS*DSSDT)
1301	T1RH=TZ(ND)+DTMP
1302	ES=RHBCALIQEXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
1303	QSRH=EP2*ES/(P0(ND)-ES)
1304	!
1305	!...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
1306	!...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
1307	!
1308	IF(QSRH.LT.QD(ND))THEN
1309	QSRH=QD(ND)
1310	! T1RH=T1+(QS-QSRH)*RL/CP
1311	T1RH=TZ(ND)
1312	ENDIF
1313	TZ(ND)=T1RH
1314	QS=QSRH
1315	TDER=TDER+(QS-QD(ND))*DDR(ND)
1316	QD(ND)=QS
1317	135 THTAD(ND)=TZ(ND)(P00/P0(ND))(0.2854(1.-0.28*QD(ND)))
1318	!
1319	!...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
1320	!...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
1321	!
1322	141 IF(TDER.LT.1.)THEN
1323	! WRITE(98,3004)I,J
1324	3004 FORMAT(' ','I=',I3,2X,'J=',I3)
1325	PPTFLX=TRPPT
1326	CPR=TRPPT
1327	TDER=0.
1328	CNDTNF=0.
1329	UPDINC=1.
1330	LDB=LFS
1331	DO 117 NDK=1,LTOP
1332	DMF(NDK)=0.
1333	DER(NDK)=0.
1334	DDR(NDK)=0.
1335	THTAD(NDK)=0.
1336	WD(NDK)=0.
1337	TZ(NDK)=0.
1338	117 QD(NDK)=0.
1339	AINCM2=100.
1340	GOTO 165
1341	ENDIF
1342	!
1343	!...ADJUST DOWNDRAFT MASS FLUX SO THAT EVAPORATION RATE IN DOWNDRAFT IS
1344	!...CONSISTENT WITH PRECIPITATION EFFICIENCY RELATIONSHIP...
1345	!
1346	DEVDMF=TDER/DMF(LFS)
1347	PPR=0.
1348	PPTFLX=PEFF*USR
1349	RCED=TRPPT-PPTFLX
1350	!
1351	!...PPR IS THE TOTAL AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE
1352	!...UPDRAFT FROM CLOUD BASE TO THE LFS...UPDRAFT MASS FLUX WILL BE
1353	!...INCREASED UP TO THE LFS TO ACCOUNT FOR UPDRAFT AIR MIXING WITH
1354	!...ENVIRONMENTAL AIR TO THE UPDRAFT, SO PPR WILL INCREASE
1355	!...PROPORTIONATELY...
1356	!
1357	DO 132 NM=KLCL,LFS
1358	132 PPR=PPR+PPTLIQ(NM)+PPTICE(NM)
1359	IF(LFS.GE.KLCL)THEN
1360	DPPTDF=(1.-PEFF)PPR(1.-EQFRC(LFS))/UMF(LFS)
1361	ELSE
1362	DPPTDF=0.
1363	ENDIF
1364	!
1365	!...CNDTNF IS THE AMOUNT OF CONDENSATE TRANSFERRED ALONG WITH UPDRAFT
1366	!...MASS THE DOWNDRAFT AT THE LFS...
1367	!
1368	CNDTNF=(QLIQ(LFS)+QICE(LFS))*(1.-EQFRC(LFS))
1369	DMFLFS=RCED/(DEVDMF+DPPTDF+CNDTNF)
1370	IF(DMFLFS.GT.0.)THEN
1371	TDER=0.
1372	GOTO 141
1373	ENDIF
1374	!
1375	!...DDINC IS THE FACTOR BY WHICH TO INCREASE THE FIRST-GUESS DOWNDRAFT
1376	!...MASS FLUX TO SATISFY THE PRECIP EFFICIENCY RELATIONSHIP, UPDINC IS T
1377	!...WHICH TO INCREASE THE UPDRAFT MASS FLUX BELOW THE LFS TO ACCOUNT FOR
1378	!...TRANSFER OF MASS FROM UPDRAFT TO DOWNDRAFT...
1379	!
1380	! DDINC=DMFLFS/DMF(LFS)
1381	IF(LFS.GE.KLCL)THEN
1382	UPDINC=(UMF(LFS)-(1.-EQFRC(LFS))*DMFLFS)/UMF(LFS)
1383	!
1384	!...LIMIT UPDINC TO LESS THAN OR EQUAL TO 1.5...
1385	!
1386	IF(UPDINC.GT.1.5)THEN
1387	UPDINC=1.5
1388	DMFLFS2=UMF(LFS)*(UPDINC-1.)/(EQFRC(LFS)-1.)
1389	RCED2=DMFLFS2*(DEVDMF+DPPTDF+CNDTNF)
1390	PPTFLX=PPTFLX+(RCED-RCED2)
1391	PEFF2=PPTFLX/USR
1392	RCED=RCED2
1393	DMFLFS=DMFLFS2
1394	ENDIF
1395	ELSE
1396	UPDINC=1.
1397	ENDIF
1398	DDINC=DMFLFS/DMF(LFS)
1399	DO 149 NK=LDB,LFS
1400	DMF(NK)=DMF(NK)*DDINC
1401	DER(NK)=DER(NK)*DDINC
1402	DDR(NK)=DDR(NK)*DDINC
1403	149 CONTINUE
1404	CPR=TRPPT+PPR*(UPDINC-1.)
1405	PPTFLX=PPTFLX+PEFFPPR(UPDINC-1.)
1406	PEFF=PEFF2
1407	TDER=TDER*DDINC
1408	!
1409	!...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
1410	! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
1411	! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
1412	!
1413	DO 155 NK=LC,LFS
1414	UMF(NK)=UMF(NK)*UPDINC
1415	UDR(NK)=UDR(NK)*UPDINC
1416	UER(NK)=UER(NK)*UPDINC
1417	PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
1418	PPTICE(NK)=PPTICE(NK)*UPDINC
1419	DETLQ(NK)=DETLQ(NK)*UPDINC
1420	155 DETIC(NK)=DETIC(NK)*UPDINC
1421	!
1422	!...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
1423	!...DOWNDRAFT...
1424	!
1425	IF(LDB.GT.1)THEN
1426	DO 156 NK=1,LDB-1
1427	DMF(NK)=0.
1428	DER(NK)=0.
1429	DDR(NK)=0.
1430	WD(NK)=0.
1431	TZ(NK)=0.
1432	QD(NK)=0.
1433	THTAD(NK)=0.
1434	156 CONTINUE
1435	ENDIF
1436	DO 157 NK=LFS+1,KX
1437	DMF(NK)=0.
1438	DER(NK)=0.
1439	DDR(NK)=0.
1440	WD(NK)=0.
1441	TZ(NK)=0.
1442	QD(NK)=0.
1443	THTAD(NK)=0.
1444	157 CONTINUE
1445	DO 158 NK=LDT+1,LFS-1
1446	TZ(NK)=0.
1447	QD(NK)=0.
1448	158 CONTINUE
1449	!
1450	!
1451	!...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE
1452	! INFLOW INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN
1453	! IS AVAILABLE IN THAT LAYER INITIALLY...
1454	!
1455	165 AINCMX=1000.
1456	LMAX=MAX0(KLCL,LFS)
1457	DO 166 NK=LC,LMAX
1458	IF((UER(NK)-DER(NK)).GT.0.)AINCM1=EMS(NK)/((UER(NK)-DER(NK))* &
1459	TIMEC)
1460	AINCMX=AMIN1(AINCMX,AINCM1)
1461	166 CONTINUE
1462	AINC=1.
1463	IF(AINCMX.LT.AINC)AINC=AINCMX
1464	!
1465	!...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRFT AND DOWNDRFT...THEY
1466	!...WILL ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE
1467	!...STABILIZATION CLOSURE...
1468	!
1469	NCOUNT=0
1470	TDER2=TDER
1471	PPTFL2=PPTFLX
1472	DO 170 NK=1,LTOP
1473	DETLQ2(NK)=DETLQ(NK)
1474	DETIC2(NK)=DETIC(NK)
1475	UDR2(NK)=UDR(NK)
1476	UER2(NK)=UER(NK)
1477	DDR2(NK)=DDR(NK)
1478	DER2(NK)=DER(NK)
1479	UMF2(NK)=UMF(NK)
1480	DMF2(NK)=DMF(NK)
1481	170 CONTINUE
1482	FABE=1.
1483	STAB=0.95
1484	NOITR=0
1485	IF(AINC/AINCMX.GT.0.999)THEN
1486	NCOUNT=0
1487	GOTO 255
1488	ENDIF
1489	ISTOP=0
1490	175 NCOUNT=NCOUNT+1
1491	!
1492	!*****************************************************************
1493	! *
1494	! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
1495	! *
1496	!*****************************************************************
1497	!
1498	!...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
1499	!...SATISFY MASS CONTINUITY...
1500	!
1501	185 CONTINUE
1502	DTT=TIMEC
1503	DO 200 NK=1,LTOP
1504	DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
1505	IF(NK.GT.1)THEN
1506	OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
1507	DTT1=0.75*DP(NK-1)/(ABS(OMG(NK))+1.E-10)
1508	DTT=AMIN1(DTT,DTT1)
1509	ENDIF
1510	200 CONTINUE
1511	DO 488 NK=1,LTOP
1512	THPA(NK)=THTA0(NK)
1513	QPA(NK)=Q0(NK)
1514	NSTEP=NINT(TIMEC/DTT+1)
1515	DTIME=TIMEC/FLOAT(NSTEP)
1516	FXM(NK)=OMG(NK)*DXSQ/G
1517	488 CONTINUE
1518	!
1519	!...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
1520	!
1521	DO 495 NTC=1,NSTEP
1522	!
1523	!...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED
1524	!...SIGN OF OMEGA...
1525	!
1526	DO 493 NK=1,LTOP
1527	THFXTOP(NK)=0.
1528	THFXBOT(NK)=0.
1529	QFXTOP(NK)=0.
1530	493 QFXBOT(NK)=0.
1531	DO 494 NK=2,LTOP
1532	IF(OMG(NK).LE.0.)THEN
1533	THFXBOT(NK)=-FXM(NK)*THPA(NK-1)
1534	QFXBOT(NK)=-FXM(NK)*QPA(NK-1)
1535	THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
1536	QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
1537	ELSE
1538	THFXBOT(NK)=-FXM(NK)*THPA(NK)
1539	QFXBOT(NK)=-FXM(NK)*QPA(NK)
1540	THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
1541	QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
1542	ENDIF
1543	494 CONTINUE
1544	!
1545	!...UPDATE THE THETA AND QV VALUES AT EACH LEVEL..
1546	!
1547	DO 492 NK=1,LTOP
1548	THPA(NK)=THPA(NK)+(THFXBOT(NK)+UDR(NK)THTAU(NK)+DDR(NK) &
1549	THTAD(NK)+THFXTOP(NK)-(UER(NK)-DER(NK))THTA0(NK)) &
1550	DTIME*EMSD(NK)
1551	QPA(NK)=QPA(NK)+(QFXBOT(NK)+UDR(NK)QDT(NK)+DDR(NK)QD(NK)+ &
1552	QFXTOP(NK)-(UER(NK)-DER(NK))Q0(NK))DTIME*EMSD(NK)
1553
1554	492 CONTINUE
1555	495 CONTINUE
1556	DO 498 NK=1,LTOP
1557	THTAG(NK)=THPA(NK)
1558	QG(NK)=QPA(NK)
1559	498 CONTINUE
1560	!
1561	!...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO,
1562	!...BORROW MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO.
1563	!
1564	DO 499 NK=1,LTOP
1565	IF(QG(NK).LT.0.)THEN
1566	IF(NK.EQ.1)THEN
1567	CALL wrf_error_fatal ( 'module_cu_kf.F: problem with kf scheme: qg = 0 at the surface' )
1568	ENDIF
1569	NK1=NK+1
1570	IF(NK.EQ.LTOP)NK1=KLCL
1571	TMA=QG(NK1)*EMS(NK1)
1572	TMB=QG(NK-1)*EMS(NK-1)
1573	TMM=(QG(NK)-1.E-9)*EMS(NK)
1574	BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
1575	ACOEFF=BCOEFF*TMA/TMB
1576	TMB=TMB*(1.-BCOEFF)
1577	TMA=TMA*(1.-ACOEFF)
1578	IF(NK.EQ.LTOP)THEN
1579	QVDIFF=(QG(NK1)-TMAEMSD(NK1))100./QG(NK1)
1580	IF(ABS(QVDIFF).GT.1.)THEN
1581	PRINT *,'--WARNING-- CLOUD BASE WATER VAPOR CHANGES BY ', &
1582	QVDIFF, &
1583	' PERCENT WHEN MOISTURE IS BORROWED TO PREVENT NEG VALUES', &
1584	' IN KAIN-FRITSCH'
1585	ENDIF
1586	ENDIF
1587	QG(NK)=1.E-9
1588	QG(NK1)=TMA*EMSD(NK1)
1589	QG(NK-1)=TMB*EMSD(NK-1)
1590	ENDIF
1591	499 CONTINUE
1592	TOPOMG=(UDR(LTOP)-UER(LTOP))DP(LTOP)EMSD(LTOP)
1593	IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
1594	! WRITE(98,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;'
1595	! * ,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1596	WRITE(6,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;' &
1597	,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1598	ISTOP=1
1599	GOTO 265
1600	ENDIF
1601	!
1602	!...CONVERT THETA TO T...
1603	!
1604	! PAY ATTENTION ...
1605	!
1606	DO 230 NK=1,LTOP
1607	EXN(NK)=(P00/P0(NK))*(0.2854(1.-0.28*QG(NK)))
1608	TG(NK)=THTAG(NK)/EXN(NK)
1609	TVG(NK)=TG(NK)(1.+0.608QG(NK))
1610	230 CONTINUE
1611	!
1612	!*******************************************************************
1613	! *
1614	! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
1615	! *
1616	!*******************************************************************
1617	!
1618	!...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
1619	!
1620	THMIX=0.
1621	QMIX=0.
1622	PMIX=0.
1623	DO 217 NK=LC,KPBL
1624	ROCPQ=0.2854(1.-0.28QG(NK))
1625	THMIX=THMIX+DP(NK)TG(NK)(P00/P0(NK))**ROCPQ
1626	QMIX=QMIX+DP(NK)*QG(NK)
1627	217 PMIX=PMIX+DP(NK)*P0(NK)
1628	THMIX=THMIX/DPTHMX
1629	QMIX=QMIX/DPTHMX
1630	PMIX=PMIX/DPTHMX
1631	ROCPQ=0.2854(1.-0.28QMIX)
1632	TMIX=THMIX(PMIX/P00)*ROCPQ
1633	ES=ALIQEXP((TMIXBLIQ-CLIQ)/(TMIX-DLIQ))
1634	QS=EP2*ES/(PMIX-ES)
1635	!
1636	!...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
1637	!
1638	IF(QMIX.GT.QS)THEN
1639	RL=XLV0-XLV1*TMIX
1640	CPM=CP(1.+0.887QMIX)
1641	DSSDT=QS(CLIQ-BLIQDLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
1642	DQ=(QMIX-QS)/(1.+RL*DSSDT/CPM)
1643	TMIX=TMIX+RL/CP*DQ
1644	QMIX=QMIX-DQ
1645	ROCPQ=0.2854(1.-0.28QMIX)
1646	THMIX=TMIX(P00/PMIX)*ROCPQ
1647	TLCL=TMIX
1648	PLCL=PMIX
1649	ELSE
1650	QMIX=AMAX1(QMIX,0.)
1651	EMIX=QMIX*PMIX/(EP2+QMIX)
1652	TLOG=ALOG(EMIX/ALIQ)
1653	TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1654	TLCL=TDPT-(.212+1.571E-3(TDPT-T00)-4.36E-4(TMIX-T00))*(TMIX- &
1655	TDPT)
1656	TLCL=AMIN1(TLCL,TMIX)
1657	CPORQ=1./ROCPQ
1658	PLCL=P00(TLCL/THMIX)*CPORQ
1659	ENDIF
1660	TVLCL=TLCL(1.+0.608QMIX)
1661	DO 235 NK=LC,KL
1662	KLCL=NK
1663	235 IF(PLCL.GE.P0(NK))GOTO 240
1664	240 K=KLCL-1
1665	DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
1666	!
1667	!...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
1668	!
1669	TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
1670	QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
1671	TVEN=TENV(1.+0.608QENV)
1672	TVBAR=0.5*(TVG(K)+TVEN)
1673	! ZLCL=Z0(K)+RTVBARALOG(P0(K)/PLCL)/G
1674	ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
1675	TVAVG=0.5(TVEN+TG(KLCL)(1.+0.608*QG(KLCL)))
1676	PLCL=P0(KLCL)EXP(G/(RTVAVG)*(Z0(KLCL)-ZLCL))
1677	THETEU(K)=TMIX(1.E5/PMIX)(0.2854(1.-0.28QMIX)) &
1678	EXP((3374.6525/TLCL-2.5403)QMIX(1.+0.81*QMIX))
1679	ES=ALIQEXP((TENVBLIQ-CLIQ)/(TENV-DLIQ))
1680	QESE=EP2*ES/(PLCL-ES)
1681	THTESG(K)=TENV(1.E5/PLCL)(0.2854(1.-0.28QESE)) &
1682	EXP((3374.6525/TENV-2.5403)QESE(1.+0.81*QESE))
1683	!
1684	!...COMPUTE ADJUSTED ABE(ABEG).
1685	!
1686	ABEG=0.
1687	THTUDL=THETEU(K)
1688	DO 245 NK=K,LTOPM1
1689	NK1=NK+1
1690	ES=ALIQEXP((TG(NK1)BLIQ-CLIQ)/(TG(NK1)-DLIQ))
1691	QESE=EP2*ES/(P0(NK1)-ES)
1692	THTESG(NK1)=TG(NK1)(1.E5/P0(NK1))(0.2854(1.-0.28QESE)) &
1693	EXP((3374.6525/TG(NK1)-2.5403)QESE(1.+0.81*QESE) &
1694	)
1695	! DZZ=CVMGT(Z0(KLCL)-ZLCL,DZA(NK),NK.EQ.K)
1696	IF(NK.EQ.K)THEN
1697	DZZ=Z0(KLCL)-ZLCL
1698	ELSE
1699	DZZ=DZA(NK)
1700	ENDIF
1701	BE=((2.THTUDL)/(THTESG(NK1)+THTESG(NK))-1.)DZZ
1702	245 IF(BE.GT.0.)ABEG=ABEG+BE*G
1703	!
1704	!...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
1705	!...THE PERIOD TIMEC...
1706	!
1707	IF(NOITR.EQ.1)THEN
1708	! WRITE(98,1060)FABE
1709	GOTO 265
1710	ENDIF
1711	DABE=AMAX1(ABE-ABEG,0.1*ABE)
1712	FABE=ABEG/(ABE+1.E-8)
1713	IF(FABE.GT.1.)THEN
1714	! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS '
1715	! *,'GRID POINT; NO CONVECTION ALLOWED!'
1716	GOTO 325
1717	ENDIF
1718	IF(NCOUNT.NE.1)THEN
1719	DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
1720	IF(DFDA.GT.0.)THEN
1721	NOITR=1
1722	AINC=AINCOLD
1723	GOTO 255
1724	ENDIF
1725	ENDIF
1726	AINCOLD=AINC
1727	FABEOLD=FABE
1728	IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
1729	! WRITE(98,1055)FABE
1730	GOTO 265
1731	ENDIF
1732	IF(FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB)GOTO 265
1733	IF(NCOUNT.GT.10)THEN
1734	! WRITE(98,1060)FABE
1735	GOTO 265
1736	ENDIF
1737	!
1738	!...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE
1739	!...CONVECTIVE MASS FLUX BY THE FACTOR AINC:
1740	!
1741	IF(FABE.EQ.0.)THEN
1742	AINC=AINC*0.5
1743	ELSE
1744	AINC=AINCSTABABE/(DABE+1.E-8)
1745	ENDIF
1746	255 AINC=AMIN1(AINCMX,AINC)
1747	!...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION
1748	!...WILL BE MINIMAL SO JUST IGNORE IT...
1749	IF(AINC.LT.0.05)GOTO 325
1750	! AINC=AMAX1(AINC,0.05)
1751	TDER=TDER2*AINC
1752	PPTFLX=PPTFL2*AINC
1753	! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABEOLD,AINCOLD
1754	DO 260 NK=1,LTOP
1755	UMF(NK)=UMF2(NK)*AINC
1756	DMF(NK)=DMF2(NK)*AINC
1757	DETLQ(NK)=DETLQ2(NK)*AINC
1758	DETIC(NK)=DETIC2(NK)*AINC
1759	UDR(NK)=UDR2(NK)*AINC
1760	UER(NK)=UER2(NK)*AINC
1761	DER(NK)=DER2(NK)*AINC
1762	DDR(NK)=DDR2(NK)*AINC
1763	260 CONTINUE
1764	!
1765	!...GO BACK UP FOR ANOTHER ITERATION...
1766	!
1767	GOTO 175
1768	265 CONTINUE
1769	!
1770	!...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
1771	!...GRID POINT...
1772	!
1773	!...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
1774	!
1775	!...FRC2 IS THE FRACTION OF TOTAL CONDENSATE
1776	!...GENERATED THAT GOES INTO PRECIPITIATION
1777	FRC2=PPTFLX/(CPR*AINC)
1778	DO 270 NK=1,LTOP
1779	QLPA(NK)=QL0(NK)
1780	QIPA(NK)=QI0(NK)
1781	QRPA(NK)=QR0(NK)
1782	QSPA(NK)=QS0(NK)
1783	RAINFB(NK)=PPTLIQ(NK)AINCFBFRC*FRC2
1784	SNOWFB(NK)=PPTICE(NK)AINCFBFRC*FRC2
1785	270 CONTINUE
1786	DO 290 NTC=1,NSTEP
1787	!
1788	!...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH
1789	!...LAYER BASED ON THE SIGN OF OMEGA...
1790	!
1791	DO 275 NK=1,LTOP
1792	QLFXIN(NK)=0.
1793	QLFXOUT(NK)=0.
1794	QIFXIN(NK)=0.
1795	QIFXOUT(NK)=0.
1796	QRFXIN(NK)=0.
1797	QRFXOUT(NK)=0.
1798	QSFXIN(NK)=0.
1799	QSFXOUT(NK)=0.
1800	275 CONTINUE
1801	DO 280 NK=2,LTOP
1802	IF(OMG(NK).LE.0.)THEN
1803	QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
1804	QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
1805	QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
1806	QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
1807	QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
1808	QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
1809	QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
1810	QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
1811	ELSE
1812	QLFXOUT(NK)=FXM(NK)*QLPA(NK)
1813	QIFXOUT(NK)=FXM(NK)*QIPA(NK)
1814	QRFXOUT(NK)=FXM(NK)*QRPA(NK)
1815	QSFXOUT(NK)=FXM(NK)*QSPA(NK)
1816	QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
1817	QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
1818	QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
1819	QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
1820	ENDIF
1821	280 CONTINUE
1822	!
1823	!...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
1824	!
1825	DO 285 NK=1,LTOP
1826	QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))DTIME &
1827	EMSD(NK)
1828	QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))DTIME &
1829	EMSD(NK)
1830	QRPA(NK)=QRPA(NK)+(QRFXIN(NK)+QLQOUT(NK)*UDR(NK)-QRFXOUT(NK) &
1831	+RAINFB(NK))DTIMEEMSD(NK)
1832	QSPA(NK)=QSPA(NK)+(QSFXIN(NK)+QICOUT(NK)*UDR(NK)-QSFXOUT(NK) &
1833	+SNOWFB(NK))DTIMEEMSD(NK)
1834	285 CONTINUE
1835	290 CONTINUE
1836	DO 295 NK=1,LTOP
1837	QLG(NK)=QLPA(NK)
1838	QIG(NK)=QIPA(NK)
1839	QRG(NK)=QRPA(NK)
1840	QSG(NK)=QSPA(NK)
1841	295 CONTINUE
1842	! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABE,AINC
1843	!
1844	!...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
1845	!
1846	IF(ISTOP.EQ.1)THEN
1847	WRITE(6,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
1848	' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
1849	,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
1850	DO 300 K=LTOP,1,-1
1851	DTT=(TG(K)-T0(K))*86400./TIMEC
1852	RL=XLV0-XLV1*TG(K)
1853	DR=-(QG(K)-Q0(K))RL86400./(TIMEC*CP)
1854	UDFRC=UDR(K)TIMECEMSD(K)
1855	UEFRC=UER(K)TIMECEMSD(K)
1856	DDFRC=DDR(K)TIMECEMSD(K)
1857	DEFRC=-DER(K)TIMECEMSD(K)
1858	WRITE (6,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)* &
1859	1.E4,UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC &
1860	,DDFRC,EMS(K)/1.E11,W0AVG1D(K)*1.E2,DETLQ(K) &
1861	TIMECEMSD(K)1.E3,DETIC(K)TIMECEMSD(K) &
1862	1.E3
1863	300 CONTINUE
1864	WRITE(6,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
1865	'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
1866	DO 305 K=KX,1,-1
1867	DTT=TG(K)-T0(K)
1868	TUC=TU(K)-T00
1869	IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
1870	TDC=TZ(K)-T00
1871	IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
1872	ES=ALIQEXP((BLIQTG(K)-CLIQ)/(TG(K)-DLIQ))
1873	QGS=ES*EP2/(P0(K)-ES)
1874	RH0=Q0(K)/QES(K)
1875	RHG=QG(K)/QGS
1876	WRITE (6,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC &
1877	,TDC,Q0(K)1000.,QG(K)1000.,(QG(K)-Q0(K))* &
1878	1000.,QU(K)1000.,QD(K)1000.,QLG(K)*1000., &
1879	QIG(K)1000.,QRG(K)1000.,QSG(K)*1000.,RH0,RHG
1880	305 CONTINUE
1881	!
1882	!...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
1883	!...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
1884	!
1885	IF(ISTOP.EQ.1)THEN
1886	DO 310 K=1,KX
1887	WRITE ( wrf_err_message , 1115 ) &
1888	Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
1889	U0(K),V0(K),DP(K)/100.,W0AVG1D(K)
1890	CALL wrf_message ( TRIM( wrf_err_message ) )
1891	310 CONTINUE
1892	CALL wrf_error_fatal ( 'module_cu_kf.F: KAIN-FRITSCH' )
1893	ENDIF
1894	ENDIF
1895	CNDTNF=(1.-EQFRC(LFS))(QLIQ(LFS)+QICE(LFS))DMF(LFS)
1896	! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
1897	!
1898	! EVALUATE MOISTURE BUDGET...
1899	!
1900	QINIT=0.
1901	QFNL=0.
1902	DPT=0.
1903	DO 315 NK=1,LTOP
1904	DPT=DPT+DP(NK)
1905	QINIT=QINIT+Q0(NK)*EMS(NK)
1906	QFNL=QFNL+QG(NK)*EMS(NK)
1907	QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
1908	315 CONTINUE
1909	QFNL=QFNL+PPTFLXTIMEC(1.-FBFRC)
1910	ERR2=(QFNL-QINIT)*100./QINIT
1911	! WRITE(98,1110)QINIT,QFNL,ERR2
1912	! IF(ABS(ERR2).GT.0.05)STOP 'QVERR'
1913	IF(ABS(ERR2).GT.0.05)CALL wrf_error_fatal( 'module_cu_kf.F: QVERR' )
1914	RELERR=ERR2QINIT/(PPTFLXTIMEC+1.E-10)
1915	! WRITE(98,1120)RELERR
1916	! WRITE(98,*)'TDER, CPR, USR, TRPPT =',
1917	! TDER,CPRAINC,USRAINC,TRPPTAINC
1918	!
1919	!...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
1920	!
1921	!...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
1922	!...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
1923	!
1924	IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
1925	NCA(I,J)=FLOAT(NIC)*DT
1926	DO 320 K=1,KX
1927	! IF(IMOIST.NE.2)THEN
1928	!
1929	!...IF HYDROMETEORS ARE NOT ALLOWED, THEY MUST BE EVAPORATED OR SUBLIMAT
1930	!...AND FED BACK AS VAPOR, ALONG WITH ASSOCIATED CHANGES IN TEMPERATURE.
1931	!...NOTE: THIS WILL INTRODUCE CHANGES IN THE CONVECTIVE TEMPERATURE AND
1932	!...WATER VAPOR FEEDBACK TENDENCIES AND MAY LEAD TO SUPERSATURATED VALUE
1933	!...OF QG...
1934	!
1935	! RLC=XLV0-XLV1*TG(K)
1936	! RLS=XLS0-XLS1*TG(K)
1937	! CPM=CP(1.+0.887QG(K))
1938	! TG(K)=TG(K)-(RLC(QLG(K)+QRG(K))+RLS(QIG(K)+QSG(K)))/CPM
1939	! QG(K)=QG(K)+(QLG(K)+QRG(K)+QIG(K)+QSG(K))
1940	! DQCDT(K)=0.
1941	! DQIDT(K)=0.
1942	! DQRDT(K)=0.
1943	! DQSDT(K)=0.
1944	! ELSE
1945	IF(.NOT. qi_flag .and. warm_rain)THEN
1946	!
1947	!...IF ICE PHASE IS NOT ALLOWED, MELT ALL FROZEN HYDROMETEORS...
1948	!
1949	CPM=CP(1.+0.887QG(K))
1950	TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
1951	DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
1952	DQIDT(K)=0.
1953	DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
1954	DQSDT(K)=0.
1955	ELSEIF(.NOT. qi_flag .and. .not. warm_rain)THEN
1956	!
1957	!...IF ICE PHASE IS ALLOWED, BUT MIXED PHASE IS NOT, MELT FROZEN HYDROME
1958	!...BELOW THE MELTING LEVEL, FREEZE LIQUID WATER ABOVE THE MELTING LEVEL
1959	!
1960	CPM=CP(1.+0.887QG(K))
1961	IF(K.LE.ML)THEN
1962	TG(K)=TG(K)-(QIG(K)+QSG(K))*RLF/CPM
1963	ELSEIF(K.GT.ML)THEN
1964	TG(K)=TG(K)+(QLG(K)+QRG(K))*RLF/CPM
1965	ENDIF
1966	DQCDT(K)=(QLG(K)+QIG(K)-QL0(K)-QI0(K))/TIMEC
1967	DQIDT(K)=0.
1968	DQRDT(K)=(QRG(K)+QSG(K)-QR0(K)-QS0(K))/TIMEC
1969	DQSDT(K)=0.
1970	ELSEIF(qi_flag) THEN
1971	!
1972	!...IF MIXED PHASE HYDROMETEORS ARE ALLOWED, FEED BACK CONVECTIVE
1973	!...TENDENCY OF HYDROMETEORS DIRECTLY...
1974	!
1975	DQCDT(K)=(QLG(K)-QL0(K))/TIMEC
1976	DQIDT(K)=(QIG(K)-QI0(K))/TIMEC
1977	DQRDT(K)=(QRG(K)-QR0(K))/TIMEC
1978	IF (qs_flag ) THEN
1979	DQSDT(K)=(QSG(K)-QS0(K))/TIMEC
1980	ELSE
1981	DQIDT(K)=DQIDT(K)+(QSG(K)-QS0(K))/TIMEC
1982	ENDIF
1983	ELSE
1984	CALL wrf_error_fatal ( 'module_cu_kf: THIS COMBINATION OF IMOIST, IICE NOT ALLOWED' )
1985	ENDIF
1986	! ENDIF
1987	DTDT(K)=(TG(K)-T0(K))/TIMEC
1988	DQDT(K)=(QG(K)-Q0(K))/TIMEC
1989	320 CONTINUE
1990
1991	! RAINCV is in the unit of mm
1992
1993	PRATEC(I,J)=PPTFLX*(1.-FBFRC)/DXSQ
1994	RAINCV(I,J)=DT*PRATEC(I,J)
1995	RNC=RAINCV(I,J)*NIC
1996	! WRITE(98,909)RNC
1997	909 FORMAT(' CONVECTIVE RAINFALL =',F8.4,' CM')
1998
1999	325 CONTINUE
2000
2001	1000 FORMAT(' ',10A8)
2002	1005 FORMAT(' ',F6.0,2X,F6.4,2X,F7.3,1X,F6.4,2X,4(F6.3,2X),2(F7.3,1X))
2003	1010 FORMAT(' ',' VERTICAL VELOCITY IS NEGATIVE AT ',F4.0,' MB')
2004	1015 FORMAT(' ','ALL REMAINING MASS DETRAINS BELOW ',F4.0,' MB')
2005	1025 FORMAT(5X,' KLCL=',I2,' ZLCL=',F7.1,'M', &
2006	' DTLCL=',F5.2,' LTOP=',I2,' P0(LTOP)=',-2PF5.1,'MB FRZ LV=', &
2007	I2,' TMIX=',0PF4.1,1X,'PMIX=',-2PF6.1,' QMIX=',3PF5.1, &
2008	' CAPE=',0PF7.1)
2009	1030 FORMAT(' ',' P0(LET) = ',F6.1,' P0(LTOP) = ',F6.1,' VMFLCL =', &
2010	E12.3,' PLCL =',F6.1,' WLCL =',F6.3,' CLDHGT =', &
2011	F8.1)
2012	1035 FORMAT(1X,'PEF(WS)=',F4.2,'(CB)=',F4.2,'LC,LET=',2I3,'WKL=' &
2013	,F6.3,'VWS=',F5.2)
2014	1040 FORMAT(' ','PRECIP EFF = 100%, ENVIR CANNOT SUPPORT DOWND' &
2015	,'RAFTS')
2016	!1045 FORMAT('NUMBER OF DOWNDRAFT ITERATIONS EXCEEDS 10...PPTFLX' &
2017	! ' IS DIFFERENT FROM THAT GIVEN BY PRECIP EFF RELATION')
2018	! FLIC HAS TROUBLE WITH THIS ONE.
2019	1045 FORMAT('NUMBER OF DOWNDRAFT ITERATIONS EXCEEDS 10')
2020	1050 FORMAT(' ','LCOUNT= ',I3,' PPTFLX/CPR, PEFF= ',F5.3,1X,F5.3, &
2021	'DMF(LFS)/UMF(LCL)= ',F5.3)
2022	1055 FORMAT(/'*** DEGREE OF STABILIZATION =',F5.3,', NO MORE MASS F' &
2023	,'LUX IS ALLOWED')
2024	!1060 FORMAT(/' ITERATION DOES NOT CONVERGE TO GIVE THE SPECIFIED ' &
2025	! 'DEGREE OF STABILIZATION! FABE= ',F6.4)
2026	1060 FORMAT(/' ITERATION DOES NOT CONVERGE. FABE= ',F6.4)
2027	1070 FORMAT (16A8)
2028	1075 FORMAT (F8.2,3(F8.2),2(F8.3),F8.2,2F8.3,F8.2,6F8.3)
2029	1080 FORMAT(2X,'LFS,LDB,LDT =',3I3,' TIMEC, NSTEP=',F5.0,I3, &
2030	'NCOUNT, FABE, AINC=',I2,1X,F5.3,F6.2)
2031	1085 FORMAT (A3,16A7,2A8)
2032	1090 FORMAT (I3,F7.2,F7.0,10F7.2,4F7.3,2F8.3)
2033	1095 FORMAT(' ',' PPT PRODUCTION RATE= ',F10.0,' TOTAL EVAP+PPT= ', &
2034	F10.0)
2035	1105 FORMAT(' ','NET LATENT HEAT RELEASE =',E12.5,' ACTUAL HEATING =', &
2036	E12.5,' J/KG-S, DIFFERENCE = ',F9.3,'PERCENT')
2037	1110 FORMAT(' ','INITIAL WATER =',E12.5,' FINAL WATER =',E12.5, &
2038	' TOTAL WATER CHANGE =',F8.2,'PERCENT')
2039	1115 FORMAT (2X,F6.0,2X,F7.2,2X,F5.1,2X,F6.3,2(2X,F5.1),2X,F7.2,2X,F7.4 &
2040	)
2041	1120 FORMAT(' ','MOISTURE ERROR AS FUNCTION OF TOTAL PPT =',F9.3, &
2042	'PERCENT')
2043
2044	END SUBROUTINE KFPARA
2045
2046	!-----------------------------------------------------------------------
2047	SUBROUTINE CONDLOAD(QLIQ,QICE,WTW,DZ,BOTERM,ENTERM,RATE,QNEWLQ, &
2048	QNEWIC,QLQOUT,QICOUT,G)
2049	!-----------------------------------------------------------------------
2050	IMPLICIT NONE
2051	!-----------------------------------------------------------------------
2052	! 9/18/88...THIS PRECIPITATION FALLOUT SCHEME IS BASED ON THE SCHEME US
2053	! BY OGURA AND CHO (1973). LIQUID WATER FALLOUT FROM A PARCEL IS CAL-
2054	! CULATED USING THE EQUATION DQ=-RATEQDT, BUT TO SIMULATE A QUASI-
2055	! CONTINUOUS PROCESS, AND TO ELIMINATE A DEPENDENCY ON VERTICAL
2056	! RESOLUTION THIS IS EXPRESSED AS Q=QEXP(-RATEDZ).
2057
2058	REAL, INTENT(IN ) :: G
2059	REAL, INTENT(IN ) :: DZ,BOTERM,ENTERM,RATE
2060	REAL, INTENT(INOUT) :: QLQOUT,QICOUT,WTW,QLIQ,QICE,QNEWLQ,QNEWIC
2061
2062	REAL :: QTOT,QNEW,QEST,G1,WAVG,CONV,RATIO3,OLDQ,RATIO4,DQ,PPTDRG
2063
2064	QTOT=QLIQ+QICE
2065	QNEW=QNEWLQ+QNEWIC
2066	!
2067	! ESTIMATE THE VERTICAL VELOCITY SO THAT AN AVERAGE VERTICAL VELOCITY C
2068	! BE CALCULATED TO ESTIMATE THE TIME REQUIRED FOR ASCENT BETWEEN MODEL
2069	! LEVELS...
2070	!
2071	QEST=0.5*(QTOT+QNEW)
2072	G1=WTW+BOTERM-ENTERM-2.GDZ*QEST/1.5
2073	IF(G1.LT.0.0)G1=0.
2074	WAVG=(SQRT(WTW)+SQRT(G1))/2.
2075	CONV=RATE*DZ/WAVG
2076	!
2077	! RATIO3 IS THE FRACTION OF LIQUID WATER IN FRESH CONDENSATE, RATIO4 IS
2078	! THE FRACTION OF LIQUID WATER IN THE TOTAL AMOUNT OF CONDENSATE INVOLV
2079	! IN THE PRECIPITATION PROCESS - NOTE THAT ONLY 60% OF THE FRESH CONDEN
2080	! SATE IS IS ALLOWED TO PARTICIPATE IN THE CONVERSION PROCESS...
2081	!
2082	RATIO3=QNEWLQ/(QNEW+1.E-10)
2083	! OLDQ=QTOT
2084	QTOT=QTOT+0.6*QNEW
2085	OLDQ=QTOT
2086	RATIO4=(0.6*QNEWLQ+QLIQ)/(QTOT+1.E-10)
2087	QTOT=QTOT*EXP(-CONV)
2088	!
2089	! DETERMINE THE AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE UPDRAFT
2090	! PARCEL AT THIS LEVEL...
2091	!
2092	DQ=OLDQ-QTOT
2093	QLQOUT=RATIO4*DQ
2094	QICOUT=(1.-RATIO4)*DQ
2095	!
2096	! ESTIMATE THE MEAN LOAD OF CONDENSATE ON THE UPDRAFT IN THE LAYER, CAL
2097	! LATE VERTICAL VELOCITY
2098	!
2099	PPTDRG=0.5(OLDQ+QTOT-0.2QNEW)
2100	WTW=WTW+BOTERM-ENTERM-2.GDZ*PPTDRG/1.5
2101	!
2102	! DETERMINE THE NEW LIQUID WATER AND ICE CONCENTRATIONS INCLUDING LOSSE
2103	! DUE TO PRECIPITATION AND GAINS FROM CONDENSATION...
2104	!
2105	QLIQ=RATIO4QTOT+RATIO30.4*QNEW
2106	QICE=(1.-RATIO4)QTOT+(1.-RATIO3)0.4*QNEW
2107	QNEWLQ=0.
2108	QNEWIC=0.
2109
2110	END SUBROUTINE CONDLOAD
2111
2112	!-----------------------------------------------------------------------
2113	SUBROUTINE DTFRZNEW(TU,P,THTEU,QVAP,QLIQ,QICE,RATIO2,TTFRZ,TBFRZ, &
2114	QNWFRZ,RL,FRC1,EFFQ,IFLAG,XLV0,XLV1,XLS0,XLS1, &
2115	EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2116	!-----------------------------------------------------------------------
2117	IMPLICIT NONE
2118	!-----------------------------------------------------------------------
2119	REAL, INTENT(IN ) :: XLV0,XLV1
2120	REAL, INTENT(IN ) :: P,TTFRZ,TBFRZ,EFFQ,XLS0,XLS1,EP2,ALIQ, &
2121	BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE
2122	REAL, INTENT(INOUT) :: TU,THTEU,QVAP,QLIQ,QICE,RATIO2, &
2123	FRC1,RL,QNWFRZ
2124	INTEGER, INTENT(INOUT) :: IFLAG
2125
2126	REAL :: CCP,RV,C5,QLQFRZ,QNEW,ESLIQ,ESICE,RLC,RLS,PI,ES,RLF,A, &
2127	B,C,DQVAP,DTFRZ,TU1,QVAP1
2128	!-----------------------------------------------------------------------
2129	!
2130	!...ALLOW GLACIATION OF THE UPDRAFT TO OCCUR AS AN APPROXIMATELY LINEAR
2131	! FUNCTION OF TEMPERATURE IN THE TEMPERATURE RANGE TTFRZ TO TBFRZ...
2132	!
2133
2134	RV=461.5
2135	C5=1.0723E-3
2136	!
2137	!...ADJUST THE LIQUID WATER CONCENTRATIONS FROM FRESH CONDENSATE AND THA
2138	! BROUGHT UP FROM LOWER LEVELS TO AN AMOUNT THAT WOULD BE PRESENT IF N
2139	! LIQUID WATER HAD FROZEN THUS FAR...THIS IS NECESSARY BECAUSE THE
2140	! EXPRESSION FOR TEMP CHANGE IS MULTIPLIED BY THE FRACTION EQUAL TO TH
2141	! PARCEL TEMP DECREASE SINCE THE LAST MODEL LEVEL DIVIDED BY THE TOTAL
2142	! GLACIATION INTERVAL, SO THAT EFFECTIVELY THIS APPROXIMATELY ALLOWS A
2143	! AMOUNT OF LIQUID WATER TO FREEZE WHICH IS EQUAL TO THIS SAME FRACTIO
2144	! OF THE LIQUID WATER THAT WAS PRESENT BEFORE THE GLACIATION PROCESS W
2145	! INITIATED...ALSO, TO ALLOW THETAU TO CONVERT APPROXIMATELY LINEARLY
2146	! ITS VALUE WITH RESPECT TO ICE, WE NEED TO ALLOW A PORTION OF THE FRE
2147	! CONDENSATE TO CONTRIBUTE TO THE GLACIATION PROCESS; THE FRACTIONAL
2148	! AMOUNT THAT APPLIES TO THIS PORTION IS 1/2 OF THE FRACTIONAL AMOUNT
2149	! FROZEN OF THE "OLD" CONDENSATE BECAUSE THIS FRESH CONDENSATE IS ONLY
2150	! PRODUCED GRADUALLY OVER THE LAYER...NOTE THAT IN TERMS OF THE DYNAMI
2151	! OF THE PRECIPITATION PROCESS, IE. PRECIPITATION FALLOUT, THIS FRACTI
2152	! AMNT OF FRESH CONDENSATE HAS ALREADY BEEN INCLUDED IN THE ICE CATEGO
2153	!
2154	QLQFRZ=QLIQ*EFFQ
2155	QNEW=QNWFRZEFFQ0.5
2156	ESLIQ=ALIQEXP((BLIQTU-CLIQ)/(TU-DLIQ))
2157	ESICE=AICEEXP((BICETU-CICE)/(TU-DICE))
2158	RLC=2.5E6-2369.276*(TU-273.16)
2159	RLS=2833922.-259.532*(TU-273.16)
2160	RLF=RLS-RLC
2161	CCP=1005.7(1.+0.89QVAP)
2162	!
2163	! A = D(ES)/DT IS THAT CALCULATED FROM BUCK`S (1981) EMPIRICAL FORMULAS
2164	! FOR SATURATION VAPOR PRESSURE...
2165	!
2166	A=(CICE-BICEDICE)/((TU-DICE)(TU-DICE))
2167	B=RLS*EP2/P
2168	C=ABESICE/CCP
2169	DQVAP=B(ESLIQ-ESICE)/(RLS+RLSC)-RLF*(QLQFRZ+QNEW)/(RLS+RLS/C)
2170	DTFRZ=(RLF(QLQFRZ+QNEW)+B(ESLIQ-ESICE))/(CCP+ABESICE)
2171	TU1=TU
2172	QVAP1=QVAP
2173	TU=TU+FRC1*DTFRZ
2174	QVAP=QVAP-FRC1*DQVAP
2175	ES=QVAP*P/(EP2+QVAP)
2176	ESLIQ=ALIQEXP((BLIQTU-CLIQ)/(TU-DLIQ))
2177	ESICE=AICEEXP((BICETU-CICE)/(TU-DICE))
2178	RATIO2=(ESLIQ-ES)/(ESLIQ-ESICE)
2179	!
2180	! TYPICALLY, RATIO2 IS VERY CLOSE TO (TTFRZ-TU)/(TTFRZ-TBFRZ), USUALLY
2181	! WITHIN 1% (USING TU BEFORE GALCIATION EFFECTS ARE APPLIED); IF THE
2182	! INITIAL UPDRAFT TEMP IS BELOW TBFRZ AND RATIO2 IS STILL LESS THAN 1,
2183	! AN ADJUSTMENT TO FRC1 AND RATIO2 IS INTRODUCED SO THAT GLACIATION
2184	! EFFECTS ARE NOT UNDERESTIMATED; CONVERSELY, IF RATIO2 IS GREATER THAN
2185	! FRC1 IS ADJUSTED SO THAT GLACIATION EFFECTS ARE NOT OVERESTIMATED...
2186	!
2187	IF(IFLAG.GT.0.AND.RATIO2.LT.1)THEN
2188	FRC1=FRC1+(1.-RATIO2)
2189	TU=TU1+FRC1*DTFRZ
2190	QVAP=QVAP1-FRC1*DQVAP
2191	RATIO2=1.
2192	IFLAG=1
2193	GOTO 20
2194	ENDIF
2195	IF(RATIO2.GT.1.)THEN
2196	FRC1=FRC1-(RATIO2-1.)
2197	FRC1=AMAX1(0.0,FRC1)
2198	TU=TU1+FRC1*DTFRZ
2199	QVAP=QVAP1-FRC1*DQVAP
2200	RATIO2=1.
2201	IFLAG=1
2202	ENDIF
2203	!
2204	! CALCULATE A HYBRID VALUE OF THETAU, ASSUMING THAT THE LATENT HEAT OF
2205	! VAPORIZATION/SUBLIMATION CAN BE ESTIMATED USING THE SAME WEIGHTING
2206	! FUNCTION AS THAT USED TO CALCULATE SATURATION VAPOR PRESSURE, CALCU-
2207	! LATE NEW LIQUID WATER AND ICE CONCENTRATIONS...
2208	!
2209	20 RLC=XLV0-XLV1*TU
2210	RLS=XLS0-XLS1*TU
2211	RL=RATIO2RLS+(1.-RATIO2)RLC
2212	PI=(1.E5/P)*(0.2854(1.-0.28*QVAP))
2213	THTEU=TUPIEXP(RLQVAPC5/TU(1.+0.81QVAP))
2214	IF(IFLAG.EQ.1)THEN
2215	QICE=QICE+FRC1*DQVAP+QLIQ
2216	QLIQ=0.
2217	ELSE
2218	QICE=QICE+FRC1*(DQVAP+QLQFRZ)
2219	QLIQ=QLIQ-FRC1*QLQFRZ
2220	ENDIF
2221	QNWFRZ=0.
2222
2223	END SUBROUTINE DTFRZNEW
2224
2225	!-----------------------------------------------------------------------
2226	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2227	! THIS SUBROUTINE INTEGRATES THE AREA UNDER THE CURVE IN THE GAUSSIAN
2228	! DISTRIBUTION...THE NUMERICAL APPROXIMATION TO THE INTEGRAL IS TAKEN F
2229	! HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMA
2230	! TABLES ED. BY ABRAMOWITZ AND STEGUN, NAT L BUREAU OF STANDARDS APPLI
2231	! MATHEMATICS SERIES. JUNE, 1964., MAY, 1968.
2232	! JACK KAIN
2233	! 7/6/89
2234	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2235	!***********************************************************************
2236	!*** GAUSSIAN TYPE MIXING PROFILE....****************************
2237	SUBROUTINE PROF5(EQ,EE,UD)
2238	!-----------------------------------------------------------------------
2239	IMPLICIT NONE
2240	!-----------------------------------------------------------------------
2241	REAL, INTENT(IN ) :: EQ
2242	REAL, INTENT(INOUT) :: EE,UD
2243	REAL :: SQRT2P,A1,A2,A3,P,SIGMA,FE,X,Y,EY,E45,T1,T2,C1,C2
2244
2245	DATA SQRT2P,A1,A2,A3,P,SIGMA,FE/2.506628,0.4361836,-0.1201676, &
2246	0.9372980,0.33267,0.166666667,0.202765151/
2247	X=(EQ-0.5)/SIGMA
2248	Y=6.*EQ-3.
2249	EY=EXP(Y*Y/(-2))
2250	E45=EXP(-4.5)
2251	T2=1./(1.+P*ABS(Y))
2252	T1=0.500498
2253	C1=A1T1+A2T1T1+A3T1T1T1
2254	C2=A1T2+A2T2T2+A3T2T2T2
2255	IF(Y.GE.0.)THEN
2256	EE=SIGMA(0.5(SQRT2P-E45C1-EYC2)+SIGMA(E45-EY))-E45EQ*EQ/2.
2257	UD=SIGMA(0.5(EYC2-E45C1)+SIGMA(E45-EY))-E45(0.5+EQ*EQ/2.- &
2258	EQ)
2259	ELSE
2260	EE=SIGMA(0.5(EYC2-E45C1)+SIGMA(E45-EY))-E45EQ*EQ/2.
2261	UD=SIGMA(0.5(SQRT2P-E45C1-EYC2)+SIGMA(E45-EY))-E45(0.5+EQ* &
2262	EQ/2.-EQ)
2263	ENDIF
2264	EE=EE/FE
2265	UD=UD/FE
2266
2267	END SUBROUTINE PROF5
2268
2269	!-----------------------------------------------------------------------
2270	SUBROUTINE TPMIX(P,THTU,TU,QU,QLIQ,QICE,QNEWLQ,QNEWIC,RATIO2,RL, &
2271	XLV0,XLV1,XLS0,XLS1, &
2272	EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2273	!-----------------------------------------------------------------------
2274	IMPLICIT NONE
2275	!-----------------------------------------------------------------------
2276	REAL, INTENT(IN ) :: XLV0,XLV1
2277	REAL, INTENT(IN ) :: P,THTU,RATIO2,RL,XLS0, &
2278	XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,&
2279	CICE,DICE
2280	REAL, INTENT(INOUT) :: QU,QLIQ,QICE,TU,QNEWLQ,QNEWIC
2281	REAL :: ES,QS,PI,THTGS,F0,T1,T0,C5,RV,ESLIQ,ESICE,F1,DT,QNEW, &
2282	DQ, QTOT,DQICE,DQLIQ,RLL,CCP
2283	INTEGER :: ITCNT
2284	!-----------------------------------------------------------------------
2285	!
2286	!...THIS SUBROUTINE ITERATIVELY EXTRACTS WET-BULB TEMPERATURE FROM EQUIV
2287	! POTENTIAL TEMPERATURE, THEN CHECKS TO SEE IF SUFFICIENT MOISTURE IS
2288	! AVAILABLE TO ACHIEVE SATURATION...IF NOT, TEMPERATURE IS ADJUSTED
2289	! ACCORDINGLY, IF SO, THE RESIDUAL LIQUID WATER/ICE CONCENTRATION IS
2290	! DETERMINED...
2291	!
2292	C5=1.0723E-3
2293	RV=461.5
2294	!
2295	! ITERATE TO FIND WET BULB TEMPERATURE AS A FUNCTION OF EQUIVALENT POT
2296	! TEMP AND PRS, ASSUMING SATURATION VAPOR PRESSURE...RATIO2 IS THE DEG
2297	! OF GLACIATION...
2298	!
2299	IF(RATIO2.LT.1.E-6)THEN
2300	ES=ALIQEXP((BLIQTU-CLIQ)/(TU-DLIQ))
2301	QS=EP2*ES/(P-ES)
2302	PI=(1.E5/P)*(0.2854(1.-0.28*QS))
2303	THTGS=TUPIEXP((3374.6525/TU-2.5403)QS(1.+0.81*QS))
2304	ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2305	ES=AICEEXP((BICETU-CICE)/(TU-DICE))
2306	QS=EP2*ES/(P-ES)
2307	PI=(1.E5/P)*(0.2854(1.-0.28*QS))
2308	THTGS=TUPIEXP((3114.834/TU-0.278296)QS(1.+0.81*QS))
2309	ELSE
2310	ESLIQ=ALIQEXP((BLIQTU-CLIQ)/(TU-DLIQ))
2311	ESICE=AICEEXP((BICETU-CICE)/(TU-DICE))
2312	ES=(1.-RATIO2)ESLIQ+RATIO2ESICE
2313	QS=EP2*ES/(P-ES)
2314	PI=(1.E5/P)*(0.2854(1.-0.28*QS))
2315	THTGS=TUPIEXP(RLQSC5/TU(1.+0.81QS))
2316	ENDIF
2317	F0=THTGS-THTU
2318	T1=TU-0.5*F0
2319	T0=TU
2320	ITCNT=0
2321	90 IF(RATIO2.LT.1.E-6)THEN
2322	ES=ALIQEXP((BLIQT1-CLIQ)/(T1-DLIQ))
2323	QS=EP2*ES/(P-ES)
2324	PI=(1.E5/P)*(0.2854(1.-0.28*QS))
2325	THTGS=T1PIEXP((3374.6525/T1-2.5403)QS(1.+0.81*QS))
2326	ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2327	ES=AICEEXP((BICET1-CICE)/(T1-DICE))
2328	QS=EP2*ES/(P-ES)
2329	PI=(1.E5/P)*(0.2854(1.-0.28*QS))
2330	THTGS=T1PIEXP((3114.834/T1-0.278296)QS(1.+0.81*QS))
2331	ELSE
2332	ESLIQ=ALIQEXP((BLIQT1-CLIQ)/(T1-DLIQ))
2333	ESICE=AICEEXP((BICET1-CICE)/(T1-DICE))
2334	ES=(1.-RATIO2)ESLIQ+RATIO2ESICE
2335	QS=EP2*ES/(P-ES)
2336	PI=(1.E5/P)*(0.2854(1.-0.28*QS))
2337	THTGS=T1PIEXP(RLQSC5/T1(1.+0.81QS))
2338	ENDIF
2339	F1=THTGS-THTU
2340	IF(ABS(F1).LT.0.01)GOTO 50
2341	ITCNT=ITCNT+1
2342	IF(ITCNT.GT.10)GOTO 50
2343	DT=F1*(T1-T0)/(F1-F0)
2344	T0=T1
2345	F0=F1
2346	T1=T1-DT
2347	GOTO 90
2348	!
2349	! IF THE PARCEL IS SUPERSATURATED, CALCULATE CONCENTRATION OF FRESH
2350	! CONDENSATE...
2351	!
2352	50 IF(QS.LE.QU)THEN
2353	QNEW=QU-QS
2354	QU=QS
2355	GOTO 96
2356	ENDIF
2357	!
2358	! IF THE PARCEL IS SUBSATURATED, TEMPERATURE AND MIXING RATIO MUST BE
2359	! ADJUSTED...IF LIQUID WATER OR ICE IS PRESENT, IT IS ALLOWED TO EVAPO
2360	! SUBLIMATE.
2361	!
2362	QNEW=0.
2363	DQ=QS-QU
2364	QTOT=QLIQ+QICE
2365	!
2366	! IF THERE IS ENOUGH LIQUID OR ICE TO SATURATE THE PARCEL, TEMP STAYS
2367	! WET BULB VALUE, VAPOR MIXING RATIO IS AT SATURATED LEVEL, AND THE MI
2368	! RATIOS OF LIQUID AND ICE ARE ADJUSTED TO MAKE UP THE ORIGINAL SATURA
2369	! DEFICIT... OTHERWISE, ANY AVAILABLE LIQ OR ICE VAPORIZES AND APPROPR
2370	! ADJUSTMENTS TO PARCEL TEMP; VAPOR, LIQUID, AND ICE MIXING RATIOS ARE
2371	!
2372	!...NOTE THAT THE LIQ AND ICE MAY BE PRESENT IN PROPORTIONS SLIGHTLY DIF
2373	! THAN SUGGESTED BY THE VALUE OF RATIO2...CHECK TO MAKE SURE THAT LIQ
2374	! ICE CONCENTRATIONS ARE NOT REDUCED TO BELOW ZERO WHEN EVAPORATION/
2375	! SUBLIMATION OCCURS...
2376	!
2377	IF(QTOT.GE.DQ)THEN
2378	DQICE=0.0
2379	DQLIQ=0.0
2380	QLIQ=QLIQ-(1.-RATIO2)*DQ
2381	IF(QLIQ.LT.0.)THEN
2382	DQICE=0.0-QLIQ
2383	QLIQ=0.0
2384	ENDIF
2385	QICE=QICE-RATIO2*DQ+DQICE
2386	IF(QICE.LT.0.)THEN
2387	DQLIQ=0.0-QICE
2388	QICE=0.0
2389	ENDIF
2390	QLIQ=QLIQ+DQLIQ
2391	QU=QS
2392	GOTO 96
2393	ELSE
2394	IF(RATIO2.LT.1.E-6)THEN
2395	RLL=XLV0-XLV1*T1
2396	ELSEIF(ABS(RATIO2-1.).LT.1.E-6)THEN
2397	RLL=XLS0-XLS1*T1
2398	ELSE
2399	RLL=RL
2400	ENDIF
2401	CCP=1005.7(1.+0.89QU)
2402	IF(QTOT.LT.1.E-10)THEN
2403	!
2404	!...IF NO LIQUID WATER OR ICE IS AVAILABLE, TEMPERATURE IS GIVEN BY:
2405	T1=T1+RLL*(DQ/(1.+DQ))/CCP
2406	GOTO 96
2407	ELSE
2408	!
2409	!...IF SOME LIQ WATER/ICE IS AVAILABLE, BUT NOT ENOUGH TO ACHIEVE SATURA
2410	! THE TEMPERATURE IS GIVEN BY:
2411	T1=T1+RLL*((DQ-QTOT)/(1+DQ-QTOT))/CCP
2412	QU=QU+QTOT
2413	QTOT=0.
2414	ENDIF
2415	QLIQ=0
2416	QICE=0.
2417	ENDIF
2418	96 TU=T1
2419	QNEWLQ=(1.-RATIO2)*QNEW
2420	QNEWIC=RATIO2*QNEW
2421	IF(ITCNT.GT.10)PRINT,'**** NUMBER OF ITERATIONS IN TPMIX =', &
2422	ITCNT
2423
2424	END SUBROUTINE TPMIX
2425	!-----------------------------------------------------------------------
2426	SUBROUTINE ENVIRTHT(P1,T1,Q1,THT1,R1,RL, &
2427	EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2428	!-----------------------------------------------------------------------
2429	IMPLICIT NONE
2430	!-----------------------------------------------------------------------
2431	REAL, INTENT(IN ) :: P1,T1,Q1,R1,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,&
2432	BICE,CICE,DICE
2433	REAL, INTENT(INOUT) :: THT1
2434	REAL:: T00,P00,C1,C2,C3,C4,C5,EE,TLOG,TDPT,TSAT,THT,TFPT,TLOGIC, &
2435	TSATLQ,TSATIC
2436
2437	DATA T00,P00,C1,C2,C3,C4,C5/273.16,1.E5,3374.6525,2.5403,3114.834,&
2438	0.278296,1.0723E-3/
2439	!
2440	! CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMPERATURE...
2441	!
2442
2443	IF(R1.LT.1.E-6)THEN
2444	EE=Q1*P1/(EP2+Q1)
2445	TLOG=ALOG(EE/ALIQ)
2446	TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2447	TSAT=TDPT-(.212+1.571E-3(TDPT-T00)-4.36E-4(T1-T00))*(T1-TDPT)
2448	THT=T1(P00/P1)(0.2854(1.-0.28*Q1))
2449	THT1=THTEXP((C1/TSAT-C2)Q1(1.+0.81Q1))
2450	ELSEIF(ABS(R1-1.).LT.1.E-6)THEN
2451	EE=Q1*P1/(EP2+Q1)
2452	TLOG=ALOG(EE/AICE)
2453	TFPT=(CICE-DICE*TLOG)/(BICE-TLOG)
2454	THT=T1(P00/P1)(0.2854(1.-0.28*Q1))
2455	TSAT=TFPT-(.182+1.13E-3(TFPT-T00)-3.58E-4(T1-T00))*(T1-TFPT)
2456	THT1=THTEXP((C3/TSAT-C4)Q1(1.+0.81Q1))
2457	ELSE
2458	EE=Q1*P1/(EP2+Q1)
2459	TLOG=ALOG(EE/ALIQ)
2460	TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
2461	TLOGIC=ALOG(EE/AICE)
2462	TFPT=(CICE-DICE*TLOGIC)/(BICE-TLOGIC)
2463	THT=T1(P00/P1)(0.2854(1.-0.28*Q1))
2464	TSATLQ=TDPT-(.212+1.571E-3(TDPT-T00)-4.36E-4(T1-T00))*(T1-TDPT)
2465	TSATIC=TFPT-(.182+1.13E-3(TFPT-T00)-3.58E-4(T1-T00))*(T1-TFPT)
2466	TSAT=R1TSATIC+(1.-R1)TSATLQ
2467	THT1=THTEXP(RLQ1C5/TSAT(1.+0.81*Q1))
2468	ENDIF
2469
2470	END SUBROUTINE ENVIRTHT
2471
2472	!-----------------------------------------------------------------------
2473	!*********************** TPDD.FOR **********************************
2474	! THIS SUBROUTINE ITERATIVELY EXTRACTS TEMPERATURE FROM EQUIVALENT *
2475	! POTENTIAL TEMP. IT IS DESIGNED FOR USE WITH DOWNDRAFT CALCULATIONS.
2476	! IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, PARCEL *
2477	! TEMP, SPECIFIC HUMIDITY, AND LIQUID WATER CONTENT ARE ITERATIVELY *
2478	! CALCULATED. *
2479	!***********************************************************************
2480	FUNCTION TPDD(P,THTED,TGS,RS,RD,RH,XLV0,XLV1, &
2481	EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE )
2482	!-----------------------------------------------------------------------
2483	IMPLICIT NONE
2484	!-----------------------------------------------------------------------
2485	REAL, INTENT(IN ) :: XLV0,XLV1
2486	REAL, INTENT(IN ) :: P,THTED,TGS,RD,RH,EP2,ALIQ,BLIQ, &
2487	CLIQ,DLIQ,AICE,BICE,CICE,DICE
2488	REAL, INTENT(INOUT) :: RS
2489	REAL :: TPDD,ES,PI,THTGS,F0,T1,T0,CCP,F1,DT,RL,DSSDT,T1RH,RSRH
2490	INTEGER :: ITCNT
2491	!-----------------------------------------------------------------------
2492	ES=ALIQEXP((BLIQTGS-CLIQ)/(TGS-DLIQ))
2493	RS=EP2*ES/(P-ES)
2494	PI=(1.E5/P)*(0.2854(1.-0.28*RS))
2495	THTGS=TGSPIEXP((3374.6525/TGS-2.5403)RS(1.+0.81*RS))
2496	F0=THTGS-THTED
2497	T1=TGS-0.5*F0
2498	T0=TGS
2499	CCP=1005.7
2500	!
2501	!...ITERATE TO FIND WET-BULB TEMPERATURE...
2502	!
2503	ITCNT=0
2504	90 ES=ALIQEXP((BLIQT1-CLIQ)/(T1-DLIQ))
2505	RS=EP2*ES/(P-ES)
2506	PI=(1.E5/P)*(0.2854(1.-0.28*RS))
2507	THTGS=T1PIEXP((3374.6525/T1-2.5403)RS(1.+0.81*RS))
2508	F1=THTGS-THTED
2509	IF(ABS(F1).LT.0.05)GOTO 50
2510	ITCNT=ITCNT+1
2511	IF(ITCNT.GT.10)GOTO 50
2512	DT=F1*(T1-T0)/(F1-F0)
2513	T0=T1
2514	F0=F1
2515	T1=T1-DT
2516	GOTO 90
2517	50 RL=XLV0-XLV1*T1
2518	!
2519	!...IF RELATIVE HUMIDITY IS SPECIFIED TO BE LESS THAN 100%, ESTIMATE THE
2520	! TEMPERATURE AND MIXING RATIO WHICH WILL YIELD THE APPROPRIATE VALUE.
2521	!
2522	IF(RH.EQ.1.)GOTO 110
2523	DSSDT=(CLIQ-BLIQDLIQ)/((T1-DLIQ)(T1-DLIQ))
2524	DT=RLRS(1.-RH)/(CCP+RLRHRS*DSSDT)
2525	T1RH=T1+DT
2526	ES=RHALIQEXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
2527	RSRH=EP2*ES/(P-ES)
2528	!
2529	!...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
2530	!...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
2531	!
2532	IF(RSRH.LT.RD)THEN
2533	RSRH=RD
2534	T1RH=T1+(RS-RSRH)*RL/CCP
2535	ENDIF
2536	T1=T1RH
2537	RS=RSRH
2538	110 TPDD=T1
2539	IF(ITCNT.GT.10)PRINT,'**** NUMBER OF ITERATIONS IN TPDD = ', &
2540	ITCNT
2541
2542	END FUNCTION TPDD
2543
2544	!====================================================================
2545	SUBROUTINE kfinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN, &
2546	RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QI,P_QS, &
2547	P_FIRST_SCALAR,restart,allowed_to_read, &
2548	ids, ide, jds, jde, kds, kde, &
2549	ims, ime, jms, jme, kms, kme, &
2550	its, ite, jts, jte, kts, kte )
2551	!--------------------------------------------------------------------
2552	IMPLICIT NONE
2553	!--------------------------------------------------------------------
2554	LOGICAL , INTENT(IN) :: restart, allowed_to_read
2555	INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
2556	ims, ime, jms, jme, kms, kme, &
2557	its, ite, jts, jte, kts, kte
2558	INTEGER , INTENT(IN) :: P_QI,P_QS,P_FIRST_SCALAR
2559
2560	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
2561	RTHCUTEN, &
2562	RQVCUTEN, &
2563	RQCCUTEN, &
2564	RQRCUTEN, &
2565	RQICUTEN, &
2566	RQSCUTEN
2567
2568	REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
2569
2570	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
2571
2572	INTEGER :: i, j, k, itf, jtf, ktf
2573
2574	jtf=min0(jte,jde-1)
2575	ktf=min0(kte,kde-1)
2576	itf=min0(ite,ide-1)
2577
2578	IF(.not.restart)THEN
2579	DO j=jts,jtf
2580	DO k=kts,ktf
2581	DO i=its,itf
2582	RTHCUTEN(i,k,j)=0.
2583	RQVCUTEN(i,k,j)=0.
2584	RQCCUTEN(i,k,j)=0.
2585	RQRCUTEN(i,k,j)=0.
2586	ENDDO
2587	ENDDO
2588	ENDDO
2589
2590	IF (P_QI .ge. P_FIRST_SCALAR) THEN
2591	DO j=jts,jtf
2592	DO k=kts,ktf
2593	DO i=its,itf
2594	RQICUTEN(i,k,j)=0.
2595	ENDDO
2596	ENDDO
2597	ENDDO
2598	ENDIF
2599
2600	IF (P_QS .ge. P_FIRST_SCALAR) THEN
2601	DO j=jts,jtf
2602	DO k=kts,ktf
2603	DO i=its,itf
2604	RQSCUTEN(i,k,j)=0.
2605	ENDDO
2606	ENDDO
2607	ENDDO
2608	ENDIF
2609
2610	DO j=jts,jtf
2611	DO i=its,itf
2612	NCA(i,j)=-100.
2613	ENDDO
2614	ENDDO
2615
2616	DO j=jts,jtf
2617	DO k=kts,ktf
2618	DO i=its,itf
2619	W0AVG(i,k,j)=0.
2620	ENDDO
2621	ENDDO
2622	ENDDO
2623
2624	ENDIF
2625
2626	END SUBROUTINE kfinit
2627
2628	!-------------------------------------------------------
2629
2630	END MODULE module_cu_kf
2631

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

Download in other formats:

Original Format