Context Navigation

module_cu_kfeta.F

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

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

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

Download in other formats:

Original Format