Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_kf.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: 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
195	! Old:
196	!
197	! W0AVG(I,K,J)=(W0AVG(I,K,J)*(TST-1.)+W0)/TST
198	! New, to support adaptive time step:
199	!
200	W0AVG(I,K,J) = ( W0AVG(I,K,J) * W0AVGfctr + W0 * W0fctr ) / W0den
201
202	ENDDO
203	ENDDO
204	ENDDO
205	lastdt = dt
206	!
207	!...CHECK FOR CONVECTIVE INITIATION EVERY 5 MINUTES (OR NTST/2)...
208	!
209
210	!
211	! Modified for adaptive time step
212	!
213	if (ADAPT_STEP_FLAG) then
214	if ( (KTAU .eq. 1) .or. (cudt .eq. 0) .or. &
215	( CURR_SECS + dt >= &
216	( int( CURR_SECS / ( cudt * 60 ) ) + 1 ) * cudt * 60 ) ) then
217	run_param = .TRUE.
218	else
219	run_param = .FALSE.
220	endif
221
222	else
223	if (MOD(KTAU,NTST) .EQ. 0 .or. KTAU .eq. 1) then
224	run_param = .TRUE.
225	else
226	run_param = .FALSE.
227	endif
228	endif
229
230	IF (run_param) then
231	DO J = jts,jte
232	DO I= its,ite
233	CU_ACT_FLAG(i,j) = .true.
234	ENDDO
235	ENDDO
236
237	DO J = jts,jte
238	DO I=its,ite
239	! if (i.eq. 110 .and. j .eq. 59 ) then
240	! write(0,*) 'nca = ',nca(i,j),' CU_ACT_FLAG = ',CU_ACT_FLAG(i,j)
241	! write(0,*) 'dt = ',dt,' ADAPT_STEP_FLAG = ',ADAPT_STEP_FLAG
242	! endif
243	! IF ( NINT(NCA(I,J)) .gt. 0 ) then
244	IF ( NCA(I,J) .gt. 0.5*DT ) then
245	CU_ACT_FLAG(i,j) = .false.
246	ELSE
247
248	DO k=kts,kte
249	DQDT(k)=0.
250	DQIDT(k)=0.
251	DQCDT(k)=0.
252	DQRDT(k)=0.
253	DQSDT(k)=0.
254	DTDT(k)=0.
255	ENDDO
256	RAINCV(I,J)=0.
257	PRATEC(I,J)=0.
258	!
259	! assign vars from 3D to 1D
260
261	DO K=kts,kte
262	U1D(K) =U(I,K,J)
263	V1D(K) =V(I,K,J)
264	T1D(K) =T(I,K,J)
265	RHO1D(K) =rho(I,K,J)
266	QV1D(K)=QV(I,K,J)
267	P1D(K) =Pcps(I,K,J)
268	W0AVG1D(K) =W0AVG(I,K,J)
269	DZ1D(k)=dz8w(I,K,J)
270	ENDDO
271
272	!
273	CALL KFPARA(I, J, &
274	U1D,V1D,T1D,QV1D,P1D,DZ1D, &
275	W0AVG1D,DT,DX,DXSQ,RHO1D, &
276	XLV0,XLV1,XLS0,XLS1,CP,R,G, &
277	EP2,SVP1,SVP2,SVP3,SVPT0, &
278	DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
279	RAINCV,PRATEC,NCA, &
280	warm_rain,qi_flag,qs_flag, &
281	ids,ide, jds,jde, kds,kde, &
282	ims,ime, jms,jme, kms,kme, &
283	its,ite, jts,jte, kts,kte )
284
285	IF ( PRESENT( RTHCUTEN ) .AND. PRESENT( RQVCUTEN ) ) THEN
286	DO K=kts,kte
287	RTHCUTEN(I,K,J)=DTDT(K)/pi(I,K,J)
288	RQVCUTEN(I,K,J)=DQDT(K)
289	ENDDO
290	ENDIF
291
292	IF( PRESENT(RQRCUTEN) .AND. PRESENT(RQCCUTEN) .AND. &
293	PRESENT(F_QR) ) THEN
294	IF ( F_QR ) THEN
295	DO K=kts,kte
296	RQRCUTEN(I,K,J)=DQRDT(K)
297	RQCCUTEN(I,K,J)=DQCDT(K)
298	ENDDO
299	ELSE
300	! This is the case for Eta microphysics without 3d rain field
301	DO K=kts,kte
302	RQRCUTEN(I,K,J)=0.
303	RQCCUTEN(I,K,J)=DQRDT(K)+DQCDT(K)
304	ENDDO
305	ENDIF
306	ENDIF
307
308	!...... QSTEN STORES GRAUPEL TENDENCY IF IT EXISTS, OTHERISE SNOW (V2)
309
310	IF( PRESENT( RQICUTEN ) .AND. qi_flag )THEN
311	DO K=kts,kte
312	RQICUTEN(I,K,J)=DQIDT(K)
313	ENDDO
314	ENDIF
315
316	IF( PRESENT ( RQSCUTEN ) .AND. qs_flag )THEN
317	DO K=kts,kte
318	RQSCUTEN(I,K,J)=DQSDT(K)
319	ENDDO
320	ENDIF
321	!
322	ENDIF
323	ENDDO
324	ENDDO
325
326	ENDIF
327
328	END SUBROUTINE KFCPS
329
330	!-----------------------------------------------------------
331	SUBROUTINE KFPARA (I, J, &
332	U0,V0,T0,QV0,P0,DZQ,W0AVG1D, &
333	DT,DX,DXSQ,rho, &
334	XLV0,XLV1,XLS0,XLS1,CP,R,G, &
335	EP2,SVP1,SVP2,SVP3,SVPT0, &
336	DQDT,DQIDT,DQCDT,DQRDT,DQSDT,DTDT, &
337	RAINCV,PRATEC,NCA, &
338	warm_rain,qi_flag,qs_flag, &
339	ids,ide, jds,jde, kds,kde, &
340	ims,ime, jms,jme, kms,kme, &
341	its,ite, jts,jte, kts,kte )
342	!-----------------------------------------------------------
343	IMPLICIT NONE
344	!-----------------------------------------------------------
345	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
346	ims,ime, jms,jme, kms,kme, &
347	its,ite, jts,jte, kts,kte, &
348	I,J
349	LOGICAL, INTENT(IN ) :: warm_rain
350	LOGICAL :: qi_flag, qs_flag
351
352	!
353	REAL, DIMENSION( kts:kte ), &
354	INTENT(IN ) :: U0, &
355	V0, &
356	T0, &
357	QV0, &
358	P0, &
359	rho, &
360	DZQ, &
361	W0AVG1D
362	!
363	REAL, INTENT(IN ) :: DT,DX,DXSQ
364	!
365
366	REAL, INTENT(IN ) :: XLV0,XLV1,XLS0,XLS1,CP,R,G
367	REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
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) :: RAINCV, &
379	PRATEC, &
380	NCA
381	!
382	!...DEFINE LOCAL VARIABLES...
383	!
384	REAL, DIMENSION( kts:kte ) :: &
385	Q0,Z0,TV0,TU,TVU,QU,TZ,TVD, &
386	QD,QES,THTES,TG,TVG,QG,WU,WD,W0,EMS,EMSD, &
387	UMF,UER,UDR,DMF,DER,DDR,UMF2,UER2, &
388	UDR2,DMF2,DER2,DDR2,DZA,THTA0,THETEE, &
389	THTAU,THETEU,THTAD,THETED,QLIQ,QICE, &
390	QLQOUT,QICOUT,PPTLIQ,PPTICE,DETLQ,DETIC, &
391	DETLQ2,DETIC2,RATIO,RATIO2
392
393	REAL, DIMENSION( kts:kte ) :: &
394	DOMGDP,EXN,RHOE,TVQU,DP,RH,EQFRC,WSPD, &
395	QDT,FXM,THTAG,THTESG,THPA,THFXTOP, &
396	THFXBOT,QPA,QFXTOP,QFXBOT,QLPA,QLFXIN, &
397	QLFXOUT,QIPA,QIFXIN,QIFXOUT,QRPA, &
398	QRFXIN,QRFXOUT,QSPA,QSFXIN,QSFXOUT, &
399	QL0,QLG,QI0,QIG,QR0,QRG,QS0,QSG
400
401	REAL, DIMENSION( kts:kte+1 ) :: OMG
402	REAL, DIMENSION( kts:kte ) :: RAINFB,SNOWFB
403
404	! LOCAL VARS
405
406	REAL :: P00,T00,CV,B61,RLF,RHIC,RHBC,PIE, &
407	TTFRZ,TBFRZ,C5,RATE
408	REAL :: GDRY,ROCP,ALIQ,BLIQ, &
409	CLIQ,DLIQ,AICE,BICE,CICE,DICE
410	REAL :: FBFRC,P300,DPTHMX,THMIX,QMIX,ZMIX,PMIX, &
411	ROCPQ,TMIX,EMIX,TLOG,TDPT,TLCL,TVLCL, &
412	CPORQ,PLCL,ES,DLP,TENV,QENV,TVEN,TVBAR, &
413	ZLCL,WKL,WABS,TRPPT,WSIGNE,DTLCL,GDT,WLCL,&
414	TVAVG,QESE,WTW,RHOLCL,AU0,VMFLCL,UPOLD, &
415	UPNEW,ABE,WKLCL,THTUDL,TUDL,TTEMP,FRC1, &
416	QNEWIC,RL,R1,QNWFRZ,EFFQ,BE,BOTERM,ENTERM,&
417	DZZ,WSQ,UDLBE,REI,EE2,UD2,TTMP,F1,F2, &
418	THTTMP,QTMP,TMPLIQ,TMPICE,TU95,TU10,EE1, &
419	UD1,CLDHGT,DPTT,QNEWLQ,DUMFDP,EE,TSAT, &
420	THTA,P150,USR,VCONV,TIMEC,SHSIGN,VWS,PEF, &
421	CBH,RCBH,PEFCBH,PEFF,PEFF2,TDER,THTMIN, &
422	DTMLTD,QS,TADVEC,DPDD,FRC,DPT,RDD,A1, &
423	DSSDT,DTMP,T1RH,QSRH,PPTFLX,CPR,CNDTNF, &
424	UPDINC,AINCM2,DEVDMF,PPR,RCED,DPPTDF, &
425	DMFLFS,DMFLFS2,RCED2,DDINC,AINCMX,AINCM1, &
426	AINC,TDER2,PPTFL2,FABE,STAB,DTT,DTT1, &
427	DTIME,TMA,TMB,TMM,BCOEFF,ACOEFF,QVDIFF, &
428	TOPOMG,CPM,DQ,ABEG,DABE,DFDA,FRC2,DR, &
429	UDFRC,TUC,QGS,RH0,RHG,QINIT,QFNL,ERR2, &
430	RELERR,RLC,RLS,RNC,FABEOLD,AINCOLD,UEFRC, &
431	DDFRC,TDC,DEFRC
432
433	INTEGER :: KX,K,KL
434	!
435	INTEGER :: ISTOP,ML,L5,L4,KMIX,LOW, &
436	LC,MXLAYR,LLFC,NLAYRS,NK, &
437	KPBL,KLCL,LCL,LET,IFLAG, &
438	KFRZ,NK1,LTOP,NJ,LTOP1, &
439	LTOPM1,LVF,KSTART,KMIN,LFS, &
440	ND,NIC,LDB,LDT,ND1,NDK, &
441	NM,LMAX,NCOUNT,NOITR, &
442	NSTEP,NTC
443	!
444	DATA P00,T00/1.E5,273.16/
445	DATA CV,B61,RLF/717.,0.608,3.339E5/
446	DATA RHIC,RHBC/1.,0.90/
447	DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
448	DATA RATE/0.01/
449	!-----------------------------------------------------------
450	GDRY=-G/CP
451	ROCP=R/CP
452	KL=kte
453	KX=kte
454	!
455	! ALIQ = 613.3
456	! BLIQ = 17.502
457	! CLIQ = 4780.8
458	! DLIQ = 32.19
459	ALIQ = SVP1*1000.
460	BLIQ = SVP2
461	CLIQ = SVP2*SVPT0
462	DLIQ = SVP3
463	AICE = 613.2
464	BICE = 22.452
465	CICE = 6133.0
466	DICE = 0.61
467	!
468
469	!...OPTION TO FEED CONVECTIVELY GENERATED RAINWATER
470	!...INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL)
471	!...FIELD. 'FBFRC' IS THE FRACTION OF AVAILABLE
472	!...PRECIPITATION TO BE FED BACK (0.0 - 1.0)...
473	!
474	FBFRC=0.0
475	!
476	!...SCHEME IS CALLED ONCE ON EACH NORTH-SOUTH SLICE, THE LOOP BELOW
477	!...CHECKS FOR THE POSSIBILITY OF INITIATING PARAMETERIZED
478	!...CONVECTION AT EACH POINT WITHIN THE SLICE
479	!
480	!...SEE IF IT IS NECESSARY TO CHECK FOR CONVECTIVE TRIGGERING AT THIS
481	!...GRID POINT. IF NCA>0, CONVECTION IS ALREADY ACTIVE AT THIS POINT,
482	!...JUST FEED BACK THE TENDENCIES SAVED FROM THE TIME WHEN CONVECTION
483	!...WAS INITIATED. IF NCA<0, CONVECTION IS NOT ACTIVE
484	!...AND YOU MAY WANT TO CHECK TO SEE IF IT CAN BE ACTIVATED FOR THE
485	!...CURRENT CONDITIONS. IN PREVIOUS APLICATIONS OF THIS SCHEME,
486	!...THE VARIABLE ICLDCK WAS USED BELOW TO SAVE TIME BY ONLY CHECKING
487	!...FOR THE POSSIBILITY OF CONVECTIVE INITIATION EVERY 5 OR 10
488	!...MINUTES...
489	!
490
491	! 10 CONTINUE
492	!SUE P300=1000.(PSB(I,J)A(KL)+PTOP-30.)+PP3D(I,J,KL)
493
494	P300=P0(1)-30000.
495	!
496	!...PRESSURE PERTURBATION TERM IS ONLY DEFINED AT MID-POINT OF
497	!...VERTICAL LAYERS...SINCE TOTAL PRESSURE IS NEEDED AT THE TOP AND
498	!...BOTTOM OF LAYERS BELOW, DO AN INTERPOLATION...
499	!
500	!...INPUT A VERTICAL SOUNDING ... NOTE THAT MODEL LAYERS ARE NUMBERED
501	!...FROM BOTTOM-UP IN THE KF SCHEME...
502	!
503	ML=0
504	!SUE tmprpsb=1./PSB(I,J)
505	!SUE CELL=PTOP*tmprpsb
506
507	DO 15 K=1,KX
508	!SUE P0(K)=1.E3(A(NK)PSB(I,J)+PTOP)+PP3D(I,J,NK)
509	!
510	!...IF Q0 IS ABOVE SATURATION VALUE, REDUCE IT TO SATURATION LEVEL...
511	!
512	ES=ALIQEXP((BLIQT0(K)-CLIQ)/(T0(K)-DLIQ))
513	QES(K)=EP2*ES/(P0(K)-ES)
514	Q0(K)=AMIN1(QES(K),QV0(K))
515	Q0(K)=AMAX1(0.000001,Q0(K))
516	QL0(K)=0.
517	QI0(K)=0.
518	QR0(K)=0.
519	QS0(K)=0.
520
521	TV0(K)=T0(K)(1.+B61Q0(K))
522	RHOE(K)=P0(K)/(R*TV0(K))
523
524	DP(K)=rho(k)gDZQ(k)
525	!
526	!...DZQ IS DZ BETWEEN SIGMA SURFACES, DZA IS DZ BETWEEN MODEL HALF LEVEL
527	! DP IS THE PRESSURE INTERVAL BETWEEN FULL SIGMA LEVELS...
528	!
529	IF(P0(K).GE.500E2)L5=K
530	IF(P0(K).GE.400E2)L4=K
531	IF(P0(K).GE.P300)LLFC=K
532	IF(T0(K).GT.T00)ML=K
533	15 CONTINUE
534
535	Z0(1)=.5*DZQ(1)
536	DO 20 K=2,KL
537	Z0(K)=Z0(K-1)+.5*(DZQ(K)+DZQ(K-1))
538	DZA(K-1)=Z0(K)-Z0(K-1)
539	20 CONTINUE
540	DZA(KL)=0.
541	KMIX=1
542	25 LOW=KMIX
543
544	IF(LOW.GT.LLFC)GOTO 325
545
546	LC=LOW
547	MXLAYR=0
548	!
549	!...ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
550	!...UNSTABLE AIR 50 TO 100 mb DEEP...TO APPROXIMATE THIS, ISOLATE A
551	!...GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
552	!...LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 60 mb..
553	!
554	NLAYRS=0
555	DPTHMX=0.
556	DO 63 NK=LC,KX
557	DPTHMX=DPTHMX+DP(NK)
558	NLAYRS=NLAYRS+1
559	63 IF(DPTHMX.GT.6.E3)GOTO 64
560	GOTO 325
561	64 KPBL=LC+NLAYRS-1
562	KMIX=LC+1
563	18 THMIX=0.
564	QMIX=0.
565	ZMIX=0.
566	PMIX=0.
567	DPTHMX=0.
568	!
569	!...FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
570	!...MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
571	!...LAYERS...
572	!
573	DO 17 NK=LC,KPBL
574	DPTHMX=DPTHMX+DP(NK)
575	ROCPQ=0.2854(1.-0.28Q0(NK))
576	THMIX=THMIX+DP(NK)T0(NK)(P00/P0(NK))**ROCPQ
577	QMIX=QMIX+DP(NK)*Q0(NK)
578	ZMIX=ZMIX+DP(NK)*Z0(NK)
579	17 PMIX=PMIX+DP(NK)*P0(NK)
580	THMIX=THMIX/DPTHMX
581	QMIX=QMIX/DPTHMX
582	ZMIX=ZMIX/DPTHMX
583	PMIX=PMIX/DPTHMX
584	ROCPQ=0.2854(1.-0.28QMIX)
585	TMIX=THMIX(PMIX/P00)*ROCPQ
586	EMIX=QMIX*PMIX/(EP2+QMIX)
587	!
588	!...FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL, PRESSURE
589	!...LEVEL OF LCL...
590	!
591	TLOG=ALOG(EMIX/ALIQ)
592	TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
593	TLCL=TDPT-(.212+1.571E-3(TDPT-T00)-4.36E-4(TMIX-T00))*(TMIX- &
594	TDPT)
595	TLCL=AMIN1(TLCL,TMIX)
596	TVLCL=TLCL(1.+0.608QMIX)
597	CPORQ=1./ROCPQ
598	PLCL=P00(TLCL/THMIX)*CPORQ
599	DO 29 NK=LC,KL
600	KLCL=NK
601	IF(PLCL.GE.P0(NK))GOTO 35
602	29 CONTINUE
603	GOTO 325
604	35 K=KLCL-1
605	DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
606	!
607	!...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
608	!
609	TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
610	QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
611	TVEN=TENV(1.+0.608QENV)
612	TVBAR=0.5*(TV0(K)+TVEN)
613	! ZLCL=Z0(K)+RTVBARALOG(P0(K)/PLCL)/G
614	ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
615	!
616	!...CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
617	!...FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992)...W0AVG IS AN
618	!...APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
619	!...VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
620	!...THAN THE INSTANTANEOUS VALUE...FORMULA RELATING TEMPERATURE
621	!...PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
622	!...SUCCESS AT GRID LENGTHS NEAR 25 km. FOR DIFFERENT GRID-LENGTHS,
623	!...ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
624	!...LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH...
625	!
626	WKLCL=0.02*ZLCL/2.5E3
627	WKL=(W0AVG1D(K)+(W0AVG1D(KLCL)-W0AVG1D(K))DLP)DX/25.E3- &
628	WKLCL
629	WABS=ABS(WKL)+1.E-10
630	WSIGNE=WKL/WABS
631	DTLCL=4.64WSIGNEWABS**0.33
632	GDT=GDTLCL(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
633	WLCL=1.+.5WSIGNESQRT(ABS(GDT)+1.E-10)
634	IF(TLCL+DTLCL.GT.TENV)GOTO 45
635	IF(KPBL.GE.LLFC)GOTO 325
636	GOTO 25
637	!
638	!...CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED...COMPUTE
639	!...EQUIVALENT POTENTIAL TEMPERATURE
640	!...(THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL...
641	!
642	45 THETEU(K)=TMIX(1.E5/PMIX)(0.2854(1.-0.28QMIX)) &
643	EXP((3374.6525/TLCL-2.5403)QMIX(1.+0.81*QMIX))
644	ES=ALIQEXP((TENVBLIQ-CLIQ)/(TENV-DLIQ))
645	TVAVG=0.5(TV0(KLCL)+TENV(1.+0.608*QENV))
646	PLCL=P0(KLCL)EXP(G/(RTVAVG)*(Z0(KLCL)-ZLCL))
647	QESE=EP2*ES/(PLCL-ES)
648	GDT=GDTLCL(ZLCL-Z0(LC))/(TV0(LC)+TVEN)
649	WLCL=1.+.5WSIGNESQRT(ABS(GDT)+1.E-10)
650	THTES(K)=TENV(1.E5/PLCL)(0.2854(1.-0.28QESE)) &
651	EXP((3374.6525/TENV-2.5403)QESE(1.+0.81*QESE))
652	WTW=WLCL*WLCL
653	IF(WLCL.LT.0.)GOTO 25
654	TVLCL=TLCL(1.+0.608QMIX)
655	RHOLCL=PLCL/(R*TVLCL)
656	!
657	LCL=KLCL
658	LET=LCL
659	!
660	!*******************************************************************
661	! *
662	! COMPUTE UPDRAFT PROPERTIES *
663	! *
664	!*******************************************************************
665	!
666	!
667	!...ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))...
668	!
669	WU(K)=WLCL
670	AU0=PIERADRAD
671	UMF(K)=RHOLCL*AU0
672	VMFLCL=UMF(K)
673	UPOLD=VMFLCL
674	UPNEW=UPOLD
675	!
676	!...RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
677	!...UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE BUOYANT ENERGY,
678	! TRPPT IS THE TOTAL RATE OF PRECIPITATION PRODUCTION...
679	!
680	RATIO2(K)=0.
681	UER(K)=0.
682	ABE=0.
683	TRPPT=0.
684	TU(K)=TLCL
685	TVU(K)=TVLCL
686	QU(K)=QMIX
687	EQFRC(K)=1.
688	QLIQ(K)=0.
689	QICE(K)=0.
690	QLQOUT(K)=0.
691	QICOUT(K)=0.
692	DETLQ(K)=0.
693	DETIC(K)=0.
694	PPTLIQ(K)=0.
695	PPTICE(K)=0.
696	IFLAG=0
697	KFRZ=LC
698	!
699	!...THE AMOUNT OF CONV AVAIL POT ENERGY (CAPE) IS CALCULATED WITH
700	! RESPECT TO UNDILUTE PARCEL ASCENT; EQ POT TEMP OF UNDILUTE
701	! PARCEL IS THTUDL, UNDILUTE TEMPERATURE IS GIVEN BY TUDL...
702	!
703	THTUDL=THETEU(K)
704	TUDL=TLCL
705	!
706	!...TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
707	! PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
708	! FREEZING IS SPECIFIED TO BEGIN. WITHIN THE GLACIATION
709	! INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
710	! PREVIOUS MODEL LEVEL...
711	!
712	TTEMP=TTFRZ
713	!
714	!...ENTER THE LOOP FOR UPDRAFT CALCULATIONS...CALCULATE UPDRAFT TEMP,
715	! MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
716	! MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL...
717	!
718	DO 60 NK=K,KL-1
719	NK1=NK+1
720	RATIO2(NK1)=RATIO2(NK)
721	!
722	!...UPDATE UPDRAFT PROPERTIES AT THE NEXT MODEL LVL TO REFLECT
723	! ENTRAINMENT OF ENVIRONMENTAL AIR...
724	!
725	FRC1=0.
726	TU(NK1)=T0(NK1)
727	THETEU(NK1)=THETEU(NK)
728	QU(NK1)=QU(NK)
729	QLIQ(NK1)=QLIQ(NK)
730	QICE(NK1)=QICE(NK)
731
732	CALL TPMIX(P0(NK1),THETEU(NK1),TU(NK1),QU(NK1),QLIQ(NK1), &
733	QICE(NK1),QNEWLQ,QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0, &
734	XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
735	TVU(NK1)=TU(NK1)(1.+0.608QU(NK1))
736	!
737	!...CHECK TO SEE IF UPDRAFT TEMP IS WITHIN THE FREEZING INTERVAL,
738	! IF IT IS, CALCULATE THE FRACTIONAL CONVERSION TO GLACIATION
739	! AND ADJUST QNEWLQ TO REFLECT THE GRADUAL CHANGE IN THETAU
740	! SINCE THE LAST MODEL LEVEL...THE GLACIATION EFFECTS WILL BE
741	! DETERMINED AFTER THE AMOUNT OF CONDENSATE AVAILABLE AFTER
742	! PRECIP FALLOUT IS DETERMINED...TTFRZ IS THE TEMP AT WHICH
743	! GLACIATION BEGINS, TBFRZ THE TEMP AT WHICH IT ENDS...
744	!
745	IF(TU(NK1).LE.TTFRZ.AND.IFLAG.LT.1)THEN
746	IF(TU(NK1).GT.TBFRZ)THEN
747	IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
748	FRC1=(TTEMP-TU(NK1))/(TTFRZ-TBFRZ)
749	R1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
750	ELSE
751	FRC1=(TTEMP-TBFRZ)/(TTFRZ-TBFRZ)
752	R1=1.
753	IFLAG=1
754	ENDIF
755	QNWFRZ=QNEWLQ
756	QNEWIC=QNEWIC+QNEWLQR10.5
757	QNEWLQ=QNEWLQ-QNEWLQR10.5
758	EFFQ=(TTFRZ-TBFRZ)/(TTEMP-TBFRZ)
759	TTEMP=TU(NK1)
760	ENDIF
761	!
762	! CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
763	!
764	IF(NK.EQ.K)THEN
765	BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
766	BOTERM=2.(Z0(NK1)-ZLCL)G*BE/1.5
767	ENTERM=0.
768	DZZ=Z0(NK1)-ZLCL
769	ELSE
770	BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
771	BOTERM=2.DZA(NK)G*BE/1.5
772	ENTERM=2.UER(NK)WTW/UPOLD
773	DZZ=DZA(NK)
774	ENDIF
775	WSQ=WTW
776	CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM,RATE, &
777	QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1), G)
778
779	!...IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
780	! IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS...
781	!
782	IF(WTW.LE.0.)GOTO 65
783	WABS=SQRT(ABS(WTW))
784	WU(NK1)=WTW/WABS
785	!
786	! UPDATE THE ABE FOR UNDILUTE ASCENT...
787	!
788	THTES(NK1)=T0(NK1)(1.E5/P0(NK1))(0.2854(1.-0.28*QES(NK1))) &
789	* &
790	EXP((3374.6525/T0(NK1)-2.5403)QES(NK1)(1.+0.81* &
791	QES(NK1)))
792	UDLBE=((2.THTUDL)/(THTES(NK)+THTES(NK1))-1.)DZZ
793	IF(UDLBE.GT.0.)ABE=ABE+UDLBE*G
794	!
795	! DETERMINE THE EFFECTS OF CLOUD GLACIATION IF WITHIN THE SPECIFIED
796	! TEMP INTERVAL...
797	!
798	IF(FRC1.GT.1.E-6)THEN
799	CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QLIQ(NK1), &
800	QICE(NK1),RATIO2(NK1),TTFRZ,TBFRZ,QNWFRZ,RL,FRC1,EFFQ, &
801	IFLAG,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE &
802	,CICE,DICE)
803	ENDIF
804	!
805	! CALL SUBROUTINE TO CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL TEMP.
806	! WITHIN GLACIATION INTERVAL, THETAE MUST BE CALCULATED WITH RESPECT TO
807	! SAME DEGREE OF GLACIATION FOR ALL ENTRAINING AIR...
808	!
809	CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),RATIO2(NK1), &
810	RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
811
812	!...REI IS THE RATE OF ENVIRONMENTAL INFLOW...
813	!
814	REI=VMFLCLDP(NK1)0.03/RAD
815	TVQU(NK1)=TU(NK1)(1.+0.608QU(NK1)-QLIQ(NK1)-QICE(NK1))
816	!
817	!...IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, NO
818	! ENTRAINMENT IS ALLOWED AT THIS LEVEL...
819	!
820	IF(TVQU(NK1).LE.TV0(NK1))THEN
821	UER(NK1)=0.0
822	UDR(NK1)=REI
823	EE2=0.
824	UD2=1.
825	EQFRC(NK1)=0.
826	GOTO 55
827	ENDIF
828	LET=NK1
829	TTMP=TVQU(NK1)
830	!
831	!...DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT AND ENVIRONMENTAL
832	! AIR FOR ESTIMATION OF ENTRAINMENT AND DETRAINMENT RATES...
833	!
834	F1=0.95
835	F2=1.-F1
836	THTTMP=F1THETEE(NK1)+F2THETEU(NK1)
837	QTMP=F1Q0(NK1)+F2QU(NK1)
838	TMPLIQ=F2*QLIQ(NK1)
839	TMPICE=F2*QICE(NK1)
840	CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
841	QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
842	DLIQ,AICE,BICE,CICE,DICE)
843	TU95=TTMP(1.+0.608QTMP-TMPLIQ-TMPICE)
844	IF(TU95.GT.TV0(NK1))THEN
845	EE2=1.
846	UD2=0.
847	EQFRC(NK1)=1.0
848	GOTO 50
849	ENDIF
850	F1=0.10
851	F2=1.-F1
852	THTTMP=F1THETEE(NK1)+F2THETEU(NK1)
853	QTMP=F1Q0(NK1)+F2QU(NK1)
854	TMPLIQ=F2*QLIQ(NK1)
855	TMPICE=F2*QICE(NK1)
856	CALL TPMIX(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,QNEWLQ, &
857	QNEWIC,RATIO2(NK1),RL,XLV0,XLV1,XLS0,XLS1,EP2,ALIQ,BLIQ,CLIQ, &
858	DLIQ,AICE,BICE,CICE,DICE)
859	TU10=TTMP(1.+0.608QTMP-TMPLIQ-TMPICE)
860	IF(TU10.EQ.TVQU(NK1))THEN
861	EE2=1.
862	UD2=0.
863	EQFRC(NK1)=1.0
864	GOTO 50
865	ENDIF
866	EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
867	EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
868	EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
869	IF(EQFRC(NK1).EQ.1)THEN
870	EE2=1.
871	UD2=0.
872	GOTO 50
873	ELSEIF(EQFRC(NK1).EQ.0.)THEN
874	EE2=0.
875	UD2=1.
876	GOTO 50
877	ELSE
878	!
879	!...SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
880	! FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES...
881	!
882	CALL PROF5(EQFRC(NK1),EE2,UD2)
883	ENDIF
884	!
885	50 IF(NK.EQ.K)THEN
886	EE1=1.
887	UD1=0.
888	ENDIF
889	!
890	!...NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE
891	! FRACTIONAL VALUES IN THE LAYER...
892	!
893	UER(NK1)=0.5REI(EE1+EE2)
894	UDR(NK1)=0.5REI(UD1+UD2)
895	!
896	!...IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
897	! UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATION
898	!
899	55 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
900	!
901	!...IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL
902	! UPDRAFT FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT THE
903	! PREVIOUS MODEL
904	!
905	IF(UDLBE.GT.0.)ABE=ABE-UDLBE*G
906	LET=NK
907	! WRITE(98,1015)P0(NK1)/100.
908	GOTO 65
909	ENDIF
910	EE1=EE2
911	UD1=UD2
912	UPOLD=UMF(NK)-UDR(NK1)
913	UPNEW=UPOLD+UER(NK1)
914	UMF(NK1)=UPNEW
915	!
916	!...DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND ICE IN
917	! THE DETRAINING UPDRAFT MASS...
918	!
919	DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
920	DETIC(NK1)=QICE(NK1)*UDR(NK1)
921	QDT(NK1)=QU(NK1)
922	QU(NK1)=(UPOLDQU(NK1)+UER(NK1)Q0(NK1))/UPNEW
923	THETEU(NK1)=(THETEU(NK1)UPOLD+THETEE(NK1)UER(NK1))/UPNEW
924	QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
925	QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
926	!
927	!...KFRZ IS THE HIGHEST MODEL LEVEL AT WHICH LIQUID CONDENSATE IS
928	! GENERATING PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF LIQUID
929	! PRECIP AT A GIVING MODEL LVL, PPTICE THE SAME FOR ICE, TRPPT IS
930	! THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE CURRENT MODEL LEVEL
931	!
932	IF(ABS(RATIO2(NK1)-1.).GT.1.E-6)KFRZ=NK1
933	PPTLIQ(NK1)=QLQOUT(NK1)*(UMF(NK)-UDR(NK1))
934	PPTICE(NK1)=QICOUT(NK1)*(UMF(NK)-UDR(NK1))
935	TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
936	IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
937	60 CONTINUE
938	!
939	!...CHECK CLOUD DEPTH...IF CLOUD IS TALL ENOUGH, ESTIMATE THE EQUILIBRIU
940	! TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE AT CLOUD TOP SO
941	! THAT MASS FLUX DECREASES TO ZERO AS A LINEAR FUNCTION OF PRESSURE
942	! BETWEEN THE LET AND CLOUD TOP...
943	!
944	!...LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL
945	! VELOCITY FIRST BECOMES NEGATIVE...
946	!
947	65 LTOP=NK
948	CLDHGT=Z0(LTOP)-ZLCL
949	!
950	!...IF CLOUD TOP HGT IS LESS THAN SPECIFIED MINIMUM HEIGHT, GO BACK AND
951	! THE NEXT HIGHEST 60MB LAYER TO SEE IF A BIGGER CLOUD CAN BE OBTAINED
952	! THAT SOURCE AIR...
953	!
954	! IF(CLDHGT.LT.4.E3.OR.ABE.LT.1.)THEN
955	IF(CLDHGT.LT.3.E3.OR.ABE.LT.1.)THEN
956	DO 70 NK=K,LTOP
957	UMF(NK)=0.
958	UDR(NK)=0.
959	UER(NK)=0.
960	DETLQ(NK)=0.
961	DETIC(NK)=0.
962	PPTLIQ(NK)=0.
963	70 PPTICE(NK)=0.
964	GOTO 25
965	ENDIF
966	!
967	!...IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS
968	! FLUX THIS LEVEL...
969	!
970	IF(LET.EQ.LTOP)THEN
971	UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
972	DETLQ(LTOP)=QLIQ(LTOP)UDR(LTOP)UPNEW/UPOLD
973	DETIC(LTOP)=QICE(LTOP)UDR(LTOP)UPNEW/UPOLD
974	TRPPT=TRPPT-(PPTLIQ(LTOP)+PPTICE(LTOP))
975	UER(LTOP)=0.
976	UMF(LTOP)=0.
977	GOTO 85
978	ENDIF
979	!
980	! BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET...
981	!
982	DPTT=0.
983	DO 71 NJ=LET+1,LTOP
984	71 DPTT=DPTT+DP(NJ)
985	DUMFDP=UMF(LET)/DPTT
986	!
987	!...ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
988	! RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
989	! PTOP
990	!
991	DO 75 NK=LET+1,LTOP
992	UDR(NK)=DP(NK)*DUMFDP
993	UMF(NK)=UMF(NK-1)-UDR(NK)
994	DETLQ(NK)=QLIQ(NK)*UDR(NK)
995	DETIC(NK)=QICE(NK)*UDR(NK)
996	TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
997	PPTLIQ(NK)=(UMF(NK-1)-UDR(NK))*QLQOUT(NK)
998	PPTICE(NK)=(UMF(NK-1)-UDR(NK))*QICOUT(NK)
999	TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
1000	75 CONTINUE
1001	!
1002	!...SEND UPDRAFT CHARACTERISTICS TO OUTPUT FILES...
1003	!
1004	85 CONTINUE
1005	!
1006	!...EXTEND THE UPDRAFT MASS FLUX PROFILE DOWN TO THE SOURCE LAYER FOR
1007	! THE UPDRAFT AIR...ALSO, DEFINE THETAE FOR LEVELS BELOW THE LCL...
1008	!
1009	DO 90 NK=1,K
1010	IF(NK.GE.LC)THEN
1011	IF(NK.EQ.LC)THEN
1012	UMF(NK)=VMFLCL*DP(NK)/DPTHMX
1013	UER(NK)=VMFLCL*DP(NK)/DPTHMX
1014	ELSEIF(NK.LE.KPBL)THEN
1015	UER(NK)=VMFLCL*DP(NK)/DPTHMX
1016	UMF(NK)=UMF(NK-1)+UER(NK)
1017	ELSE
1018	UMF(NK)=VMFLCL
1019	UER(NK)=0.
1020	ENDIF
1021	TU(NK)=TMIX+(Z0(NK)-ZMIX)*GDRY
1022	QU(NK)=QMIX
1023	WU(NK)=WLCL
1024	ELSE
1025	TU(NK)=0.
1026	QU(NK)=0.
1027	UMF(NK)=0.
1028	WU(NK)=0.
1029	UER(NK)=0.
1030	ENDIF
1031	UDR(NK)=0.
1032	QDT(NK)=0.
1033	QLIQ(NK)=0.
1034	QICE(NK)=0.
1035	QLQOUT(NK)=0.
1036	QICOUT(NK)=0.
1037	PPTLIQ(NK)=0.
1038	PPTICE(NK)=0.
1039	DETLQ(NK)=0.
1040	DETIC(NK)=0.
1041	RATIO2(NK)=0.
1042	EE=Q0(NK)*P0(NK)/(EP2+Q0(NK))
1043	TLOG=ALOG(EE/ALIQ)
1044	TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1045	TSAT=TDPT-(.212+1.571E-3(TDPT-T00)-4.36E-4(T0(NK)-T00))*( &
1046	T0(NK)-TDPT)
1047	THTA=T0(NK)(1.E5/P0(NK))(0.2854(1.-0.28*Q0(NK)))
1048	THETEE(NK)=THTA* &
1049	EXP((3374.6525/TSAT-2.5403)Q0(NK)(1.+0.81*Q0(NK)) &
1050	)
1051	THTES(NK)=THTA* &
1052	EXP((3374.6525/T0(NK)-2.5403)QES(NK)(1.+0.81* &
1053	QES(NK)))
1054	EQFRC(NK)=1.0
1055	90 CONTINUE
1056	!
1057	LTOP1=LTOP+1
1058	LTOPM1=LTOP-1
1059	!
1060	!...DEFINE VARIABLES ABOVE CLOUD TOP...
1061	!
1062	DO 95 NK=LTOP1,KX
1063	UMF(NK)=0.
1064	UDR(NK)=0.
1065	UER(NK)=0.
1066	QDT(NK)=0.
1067	QLIQ(NK)=0.
1068	QICE(NK)=0.
1069	QLQOUT(NK)=0.
1070	QICOUT(NK)=0.
1071	DETLQ(NK)=0.
1072	DETIC(NK)=0.
1073	PPTLIQ(NK)=0.
1074	PPTICE(NK)=0.
1075	IF(NK.GT.LTOP1)THEN
1076	TU(NK)=0.
1077	QU(NK)=0.
1078	WU(NK)=0.
1079	ENDIF
1080	THTA0(NK)=0.
1081	THTAU(NK)=0.
1082	EMS(NK)=DP(NK)*DXSQ/G
1083	EMSD(NK)=1./EMS(NK)
1084	TG(NK)=T0(NK)
1085	QG(NK)=Q0(NK)
1086	QLG(NK)=0.
1087	QIG(NK)=0.
1088	QRG(NK)=0.
1089	QSG(NK)=0.
1090	95 OMG(NK)=0.
1091	OMG(KL+1)=0.
1092	P150=P0(KLCL)-1.50E4
1093	DO 100 NK=1,LTOP
1094	THTAD(NK)=0.
1095	EMS(NK)=DP(NK)*DXSQ/G
1096	EMSD(NK)=1./EMS(NK)
1097	!
1098	!...INITIALIZE SOME VARIABLES TO BE USED LATER IN THE VERT ADVECTION
1099	! SCHEME
1100	!
1101	EXN(NK)=(P00/P0(NK))*(0.2854(1.-0.28*QDT(NK)))
1102	THTAU(NK)=TU(NK)*EXN(NK)
1103	EXN(NK)=(P00/P0(NK))*(0.2854(1.-0.28*Q0(NK)))
1104	THTA0(NK)=T0(NK)*EXN(NK)
1105	!
1106	!...LVF IS THE LEVEL AT WHICH MOISTURE FLUX IS ESTIMATED AS THE BASIS
1107	!...FOR PRECIPITATION EFFICIENCY CALCULATIONS...
1108	!
1109	IF(P0(NK).GT.P150)LVF=NK
1110	100 OMG(NK)=0.
1111	LVF=MIN0(LVF,LET)
1112	USR=UMF(LVF+1)*(QU(LVF+1)+QLIQ(LVF+1)+QICE(LVF+1))
1113	USR=AMIN1(USR,TRPPT)
1114	IF(USR.LT.1.E-8)USR=TRPPT
1115	!
1116	! WRITE(98,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
1117	! * TMIX-T00,PMIX,QMIX,ABE
1118	! WRITE(98,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
1119	! * WLCL,CLDHGT
1120	!
1121	!...COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
1122	!...AND MIDTROPOSPHERE IS USED.
1123	!
1124	WSPD(KLCL)=SQRT(U0(KLCL)U0(KLCL)+V0(KLCL)V0(KLCL))
1125	WSPD(L5)=SQRT(U0(L5)U0(L5)+V0(L5)V0(L5))
1126	WSPD(LTOP)=SQRT(U0(LTOP)U0(LTOP)+V0(LTOP)V0(LTOP))
1127	VCONV=.5*(WSPD(KLCL)+WSPD(L5))
1128	if (VCONV .gt. 0.) then
1129	TIMEC=DX/VCONV
1130	else
1131	TIMEC=3600.
1132	endif
1133	! TIMEC=DX/VCONV
1134	TADVEC=TIMEC
1135	TIMEC=AMAX1(1800.,TIMEC)
1136	TIMEC=AMIN1(3600.,TIMEC)
1137	NIC=NINT(TIMEC/DT)
1138	TIMEC=FLOAT(NIC)*DT
1139	!
1140	!...COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
1141	!
1142	! SHSIGN = CVMGT(1.,-1.,WSPD(LTOP).GT.WSPD(KLCL))
1143	IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
1144	SHSIGN=1.
1145	ELSE
1146	SHSIGN=-1.
1147	ENDIF
1148	VWS=(U0(LTOP)-U0(KLCL))(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL)) &
1149	(V0(LTOP)-V0(KLCL))
1150	VWS=1.E3SHSIGNSQRT(VWS)/(Z0(LTOP)-Z0(LCL))
1151	PEF=1.591+VWS(-.639+VWS(9.53E-2-VWS*4.96E-3))
1152	PEF=AMAX1(PEF,.2)
1153	PEF=AMIN1(PEF,.9)
1154	!
1155	!...PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE.
1156	!
1157	CBH=(ZLCL-Z0(1))*3.281E-3
1158	IF(CBH.LT.3.)THEN
1159	RCBH=.02
1160	ELSE
1161	RCBH=.96729352+CBH(-.70034167+CBH(.162179896+CBH*(- &
1162	1.2569798E-2+CBH(4.2772E-4-CBH5.44E-6))))
1163	ENDIF
1164	IF(CBH.GT.25)RCBH=2.4
1165	PEFCBH=1./(1.+RCBH)
1166	PEFCBH=AMIN1(PEFCBH,.9)
1167	!
1168	!... MEAN PEF. IS USED TO COMPUTE RAINFALL.
1169	!
1170	PEFF=.5*(PEF+PEFCBH)
1171	PEFF2=PEFF
1172	! WRITE(98,1035)PEF,PEFCBH,LC,LET,WKL,VWS
1173	!
1174	!*****************************************************************
1175	! *
1176	! COMPUTE DOWNDRAFT PROPERTIES *
1177	! *
1178	!*****************************************************************
1179	!
1180	!...LET DOWNDRAFT ORIGINATE AT THE LEVEL OF MINIMUM SATURATION EQUIVALEN
1181	!...POTENTIAL TEMPERATURE (SEQT) IN THE CLOUD LAYER, EXTEND DOWNWARD TO
1182	!...SURFACE, OR TO THE LAYER BELOW CLOUD BASE AT WHICH ENVIR SEQT IS LES
1183	!...THAN MIN SEQT IN THE CLOUD LAYER...LET DOWNDRAFT DETRAIN OVER A LAYE
1184	!...OF SPECIFIED PRESSURE-DEPTH (DPDD)...
1185	!
1186	TDER=0.
1187	KSTART=MAX0(KPBL,KLCL)
1188	THTMIN=THTES(KSTART+1)
1189	KMIN=KSTART+1
1190	DO 104 NK=KSTART+2,LTOP-1
1191	THTMIN=AMIN1(THTMIN,THTES(NK))
1192	IF(THTMIN.EQ.THTES(NK))KMIN=NK
1193	104 CONTINUE
1194	LFS=KMIN
1195	IF(RATIO2(LFS).GT.0.)CALL ENVIRTHT(P0(LFS),T0(LFS),Q0(LFS), &
1196	THETEE(LFS),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1197	EQFRC(LFS)=(THTES(LFS)-THETEU(LFS))/(THETEE(LFS)-THETEU(LFS))
1198	EQFRC(LFS)=AMAX1(EQFRC(LFS),0.)
1199	EQFRC(LFS)=AMIN1(EQFRC(LFS),1.)
1200	THETED(LFS)=THTES(LFS)
1201	!
1202	!...ESTIMATE THE EFFECT OF MELTING PRECIPITATION IN THE DOWNDRAFT...
1203	!
1204	IF(ML.GT.0)THEN
1205	DTMLTD=0.5(QU(KLCL)-QU(LTOP))RLF/CP
1206	ELSE
1207	DTMLTD=0.
1208	ENDIF
1209	TZ(LFS)=T0(LFS)-DTMLTD
1210	ES=ALIQEXP((TZ(LFS)BLIQ-CLIQ)/(TZ(LFS)-DLIQ))
1211	QS=EP2*ES/(P0(LFS)-ES)
1212	QD(LFS)=EQFRC(LFS)Q0(LFS)+(1.-EQFRC(LFS))QU(LFS)
1213	THTAD(LFS)=TZ(LFS)(P00/P0(LFS))(0.2854(1.-0.28*QD(LFS)))
1214	IF(QD(LFS).GE.QS)THEN
1215	THETED(LFS)=THTAD(LFS)* &
1216	EXP((3374.6525/TZ(LFS)-2.5403)QS(1.+0.81*QS))
1217	ELSE
1218	CALL ENVIRTHT(P0(LFS),TZ(LFS),QD(LFS),THETED(LFS),0.,RL,EP2,ALIQ, &
1219	BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1220	ENDIF
1221	DO 107 NK=1,LFS
1222	ND=LFS-NK
1223	IF(THETED(LFS).GT.THTES(ND).OR.ND.EQ.1)THEN
1224	LDB=ND
1225	!
1226	!...IF DOWNDRAFT NEVER BECOMES NEGATIVELY BUOYANT OR IF IT
1227	!...IS SHALLOWER 50 mb, DON'T ALLOW IT TO OCCUR AT ALL...
1228	!
1229	IF(NK.EQ.1.OR.(P0(LDB)-P0(LFS)).LT.50.E2)GOTO 141
1230	! testing ---- no downdraft
1231	! GOTO 141
1232	GOTO 110
1233	ENDIF
1234	107 CONTINUE
1235	!
1236	!...ALLOW DOWNDRAFT TO DETRAIN IN A SINGLE LAYER, BUT WITH DOWNDRAFT AIR
1237	!...TYPICALLY FLUSHED UP INTO HIGHER LAYERS AS ALLOWED IN THE TOTAL
1238	!...VERTICAL ADVECTION CALCULATIONS FARTHER DOWN IN THE CODE...
1239	!
1240	110 DPDD=DP(LDB)
1241	LDT=LDB
1242	FRC=1.
1243	DPT=0.
1244	! DO 115 NK=LDB,LFS
1245	! DPT=DPT+DP(NK)
1246	! IF(DPT.GT.DPDD)THEN
1247	! LDT=NK
1248	! FRC=(DPDD+DP(NK)-DPT)/DP(NK)
1249	! GOTO 120
1250	! ENDIF
1251	! IF(NK.EQ.LFS-1)THEN
1252	! LDT=NK
1253	! FRC=1.
1254	! DPDD=DPT
1255	! GOTO 120
1256	! ENDIF
1257	!115 CONTINUE
1258	120 CONTINUE
1259	!
1260	!...TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX..
1261	!
1262	TVD(LFS)=T0(LFS)(1.+0.608QES(LFS))
1263	RDD=P0(LFS)/(R*TVD(LFS))
1264	A1=(1.-PEFF)*AU0
1265	DMF(LFS)=-A1*RDD
1266	DER(LFS)=EQFRC(LFS)*DMF(LFS)
1267	DDR(LFS)=0.
1268	DO 140 ND=LFS-1,LDB,-1
1269	ND1=ND+1
1270	IF(ND.LE.LDT)THEN
1271	DER(ND)=0.
1272	DDR(ND)=-DMF(LDT+1)DP(ND)FRC/DPDD
1273	DMF(ND)=DMF(ND1)+DDR(ND)
1274	FRC=1.
1275	THETED(ND)=THETED(ND1)
1276	QD(ND)=QD(ND1)
1277	ELSE
1278	DER(ND)=DMF(LFS)0.03DP(ND)/RAD
1279	DDR(ND)=0.
1280	DMF(ND)=DMF(ND1)+DER(ND)
1281	IF(RATIO2(ND).GT.0.)CALL ENVIRTHT(P0(ND),T0(ND),Q0(ND), &
1282	THETEE(ND),0.,RL,EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1283	THETED(ND)=(THETED(ND1)DMF(ND1)+THETEE(ND)DER(ND))/DMF(ND)
1284	QD(ND)=(QD(ND1)DMF(ND1)+Q0(ND)DER(ND))/DMF(ND)
1285	ENDIF
1286	140 CONTINUE
1287	TDER=0.
1288	!
1289	!...CALCULATION AN EVAPORATION RATE FOR GIVEN MASS FLUX...
1290	!
1291	DO 135 ND=LDB,LDT
1292	TZ(ND)= &
1293	TPDD(P0(ND),THETED(LDT),T0(ND),QS,QD(ND),1.0,XLV0,XLV1, &
1294	EP2,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
1295	ES=ALIQEXP((TZ(ND)BLIQ-CLIQ)/(TZ(ND)-DLIQ))
1296	QS=EP2*ES/(P0(ND)-ES)
1297	DSSDT=(CLIQ-BLIQDLIQ)/((TZ(ND)-DLIQ)(TZ(ND)-DLIQ))
1298	RL=XLV0-XLV1*TZ(ND)
1299	DTMP=RLQS(1.-RHBC)/(CP+RLRHBCQS*DSSDT)
1300	T1RH=TZ(ND)+DTMP
1301	ES=RHBCALIQEXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
1302	QSRH=EP2*ES/(P0(ND)-ES)
1303	!
1304	!...CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN ACTUAL
1305	!...MIXING RATIO...IF SO, ADJUST TO GIVE ZERO EVAPORATION...
1306	!
1307	IF(QSRH.LT.QD(ND))THEN
1308	QSRH=QD(ND)
1309	! T1RH=T1+(QS-QSRH)*RL/CP
1310	T1RH=TZ(ND)
1311	ENDIF
1312	TZ(ND)=T1RH
1313	QS=QSRH
1314	TDER=TDER+(QS-QD(ND))*DDR(ND)
1315	QD(ND)=QS
1316	135 THTAD(ND)=TZ(ND)(P00/P0(ND))(0.2854(1.-0.28*QD(ND)))
1317	!
1318	!...IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
1319	!...HUMIDITY, NO DOWNDRAFT IS ALLOWED...
1320	!
1321	141 IF(TDER.LT.1.)THEN
1322	! WRITE(98,3004)I,J
1323	3004 FORMAT(' ','I=',I3,2X,'J=',I3)
1324	PPTFLX=TRPPT
1325	CPR=TRPPT
1326	TDER=0.
1327	CNDTNF=0.
1328	UPDINC=1.
1329	LDB=LFS
1330	DO 117 NDK=1,LTOP
1331	DMF(NDK)=0.
1332	DER(NDK)=0.
1333	DDR(NDK)=0.
1334	THTAD(NDK)=0.
1335	WD(NDK)=0.
1336	TZ(NDK)=0.
1337	117 QD(NDK)=0.
1338	AINCM2=100.
1339	GOTO 165
1340	ENDIF
1341	!
1342	!...ADJUST DOWNDRAFT MASS FLUX SO THAT EVAPORATION RATE IN DOWNDRAFT IS
1343	!...CONSISTENT WITH PRECIPITATION EFFICIENCY RELATIONSHIP...
1344	!
1345	DEVDMF=TDER/DMF(LFS)
1346	PPR=0.
1347	PPTFLX=PEFF*USR
1348	RCED=TRPPT-PPTFLX
1349	!
1350	!...PPR IS THE TOTAL AMOUNT OF PRECIPITATION THAT FALLS OUT OF THE
1351	!...UPDRAFT FROM CLOUD BASE TO THE LFS...UPDRAFT MASS FLUX WILL BE
1352	!...INCREASED UP TO THE LFS TO ACCOUNT FOR UPDRAFT AIR MIXING WITH
1353	!...ENVIRONMENTAL AIR TO THE UPDRAFT, SO PPR WILL INCREASE
1354	!...PROPORTIONATELY...
1355	!
1356	DO 132 NM=KLCL,LFS
1357	132 PPR=PPR+PPTLIQ(NM)+PPTICE(NM)
1358	IF(LFS.GE.KLCL)THEN
1359	DPPTDF=(1.-PEFF)PPR(1.-EQFRC(LFS))/UMF(LFS)
1360	ELSE
1361	DPPTDF=0.
1362	ENDIF
1363	!
1364	!...CNDTNF IS THE AMOUNT OF CONDENSATE TRANSFERRED ALONG WITH UPDRAFT
1365	!...MASS THE DOWNDRAFT AT THE LFS...
1366	!
1367	CNDTNF=(QLIQ(LFS)+QICE(LFS))*(1.-EQFRC(LFS))
1368	DMFLFS=RCED/(DEVDMF+DPPTDF+CNDTNF)
1369	IF(DMFLFS.GT.0.)THEN
1370	TDER=0.
1371	GOTO 141
1372	ENDIF
1373	!
1374	!...DDINC IS THE FACTOR BY WHICH TO INCREASE THE FIRST-GUESS DOWNDRAFT
1375	!...MASS FLUX TO SATISFY THE PRECIP EFFICIENCY RELATIONSHIP, UPDINC IS T
1376	!...WHICH TO INCREASE THE UPDRAFT MASS FLUX BELOW THE LFS TO ACCOUNT FOR
1377	!...TRANSFER OF MASS FROM UPDRAFT TO DOWNDRAFT...
1378	!
1379	! DDINC=DMFLFS/DMF(LFS)
1380	IF(LFS.GE.KLCL)THEN
1381	UPDINC=(UMF(LFS)-(1.-EQFRC(LFS))*DMFLFS)/UMF(LFS)
1382	!
1383	!...LIMIT UPDINC TO LESS THAN OR EQUAL TO 1.5...
1384	!
1385	IF(UPDINC.GT.1.5)THEN
1386	UPDINC=1.5
1387	DMFLFS2=UMF(LFS)*(UPDINC-1.)/(EQFRC(LFS)-1.)
1388	RCED2=DMFLFS2*(DEVDMF+DPPTDF+CNDTNF)
1389	PPTFLX=PPTFLX+(RCED-RCED2)
1390	PEFF2=PPTFLX/USR
1391	RCED=RCED2
1392	DMFLFS=DMFLFS2
1393	ENDIF
1394	ELSE
1395	UPDINC=1.
1396	ENDIF
1397	DDINC=DMFLFS/DMF(LFS)
1398	DO 149 NK=LDB,LFS
1399	DMF(NK)=DMF(NK)*DDINC
1400	DER(NK)=DER(NK)*DDINC
1401	DDR(NK)=DDR(NK)*DDINC
1402	149 CONTINUE
1403	CPR=TRPPT+PPR*(UPDINC-1.)
1404	PPTFLX=PPTFLX+PEFFPPR(UPDINC-1.)
1405	PEFF=PEFF2
1406	TDER=TDER*DDINC
1407	!
1408	!...ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER AN
1409	! DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE ESTIMATE
1410	! FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS...
1411	!
1412	DO 155 NK=LC,LFS
1413	UMF(NK)=UMF(NK)*UPDINC
1414	UDR(NK)=UDR(NK)*UPDINC
1415	UER(NK)=UER(NK)*UPDINC
1416	PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
1417	PPTICE(NK)=PPTICE(NK)*UPDINC
1418	DETLQ(NK)=DETLQ(NK)*UPDINC
1419	155 DETIC(NK)=DETIC(NK)*UPDINC
1420	!
1421	!...ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW THE
1422	!...DOWNDRAFT...
1423	!
1424	IF(LDB.GT.1)THEN
1425	DO 156 NK=1,LDB-1
1426	DMF(NK)=0.
1427	DER(NK)=0.
1428	DDR(NK)=0.
1429	WD(NK)=0.
1430	TZ(NK)=0.
1431	QD(NK)=0.
1432	THTAD(NK)=0.
1433	156 CONTINUE
1434	ENDIF
1435	DO 157 NK=LFS+1,KX
1436	DMF(NK)=0.
1437	DER(NK)=0.
1438	DDR(NK)=0.
1439	WD(NK)=0.
1440	TZ(NK)=0.
1441	QD(NK)=0.
1442	THTAD(NK)=0.
1443	157 CONTINUE
1444	DO 158 NK=LDT+1,LFS-1
1445	TZ(NK)=0.
1446	QD(NK)=0.
1447	158 CONTINUE
1448	!
1449	!
1450	!...SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE
1451	! INFLOW INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN
1452	! IS AVAILABLE IN THAT LAYER INITIALLY...
1453	!
1454	165 AINCMX=1000.
1455	LMAX=MAX0(KLCL,LFS)
1456	DO 166 NK=LC,LMAX
1457	IF((UER(NK)-DER(NK)).GT.0.)AINCM1=EMS(NK)/((UER(NK)-DER(NK))* &
1458	TIMEC)
1459	AINCMX=AMIN1(AINCMX,AINCM1)
1460	166 CONTINUE
1461	AINC=1.
1462	IF(AINCMX.LT.AINC)AINC=AINCMX
1463	!
1464	!...SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRFT AND DOWNDRFT...THEY
1465	!...WILL ITERATIVELY ADJUSTED BY THE FACTOR AINC TO SATISFY THE
1466	!...STABILIZATION CLOSURE...
1467	!
1468	NCOUNT=0
1469	TDER2=TDER
1470	PPTFL2=PPTFLX
1471	DO 170 NK=1,LTOP
1472	DETLQ2(NK)=DETLQ(NK)
1473	DETIC2(NK)=DETIC(NK)
1474	UDR2(NK)=UDR(NK)
1475	UER2(NK)=UER(NK)
1476	DDR2(NK)=DDR(NK)
1477	DER2(NK)=DER(NK)
1478	UMF2(NK)=UMF(NK)
1479	DMF2(NK)=DMF(NK)
1480	170 CONTINUE
1481	FABE=1.
1482	STAB=0.95
1483	NOITR=0
1484	IF(AINC/AINCMX.GT.0.999)THEN
1485	NCOUNT=0
1486	GOTO 255
1487	ENDIF
1488	ISTOP=0
1489	175 NCOUNT=NCOUNT+1
1490	!
1491	!*****************************************************************
1492	! *
1493	! COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE *
1494	! *
1495	!*****************************************************************
1496	!
1497	!...DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER TO
1498	!...SATISFY MASS CONTINUITY...
1499	!
1500	185 CONTINUE
1501	DTT=TIMEC
1502	DO 200 NK=1,LTOP
1503	DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
1504	IF(NK.GT.1)THEN
1505	OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
1506	DTT1=0.75*DP(NK-1)/(ABS(OMG(NK))+1.E-10)
1507	DTT=AMIN1(DTT,DTT1)
1508	ENDIF
1509	200 CONTINUE
1510	DO 488 NK=1,LTOP
1511	THPA(NK)=THTA0(NK)
1512	QPA(NK)=Q0(NK)
1513	NSTEP=NINT(TIMEC/DTT+1)
1514	DTIME=TIMEC/FLOAT(NSTEP)
1515	FXM(NK)=OMG(NK)*DXSQ/G
1516	488 CONTINUE
1517	!
1518	!...DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV...
1519	!
1520	DO 495 NTC=1,NSTEP
1521	!
1522	!...ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED
1523	!...SIGN OF OMEGA...
1524	!
1525	DO 493 NK=1,LTOP
1526	THFXTOP(NK)=0.
1527	THFXBOT(NK)=0.
1528	QFXTOP(NK)=0.
1529	493 QFXBOT(NK)=0.
1530	DO 494 NK=2,LTOP
1531	IF(OMG(NK).LE.0.)THEN
1532	THFXBOT(NK)=-FXM(NK)*THPA(NK-1)
1533	QFXBOT(NK)=-FXM(NK)*QPA(NK-1)
1534	THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
1535	QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
1536	ELSE
1537	THFXBOT(NK)=-FXM(NK)*THPA(NK)
1538	QFXBOT(NK)=-FXM(NK)*QPA(NK)
1539	THFXTOP(NK-1)=THFXTOP(NK-1)-THFXBOT(NK)
1540	QFXTOP(NK-1)=QFXTOP(NK-1)-QFXBOT(NK)
1541	ENDIF
1542	494 CONTINUE
1543	!
1544	!...UPDATE THE THETA AND QV VALUES AT EACH LEVEL..
1545	!
1546	DO 492 NK=1,LTOP
1547	THPA(NK)=THPA(NK)+(THFXBOT(NK)+UDR(NK)THTAU(NK)+DDR(NK) &
1548	THTAD(NK)+THFXTOP(NK)-(UER(NK)-DER(NK))THTA0(NK)) &
1549	DTIME*EMSD(NK)
1550	QPA(NK)=QPA(NK)+(QFXBOT(NK)+UDR(NK)QDT(NK)+DDR(NK)QD(NK)+ &
1551	QFXTOP(NK)-(UER(NK)-DER(NK))Q0(NK))DTIME*EMSD(NK)
1552
1553	492 CONTINUE
1554	495 CONTINUE
1555	DO 498 NK=1,LTOP
1556	THTAG(NK)=THPA(NK)
1557	QG(NK)=QPA(NK)
1558	498 CONTINUE
1559	!
1560	!...CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE; IF SO,
1561	!...BORROW MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE ZERO.
1562	!
1563	DO 499 NK=1,LTOP
1564	IF(QG(NK).LT.0.)THEN
1565	IF(NK.EQ.1)THEN
1566	CALL wrf_error_fatal ( 'module_cu_kf.F: problem with kf scheme: qg = 0 at the surface' )
1567	ENDIF
1568	NK1=NK+1
1569	IF(NK.EQ.LTOP)NK1=KLCL
1570	TMA=QG(NK1)*EMS(NK1)
1571	TMB=QG(NK-1)*EMS(NK-1)
1572	TMM=(QG(NK)-1.E-9)*EMS(NK)
1573	BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
1574	ACOEFF=BCOEFF*TMA/TMB
1575	TMB=TMB*(1.-BCOEFF)
1576	TMA=TMA*(1.-ACOEFF)
1577	IF(NK.EQ.LTOP)THEN
1578	QVDIFF=(QG(NK1)-TMAEMSD(NK1))100./QG(NK1)
1579	IF(ABS(QVDIFF).GT.1.)THEN
1580	PRINT *,'--WARNING-- CLOUD BASE WATER VAPOR CHANGES BY ', &
1581	QVDIFF, &
1582	' PERCENT WHEN MOISTURE IS BORROWED TO PREVENT NEG VALUES', &
1583	' IN KAIN-FRITSCH'
1584	ENDIF
1585	ENDIF
1586	QG(NK)=1.E-9
1587	QG(NK1)=TMA*EMSD(NK1)
1588	QG(NK-1)=TMB*EMSD(NK-1)
1589	ENDIF
1590	499 CONTINUE
1591	TOPOMG=(UDR(LTOP)-UER(LTOP))DP(LTOP)EMSD(LTOP)
1592	IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
1593	! WRITE(98,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;'
1594	! * ,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1595	WRITE(6,*)'ERROR: MASS DOES NOT BALANCE IN KF SCHEME;' &
1596	,'TOPOMG, OMG =',TOPOMG,OMG(LTOP)
1597	ISTOP=1
1598	GOTO 265
1599	ENDIF
1600	!
1601	!...CONVERT THETA TO T...
1602	!
1603	! PAY ATTENTION ...
1604	!
1605	DO 230 NK=1,LTOP
1606	EXN(NK)=(P00/P0(NK))*(0.2854(1.-0.28*QG(NK)))
1607	TG(NK)=THTAG(NK)/EXN(NK)
1608	TVG(NK)=TG(NK)(1.+0.608QG(NK))
1609	230 CONTINUE
1610	!
1611	!*******************************************************************
1612	! *
1613	! COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY. *
1614	! *
1615	!*******************************************************************
1616	!
1617	!...THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
1618	!
1619	THMIX=0.
1620	QMIX=0.
1621	PMIX=0.
1622	DO 217 NK=LC,KPBL
1623	ROCPQ=0.2854(1.-0.28QG(NK))
1624	THMIX=THMIX+DP(NK)TG(NK)(P00/P0(NK))**ROCPQ
1625	QMIX=QMIX+DP(NK)*QG(NK)
1626	217 PMIX=PMIX+DP(NK)*P0(NK)
1627	THMIX=THMIX/DPTHMX
1628	QMIX=QMIX/DPTHMX
1629	PMIX=PMIX/DPTHMX
1630	ROCPQ=0.2854(1.-0.28QMIX)
1631	TMIX=THMIX(PMIX/P00)*ROCPQ
1632	ES=ALIQEXP((TMIXBLIQ-CLIQ)/(TMIX-DLIQ))
1633	QS=EP2*ES/(PMIX-ES)
1634	!
1635	!...REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY...
1636	!
1637	IF(QMIX.GT.QS)THEN
1638	RL=XLV0-XLV1*TMIX
1639	CPM=CP(1.+0.887QMIX)
1640	DSSDT=QS(CLIQ-BLIQDLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
1641	DQ=(QMIX-QS)/(1.+RL*DSSDT/CPM)
1642	TMIX=TMIX+RL/CP*DQ
1643	QMIX=QMIX-DQ
1644	ROCPQ=0.2854(1.-0.28QMIX)
1645	THMIX=TMIX(P00/PMIX)*ROCPQ
1646	TLCL=TMIX
1647	PLCL=PMIX
1648	ELSE
1649	QMIX=AMAX1(QMIX,0.)
1650	EMIX=QMIX*PMIX/(EP2+QMIX)
1651	TLOG=ALOG(EMIX/ALIQ)
1652	TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
1653	TLCL=TDPT-(.212+1.571E-3(TDPT-T00)-4.36E-4(TMIX-T00))*(TMIX- &
1654	TDPT)
1655	TLCL=AMIN1(TLCL,TMIX)
1656	CPORQ=1./ROCPQ
1657	PLCL=P00(TLCL/THMIX)*CPORQ
1658	ENDIF
1659	TVLCL=TLCL(1.+0.608QMIX)
1660	DO 235 NK=LC,KL
1661	KLCL=NK
1662	235 IF(PLCL.GE.P0(NK))GOTO 240
1663	240 K=KLCL-1
1664	DLP=ALOG(PLCL/P0(K))/ALOG(P0(KLCL)/P0(K))
1665	!
1666	!...ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL...
1667	!
1668	TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
1669	QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
1670	TVEN=TENV(1.+0.608QENV)
1671	TVBAR=0.5*(TVG(K)+TVEN)
1672	! ZLCL=Z0(K)+RTVBARALOG(P0(K)/PLCL)/G
1673	ZLCL=Z0(K)+(Z0(KLCL)-Z0(K))*DLP
1674	TVAVG=0.5(TVEN+TG(KLCL)(1.+0.608*QG(KLCL)))
1675	PLCL=P0(KLCL)EXP(G/(RTVAVG)*(Z0(KLCL)-ZLCL))
1676	THETEU(K)=TMIX(1.E5/PMIX)(0.2854(1.-0.28QMIX)) &
1677	EXP((3374.6525/TLCL-2.5403)QMIX(1.+0.81*QMIX))
1678	ES=ALIQEXP((TENVBLIQ-CLIQ)/(TENV-DLIQ))
1679	QESE=EP2*ES/(PLCL-ES)
1680	THTESG(K)=TENV(1.E5/PLCL)(0.2854(1.-0.28QESE)) &
1681	EXP((3374.6525/TENV-2.5403)QESE(1.+0.81*QESE))
1682	!
1683	!...COMPUTE ADJUSTED ABE(ABEG).
1684	!
1685	ABEG=0.
1686	THTUDL=THETEU(K)
1687	DO 245 NK=K,LTOPM1
1688	NK1=NK+1
1689	ES=ALIQEXP((TG(NK1)BLIQ-CLIQ)/(TG(NK1)-DLIQ))
1690	QESE=EP2*ES/(P0(NK1)-ES)
1691	THTESG(NK1)=TG(NK1)(1.E5/P0(NK1))(0.2854(1.-0.28QESE)) &
1692	EXP((3374.6525/TG(NK1)-2.5403)QESE(1.+0.81*QESE) &
1693	)
1694	! DZZ=CVMGT(Z0(KLCL)-ZLCL,DZA(NK),NK.EQ.K)
1695	IF(NK.EQ.K)THEN
1696	DZZ=Z0(KLCL)-ZLCL
1697	ELSE
1698	DZZ=DZA(NK)
1699	ENDIF
1700	BE=((2.THTUDL)/(THTESG(NK1)+THTESG(NK))-1.)DZZ
1701	245 IF(BE.GT.0.)ABEG=ABEG+BE*G
1702	!
1703	!...ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
1704	!...THE PERIOD TIMEC...
1705	!
1706	IF(NOITR.EQ.1)THEN
1707	! WRITE(98,1060)FABE
1708	GOTO 265
1709	ENDIF
1710	DABE=AMAX1(ABE-ABEG,0.1*ABE)
1711	FABE=ABEG/(ABE+1.E-8)
1712	IF(FABE.GT.1.)THEN
1713	! WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS '
1714	! *,'GRID POINT; NO CONVECTION ALLOWED!'
1715	GOTO 325
1716	ENDIF
1717	IF(NCOUNT.NE.1)THEN
1718	DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
1719	IF(DFDA.GT.0.)THEN
1720	NOITR=1
1721	AINC=AINCOLD
1722	GOTO 255
1723	ENDIF
1724	ENDIF
1725	AINCOLD=AINC
1726	FABEOLD=FABE
1727	IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
1728	! WRITE(98,1055)FABE
1729	GOTO 265
1730	ENDIF
1731	IF(FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB)GOTO 265
1732	IF(NCOUNT.GT.10)THEN
1733	! WRITE(98,1060)FABE
1734	GOTO 265
1735	ENDIF
1736	!
1737	!...IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE
1738	!...CONVECTIVE MASS FLUX BY THE FACTOR AINC:
1739	!
1740	IF(FABE.EQ.0.)THEN
1741	AINC=AINC*0.5
1742	ELSE
1743	AINC=AINCSTABABE/(DABE+1.E-8)
1744	ENDIF
1745	255 AINC=AMIN1(AINCMX,AINC)
1746	!...IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION
1747	!...WILL BE MINIMAL SO JUST IGNORE IT...
1748	IF(AINC.LT.0.05)GOTO 325
1749	! AINC=AMAX1(AINC,0.05)
1750	TDER=TDER2*AINC
1751	PPTFLX=PPTFL2*AINC
1752	! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABEOLD,AINCOLD
1753	DO 260 NK=1,LTOP
1754	UMF(NK)=UMF2(NK)*AINC
1755	DMF(NK)=DMF2(NK)*AINC
1756	DETLQ(NK)=DETLQ2(NK)*AINC
1757	DETIC(NK)=DETIC2(NK)*AINC
1758	UDR(NK)=UDR2(NK)*AINC
1759	UER(NK)=UER2(NK)*AINC
1760	DER(NK)=DER2(NK)*AINC
1761	DDR(NK)=DDR2(NK)*AINC
1762	260 CONTINUE
1763	!
1764	!...GO BACK UP FOR ANOTHER ITERATION...
1765	!
1766	GOTO 175
1767	265 CONTINUE
1768	!
1769	!...CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR THIS
1770	!...GRID POINT...
1771	!
1772	!...COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV...
1773	!
1774	!...FRC2 IS THE FRACTION OF TOTAL CONDENSATE
1775	!...GENERATED THAT GOES INTO PRECIPITIATION
1776	FRC2=PPTFLX/(CPR*AINC)
1777	DO 270 NK=1,LTOP
1778	QLPA(NK)=QL0(NK)
1779	QIPA(NK)=QI0(NK)
1780	QRPA(NK)=QR0(NK)
1781	QSPA(NK)=QS0(NK)
1782	RAINFB(NK)=PPTLIQ(NK)AINCFBFRC*FRC2
1783	SNOWFB(NK)=PPTICE(NK)AINCFBFRC*FRC2
1784	270 CONTINUE
1785	DO 290 NTC=1,NSTEP
1786	!
1787	!...ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH
1788	!...LAYER BASED ON THE SIGN OF OMEGA...
1789	!
1790	DO 275 NK=1,LTOP
1791	QLFXIN(NK)=0.
1792	QLFXOUT(NK)=0.
1793	QIFXIN(NK)=0.
1794	QIFXOUT(NK)=0.
1795	QRFXIN(NK)=0.
1796	QRFXOUT(NK)=0.
1797	QSFXIN(NK)=0.
1798	QSFXOUT(NK)=0.
1799	275 CONTINUE
1800	DO 280 NK=2,LTOP
1801	IF(OMG(NK).LE.0.)THEN
1802	QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
1803	QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
1804	QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
1805	QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
1806	QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
1807	QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
1808	QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
1809	QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
1810	ELSE
1811	QLFXOUT(NK)=FXM(NK)*QLPA(NK)
1812	QIFXOUT(NK)=FXM(NK)*QIPA(NK)
1813	QRFXOUT(NK)=FXM(NK)*QRPA(NK)
1814	QSFXOUT(NK)=FXM(NK)*QSPA(NK)
1815	QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
1816	QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
1817	QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
1818	QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
1819	ENDIF
1820	280 CONTINUE
1821	!
1822	!...UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL...
1823	!
1824	DO 285 NK=1,LTOP
1825	QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))DTIME &
1826	EMSD(NK)
1827	QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))DTIME &
1828	EMSD(NK)
1829	QRPA(NK)=QRPA(NK)+(QRFXIN(NK)+QLQOUT(NK)*UDR(NK)-QRFXOUT(NK) &
1830	+RAINFB(NK))DTIMEEMSD(NK)
1831	QSPA(NK)=QSPA(NK)+(QSFXIN(NK)+QICOUT(NK)*UDR(NK)-QSFXOUT(NK) &
1832	+SNOWFB(NK))DTIMEEMSD(NK)
1833	285 CONTINUE
1834	290 CONTINUE
1835	DO 295 NK=1,LTOP
1836	QLG(NK)=QLPA(NK)
1837	QIG(NK)=QIPA(NK)
1838	QRG(NK)=QRPA(NK)
1839	QSG(NK)=QSPA(NK)
1840	295 CONTINUE
1841	! WRITE(98,1080)LFS,LDB,LDT,TIMEC,NSTEP,NCOUNT,FABE,AINC
1842	!
1843	!...SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES...
1844	!
1845	IF(ISTOP.EQ.1)THEN
1846	WRITE(6,1070)' P ',' DP ',' DT K/D ',' DR K/D ',' OMG ', &
1847	' DOMGDP ',' UMF ',' UER ',' UDR ',' DMF ',' DER ' &
1848	,' DDR ',' EMS ',' W0 ',' DETLQ ',' DETIC '
1849	DO 300 K=LTOP,1,-1
1850	DTT=(TG(K)-T0(K))*86400./TIMEC
1851	RL=XLV0-XLV1*TG(K)
1852	DR=-(QG(K)-Q0(K))RL86400./(TIMEC*CP)
1853	UDFRC=UDR(K)TIMECEMSD(K)
1854	UEFRC=UER(K)TIMECEMSD(K)
1855	DDFRC=DDR(K)TIMECEMSD(K)
1856	DEFRC=-DER(K)TIMECEMSD(K)
1857	WRITE (6,1075)P0(K)/100.,DP(K)/100.,DTT,DR,OMG(K),DOMGDP(K)* &
1858	1.E4,UMF(K)/1.E6,UEFRC,UDFRC,DMF(K)/1.E6,DEFRC &
1859	,DDFRC,EMS(K)/1.E11,W0AVG1D(K)*1.E2,DETLQ(K) &
1860	TIMECEMSD(K)1.E3,DETIC(K)TIMECEMSD(K) &
1861	1.E3
1862	300 CONTINUE
1863	WRITE(6,1085)'K','P','Z','T0','TG','DT','TU','TD','Q0','QG', &
1864	'DQ','QU','QD','QLG','QIG','QRG','QSG','RH0','RHG'
1865	DO 305 K=KX,1,-1
1866	DTT=TG(K)-T0(K)
1867	TUC=TU(K)-T00
1868	IF(K.LT.LC.OR.K.GT.LTOP)TUC=0.
1869	TDC=TZ(K)-T00
1870	IF((K.LT.LDB.OR.K.GT.LDT).AND.K.NE.LFS)TDC=0.
1871	ES=ALIQEXP((BLIQTG(K)-CLIQ)/(TG(K)-DLIQ))
1872	QGS=ES*EP2/(P0(K)-ES)
1873	RH0=Q0(K)/QES(K)
1874	RHG=QG(K)/QGS
1875	WRITE (6,1090)K,P0(K)/100.,Z0(K),T0(K)-T00,TG(K)-T00,DTT,TUC &
1876	,TDC,Q0(K)1000.,QG(K)1000.,(QG(K)-Q0(K))* &
1877	1000.,QU(K)1000.,QD(K)1000.,QLG(K)*1000., &
1878	QIG(K)1000.,QRG(K)1000.,QSG(K)*1000.,RH0,RHG
1879	305 CONTINUE
1880	!
1881	!...IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT A
1882	!...TO BE USED LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN...
1883	!
1884	IF(ISTOP.EQ.1)THEN
1885	DO 310 K=1,KX
1886	WRITE ( wrf_err_message , 1115 ) &
1887	Z0(K),P0(K)/100.,T0(K)-273.16,Q0(K)*1000., &
1888	U0(K),V0(K),DP(K)/100.,W0AVG1D(K)
1889	CALL wrf_message ( TRIM( wrf_err_message ) )
1890	310 CONTINUE
1891	CALL wrf_error_fatal ( 'module_cu_kf.F: KAIN-FRITSCH' )
1892	ENDIF
1893	ENDIF
1894	CNDTNF=(1.-EQFRC(LFS))(QLIQ(LFS)+QICE(LFS))DMF(LFS)
1895	! WRITE(98,1095)CPR*AINC,TDER+PPTFLX+CNDTNF
1896	!
1897	! EVALUATE MOISTURE BUDGET...
1898	!
1899	QINIT=0.
1900	QFNL=0.
1901	DPT=0.
1902	DO 315 NK=1,LTOP
1903	DPT=DPT+DP(NK)
1904	QINIT=QINIT+Q0(NK)*EMS(NK)
1905	QFNL=QFNL+QG(NK)*EMS(NK)
1906	QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
1907	315 CONTINUE
1908	QFNL=QFNL+PPTFLXTIMEC(1.-FBFRC)
1909	ERR2=(QFNL-QINIT)*100./QINIT
1910	! WRITE(98,1110)QINIT,QFNL,ERR2
1911	! IF(ABS(ERR2).GT.0.05)STOP 'QVERR'
1912	IF(ABS(ERR2).GT.0.05)CALL wrf_error_fatal( 'module_cu_kf.F: QVERR' )
1913	RELERR=ERR2QINIT/(PPTFLXTIMEC+1.E-10)
1914	! WRITE(98,1120)RELERR
1915	! WRITE(98,*)'TDER, CPR, USR, TRPPT =',
1916	! TDER,CPRAINC,USRAINC,TRPPTAINC
1917	!
1918	!...FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
1919	!
1920	!...IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
1921	!...TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC...
1922	!
1923	IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/DT)
1924	! NCA(I,J)=FLOAT(NIC)
1925	NCA(I,J)=TADVEC
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