Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_etanew.F @ 3094

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

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

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

Download in other formats:

Original Format