Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_etanew.F @ 3547

Last change on this file since 3547 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 107.4 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),PRIVATE,SAVE :: TBPVS,TBPVS0
33	REAL, PRIVATE,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=-30. &
69	& ,T_ICEK=T0C+T_ICE &
70	& ,T_ICE_init=-5. &
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	!
1097	!--- Reduce excessive accumulation of ice at upper levels
1098	! associated with strong grid-resolved ascent
1099	!
1100	!--- Force NLICE to be between NLImin and NLImax
1101	!
1102	DUM=MAX(NLImin, MIN(NLImax, NLICE) )
1103	XLI=RHO*(QTICE/DUM-XSIMASS)/RimeF1
1104	IF (XLI .LE. MASSI(MDImin) ) THEN
1105	INDEXS=MDImin
1106	ELSE IF (XLI .LE. MASSI(450) ) THEN
1107	DLI=9.5885E5XLI*.42066 ! DLI in microns
1108	INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
1109	ELSE IF (XLI .LE. MASSI(MDImax) ) THEN
1110	DLI=3.9751E6XLI*.49870 ! DLI in microns
1111	INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
1112	ELSE
1113	INDEXS=MDImax
1114	!
1115	!--- 8/22/01: Increase density of large ice if maximum limits
1116	! are reached for number concentration (NLImax) and mean size
1117	! (MDImax). Done to increase fall out of ice.
1118	!
1119	IF (DUM .GE. NLImax) &
1120	& RimeF1=RHO*(QTICE/NLImax-XSIMASS)/MASSI(INDEXS)
1121	ENDIF ! End IF (XLI .LE. MASSI(MDImin) )
1122	! WRITE(6,"(4(a12,g11.4,1x))")
1123	! & '{$ TC=',TC,'P=',.01*PP,'NLICE=',NLICE,'DUM=',DUM,
1124	! & '{$ XLI=',XLI,'INDEXS=',FLOAT(INDEXS),'RHO=',RHO,'QTICE=',QTICE,
1125	! & '{$ XSIMASS=',XSIMASS,'RimeF1=',RimeF1
1126	ENDIF ! End IF ( (NLICE.GE.NLImin .AND. NLICE.LE.NLImax) ...
1127	ENDDO ! End DO IPASS=0,1
1128	ENDIF ! End IF (QI.LE.EPSQ .AND. ASNOW.LE.CLIMIT)
1129	ENDIF ! End IF (ICE_logical)
1130	!
1131	!----------------------------------------------------------------------
1132	!--------------- Calculate individual processes -----------------------
1133	!----------------------------------------------------------------------
1134	!
1135	!--- Cloud water autoconversion to rain and collection by rain
1136	!
1137	IF (QW.GT.EPSQ .AND. TC.GE.T_ICE) THEN
1138	!
1139	!--- QW0 could be modified based on land/sea properties,
1140	! presence of convection, etc. This is why QAUT0 and CRAUT
1141	! are passed into the subroutine as externally determined
1142	! parameters. Can be changed in the future if desired.
1143	!
1144	QW0=QAUT0*RRHO
1145	PRAUT=MAX(0., QW-QW0)*CRAUT
1146	IF (QLICE .GT. EPSQ) THEN
1147	!
1148	!--- Collection of cloud water by large ice particles ("snow")
1149	! PIACWI=PIACW for riming, PIACWI=0 for shedding
1150	!
1151	FWS=MIN(1., CIACWVEL_INCNLICEACCRI(INDEXS)/PP*C1)
1152	PIACW=FWS*QW
1153	IF (TC .LT. 0.) PIACWI=PIACW ! Large ice riming
1154	ENDIF ! End IF (QLICE .GT. EPSQ)
1155	ENDIF ! End IF (QW.GT.EPSQ .AND. TC.GE.T_ICE)
1156	!
1157	!----------------------------------------------------------------------
1158	!--- Loop around some of the ice-phase processes if no ice should be present
1159	!----------------------------------------------------------------------
1160	!
1161	IF (ICE_logical .EQV. .FALSE.) GO TO 20
1162	!
1163	!--- Now the pretzel logic of calculating ice deposition
1164	!
1165	IF (TC.LT.T_ICE .AND. (WV.GT.QSIgrd .OR. QW.GT.EPSQ)) THEN
1166	!
1167	!--- Adjust to ice saturation at T<T_ICE (-10C) if supersaturated.
1168	! Sources of ice due to nucleation and convective detrainment are
1169	! either poorly understood, poorly resolved at typical NWP
1170	! resolutions, or are not represented (e.g., no detrained
1171	! condensate in BMJ Cu scheme).
1172	!
1173	PCOND=-QW
1174	DUM1=TK+XLV1*PCOND ! Updated (dummy) temperature (deg K)
1175	DUM2=WV+QW ! Updated (dummy) water vapor mixing ratio
1176	DUM=1000.*FPVS(DUM1) ! Updated (dummy) saturation vapor pressure w/r/t ice
1177	DUM=RHgrdEPSDUM/(PP-DUM) ! Updated (dummy) saturation mixing ratio w/r/t ice
1178	IF (DUM2 .GT. DUM) PIDEP=DEPOSIT (PP, DUM1, DUM2)
1179	DWVi=0. ! Used only for debugging
1180	!
1181	ELSE IF (TC .LT. 0.) THEN
1182	!
1183	!--- These quantities are handy for ice deposition/sublimation
1184	! PIDEP_max - max deposition or minimum sublimation to ice saturation
1185	!
1186	DENOMI=1.+XLS2QSITK2
1187	DWVi=MIN(WVQW,QSW)-QSI
1188	PIDEP_max=MAX(PILOSS, DWVi/DENOMI)
1189	IF (QTICE .GT. 0.) THEN
1190	!
1191	!--- Calculate ice deposition/sublimation
1192	! * SFACTOR - [VEL_INC*.5][Schmidt*(1./3.)][(RHO/DYNVIS)**.5],
1193	! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
1194	! * Units: SFACTOR - s.5/m ; ABI - m2/s ; NLICE - m**-3 ;
1195	! VENTIL, VENTIS - m**-2 ; VENTI1 - m ;
1196	! VENTI2 - m2/s.5 ; DIDEP - unitless
1197	!
1198	SFACTOR=VEL_INC*.5(RHO/(DIFFUSDIFFUSDYNVIS))**C2
1199	ABI=1./(RHOXLS3QSI*TK2/THERM_COND+1./DIFFUS)
1200	!
1201	!--- VENTIL - Number concentration * ventilation factors for large ice
1202	!--- VENTIS - Number concentration * ventilation factors for small ice
1203	!
1204	!--- Variation in the number concentration of ice with time is not
1205	! accounted for in these calculations (could be in the future).
1206	!
1207	VENTIL=(VENTI1(INDEXS)+SFACTORVENTI2(INDEXS))NLICE
1208	VENTIS=(VENTI1(MDImin)+SFACTORVENTI2(MDImin))NSmICE
1209	DIDEP=ABI(VENTIL+VENTIS)DTPH
1210	!
1211	!--- Account for change in water vapor supply w/ time
1212	!
1213	IF (DIDEP .GE. Xratio)then
1214	DIDEP=(1.-EXP(-DIDEP*DENOMI))/DENOMI
1215	endif
1216	IF (DWVi .GT. 0.) THEN
1217	PIDEP=MIN(DWVi*DIDEP, PIDEP_max)
1218	ELSE IF (DWVi .LT. 0.) THEN
1219	PIDEP=MAX(DWVi*DIDEP, PIDEP_max)
1220	ENDIF
1221	!
1222	ELSE IF (WVQW.GT.QSI .AND. TC.LE.T_ICE_init) THEN
1223	!
1224	!--- Ice nucleation in near water-saturated conditions. Ice crystal
1225	! growth during time step calculated using Miller & Young (1979, JAS).
1226	!--- These deposition rates could drive conditions below water saturation,
1227	! which is the basis of these calculations. Intended to approximate
1228	! more complex & computationally intensive calculations.
1229	!
1230	INDEX_MY=MAX(MY_T1, MIN( INT(.5-TC), MY_T2 ) )
1231	!
1232	!--- DUM1 is the supersaturation w/r/t ice at water-saturated conditions
1233	!
1234	!--- DUM2 is the number of ice crystals nucleated at water-saturated
1235	! conditions based on Meyers et al. (JAM, 1992).
1236	!
1237	!--- Prevent unrealistically large ice initiation (limited by PIDEP_max)
1238	! if DUM2 values are increased in future experiments
1239	!
1240	DUM1=QSW/QSI-1.
1241	DUM2=1.E3EXP(12.96DUM1-.639)
1242	PIDEP=MIN(PIDEP_max, DUM2MY_GROWTH(INDEX_MY)RRHO)
1243	!
1244	ENDIF ! End IF (QTICE .GT. 0.)
1245	!
1246	ENDIF ! End IF (TC.LT.T_ICE .AND. (WV.GT.QSIgrd .OR. QW.GT.EPSQ))
1247	!
1248	!------------------------------------------------------------------------
1249	!
1250	20 CONTINUE ! Jump here if conditions for ice are not present
1251
1252
1253	!
1254	!------------------------------------------------------------------------
1255	!
1256	!--- Cloud water condensation
1257	!
1258	IF (TC.GE.T_ICE .AND. (QW.GT.EPSQ .OR. WV.GT.QSWgrd)) THEN
1259	IF (PIACWI.EQ.0. .AND. PIDEP.EQ.0.) THEN
1260	PCOND=CONDENSE (PP, QW, TK, WV)
1261	ELSE
1262	!
1263	!--- Modify cloud condensation in response to ice processes
1264	!
1265	DUM=XLVQSWgrdRCPRV*TK2
1266	DENOMWI=1.+XLS*DUM
1267	DENOMF=XLF*DUM
1268	DUM=MAX(0., PIDEP)
1269	PCOND=(WV-QSWgrd-DENOMWIDUM-DENOMFPIACWI)/DENOMW
1270	DUM1=-QW
1271	DUM2=PCOND-PIACW
1272	IF (DUM2 .LT. DUM1) THEN
1273	!
1274	!--- Limit cloud water sinks
1275	!
1276	DUM=DUM1/DUM2
1277	PCOND=DUM*PCOND
1278	PIACW=DUM*PIACW
1279	PIACWI=DUM*PIACWI
1280	ENDIF ! End IF (DUM2 .LT. DUM1)
1281	ENDIF ! End IF (PIACWI.EQ.0. .AND. PIDEP.EQ.0.)
1282	ENDIF ! End IF (TC.GE.T_ICE .AND. (QW.GT.EPSQ .OR. WV.GT.QSWgrd))
1283	!
1284	!--- Limit freezing of accreted rime to prevent temperature oscillations,
1285	! a crude Schumann-Ludlam limit (p. 209 of Young, 1993).
1286	!
1287	TCC=TC+XLV1PCOND+XLS1PIDEP+XLF1*PIACWI
1288	IF (TCC .GT. 0.) THEN
1289	PIACWI=0.
1290	TCC=TC+XLV1PCOND+XLS1PIDEP
1291	ENDIF
1292	IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical) THEN
1293	!
1294	!--- Calculate melting and evaporation/condensation
1295	! * Units: SFACTOR - s.5/m ; ABI - m2/s ; NLICE - m**-3 ;
1296	! VENTIL - m**-2 ; VENTI1 - m ;
1297	! VENTI2 - m2/s.5 ; CIEVP - /s
1298	!
1299	SFACTOR=VEL_INC*.5(RHO/(DIFFUSDIFFUSDYNVIS))**C2
1300	VENTIL=NLICE(VENTI1(INDEXS)+SFACTORVENTI2(INDEXS))
1301	AIEVP=VENTILDIFFUSDTPH
1302	IF (AIEVP .LT. Xratio) THEN
1303	DIEVP=AIEVP
1304	ELSE
1305	DIEVP=1.-EXP(-AIEVP)
1306	ENDIF
1307	QSW0=EPS*ESW0/(PP-ESW0)
1308	DWV0=MIN(WV,QSW)-QSW0
1309	DUM=QW+PCOND
1310	IF (WV.LT.QSW .AND. DUM.LE.EPSQ) THEN
1311	!
1312	!--- Evaporation from melting snow (sink of snow) or shedding
1313	! of water condensed onto melting snow (source of rain)
1314	!
1315	DUM=DWV0*DIEVP
1316	PIEVP=MAX( MIN(0., DUM), PILOSS)
1317	PICND=MAX(0., DUM)
1318	ENDIF ! End IF (WV.LT.QSW .AND. DUM.LE.EPSQ)
1319	PIMLT=THERM_CONDTCCVENTILRRHODTPH/XLF
1320	!
1321	!--- Limit melting to prevent temperature oscillations across 0C
1322	!
1323	DUM1=MAX( 0., (TCC+XLV1*PIEVP)/XLF1 )
1324	PIMLT=MIN(PIMLT, DUM1)
1325	!
1326	!--- Limit loss of snow by melting (>0) and evaporation
1327	!
1328	DUM=PIEVP-PIMLT
1329	IF (DUM .LT. PILOSS) THEN
1330	DUM1=PILOSS/DUM
1331	PIMLT=PIMLT*DUM1
1332	PIEVP=PIEVP*DUM1
1333	ENDIF ! End IF (DUM .GT. QTICE)
1334	ENDIF ! End IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical)
1335	!
1336	!--- IMPORTANT: Estimate time-averaged properties.
1337	!
1338	! * TOT_RAIN - total mass of rain before microphysics, which is the sum of
1339	! the total mass of rain in the current layer and the input
1340	! flux of rain from above
1341	! * VRAIN1 - fall speed of rain into grid from above (with air resistance correction)
1342	! * QTRAIN - time-averaged mixing ratio of rain (kg/kg)
1343	! * PRLOSS - greatest loss (<0) of rain, removing all rain falling from
1344	! above and the rain within the layer
1345	! * RQR - rain content (kg/m**3)
1346	! * INDEXR - mean size of rain drops to the nearest 1 micron in size
1347	! * N0r - intercept of rain size distribution (typically 106 m-4)
1348	!
1349	TOT_RAIN=0.
1350	VRAIN1=0.
1351	QTRAIN=0.
1352	PRLOSS=0.
1353	RQR=0.
1354	N0r=0.
1355	INDEXR1=0
1356	INDEXR=MDRmin
1357	IF (RAIN_logical) THEN
1358	IF (ARAIN .LE. 0.) THEN
1359	INDEXR=MDRmin
1360	VRAIN1=0.
1361	ELSE
1362	!
1363	!--- INDEXR (related to mean diameter) & N0r could be modified
1364	! by land/sea properties, presence of convection, etc.
1365	!
1366	!--- Rain rate normalized to a density of 1.194 kg/m**3
1367	!
1368	RR=ARAIN/(DTPH*GAMMAR)
1369	!
1370	IF (RR .LE. RR_DRmin) THEN
1371	!
1372	!--- Assume fixed mean diameter of rain (0.2 mm) for low rain rates,
1373	! instead vary N0r with rain rate
1374	!
1375	INDEXR=MDRmin
1376	ELSE IF (RR .LE. RR_DR1) THEN
1377	!
1378	!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1379	! for mean diameters (Dr) between 0.05 and 0.10 mm:
1380	! V(Dr)=5.6023e4Dr*1.136, V in m/s and Dr in m
1381	! RR = PI1000.N0r05.6023e4Dr*(4+1.136) = 1.408e15Dr**5.136,
1382	! RR in kg/(m*2s)
1383	! Dr (m) = 1.123e-3RR.1947 -> Dr (microns) = 1.123e3RR**.1947
1384	!
1385	INDEXR=INT( 1.123E3RR*.1947 + .5 )
1386	INDEXR=MAX( MDRmin, MIN(INDEXR, MDR1) )
1387	ELSE IF (RR .LE. RR_DR2) THEN
1388	!
1389	!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1390	! for mean diameters (Dr) between 0.10 and 0.20 mm:
1391	! V(Dr)=1.0867e4Dr*.958, V in m/s and Dr in m
1392	! RR = PI1000.N0r01.0867e4Dr*(4+.958) = 2.731e14Dr**4.958,
1393	! RR in kg/(m*2s)
1394	! Dr (m) = 1.225e-3RR.2017 -> Dr (microns) = 1.225e3RR**.2017
1395	!
1396	INDEXR=INT( 1.225E3RR*.2017 + .5 )
1397	INDEXR=MAX( MDR1, MIN(INDEXR, MDR2) )
1398	ELSE IF (RR .LE. RR_DR3) THEN
1399	!
1400	!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1401	! for mean diameters (Dr) between 0.20 and 0.32 mm:
1402	! V(Dr)=2831.Dr*.80, V in m/s and Dr in m
1403	! RR = PI1000.N0r02831.Dr*(4+.80) = 7.115e13Dr**4.80,
1404	! RR in kg/(m*2s)
1405	! Dr (m) = 1.3006e-3RR.2083 -> Dr (microns) = 1.3006e3RR**.2083
1406	!
1407	INDEXR=INT( 1.3006E3RR*.2083 + .5 )
1408	INDEXR=MAX( MDR2, MIN(INDEXR, MDR3) )
1409	ELSE IF (RR .LE. RR_DRmax) THEN
1410	!
1411	!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1412	! for mean diameters (Dr) between 0.32 and 0.45 mm:
1413	! V(Dr)=944.8Dr*.6636, V in m/s and Dr in m
1414	! RR = PI1000.N0r0944.8Dr*(4+.6636) = 2.3745e13Dr**4.6636,
1415	! RR in kg/(m*2s)
1416	! Dr (m) = 1.355e-3RR.2144 -> Dr (microns) = 1.355e3RR**.2144
1417	!
1418	INDEXR=INT( 1.355E3RR*.2144 + .5 )
1419	INDEXR=MAX( MDR3, MIN(INDEXR, MDRmax) )
1420	ELSE
1421	!
1422	!--- Assume fixed mean diameter of rain (0.45 mm) for high rain rates,
1423	! instead vary N0r with rain rate
1424	!
1425	INDEXR=MDRmax
1426	ENDIF ! End IF (RR .LE. RR_DRmin) etc.
1427	VRAIN1=GAMMAR*VRAIN(INDEXR)
1428	ENDIF ! End IF (ARAIN .LE. 0.)
1429	INDEXR1=INDEXR ! For debugging only
1430	TOT_RAIN=THICKQR+BLENDARAIN
1431	QTRAIN=TOT_RAIN/(THICK+BLDTRH*VRAIN1)
1432	PRLOSS=-TOT_RAIN/THICK
1433	RQR=RHO*QTRAIN
1434	!
1435	!--- RQR - time-averaged rain content (kg/m**3)
1436	!
1437	IF (RQR .LE. RQR_DRmin) THEN
1438	N0r=MAX(N0rmin, CN0r_DMRmin*RQR)
1439	INDEXR=MDRmin
1440	ELSE IF (RQR .GE. RQR_DRmax) THEN
1441	N0r=CN0r_DMRmax*RQR
1442	INDEXR=MDRmax
1443	ELSE
1444	N0r=N0r0
1445	INDEXR=MAX( XMRmin, MIN(CN0r0RQR*.25, XMRmax) )
1446	ENDIF
1447	!
1448	IF (TC .LT. T_ICE) THEN
1449	PIACR=-PRLOSS
1450	ELSE
1451	DWVr=WV-PCOND-QSW
1452	DUM=QW+PCOND
1453	IF (DWVr.LT.0. .AND. DUM.LE.EPSQ) THEN
1454	!
1455	!--- Rain evaporation
1456	!
1457	! * RFACTOR - [GAMMAR*.5][Schmidt*(1./3.)][(RHO/DYNVIS)**.5],
1458	! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
1459	!
1460	! * Units: RFACTOR - s.5/m ; ABW - m2/s ; VENTR - m**-2 ;
1461	! N0r - m-4 ; VENTR1 - m2 ; VENTR2 - m3/s.5 ;
1462	! CREVP - unitless
1463	!
1464	RFACTOR=GAMMAR*.5(RHO/(DIFFUSDIFFUSDYNVIS))**C2
1465	ABW=1./(RHO*XLV2/THERM_COND+1./DIFFUS)
1466	!
1467	!--- Note that VENTR1, VENTR2 lookup tables do not include the
1468	! 1/Davg multiplier as in the ice tables
1469	!
1470	VENTR=N0r(VENTR1(INDEXR)+RFACTORVENTR2(INDEXR))
1471	CREVP=ABWVENTRDTPH
1472	IF (CREVP .LT. Xratio) THEN
1473	DUM=DWVr*CREVP
1474	ELSE
1475	DUM=DWVr(1.-EXP(-CREVPDENOMW))/DENOMW
1476	ENDIF
1477	PREVP=MAX(DUM, PRLOSS)
1478	ELSE IF (QW .GT. EPSQ) THEN
1479	FWR=CRACWGAMMARN0r*ACCRR(INDEXR)
1480	PRACW=MIN(1.,FWR)*QW
1481	ENDIF ! End IF (DWVr.LT.0. .AND. DUM.LE.EPSQ)
1482	!
1483	IF (TC.LT.0. .AND. TCC.LT.0.) THEN
1484	!
1485	!--- Biggs (1953) heteorogeneous freezing (e.g., Lin et al., 1983)
1486	! - Rescaled mean drop diameter from microns (INDEXR) to mm (DUM) to prevent underflow
1487	!
1488	DUM=.001*FLOAT(INDEXR)
1489	DUM=(EXP(ABFRTC)-1.)DUMDUMDUMDUMDUMDUMDUM
1490	PIACR=MIN(CBFRN0rRRHO*DUM, QTRAIN)
1491	IF (QLICE .GT. EPSQ) THEN
1492	!
1493	!--- Freezing of rain by collisions w/ large ice
1494	!
1495	DUM=GAMMAR*VRAIN(INDEXR)
1496	DUM1=DUM-VSNOW
1497	!
1498	!--- DUM2 - Difference in spectral fall speeds of rain and
1499	! large ice, parameterized following eq. (48) on p. 112 of
1500	! Murakami (J. Meteor. Soc. Japan, 1990)
1501	!
1502	DUM2=(DUM1DUM1+.04DUMVSNOW)*.5
1503	DUM1=5.E-12INDEXRINDEXR+2.E-12INDEXRINDEXS &
1504	& +.5E-12INDEXSINDEXS
1505	FIR=MIN(1., CIACRNLICEDUM1*DUM2)
1506	!
1507	!--- Future? Should COLLECTION BY SMALL ICE SHOULD BE INCLUDED???
1508	!
1509	PIACR=MIN(PIACR+FIR*QTRAIN, QTRAIN)
1510	ENDIF ! End IF (QLICE .GT. EPSQ)
1511	DUM=PREVP-PIACR
1512	If (DUM .LT. PRLOSS) THEN
1513	DUM1=PRLOSS/DUM
1514	PREVP=DUM1*PREVP
1515	PIACR=DUM1*PIACR
1516	ENDIF ! End If (DUM .LT. PRLOSS)
1517	ENDIF ! End IF (TC.LT.0. .AND. TCC.LT.0.)
1518	ENDIF ! End IF (TC .LT. T_ICE)
1519	ENDIF ! End IF (RAIN_logical)
1520	!
1521	!----------------------------------------------------------------------
1522	!---------------------- Main Budget Equations -------------------------
1523	!----------------------------------------------------------------------
1524	!
1525	!
1526	!-----------------------------------------------------------------------
1527	!--- Update fields, determine characteristics for next lower layer ----
1528	!-----------------------------------------------------------------------
1529	!
1530	!--- Carefully limit sinks of cloud water
1531	!
1532	DUM1=PIACW+PRAUT+PRACW-MIN(0.,PCOND)
1533	IF (DUM1 .GT. QW) THEN
1534	DUM=QW/DUM1
1535	PIACW=DUM*PIACW
1536	PIACWI=DUM*PIACWI
1537	PRAUT=DUM*PRAUT
1538	PRACW=DUM*PRACW
1539	IF (PCOND .LT. 0.) PCOND=DUM*PCOND
1540	ENDIF
1541	PIACWR=PIACW-PIACWI ! TC >= 0C
1542	!
1543	!--- QWnew - updated cloud water mixing ratio
1544	!
1545	DELW=PCOND-PIACW-PRAUT-PRACW
1546	QWnew=QW+DELW
1547	IF (QWnew .LE. EPSQ) QWnew=0.
1548	IF (QW.GT.0. .AND. QWnew.NE.0.) THEN
1549	DUM=QWnew/QW
1550	IF (DUM .LT. TOLER) QWnew=0.
1551	ENDIF
1552	!
1553	!--- Update temperature and water vapor mixing ratios
1554	!
1555	DELT= XLV1*(PCOND+PIEVP+PICND+PREVP) &
1556	& +XLS1PIDEP+XLF1(PIACWI+PIACR-PIMLT)
1557	Tnew=TK+DELT
1558	!
1559	DELV=-PCOND-PIDEP-PIEVP-PICND-PREVP
1560	WVnew=WV+DELV
1561	!
1562	!--- Update ice mixing ratios
1563	!
1564	!---
1565	! * TOT_ICEnew - total mass (small & large) ice after microphysics,
1566	! which is the sum of the total mass of large ice in the
1567	! current layer and the flux of ice out of the grid box below
1568	! * RimeF - Rime Factor, which is the mass ratio of total (unrimed &
1569	! rimed) ice mass to the unrimed ice mass (>=1)
1570	! * QInew - updated mixing ratio of total (large & small) ice in layer
1571	! -> TOT_ICEnew=QInewTHICK+BLDTRHQLICEnew*VSNOW
1572	! -> But QLICEnew=QInew*FLIMASS, so
1573	! -> TOT_ICEnew=QInew(THICK+BLDTRHFLIMASS*VSNOW)
1574	! * ASNOWnew - updated accumulation of snow at bottom of grid cell
1575	!---
1576	!
1577	DELI=0.
1578	RimeF=1.
1579	IF (ICE_logical) THEN
1580	DELI=PIDEP+PIEVP+PIACWI+PIACR-PIMLT
1581	TOT_ICEnew=TOT_ICE+THICK*DELI
1582	IF (TOT_ICE.GT.0. .AND. TOT_ICEnew.NE.0.) THEN
1583	DUM=TOT_ICEnew/TOT_ICE
1584	IF (DUM .LT. TOLER) TOT_ICEnew=0.
1585	ENDIF
1586	IF (TOT_ICEnew .LE. CLIMIT) THEN
1587	TOT_ICEnew=0.
1588	RimeF=1.
1589	QInew=0.
1590	ASNOWnew=0.
1591	ELSE
1592	!
1593	!--- Update rime factor if appropriate
1594	!
1595	DUM=PIACWI+PIACR
1596	IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ) THEN
1597	RimeF=RimeF1
1598	ELSE
1599	!
1600	!--- Rime Factor, RimeF = (Total ice mass)/(Total unrimed ice mass)
1601	! DUM1 - Total ice mass, rimed & unrimed
1602	! DUM2 - Estimated mass of unrimed ice
1603	!
1604	DUM1=TOT_ICE+THICK*(PIDEP+DUM)
1605	DUM2=TOT_ICE/RimeF1+THICK*PIDEP
1606	IF (DUM2 .LE. 0.) THEN
1607	RimeF=RFmax
1608	ELSE
1609	RimeF=MIN(RFmax, MAX(1., DUM1/DUM2) )
1610	ENDIF
1611	ENDIF ! End IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ)
1612	QInew=TOT_ICEnew/(THICK+BLDTRHFLIMASSVSNOW)
1613	IF (QInew .LE. EPSQ) QInew=0.
1614	IF (QI.GT.0. .AND. QInew.NE.0.) THEN
1615	DUM=QInew/QI
1616	IF (DUM .LT. TOLER) QInew=0.
1617	ENDIF
1618	ASNOWnew=BLDTRHFLIMASSVSNOW*QInew
1619	IF (ASNOW.GT.0. .AND. ASNOWnew.NE.0.) THEN
1620	DUM=ASNOWnew/ASNOW
1621	IF (DUM .LT. TOLER) ASNOWnew=0.
1622	ENDIF
1623	ENDIF ! End IF (TOT_ICEnew .LE. CLIMIT)
1624	ENDIF ! End IF (ICE_logical)
1625
1626
1627	!
1628	!--- Update rain mixing ratios
1629	!
1630	!---
1631	! * TOT_RAINnew - total mass of rain after microphysics
1632	! current layer and the input flux of ice from above
1633	! * VRAIN2 - time-averaged fall speed of rain in grid and below
1634	! (with air resistance correction)
1635	! * QRnew - updated rain mixing ratio in layer
1636	! -> TOT_RAINnew=QRnew(THICK+BLDTRHVRAIN2)
1637	! * ARAINnew - updated accumulation of rain at bottom of grid cell
1638	!---
1639	!
1640	DELR=PRAUT+PRACW+PIACWR-PIACR+PIMLT+PREVP+PICND
1641	TOT_RAINnew=TOT_RAIN+THICK*DELR
1642	IF (TOT_RAIN.GT.0. .AND. TOT_RAINnew.NE.0.) THEN
1643	DUM=TOT_RAINnew/TOT_RAIN
1644	IF (DUM .LT. TOLER) TOT_RAINnew=0.
1645	ENDIF
1646	IF (TOT_RAINnew .LE. CLIMIT) THEN
1647	TOT_RAINnew=0.
1648	VRAIN2=0.
1649	QRnew=0.
1650	ARAINnew=0.
1651	ELSE
1652	!
1653	!--- 1st guess time-averaged rain rate at bottom of grid box
1654	!
1655	RR=TOT_RAINnew/(DTPH*GAMMAR)
1656	!
1657	!--- Use same algorithm as above for calculating mean drop diameter
1658	! (IDR, in microns), which is used to estimate the time-averaged
1659	! fall speed of rain drops at the bottom of the grid layer. This
1660	! isn't perfect, but the alternative is solving a transcendental
1661	! equation that is numerically inefficient and nasty to program
1662	! (coded in earlier versions of GSMCOLUMN prior to 8-22-01).
1663	!
1664	IF (RR .LE. RR_DRmin) THEN
1665	IDR=MDRmin
1666	ELSE IF (RR .LE. RR_DR1) THEN
1667	IDR=INT( 1.123E3RR*.1947 + .5 )
1668	IDR=MAX( MDRmin, MIN(IDR, MDR1) )
1669	ELSE IF (RR .LE. RR_DR2) THEN
1670	IDR=INT( 1.225E3RR*.2017 + .5 )
1671	IDR=MAX( MDR1, MIN(IDR, MDR2) )
1672	ELSE IF (RR .LE. RR_DR3) THEN
1673	IDR=INT( 1.3006E3RR*.2083 + .5 )
1674	IDR=MAX( MDR2, MIN(IDR, MDR3) )
1675	ELSE IF (RR .LE. RR_DRmax) THEN
1676	IDR=INT( 1.355E3RR*.2144 + .5 )
1677	IDR=MAX( MDR3, MIN(IDR, MDRmax) )
1678	ELSE
1679	IDR=MDRmax
1680	ENDIF ! End IF (RR .LE. RR_DRmin)
1681	VRAIN2=GAMMAR*VRAIN(IDR)
1682	QRnew=TOT_RAINnew/(THICK+BLDTRH*VRAIN2)
1683	IF (QRnew .LE. EPSQ) QRnew=0.
1684	IF (QR.GT.0. .AND. QRnew.NE.0.) THEN
1685	DUM=QRnew/QR
1686	IF (DUM .LT. TOLER) QRnew=0.
1687	ENDIF
1688	ARAINnew=BLDTRHVRAIN2QRnew
1689	IF (ARAIN.GT.0. .AND. ARAINnew.NE.0.) THEN
1690	DUM=ARAINnew/ARAIN
1691	IF (DUM .LT. TOLER) ARAINnew=0.
1692	ENDIF
1693	ENDIF
1694	!
1695	WCnew=QWnew+QRnew+QInew
1696	!
1697	!----------------------------------------------------------------------
1698	!-------------- Begin debugging & verification ------------------------
1699	!----------------------------------------------------------------------
1700	!
1701	!--- QT, QTnew - total water (vapor & condensate) before & after microphysics, resp.
1702	!
1703
1704
1705	QT=THICK*(WV+WC)+ARAIN+ASNOW
1706	QTnew=THICK*(WVnew+WCnew)+ARAINnew+ASNOWnew
1707	BUDGET=QT-QTnew
1708	!
1709	!--- Additional check on budget preservation, accounting for truncation effects
1710	!
1711	DBG_logical=.FALSE.
1712	! DUM=ABS(BUDGET)
1713	! IF (DUM .GT. TOLER) THEN
1714	! DUM=DUM/MIN(QT, QTnew)
1715	! IF (DUM .GT. TOLER) DBG_logical=.TRUE.
1716	! ENDIF
1717	!!
1718	! DUM=(RHgrd+.001)*QSInew
1719	! IF ( (QWnew.GT.EPSQ) .OR. QRnew.GT.EPSQ .OR. WVnew.GT.DUM)
1720	! & .AND. TC.LT.T_ICE ) DBG_logical=.TRUE.
1721	!
1722	! IF (TC.GT.5. .AND. QInew.GT.EPSQ) DBG_logical=.TRUE.
1723	!
1724	IF ((WVnew.LT.EPSQ .OR. DBG_logical) .AND. PRINT_diag) THEN
1725	!
1726	WRITE(6,"(/2(a,i4),2(a,i2))") '{} i=',I_index,' j=',J_index,&
1727	& ' L=',L,' LSFC=',LSFC
1728	!
1729	ESW=1000.*FPVS0(Tnew)
1730	QSWnew=EPS*ESW/(PP-ESW)
1731	IF (TC.LT.0. .OR. Tnew .LT. 0.) THEN
1732	ESI=1000.*FPVS(Tnew)
1733	QSInew=EPS*ESI/(PP-ESI)
1734	ELSE
1735	QSI=QSW
1736	QSInew=QSWnew
1737	ENDIF
1738	WSnew=QSInew
1739	WRITE(6,"(4(a12,g11.4,1x))") &
1740	& '{} TCold=',TC,'TCnew=',Tnew-T0C,'P=',.01*PP,'RHO=',RHO, &
1741	& '{} THICK=',THICK,'RHold=',WV/WS,'RHnew=',WVnew/WSnew, &
1742	& 'RHgrd=',RHgrd, &
1743	& '{} RHWold=',WV/QSW,'RHWnew=',WVnew/QSWnew,'RHIold=',WV/QSI, &
1744	& 'RHInew=',WVnew/QSInew, &
1745	& '{} QSWold=',QSW,'QSWnew=',QSWnew,'QSIold=',QSI,'QSInew=',QSInew, &
1746	& '{} WSold=',WS,'WSnew=',WSnew,'WVold=',WV,'WVnew=',WVnew, &
1747	& '{} WCold=',WC,'WCnew=',WCnew,'QWold=',QW,'QWnew=',QWnew, &
1748	& '{} QIold=',QI,'QInew=',QInew,'QRold=',QR,'QRnew=',QRnew, &
1749	& '{} ARAINold=',ARAIN,'ARAINnew=',ARAINnew,'ASNOWold=',ASNOW, &
1750	& 'ASNOWnew=',ASNOWnew, &
1751	& '{} TOT_RAIN=',TOT_RAIN,'TOT_RAINnew=',TOT_RAINnew, &
1752	& 'TOT_ICE=',TOT_ICE,'TOT_ICEnew=',TOT_ICEnew, &
1753	& '{} BUDGET=',BUDGET,'QTold=',QT,'QTnew=',QTnew
1754	!
1755	WRITE(6,"(4(a12,g11.4,1x))") &
1756	& '{} DELT=',DELT,'DELV=',DELV,'DELW=',DELW,'DELI=',DELI, &
1757	& '{} DELR=',DELR,'PCOND=',PCOND,'PIDEP=',PIDEP,'PIEVP=',PIEVP, &
1758	& '{} PICND=',PICND,'PREVP=',PREVP,'PRAUT=',PRAUT,'PRACW=',PRACW, &
1759	& '{} PIACW=',PIACW,'PIACWI=',PIACWI,'PIACWR=',PIACWR,'PIMLT=', &
1760	& PIMLT, &
1761	& '{} PIACR=',PIACR
1762	!
1763	IF (ICE_logical) WRITE(6,"(4(a12,g11.4,1x))") &
1764	& '{} RimeF1=',RimeF1,'GAMMAS=',GAMMAS,'VrimeF=',VrimeF, &
1765	& 'VSNOW=',VSNOW, &
1766	& '{} INDEXS=',FLOAT(INDEXS),'FLARGE=',FLARGE,'FSMALL=',FSMALL, &
1767	& 'FLIMASS=',FLIMASS, &
1768	& '{} XSIMASS=',XSIMASS,'XLIMASS=',XLIMASS,'QLICE=',QLICE, &
1769	& 'QTICE=',QTICE, &
1770	& '{} NLICE=',NLICE,'NSmICE=',NSmICE,'PILOSS=',PILOSS, &
1771	& 'EMAIRI=',EMAIRI, &
1772	& '{} RimeF=',RimeF
1773	!
1774	IF (TOT_RAIN.GT.0. .OR. TOT_RAINnew.GT.0.) &
1775	& WRITE(6,"(4(a12,g11.4,1x))") &
1776	& '{} INDEXR1=',FLOAT(INDEXR1),'INDEXR=',FLOAT(INDEXR), &
1777	& 'GAMMAR=',GAMMAR,'N0r=',N0r, &
1778	& '{} VRAIN1=',VRAIN1,'VRAIN2=',VRAIN2,'QTRAIN=',QTRAIN,'RQR=',RQR, &
1779	& '{} PRLOSS=',PRLOSS,'VOLR1=',THICK+BLDTRH*VRAIN1, &
1780	& 'VOLR2=',THICK+BLDTRH*VRAIN2
1781	!
1782	IF (PRAUT .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} QW0=',QW0
1783	!
1784	IF (PRACW .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} FWR=',FWR
1785	!
1786	IF (PIACR .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} FIR=',FIR
1787	!
1788	DUM=PIMLT+PICND-PREVP-PIEVP
1789	IF (DUM.GT.0. .or. DWVi.NE.0.) &
1790	& WRITE(6,"(4(a12,g11.4,1x))") &
1791	& '{} TFACTOR=',TFACTOR,'DYNVIS=',DYNVIS, &
1792	& 'THERM_CON=',THERM_COND,'DIFFUS=',DIFFUS
1793	!
1794	IF (PREVP .LT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
1795	& '{} RFACTOR=',RFACTOR,'ABW=',ABW,'VENTR=',VENTR,'CREVP=',CREVP, &
1796	& '{} DWVr=',DWVr,'DENOMW=',DENOMW
1797	!
1798	IF (PIDEP.NE.0. .AND. DWVi.NE.0.) &
1799	& WRITE(6,"(4(a12,g11.4,1x))") &
1800	& '{} DWVi=',DWVi,'DENOMI=',DENOMI,'PIDEP_max=',PIDEP_max, &
1801	& 'SFACTOR=',SFACTOR, &
1802	& '{} ABI=',ABI,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
1803	& 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
1804	& '{} VENTIS=',VENTIS,'DIDEP=',DIDEP
1805	!
1806	IF (PIDEP.GT.0. .AND. PCOND.NE.0.) &
1807	& WRITE(6,"(4(a12,g11.4,1x))") &
1808	& '{} DENOMW=',DENOMW,'DENOMWI=',DENOMWI,'DENOMF=',DENOMF, &
1809	& 'DUM2=',PCOND-PIACW
1810	!
1811	IF (FWS .GT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
1812	& '{} FWS=',FWS
1813	!
1814	DUM=PIMLT+PICND-PIEVP
1815	IF (DUM.GT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
1816	& '{} SFACTOR=',SFACTOR,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
1817	& 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
1818	& '{} AIEVP=',AIEVP,'DIEVP=',DIEVP,'QSW0=',QSW0,'DWV0=',DWV0
1819	!
1820	ENDIF
1821
1822
1823	!
1824	!-----------------------------------------------------------------------
1825	!--------------- Water budget statistics & maximum values --------------
1826	!-----------------------------------------------------------------------
1827	!
1828	IF (PRINT_diag) THEN
1829	ITdx=MAX( ITLO, MIN( INT(Tnew-T0C), ITHI ) )
1830	IF (QInew .GT. EPSQ) NSTATS(ITdx,1)=NSTATS(ITdx,1)+1
1831	IF (QInew.GT.EPSQ .AND. QRnew+QWnew.GT.EPSQ) &
1832	& NSTATS(ITdx,2)=NSTATS(ITdx,2)+1
1833	IF (QWnew .GT. EPSQ) NSTATS(ITdx,3)=NSTATS(ITdx,3)+1
1834	IF (QRnew .GT. EPSQ) NSTATS(ITdx,4)=NSTATS(ITdx,4)+1
1835	!
1836	QMAX(ITdx,1)=MAX(QMAX(ITdx,1), QInew)
1837	QMAX(ITdx,2)=MAX(QMAX(ITdx,2), QWnew)
1838	QMAX(ITdx,3)=MAX(QMAX(ITdx,3), QRnew)
1839	QMAX(ITdx,4)=MAX(QMAX(ITdx,4), ASNOWnew)
1840	QMAX(ITdx,5)=MAX(QMAX(ITdx,5), ARAINnew)
1841	QTOT(ITdx,1)=QTOT(ITdx,1)+QInew*THICK
1842	QTOT(ITdx,2)=QTOT(ITdx,2)+QWnew*THICK
1843	QTOT(ITdx,3)=QTOT(ITdx,3)+QRnew*THICK
1844	!
1845	QTOT(ITdx,4)=QTOT(ITdx,4)+PCOND*THICK
1846	QTOT(ITdx,5)=QTOT(ITdx,5)+PICND*THICK
1847	QTOT(ITdx,6)=QTOT(ITdx,6)+PIEVP*THICK
1848	QTOT(ITdx,7)=QTOT(ITdx,7)+PIDEP*THICK
1849	QTOT(ITdx,8)=QTOT(ITdx,8)+PREVP*THICK
1850	QTOT(ITdx,9)=QTOT(ITdx,9)+PRAUT*THICK
1851	QTOT(ITdx,10)=QTOT(ITdx,10)+PRACW*THICK
1852	QTOT(ITdx,11)=QTOT(ITdx,11)+PIMLT*THICK
1853	QTOT(ITdx,12)=QTOT(ITdx,12)+PIACW*THICK
1854	QTOT(ITdx,13)=QTOT(ITdx,13)+PIACWI*THICK
1855	QTOT(ITdx,14)=QTOT(ITdx,14)+PIACWR*THICK
1856	QTOT(ITdx,15)=QTOT(ITdx,15)+PIACR*THICK
1857	!
1858	QTOT(ITdx,16)=QTOT(ITdx,16)+(WVnew-WV)*THICK
1859	QTOT(ITdx,17)=QTOT(ITdx,17)+(QWnew-QW)*THICK
1860	QTOT(ITdx,18)=QTOT(ITdx,18)+(QInew-QI)*THICK
1861	QTOT(ITdx,19)=QTOT(ITdx,19)+(QRnew-QR)*THICK
1862	QTOT(ITdx,20)=QTOT(ITdx,20)+(ARAINnew-ARAIN)
1863	QTOT(ITdx,21)=QTOT(ITdx,21)+(ASNOWnew-ASNOW)
1864	IF (QInew .GT. 0.) &
1865	& QTOT(ITdx,22)=QTOT(ITdx,22)+QInew*THICK/RimeF
1866	!
1867	ENDIF
1868	!
1869	!----------------------------------------------------------------------
1870	!------------------------- Update arrays ------------------------------
1871	!----------------------------------------------------------------------
1872	!
1873
1874
1875	T_col(L)=Tnew ! Updated temperature
1876	!
1877	QV_col(L)=max(EPSQ, WVnew/(1.+WVnew)) ! Updated specific humidity
1878	WC_col(L)=max(EPSQ, WCnew) ! Updated total condensate mixing ratio
1879	QI_col(L)=max(EPSQ, QInew) ! Updated ice mixing ratio
1880	QR_col(L)=max(EPSQ, QRnew) ! Updated rain mixing ratio
1881	QW_col(L)=max(EPSQ, QWnew) ! Updated cloud water mixing ratio
1882	RimeF_col(L)=RimeF ! Updated rime factor
1883	ASNOW=ASNOWnew ! Updated accumulated snow
1884	ARAIN=ARAINnew ! Updated accumulated rain
1885	!
1886	!#######################################################################
1887	!
1888	10 CONTINUE ! ##### End "L" loop through model levels #####
1889
1890
1891	!
1892	!#######################################################################
1893	!
1894	!-----------------------------------------------------------------------
1895	!--------------------------- Return to GSMDRIVE -----------------------
1896	!-----------------------------------------------------------------------
1897	!
1898	CONTAINS
1899	!#######################################################################
1900	!--------- Produces accurate calculation of cloud condensation ---------
1901	!#######################################################################
1902	!
1903	REAL FUNCTION CONDENSE (PP, QW, TK, WV)
1904	!
1905	!---------------------------------------------------------------------------------
1906	!------ The Asai (1965) algorithm takes into consideration the release of ------
1907	!------ latent heat in increasing the temperature & in increasing the ------
1908	!------ saturation mixing ratio (following the Clausius-Clapeyron eqn.). ------
1909	!---------------------------------------------------------------------------------
1910	!
1911	IMPLICIT NONE
1912	!
1913	INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
1914	REAL (KIND=HIGH_PRES), PARAMETER :: &
1915	& RHLIMIT=.001, RHLIMIT1=-RHLIMIT
1916	REAL (KIND=HIGH_PRES) :: COND, SSAT, WCdum
1917	!
1918	REAL,INTENT(IN) :: QW,PP,WV,TK
1919	REAL WVdum,Tdum,XLV2,DWV,WS,ESW,XLV1,XLV
1920	integer nsteps
1921	!
1922	!-----------------------------------------------------------------------
1923	!
1924	!--- LV (T) is from Bolton (JAS, 1980)
1925	!
1926	XLV=3.148E6-2370.*TK
1927	XLV1=XLV*RCP
1928	XLV2=XLVXLVRCPRV
1929	Tdum=TK
1930	WVdum=WV
1931	WCdum=QW
1932	ESW=1000.*FPVS0(Tdum) ! Saturation vapor press w/r/t water
1933	WS=RHgrdEPSESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
1934	DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1935	SSAT=DWV/WS ! Supersaturation ratio
1936	CONDENSE=0.
1937	nsteps = 0
1938	DO WHILE ((SSAT.LT.RHLIMIT1 .AND. WCdum.GT.EPSQ) &
1939	& .OR. SSAT.GT.RHLIMIT)
1940	nsteps = nsteps + 1
1941	COND=DWV/(1.+XLV2WS/(TdumTdum)) ! Asai (1965, J. Japan)
1942	COND=MAX(COND, -WCdum) ! Limit cloud water evaporation
1943	Tdum=Tdum+XLV1*COND ! Updated temperature
1944	WVdum=WVdum-COND ! Updated water vapor mixing ratio
1945	WCdum=WCdum+COND ! Updated cloud water mixing ratio
1946	CONDENSE=CONDENSE+COND ! Total cloud water condensation
1947	ESW=1000.*FPVS0(Tdum) ! Updated saturation vapor press w/r/t water
1948	WS=RHgrdEPSESW/(PP-ESW) ! Updated saturation mixing ratio w/r/t water
1949	DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1950	SSAT=DWV/WS ! Grid-scale supersaturation ratio
1951	ENDDO
1952	!
1953	END FUNCTION CONDENSE
1954	!
1955	!#######################################################################
1956	!---------------- Calculate ice deposition at T<T_ICE ------------------
1957	!#######################################################################
1958	!
1959	REAL FUNCTION DEPOSIT (PP, Tdum, WVdum)
1960	!
1961	!--- Also uses the Asai (1965) algorithm, but uses a different target
1962	! vapor pressure for the adjustment
1963	!
1964	IMPLICIT NONE
1965	!
1966	INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
1967	REAL (KIND=HIGH_PRES), PARAMETER :: RHLIMIT=.001, &
1968	& RHLIMIT1=-RHLIMIT
1969	REAL (KIND=HIGH_PRES) :: DEP, SSAT
1970	!
1971	real,INTENT(IN) :: PP
1972	real,INTENT(INOUT) :: WVdum,Tdum
1973	real ESI,WS,DWV
1974	!
1975	!-----------------------------------------------------------------------
1976	!
1977	ESI=1000.*FPVS(Tdum) ! Saturation vapor press w/r/t ice
1978	WS=RHgrdEPSESI/(PP-ESI) ! Saturation mixing ratio
1979	DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1980	SSAT=DWV/WS ! Supersaturation ratio
1981	DEPOSIT=0.
1982	DO WHILE (SSAT.GT.RHLIMIT .OR. SSAT.LT.RHLIMIT1)
1983	!
1984	!--- Note that XLVS2=LSLV/(CPRV)=LVWS/(RVTT)(LS/CP*DEP1),
1985	! where WS is the saturation mixing ratio following Clausius-
1986	! Clapeyron (see Asai,1965; Young,1993,p.405)
1987	!
1988	DEP=DWV/(1.+XLS2WS/(TdumTdum)) ! Asai (1965, J. Japan)
1989	Tdum=Tdum+XLS1*DEP ! Updated temperature
1990	WVdum=WVdum-DEP ! Updated ice mixing ratio
1991	DEPOSIT=DEPOSIT+DEP ! Total ice deposition
1992	ESI=1000.*FPVS(Tdum) ! Updated saturation vapor press w/r/t ice
1993	WS=RHgrdEPSESI/(PP-ESI) ! Updated saturation mixing ratio w/r/t ice
1994	DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1995	SSAT=DWV/WS ! Grid-scale supersaturation ratio
1996	ENDDO
1997	!
1998	END FUNCTION DEPOSIT
1999	!
2000	END SUBROUTINE EGCP01COLUMN
2001
2002	!#######################################################################
2003	!------- Initialize constants & lookup tables for microphysics ---------
2004	!#######################################################################
2005	!
2006
2007	! SH 0211/2002
2008
2009	!-----------------------------------------------------------------------
2010	SUBROUTINE etanewinit (GSMDT,DT,DELX,DELY,LOWLYR,restart, &
2011	& F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
2012	& MP_RESTART_STATE,TBPVS_STATE,TBPVS0_STATE, &
2013	& ALLOWED_TO_READ, &
2014	& IDS,IDE,JDS,JDE,KDS,KDE, &
2015	& IMS,IME,JMS,JME,KMS,KME, &
2016	& ITS,ITE,JTS,JTE,KTS,KTE )
2017	!-----------------------------------------------------------------------
2018	!-------------------------------------------------------------------------------
2019	!--- SUBPROGRAM DOCUMENTATION BLOCK
2020	! PRGRMMR: Ferrier ORG: W/NP22 DATE: February 2001
2021	! Jin ORG: W/NP22 DATE: January 2002
2022	! (Modification for WRF structure)
2023	!
2024	!-------------------------------------------------------------------------------
2025	! ABSTRACT:
2026	! * Reads various microphysical lookup tables used in COLUMN_MICRO
2027	! * Lookup tables were created "offline" and are read in during execution
2028	! * Creates lookup tables for saturation vapor pressure w/r/t water & ice
2029	!-------------------------------------------------------------------------------
2030	!
2031	! USAGE: CALL etanewinit FROM SUBROUTINE GSMDRIVE AT MODEL START TIME
2032	!
2033	! INPUT ARGUMENT LIST:
2034	! DTPH - physics time step (s)
2035	!
2036	! OUTPUT ARGUMENT LIST:
2037	! NONE
2038	!
2039	! OUTPUT FILES:
2040	! NONE
2041	!
2042	! SUBROUTINES:
2043	! MY_GROWTH_RATES - lookup table for growth of nucleated ice
2044	! GPVS - lookup tables for saturation vapor pressure (water, ice)
2045	!
2046	! UNIQUE: NONE
2047	!
2048	! LIBRARY: NONE
2049	!
2050	! COMMON BLOCKS:
2051	! CMICRO_CONS - constants used in GSMCOLUMN
2052	! CMY600 - lookup table for growth of ice crystals in
2053	! water saturated conditions (Miller & Young, 1979)
2054	! IVENT_TABLES - lookup tables for ventilation effects of ice
2055	! IACCR_TABLES - lookup tables for accretion rates of ice
2056	! IMASS_TABLES - lookup tables for mass content of ice
2057	! IRATE_TABLES - lookup tables for precipitation rates of ice
2058	! IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
2059	! MAPOT - Need lat/lon grid resolution
2060	! RVENT_TABLES - lookup tables for ventilation effects of rain
2061	! RACCR_TABLES - lookup tables for accretion rates of rain
2062	! RMASS_TABLES - lookup tables for mass content of rain
2063	! RVELR_TABLES - lookup tables for fall speeds of rain
2064	! RRATE_TABLES - lookup tables for precipitation rates of rain
2065	!
2066	! ATTRIBUTES:
2067	! LANGUAGE: FORTRAN 90
2068	! MACHINE : IBM SP
2069	!
2070	!-----------------------------------------------------------------------
2071	!
2072	!
2073	!-----------------------------------------------------------------------
2074	IMPLICIT NONE
2075	!-----------------------------------------------------------------------
2076	!-------------------------------------------------------------------------
2077	!-------------- Parameters & arrays for lookup tables --------------------
2078	!-------------------------------------------------------------------------
2079	!
2080	!--- Common block of constants used in column microphysics
2081	!
2082	!WRF
2083	! real DLMD,DPHD
2084	!WRF
2085	!
2086	!-----------------------------------------------------------------------
2087	!--- Parameters & data statement for local calculations
2088	!-----------------------------------------------------------------------
2089	!
2090	INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3
2091	!
2092	! VARIABLES PASSED IN
2093	integer,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
2094	& ,IMS,IME,JMS,JME,KMS,KME &
2095	& ,ITS,ITE,JTS,JTE,KTS,KTE
2096	!WRF
2097	INTEGER, DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: LOWLYR
2098	!
2099	real, INTENT(IN) :: DELX,DELY
2100	real,DIMENSION(*), INTENT(INOUT) :: MP_RESTART_STATE
2101	real,DIMENSION(NX), INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
2102	real,DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(OUT) :: &
2103	& F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
2104	INTEGER, PARAMETER :: ITLO=-60, ITHI=40
2105	! integer,DIMENSION(ITLO:ITHI,4),INTENT(INOUT) :: NSTATS
2106	! real,DIMENSION(ITLO:ITHI,5),INTENT(INOUT) :: QMAX
2107	! real,DIMENSION(ITLO:ITHI,22),INTENT(INOUT) :: QTOT
2108	! real,INTENT(INOUT) :: PRECtot(2),PRECmax(2)
2109	real,INTENT(IN) :: DT,GSMDT
2110	LOGICAL,INTENT(IN) :: allowed_to_read,restart
2111	!
2112	!-----------------------------------------------------------------------
2113	! LOCAL VARIABLES
2114	!-----------------------------------------------------------------------
2115	REAL :: BBFR,DTPH,PI,DX,Thour_print
2116	INTEGER :: I,IM,J,L,K,JTF,KTF,ITF
2117	INTEGER :: etampnew_unit1
2118	LOGICAL, PARAMETER :: PRINT_diag=.FALSE.
2119	LOGICAL :: opened
2120	LOGICAL , EXTERNAL :: wrf_dm_on_monitor
2121	CHARACTER*80 errmess
2122	!
2123	!-----------------------------------------------------------------------
2124	!
2125	JTF=MIN0(JTE,JDE-1)
2126	KTF=MIN0(KTE,KDE-1)
2127	ITF=MIN0(ITE,IDE-1)
2128	!
2129	DO J=JTS,JTF
2130	DO I=ITS,ITF
2131	LOWLYR(I,J)=1
2132	ENDDO
2133	ENDDO
2134	!
2135	IF(.NOT.RESTART)THEN
2136	DO J = jts,jte
2137	DO K = kts,kte
2138	DO I= its,ite
2139	F_ICE_PHY(i,k,j)=0.
2140	F_RAIN_PHY(i,k,j)=0.
2141	F_RIMEF_PHY(i,k,j)=1.
2142	ENDDO
2143	ENDDO
2144	ENDDO
2145	ENDIF
2146	!
2147	!-----------------------------------------------------------------------
2148	IF(ALLOWED_TO_READ)THEN
2149	!-----------------------------------------------------------------------
2150	!
2151	DX=((DELX)2+(DELY)2)**.5/1000. ! Model resolution at equator (km)
2152	DX=MIN(100., MAX(5., DX) )
2153	!
2154	!-- Relative humidity threshold for the onset of grid-scale condensation
2155	!!-- 9/1/01: Assume the following functional dependence for 5 - 100 km resolution:
2156	!! RHgrd=0.90 for dx=100 km, 0.98 for dx=5 km, where
2157	! RHgrd=0.90+.08((100.-DX)/95.)*.5
2158	!
2159	DTPH=MAX(GSMDT*60.,DT)
2160	DTPH=NINT(DTPH/DT)*DT
2161	!
2162	!--- Create lookup tables for saturation vapor pressure w/r/t water & ice
2163	!
2164	CALL GPVS
2165	!
2166	!--- Read in various lookup tables
2167	!
2168	IF ( wrf_dm_on_monitor() ) THEN
2169	DO i = 31,99
2170	INQUIRE ( i , OPENED = opened )
2171	IF ( .NOT. opened ) THEN
2172	etampnew_unit1 = i
2173	GOTO 2061
2174	ENDIF
2175	ENDDO
2176	etampnew_unit1 = -1
2177	2061 CONTINUE
2178	ENDIF
2179	!
2180	CALL wrf_dm_bcast_bytes ( etampnew_unit1 , IWORDSIZE )
2181	!
2182	IF ( etampnew_unit1 < 0 ) THEN
2183	CALL wrf_error_fatal ( 'module_mp_etanew: etanewinit: Can not find '// &
2184	'unused fortran unit to read in lookup table.' )
2185	ENDIF
2186	!
2187	IF ( wrf_dm_on_monitor() ) THEN
2188	!!was OPEN (UNIT=1,FILE="eta_micro_lookup.dat",FORM="UNFORMATTED")
2189	OPEN(UNIT=etampnew_unit1,FILE="ETAMPNEW_DATA", &
2190	& FORM="UNFORMATTED",STATUS="OLD",ERR=9061)
2191	!
2192	READ(etampnew_unit1) VENTR1
2193	READ(etampnew_unit1) VENTR2
2194	READ(etampnew_unit1) ACCRR
2195	READ(etampnew_unit1) MASSR
2196	READ(etampnew_unit1) VRAIN
2197	READ(etampnew_unit1) RRATE
2198	READ(etampnew_unit1) VENTI1
2199	READ(etampnew_unit1) VENTI2
2200	READ(etampnew_unit1) ACCRI
2201	READ(etampnew_unit1) MASSI
2202	READ(etampnew_unit1) VSNOWI
2203	READ(etampnew_unit1) VEL_RF
2204	! read(etampnew_unit1) my_growth ! Applicable only for DTPH=180 s
2205	CLOSE (etampnew_unit1)
2206	ENDIF
2207	!
2208	CALL wrf_dm_bcast_bytes ( VENTR1 , size ( VENTR1 ) * RWORDSIZE )
2209	CALL wrf_dm_bcast_bytes ( VENTR2 , size ( VENTR2 ) * RWORDSIZE )
2210	CALL wrf_dm_bcast_bytes ( ACCRR , size ( ACCRR ) * RWORDSIZE )
2211	CALL wrf_dm_bcast_bytes ( MASSR , size ( MASSR ) * RWORDSIZE )
2212	CALL wrf_dm_bcast_bytes ( VRAIN , size ( VRAIN ) * RWORDSIZE )
2213	CALL wrf_dm_bcast_bytes ( RRATE , size ( RRATE ) * RWORDSIZE )
2214	CALL wrf_dm_bcast_bytes ( VENTI1 , size ( VENTI1 ) * RWORDSIZE )
2215	CALL wrf_dm_bcast_bytes ( VENTI2 , size ( VENTI2 ) * RWORDSIZE )
2216	CALL wrf_dm_bcast_bytes ( ACCRI , size ( ACCRI ) * RWORDSIZE )
2217	CALL wrf_dm_bcast_bytes ( MASSI , size ( MASSI ) * RWORDSIZE )
2218	CALL wrf_dm_bcast_bytes ( VSNOWI , size ( VSNOWI ) * RWORDSIZE )
2219	CALL wrf_dm_bcast_bytes ( VEL_RF , size ( VEL_RF ) * RWORDSIZE )
2220	!
2221	!--- Calculates coefficients for growth rates of ice nucleated in water
2222	! saturated conditions, scaled by physics time step (lookup table)
2223	!
2224	CALL MY_GROWTH_RATES (DTPH)
2225	! CALL MY_GROWTH_RATES (DTPH,MY_GROWTH)
2226	!
2227	PI=ACOS(-1.)
2228	!
2229	!--- Constants associated with Biggs (1953) freezing of rain, as parameterized
2230	! following Lin et al. (JCAM, 1983) & Reisner et al. (1998, QJRMS).
2231	!
2232	ABFR=-0.66
2233	BBFR=100.
2234	CBFR=20.PIPIBBFRRHOL*1.E-21
2235	!
2236	!--- CIACW is used in calculating riming rates
2237	! The assumed effective collection efficiency of cloud water rimed onto
2238	! ice is =0.5 below:
2239	!
2240	CIACW=DTPH0.25PI0.5(1.E5)**C1
2241	!
2242	!--- CIACR is used in calculating freezing of rain colliding with large ice
2243	! The assumed collection efficiency is 1.0
2244	!
2245	CIACR=PI*DTPH
2246	!
2247	!--- Based on rain lookup tables for mean diameters from 0.05 to 0.45 mm
2248	! * Four different functional relationships of mean drop diameter as
2249	! a function of rain rate (RR), derived based on simple fits to
2250	! mass-weighted fall speeds of rain as functions of mean diameter
2251	! from the lookup tables.
2252	!
2253	RR_DRmin=N0r0*RRATE(MDRmin) ! RR for mean drop diameter of .05 mm
2254	RR_DR1=N0r0*RRATE(MDR1) ! RR for mean drop diameter of .10 mm
2255	RR_DR2=N0r0*RRATE(MDR2) ! RR for mean drop diameter of .20 mm
2256	RR_DR3=N0r0*RRATE(MDR3) ! RR for mean drop diameter of .32 mm
2257	RR_DRmax=N0r0*RRATE(MDRmax) ! RR for mean drop diameter of .45 mm
2258	!
2259	RQR_DRmin=N0r0*MASSR(MDRmin) ! Rain content for mean drop diameter of .05 mm
2260	RQR_DR1=N0r0*MASSR(MDR1) ! Rain content for mean drop diameter of .10 mm
2261	RQR_DR2=N0r0*MASSR(MDR2) ! Rain content for mean drop diameter of .20 mm
2262	RQR_DR3=N0r0*MASSR(MDR3) ! Rain content for mean drop diameter of .32 mm
2263	RQR_DRmax=N0r0*MASSR(MDRmax) ! Rain content for mean drop diameter of .45 mm
2264	C_N0r0=PIRHOLN0r0
2265	CN0r0=1.E6/C_N0r0**.25
2266	CN0r_DMRmin=1./(PIRHOLDMRmin**4)
2267	CN0r_DMRmax=1./(PIRHOLDMRmax**4)
2268	!
2269	!--- CRACW is used in calculating collection of cloud water by rain (an
2270	! assumed collection efficiency of 1.0)
2271	!
2272	CRACW=DTPH0.25PI*1.0
2273	!
2274	ESW0=1000.*FPVS0(T0C) ! Saturation vapor pressure at 0C
2275	RFmax=1.1**Nrime ! Maximum rime factor allowed
2276	!
2277	!------------------------------------------------------------------------
2278	!--------------- Constants passed through argument list -----------------
2279	!------------------------------------------------------------------------
2280	!
2281	!--- Important parameters for self collection (autoconversion) of
2282	! cloud water to rain.
2283	!
2284	!--- CRAUT is proportional to the rate that cloud water is converted by
2285	! self collection to rain (autoconversion rate)
2286	!
2287	CRAUT=1.-EXP(-1.E-3*DTPH)
2288	!
2289	!--- QAUT0 is the threshold cloud content for autoconversion to rain
2290	! needed for droplets to reach a diameter of 20 microns (following
2291	! Manton and Cotton, 1977; Banta and Hanson, 1987, JCAM)
2292	!--- QAUT0=1.2567, 0.8378, or 0.4189 g/m**3 for droplet number concentrations
2293	! of 300, 200, and 100 cm**-3, respectively
2294	!
2295	QAUT0=PIRHOLNCW(20.E-6)*3/6.
2296	!
2297	!--- For calculating snow optical depths by considering bulk density of
2298	! snow based on emails from Q. Fu (6/27-28/01), where optical
2299	! depth (T) = 1.5SWP/(ReffDENS), SWP is snow water path, Reff
2300	! is effective radius, and DENS is the bulk density of snow.
2301	!
2302	! SWP (kg/m*2)=(1.E-3 kg/g)SWPrad, SWPrad in g/m**2 used in radiation
2303	! T = 1.51.E3SWPrad/(Reff*DENS)
2304	!
2305	! See derivation for MASSI(INDEXS), note equal to RHO*QSNOW/NSNOW
2306	!
2307	! SDENS=1.5e3/DENS, DENS=MASSI(INDEXS)/[PI(1.E-6INDEXS)**3]
2308	!
2309	DO I=MDImin,MDImax
2310	SDENS(I)=PI1.5E-15FLOAT(III)/MASSI(I)
2311	ENDDO
2312	!
2313	Thour_print=-DTPH/3600.
2314
2315	! SH 0211/2002
2316	! IF (PRINT_diag) THEN
2317
2318	!-------- Total and maximum quantities
2319	!
2320	! NSTATS=0 ! Microphysical statistics dealing w/ grid-point counts
2321	! QMAX=0. ! Microphysical statistics dealing w/ hydrometeor mass
2322	! QTOT=0. ! Microphysical statistics dealing w/ hydrometeor mass
2323	! PRECmax=0. ! Maximum precip rates (rain, snow) at surface (mm/h)
2324	! PRECtot=0. ! Total precipitation (rain, snow) accumulation at surface
2325	! ENDIF
2326
2327	!wrf
2328	IF(.NOT.RESTART)THEN
2329	MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
2330	MP_RESTART_STATE(MY_T2+1)=C1XPVS0
2331	MP_RESTART_STATE(MY_T2+2)=C2XPVS0
2332	MP_RESTART_STATE(MY_T2+3)=C1XPVS
2333	MP_RESTART_STATE(MY_T2+4)=C2XPVS
2334	MP_RESTART_STATE(MY_T2+5)=CIACW
2335	MP_RESTART_STATE(MY_T2+6)=CIACR
2336	MP_RESTART_STATE(MY_T2+7)=CRACW
2337	MP_RESTART_STATE(MY_T2+8)=CRAUT
2338	TBPVS_STATE(1:NX) =TBPVS(1:NX)
2339	TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
2340	ENDIF
2341
2342	ENDIF ! Allowed_to_read
2343
2344	RETURN
2345	!
2346	!-----------------------------------------------------------------------
2347	!
2348	9061 CONTINUE
2349	WRITE( errmess , '(A,I4)' ) &
2350	'module_mp_etanew: error opening ETAMPNEW_DATA on unit ' &
2351	&, etampnew_unit1
2352	CALL wrf_error_fatal(errmess)
2353	!
2354	!-----------------------------------------------------------------------
2355	END SUBROUTINE etanewinit
2356	!
2357	SUBROUTINE MY_GROWTH_RATES (DTPH)
2358	! SUBROUTINE MY_GROWTH_RATES (DTPH,MY_GROWTH)
2359	!
2360	!--- Below are tabulated values for the predicted mass of ice crystals
2361	! after 600 s of growth in water saturated conditions, based on
2362	! calculations from Miller and Young (JAS, 1979). These values are
2363	! crudely estimated from tabulated curves at 600 s from Fig. 6.9 of
2364	! Young (1993). Values at temperatures colder than -27C were
2365	! assumed to be invariant with temperature.
2366	!
2367	!--- Used to normalize Miller & Young (1979) calculations of ice growth
2368	! over large time steps using their tabulated values at 600 s.
2369	! Assumes 3D growth with time**1.5 following eq. (6.3) in Young (1993).
2370	!
2371	IMPLICIT NONE
2372	!
2373	REAL,INTENT(IN) :: DTPH
2374	!
2375	REAL DT_ICE
2376	REAL,DIMENSION(35) :: MY_600
2377	!WRF
2378	!
2379	!-----------------------------------------------------------------------
2380	DATA MY_600 / &
2381	& 5.5e-8, 1.4E-7, 2.8E-7, 6.E-7, 3.3E-6, &
2382	& 2.E-6, 9.E-7, 8.8E-7, 8.2E-7, 9.4e-7, &
2383	& 1.2E-6, 1.85E-6, 5.5E-6, 1.5E-5, 1.7E-5, &
2384	& 1.5E-5, 1.E-5, 3.4E-6, 1.85E-6, 1.35E-6, &
2385	& 1.05E-6, 1.E-6, 9.5E-7, 9.0E-7, 9.5E-7, &
2386	& 9.5E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7, &
2387	& 9.E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7 / ! -31 to -35 deg C
2388	!
2389	!-----------------------------------------------------------------------
2390	!
2391	DT_ICE=(DTPH/600.)**1.5
2392	MY_GROWTH=DT_ICE*MY_600
2393	!
2394	!-----------------------------------------------------------------------
2395	!
2396	END SUBROUTINE MY_GROWTH_RATES
2397	!
2398	!-----------------------------------------------------------------------
2399	!--------- Old GFS saturation vapor pressure lookup tables -----------
2400	!-----------------------------------------------------------------------
2401	!
2402	SUBROUTINE GPVS
2403	! ******************************************************************
2404	!$$$ SUBPROGRAM DOCUMENTATION BLOCK
2405	! . . .
2406	! SUBPROGRAM: GPVS COMPUTE SATURATION VAPOR PRESSURE TABLE
2407	! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2408	!
2409	! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE TABLE AS A FUNCTION OF
2410	! TEMPERATURE FOR THE TABLE LOOKUP FUNCTION FPVS.
2411	! EXACT SATURATION VAPOR PRESSURES ARE CALCULATED IN SUBPROGRAM FPVSX.
2412	! THE CURRENT IMPLEMENTATION COMPUTES A TABLE WITH A LENGTH
2413	! OF 7501 FOR TEMPERATURES RANGING FROM 180.0 TO 330.0 KELVIN.
2414	!
2415	! PROGRAM HISTORY LOG:
2416	! 91-05-07 IREDELL
2417	! 94-12-30 IREDELL EXPAND TABLE
2418	! 96-02-19 HONG ICE EFFECT
2419	! 01-11-29 JIN MODIFIED FOR WRF
2420	!
2421	! USAGE: CALL GPVS
2422	!
2423	! SUBPROGRAMS CALLED:
2424	! (FPVSX) - INLINABLE FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE
2425	!
2426	! COMMON BLOCKS:
2427	! COMPVS - SCALING PARAMETERS AND TABLE FOR FUNCTION FPVS.
2428	!
2429	! ATTRIBUTES:
2430	! LANGUAGE: FORTRAN 90
2431	!
2432	!$$$
2433	IMPLICIT NONE
2434	real :: X,XINC,T
2435	integer :: JX
2436	!----------------------------------------------------------------------
2437	XINC=(XMAX-XMIN)/(NX-1)
2438	C1XPVS=1.-XMIN/XINC
2439	C2XPVS=1./XINC
2440	C1XPVS0=1.-XMIN/XINC
2441	C2XPVS0=1./XINC
2442	!
2443	DO JX=1,NX
2444	X=XMIN+(JX-1)*XINC
2445	T=X
2446	TBPVS(JX)=FPVSX(T)
2447	TBPVS0(JX)=FPVSX0(T)
2448	ENDDO
2449	!
2450	END SUBROUTINE GPVS
2451	!-----------------------------------------------------------------------
2452	!***********************************************************************
2453	!-----------------------------------------------------------------------
2454	REAL FUNCTION FPVS(T)
2455	!-----------------------------------------------------------------------
2456	!$$$ SUBPROGRAM DOCUMENTATION BLOCK
2457	! . . .
2458	! SUBPROGRAM: FPVS COMPUTE SATURATION VAPOR PRESSURE
2459	! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2460	!
2461	! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE FROM THE TEMPERATURE.
2462	! A LINEAR INTERPOLATION IS DONE BETWEEN VALUES IN A LOOKUP TABLE
2463	! COMPUTED IN GPVS. SEE DOCUMENTATION FOR FPVSX FOR DETAILS.
2464	! INPUT VALUES OUTSIDE TABLE RANGE ARE RESET TO TABLE EXTREMA.
2465	! THE INTERPOLATION ACCURACY IS ALMOST 6 DECIMAL PLACES.
2466	! ON THE CRAY, FPVS IS ABOUT 4 TIMES FASTER THAN EXACT CALCULATION.
2467	! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
2468	!
2469	! PROGRAM HISTORY LOG:
2470	! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
2471	! 94-12-30 IREDELL EXPAND TABLE
2472	! 96-02-19 HONG ICE EFFECT
2473	! 01-11-29 JIN MODIFIED FOR WRF
2474	!
2475	! USAGE: PVS=FPVS(T)
2476	!
2477	! INPUT ARGUMENT LIST:
2478	! T - REAL TEMPERATURE IN KELVIN
2479	!
2480	! OUTPUT ARGUMENT LIST:
2481	! FPVS - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
2482	!
2483	! ATTRIBUTES:
2484	! LANGUAGE: FORTRAN 90
2485	!$$$
2486	IMPLICIT NONE
2487	real,INTENT(IN) :: T
2488	real XJ
2489	integer :: JX
2490	!-----------------------------------------------------------------------
2491	XJ=MIN(MAX(C1XPVS+C2XPVS*T,1.),FLOAT(NX))
2492	JX=MIN(XJ,NX-1.)
2493	FPVS=TBPVS(JX)+(XJ-JX)*(TBPVS(JX+1)-TBPVS(JX))
2494	!
2495	END FUNCTION FPVS
2496	!-----------------------------------------------------------------------
2497	!-----------------------------------------------------------------------
2498	REAL FUNCTION FPVS0(T)
2499	!-----------------------------------------------------------------------
2500	IMPLICIT NONE
2501	real,INTENT(IN) :: T
2502	real :: XJ1
2503	integer :: JX1
2504	!-----------------------------------------------------------------------
2505	XJ1=MIN(MAX(C1XPVS0+C2XPVS0*T,1.),FLOAT(NX))
2506	JX1=MIN(XJ1,NX-1.)
2507	FPVS0=TBPVS0(JX1)+(XJ1-JX1)*(TBPVS0(JX1+1)-TBPVS0(JX1))
2508	!
2509	END FUNCTION FPVS0
2510	!-----------------------------------------------------------------------
2511	!***********************************************************************
2512	!-----------------------------------------------------------------------
2513	REAL FUNCTION FPVSX(T)
2514	!-----------------------------------------------------------------------
2515	!$$$ SUBPROGRAM DOCUMENTATION BLOCK
2516	! . . .
2517	! SUBPROGRAM: FPVSX COMPUTE SATURATION VAPOR PRESSURE
2518	! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2519	!
2520	! ABSTRACT: EXACTLY COMPUTE SATURATION VAPOR PRESSURE FROM TEMPERATURE.
2521	! THE WATER MODEL ASSUMES A PERFECT GAS, CONSTANT SPECIFIC HEATS
2522	! FOR GAS AND LIQUID, AND NEGLECTS THE VOLUME OF THE LIQUID.
2523	! THE MODEL DOES ACCOUNT FOR THE VARIATION OF THE LATENT HEAT
2524	! OF CONDENSATION WITH TEMPERATURE. THE ICE OPTION IS NOT INCLUDED.
2525	! THE CLAUSIUS-CLAPEYRON EQUATION IS INTEGRATED FROM THE TRIPLE POINT
2526	! TO GET THE FORMULA
2527	! PVS=PSATK(TRXA)EXP(XB*(1.-TR))
2528	! WHERE TR IS TTP/T AND OTHER VALUES ARE PHYSICAL CONSTANTS
2529	! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
2530	!
2531	! PROGRAM HISTORY LOG:
2532	! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
2533	! 94-12-30 IREDELL EXACT COMPUTATION
2534	! 96-02-19 HONG ICE EFFECT
2535	! 01-11-29 JIN MODIFIED FOR WRF
2536	!
2537	! USAGE: PVS=FPVSX(T)
2538	! REFERENCE: EMANUEL(1994),116-117
2539	!
2540	! INPUT ARGUMENT LIST:
2541	! T - REAL TEMPERATURE IN KELVIN
2542	!
2543	! OUTPUT ARGUMENT LIST:
2544	! FPVSX - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
2545	!
2546	! ATTRIBUTES:
2547	! LANGUAGE: FORTRAN 90
2548	!$$$
2549	IMPLICIT NONE
2550	!-----------------------------------------------------------------------
2551	real, parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
2552	, CLIQ=4.1855E+3,CVAP= 1.8460E+3 &
2553	, CICE=2.1060E+3,HSUB=2.8340E+6
2554	!
2555	real, parameter :: PSATK=PSAT*1.E-3
2556	real, parameter :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
2557	real, parameter :: DLDTI=CVAP-CICE &
2558	, XAI=-DLDTI/RV,XBI=XAI+HSUB/(RV*TTP)
2559	real T,TR
2560	!-----------------------------------------------------------------------
2561	TR=TTP/T
2562	!
2563	IF(T.GE.TTP)THEN
2564	FPVSX=PSATK(TRXA)EXP(XB*(1.-TR))
2565	ELSE
2566	FPVSX=PSATK(TRXAI)EXP(XBI*(1.-TR))
2567	ENDIF
2568	!
2569	END FUNCTION FPVSX
2570	!-----------------------------------------------------------------------
2571	!-----------------------------------------------------------------------
2572	REAL FUNCTION FPVSX0(T)
2573	!-----------------------------------------------------------------------
2574	IMPLICIT NONE
2575	real,parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
2576	, CLIQ=4.1855E+3,CVAP=1.8460E+3 &
2577	, CICE=2.1060E+3,HSUB=2.8340E+6
2578	real,PARAMETER :: PSATK=PSAT*1.E-3
2579	real,PARAMETER :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
2580	real,PARAMETER :: DLDTI=CVAP-CICE &
2581	, XAI=-DLDT/RV,XBI=XA+HSUB/(RV*TTP)
2582	real :: T,TR
2583	!-----------------------------------------------------------------------
2584	TR=TTP/T
2585	FPVSX0=PSATK(TRXA)EXP(XB*(1.-TR))
2586	!
2587	END FUNCTION FPVSX0
2588	!
2589	END MODULE module_mp_etanew

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

Download in other formats:

Original Format