Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_etanew.F @ 2295

Last change on this file since 2295 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 108.6 KB

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

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

Download in other formats:

Original Format