Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_HWRF.F @ 1372

Last change on this file since 1372 was 1, checked in by lfita, 11 years ago

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

WRF: version v3.3
LMDZ: version v1818

More details in:

Property svn:executable set to *

File size: 110.9 KB

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

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

Download in other formats:

Original Format