Context Navigation

module_cu_kf.F

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

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

Download in other formats:

Original Format