Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_mp_etanew.F @ 3094

Last change on this file since 3094 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: 107.5 KB

Rev	Line
[2759]	1	!WRF:MODEL_MP:PHYSICS
	2	!
	3	MODULE module_mp_etanew
	4	!
	5	!-----------------------------------------------------------------------
	6	REAL,PRIVATE,SAVE :: ABFR, CBFR, CIACW, CIACR, C_N0r0, &
	7	& CN0r0, CN0r_DMRmin, CN0r_DMRmax, CRACW, CRAUT, ESW0, &
	8	& RFmax, RQR_DR1, RQR_DR2, RQR_DR3, RQR_DRmin, &
	9	& RQR_DRmax, RR_DRmin, RR_DR1, RR_DR2, RR_DR3, RR_DRmax
	10	!
	11	INTEGER, PRIVATE,PARAMETER :: MY_T1=1, MY_T2=35
	12	REAL,PRIVATE,DIMENSION(MY_T1:MY_T2),SAVE :: MY_GROWTH
	13	!
	14	REAL, PRIVATE,PARAMETER :: DMImin=.05e-3, DMImax=1.e-3, &
	15	& DelDMI=1.e-6,XMImin=1.e6*DMImin
	16	INTEGER, PUBLIC,PARAMETER :: XMImax=1.e6*DMImax, XMIexp=.0536, &
	17	& MDImin=XMImin, MDImax=XMImax
	18	REAL, PRIVATE,DIMENSION(MDImin:MDImax) :: &
	19	& ACCRI,SDENS,VSNOWI,VENTI1,VENTI2
	20	!
	21	REAL, PRIVATE,PARAMETER :: DMRmin=.05e-3, DMRmax=.45e-3, &
	22	& DelDMR=1.e-6,XMRmin=1.e6DMRmin, XMRmax=1.e6DMRmax
	23	INTEGER, PRIVATE,PARAMETER :: MDRmin=XMRmin, MDRmax=XMRmax
	24	REAL, PRIVATE,DIMENSION(MDRmin:MDRmax):: &
	25	& ACCRR,MASSR,RRATE,VRAIN,VENTR1,VENTR2
	26	!
	27	INTEGER, PRIVATE,PARAMETER :: Nrime=40
	28	REAL, DIMENSION(2:9,0:Nrime),PRIVATE,SAVE :: VEL_RF
	29	!
	30	INTEGER,PARAMETER :: NX=7501
	31	REAL, PARAMETER :: XMIN=180.0,XMAX=330.0
	32	REAL, DIMENSION(NX),SAVE :: TBPVS,TBPVS0
	33	REAL, SAVE :: C1XPVS0,C2XPVS0,C1XPVS,C2XPVS
	34	!
	35	REAL, PRIVATE,PARAMETER :: &
	36	!--- Physical constants follow:
	37	& CP=1004.6, EPSQ=1.E-12, GRAV=9.806, RHOL=1000., RD=287.04 &
	38	& ,RV=461.5, T0C=273.15, XLS=2.834E6 &
	39	!--- Derived physical constants follow:
	40	& ,EPS=RD/RV, EPS1=RV/RD-1., EPSQ1=1.001*EPSQ &
	41	& ,RCP=1./CP, RCPRV=RCP/RV, RGRAV=1./GRAV, RRHOL=1./RHOL &
	42	& ,XLS1=XLSRCP, XLS2=XLSXLSRCPRV, XLS3=XLSXLS/RV &
	43	!--- Constants specific to the parameterization follow:
	44	!--- CLIMIT/CLIMIT1 are lower limits for treating accumulated precipitation
	45	& ,CLIMIT=10.*EPSQ, CLIMIT1=-CLIMIT &
	46	& ,C1=1./3. &
	47	& ,DMR1=.1E-3, DMR2=.2E-3, DMR3=.32E-3 &
	48	& ,XMR1=1.e6DMR1, XMR2=1.e6DMR2, XMR3=1.e6*DMR3
	49	INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3
	50	!
	51	! ======================================================================
	52	!--- Important tunable parameters that are exported to other modules
	53	! * RHgrd - threshold relative humidity for onset of condensation
	54	! * T_ICE - temperature (C) threshold at which all remaining liquid water
	55	! is glaciated to ice
	56	! * T_ICE_init - maximum temperature (C) at which ice nucleation occurs
	57	! * NLImax - maximum number concentrations (m**-3) of large ice (snow/graupel/sleet)
	58	! * NLImin - minimum number concentrations (m**-3) of large ice (snow/graupel/sleet)
	59	! * N0r0 - assumed intercept (m**-4) of rain drops if drop diameters are between 0.2 and 0.45 mm
	60	! * N0rmin - minimum intercept (m**-4) for rain drops
	61	! * NCW - number concentrations of cloud droplets (m**-3)
	62	! * FLARGE1, FLARGE2 - number fraction of large ice to total (large+snow) ice
	63	! at T>0C and in presence of sublimation (FLARGE1), otherwise in
	64	! presence of ice saturated/supersaturated conditions
	65	! ======================================================================
	66	REAL, PUBLIC,PARAMETER :: &
	67	& RHgrd=1. &
	68	& ,T_ICE=-40. &
	69	& ,T_ICEK=T0C+T_ICE &
	70	& ,T_ICE_init=-15. &
	71	& ,NLImax=5.E3 &
	72	& ,NLImin=1.E3 &
	73	& ,N0r0=8.E6 &
	74	& ,N0rmin=1.E4 &
	75	& ,NCW=100.E6 &
	76	& ,FLARGE1=1. &
	77	& ,FLARGE2=.2
	78	!--- Other public variables passed to other routines:
	79	REAL,PUBLIC,SAVE :: QAUT0
	80	REAL, PUBLIC,DIMENSION(MDImin:MDImax) :: MASSI
	81	!
	82	!
	83	CONTAINS
	84
	85	!-----------------------------------------------------------------------
	86	!-----------------------------------------------------------------------
	87	SUBROUTINE ETAMP_NEW (itimestep,DT,DX,DY, &
	88	& dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qt, &
	89	& LOWLYR,SR, &
	90	& F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
	91	& QC,QR,QS, &
	92	& mp_restart_state,tbpvs_state,tbpvs0_state, &
	93	& RAINNC,RAINNCV, &
	94	& ids,ide, jds,jde, kds,kde, &
	95	& ims,ime, jms,jme, kms,kme, &
	96	& its,ite, jts,jte, kts,kte )
	97	!-----------------------------------------------------------------------
	98	IMPLICIT NONE
	99	!-----------------------------------------------------------------------
	100	INTEGER, PARAMETER :: ITLO=-60, ITHI=40
	101
	102	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
	103	& ,IMS,IME,JMS,JME,KMS,KME &
	104	& ,ITS,ITE,JTS,JTE,KTS,KTE &
	105	& ,ITIMESTEP
	106
	107	REAL, INTENT(IN) :: DT,DX,DY
	108	REAL, INTENT(IN), DIMENSION(ims:ime, kms:kme, jms:jme):: &
	109	& dz8w,p_phy,pi_phy,rho_phy
	110	REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme):: &
	111	& th_phy,qv,qt
	112	REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
	113	& qc,qr,qs
	114	REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
	115	& F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
	116	REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: &
	117	& RAINNC,RAINNCV
	118	REAL, INTENT(OUT), DIMENSION(ims:ime,jms:jme):: SR
	119	!
	120	REAL,DIMENSION(*),INTENT(INOUT) :: MP_RESTART_STATE
	121	!
	122	REAL,DIMENSION(nx),INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
	123	!
	124	INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
	125
	126	!-----------------------------------------------------------------------
	127	! LOCAL VARS
	128	!-----------------------------------------------------------------------
	129
	130	! NSTATS,QMAX,QTOT are diagnostic vars
	131
	132	INTEGER,DIMENSION(ITLO:ITHI,4) :: NSTATS
	133	REAL, DIMENSION(ITLO:ITHI,5) :: QMAX
	134	REAL, DIMENSION(ITLO:ITHI,22):: QTOT
	135
	136	! SOME VARS WILL BE USED FOR DATA ASSIMILATION (DON'T NEED THEM NOW).
	137	! THEY ARE TREATED AS LOCAL VARS, BUT WILL BECOME STATE VARS IN THE
	138	! FUTURE. SO, WE DECLARED THEM AS MEMORY SIZES FOR THE FUTURE USE
	139
	140	! TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related
	141	! the microphysics scheme. Instead, they will be used by Eta precip
	142	! assimilation.
	143
	144	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
	145	& TLATGS_PHY,TRAIN_PHY
	146	REAL, DIMENSION(ims:ime,jms:jme):: APREC,PREC,ACPREC
	147	REAL, DIMENSION(its:ite, kts:kte, jts:jte):: t_phy
	148
	149	INTEGER :: I,J,K,KFLIP
	150	REAL :: WC
	151	!
	152	!-----------------------------------------------------------------------
	153	!**********************************************************************
	154	!-----------------------------------------------------------------------
	155	!
	156	MY_GROWTH(MY_T1:MY_T2)=MP_RESTART_STATE(MY_T1:MY_T2)
	157	!
	158	C1XPVS0=MP_RESTART_STATE(MY_T2+1)
	159	C2XPVS0=MP_RESTART_STATE(MY_T2+2)
	160	C1XPVS =MP_RESTART_STATE(MY_T2+3)
	161	C2XPVS =MP_RESTART_STATE(MY_T2+4)
	162	CIACW =MP_RESTART_STATE(MY_T2+5)
	163	CIACR =MP_RESTART_STATE(MY_T2+6)
	164	CRACW =MP_RESTART_STATE(MY_T2+7)
	165	CRAUT =MP_RESTART_STATE(MY_T2+8)
	166	!
	167	TBPVS(1:NX) =TBPVS_STATE(1:NX)
	168	TBPVS0(1:NX)=TBPVS0_STATE(1:NX)
	169	!
	170	DO j = jts,jte
	171	DO k = kts,kte
	172	DO i = its,ite
	173	t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
	174	qv(i,k,j)=qv(i,k,j)/(1.+qv(i,k,j)) !Convert to specific humidity
	175	ENDDO
	176	ENDDO
	177	ENDDO
	178
	179	! initial diagnostic variables and data assimilation vars
	180	! (will need to delete this part in the future)
	181
	182	DO k = 1,4
	183	DO i = ITLO,ITHI
	184	NSTATS(i,k)=0.
	185	ENDDO
	186	ENDDO
	187
	188	DO k = 1,5
	189	DO i = ITLO,ITHI
	190	QMAX(i,k)=0.
	191	ENDDO
	192	ENDDO
	193
	194	DO k = 1,22
	195	DO i = ITLO,ITHI
	196	QTOT(i,k)=0.
	197	ENDDO
	198	ENDDO
	199
	200	! initial data assimilation vars (will need to delete this part in the future)
	201
	202	DO j = jts,jte
	203	DO k = kts,kte
	204	DO i = its,ite
	205	TLATGS_PHY (i,k,j)=0.
	206	TRAIN_PHY (i,k,j)=0.
	207	ENDDO
	208	ENDDO
	209	ENDDO
	210
	211	DO j = jts,jte
	212	DO i = its,ite
	213	ACPREC(i,j)=0.
	214	APREC (i,j)=0.
	215	PREC (i,j)=0.
	216	SR (i,j)=0.
	217	ENDDO
	218	ENDDO
	219
	220	!-----------------------------------------------------------------------
	221
	222	CALL EGCP01DRV(DT,LOWLYR, &
	223	& APREC,PREC,ACPREC,SR,NSTATS,QMAX,QTOT, &
	224	& dz8w,rho_phy,qt,t_phy,qv,F_ICE_PHY,P_PHY, &
	225	& F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
	226	& ids,ide, jds,jde, kds,kde, &
	227	& ims,ime, jms,jme, kms,kme, &
	228	& its,ite, jts,jte, kts,kte )
	229	!-----------------------------------------------------------------------
	230
	231	DO j = jts,jte
	232	DO k = kts,kte
	233	DO i = its,ite
	234	th_phy(i,k,j) = t_phy(i,k,j)/pi_phy(i,k,j)
	235	qv(i,k,j)=qv(i,k,j)/(1.-qv(i,k,j)) !Convert to mixing ratio
	236	WC=qt(I,K,J)
	237	QS(I,K,J)=0.
	238	QR(I,K,J)=0.
	239	QC(I,K,J)=0.
	240	IF(F_ICE_PHY(I,K,J)>=1.)THEN
	241	QS(I,K,J)=WC
	242	ELSEIF(F_ICE_PHY(I,K,J)<=0.)THEN
	243	QC(I,K,J)=WC
	244	ELSE
	245	QS(I,K,J)=F_ICE_PHY(I,K,J)*WC
	246	QC(I,K,J)=WC-QS(I,K,J)
	247	ENDIF
	248	!
	249	IF(QC(I,K,J)>0..AND.F_RAIN_PHY(I,K,J)>0.)THEN
	250	IF(F_RAIN_PHY(I,K,J).GE.1.)THEN
	251	QR(I,K,J)=QC(I,K,J)
	252	QC(I,K,J)=0.
	253	ELSE
	254	QR(I,K,J)=F_RAIN_PHY(I,K,J)*QC(I,K,J)
	255	QC(I,K,J)=QC(I,K,J)-QR(I,K,J)
	256	ENDIF
	257	ENDIF
	258	ENDDO
	259	ENDDO
	260	ENDDO
	261	!
	262	! update rain (from m to mm)
	263
	264	DO j=jts,jte
	265	DO i=its,ite
	266	RAINNC(i,j)=APREC(i,j)*1000.+RAINNC(i,j)
	267	RAINNCV(i,j)=APREC(i,j)*1000.
	268	ENDDO
	269	ENDDO
	270	!
	271	MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
	272	MP_RESTART_STATE(MY_T2+1)=C1XPVS0
	273	MP_RESTART_STATE(MY_T2+2)=C2XPVS0
	274	MP_RESTART_STATE(MY_T2+3)=C1XPVS
	275	MP_RESTART_STATE(MY_T2+4)=C2XPVS
	276	MP_RESTART_STATE(MY_T2+5)=CIACW
	277	MP_RESTART_STATE(MY_T2+6)=CIACR
	278	MP_RESTART_STATE(MY_T2+7)=CRACW
	279	MP_RESTART_STATE(MY_T2+8)=CRAUT
	280	!
	281	TBPVS_STATE(1:NX) =TBPVS(1:NX)
	282	TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
	283
	284	!-----------------------------------------------------------------------
	285
	286	END SUBROUTINE ETAMP_NEW
	287
	288	!-----------------------------------------------------------------------
	289
	290	SUBROUTINE EGCP01DRV( &
	291	& DTPH,LOWLYR,APREC,PREC,ACPREC,SR, &
	292	& NSTATS,QMAX,QTOT, &
	293	& dz8w,RHO_PHY,CWM_PHY,T_PHY,Q_PHY,F_ICE_PHY,P_PHY, &
	294	& F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
	295	& ids,ide, jds,jde, kds,kde, &
	296	& ims,ime, jms,jme, kms,kme, &
	297	& its,ite, jts,jte, kts,kte)
	298	!-----------------------------------------------------------------------
	299	! DTPH Physics time step (s)
	300	! CWM_PHY (qt) Mixing ratio of the total condensate. kg/kg
	301	! Q_PHY Mixing ratio of water vapor. kg/kg
	302	! F_RAIN_PHY Fraction of rain.
	303	! F_ICE_PHY Fraction of ice.
	304	! F_RIMEF_PHY Mass ratio of rimed ice (rime factor).
	305	!
	306	!TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related the
	307	!micrphysics sechme. Instead, they will be used by Eta precip assimilation.
	308	!
	309	!NSTATS,QMAX,QTOT are used for diagnosis purposes.
	310	!
	311	!-----------------------------------------------------------------------
	312	!--- Variables APREC,PREC,ACPREC,SR are calculated for precip assimilation
	313	! and/or ZHAO's scheme in Eta and are not required by this microphysics
	314	! scheme itself.
	315	!--- NSTATS,QMAX,QTOT are used for diagnosis purposes only. They will be
	316	! printed out when PRINT_diag is true.
	317	!
	318	!-----------------------------------------------------------------------
	319	IMPLICIT NONE
	320	!-----------------------------------------------------------------------
	321	!
	322	INTEGER, PARAMETER :: ITLO=-60, ITHI=40
	323	LOGICAL, PARAMETER :: PRINT_diag=.FALSE.
	324	! VARIABLES PASSED IN/OUT
	325	INTEGER,INTENT(IN ) :: ids,ide, jds,jde, kds,kde &
	326	& ,ims,ime, jms,jme, kms,kme &
	327	& ,its,ite, jts,jte, kts,kte
	328	REAL,INTENT(IN) :: DTPH
	329	INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
	330	INTEGER,DIMENSION(ITLO:ITHI,4),INTENT(INOUT) :: NSTATS
	331	REAL,DIMENSION(ITLO:ITHI,5),INTENT(INOUT) :: QMAX
	332	REAL,DIMENSION(ITLO:ITHI,22),INTENT(INOUT) :: QTOT
	333	REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: &
	334	& APREC,PREC,ACPREC,SR
	335	REAL,DIMENSION( its:ite, kts:kte, jts:jte ),INTENT(INOUT) :: t_phy
	336	REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: &
	337	& dz8w,P_PHY,RHO_PHY
	338	REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(INOUT) :: &
	339	& CWM_PHY, F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY &
	340	& ,Q_PHY,TRAIN_PHY
	341	!
	342	!-----------------------------------------------------------------------
	343	!LOCAL VARIABLES
	344	!-----------------------------------------------------------------------
	345	!
	346	#define CACHE_FRIENDLY_MP_ETANEW
	347	#ifdef CACHE_FRIENDLY_MP_ETANEW
	348	# define TEMP_DIMS kts:kte,its:ite,jts:jte
	349	# define TEMP_DEX L,I,J
	350	#else
	351	# define TEMP_DIMS its:ite,jts:jte,kts:kte
	352	# define TEMP_DEX I,J,L
	353	#endif
	354	!
	355	INTEGER :: LSFC,I,J,I_index,J_index,L,K,KFLIP
	356	REAL,DIMENSION(TEMP_DIMS) :: CWM,T,Q,TRAIN,TLATGS,P
	357	REAL,DIMENSION(kts:kte,its:ite,jts:jte) :: F_ice,F_rain,F_RimeF
	358	INTEGER,DIMENSION(its:ite,jts:jte) :: LMH
	359	REAL :: TC,WC,QI,QR,QW,Fice,Frain,DUM,ASNOW,ARAIN
	360	REAL,DIMENSION(kts:kte) :: P_col,Q_col,T_col,QV_col,WC_col, &
	361	RimeF_col,QI_col,QR_col,QW_col, THICK_col,DPCOL
	362	REAL,DIMENSION(2) :: PRECtot,PRECmax
	363	!-----------------------------------------------------------------------
	364	!
	365	DO J=JTS,JTE
	366	DO I=ITS,ITE
	367	LMH(I,J) = KTE-LOWLYR(I,J)+1
	368	ENDDO
	369	ENDDO
	370
	371	DO 98 J=JTS,JTE
	372	DO 98 I=ITS,ITE
	373	DO L=KTS,KTE
	374	KFLIP=KTE+1-L
	375	CWM(TEMP_DEX)=CWM_PHY(I,KFLIP,J)
	376	T(TEMP_DEX)=T_PHY(I,KFLIP,J)
	377	Q(TEMP_DEX)=Q_PHY(I,KFLIP,J)
	378	P(TEMP_DEX)=P_PHY(I,KFLIP,J)
	379	TLATGS(TEMP_DEX)=TLATGS_PHY(I,KFLIP,J)
	380	TRAIN(TEMP_DEX)=TRAIN_PHY(I,KFLIP,J)
	381	F_ice(L,I,J)=F_ice_PHY(I,KFLIP,J)
	382	F_rain(L,I,J)=F_rain_PHY(I,KFLIP,J)
	383	F_RimeF(L,I,J)=F_RimeF_PHY(I,KFLIP,J)
	384	ENDDO
	385	98 CONTINUE
	386
	387	DO 100 J=JTS,JTE
	388	DO 100 I=ITS,ITE
	389	LSFC=LMH(I,J) ! "L" of surface
	390	!
	391	DO K=KTS,KTE
	392	KFLIP=KTE+1-K
	393	DPCOL(K)=RHO_PHY(I,KFLIP,J)GRAVdz8w(I,KFLIP,J)
	394	ENDDO
	395	!
	396	!
	397	!--- Initialize column data (1D arrays)
	398	!
	399	L=1
	400	IF (CWM(TEMP_DEX) .LE. EPSQ) CWM(TEMP_DEX)=EPSQ
	401	F_ice(1,I,J)=1.
	402	F_rain(1,I,J)=0.
	403	F_RimeF(1,I,J)=1.
	404	DO L=1,LSFC
	405	!
	406	!--- Pressure (Pa) = (Psfc-Ptop)*(ETA/ETA_sfc)+Ptop
	407	!
	408	P_col(L)=P(TEMP_DEX)
	409	!
	410	!--- Layer thickness = RHODZ = -DP/G = (Psfc-Ptop)D_ETA/(G*ETA_sfc)
	411	!
	412	THICK_col(L)=DPCOL(L)*RGRAV
	413	T_col(L)=T(TEMP_DEX)
	414	TC=T_col(L)-T0C
	415	QV_col(L)=max(EPSQ, Q(TEMP_DEX))
	416	IF (CWM(TEMP_DEX) .LE. EPSQ1) THEN
	417	WC_col(L)=0.
	418	IF (TC .LT. T_ICE) THEN
	419	F_ice(L,I,J)=1.
	420	ELSE
	421	F_ice(L,I,J)=0.
	422	ENDIF
	423	F_rain(L,I,J)=0.
	424	F_RimeF(L,I,J)=1.
	425	ELSE
	426	WC_col(L)=CWM(TEMP_DEX)
	427	ENDIF
	428	!
	429	!--- Determine composition of condensate in terms of
	430	! cloud water, ice, & rain
	431	!
	432	WC=WC_col(L)
	433	QI=0.
	434	QR=0.
	435	QW=0.
	436	Fice=F_ice(L,I,J)
	437	Frain=F_rain(L,I,J)
	438	IF (Fice .GE. 1.) THEN
	439	QI=WC
	440	ELSE IF (Fice .LE. 0.) THEN
	441	QW=WC
	442	ELSE
	443	QI=Fice*WC
	444	QW=WC-QI
	445	ENDIF
	446	IF (QW.GT.0. .AND. Frain.GT.0.) THEN
	447	IF (Frain .GE. 1.) THEN
	448	QR=QW
	449	QW=0.
	450	ELSE
	451	QR=Frain*QW
	452	QW=QW-QR
	453	ENDIF
	454	ENDIF
	455	RimeF_col(L)=F_RimeF(L,I,J) ! (real)
	456	QI_col(L)=QI
	457	QR_col(L)=QR
	458	QW_col(L)=QW
	459	ENDDO
	460	!
	461	!#######################################################################
	462	!
	463	!--- Perform the microphysical calculations in this column
	464	!
	465	I_index=I
	466	J_index=J
	467	CALL EGCP01COLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, LSFC, &
	468	& P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
	469	& THICK_col, WC_col,KTS,KTE,NSTATS,QMAX,QTOT )
	470
	471
	472	!
	473	!#######################################################################
	474	!
	475	!
	476	!--- Update storage arrays
	477	!
	478	DO L=1,LSFC
	479	TRAIN(TEMP_DEX)=(T_col(L)-T(TEMP_DEX))/DTPH
	480	TLATGS(TEMP_DEX)=T_col(L)-T(TEMP_DEX)
	481	T(TEMP_DEX)=T_col(L)
	482	Q(TEMP_DEX)=QV_col(L)
	483	CWM(TEMP_DEX)=WC_col(L)
	484	!
	485	!--- REAL*4 array storage
	486	!
	487	F_RimeF(L,I,J)=MAX(1., RimeF_col(L))
	488	IF (QI_col(L) .LE. EPSQ) THEN
	489	F_ice(L,I,J)=0.
	490	IF (T_col(L) .LT. T_ICEK) F_ice(L,I,J)=1.
	491	ELSE
	492	F_ice(L,I,J)=MAX( 0., MIN(1., QI_col(L)/WC_col(L)) )
	493	ENDIF
	494	IF (QR_col(L) .LE. EPSQ) THEN
	495	DUM=0
	496	ELSE
	497	DUM=QR_col(L)/(QR_col(L)+QW_col(L))
	498	ENDIF
	499	F_rain(L,I,J)=DUM
	500	!
	501	ENDDO
	502	!
	503	!--- Update accumulated precipitation statistics
	504	!
	505	!--- Surface precipitation statistics; SR is fraction of surface
	506	! precipitation (if >0) associated with snow
	507	!
	508	APREC(I,J)=(ARAIN+ASNOW)*RRHOL ! Accumulated surface precip (depth in m) !<--- Ying
	509	PREC(I,J)=PREC(I,J)+APREC(I,J)
	510	ACPREC(I,J)=ACPREC(I,J)+APREC(I,J)
	511	IF(APREC(I,J) .LT. 1.E-8) THEN
	512	SR(I,J)=0.
	513	ELSE
	514	SR(I,J)=RRHOL*ASNOW/APREC(I,J)
	515	ENDIF
	516	!
	517	!--- Debug statistics
	518	!
	519	IF (PRINT_diag) THEN
	520	PRECtot(1)=PRECtot(1)+ARAIN
	521	PRECtot(2)=PRECtot(2)+ASNOW
	522	PRECmax(1)=MAX(PRECmax(1), ARAIN)
	523	PRECmax(2)=MAX(PRECmax(2), ASNOW)
	524	ENDIF
	525
	526
	527	!#######################################################################
	528	!#######################################################################
	529	!
	530	100 CONTINUE ! End "I" & "J" loops
	531	DO 101 J=JTS,JTE
	532	DO 101 I=ITS,ITE
	533	DO L=KTS,KTE
	534	KFLIP=KTE+1-L
	535	CWM_PHY(I,KFLIP,J)=CWM(TEMP_DEX)
	536	T_PHY(I,KFLIP,J)=T(TEMP_DEX)
	537	Q_PHY(I,KFLIP,J)=Q(TEMP_DEX)
	538	TLATGS_PHY(I,KFLIP,J)=TLATGS(TEMP_DEX)
	539	TRAIN_PHY(I,KFLIP,J)=TRAIN(TEMP_DEX)
	540	F_ice_PHY(I,KFLIP,J)=F_ice(L,I,J)
	541	F_rain_PHY(I,KFLIP,J)=F_rain(L,I,J)
	542	F_RimeF_PHY(I,KFLIP,J)=F_RimeF(L,I,J)
	543	ENDDO
	544	101 CONTINUE
	545	END SUBROUTINE EGCP01DRV
	546	!
	547	!
	548	!###############################################################################
	549	! ***** VERSION OF MICROPHYSICS DESIGNED FOR HIGHER RESOLUTION MESO ETA MODEL
	550	! (1) Represents sedimentation by preserving a portion of the precipitation
	551	! through top-down integration from cloud-top. Modified procedure to
	552	! Zhao and Carr (1997).
	553	! (2) Microphysical equations are modified to be less sensitive to time
	554	! steps by use of Clausius-Clapeyron equation to account for changes in
	555	! saturation mixing ratios in response to latent heating/cooling.
	556	! (3) Prevent spurious temperature oscillations across 0C due to
	557	! microphysics.
	558	! (4) Uses lookup tables for: calculating two different ventilation
	559	! coefficients in condensation and deposition processes; accretion of
	560	! cloud water by precipitation; precipitation mass; precipitation rate
	561	! (and mass-weighted precipitation fall speeds).
	562	! (5) Assumes temperature-dependent variation in mean diameter of large ice
	563	! (Houze et al., 1979; Ryan et al., 1996).
	564	! -> 8/22/01: This relationship has been extended to colder temperatures
	565	! to parameterize smaller large-ice particles down to mean sizes of MDImin,
	566	! which is 50 microns reached at -55.9C.
	567	! (6) Attempts to differentiate growth of large and small ice, mainly for
	568	! improved transition from thin cirrus to thick, precipitating ice
	569	! anvils.
	570	! -> 8/22/01: This feature has been diminished by effectively adjusting to
	571	! ice saturation during depositional growth at temperatures colder than
	572	! -10C. Ice sublimation is calculated more explicitly. The logic is
	573	! that sources of are either poorly understood (e.g., nucleation for NWP)
	574	! or are not represented in the Eta model (e.g., detrainment of ice from
	575	! convection). Otherwise the model is too wet compared to the radiosonde
	576	! observations based on 1 Feb - 18 March 2001 retrospective runs.
	577	! (7) Top-down integration also attempts to treat mixed-phase processes,
	578	! allowing a mixture of ice and water. Based on numerous observational
	579	! studies, ice growth is based on nucleation at cloud top &
	580	! subsequent growth by vapor deposition and riming as the ice particles
	581	! fall through the cloud. Effective nucleation rates are a function
	582	! of ice supersaturation following Meyers et al. (JAM, 1992).
	583	! -> 8/22/01: The simulated relative humidities were far too moist compared
	584	! to the rawinsonde observations. This feature has been substantially
	585	! diminished, limited to a much narrower temperature range of 0 to -10C.
	586	! (8) Depositional growth of newly nucleated ice is calculated for large time
	587	! steps using Fig. 8 of Miller and Young (JAS, 1979), at 1 deg intervals
	588	! using their ice crystal masses calculated after 600 s of growth in water
	589	! saturated conditions. The growth rates are normalized by time step
	590	! assuming 3D growth with time**1.5 following eq. (6.3) in Young (1993).
	591	! -> 8/22/01: This feature has been effectively limited to 0 to -10C.
	592	! (9) Ice precipitation rates can increase due to increase in response to
	593	! cloud water riming due to (a) increased density & mass of the rimed
	594	! ice, and (b) increased fall speeds of rimed ice.
	595	! -> 8/22/01: This feature has been effectively limited to 0 to -10C.
	596	!###############################################################################
	597	!###############################################################################
	598	!
	599	SUBROUTINE EGCP01COLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, &
	600	& LSFC, P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
	601	& THICK_col, WC_col ,KTS,KTE,NSTATS,QMAX,QTOT)
	602	!
	603	!###############################################################################
	604	!###############################################################################
	605	!
	606	!-------------------------------------------------------------------------------
	607	!----- NOTE: Code is currently set up w/o threading!
	608	!-------------------------------------------------------------------------------
	609	!$$$ SUBPROGRAM DOCUMENTATION BLOCK
	610	! . . .
	611	! SUBPROGRAM: Grid-scale microphysical processes - condensation & precipitation
	612	! PRGRMMR: Ferrier ORG: W/NP22 DATE: 08-2001
	613	! PRGRMMR: Jin (Modification for WRF structure)
	614	!-------------------------------------------------------------------------------
	615	! ABSTRACT:
	616	! * Merges original GSCOND & PRECPD subroutines.
	617	! * Code has been substantially streamlined and restructured.
	618	! * Exchange between water vapor & small cloud condensate is calculated using
	619	! the original Asai (1965, J. Japan) algorithm. See also references to
	620	! Yau and Austin (1979, JAS), Rutledge and Hobbs (1983, JAS), and Tao et al.
	621	! (1989, MWR). This algorithm replaces the Sundqvist et al. (1989, MWR)
	622	! parameterization.
	623	!-------------------------------------------------------------------------------
	624	!
	625	! USAGE:
	626	! * CALL EGCP01COLUMN FROM SUBROUTINE EGCP01DRV
	627	!
	628	! INPUT ARGUMENT LIST:
	629	! DTPH - physics time step (s)
	630	! I_index - I index
	631	! J_index - J index
	632	! LSFC - Eta level of level above surface, ground
	633	! P_col - vertical column of model pressure (Pa)
	634	! QI_col - vertical column of model ice mixing ratio (kg/kg)
	635	! QR_col - vertical column of model rain ratio (kg/kg)
	636	! QV_col - vertical column of model water vapor specific humidity (kg/kg)
	637	! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
	638	! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
	639	! T_col - vertical column of model temperature (deg K)
	640	! THICK_col - vertical column of model mass thickness (density*height increment)
	641	! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
	642	!
	643	!
	644	! OUTPUT ARGUMENT LIST:
	645	! ARAIN - accumulated rainfall at the surface (kg)
	646	! ASNOW - accumulated snowfall at the surface (kg)
	647	! QV_col - vertical column of model water vapor specific humidity (kg/kg)
	648	! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
	649	! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
	650	! QI_col - vertical column of model ice mixing ratio (kg/kg)
	651	! QR_col - vertical column of model rain ratio (kg/kg)
	652	! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
	653	! T_col - vertical column of model temperature (deg K)
	654	!
	655	! OUTPUT FILES:
	656	! NONE
	657	!
	658	! Subprograms & Functions called:
	659	! * Real Function CONDENSE - cloud water condensation
	660	! * Real Function DEPOSIT - ice deposition (not sublimation)
	661	!
	662	! UNIQUE: NONE
	663	!
	664	! LIBRARY: NONE
	665	!
	666	! COMMON BLOCKS:
	667	! CMICRO_CONS - key constants initialized in GSMCONST
	668	! CMICRO_STATS - accumulated and maximum statistics
	669	! CMY_GROWTH - lookup table for growth of ice crystals in
	670	! water saturated conditions (Miller & Young, 1979)
	671	! IVENT_TABLES - lookup tables for ventilation effects of ice
	672	! IACCR_TABLES - lookup tables for accretion rates of ice
	673	! IMASS_TABLES - lookup tables for mass content of ice
	674	! IRATE_TABLES - lookup tables for precipitation rates of ice
	675	! IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
	676	! RVENT_TABLES - lookup tables for ventilation effects of rain
	677	! RACCR_TABLES - lookup tables for accretion rates of rain
	678	! RMASS_TABLES - lookup tables for mass content of rain
	679	! RVELR_TABLES - lookup tables for fall speeds of rain
	680	! RRATE_TABLES - lookup tables for precipitation rates of rain
	681	!
	682	! ATTRIBUTES:
	683	! LANGUAGE: FORTRAN 90
	684	! MACHINE : IBM SP
	685	!
	686	!
	687	!-------------------------------------------------------------------------
	688	!--------------- Arrays & constants in argument list ---------------------
	689	!-------------------------------------------------------------------------
	690	!
	691	IMPLICIT NONE
	692	!
	693	INTEGER,INTENT(IN) :: KTS,KTE,I_index, J_index, LSFC
	694	REAL,INTENT(INOUT) :: ARAIN, ASNOW
	695	REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: P_col, QI_col,QR_col &
	696	& ,QV_col ,QW_col, RimeF_col, T_col, THICK_col,WC_col
	697	!
	698	!-------------------------------------------------------------------------
	699	!-------------- Common blocks for microphysical statistics ---------------
	700	!-------------------------------------------------------------------------
	701	!
	702	!-------------------------------------------------------------------------
	703	!--------- Common blocks for constants initialized in GSMCONST ----------
	704	!
	705	INTEGER, PARAMETER :: ITLO=-60, ITHI=40
	706	INTEGER,INTENT(INOUT) :: NSTATS(ITLO:ITHI,4)
	707	REAL,INTENT(INOUT) :: QMAX(ITLO:ITHI,5),QTOT(ITLO:ITHI,22)
	708	!
	709	!-------------------------------------------------------------------------
	710	!--------------- Common blocks for various lookup tables -----------------
	711	!
	712	!--- Discretized growth rates of small ice crystals after their nucleation
	713	! at 1 C intervals from -1 C to -35 C, based on calculations by Miller
	714	! and Young (1979, JAS) after 600 s of growth. Resultant growth rates
	715	! are multiplied by physics time step in GSMCONST.
	716	!
	717	!-------------------------------------------------------------------------
	718	!
	719	!--- Mean ice-particle diameters varying from 50 microns to 1000 microns
	720	! (1 mm), assuming an exponential size distribution.
	721	!
	722	!---- Meaning of the following arrays:
	723	! - mdiam - mean diameter (m)
	724	! - VENTI1 - integrated quantity associated w/ ventilation effects
	725	! (capacitance only) for calculating vapor deposition onto ice
	726	! - VENTI2 - integrated quantity associated w/ ventilation effects
	727	! (with fall speed) for calculating vapor deposition onto ice
	728	! - ACCRI - integrated quantity associated w/ cloud water collection by ice
	729	! - MASSI - integrated quantity associated w/ ice mass
	730	! - VSNOWI - mass-weighted fall speed of snow (large ice), used to calculate
	731	! precipitation rates
	732	!
	733	!
	734	!-------------------------------------------------------------------------
	735	!
	736	!--- VEL_RF - velocity increase of rimed particles as functions of crude
	737	! particle size categories (at 0.1 mm intervals of mean ice particle
	738	! sizes) and rime factor (different values of Rime Factor of 1.1**N,
	739	! where N=0 to Nrime).
	740	!
	741	!-------------------------------------------------------------------------
	742	!
	743	!--- Mean rain drop diameters varying from 50 microns (0.05 mm) to 450 microns
	744	! (0.45 mm), assuming an exponential size distribution.
	745	!
	746	!-------------------------------------------------------------------------
	747	!------- Key parameters, local variables, & important comments ---------
	748	!-----------------------------------------------------------------------
	749	!
	750	!--- TOLER => Tolerance or precision for accumulated precipitation
	751	!
	752	REAL, PARAMETER :: TOLER=5.E-7, C2=1./6., RHO0=1.194, Xratio=.025
	753	!
	754	!--- If BLEND=1:
	755	! precipitation (large) ice amounts are estimated at each level as a
	756	! blend of ice falling from the grid point above and the precip ice
	757	! present at the start of the time step (see TOT_ICE below).
	758	!--- If BLEND=0:
	759	! precipitation (large) ice amounts are estimated to be the precip
	760	! ice present at the start of the time step.
	761	!
	762	!--- Extended to include sedimentation of rain on 2/5/01
	763	!
	764	REAL, PARAMETER :: BLEND=1.
	765	!
	766	!--- This variable is for debugging purposes (if .true.)
	767	!
	768	LOGICAL, PARAMETER :: PRINT_diag=.FALSE.
	769	!
	770	!-----------------------------------------------------------------------
	771	!--- Local variables
	772	!-----------------------------------------------------------------------
	773	!
	774	REAL EMAIRI, N0r, NLICE, NSmICE
	775	LOGICAL CLEAR, ICE_logical, DBG_logical, RAIN_logical
	776	INTEGER :: IDR,INDEX_MY,INDEXR,INDEXR1,INDEXS,IPASS,ITDX,IXRF, &
	777	& IXS,LBEF,L
	778	!
	779	REAL :: ABI,ABW,AIEVP,ARAINnew,ASNOWnew,BLDTRH,BUDGET, &
	780	& CREVP,DELI,DELR,DELT,DELV,DELW,DENOMF, &
	781	& DENOMI,DENOMW,DENOMWI,DIDEP, &
	782	& DIEVP,DIFFUS,DLI,DTPH,DTRHO,DUM,DUM1, &
	783	& DUM2,DWV0,DWVI,DWVR,DYNVIS,ESI,ESW,FIR,FLARGE,FLIMASS, &
	784	& FSMALL,FWR,FWS,GAMMAR,GAMMAS, &
	785	& PCOND,PIACR,PIACW,PIACWI,PIACWR,PICND,PIDEP,PIDEP_max, &
	786	& PIEVP,PILOSS,PIMLT,PP,PRACW,PRAUT,PREVP,PRLOSS, &
	787	& QI,QInew,QLICE,QR,QRnew,QSI,QSIgrd,QSInew,QSW,QSW0, &
	788	& QSWgrd,QSWnew,QT,QTICE,QTnew,QTRAIN,QV,QW,QW0,QWnew, &
	789	& RFACTOR,RHO,RIMEF,RIMEF1,RQR,RR,RRHO,SFACTOR, &
	790	& TC,TCC,TFACTOR,THERM_COND,THICK,TK,TK2,TNEW, &
	791	& TOT_ICE,TOT_ICEnew,TOT_RAIN,TOT_RAINnew, &
	792	& VEL_INC,VENTR,VENTIL,VENTIS,VRAIN1,VRAIN2,VRIMEF,VSNOW, &
	793	& WC,WCnew,WSgrd,WS,WSnew,WV,WVnew,WVQW, &
	794	& XLF,XLF1,XLI,XLV,XLV1,XLV2,XLIMASS,XRF,XSIMASS
	795	!
	796	!#######################################################################
	797	!########################## Begin Execution ############################
	798	!#######################################################################
	799	!
	800	!
	801	ARAIN=0. ! Accumulated rainfall into grid box from above (kg/m**2)
	802	ASNOW=0. ! Accumulated snowfall into grid box from above (kg/m**2)
	803	!
	804	!-----------------------------------------------------------------------
	805	!------------ Loop from top (L=1) to surface (L=LSFC) ------------------
	806	!-----------------------------------------------------------------------
	807	!
	808
	809	DO 10 L=1,LSFC
	810
	811	!--- Skip this level and go to the next lower level if no condensate
	812	! and very low specific humidities
	813	!
	814	IF (QV_col(L).LE.EPSQ .AND. WC_col(L).LE.EPSQ) GO TO 10
	815	!
	816	!-----------------------------------------------------------------------
	817	!------------ Proceed with cloud microphysics calculations -------------
	818	!-----------------------------------------------------------------------
	819	!
	820	TK=T_col(L) ! Temperature (deg K)
	821	TC=TK-T0C ! Temperature (deg C)
	822	PP=P_col(L) ! Pressure (Pa)
	823	QV=QV_col(L) ! Specific humidity of water vapor (kg/kg)
	824	WV=QV/(1.-QV) ! Water vapor mixing ratio (kg/kg)
	825	WC=WC_col(L) ! Grid-scale mixing ratio of total condensate (water or ice; kg/kg)
	826	!
	827	!-----------------------------------------------------------------------
	828	!--- Moisture variables below are mixing ratios & not specifc humidities
	829	!-----------------------------------------------------------------------
	830	!
	831	CLEAR=.TRUE.
	832	!
	833	!--- This check is to determine grid-scale saturation when no condensate is present
	834	!
	835	ESW=1000.*FPVS0(TK) ! Saturation vapor pressure w/r/t water
	836	QSW=EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
	837	WS=QSW ! General saturation mixing ratio (water/ice)
	838	IF (TC .LT. 0.) THEN
	839	ESI=1000.*FPVS(TK) ! Saturation vapor pressure w/r/t ice
	840	QSI=EPS*ESI/(PP-ESI) ! Saturation mixing ratio w/r/t water
	841	WS=QSI ! General saturation mixing ratio (water/ice)
	842	ENDIF
	843	!
	844	!--- Effective grid-scale Saturation mixing ratios
	845	!
	846	QSWgrd=RHgrd*QSW
	847	QSIgrd=RHgrd*QSI
	848	WSgrd=RHgrd*WS
	849	!
	850	!--- Check if air is subsaturated and w/o condensate
	851	!
	852	IF (WV.GT.WSgrd .OR. WC.GT.EPSQ) CLEAR=.FALSE.
	853	!
	854	!--- Check if any rain is falling into layer from above
	855	!
	856	IF (ARAIN .GT. CLIMIT) THEN
	857	CLEAR=.FALSE.
	858	ELSE
	859	ARAIN=0.
	860	ENDIF
	861	!
	862	!--- Check if any ice is falling into layer from above
	863	!
	864	!--- NOTE that "SNOW" in variable names is synonomous with
	865	! large, precipitation ice particles
	866	!
	867	IF (ASNOW .GT. CLIMIT) THEN
	868	CLEAR=.FALSE.
	869	ELSE
	870	ASNOW=0.
	871	ENDIF
	872	!
	873	!-----------------------------------------------------------------------
	874	!-- Loop to the end if in clear, subsaturated air free of condensate ---
	875	!-----------------------------------------------------------------------
	876	!
	877	IF (CLEAR) GO TO 10
	878	!
	879	!-----------------------------------------------------------------------
	880	!--------- Initialize RHO, THICK & microphysical processes -------------
	881	!-----------------------------------------------------------------------
	882	!
	883	!
	884	!--- Virtual temperature, TV=T(1./EPS-1)Q, Q is specific humidity;
	885	! (see pp. 63-65 in Fleagle & Businger, 1963)
	886	!
	887	RHO=PP/(RDTK(1.+EPS1QV)) ! Air density (kg/m*3)
	888	RRHO=1./RHO ! Reciprocal of air density
	889	DTRHO=DTPHRHO ! Time step air density
	890	BLDTRH=BLENDDTRHO ! Blend parameter time step * air density
	891	THICK=THICK_col(L) ! Layer thickness = RHODZ = -DP/G = (Psfc-Ptop)D_ETA/(G*ETA_sfc)
	892	!
	893	ARAINnew=0. ! Updated accumulated rainfall
	894	ASNOWnew=0. ! Updated accumulated snowfall
	895	QI=QI_col(L) ! Ice mixing ratio
	896	QInew=0. ! Updated ice mixing ratio
	897	QR=QR_col(L) ! Rain mixing ratio
	898	QRnew=0. ! Updated rain ratio
	899	QW=QW_col(L) ! Cloud water mixing ratio
	900	QWnew=0. ! Updated cloud water ratio
	901	!
	902	PCOND=0. ! Condensation (>0) or evaporation (<0) of cloud water (kg/kg)
	903	PIDEP=0. ! Deposition (>0) or sublimation (<0) of ice crystals (kg/kg)
	904	PIACW=0. ! Cloud water collection (riming) by precipitation ice (kg/kg; >0)
	905	PIACWI=0. ! Growth of precip ice by riming (kg/kg; >0)
	906	PIACWR=0. ! Shedding of accreted cloud water to form rain (kg/kg; >0)
	907	PIACR=0. ! Freezing of rain onto large ice at supercooled temps (kg/kg; >0)
	908	PICND=0. ! Condensation (>0) onto wet, melting ice (kg/kg)
	909	PIEVP=0. ! Evaporation (<0) from wet, melting ice (kg/kg)
	910	PIMLT=0. ! Melting ice (kg/kg; >0)
	911	PRAUT=0. ! Cloud water autoconversion to rain (kg/kg; >0)
	912	PRACW=0. ! Cloud water collection (accretion) by rain (kg/kg; >0)
	913	PREVP=0. ! Rain evaporation (kg/kg; <0)
	914	!
	915	!--- Double check input hydrometeor mixing ratios
	916	!
	917	! DUM=WC-(QI+QW+QR)
	918	! DUM1=ABS(DUM)
	919	! DUM2=TOLER*MIN(WC, QI+QW+QR)
	920	! IF (DUM1 .GT. DUM2) THEN
	921	! WRITE(6,"(/2(a,i4),a,i2)") '{@ i=',I_index,' j=',J_index,
	922	! & ' L=',L
	923	! WRITE(6,"(4(a12,g11.4,1x))")
	924	! & '{@ TCold=',TC,'P=',.01*PP,'DIFF=',DUM,'WCold=',WC,
	925	! & '{@ QIold=',QI,'QWold=',QW,'QRold=',QR
	926	! ENDIF
	927	!
	928	!***********************************************************************
	929	!********* MAIN MICROPHYSICS CALCULATIONS NOW FOLLOW! **************
	930	!***********************************************************************
	931	!
	932	!--- Calculate a few variables, which are used more than once below
	933	!
	934	!--- Latent heat of vaporization as a function of temperature from
	935	! Bolton (1980, JAS)
	936	!
	937	XLV=3.148E6-2370*TK ! Latent heat of vaporization (Lv)
	938	XLF=XLS-XLV ! Latent heat of fusion (Lf)
	939	XLV1=XLV*RCP ! Lv/Cp
	940	XLF1=XLF*RCP ! Lf/Cp
	941	TK2=1./(TKTK) ! 1./TK*2
	942	XLV2=XLVXLVQSWTK2/RV ! Lv2Qsw/(RvTK*2)
	943	DENOMW=1.+XLV2*RCP ! Denominator term, Clausius-Clapeyron correction
	944	!
	945	!--- Basic thermodynamic quantities
	946	! * DYNVIS - dynamic viscosity [ kg/(m*s) ]
	947	! * THERM_COND - thermal conductivity [ J/(msK) ]
	948	! * DIFFUS - diffusivity of water vapor [ m**2/s ]
	949	!
	950	TFACTOR=TK**1.5/(TK+120.)
	951	DYNVIS=1.496E-6*TFACTOR
	952	THERM_COND=2.116E-3*TFACTOR
	953	DIFFUS=8.794E-5TK*1.81/PP
	954	!
	955	!--- Air resistance term for the fall speed of ice following the
	956	! basic research by Heymsfield, Kajikawa, others
	957	!
	958	GAMMAS=(1.E5/PP)**C1
	959	!
	960	!--- Air resistance for rain fall speed (Beard, 1985, JAS, p.470)
	961	!
	962	GAMMAR=(RHO0/RHO)**.4
	963	!
	964	!----------------------------------------------------------------------
	965	!------------- IMPORTANT MICROPHYSICS DECISION TREE -----------------
	966	!----------------------------------------------------------------------
	967	!
	968	!--- Determine if conditions supporting ice are present
	969	!
	970	IF (TC.LT.0. .OR. QI.GT. EPSQ .OR. ASNOW.GT.CLIMIT) THEN
	971	ICE_logical=.TRUE.
	972	ELSE
	973	ICE_logical=.FALSE.
	974	QLICE=0.
	975	QTICE=0.
	976	ENDIF
	977	!
	978	!--- Determine if rain is present
	979	!
	980	RAIN_logical=.FALSE.
	981	IF (ARAIN.GT.CLIMIT .OR. QR.GT.EPSQ) RAIN_logical=.TRUE.
	982	!
	983	IF (ICE_logical) THEN
	984	!
	985	!--- IMPORTANT: Estimate time-averaged properties.
	986	!
	987	!---
	988	! * FLARGE - ratio of number of large ice to total (large & small) ice
	989	! * FSMALL - ratio of number of small ice crystals to large ice particles
	990	! -> Small ice particles are assumed to have a mean diameter of 50 microns.
	991	! * XSIMASS - used for calculating small ice mixing ratio
	992	!---
	993	! * TOT_ICE - total mass (small & large) ice before microphysics,
	994	! which is the sum of the total mass of large ice in the
	995	! current layer and the input flux of ice from above
	996	! * PILOSS - greatest loss (<0) of total (small & large) ice by
	997	! sublimation, removing all of the ice falling from above
	998	! and the ice within the layer
	999	! * RimeF1 - Rime Factor, which is the mass ratio of total (unrimed & rimed)
	1000	! ice mass to the unrimed ice mass (>=1)
	1001	! * VrimeF - the velocity increase due to rime factor or melting (ratio, >=1)
	1002	! * VSNOW - Fall speed of rimed snow w/ air resistance correction
	1003	! * EMAIRI - equivalent mass of air associated layer and with fall of snow into layer
	1004	! * XLIMASS - used for calculating large ice mixing ratio
	1005	! * FLIMASS - mass fraction of large ice
	1006	! * QTICE - time-averaged mixing ratio of total ice
	1007	! * QLICE - time-averaged mixing ratio of large ice
	1008	! * NLICE - time-averaged number concentration of large ice
	1009	! * NSmICE - number concentration of small ice crystals at current level
	1010	!---
	1011	!--- Assumed number fraction of large ice particles to total (large & small)
	1012	! ice particles, which is based on a general impression of the literature.
	1013	!
	1014	WVQW=WV+QW ! Water vapor & cloud water
	1015	!
	1016
	1017
	1018	IF (TC.GE.0. .OR. WVQW.LT.QSIgrd) THEN
	1019	!
	1020	!--- Eliminate small ice particle contributions for melting & sublimation
	1021	!
	1022	FLARGE=FLARGE1
	1023	ELSE
	1024	!
	1025	!--- Enhanced number of small ice particles during depositional growth
	1026	! (effective only when 0C > T >= T_ice [-10C] )
	1027	!
	1028	FLARGE=FLARGE2
	1029	!
	1030	!--- Larger number of small ice particles due to rime splintering
	1031	!
	1032	IF (TC.GE.-8. .AND. TC.LE.-3.) FLARGE=.5*FLARGE
	1033	!
	1034	ENDIF ! End IF (TC.GE.0. .OR. WVQW.LT.QSIgrd)
	1035	FSMALL=(1.-FLARGE)/FLARGE
	1036	XSIMASS=RRHOMASSI(MDImin)FSMALL
	1037	IF (QI.LE.EPSQ .AND. ASNOW.LE.CLIMIT) THEN
	1038	INDEXS=MDImin
	1039	TOT_ICE=0.
	1040	PILOSS=0.
	1041	RimeF1=1.
	1042	VrimeF=1.
	1043	VEL_INC=GAMMAS
	1044	VSNOW=0.
	1045	EMAIRI=THICK
	1046	XLIMASS=RRHORimeF1MASSI(INDEXS)
	1047	FLIMASS=XLIMASS/(XLIMASS+XSIMASS)
	1048	QLICE=0.
	1049	QTICE=0.
	1050	NLICE=0.
	1051	NSmICE=0.
	1052	ELSE
	1053	!
	1054	!--- For T<0C mean particle size follows Houze et al. (JAS, 1979, p. 160),
	1055	! converted from Fig. 5 plot of LAMDAs. Similar set of relationships
	1056	! also shown in Fig. 8 of Ryan (BAMS, 1996, p. 66).
	1057	!
	1058	DUM=XMImaxEXP(.0536TC)
	1059	INDEXS=MIN(MDImax, MAX(MDImin, INT(DUM) ) )
	1060	TOT_ICE=THICKQI+BLENDASNOW
	1061	PILOSS=-TOT_ICE/THICK
	1062	LBEF=MAX(1,L-1)
	1063	DUM1=RimeF_col(LBEF)
	1064	DUM2=RimeF_col(L)
	1065	RimeF1=(DUM2THICKQI+DUM1BLENDASNOW)/TOT_ICE
	1066	RimeF1=MIN(RimeF1, RFmax)
	1067	DO IPASS=0,1
	1068	IF (RimeF1 .LE. 1.) THEN
	1069	RimeF1=1.
	1070	VrimeF=1.
	1071	ELSE
	1072	IXS=MAX(2, MIN(INDEXS/100, 9))
	1073	XRF=10.492*ALOG(RimeF1)
	1074	IXRF=MAX(0, MIN(INT(XRF), Nrime))
	1075	IF (IXRF .GE. Nrime) THEN
	1076	VrimeF=VEL_RF(IXS,Nrime)
	1077	ELSE
	1078	VrimeF=VEL_RF(IXS,IXRF)+(XRF-FLOAT(IXRF))* &
	1079	& (VEL_RF(IXS,IXRF+1)-VEL_RF(IXS,IXRF))
	1080	ENDIF
	1081	ENDIF ! End IF (RimeF1 .LE. 1.)
	1082	VEL_INC=GAMMAS*VrimeF
	1083	VSNOW=VEL_INC*VSNOWI(INDEXS)
	1084	EMAIRI=THICK+BLDTRH*VSNOW
	1085	XLIMASS=RRHORimeF1MASSI(INDEXS)
	1086	FLIMASS=XLIMASS/(XLIMASS+XSIMASS)
	1087	QTICE=TOT_ICE/EMAIRI
	1088	QLICE=FLIMASS*QTICE
	1089	NLICE=QLICE/XLIMASS
	1090	NSmICE=Fsmall*NLICE
	1091	!
	1092	IF ( (NLICE.GE.NLImin .AND. NLICE.LE.NLImax) &
	1093	& .OR. IPASS.EQ.1) THEN
	1094	EXIT
	1095	ELSE
	1096	IF (TC < 0) THEN
	1097	XLI=RHO*(QTICE/DUM-XSIMASS)/RimeF1
	1098	IF (XLI .LE. MASSI(MDImin) ) THEN
	1099	INDEXS=MDImin
	1100	ELSE IF (XLI .LE. MASSI(450) ) THEN
	1101	DLI=9.5885E5XLI*.42066 ! DLI in microns
	1102	INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
	1103	ELSE IF (XLI .LE. MASSI(MDImax) ) THEN
	1104	DLI=3.9751E6XLI*.49870 ! DLI in microns
	1105	INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
	1106	ELSE
	1107	INDEXS=MDImax
	1108	ENDIF ! End IF (XLI .LE. MASSI(MDImin) )
	1109	ENDIF ! End IF (TC < 0)
	1110	!
	1111	!--- Reduce excessive accumulation of ice at upper levels
	1112	! associated with strong grid-resolved ascent
	1113	!
	1114	!--- Force NLICE to be between NLImin and NLImax
	1115	!
	1116	!
	1117	!--- 8/22/01: Increase density of large ice if maximum limits
	1118	! are reached for number concentration (NLImax) and mean size
	1119	! (MDImax). Done to increase fall out of ice.
	1120	!
	1121	DUM=MAX(NLImin, MIN(NLImax, NLICE) )
	1122	IF (DUM.GE.NLImax .AND. INDEXS.GE.MDImax) &
	1123	& RimeF1=RHO*(QTICE/NLImax-XSIMASS)/MASSI(INDEXS)
	1124	! WRITE(6,"(4(a12,g11.4,1x))")
	1125	! & '{$ TC=',TC,'P=',.01*PP,'NLICE=',NLICE,'DUM=',DUM,
	1126	! & '{$ XLI=',XLI,'INDEXS=',FLOAT(INDEXS),'RHO=',RHO,'QTICE=',QTICE,
	1127	! & '{$ XSIMASS=',XSIMASS,'RimeF1=',RimeF1
	1128	ENDIF ! End IF ( (NLICE.GE.NLImin .AND. NLICE.LE.NLImax) ...
	1129	ENDDO ! End DO IPASS=0,1
	1130	ENDIF ! End IF (QI.LE.EPSQ .AND. ASNOW.LE.CLIMIT)
	1131	ENDIF ! End IF (ICE_logical)
	1132	!
	1133	!----------------------------------------------------------------------
	1134	!--------------- Calculate individual processes -----------------------
	1135	!----------------------------------------------------------------------
	1136	!
	1137	!--- Cloud water autoconversion to rain and collection by rain
	1138	!
	1139	IF (QW.GT.EPSQ .AND. TC.GE.T_ICE) THEN
	1140	!
	1141	!--- QW0 could be modified based on land/sea properties,
	1142	! presence of convection, etc. This is why QAUT0 and CRAUT
	1143	! are passed into the subroutine as externally determined
	1144	! parameters. Can be changed in the future if desired.
	1145	!
	1146	QW0=QAUT0*RRHO
	1147	PRAUT=MAX(0., QW-QW0)*CRAUT
	1148	IF (QLICE .GT. EPSQ) THEN
	1149	!
	1150	!--- Collection of cloud water by large ice particles ("snow")
	1151	! PIACWI=PIACW for riming, PIACWI=0 for shedding
	1152	!
	1153	FWS=MIN(1., CIACWVEL_INCNLICEACCRI(INDEXS)/PP*C1)
	1154	PIACW=FWS*QW
	1155	IF (TC .LT. 0.) PIACWI=PIACW ! Large ice riming
	1156	ENDIF ! End IF (QLICE .GT. EPSQ)
	1157	ENDIF ! End IF (QW.GT.EPSQ .AND. TC.GE.T_ICE)
	1158	!
	1159	!----------------------------------------------------------------------
	1160	!--- Loop around some of the ice-phase processes if no ice should be present
	1161	!----------------------------------------------------------------------
	1162	!
	1163	IF (ICE_logical .EQV. .FALSE.) GO TO 20
	1164	!
	1165	!--- Now the pretzel logic of calculating ice deposition
	1166	!
	1167	IF (TC.LT.T_ICE .AND. (WV.GT.QSIgrd .OR. QW.GT.EPSQ)) THEN
	1168	!
	1169	!--- Adjust to ice saturation at T<T_ICE (-10C) if supersaturated.
	1170	! Sources of ice due to nucleation and convective detrainment are
	1171	! either poorly understood, poorly resolved at typical NWP
	1172	! resolutions, or are not represented (e.g., no detrained
	1173	! condensate in BMJ Cu scheme).
	1174	!
	1175	PCOND=-QW
	1176	DUM1=TK+XLV1*PCOND ! Updated (dummy) temperature (deg K)
	1177	DUM2=WV+QW ! Updated (dummy) water vapor mixing ratio
	1178	DUM=1000.*FPVS(DUM1) ! Updated (dummy) saturation vapor pressure w/r/t ice
	1179	DUM=RHgrdEPSDUM/(PP-DUM) ! Updated (dummy) saturation mixing ratio w/r/t ice
	1180	IF (DUM2 .GT. DUM) PIDEP=DEPOSIT (PP, DUM1, DUM2)
	1181	DWVi=0. ! Used only for debugging
	1182	!
	1183	ELSE IF (TC .LT. 0.) THEN
	1184	!
	1185	!--- These quantities are handy for ice deposition/sublimation
	1186	! PIDEP_max - max deposition or minimum sublimation to ice saturation
	1187	!
	1188	DENOMI=1.+XLS2QSITK2
	1189	DWVi=MIN(WVQW,QSW)-QSI
	1190	PIDEP_max=MAX(PILOSS, DWVi/DENOMI)
	1191	IF (QTICE .GT. 0.) THEN
	1192	!
	1193	!--- Calculate ice deposition/sublimation
	1194	! * SFACTOR - [VEL_INC*.5][Schmidt*(1./3.)][(RHO/DYNVIS)**.5],
	1195	! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
	1196	! * Units: SFACTOR - s.5/m ; ABI - m2/s ; NLICE - m**-3 ;
	1197	! VENTIL, VENTIS - m**-2 ; VENTI1 - m ;
	1198	! VENTI2 - m2/s.5 ; DIDEP - unitless
	1199	!
	1200	SFACTOR=VEL_INC*.5(RHO/(DIFFUSDIFFUSDYNVIS))**C2
	1201	ABI=1./(RHOXLS3QSI*TK2/THERM_COND+1./DIFFUS)
	1202	!
	1203	!--- VENTIL - Number concentration * ventilation factors for large ice
	1204	!--- VENTIS - Number concentration * ventilation factors for small ice
	1205	!
	1206	!--- Variation in the number concentration of ice with time is not
	1207	! accounted for in these calculations (could be in the future).
	1208	!
	1209	VENTIL=(VENTI1(INDEXS)+SFACTORVENTI2(INDEXS))NLICE
	1210	VENTIS=(VENTI1(MDImin)+SFACTORVENTI2(MDImin))NSmICE
	1211	DIDEP=ABI(VENTIL+VENTIS)DTPH
	1212	!
	1213	!--- Account for change in water vapor supply w/ time
	1214	!
	1215	IF (DIDEP .GE. Xratio)then
	1216	DIDEP=(1.-EXP(-DIDEP*DENOMI))/DENOMI
	1217	endif
	1218	IF (DWVi .GT. 0.) THEN
	1219	PIDEP=MIN(DWVi*DIDEP, PIDEP_max)
	1220	ELSE IF (DWVi .LT. 0.) THEN
	1221	PIDEP=MAX(DWVi*DIDEP, PIDEP_max)
	1222	ENDIF
	1223	!
	1224	ELSE IF (WVQW.GT.QSI .AND. TC.LE.T_ICE_init) THEN
	1225	!
	1226	!--- Ice nucleation in near water-saturated conditions. Ice crystal
	1227	! growth during time step calculated using Miller & Young (1979, JAS).
	1228	!--- These deposition rates could drive conditions below water saturation,
	1229	! which is the basis of these calculations. Intended to approximate
	1230	! more complex & computationally intensive calculations.
	1231	!
	1232	INDEX_MY=MAX(MY_T1, MIN( INT(.5-TC), MY_T2 ) )
	1233	!
	1234	!--- DUM1 is the supersaturation w/r/t ice at water-saturated conditions
	1235	!
	1236	!--- DUM2 is the number of ice crystals nucleated at water-saturated
	1237	! conditions based on Meyers et al. (JAM, 1992).
	1238	!
	1239	!--- Prevent unrealistically large ice initiation (limited by PIDEP_max)
	1240	! if DUM2 values are increased in future experiments
	1241	!
	1242	DUM1=QSW/QSI-1.
	1243	DUM2=1.E3EXP(12.96DUM1-.639)
	1244	PIDEP=MIN(PIDEP_max, DUM2MY_GROWTH(INDEX_MY)RRHO)
	1245	!
	1246	ENDIF ! End IF (QTICE .GT. 0.)
	1247	!
	1248	ENDIF ! End IF (TC.LT.T_ICE .AND. (WV.GT.QSIgrd .OR. QW.GT.EPSQ))
	1249	!
	1250	!------------------------------------------------------------------------
	1251	!
	1252	20 CONTINUE ! Jump here if conditions for ice are not present
	1253
	1254
	1255	!
	1256	!------------------------------------------------------------------------
	1257	!
	1258	!--- Cloud water condensation
	1259	!
	1260	IF (TC.GE.T_ICE .AND. (QW.GT.EPSQ .OR. WV.GT.QSWgrd)) THEN
	1261	IF (PIACWI.EQ.0. .AND. PIDEP.EQ.0.) THEN
	1262	PCOND=CONDENSE (PP, QW, TK, WV)
	1263	ELSE
	1264	!
	1265	!--- Modify cloud condensation in response to ice processes
	1266	!
	1267	DUM=XLVQSWgrdRCPRV*TK2
	1268	DENOMWI=1.+XLS*DUM
	1269	DENOMF=XLF*DUM
	1270	DUM=MAX(0., PIDEP)
	1271	PCOND=(WV-QSWgrd-DENOMWIDUM-DENOMFPIACWI)/DENOMW
	1272	DUM1=-QW
	1273	DUM2=PCOND-PIACW
	1274	IF (DUM2 .LT. DUM1) THEN
	1275	!
	1276	!--- Limit cloud water sinks
	1277	!
	1278	DUM=DUM1/DUM2
	1279	PCOND=DUM*PCOND
	1280	PIACW=DUM*PIACW
	1281	PIACWI=DUM*PIACWI
	1282	ENDIF ! End IF (DUM2 .LT. DUM1)
	1283	ENDIF ! End IF (PIACWI.EQ.0. .AND. PIDEP.EQ.0.)
	1284	ENDIF ! End IF (TC.GE.T_ICE .AND. (QW.GT.EPSQ .OR. WV.GT.QSWgrd))
	1285	!
	1286	!--- Limit freezing of accreted rime to prevent temperature oscillations,
	1287	! a crude Schumann-Ludlam limit (p. 209 of Young, 1993).
	1288	!
	1289	TCC=TC+XLV1PCOND+XLS1PIDEP+XLF1*PIACWI
	1290	IF (TCC .GT. 0.) THEN
	1291	PIACWI=0.
	1292	TCC=TC+XLV1PCOND+XLS1PIDEP
	1293	ENDIF
	1294	IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical) THEN
	1295	!
	1296	!--- Calculate melting and evaporation/condensation
	1297	! * Units: SFACTOR - s.5/m ; ABI - m2/s ; NLICE - m**-3 ;
	1298	! VENTIL - m**-2 ; VENTI1 - m ;
	1299	! VENTI2 - m2/s.5 ; CIEVP - /s
	1300	!
	1301	SFACTOR=VEL_INC*.5(RHO/(DIFFUSDIFFUSDYNVIS))**C2
	1302	VENTIL=NLICE(VENTI1(INDEXS)+SFACTORVENTI2(INDEXS))
	1303	AIEVP=VENTILDIFFUSDTPH
	1304	IF (AIEVP .LT. Xratio) THEN
	1305	DIEVP=AIEVP
	1306	ELSE
	1307	DIEVP=1.-EXP(-AIEVP)
	1308	ENDIF
	1309	QSW0=EPS*ESW0/(PP-ESW0)
	1310	DWV0=MIN(WV,QSW)-QSW0
	1311	DUM=QW+PCOND
	1312	IF (WV.LT.QSW .AND. DUM.LE.EPSQ) THEN
	1313	!
	1314	!--- Evaporation from melting snow (sink of snow) or shedding
	1315	! of water condensed onto melting snow (source of rain)
	1316	!
	1317	DUM=DWV0*DIEVP
	1318	PIEVP=MAX( MIN(0., DUM), PILOSS)
	1319	PICND=MAX(0., DUM)
	1320	ENDIF ! End IF (WV.LT.QSW .AND. DUM.LE.EPSQ)
	1321	PIMLT=THERM_CONDTCCVENTILRRHODTPH/XLF
	1322	!
	1323	!--- Limit melting to prevent temperature oscillations across 0C
	1324	!
	1325	DUM1=MAX( 0., (TCC+XLV1*PIEVP)/XLF1 )
	1326	PIMLT=MIN(PIMLT, DUM1)
	1327	!
	1328	!--- Limit loss of snow by melting (>0) and evaporation
	1329	!
	1330	DUM=PIEVP-PIMLT
	1331	IF (DUM .LT. PILOSS) THEN
	1332	DUM1=PILOSS/DUM
	1333	PIMLT=PIMLT*DUM1
	1334	PIEVP=PIEVP*DUM1
	1335	ENDIF ! End IF (DUM .GT. QTICE)
	1336	ENDIF ! End IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical)
	1337	!
	1338	!--- IMPORTANT: Estimate time-averaged properties.
	1339	!
	1340	! * TOT_RAIN - total mass of rain before microphysics, which is the sum of
	1341	! the total mass of rain in the current layer and the input
	1342	! flux of rain from above
	1343	! * VRAIN1 - fall speed of rain into grid from above (with air resistance correction)
	1344	! * QTRAIN - time-averaged mixing ratio of rain (kg/kg)
	1345	! * PRLOSS - greatest loss (<0) of rain, removing all rain falling from
	1346	! above and the rain within the layer
	1347	! * RQR - rain content (kg/m**3)
	1348	! * INDEXR - mean size of rain drops to the nearest 1 micron in size
	1349	! * N0r - intercept of rain size distribution (typically 106 m-4)
	1350	!
	1351	TOT_RAIN=0.
	1352	VRAIN1=0.
	1353	QTRAIN=0.
	1354	PRLOSS=0.
	1355	RQR=0.
	1356	N0r=0.
	1357	INDEXR=MDRmin
	1358	INDEXR1=INDEXR !-- For debugging only
	1359	IF (RAIN_logical) THEN
	1360	IF (ARAIN .LE. 0.) THEN
	1361	INDEXR=MDRmin
	1362	VRAIN1=0.
	1363	ELSE
	1364	!
	1365	!--- INDEXR (related to mean diameter) & N0r could be modified
	1366	! by land/sea properties, presence of convection, etc.
	1367	!
	1368	!--- Rain rate normalized to a density of 1.194 kg/m**3
	1369	!
	1370	RR=ARAIN/(DTPH*GAMMAR)
	1371	!
	1372	IF (RR .LE. RR_DRmin) THEN
	1373	!
	1374	!--- Assume fixed mean diameter of rain (0.2 mm) for low rain rates,
	1375	! instead vary N0r with rain rate
	1376	!
	1377	INDEXR=MDRmin
	1378	ELSE IF (RR .LE. RR_DR1) THEN
	1379	!
	1380	!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
	1381	! for mean diameters (Dr) between 0.05 and 0.10 mm:
	1382	! V(Dr)=5.6023e4Dr*1.136, V in m/s and Dr in m
	1383	! RR = PI1000.N0r05.6023e4Dr*(4+1.136) = 1.408e15Dr**5.136,
	1384	! RR in kg/(m*2s)
	1385	! Dr (m) = 1.123e-3RR.1947 -> Dr (microns) = 1.123e3RR**.1947
	1386	!
	1387	INDEXR=INT( 1.123E3RR*.1947 + .5 )
	1388	INDEXR=MAX( MDRmin, MIN(INDEXR, MDR1) )
	1389	ELSE IF (RR .LE. RR_DR2) THEN
	1390	!
	1391	!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
	1392	! for mean diameters (Dr) between 0.10 and 0.20 mm:
	1393	! V(Dr)=1.0867e4Dr*.958, V in m/s and Dr in m
	1394	! RR = PI1000.N0r01.0867e4Dr*(4+.958) = 2.731e14Dr**4.958,
	1395	! RR in kg/(m*2s)
	1396	! Dr (m) = 1.225e-3RR.2017 -> Dr (microns) = 1.225e3RR**.2017
	1397	!
	1398	INDEXR=INT( 1.225E3RR*.2017 + .5 )
	1399	INDEXR=MAX( MDR1, MIN(INDEXR, MDR2) )
	1400	ELSE IF (RR .LE. RR_DR3) THEN
	1401	!
	1402	!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
	1403	! for mean diameters (Dr) between 0.20 and 0.32 mm:
	1404	! V(Dr)=2831.Dr*.80, V in m/s and Dr in m
	1405	! RR = PI1000.N0r02831.Dr*(4+.80) = 7.115e13Dr**4.80,
	1406	! RR in kg/(m*2s)
	1407	! Dr (m) = 1.3006e-3RR.2083 -> Dr (microns) = 1.3006e3RR**.2083
	1408	!
	1409	INDEXR=INT( 1.3006E3RR*.2083 + .5 )
	1410	INDEXR=MAX( MDR2, MIN(INDEXR, MDR3) )
	1411	ELSE IF (RR .LE. RR_DRmax) THEN
	1412	!
	1413	!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
	1414	! for mean diameters (Dr) between 0.32 and 0.45 mm:
	1415	! V(Dr)=944.8Dr*.6636, V in m/s and Dr in m
	1416	! RR = PI1000.N0r0944.8Dr*(4+.6636) = 2.3745e13Dr**4.6636,
	1417	! RR in kg/(m*2s)
	1418	! Dr (m) = 1.355e-3RR.2144 -> Dr (microns) = 1.355e3RR**.2144
	1419	!
	1420	INDEXR=INT( 1.355E3RR*.2144 + .5 )
	1421	INDEXR=MAX( MDR3, MIN(INDEXR, MDRmax) )
	1422	ELSE
	1423	!
	1424	!--- Assume fixed mean diameter of rain (0.45 mm) for high rain rates,
	1425	! instead vary N0r with rain rate
	1426	!
	1427	INDEXR=MDRmax
	1428	ENDIF ! End IF (RR .LE. RR_DRmin) etc.
	1429	VRAIN1=GAMMAR*VRAIN(INDEXR)
	1430	ENDIF ! End IF (ARAIN .LE. 0.)
	1431	INDEXR1=INDEXR ! For debugging only
	1432	TOT_RAIN=THICKQR+BLENDARAIN
	1433	QTRAIN=TOT_RAIN/(THICK+BLDTRH*VRAIN1)
	1434	PRLOSS=-TOT_RAIN/THICK
	1435	RQR=RHO*QTRAIN
	1436	!
	1437	!--- RQR - time-averaged rain content (kg/m**3)
	1438	!
	1439	IF (RQR .LE. RQR_DRmin) THEN
	1440	N0r=MAX(N0rmin, CN0r_DMRmin*RQR)
	1441	INDEXR=MDRmin
	1442	ELSE IF (RQR .GE. RQR_DRmax) THEN
	1443	N0r=CN0r_DMRmax*RQR
	1444	INDEXR=MDRmax
	1445	ELSE
	1446	N0r=N0r0
	1447	INDEXR=MAX( XMRmin, MIN(CN0r0RQR*.25, XMRmax) )
	1448	ENDIF
	1449	!
	1450	IF (TC .LT. T_ICE) THEN
	1451	PIACR=-PRLOSS
	1452	ELSE
	1453	DWVr=WV-PCOND-QSW
	1454	DUM=QW+PCOND
	1455	IF (DWVr.LT.0. .AND. DUM.LE.EPSQ) THEN
	1456	!
	1457	!--- Rain evaporation
	1458	!
	1459	! * RFACTOR - [GAMMAR*.5][Schmidt*(1./3.)][(RHO/DYNVIS)**.5],
	1460	! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
	1461	!
	1462	! * Units: RFACTOR - s.5/m ; ABW - m2/s ; VENTR - m**-2 ;
	1463	! N0r - m-4 ; VENTR1 - m2 ; VENTR2 - m3/s.5 ;
	1464	! CREVP - unitless
	1465	!
	1466	RFACTOR=GAMMAR*.5(RHO/(DIFFUSDIFFUSDYNVIS))**C2
	1467	ABW=1./(RHO*XLV2/THERM_COND+1./DIFFUS)
	1468	!
	1469	!--- Note that VENTR1, VENTR2 lookup tables do not include the
	1470	! 1/Davg multiplier as in the ice tables
	1471	!
	1472	VENTR=N0r(VENTR1(INDEXR)+RFACTORVENTR2(INDEXR))
	1473	CREVP=ABWVENTRDTPH
	1474	IF (CREVP .LT. Xratio) THEN
	1475	DUM=DWVr*CREVP
	1476	ELSE
	1477	DUM=DWVr(1.-EXP(-CREVPDENOMW))/DENOMW
	1478	ENDIF
	1479	PREVP=MAX(DUM, PRLOSS)
	1480	ELSE IF (QW .GT. EPSQ) THEN
	1481	FWR=CRACWGAMMARN0r*ACCRR(INDEXR)
	1482	PRACW=MIN(1.,FWR)*QW
	1483	ENDIF ! End IF (DWVr.LT.0. .AND. DUM.LE.EPSQ)
	1484	!
	1485	IF (TC.LT.0. .AND. TCC.LT.0.) THEN
	1486	!
	1487	!--- Biggs (1953) heteorogeneous freezing (e.g., Lin et al., 1983)
	1488	! - Rescaled mean drop diameter from microns (INDEXR) to mm (DUM) to prevent underflow
	1489	!
	1490	DUM=.001*FLOAT(INDEXR)
	1491	DUM=(EXP(ABFRTC)-1.)DUMDUMDUMDUMDUMDUMDUM
	1492	PIACR=MIN(CBFRN0rRRHO*DUM, QTRAIN)
	1493	IF (QLICE .GT. EPSQ) THEN
	1494	!
	1495	!--- Freezing of rain by collisions w/ large ice
	1496	!
	1497	DUM=GAMMAR*VRAIN(INDEXR)
	1498	DUM1=DUM-VSNOW
	1499	!
	1500	!--- DUM2 - Difference in spectral fall speeds of rain and
	1501	! large ice, parameterized following eq. (48) on p. 112 of
	1502	! Murakami (J. Meteor. Soc. Japan, 1990)
	1503	!
	1504	DUM2=(DUM1DUM1+.04DUMVSNOW)*.5
	1505	DUM1=5.E-12INDEXRINDEXR+2.E-12INDEXRINDEXS &
	1506	& +.5E-12INDEXSINDEXS
	1507	FIR=MIN(1., CIACRNLICEDUM1*DUM2)
	1508	!
	1509	!--- Future? Should COLLECTION BY SMALL ICE SHOULD BE INCLUDED???
	1510	!
	1511	PIACR=MIN(PIACR+FIR*QTRAIN, QTRAIN)
	1512	ENDIF ! End IF (QLICE .GT. EPSQ)
	1513	DUM=PREVP-PIACR
	1514	If (DUM .LT. PRLOSS) THEN
	1515	DUM1=PRLOSS/DUM
	1516	PREVP=DUM1*PREVP
	1517	PIACR=DUM1*PIACR
	1518	ENDIF ! End If (DUM .LT. PRLOSS)
	1519	ENDIF ! End IF (TC.LT.0. .AND. TCC.LT.0.)
	1520	ENDIF ! End IF (TC .LT. T_ICE)
	1521	ENDIF ! End IF (RAIN_logical)
	1522	!
	1523	!----------------------------------------------------------------------
	1524	!---------------------- Main Budget Equations -------------------------
	1525	!----------------------------------------------------------------------
	1526	!
	1527	!
	1528	!-----------------------------------------------------------------------
	1529	!--- Update fields, determine characteristics for next lower layer ----
	1530	!-----------------------------------------------------------------------
	1531	!
	1532	!--- Carefully limit sinks of cloud water
	1533	!
	1534	DUM1=PIACW+PRAUT+PRACW-MIN(0.,PCOND)
	1535	IF (DUM1 .GT. QW) THEN
	1536	DUM=QW/DUM1
	1537	PIACW=DUM*PIACW
	1538	PIACWI=DUM*PIACWI
	1539	PRAUT=DUM*PRAUT
	1540	PRACW=DUM*PRACW
	1541	IF (PCOND .LT. 0.) PCOND=DUM*PCOND
	1542	ENDIF
	1543	PIACWR=PIACW-PIACWI ! TC >= 0C
	1544	!
	1545	!--- QWnew - updated cloud water mixing ratio
	1546	!
	1547	DELW=PCOND-PIACW-PRAUT-PRACW
	1548	QWnew=QW+DELW
	1549	IF (QWnew .LE. EPSQ) QWnew=0.
	1550	IF (QW.GT.0. .AND. QWnew.NE.0.) THEN
	1551	DUM=QWnew/QW
	1552	IF (DUM .LT. TOLER) QWnew=0.
	1553	ENDIF
	1554	!
	1555	!--- Update temperature and water vapor mixing ratios
	1556	!
	1557	DELT= XLV1*(PCOND+PIEVP+PICND+PREVP) &
	1558	& +XLS1PIDEP+XLF1(PIACWI+PIACR-PIMLT)
	1559	Tnew=TK+DELT
	1560	!
	1561	DELV=-PCOND-PIDEP-PIEVP-PICND-PREVP
	1562	WVnew=WV+DELV
	1563	!
	1564	!--- Update ice mixing ratios
	1565	!
	1566	!---
	1567	! * TOT_ICEnew - total mass (small & large) ice after microphysics,
	1568	! which is the sum of the total mass of large ice in the
	1569	! current layer and the flux of ice out of the grid box below
	1570	! * RimeF - Rime Factor, which is the mass ratio of total (unrimed &
	1571	! rimed) ice mass to the unrimed ice mass (>=1)
	1572	! * QInew - updated mixing ratio of total (large & small) ice in layer
	1573	! -> TOT_ICEnew=QInewTHICK+BLDTRHQLICEnew*VSNOW
	1574	! -> But QLICEnew=QInew*FLIMASS, so
	1575	! -> TOT_ICEnew=QInew(THICK+BLDTRHFLIMASS*VSNOW)
	1576	! * ASNOWnew - updated accumulation of snow at bottom of grid cell
	1577	!---
	1578	!
	1579	DELI=0.
	1580	RimeF=1.
	1581	IF (ICE_logical) THEN
	1582	DELI=PIDEP+PIEVP+PIACWI+PIACR-PIMLT
	1583	TOT_ICEnew=TOT_ICE+THICK*DELI
	1584	IF (TOT_ICE.GT.0. .AND. TOT_ICEnew.NE.0.) THEN
	1585	DUM=TOT_ICEnew/TOT_ICE
	1586	IF (DUM .LT. TOLER) TOT_ICEnew=0.
	1587	ENDIF
	1588	IF (TOT_ICEnew .LE. CLIMIT) THEN
	1589	TOT_ICEnew=0.
	1590	RimeF=1.
	1591	QInew=0.
	1592	ASNOWnew=0.
	1593	ELSE
	1594	!
	1595	!--- Update rime factor if appropriate
	1596	!
	1597	DUM=PIACWI+PIACR
	1598	IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ) THEN
	1599	RimeF=RimeF1
	1600	ELSE
	1601	!
	1602	!--- Rime Factor, RimeF = (Total ice mass)/(Total unrimed ice mass)
	1603	! DUM1 - Total ice mass, rimed & unrimed
	1604	! DUM2 - Estimated mass of unrimed ice
	1605	!
	1606	DUM1=TOT_ICE+THICK*(PIDEP+DUM)
	1607	DUM2=TOT_ICE/RimeF1+THICK*PIDEP
	1608	IF (DUM2 .LE. 0.) THEN
	1609	RimeF=RFmax
	1610	ELSE
	1611	RimeF=MIN(RFmax, MAX(1., DUM1/DUM2) )
	1612	ENDIF
	1613	ENDIF ! End IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ)
	1614	QInew=TOT_ICEnew/(THICK+BLDTRHFLIMASSVSNOW)
	1615	IF (QInew .LE. EPSQ) QInew=0.
	1616	IF (QI.GT.0. .AND. QInew.NE.0.) THEN
	1617	DUM=QInew/QI
	1618	IF (DUM .LT. TOLER) QInew=0.
	1619	ENDIF
	1620	ASNOWnew=BLDTRHFLIMASSVSNOW*QInew
	1621	IF (ASNOW.GT.0. .AND. ASNOWnew.NE.0.) THEN
	1622	DUM=ASNOWnew/ASNOW
	1623	IF (DUM .LT. TOLER) ASNOWnew=0.
	1624	ENDIF
	1625	ENDIF ! End IF (TOT_ICEnew .LE. CLIMIT)
	1626	ENDIF ! End IF (ICE_logical)
	1627
	1628
	1629	!
	1630	!--- Update rain mixing ratios
	1631	!
	1632	!---
	1633	! * TOT_RAINnew - total mass of rain after microphysics
	1634	! current layer and the input flux of ice from above
	1635	! * VRAIN2 - time-averaged fall speed of rain in grid and below
	1636	! (with air resistance correction)
	1637	! * QRnew - updated rain mixing ratio in layer
	1638	! -> TOT_RAINnew=QRnew(THICK+BLDTRHVRAIN2)
	1639	! * ARAINnew - updated accumulation of rain at bottom of grid cell
	1640	!---
	1641	!
	1642	DELR=PRAUT+PRACW+PIACWR-PIACR+PIMLT+PREVP+PICND
	1643	TOT_RAINnew=TOT_RAIN+THICK*DELR
	1644	IF (TOT_RAIN.GT.0. .AND. TOT_RAINnew.NE.0.) THEN
	1645	DUM=TOT_RAINnew/TOT_RAIN
	1646	IF (DUM .LT. TOLER) TOT_RAINnew=0.
	1647	ENDIF
	1648	IF (TOT_RAINnew .LE. CLIMIT) THEN
	1649	TOT_RAINnew=0.
	1650	VRAIN2=0.
	1651	QRnew=0.
	1652	ARAINnew=0.
	1653	ELSE
	1654	!
	1655	!--- 1st guess time-averaged rain rate at bottom of grid box
	1656	!
	1657	RR=TOT_RAINnew/(DTPH*GAMMAR)
	1658	!
	1659	!--- Use same algorithm as above for calculating mean drop diameter
	1660	! (IDR, in microns), which is used to estimate the time-averaged
	1661	! fall speed of rain drops at the bottom of the grid layer. This
	1662	! isn't perfect, but the alternative is solving a transcendental
	1663	! equation that is numerically inefficient and nasty to program
	1664	! (coded in earlier versions of GSMCOLUMN prior to 8-22-01).
	1665	!
	1666	IF (RR .LE. RR_DRmin) THEN
	1667	IDR=MDRmin
	1668	ELSE IF (RR .LE. RR_DR1) THEN
	1669	IDR=INT( 1.123E3RR*.1947 + .5 )
	1670	IDR=MAX( MDRmin, MIN(IDR, MDR1) )
	1671	ELSE IF (RR .LE. RR_DR2) THEN
	1672	IDR=INT( 1.225E3RR*.2017 + .5 )
	1673	IDR=MAX( MDR1, MIN(IDR, MDR2) )
	1674	ELSE IF (RR .LE. RR_DR3) THEN
	1675	IDR=INT( 1.3006E3RR*.2083 + .5 )
	1676	IDR=MAX( MDR2, MIN(IDR, MDR3) )
	1677	ELSE IF (RR .LE. RR_DRmax) THEN
	1678	IDR=INT( 1.355E3RR*.2144 + .5 )
	1679	IDR=MAX( MDR3, MIN(IDR, MDRmax) )
	1680	ELSE
	1681	IDR=MDRmax
	1682	ENDIF ! End IF (RR .LE. RR_DRmin)
	1683	VRAIN2=GAMMAR*VRAIN(IDR)
	1684	QRnew=TOT_RAINnew/(THICK+BLDTRH*VRAIN2)
	1685	IF (QRnew .LE. EPSQ) QRnew=0.
	1686	IF (QR.GT.0. .AND. QRnew.NE.0.) THEN
	1687	DUM=QRnew/QR
	1688	IF (DUM .LT. TOLER) QRnew=0.
	1689	ENDIF
	1690	ARAINnew=BLDTRHVRAIN2QRnew
	1691	IF (ARAIN.GT.0. .AND. ARAINnew.NE.0.) THEN
	1692	DUM=ARAINnew/ARAIN
	1693	IF (DUM .LT. TOLER) ARAINnew=0.
	1694	ENDIF
	1695	ENDIF
	1696	!
	1697	WCnew=QWnew+QRnew+QInew
	1698	!
	1699	!----------------------------------------------------------------------
	1700	!-------------- Begin debugging & verification ------------------------
	1701	!----------------------------------------------------------------------
	1702	!
	1703	!--- QT, QTnew - total water (vapor & condensate) before & after microphysics, resp.
	1704	!
	1705
	1706
	1707	QT=THICK*(WV+WC)+ARAIN+ASNOW
	1708	QTnew=THICK*(WVnew+WCnew)+ARAINnew+ASNOWnew
	1709	BUDGET=QT-QTnew
	1710	!
	1711	!--- Additional check on budget preservation, accounting for truncation effects
	1712	!
	1713	DBG_logical=.FALSE.
	1714	! DUM=ABS(BUDGET)
	1715	! IF (DUM .GT. TOLER) THEN
	1716	! DUM=DUM/MIN(QT, QTnew)
	1717	! IF (DUM .GT. TOLER) DBG_logical=.TRUE.
	1718	! ENDIF
	1719	!!
	1720	! DUM=(RHgrd+.001)*QSInew
	1721	! IF ( (QWnew.GT.EPSQ) .OR. QRnew.GT.EPSQ .OR. WVnew.GT.DUM)
	1722	! & .AND. TC.LT.T_ICE ) DBG_logical=.TRUE.
	1723	!
	1724	! IF (TC.GT.5. .AND. QInew.GT.EPSQ) DBG_logical=.TRUE.
	1725	!
	1726	IF ((WVnew.LT.EPSQ .OR. DBG_logical) .AND. PRINT_diag) THEN
	1727	!
	1728	WRITE(6,"(/2(a,i4),2(a,i2))") '{} i=',I_index,' j=',J_index,&
	1729	& ' L=',L,' LSFC=',LSFC
	1730	!
	1731	ESW=1000.*FPVS0(Tnew)
	1732	QSWnew=EPS*ESW/(PP-ESW)
	1733	IF (TC.LT.0. .OR. Tnew .LT. 0.) THEN
	1734	ESI=1000.*FPVS(Tnew)
	1735	QSInew=EPS*ESI/(PP-ESI)
	1736	ELSE
	1737	QSI=QSW
	1738	QSInew=QSWnew
	1739	ENDIF
	1740	WSnew=QSInew
	1741	WRITE(6,"(4(a12,g11.4,1x))") &
	1742	& '{} TCold=',TC,'TCnew=',Tnew-T0C,'P=',.01*PP,'RHO=',RHO, &
	1743	& '{} THICK=',THICK,'RHold=',WV/WS,'RHnew=',WVnew/WSnew, &
	1744	& 'RHgrd=',RHgrd, &
	1745	& '{} RHWold=',WV/QSW,'RHWnew=',WVnew/QSWnew,'RHIold=',WV/QSI, &
	1746	& 'RHInew=',WVnew/QSInew, &
	1747	& '{} QSWold=',QSW,'QSWnew=',QSWnew,'QSIold=',QSI,'QSInew=',QSInew, &
	1748	& '{} WSold=',WS,'WSnew=',WSnew,'WVold=',WV,'WVnew=',WVnew, &
	1749	& '{} WCold=',WC,'WCnew=',WCnew,'QWold=',QW,'QWnew=',QWnew, &
	1750	& '{} QIold=',QI,'QInew=',QInew,'QRold=',QR,'QRnew=',QRnew, &
	1751	& '{} ARAINold=',ARAIN,'ARAINnew=',ARAINnew,'ASNOWold=',ASNOW, &
	1752	& 'ASNOWnew=',ASNOWnew, &
	1753	& '{} TOT_RAIN=',TOT_RAIN,'TOT_RAINnew=',TOT_RAINnew, &
	1754	& 'TOT_ICE=',TOT_ICE,'TOT_ICEnew=',TOT_ICEnew, &
	1755	& '{} BUDGET=',BUDGET,'QTold=',QT,'QTnew=',QTnew
	1756	!
	1757	WRITE(6,"(4(a12,g11.4,1x))") &
	1758	& '{} DELT=',DELT,'DELV=',DELV,'DELW=',DELW,'DELI=',DELI, &
	1759	& '{} DELR=',DELR,'PCOND=',PCOND,'PIDEP=',PIDEP,'PIEVP=',PIEVP, &
	1760	& '{} PICND=',PICND,'PREVP=',PREVP,'PRAUT=',PRAUT,'PRACW=',PRACW, &
	1761	& '{} PIACW=',PIACW,'PIACWI=',PIACWI,'PIACWR=',PIACWR,'PIMLT=', &
	1762	& PIMLT, &
	1763	& '{} PIACR=',PIACR
	1764	!
	1765	IF (ICE_logical) WRITE(6,"(4(a12,g11.4,1x))") &
	1766	& '{} RimeF1=',RimeF1,'GAMMAS=',GAMMAS,'VrimeF=',VrimeF, &
	1767	& 'VSNOW=',VSNOW, &
	1768	& '{} INDEXS=',FLOAT(INDEXS),'FLARGE=',FLARGE,'FSMALL=',FSMALL, &
	1769	& 'FLIMASS=',FLIMASS, &
	1770	& '{} XSIMASS=',XSIMASS,'XLIMASS=',XLIMASS,'QLICE=',QLICE, &
	1771	& 'QTICE=',QTICE, &
	1772	& '{} NLICE=',NLICE,'NSmICE=',NSmICE,'PILOSS=',PILOSS, &
	1773	& 'EMAIRI=',EMAIRI, &
	1774	& '{} RimeF=',RimeF
	1775	!
	1776	IF (TOT_RAIN.GT.0. .OR. TOT_RAINnew.GT.0.) &
	1777	& WRITE(6,"(4(a12,g11.4,1x))") &
	1778	& '{} INDEXR1=',FLOAT(INDEXR1),'INDEXR=',FLOAT(INDEXR), &
	1779	& 'GAMMAR=',GAMMAR,'N0r=',N0r, &
	1780	& '{} VRAIN1=',VRAIN1,'VRAIN2=',VRAIN2,'QTRAIN=',QTRAIN,'RQR=',RQR, &
	1781	& '{} PRLOSS=',PRLOSS,'VOLR1=',THICK+BLDTRH*VRAIN1, &
	1782	& 'VOLR2=',THICK+BLDTRH*VRAIN2
	1783	!
	1784	IF (PRAUT .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} QW0=',QW0
	1785	!
	1786	IF (PRACW .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} FWR=',FWR
	1787	!
	1788	IF (PIACR .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} FIR=',FIR
	1789	!
	1790	DUM=PIMLT+PICND-PREVP-PIEVP
	1791	IF (DUM.GT.0. .or. DWVi.NE.0.) &
	1792	& WRITE(6,"(4(a12,g11.4,1x))") &
	1793	& '{} TFACTOR=',TFACTOR,'DYNVIS=',DYNVIS, &
	1794	& 'THERM_CON=',THERM_COND,'DIFFUS=',DIFFUS
	1795	!
	1796	IF (PREVP .LT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
	1797	& '{} RFACTOR=',RFACTOR,'ABW=',ABW,'VENTR=',VENTR,'CREVP=',CREVP, &
	1798	& '{} DWVr=',DWVr,'DENOMW=',DENOMW
	1799	!
	1800	IF (PIDEP.NE.0. .AND. DWVi.NE.0.) &
	1801	& WRITE(6,"(4(a12,g11.4,1x))") &
	1802	& '{} DWVi=',DWVi,'DENOMI=',DENOMI,'PIDEP_max=',PIDEP_max, &
	1803	& 'SFACTOR=',SFACTOR, &
	1804	& '{} ABI=',ABI,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
	1805	& 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
	1806	& '{} VENTIS=',VENTIS,'DIDEP=',DIDEP
	1807	!
	1808	IF (PIDEP.GT.0. .AND. PCOND.NE.0.) &
	1809	& WRITE(6,"(4(a12,g11.4,1x))") &
	1810	& '{} DENOMW=',DENOMW,'DENOMWI=',DENOMWI,'DENOMF=',DENOMF, &
	1811	& 'DUM2=',PCOND-PIACW
	1812	!
	1813	IF (FWS .GT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
	1814	& '{} FWS=',FWS
	1815	!
	1816	DUM=PIMLT+PICND-PIEVP
	1817	IF (DUM.GT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
	1818	& '{} SFACTOR=',SFACTOR,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
	1819	& 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
	1820	& '{} AIEVP=',AIEVP,'DIEVP=',DIEVP,'QSW0=',QSW0,'DWV0=',DWV0
	1821	!
	1822	ENDIF
	1823
	1824
	1825	!
	1826	!-----------------------------------------------------------------------
	1827	!--------------- Water budget statistics & maximum values --------------
	1828	!-----------------------------------------------------------------------
	1829	!
	1830	IF (PRINT_diag) THEN
	1831	ITdx=MAX( ITLO, MIN( INT(Tnew-T0C), ITHI ) )
	1832	IF (QInew .GT. EPSQ) NSTATS(ITdx,1)=NSTATS(ITdx,1)+1
	1833	IF (QInew.GT.EPSQ .AND. QRnew+QWnew.GT.EPSQ) &
	1834	& NSTATS(ITdx,2)=NSTATS(ITdx,2)+1
	1835	IF (QWnew .GT. EPSQ) NSTATS(ITdx,3)=NSTATS(ITdx,3)+1
	1836	IF (QRnew .GT. EPSQ) NSTATS(ITdx,4)=NSTATS(ITdx,4)+1
	1837	!
	1838	QMAX(ITdx,1)=MAX(QMAX(ITdx,1), QInew)
	1839	QMAX(ITdx,2)=MAX(QMAX(ITdx,2), QWnew)
	1840	QMAX(ITdx,3)=MAX(QMAX(ITdx,3), QRnew)
	1841	QMAX(ITdx,4)=MAX(QMAX(ITdx,4), ASNOWnew)
	1842	QMAX(ITdx,5)=MAX(QMAX(ITdx,5), ARAINnew)
	1843	QTOT(ITdx,1)=QTOT(ITdx,1)+QInew*THICK
	1844	QTOT(ITdx,2)=QTOT(ITdx,2)+QWnew*THICK
	1845	QTOT(ITdx,3)=QTOT(ITdx,3)+QRnew*THICK
	1846	!
	1847	QTOT(ITdx,4)=QTOT(ITdx,4)+PCOND*THICK
	1848	QTOT(ITdx,5)=QTOT(ITdx,5)+PICND*THICK
	1849	QTOT(ITdx,6)=QTOT(ITdx,6)+PIEVP*THICK
	1850	QTOT(ITdx,7)=QTOT(ITdx,7)+PIDEP*THICK
	1851	QTOT(ITdx,8)=QTOT(ITdx,8)+PREVP*THICK
	1852	QTOT(ITdx,9)=QTOT(ITdx,9)+PRAUT*THICK
	1853	QTOT(ITdx,10)=QTOT(ITdx,10)+PRACW*THICK
	1854	QTOT(ITdx,11)=QTOT(ITdx,11)+PIMLT*THICK
	1855	QTOT(ITdx,12)=QTOT(ITdx,12)+PIACW*THICK
	1856	QTOT(ITdx,13)=QTOT(ITdx,13)+PIACWI*THICK
	1857	QTOT(ITdx,14)=QTOT(ITdx,14)+PIACWR*THICK
	1858	QTOT(ITdx,15)=QTOT(ITdx,15)+PIACR*THICK
	1859	!
	1860	QTOT(ITdx,16)=QTOT(ITdx,16)+(WVnew-WV)*THICK
	1861	QTOT(ITdx,17)=QTOT(ITdx,17)+(QWnew-QW)*THICK
	1862	QTOT(ITdx,18)=QTOT(ITdx,18)+(QInew-QI)*THICK
	1863	QTOT(ITdx,19)=QTOT(ITdx,19)+(QRnew-QR)*THICK
	1864	QTOT(ITdx,20)=QTOT(ITdx,20)+(ARAINnew-ARAIN)
	1865	QTOT(ITdx,21)=QTOT(ITdx,21)+(ASNOWnew-ASNOW)
	1866	IF (QInew .GT. 0.) &
	1867	& QTOT(ITdx,22)=QTOT(ITdx,22)+QInew*THICK/RimeF
	1868	!
	1869	ENDIF
	1870	!
	1871	!----------------------------------------------------------------------
	1872	!------------------------- Update arrays ------------------------------
	1873	!----------------------------------------------------------------------
	1874	!
	1875
	1876
	1877	T_col(L)=Tnew ! Updated temperature
	1878	!
	1879	QV_col(L)=max(EPSQ, WVnew/(1.+WVnew)) ! Updated specific humidity
	1880	WC_col(L)=max(EPSQ, WCnew) ! Updated total condensate mixing ratio
	1881	QI_col(L)=max(EPSQ, QInew) ! Updated ice mixing ratio
	1882	QR_col(L)=max(EPSQ, QRnew) ! Updated rain mixing ratio
	1883	QW_col(L)=max(EPSQ, QWnew) ! Updated cloud water mixing ratio
	1884	RimeF_col(L)=RimeF ! Updated rime factor
	1885	ASNOW=ASNOWnew ! Updated accumulated snow
	1886	ARAIN=ARAINnew ! Updated accumulated rain
	1887	!
	1888	!#######################################################################
	1889	!
	1890	10 CONTINUE ! ##### End "L" loop through model levels #####
	1891
	1892
	1893	!
	1894	!#######################################################################
	1895	!
	1896	!-----------------------------------------------------------------------
	1897	!--------------------------- Return to GSMDRIVE -----------------------
	1898	!-----------------------------------------------------------------------
	1899	!
	1900	CONTAINS
	1901	!#######################################################################
	1902	!--------- Produces accurate calculation of cloud condensation ---------
	1903	!#######################################################################
	1904	!
	1905	REAL FUNCTION CONDENSE (PP, QW, TK, WV)
	1906	!
	1907	!---------------------------------------------------------------------------------
	1908	!------ The Asai (1965) algorithm takes into consideration the release of ------
	1909	!------ latent heat in increasing the temperature & in increasing the ------
	1910	!------ saturation mixing ratio (following the Clausius-Clapeyron eqn.). ------
	1911	!---------------------------------------------------------------------------------
	1912	!
	1913	IMPLICIT NONE
	1914	!
	1915	INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
	1916	REAL (KIND=HIGH_PRES), PARAMETER :: &
	1917	& RHLIMIT=.001, RHLIMIT1=-RHLIMIT
	1918	REAL (KIND=HIGH_PRES) :: COND, SSAT, WCdum
	1919	!
	1920	REAL,INTENT(IN) :: QW,PP,WV,TK
	1921	REAL WVdum,Tdum,XLV2,DWV,WS,ESW,XLV1,XLV
	1922	integer nsteps
	1923	!
	1924	!-----------------------------------------------------------------------
	1925	!
	1926	!--- LV (T) is from Bolton (JAS, 1980)
	1927	!
	1928	XLV=3.148E6-2370.*TK
	1929	XLV1=XLV*RCP
	1930	XLV2=XLVXLVRCPRV
	1931	Tdum=TK
	1932	WVdum=WV
	1933	WCdum=QW
	1934	ESW=1000.*FPVS0(Tdum) ! Saturation vapor press w/r/t water
	1935	WS=RHgrdEPSESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
	1936	DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
	1937	SSAT=DWV/WS ! Supersaturation ratio
	1938	CONDENSE=0.
	1939	nsteps = 0
	1940	DO WHILE ((SSAT.LT.RHLIMIT1 .AND. WCdum.GT.EPSQ) &
	1941	& .OR. SSAT.GT.RHLIMIT)
	1942	nsteps = nsteps + 1
	1943	COND=DWV/(1.+XLV2WS/(TdumTdum)) ! Asai (1965, J. Japan)
	1944	COND=MAX(COND, -WCdum) ! Limit cloud water evaporation
	1945	Tdum=Tdum+XLV1*COND ! Updated temperature
	1946	WVdum=WVdum-COND ! Updated water vapor mixing ratio
	1947	WCdum=WCdum+COND ! Updated cloud water mixing ratio
	1948	CONDENSE=CONDENSE+COND ! Total cloud water condensation
	1949	ESW=1000.*FPVS0(Tdum) ! Updated saturation vapor press w/r/t water
	1950	WS=RHgrdEPSESW/(PP-ESW) ! Updated saturation mixing ratio w/r/t water
	1951	DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
	1952	SSAT=DWV/WS ! Grid-scale supersaturation ratio
	1953	ENDDO
	1954	!
	1955	END FUNCTION CONDENSE
	1956	!
	1957	!#######################################################################
	1958	!---------------- Calculate ice deposition at T<T_ICE ------------------
	1959	!#######################################################################
	1960	!
	1961	REAL FUNCTION DEPOSIT (PP, Tdum, WVdum)
	1962	!
	1963	!--- Also uses the Asai (1965) algorithm, but uses a different target
	1964	! vapor pressure for the adjustment
	1965	!
	1966	IMPLICIT NONE
	1967	!
	1968	INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
	1969	REAL (KIND=HIGH_PRES), PARAMETER :: RHLIMIT=.001, &
	1970	& RHLIMIT1=-RHLIMIT
	1971	REAL (KIND=HIGH_PRES) :: DEP, SSAT
	1972	!
	1973	real,INTENT(IN) :: PP
	1974	real,INTENT(INOUT) :: WVdum,Tdum
	1975	real ESI,WS,DWV
	1976	!
	1977	!-----------------------------------------------------------------------
	1978	!
	1979	ESI=1000.*FPVS(Tdum) ! Saturation vapor press w/r/t ice
	1980	WS=RHgrdEPSESI/(PP-ESI) ! Saturation mixing ratio
	1981	DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
	1982	SSAT=DWV/WS ! Supersaturation ratio
	1983	DEPOSIT=0.
	1984	DO WHILE (SSAT.GT.RHLIMIT .OR. SSAT.LT.RHLIMIT1)
	1985	!
	1986	!--- Note that XLVS2=LSLV/(CPRV)=LVWS/(RVTT)(LS/CP*DEP1),
	1987	! where WS is the saturation mixing ratio following Clausius-
	1988	! Clapeyron (see Asai,1965; Young,1993,p.405)
	1989	!
	1990	DEP=DWV/(1.+XLS2WS/(TdumTdum)) ! Asai (1965, J. Japan)
	1991	Tdum=Tdum+XLS1*DEP ! Updated temperature
	1992	WVdum=WVdum-DEP ! Updated ice mixing ratio
	1993	DEPOSIT=DEPOSIT+DEP ! Total ice deposition
	1994	ESI=1000.*FPVS(Tdum) ! Updated saturation vapor press w/r/t ice
	1995	WS=RHgrdEPSESI/(PP-ESI) ! Updated saturation mixing ratio w/r/t ice
	1996	DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
	1997	SSAT=DWV/WS ! Grid-scale supersaturation ratio
	1998	ENDDO
	1999	!
	2000	END FUNCTION DEPOSIT
	2001	!
	2002	END SUBROUTINE EGCP01COLUMN
	2003
	2004	!#######################################################################
	2005	!------- Initialize constants & lookup tables for microphysics ---------
	2006	!#######################################################################
	2007	!
	2008
	2009	! SH 0211/2002
	2010
	2011	!-----------------------------------------------------------------------
	2012	SUBROUTINE etanewinit (GSMDT,DT,DELX,DELY,LOWLYR,restart, &
	2013	& F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
	2014	& MP_RESTART_STATE,TBPVS_STATE,TBPVS0_STATE, &
	2015	& ALLOWED_TO_READ, &
	2016	& IDS,IDE,JDS,JDE,KDS,KDE, &
	2017	& IMS,IME,JMS,JME,KMS,KME, &
	2018	& ITS,ITE,JTS,JTE,KTS,KTE )
	2019	!-----------------------------------------------------------------------
	2020	!-------------------------------------------------------------------------------
	2021	!--- SUBPROGRAM DOCUMENTATION BLOCK
	2022	! PRGRMMR: Ferrier ORG: W/NP22 DATE: February 2001
	2023	! Jin ORG: W/NP22 DATE: January 2002
	2024	! (Modification for WRF structure)
	2025	!
	2026	!-------------------------------------------------------------------------------
	2027	! ABSTRACT:
	2028	! * Reads various microphysical lookup tables used in COLUMN_MICRO
	2029	! * Lookup tables were created "offline" and are read in during execution
	2030	! * Creates lookup tables for saturation vapor pressure w/r/t water & ice
	2031	!-------------------------------------------------------------------------------
	2032	!
	2033	! USAGE: CALL etanewinit FROM SUBROUTINE GSMDRIVE AT MODEL START TIME
	2034	!
	2035	! INPUT ARGUMENT LIST:
	2036	! DTPH - physics time step (s)
	2037	!
	2038	! OUTPUT ARGUMENT LIST:
	2039	! NONE
	2040	!
	2041	! OUTPUT FILES:
	2042	! NONE
	2043	!
	2044	! SUBROUTINES:
	2045	! MY_GROWTH_RATES - lookup table for growth of nucleated ice
	2046	! GPVS - lookup tables for saturation vapor pressure (water, ice)
	2047	!
	2048	! UNIQUE: NONE
	2049	!
	2050	! LIBRARY: NONE
	2051	!
	2052	! COMMON BLOCKS:
	2053	! CMICRO_CONS - constants used in GSMCOLUMN
	2054	! CMY600 - lookup table for growth of ice crystals in
	2055	! water saturated conditions (Miller & Young, 1979)
	2056	! IVENT_TABLES - lookup tables for ventilation effects of ice
	2057	! IACCR_TABLES - lookup tables for accretion rates of ice
	2058	! IMASS_TABLES - lookup tables for mass content of ice
	2059	! IRATE_TABLES - lookup tables for precipitation rates of ice
	2060	! IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
	2061	! MAPOT - Need lat/lon grid resolution
	2062	! RVENT_TABLES - lookup tables for ventilation effects of rain
	2063	! RACCR_TABLES - lookup tables for accretion rates of rain
	2064	! RMASS_TABLES - lookup tables for mass content of rain
	2065	! RVELR_TABLES - lookup tables for fall speeds of rain
	2066	! RRATE_TABLES - lookup tables for precipitation rates of rain
	2067	!
	2068	! ATTRIBUTES:
	2069	! LANGUAGE: FORTRAN 90
	2070	! MACHINE : IBM SP
	2071	!
	2072	!-----------------------------------------------------------------------
	2073	!
	2074	!
	2075	!-----------------------------------------------------------------------
	2076	IMPLICIT NONE
	2077	!-----------------------------------------------------------------------
	2078	!-------------------------------------------------------------------------
	2079	!-------------- Parameters & arrays for lookup tables --------------------
	2080	!-------------------------------------------------------------------------
	2081	!
	2082	!--- Common block of constants used in column microphysics
	2083	!
	2084	!WRF
	2085	! real DLMD,DPHD
	2086	!WRF
	2087	!
	2088	!-----------------------------------------------------------------------
	2089	!--- Parameters & data statement for local calculations
	2090	!-----------------------------------------------------------------------
	2091	!
	2092	INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3
	2093	!
	2094	! VARIABLES PASSED IN
	2095	integer,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
	2096	& ,IMS,IME,JMS,JME,KMS,KME &
	2097	& ,ITS,ITE,JTS,JTE,KTS,KTE
	2098	!WRF
	2099	INTEGER, DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: LOWLYR
	2100	!
	2101	real, INTENT(IN) :: DELX,DELY
	2102	real,DIMENSION(*), INTENT(INOUT) :: MP_RESTART_STATE
	2103	real,DIMENSION(NX), INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
	2104	real,DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(OUT) :: &
	2105	& F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
	2106	INTEGER, PARAMETER :: ITLO=-60, ITHI=40
	2107	! integer,DIMENSION(ITLO:ITHI,4),INTENT(INOUT) :: NSTATS
	2108	! real,DIMENSION(ITLO:ITHI,5),INTENT(INOUT) :: QMAX
	2109	! real,DIMENSION(ITLO:ITHI,22),INTENT(INOUT) :: QTOT
	2110	! real,INTENT(INOUT) :: PRECtot(2),PRECmax(2)
	2111	real,INTENT(IN) :: DT,GSMDT
	2112	LOGICAL,INTENT(IN) :: allowed_to_read,restart
	2113	!
	2114	!-----------------------------------------------------------------------
	2115	! LOCAL VARIABLES
	2116	!-----------------------------------------------------------------------
	2117	REAL :: BBFR,DTPH,PI,DX,Thour_print
	2118	INTEGER :: I,IM,J,L,K,JTF,KTF,ITF
	2119	INTEGER :: etampnew_unit1
	2120	LOGICAL, PARAMETER :: PRINT_diag=.FALSE.
	2121	LOGICAL :: opened
	2122	LOGICAL , EXTERNAL :: wrf_dm_on_monitor
	2123	CHARACTER*80 errmess
	2124	!
	2125	!-----------------------------------------------------------------------
	2126	!
	2127	JTF=MIN0(JTE,JDE-1)
	2128	KTF=MIN0(KTE,KDE-1)
	2129	ITF=MIN0(ITE,IDE-1)
	2130	!
	2131	DO J=JTS,JTF
	2132	DO I=ITS,ITF
	2133	LOWLYR(I,J)=1
	2134	ENDDO
	2135	ENDDO
	2136	!
	2137	IF(.NOT.RESTART)THEN
	2138	DO J = jts,jte
	2139	DO K = kts,kte
	2140	DO I= its,ite
	2141	F_ICE_PHY(i,k,j)=0.
	2142	F_RAIN_PHY(i,k,j)=0.
	2143	F_RIMEF_PHY(i,k,j)=1.
	2144	ENDDO
	2145	ENDDO
	2146	ENDDO
	2147	ENDIF
	2148	!
	2149	!-----------------------------------------------------------------------
	2150	IF(ALLOWED_TO_READ)THEN
	2151	!-----------------------------------------------------------------------
	2152	!
	2153	DX=((DELX)2+(DELY)2)**.5/1000. ! Model resolution at equator (km)
	2154	DX=MIN(100., MAX(5., DX) )
	2155	!
	2156	!-- Relative humidity threshold for the onset of grid-scale condensation
	2157	!!-- 9/1/01: Assume the following functional dependence for 5 - 100 km resolution:
	2158	!! RHgrd=0.90 for dx=100 km, 0.98 for dx=5 km, where
	2159	! RHgrd=0.90+.08((100.-DX)/95.)*.5
	2160	!
	2161	DTPH=MAX(GSMDT*60.,DT)
	2162	DTPH=NINT(DTPH/DT)*DT
	2163	!
	2164	!--- Create lookup tables for saturation vapor pressure w/r/t water & ice
	2165	!
	2166	CALL GPVS
	2167	!
	2168	!--- Read in various lookup tables
	2169	!
	2170	IF ( wrf_dm_on_monitor() ) THEN
	2171	DO i = 31,99
	2172	INQUIRE ( i , OPENED = opened )
	2173	IF ( .NOT. opened ) THEN
	2174	etampnew_unit1 = i
	2175	GOTO 2061
	2176	ENDIF
	2177	ENDDO
	2178	etampnew_unit1 = -1
	2179	2061 CONTINUE
	2180	ENDIF
	2181	!
	2182	CALL wrf_dm_bcast_bytes ( etampnew_unit1 , IWORDSIZE )
	2183	!
	2184	IF ( etampnew_unit1 < 0 ) THEN
	2185	CALL wrf_error_fatal ( 'module_mp_etanew: etanewinit: Can not find '// &
	2186	'unused fortran unit to read in lookup table.' )
	2187	ENDIF
	2188	!
	2189	IF ( wrf_dm_on_monitor() ) THEN
	2190	!!was OPEN (UNIT=1,FILE="eta_micro_lookup.dat",FORM="UNFORMATTED")
	2191	OPEN(UNIT=etampnew_unit1,FILE="ETAMPNEW_DATA", &
	2192	& FORM="UNFORMATTED",STATUS="OLD",ERR=9061)
	2193	!
	2194	READ(etampnew_unit1) VENTR1
	2195	READ(etampnew_unit1) VENTR2
	2196	READ(etampnew_unit1) ACCRR
	2197	READ(etampnew_unit1) MASSR
	2198	READ(etampnew_unit1) VRAIN
	2199	READ(etampnew_unit1) RRATE
	2200	READ(etampnew_unit1) VENTI1
	2201	READ(etampnew_unit1) VENTI2
	2202	READ(etampnew_unit1) ACCRI
	2203	READ(etampnew_unit1) MASSI
	2204	READ(etampnew_unit1) VSNOWI
	2205	READ(etampnew_unit1) VEL_RF
	2206	! read(etampnew_unit1) my_growth ! Applicable only for DTPH=180 s
	2207	CLOSE (etampnew_unit1)
	2208	ENDIF
	2209	!
	2210	CALL wrf_dm_bcast_bytes ( VENTR1 , size ( VENTR1 ) * RWORDSIZE )
	2211	CALL wrf_dm_bcast_bytes ( VENTR2 , size ( VENTR2 ) * RWORDSIZE )
	2212	CALL wrf_dm_bcast_bytes ( ACCRR , size ( ACCRR ) * RWORDSIZE )
	2213	CALL wrf_dm_bcast_bytes ( MASSR , size ( MASSR ) * RWORDSIZE )
	2214	CALL wrf_dm_bcast_bytes ( VRAIN , size ( VRAIN ) * RWORDSIZE )
	2215	CALL wrf_dm_bcast_bytes ( RRATE , size ( RRATE ) * RWORDSIZE )
	2216	CALL wrf_dm_bcast_bytes ( VENTI1 , size ( VENTI1 ) * RWORDSIZE )
	2217	CALL wrf_dm_bcast_bytes ( VENTI2 , size ( VENTI2 ) * RWORDSIZE )
	2218	CALL wrf_dm_bcast_bytes ( ACCRI , size ( ACCRI ) * RWORDSIZE )
	2219	CALL wrf_dm_bcast_bytes ( MASSI , size ( MASSI ) * RWORDSIZE )
	2220	CALL wrf_dm_bcast_bytes ( VSNOWI , size ( VSNOWI ) * RWORDSIZE )
	2221	CALL wrf_dm_bcast_bytes ( VEL_RF , size ( VEL_RF ) * RWORDSIZE )
	2222	!
	2223	!--- Calculates coefficients for growth rates of ice nucleated in water
	2224	! saturated conditions, scaled by physics time step (lookup table)
	2225	!
	2226	CALL MY_GROWTH_RATES (DTPH)
	2227	! CALL MY_GROWTH_RATES (DTPH,MY_GROWTH)
	2228	!
	2229	PI=ACOS(-1.)
	2230	!
	2231	!--- Constants associated with Biggs (1953) freezing of rain, as parameterized
	2232	! following Lin et al. (JCAM, 1983) & Reisner et al. (1998, QJRMS).
	2233	!
	2234	ABFR=-0.66
	2235	BBFR=100.
	2236	CBFR=20.PIPIBBFRRHOL*1.E-21
	2237	!
	2238	!--- CIACW is used in calculating riming rates
	2239	! The assumed effective collection efficiency of cloud water rimed onto
	2240	! ice is =0.5 below:
	2241	!
	2242	CIACW=DTPH0.25PI0.5(1.E5)**C1
	2243	!
	2244	!--- CIACR is used in calculating freezing of rain colliding with large ice
	2245	! The assumed collection efficiency is 1.0
	2246	!
	2247	CIACR=PI*DTPH
	2248	!
	2249	!--- Based on rain lookup tables for mean diameters from 0.05 to 0.45 mm
	2250	! * Four different functional relationships of mean drop diameter as
	2251	! a function of rain rate (RR), derived based on simple fits to
	2252	! mass-weighted fall speeds of rain as functions of mean diameter
	2253	! from the lookup tables.
	2254	!
	2255	RR_DRmin=N0r0*RRATE(MDRmin) ! RR for mean drop diameter of .05 mm
	2256	RR_DR1=N0r0*RRATE(MDR1) ! RR for mean drop diameter of .10 mm
	2257	RR_DR2=N0r0*RRATE(MDR2) ! RR for mean drop diameter of .20 mm
	2258	RR_DR3=N0r0*RRATE(MDR3) ! RR for mean drop diameter of .32 mm
	2259	RR_DRmax=N0r0*RRATE(MDRmax) ! RR for mean drop diameter of .45 mm
	2260	!
	2261	RQR_DRmin=N0r0*MASSR(MDRmin) ! Rain content for mean drop diameter of .05 mm
	2262	RQR_DR1=N0r0*MASSR(MDR1) ! Rain content for mean drop diameter of .10 mm
	2263	RQR_DR2=N0r0*MASSR(MDR2) ! Rain content for mean drop diameter of .20 mm
	2264	RQR_DR3=N0r0*MASSR(MDR3) ! Rain content for mean drop diameter of .32 mm
	2265	RQR_DRmax=N0r0*MASSR(MDRmax) ! Rain content for mean drop diameter of .45 mm
	2266	C_N0r0=PIRHOLN0r0
	2267	CN0r0=1.E6/C_N0r0**.25
	2268	CN0r_DMRmin=1./(PIRHOLDMRmin**4)
	2269	CN0r_DMRmax=1./(PIRHOLDMRmax**4)
	2270	!
	2271	!--- CRACW is used in calculating collection of cloud water by rain (an
	2272	! assumed collection efficiency of 1.0)
	2273	!
	2274	CRACW=DTPH0.25PI*1.0
	2275	!
	2276	ESW0=1000.*FPVS0(T0C) ! Saturation vapor pressure at 0C
	2277	RFmax=1.1**Nrime ! Maximum rime factor allowed
	2278	!
	2279	!------------------------------------------------------------------------
	2280	!--------------- Constants passed through argument list -----------------
	2281	!------------------------------------------------------------------------
	2282	!
	2283	!--- Important parameters for self collection (autoconversion) of
	2284	! cloud water to rain.
	2285	!
	2286	!--- CRAUT is proportional to the rate that cloud water is converted by
	2287	! self collection to rain (autoconversion rate)
	2288	!
	2289	CRAUT=1.-EXP(-1.E-3*DTPH)
	2290	!
	2291	!--- QAUT0 is the threshold cloud content for autoconversion to rain
	2292	! needed for droplets to reach a diameter of 20 microns (following
	2293	! Manton and Cotton, 1977; Banta and Hanson, 1987, JCAM)
	2294	!--- QAUT0=1.2567, 0.8378, or 0.4189 g/m**3 for droplet number concentrations
	2295	! of 300, 200, and 100 cm**-3, respectively
	2296	!
	2297	QAUT0=PIRHOLNCW(20.E-6)*3/6.
	2298	!
	2299	!--- For calculating snow optical depths by considering bulk density of
	2300	! snow based on emails from Q. Fu (6/27-28/01), where optical
	2301	! depth (T) = 1.5SWP/(ReffDENS), SWP is snow water path, Reff
	2302	! is effective radius, and DENS is the bulk density of snow.
	2303	!
	2304	! SWP (kg/m*2)=(1.E-3 kg/g)SWPrad, SWPrad in g/m**2 used in radiation
	2305	! T = 1.51.E3SWPrad/(Reff*DENS)
	2306	!
	2307	! See derivation for MASSI(INDEXS), note equal to RHO*QSNOW/NSNOW
	2308	!
	2309	! SDENS=1.5e3/DENS, DENS=MASSI(INDEXS)/[PI(1.E-6INDEXS)**3]
	2310	!
	2311	DO I=MDImin,MDImax
	2312	SDENS(I)=PI1.5E-15FLOAT(III)/MASSI(I)
	2313	ENDDO
	2314	!
	2315	Thour_print=-DTPH/3600.
	2316
	2317	! SH 0211/2002
	2318	! IF (PRINT_diag) THEN
	2319
	2320	!-------- Total and maximum quantities
	2321	!
	2322	! NSTATS=0 ! Microphysical statistics dealing w/ grid-point counts
	2323	! QMAX=0. ! Microphysical statistics dealing w/ hydrometeor mass
	2324	! QTOT=0. ! Microphysical statistics dealing w/ hydrometeor mass
	2325	! PRECmax=0. ! Maximum precip rates (rain, snow) at surface (mm/h)
	2326	! PRECtot=0. ! Total precipitation (rain, snow) accumulation at surface
	2327	! ENDIF
	2328
	2329	!wrf
	2330	IF(.NOT.RESTART)THEN
	2331	MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
	2332	MP_RESTART_STATE(MY_T2+1)=C1XPVS0
	2333	MP_RESTART_STATE(MY_T2+2)=C2XPVS0
	2334	MP_RESTART_STATE(MY_T2+3)=C1XPVS
	2335	MP_RESTART_STATE(MY_T2+4)=C2XPVS
	2336	MP_RESTART_STATE(MY_T2+5)=CIACW
	2337	MP_RESTART_STATE(MY_T2+6)=CIACR
	2338	MP_RESTART_STATE(MY_T2+7)=CRACW
	2339	MP_RESTART_STATE(MY_T2+8)=CRAUT
	2340	TBPVS_STATE(1:NX) =TBPVS(1:NX)
	2341	TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
	2342	ENDIF
	2343
	2344	ENDIF ! Allowed_to_read
	2345
	2346	RETURN
	2347	!
	2348	!-----------------------------------------------------------------------
	2349	!
	2350	9061 CONTINUE
	2351	WRITE( errmess , '(A,I4)' ) &
	2352	'module_mp_etanew: error opening ETAMPNEW_DATA on unit ' &
	2353	&, etampnew_unit1
	2354	CALL wrf_error_fatal(errmess)
	2355	!
	2356	!-----------------------------------------------------------------------
	2357	END SUBROUTINE etanewinit
	2358	!
	2359	SUBROUTINE MY_GROWTH_RATES (DTPH)
	2360	! SUBROUTINE MY_GROWTH_RATES (DTPH,MY_GROWTH)
	2361	!
	2362	!--- Below are tabulated values for the predicted mass of ice crystals
	2363	! after 600 s of growth in water saturated conditions, based on
	2364	! calculations from Miller and Young (JAS, 1979). These values are
	2365	! crudely estimated from tabulated curves at 600 s from Fig. 6.9 of
	2366	! Young (1993). Values at temperatures colder than -27C were
	2367	! assumed to be invariant with temperature.
	2368	!
	2369	!--- Used to normalize Miller & Young (1979) calculations of ice growth
	2370	! over large time steps using their tabulated values at 600 s.
	2371	! Assumes 3D growth with time**1.5 following eq. (6.3) in Young (1993).
	2372	!
	2373	IMPLICIT NONE
	2374	!
	2375	REAL,INTENT(IN) :: DTPH
	2376	!
	2377	REAL DT_ICE
	2378	REAL,DIMENSION(35) :: MY_600
	2379	!WRF
	2380	!
	2381	!-----------------------------------------------------------------------
	2382	DATA MY_600 / &
	2383	& 5.5e-8, 1.4E-7, 2.8E-7, 6.E-7, 3.3E-6, &
	2384	& 2.E-6, 9.E-7, 8.8E-7, 8.2E-7, 9.4e-7, &
	2385	& 1.2E-6, 1.85E-6, 5.5E-6, 1.5E-5, 1.7E-5, &
	2386	& 1.5E-5, 1.E-5, 3.4E-6, 1.85E-6, 1.35E-6, &
	2387	& 1.05E-6, 1.E-6, 9.5E-7, 9.0E-7, 9.5E-7, &
	2388	& 9.5E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7, &
	2389	& 9.E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7 / ! -31 to -35 deg C
	2390	!
	2391	!-----------------------------------------------------------------------
	2392	!
	2393	DT_ICE=(DTPH/600.)**1.5
	2394	MY_GROWTH=DT_ICE*MY_600
	2395	!
	2396	!-----------------------------------------------------------------------
	2397	!
	2398	END SUBROUTINE MY_GROWTH_RATES
	2399	!
	2400	!-----------------------------------------------------------------------
	2401	!--------- Old GFS saturation vapor pressure lookup tables -----------
	2402	!-----------------------------------------------------------------------
	2403	!
	2404	SUBROUTINE GPVS
	2405	! ******************************************************************
	2406	!$$$ SUBPROGRAM DOCUMENTATION BLOCK
	2407	! . . .
	2408	! SUBPROGRAM: GPVS COMPUTE SATURATION VAPOR PRESSURE TABLE
	2409	! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
	2410	!
	2411	! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE TABLE AS A FUNCTION OF
	2412	! TEMPERATURE FOR THE TABLE LOOKUP FUNCTION FPVS.
	2413	! EXACT SATURATION VAPOR PRESSURES ARE CALCULATED IN SUBPROGRAM FPVSX.
	2414	! THE CURRENT IMPLEMENTATION COMPUTES A TABLE WITH A LENGTH
	2415	! OF 7501 FOR TEMPERATURES RANGING FROM 180.0 TO 330.0 KELVIN.
	2416	!
	2417	! PROGRAM HISTORY LOG:
	2418	! 91-05-07 IREDELL
	2419	! 94-12-30 IREDELL EXPAND TABLE
	2420	! 96-02-19 HONG ICE EFFECT
	2421	! 01-11-29 JIN MODIFIED FOR WRF
	2422	!
	2423	! USAGE: CALL GPVS
	2424	!
	2425	! SUBPROGRAMS CALLED:
	2426	! (FPVSX) - INLINABLE FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE
	2427	!
	2428	! COMMON BLOCKS:
	2429	! COMPVS - SCALING PARAMETERS AND TABLE FOR FUNCTION FPVS.
	2430	!
	2431	! ATTRIBUTES:
	2432	! LANGUAGE: FORTRAN 90
	2433	!
	2434	!$$$
	2435	IMPLICIT NONE
	2436	real :: X,XINC,T
	2437	integer :: JX
	2438	!----------------------------------------------------------------------
	2439	XINC=(XMAX-XMIN)/(NX-1)
	2440	C1XPVS=1.-XMIN/XINC
	2441	C2XPVS=1./XINC
	2442	C1XPVS0=1.-XMIN/XINC
	2443	C2XPVS0=1./XINC
	2444	!
	2445	DO JX=1,NX
	2446	X=XMIN+(JX-1)*XINC
	2447	T=X
	2448	TBPVS(JX)=FPVSX(T)
	2449	TBPVS0(JX)=FPVSX0(T)
	2450	ENDDO
	2451	!
	2452	END SUBROUTINE GPVS
	2453	!-----------------------------------------------------------------------
	2454	!***********************************************************************
	2455	!-----------------------------------------------------------------------
	2456	REAL FUNCTION FPVS(T)
	2457	!-----------------------------------------------------------------------
	2458	!$$$ SUBPROGRAM DOCUMENTATION BLOCK
	2459	! . . .
	2460	! SUBPROGRAM: FPVS COMPUTE SATURATION VAPOR PRESSURE
	2461	! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
	2462	!
	2463	! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE FROM THE TEMPERATURE.
	2464	! A LINEAR INTERPOLATION IS DONE BETWEEN VALUES IN A LOOKUP TABLE
	2465	! COMPUTED IN GPVS. SEE DOCUMENTATION FOR FPVSX FOR DETAILS.
	2466	! INPUT VALUES OUTSIDE TABLE RANGE ARE RESET TO TABLE EXTREMA.
	2467	! THE INTERPOLATION ACCURACY IS ALMOST 6 DECIMAL PLACES.
	2468	! ON THE CRAY, FPVS IS ABOUT 4 TIMES FASTER THAN EXACT CALCULATION.
	2469	! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
	2470	!
	2471	! PROGRAM HISTORY LOG:
	2472	! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
	2473	! 94-12-30 IREDELL EXPAND TABLE
	2474	! 96-02-19 HONG ICE EFFECT
	2475	! 01-11-29 JIN MODIFIED FOR WRF
	2476	!
	2477	! USAGE: PVS=FPVS(T)
	2478	!
	2479	! INPUT ARGUMENT LIST:
	2480	! T - REAL TEMPERATURE IN KELVIN
	2481	!
	2482	! OUTPUT ARGUMENT LIST:
	2483	! FPVS - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
	2484	!
	2485	! ATTRIBUTES:
	2486	! LANGUAGE: FORTRAN 90
	2487	!$$$
	2488	IMPLICIT NONE
	2489	real,INTENT(IN) :: T
	2490	real XJ
	2491	integer :: JX
	2492	!-----------------------------------------------------------------------
	2493	XJ=MIN(MAX(C1XPVS+C2XPVS*T,1.),FLOAT(NX))
	2494	JX=MIN(XJ,NX-1.)
	2495	FPVS=TBPVS(JX)+(XJ-JX)*(TBPVS(JX+1)-TBPVS(JX))
	2496	!
	2497	END FUNCTION FPVS
	2498	!-----------------------------------------------------------------------
	2499	!-----------------------------------------------------------------------
	2500	REAL FUNCTION FPVS0(T)
	2501	!-----------------------------------------------------------------------
	2502	IMPLICIT NONE
	2503	real,INTENT(IN) :: T
	2504	real :: XJ1
	2505	integer :: JX1
	2506	!-----------------------------------------------------------------------
	2507	XJ1=MIN(MAX(C1XPVS0+C2XPVS0*T,1.),FLOAT(NX))
	2508	JX1=MIN(XJ1,NX-1.)
	2509	FPVS0=TBPVS0(JX1)+(XJ1-JX1)*(TBPVS0(JX1+1)-TBPVS0(JX1))
	2510	!
	2511	END FUNCTION FPVS0
	2512	!-----------------------------------------------------------------------
	2513	!***********************************************************************
	2514	!-----------------------------------------------------------------------
	2515	REAL FUNCTION FPVSX(T)
	2516	!-----------------------------------------------------------------------
	2517	!$$$ SUBPROGRAM DOCUMENTATION BLOCK
	2518	! . . .
	2519	! SUBPROGRAM: FPVSX COMPUTE SATURATION VAPOR PRESSURE
	2520	! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
	2521	!
	2522	! ABSTRACT: EXACTLY COMPUTE SATURATION VAPOR PRESSURE FROM TEMPERATURE.
	2523	! THE WATER MODEL ASSUMES A PERFECT GAS, CONSTANT SPECIFIC HEATS
	2524	! FOR GAS AND LIQUID, AND NEGLECTS THE VOLUME OF THE LIQUID.
	2525	! THE MODEL DOES ACCOUNT FOR THE VARIATION OF THE LATENT HEAT
	2526	! OF CONDENSATION WITH TEMPERATURE. THE ICE OPTION IS NOT INCLUDED.
	2527	! THE CLAUSIUS-CLAPEYRON EQUATION IS INTEGRATED FROM THE TRIPLE POINT
	2528	! TO GET THE FORMULA
	2529	! PVS=PSATK(TRXA)EXP(XB*(1.-TR))
	2530	! WHERE TR IS TTP/T AND OTHER VALUES ARE PHYSICAL CONSTANTS
	2531	! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
	2532	!
	2533	! PROGRAM HISTORY LOG:
	2534	! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
	2535	! 94-12-30 IREDELL EXACT COMPUTATION
	2536	! 96-02-19 HONG ICE EFFECT
	2537	! 01-11-29 JIN MODIFIED FOR WRF
	2538	!
	2539	! USAGE: PVS=FPVSX(T)
	2540	! REFERENCE: EMANUEL(1994),116-117
	2541	!
	2542	! INPUT ARGUMENT LIST:
	2543	! T - REAL TEMPERATURE IN KELVIN
	2544	!
	2545	! OUTPUT ARGUMENT LIST:
	2546	! FPVSX - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
	2547	!
	2548	! ATTRIBUTES:
	2549	! LANGUAGE: FORTRAN 90
	2550	!$$$
	2551	IMPLICIT NONE
	2552	!-----------------------------------------------------------------------
	2553	real, parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
	2554	, CLIQ=4.1855E+3,CVAP= 1.8460E+3 &
	2555	, CICE=2.1060E+3,HSUB=2.8340E+6
	2556	!
	2557	real, parameter :: PSATK=PSAT*1.E-3
	2558	real, parameter :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
	2559	real, parameter :: DLDTI=CVAP-CICE &
	2560	, XAI=-DLDTI/RV,XBI=XAI+HSUB/(RV*TTP)
	2561	real T,TR
	2562	!-----------------------------------------------------------------------
	2563	TR=TTP/T
	2564	!
	2565	IF(T.GE.TTP)THEN
	2566	FPVSX=PSATK(TRXA)EXP(XB*(1.-TR))
	2567	ELSE
	2568	FPVSX=PSATK(TRXAI)EXP(XBI*(1.-TR))
	2569	ENDIF
	2570	!
	2571	END FUNCTION FPVSX
	2572	!-----------------------------------------------------------------------
	2573	!-----------------------------------------------------------------------
	2574	REAL FUNCTION FPVSX0(T)
	2575	!-----------------------------------------------------------------------
	2576	IMPLICIT NONE
	2577	real,parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
	2578	, CLIQ=4.1855E+3,CVAP=1.8460E+3 &
	2579	, CICE=2.1060E+3,HSUB=2.8340E+6
	2580	real,PARAMETER :: PSATK=PSAT*1.E-3
	2581	real,PARAMETER :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
	2582	real,PARAMETER :: DLDTI=CVAP-CICE &
	2583	, XAI=-DLDT/RV,XBI=XA+HSUB/(RV*TTP)
	2584	real :: T,TR
	2585	!-----------------------------------------------------------------------
	2586	TR=TTP/T
	2587	FPVSX0=PSATK(TRXA)EXP(XB*(1.-TR))
	2588	!
	2589	END FUNCTION FPVSX0
	2590	!
	2591	END MODULE module_mp_etanew

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

Download in other formats:

Original Format