Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_cu_kf.F @ 3026

Last change on this file since 3026 was 2759, checked in by aslmd, 3 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 155.9 KB

Rev	Line
[2759]	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