| blob | bfb5a1e6d536daadd6bfb3e53b8ca467f8b232e1 |
1 !WRF:MODEL_MP:PHYSICS
2 !
3 !-- "Modified" fer_hires microphysics - 11 July 2016 version
4 !
5 !-- Updates based on NAM changes in 2011:
6 ! (a) Expanded rain lookup tables from 0.45 mm to 1 mm mean diameter.
7 ! (b) Allow cloud ice to fall (fall speeds based on 50 micron mean diameters).
8 ! (c) Cloud water autoconversion to rain follows Liu et al. (JAS, 2006)
9 ! (d) Fix to MY_GROWTH by multiplying estimates by 1.e-3
10 ! (e) Added integer function GET_INDEXR
11 ! (f) Added warning messages when unusual conditions occur, screened for
12 ! 5 different types of problems, such as (1) condensate in the
13 ! stratosphere, (2) temperature = NaN, (3) water supersaturation at
14 ! <180K, (4) too many iterations (>10) in the condensation function,
15 ! and (5) too many iterations (>10) in the deposition function.
16 !
17 !-- Updates based on jan19 2014 changes in the NMMB:
18 !
19 ! (1) Ice nucleation: Fletcher (1962) replaces Meyers et al. (1992)
20 ! (2) Cloud ice is a simple function of the number concentration from (1), and it
21 ! is no longer a fractional function of the large ice. Thus, the FLARGE &
22 ! FSMALL parameters are no longer used.
23 ! (3) T_ICE_init=-12 deg C provides a slight delay in the initial onset of ice.
24 ! (4) NLImax is a function of rime factor (RF) and temperature.
25 ! a) For RF>10, NLImax=1.e3. Mean ice diameters can exceed the 1 mm maximum
26 ! size in the tables so that NLICE=NLImax=1.e3.
27 ! b) Otherwise, NLImax is 10 L-1 at 0C and decreasing to 5 L-1 at <=-40C.
28 ! NLICE>NLImax at the maximum ice diameter of 1 mm.
29 ! (5) Can turn off ice processes by setting T_ICE & T_ICE_init to be < -100 deg C
30 ! (6) Modified the homogeneous freezing of cloud water when T<T_ICE.
31 ! (7) Reduce the fall speeds of rimed ice by making VEL_INC a function of
32 ! VrimeF**1.5 and not VrimeF**2.
33 ! (8) RHgrd=0.98 (98% RH) for the onset of condensation, to match what's been
34 ! tested for many months in the NAMX. Made obsolete with change in (13).
35 ! (9) Rime factor is *never* adjusted when NLICE>NLImax.
36 ! (10) Ice deposition does not change the rime factor (RF) when RF>=10 & T>T_ICE.
37 ! (11) Limit GAMMAS to <=1.5 (air resistance impact on ice fall speeds)
38 ! (12) NSImax is maximum # conc of ice crystals. At cold temperature NSImax is
39 ! calculated based on assuming 10% of total ice content is due to cloud ice.
40 !
41 !-- Further modifications starting on 23 July 2015
42 ! (13) RHgrd is passed in as an input argument so that it can vary for different
43 ! domains (RHgrd=0.98 for 12-km parent, 1.0 for 3-km nests)
44 ! (14) Use the old "PRAUT" cloud water autoconversion *threshold* (QAUT0)
46 !-- Further modifications starting on 28 July 2015
47 ! (15) Added calculations for radar reflectivity and number concentrations of
48 ! rain (Nrain) and precipitating ice (Nsnow).
49 ! (16) Removed double counting of air resistance term for riming onto ice (PIACW)
50 ! (17) The maximum rime factor (RFmx) is now a function of MASSI(INDEXS), accounting
51 ! for the increase in unrimed ice particle densities as values of INDEXS
52 ! decrease from the maximum upper limit of 1000 microns to the lower limit of
53 ! 50 microns, coinciding with the assumed size of cloud ice; see lines 1128-1134.
54 ! (18) A new closure is used for updating the rime factor, which is described in
55 ! detail near lines 1643-1682. The revised code is near lines 1683-1718.
56 ! (19) Restructured the two-pass algorithm to be more robust, removed the HAIL
57 ! & LARGE_RF logical variables so that NLICE>NLImax can occur.
58 ! (20) Increased nsimax (see !aug27 below)
59 ! (21) Modified the rain sedimentation (see two !aug27 blocks below)
60 ! (22) NInuclei is the lower of Fletcher (1962), Cooper (1986), or NSImax.
61 ! (23) NLImax is no longer used or enforced. Instead, INDEXS=MDImax when RF>20,
62 ! else INDEXS is a function of temperature. Look for !sep10 comment.
63 ! (24) An override was inserted for (18), such that the rime density is not diluted
64 ! diluted when RF>20. Look for !sep10 comment.
65 ! (25) Radar reflectivity calculations were changes to reduce radar bright bands,
66 ! limit enhanced, mixed-phase reflectivity to RF>=20. Look for !sep10 comments.
67 ! (26) NLICE is not to exceed NSI_max (250 L^-1) when RF<20. Look for !sep16 comments.
68 ! Commented out! (28) Increase hail fall speeds using Thompson et al. (2008). Look for !sep22 comments.
69 ! (29) Modify NLImax, INDEXS for RF>=20. Look for !sep22 comments.
70 ! (30) Check on NSmICE, Vci based on whether FLIMASS<1. Look for !sep22a comments.
71 ! Revised in (34)! (31) Introduced RFlag logical, which if =T enforces a lower limit of drop sizes not
72 ! to go below INDEXRmin and N0r is adjusted. Look for !nov25 comments (corrections,
73 ! refinements to sep25 & nov18 versions, includes an additional fix in nov25-fix).
74 ! Also set INDEXRmin=500 rather than 250 microns.
75 !-----------------------------------------------------------------------------
76 !--- The following changes now refer to dates when those were made in 2016.
77 !-----------------------------------------------------------------------------
78 ! (32) Convective (RF>=20, Ng~10 L^-1, RHOg~500 kg m^-3), transition (RF=10, Ng~25 L^-1,
79 ! RHOg~300 kg m^-3), & stratiform (RF<2) profiles are blended based on RF. !mar08
80 ! (33) Fixed bug in Biggs' freezing, put back in collisional drop freezing. !mar03
81 ! (34) Changes in (31) are revised so that INDEXRmin at and below 0C level is
82 ! based on a rain rate equal to the snowfall rate above the 0C level. !mar03
83 ! (35) Increase radar reflectivity when RF>10 and RQSnew > 2.5 g m^-3. !mar12
84 ! (36) !mar10 combines all elements of (32)-(35) together.
85 ! (37) Bug fixes for the changes in (34) and the RFLAG variable !apr18
86 ! (38) Revised Schumann-Ludlam limit. !apr18
87 ! (39) Simplified PCOND (cloud cond/evap) calculation !apr21
88 ! (40) Slight change in calculating RF. !apr22
89 ! (41) Reduce RF values for calculating mean sizes of snow, graupel, sleet/hail !apr22a
90 ! (42) Increase reflectivity from large, wet, high rime factor ice (graupel) by
91 ! assuming |Kw|**2/|Ki|**2 = 0.224 (Smith, 1984, JCAM).
92 ! (43) Major restructuring of code to allow N0r to vary from N0r0 !may11
93 ! (44) More major restructuring of code to use fixed XLS, XLV, XLF !may12
94 ! (45) Increased VEL_INC ~ VrimeF**2, put the enhanced graupel/hail fall speeds
95 ! from Thompson into the code but only in limited circumstances, restructured
96 ! and streamlined the INDEXS calculation, removed the upper limit for
97 ! for the vapor mixing ratio is at water saturation when calculating ice
98 ! deposition, and N0r is gradually increased for conditions supporting
99 ! drizzle when rain contents decrease below 0.25 g/m**3. !may17
100 ! (46) The may11 code changes that increase N0r0 when rain contents exceed 1 g m^-3
101 ! have been removed, limit the number of iterations calculating final rain
102 ! parameters, remove the revised N0r calculation for reflectivity. All of
103 ! the changes following those made in the may10 code. !may20
104 ! (47) Reduce the assumed # concentration of hail/sleet when RF>10 from 5 L^-1 to
105 ! 1 L^-1, and also reduce it for graupel when RF>5 from 10 L^-1 to 5 L^-1.
106 ! This is being done to try and make greater use of the Thompson graupel/hail
107 ! fallspeeds by having INDEXS==MDImax.
108 ! (48) Increased NCW from 200e6 to 300e6 for a more delayed onset of drizzle,
109 ! simplified drizzle algorithm to reduce/eliminate N0r bulls eyes and to allow
110 ! for supercooled drizzle, and set limits for 8.e6 <= N0r (m^-4) <= 1.e9 !may31
111 ! (49) Further restructuring of code to better define STRAT, DRZL logicals,
112 ! add these rain flags to mprates arrays !jun01
113 ! (50) Increase in reflectivity due to wet ice was commented out.
114 ! (51) Fixed minor bug to update INDEXR2 in the "rain_pass: do" loop. !jun13
115 ! (52) Final changes to Nsnow for boosting reflectivities from ice for
116 ! mass contents exceeding 5 g m^-3. !jun16
117 ! (53) Cosmetic changes only that do not affect the calculations. Removed old, unused
118 ! diagnostic arrays. Updated comments.
119 !
120 !-----------------------------------------------------------------------------
121 !
123 MODULE module_mp_fer_hires
125 !-----------------------------------------------------------------------
126 !-- The following changes were made on 24 July 2006.
127 ! (1) All known version 2.1 dependencies were removed from the
128 ! operational WRF NMM model code (search for "!HWRF")
129 ! (2) Incorporated code changes from the GFDL model (search for "!GFDL")
130 !-----------------------------------------------------------------------
131 !
132 REAL,PRIVATE,SAVE :: ABFR, CBFR, CIACW, CIACR, C_N0r0, &
133 & C_NR, CRAIN, &
134 & ARAUT, BRAUT, CN0r0, CN0r_DMRmin, CN0r_DMRmax, CRACW, ESW0, &
135 & RFmax, RQR_DRmin, RQR_DRmax, RR_DRmin, RR_DR1, RR_DR2, &
136 & RR_DR3, RR_DR4, RR_DR5, RR_DRmax, BETA6, PI_E, RFmx1, ARcw, & !jul31
137 & RH_NgC, RH_NgT, RQhail, AVhail, BVhail, QAUT0 !mar08 !may17
138 !
139 INTEGER,PRIVATE,PARAMETER :: INDEXRstrmax=500 !mar03, stratiform maximum
140 INTEGER, PRIVATE,PARAMETER :: MY_T1=1, MY_T2=35
141 REAL,PRIVATE,DIMENSION(MY_T1:MY_T2),SAVE :: MY_GROWTH
142 !
143 REAL, PRIVATE,PARAMETER :: DMImin=.05e-3, DMImax=1.e-3, &
144 & DelDMI=1.e-6,XMImin=1.e6*DMImin
145 REAL, PUBLIC,PARAMETER :: XMImax=1.e6*DMImax, XMIexp=.0536
146 INTEGER, PUBLIC,PARAMETER :: MDImin=XMImin, MDImax=XMImax
147 REAL, PRIVATE,DIMENSION(MDImin:MDImax) :: &
148 & ACCRI,VSNOWI,VENTI1,VENTI2
149 REAL, PUBLIC,DIMENSION(MDImin:MDImax) :: SDENS !-- For RRTM
150 !
151 REAL, PRIVATE,PARAMETER :: DMRmin=.05e-3, DMRmax=1.0e-3, &
152 & DelDMR=1.e-6, XMRmin=1.e6*DMRmin, XMRmax=1.e6*DMRmax
153 INTEGER, PUBLIC,PARAMETER :: MDRmin=XMRmin, MDRmax=XMRmax
154 !
155 REAL, PRIVATE,DIMENSION(MDRmin:MDRmax):: &
156 & ACCRR,MASSR,RRATE,VRAIN,VENTR1,VENTR2
157 !
158 INTEGER, PRIVATE,PARAMETER :: Nrime=40
159 REAL, DIMENSION(2:9,0:Nrime),PRIVATE,SAVE :: VEL_RF
160 !
161 INTEGER,PARAMETER :: NX=7501
162 REAL, PARAMETER :: XMIN=180.0,XMAX=330.0
163 REAL, DIMENSION(NX),PRIVATE,SAVE :: TBPVS,TBPVS0
164 REAL, PRIVATE,SAVE :: C1XPVS0,C2XPVS0,C1XPVS,C2XPVS
165 !
166 REAL, PRIVATE,PARAMETER :: &
167 !--- Physical constants follow:
168 & CP=1004.6, EPSQ=1.E-12, GRAV=9.806, RHOL=1000., RD=287.04 &
169 & ,RV=461.5, T0C=273.15, XLV=2.5E6, XLF=3.3358e+5 &
170 & ,pi=3.141592653589793 &
171 !--- Derived physical constants follow:
172 & ,EPS=RD/RV, EPS1=RV/RD-1., EPSQ1=1.001*EPSQ &
173 & ,RCP=1./CP, RCPRV=RCP/RV, RGRAV=1./GRAV, RRHOL=1./RHOL &
174 & ,XLS=XLV+XLF, XLV1=XLV/CP, XLF1=XLF/CP, XLS1=XLS/CP &
175 & ,XLV2=XLV*XLV/RV, XLS2=XLS*XLS/RV &
176 & ,XLV3=XLV*XLV*RCPRV, XLS3=XLS*XLS*RCPRV &
177 !--- Constants specific to the parameterization follow:
178 !--- CLIMIT/CLIMIT1 are lower limits for treating accumulated precipitation
179 & ,CLIMIT=10.*EPSQ, CLIMIT1=-CLIMIT &
180 & ,C1=1./3. &
181 & ,DMR1=.1E-3, DMR2=.2E-3, DMR3=.32E-3, DMR4=0.45E-3 &
182 & ,DMR5=0.67E-3 &
183 & ,XMR1=1.e6*DMR1, XMR2=1.e6*DMR2, XMR3=1.e6*DMR3 &
184 & ,XMR4=1.e6*DMR4, XMR5=1.e6*DMR5, RQRmix=0.05E-3, RQSmix=1.E-3 & !jul28 !apr27
185 & ,Cdry=1.634e13, Cwet=1./.224 !jul28 !apr27
186 INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3, MDR4=XMR4 &
187 & , MDR5=XMR5
189 !-- Debug 20120111
190 LOGICAL, SAVE :: WARN1=.TRUE.,WARN2=.TRUE.,WARN3=.TRUE.,WARN5=.TRUE.
191 REAL, SAVE :: Pwarn=75.E2, QTwarn=1.E-3
192 INTEGER, PARAMETER :: MAX_ITERATIONS=10
194 !
195 ! ======================================================================
196 !--- Important tunable parameters that are exported to other modules
197 !GFDL * RHgrd - generic reference to the threshold relative humidity for
198 !GFDL the onset of condensation
199 !GFDL (new) * RHgrd_in - "RHgrd" for the inner domain
200 !GFDL (new) * RHgrd_out - "RHgrd" for the outer domain
201 !HWRF 6/11/2010 mod - use lower RHgrd_out for p < 850 hPa
202 ! * T_ICE - temperature (C) threshold at which all remaining liquid water
203 ! is glaciated to ice
204 ! * T_ICE_init - maximum temperature (C) at which ice nucleation occurs
205 !
206 !-- To turn off ice processes, set T_ICE & T_ICE_init to <= -100. (i.e., -100 C)
207 !
208 ! * NLImax - maximum number concentrations (m**-3) of large ice (snow/graupel/sleet)
209 ! * NLImin - minimum number concentrations (m**-3) of large ice (snow/graupel/sleet)
210 ! * NSImax - maximum number concentrations (m**-3) of small ice crystals
211 ! * N0r0 - assumed intercept (m**-4) of rain drops if drop diameters are between 0.2 and 1.0 mm
212 ! * N0rmin - minimum intercept (m**-4) for rain drops
213 ! * NCW - number concentrations of cloud droplets (m**-3)
214 ! * PRINT_diag - for extended model diagnostics (code currently commented out)
215 ! ======================================================================
216 REAL, PUBLIC,PARAMETER :: &
217 & RHgrd_in=1. & !GFDL
218 & ,RHgrd_out=0.975 & !GFDL
219 & ,P_RHgrd_out=850.E2 & !HWRF 6/11/2010
220 & ,T_ICE=-40. &
221 & ,T_ICEK=T0C+T_ICE &
222 & ,T_ICE_init=-12. &
223 & ,NSI_max=250.E3 &
224 & ,NLImin=1.E3 &
225 & ,N0r0=8.E6 &
226 & ,N0rmin=1.E4 &
227 !!2-09-2012 & ,NCW=60.E6 & !GFDL
228 !! based on Aligo's email,NCW is changed to 250E6
229 & ,NCW=250.E6 !GFDL
230 !HWRF & ,NCW=100.E6 &
231 LOGICAL, PARAMETER :: PRINT_diag=.FALSE. !GFDL
232 !--- Other public variables passed to other routines:
233 REAL, PUBLIC,DIMENSION(MDImin:MDImax) :: MASSI
234 !
235 !
236 CONTAINS
238 !-----------------------------------------------------------------------
239 !-----------------------------------------------------------------------
240 SUBROUTINE FER_HIRES (itimestep,DT,DX,DY,GID,RAINNC,RAINNCV, & !GID
241 & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qt, & !gopal's doing
242 & LOWLYR,SR, &
243 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
244 & QC,QR,QI, &
245 & ids,ide, jds,jde, kds,kde, &
246 & ims,ime, jms,jme, kms,kme, &
247 & its,ite, jts,jte, kts,kte )
248 !HWRF SUBROUTINE ETAMP_NEW (itimestep,DT,DX,DY, &
249 !HWRF & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qc, &
250 !HWRF & LOWLYR,SR, &
251 !HWRF & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
252 !HWRF & mp_restart_state,tbpvs_state,tbpvs0_state, &
253 !HWRF & RAINNC,RAINNCV, &
254 !HWRF & ids,ide, jds,jde, kds,kde, &
255 !HWRF & ims,ime, jms,jme, kms,kme, &
256 !HWRF & its,ite, jts,jte, kts,kte )
257 !-----------------------------------------------------------------------
258 IMPLICIT NONE
259 !-----------------------------------------------------------------------
261 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
262 & ,IMS,IME,JMS,JME,KMS,KME &
263 & ,ITS,ITE,JTS,JTE,KTS,KTE &
264 & ,ITIMESTEP,GID ! GID gopal's doing
266 REAL, INTENT(IN) :: DT,DX,DY
267 REAL, INTENT(IN), DIMENSION(ims:ime, kms:kme, jms:jme):: &
268 & dz8w,p_phy,pi_phy,rho_phy
269 REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme):: &
270 & th_phy,qv,qt,qc,qr,qi
271 REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
272 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
273 REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: &
274 & RAINNC,RAINNCV
275 REAL, INTENT(OUT), DIMENSION(ims:ime,jms:jme):: SR
276 !
277 !HWRF REAL,DIMENSION(*),INTENT(INOUT) :: MP_RESTART_STATE
278 !
279 !HWRF REAL,DIMENSION(nx),INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
280 !
281 INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
283 !-----------------------------------------------------------------------
284 ! LOCAL VARS
285 !-----------------------------------------------------------------------
287 ! SOME VARS WILL BE USED FOR DATA ASSIMILATION (DON'T NEED THEM NOW).
288 ! THEY ARE TREATED AS LOCAL VARS, BUT WILL BECOME STATE VARS IN THE
289 ! FUTURE. SO, WE DECLARED THEM AS MEMORY SIZES FOR THE FUTURE USE
291 ! TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related
292 ! the microphysics scheme. Instead, they will be used by Eta precip
293 ! assimilation.
295 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
296 & TLATGS_PHY,TRAIN_PHY
297 REAL, DIMENSION(ims:ime,jms:jme):: APREC,PREC,ACPREC
298 REAL, DIMENSION(its:ite, kts:kte, jts:jte):: t_phy
300 INTEGER :: I,J,K,KFLIP
301 REAL :: WC
302 !
303 !-----------------------------------------------------------------------
304 !**********************************************************************
305 !-----------------------------------------------------------------------
306 !
307 !HWRF MY_GROWTH(MY_T1:MY_T2)=MP_RESTART_STATE(MY_T1:MY_T2)
308 !HWRF!
309 !HWRF C1XPVS0=MP_RESTART_STATE(MY_T2+1)
310 !HWRF C2XPVS0=MP_RESTART_STATE(MY_T2+2)
311 !HWRF C1XPVS =MP_RESTART_STATE(MY_T2+3)
312 !HWRF C2XPVS =MP_RESTART_STATE(MY_T2+4)
313 !HWRF CIACW =MP_RESTART_STATE(MY_T2+5)
314 !HWRF CIACR =MP_RESTART_STATE(MY_T2+6)
315 !HWRF CRACW =MP_RESTART_STATE(MY_T2+7)
316 !HWRF CRAUT =MP_RESTART_STATE(MY_T2+8)
317 !HWRF!
318 !HWRF TBPVS(1:NX) =TBPVS_STATE(1:NX)
319 !HWRF TBPVS0(1:NX)=TBPVS0_STATE(1:NX)
320 !
321 !----------
322 !2015-03-30, recalculate some constants which may depend on phy time step
323 CALL MY_GROWTH_RATES (DT)
325 !--- CIACW is used in calculating riming rates
326 ! The assumed effective collection efficiency of cloud water rimed onto
327 ! ice is =0.5 below:
328 !
329 CIACW=DT*0.25*PI_E*0.5*(1.E5)**C1
330 !
331 !--- CIACR is used in calculating freezing of rain colliding with large ice
332 ! The assumed collection efficiency is 1.0
333 !
334 CIACR=PI_E*DT
335 !
336 !--- CRACW is used in calculating collection of cloud water by rain (an
337 ! assumed collection efficiency of 1.0)
338 !
339 CRACW=DT*0.25*PI_E*1.0
340 !
341 !-- See comments in subroutine etanewhr_init starting with variable RDIS=
342 !
343 BRAUT=DT*1.1E10*BETA6/NCW
345 ! write(*,*)'dt=',dt
346 ! write(*,*)'pi_e=',pi_e
347 ! write(*,*)'ciacw=',ciacw
348 ! write(*,*)'ciacr=',ciacr
349 ! write(*,*)'cracw=',cracw
350 ! write(*,*)'araut=',araut
351 ! write(*,*)'braut=',braut
352 !! END OF adding, 2015-03-30
353 !-----------
355 DO j = jts,jte
356 DO k = kts,kte
357 DO i = its,ite
358 t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
359 qv(i,k,j)=qv(i,k,j)/(1.+qv(i,k,j)) !Convert to specific humidity
360 ENDDO
361 ENDDO
362 ENDDO
364 ! initial data assimilation vars (will need to delete this part in the future)
366 DO j = jts,jte
367 DO k = kts,kte
368 DO i = its,ite
369 TLATGS_PHY (i,k,j)=0.
370 TRAIN_PHY (i,k,j)=0.
371 ENDDO
372 ENDDO
373 ENDDO
375 DO j = jts,jte
376 DO i = its,ite
377 ACPREC(i,j)=0.
378 APREC (i,j)=0.
379 PREC (i,j)=0.
380 SR (i,j)=0.
381 ENDDO
382 ENDDO
384 !-- 6/11/2010: Update QT, F_ice, F_rain arrays
385 DO j = jts,jte
386 DO k = kts,kte
387 DO i = its,ite
388 QT(I,K,J)=QC(I,K,J)+QR(I,K,J)+QI(I,K,J)
389 IF (QI(I,K,J) <= EPSQ) THEN
390 F_ICE_PHY(I,K,J)=0.
391 F_RIMEF_PHY(I,K,J)=1.
392 IF (T_PHY(I,K,J) < T_ICEK) F_ICE_PHY(I,K,J)=1.
393 ELSE
394 F_ICE_PHY(I,K,J)=MAX( 0., MIN(1., QI(I,K,J)/QT(I,K,J) ) )
395 ENDIF
396 IF (QR(I,K,J) <= EPSQ) THEN
397 F_RAIN_PHY(I,K,J)=0.
398 ELSE
399 F_RAIN_PHY(I,K,J)=QR(I,K,J)/(QR(I,K,J)+QC(I,K,J))
400 ENDIF
401 ENDDO
402 ENDDO
403 ENDDO
405 !-----------------------------------------------------------------------
407 CALL EGCP01DRV(GID,DT,LOWLYR, &
408 & APREC,PREC,ACPREC,SR, &
409 & dz8w,rho_phy,qt,t_phy,qv,F_ICE_PHY,P_PHY, &
410 & F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
411 & ids,ide, jds,jde, kds,kde, &
412 & ims,ime, jms,jme, kms,kme, &
413 & its,ite, jts,jte, kts,kte )
414 !-----------------------------------------------------------------------
416 DO j = jts,jte
417 DO k = kts,kte
418 DO i = its,ite
419 th_phy(i,k,j) = t_phy(i,k,j)/pi_phy(i,k,j)
420 qv(i,k,j)=qv(i,k,j)/(1.-qv(i,k,j)) !Convert to mixing ratio
421 WC=qt(I,K,J)
422 QI(I,K,J)=0.
423 QR(I,K,J)=0.
424 QC(I,K,J)=0.
425 IF(F_ICE_PHY(I,K,J)>=1.)THEN
426 QI(I,K,J)=WC
427 ELSEIF(F_ICE_PHY(I,K,J)<=0.)THEN
428 QC(I,K,J)=WC
429 ELSE
430 QI(I,K,J)=F_ICE_PHY(I,K,J)*WC
431 QC(I,K,J)=WC-QI(I,K,J)
432 ENDIF
433 !
434 IF(QC(I,K,J)>0..AND.F_RAIN_PHY(I,K,J)>0.)THEN
435 IF(F_RAIN_PHY(I,K,J).GE.1.)THEN
436 QR(I,K,J)=QC(I,K,J)
437 QC(I,K,J)=0.
438 ELSE
439 QR(I,K,J)=F_RAIN_PHY(I,K,J)*QC(I,K,J)
440 QC(I,K,J)=QC(I,K,J)-QR(I,K,J)
441 ENDIF
442 endif
443 ENDDO
444 ENDDO
445 ENDDO
446 !
447 ! update rain (from m to mm)
449 DO j=jts,jte
450 DO i=its,ite
451 RAINNC(i,j)=APREC(i,j)*1000.+RAINNC(i,j)
452 RAINNCV(i,j)=APREC(i,j)*1000.
453 ENDDO
454 ENDDO
455 !
456 !HWRF MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
457 !HWRF MP_RESTART_STATE(MY_T2+1)=C1XPVS0
458 !HWRF MP_RESTART_STATE(MY_T2+2)=C2XPVS0
459 !HWRF MP_RESTART_STATE(MY_T2+3)=C1XPVS
460 !HWRF MP_RESTART_STATE(MY_T2+4)=C2XPVS
461 !HWRF MP_RESTART_STATE(MY_T2+5)=CIACW
462 !HWRF MP_RESTART_STATE(MY_T2+6)=CIACR
463 !HWRF MP_RESTART_STATE(MY_T2+7)=CRACW
464 !HWRF MP_RESTART_STATE(MY_T2+8)=CRAUT
465 !HWRF!
466 !HWRF TBPVS_STATE(1:NX) =TBPVS(1:NX)
467 !HWRF TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
469 !-----------------------------------------------------------------------
471 END SUBROUTINE FER_HIRES
473 !-----------------------------------------------------------------------
474 !-----------------------------------------------------------------------
475 ! NOTE: The only differences between FER_HIRES and FER_HIRES_ADVECT
476 ! is that the QT, and F_* are all local variables in the advected
477 ! version, and QRIMEF is only in the advected version. The innards
478 ! are all the same.
479 SUBROUTINE FER_HIRES_ADVECT (itimestep,DT,DX,DY,GID,RAINNC,RAINNCV, & !GID
480 & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv, & !gopal's doing
481 & LOWLYR,SR, &
482 & QC,QR,QI,QRIMEF, &
483 & ids,ide, jds,jde, kds,kde, &
484 & ims,ime, jms,jme, kms,kme, &
485 & its,ite, jts,jte, kts,kte )
486 !HWRF SUBROUTINE ETAMP_NEW (itimestep,DT,DX,DY, &
487 !HWRF & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qc, &
488 !HWRF & LOWLYR,SR, &
489 !HWRF & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
490 !HWRF & mp_restart_state,tbpvs_state,tbpvs0_state, &
491 !HWRF & RAINNC,RAINNCV, &
492 !HWRF & ids,ide, jds,jde, kds,kde, &
493 !HWRF & ims,ime, jms,jme, kms,kme, &
494 !HWRF & its,ite, jts,jte, kts,kte )
495 !-----------------------------------------------------------------------
496 IMPLICIT NONE
497 !-----------------------------------------------------------------------
499 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
500 & ,IMS,IME,JMS,JME,KMS,KME &
501 & ,ITS,ITE,JTS,JTE,KTS,KTE &
502 & ,ITIMESTEP,GID ! GID gopal's doing
504 REAL, INTENT(IN) :: DT,DX,DY
505 REAL, INTENT(IN), DIMENSION(ims:ime, kms:kme, jms:jme):: &
506 & dz8w,p_phy,pi_phy,rho_phy
507 REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme):: &
508 & th_phy,qv,qc,qr,qi,qrimef
509 REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: &
510 & RAINNC,RAINNCV
511 REAL, INTENT(OUT), DIMENSION(ims:ime,jms:jme):: SR
512 !
513 !HWRF REAL,DIMENSION(*),INTENT(INOUT) :: MP_RESTART_STATE
514 !
515 !HWRF REAL,DIMENSION(nx),INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
516 !
517 INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
519 !-----------------------------------------------------------------------
520 ! LOCAL VARS
521 !-----------------------------------------------------------------------
523 REAL, DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
524 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, QT
526 ! SOME VARS WILL BE USED FOR DATA ASSIMILATION (DON'T NEED THEM NOW).
527 ! THEY ARE TREATED AS LOCAL VARS, BUT WILL BECOME STATE VARS IN THE
528 ! FUTURE. SO, WE DECLARED THEM AS MEMORY SIZES FOR THE FUTURE USE
529 ! TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related
530 ! the microphysics scheme. Instead, they will be used by Eta precip
531 ! assimilation.
533 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
534 & TLATGS_PHY,TRAIN_PHY
535 REAL, DIMENSION(ims:ime,jms:jme):: APREC,PREC,ACPREC
536 REAL, DIMENSION(its:ite, kts:kte, jts:jte):: t_phy
538 INTEGER :: I,J,K,KFLIP
539 REAL :: WC
540 !
541 !-----------------------------------------------------------------------
542 !**********************************************************************
543 !-----------------------------------------------------------------------
544 !
545 !HWRF MY_GROWTH(MY_T1:MY_T2)=MP_RESTART_STATE(MY_T1:MY_T2)
546 !HWRF!
547 !HWRF C1XPVS0=MP_RESTART_STATE(MY_T2+1)
548 !HWRF C2XPVS0=MP_RESTART_STATE(MY_T2+2)
549 !HWRF C1XPVS =MP_RESTART_STATE(MY_T2+3)
550 !HWRF C2XPVS =MP_RESTART_STATE(MY_T2+4)
551 !HWRF CIACW =MP_RESTART_STATE(MY_T2+5)
552 !HWRF CIACR =MP_RESTART_STATE(MY_T2+6)
553 !HWRF CRACW =MP_RESTART_STATE(MY_T2+7)
554 !HWRF CRAUT =MP_RESTART_STATE(MY_T2+8)
555 !HWRF!
556 !HWRF TBPVS(1:NX) =TBPVS_STATE(1:NX)
557 !HWRF TBPVS0(1:NX)=TBPVS0_STATE(1:NX)
558 !
559 !----------
560 !2015-03-30, recalculate some constants which may depend on phy time step
561 CALL MY_GROWTH_RATES (DT)
563 !--- CIACW is used in calculating riming rates
564 ! The assumed effective collection efficiency of cloud water rimed onto
565 ! ice is =0.5 below:
566 !
567 CIACW=DT*0.25*PI_E*0.5*(1.E5)**C1
568 !
569 !--- CIACR is used in calculating freezing of rain colliding with large ice
570 ! The assumed collection efficiency is 1.0
571 !
572 CIACR=PI_E*DT
573 !
574 !--- CRACW is used in calculating collection of cloud water by rain (an
575 ! assumed collection efficiency of 1.0)
576 !
577 CRACW=DT*0.25*PI_E*1.0
578 !
579 !-- See comments in subroutine etanewhr_init starting with variable RDIS=
580 !
581 BRAUT=DT*1.1E10*BETA6/NCW
583 ! write(*,*)'dt=',dt
584 ! write(*,*)'pi_e=',pi_e
585 ! write(*,*)'ciacw=',ciacw
586 ! write(*,*)'ciacr=',ciacr
587 ! write(*,*)'cracw=',cracw
588 ! write(*,*)'araut=',araut
589 ! write(*,*)'braut=',braut
590 !! END OF adding, 2015-03-30
591 !-----------
593 DO j = jts,jte
594 DO k = kts,kte
595 DO i = its,ite
596 t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
597 qv(i,k,j)=qv(i,k,j)/(1.+qv(i,k,j)) !Convert to specific humidity
598 ENDDO
599 ENDDO
600 ENDDO
602 ! initial data assimilation vars (will need to delete this part in the future)
604 DO j = jts,jte
605 DO k = kts,kte
606 DO i = its,ite
607 TLATGS_PHY (i,k,j)=0.
608 TRAIN_PHY (i,k,j)=0.
609 ENDDO
610 ENDDO
611 ENDDO
613 DO j = jts,jte
614 DO i = its,ite
615 ACPREC(i,j)=0.
616 APREC (i,j)=0.
617 PREC (i,j)=0.
618 SR (i,j)=0.
619 ENDDO
620 ENDDO
622 !-- 6/11/2010: Update QT, F_ice, F_rain arrays
624 DO j = jts,jte
625 DO k = kts,kte
626 DO i = its,ite
627 QT(I,K,J)=QC(I,K,J)+QR(I,K,J)+QI(I,K,J)
628 IF (QI(I,K,J) <= EPSQ) THEN
629 F_ICE_PHY(I,K,J)=0.
630 F_RIMEF_PHY(I,K,J)=1.
631 IF (T_PHY(I,K,J) < T_ICEK) F_ICE_PHY(I,K,J)=1.
632 ELSE
633 F_ICE_PHY(I,K,J)=MAX( 0., MIN(1., QI(I,K,J)/QT(I,K,J) ) )
634 F_RIMEF_PHY(I,K,J)=QRIMEF(I,K,J)/QI(I,K,J)
635 ENDIF
636 IF (QR(I,K,J) <= EPSQ) THEN
637 F_RAIN_PHY(I,K,J)=0.
638 ELSE
639 F_RAIN_PHY(I,K,J)=QR(I,K,J)/(QR(I,K,J)+QC(I,K,J))
640 ENDIF
641 ENDDO
642 ENDDO
643 ENDDO
645 !-----------------------------------------------------------------------
647 CALL EGCP01DRV(GID,DT,LOWLYR, &
648 & APREC,PREC,ACPREC,SR, &
649 & dz8w,rho_phy,qt,t_phy,qv,F_ICE_PHY,P_PHY, &
650 & F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
651 & ids,ide, jds,jde, kds,kde, &
652 & ims,ime, jms,jme, kms,kme, &
653 & its,ite, jts,jte, kts,kte )
654 !-----------------------------------------------------------------------
656 DO j = jts,jte
657 DO k = kts,kte
658 DO i = its,ite
659 th_phy(i,k,j) = t_phy(i,k,j)/pi_phy(i,k,j)
660 qv(i,k,j)=qv(i,k,j)/(1.-qv(i,k,j)) !Convert to mixing ratio
661 WC=qt(I,K,J)
662 QI(I,K,J)=0.
663 QR(I,K,J)=0.
664 QC(I,K,J)=0.
665 IF(F_ICE_PHY(I,K,J)>=1.)THEN
666 QI(I,K,J)=WC
667 ELSEIF(F_ICE_PHY(I,K,J)<=0.)THEN
668 QC(I,K,J)=WC
669 ELSE
670 QI(I,K,J)=F_ICE_PHY(I,K,J)*WC
671 QC(I,K,J)=WC-QI(I,K,J)
672 ENDIF
673 !
674 IF(QC(I,K,J)>0..AND.F_RAIN_PHY(I,K,J)>0.)THEN
675 IF(F_RAIN_PHY(I,K,J).GE.1.)THEN
676 QR(I,K,J)=QC(I,K,J)
677 QC(I,K,J)=0.
678 ELSE
679 QR(I,K,J)=F_RAIN_PHY(I,K,J)*QC(I,K,J)
680 QC(I,K,J)=QC(I,K,J)-QR(I,K,J)
681 ENDIF
682 endif
683 QRIMEF(I,K,J)=QI(I,K,J)*F_RIMEF_PHY(I,K,J)
684 ENDDO
685 ENDDO
686 ENDDO
687 !
688 ! update rain (from m to mm)
690 DO j=jts,jte
691 DO i=its,ite
692 RAINNC(i,j)=APREC(i,j)*1000.+RAINNC(i,j)
693 RAINNCV(i,j)=APREC(i,j)*1000.
694 ENDDO
695 ENDDO
696 !
697 !HWRF MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
698 !HWRF MP_RESTART_STATE(MY_T2+1)=C1XPVS0
699 !HWRF MP_RESTART_STATE(MY_T2+2)=C2XPVS0
700 !HWRF MP_RESTART_STATE(MY_T2+3)=C1XPVS
701 !HWRF MP_RESTART_STATE(MY_T2+4)=C2XPVS
702 !HWRF MP_RESTART_STATE(MY_T2+5)=CIACW
703 !HWRF MP_RESTART_STATE(MY_T2+6)=CIACR
704 !HWRF MP_RESTART_STATE(MY_T2+7)=CRACW
705 !HWRF MP_RESTART_STATE(MY_T2+8)=CRAUT
706 !HWRF!
707 !HWRF TBPVS_STATE(1:NX) =TBPVS(1:NX)
708 !HWRF TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
710 !-----------------------------------------------------------------------
712 END SUBROUTINE FER_HIRES_ADVECT
714 !-----------------------------------------------------------------------
715 !-----------------------------------------------------------------------
717 SUBROUTINE EGCP01DRV(GID, & !GID gopal's doing
718 & DTPH,LOWLYR,APREC,PREC,ACPREC,SR, &
719 & dz8w,RHO_PHY,CWM_PHY,T_PHY,Q_PHY,F_ICE_PHY,P_PHY, &
720 & F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
721 & ids,ide, jds,jde, kds,kde, &
722 & ims,ime, jms,jme, kms,kme, &
723 & its,ite, jts,jte, kts,kte)
724 !-----------------------------------------------------------------------
725 ! DTPH Physics time step (s)
726 ! CWM_PHY (qt) Mixing ratio of the total condensate. kg/kg
727 ! Q_PHY Mixing ratio of water vapor. kg/kg
728 ! F_RAIN_PHY Fraction of rain.
729 ! F_ICE_PHY Fraction of ice.
730 ! F_RIMEF_PHY Mass ratio of rimed ice (rime factor).
731 !
732 !TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related the
733 !micrphysics sechme. Instead, they will be used by Eta precip assimilation.
734 !
735 !-----------------------------------------------------------------------
736 !--- Variables APREC,PREC,ACPREC,SR are calculated for precip assimilation
737 ! and/or ZHAO's scheme in Eta and are not required by this microphysics
738 ! scheme itself.
739 !-----------------------------------------------------------------------
740 IMPLICIT NONE
741 !-----------------------------------------------------------------------
742 !
743 ! VARIABLES PASSED IN/OUT
744 INTEGER,INTENT(IN ) :: ids,ide, jds,jde, kds,kde &
745 & ,ims,ime, jms,jme, kms,kme &
746 & ,its,ite, jts,jte, kts,kte
747 INTEGER,INTENT(IN ) :: GID ! grid%id gopal's doing
748 REAL,INTENT(IN) :: DTPH
749 INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
750 REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: &
751 & APREC,PREC,ACPREC,SR
752 REAL,DIMENSION( its:ite, kts:kte, jts:jte ),INTENT(INOUT) :: t_phy
753 REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: &
754 & dz8w,P_PHY,RHO_PHY
755 REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(INOUT) :: &
756 & CWM_PHY, F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY &
757 & ,Q_PHY,TRAIN_PHY
758 !
759 !-----------------------------------------------------------------------
760 !LOCAL VARIABLES
761 !-----------------------------------------------------------------------
762 !
763 !HWRF - Below are directives in the operational code that have been removed,
764 ! where "TEMP_DEX" has been replaced with "I,J,L" and "TEMP_DIMS" has
765 ! been replaced with "its:ite,jts:jte,kts:kte"
766 !HWRF#define CACHE_FRIENDLY_MP_ETANEW
767 !HWRF#ifdef CACHE_FRIENDLY_MP_ETANEW
768 !HWRF# define TEMP_DIMS kts:kte,its:ite,jts:jte
769 !HWRF# define TEMP_DEX L,I,J
770 !HWRF#else
771 !HWRF# define TEMP_DIMS its:ite,jts:jte,kts:kte
772 !HWRF# define TEMP_DEX I,J,L
773 !HWRF#endif
774 !HWRF!
775 INTEGER :: LSFC,I,J,I_index,J_index,L,K,KFLIP
776 !HWRF REAL,DIMENSION(TEMP_DIMS) :: CWM,T,Q,TRAIN,TLATGS,P
777 REAL,DIMENSION(its:ite,jts:jte,kts:kte) :: &
778 & CWM,T,Q,TRAIN,TLATGS,P
779 REAL,DIMENSION(kts:kte,its:ite,jts:jte) :: F_ice,F_rain,F_RimeF
780 INTEGER,DIMENSION(its:ite,jts:jte) :: LMH
781 REAL :: TC,WC,QI,QR,QW,Fice,Frain,DUM,ASNOW,ARAIN
782 REAL,DIMENSION(kts:kte) :: P_col,Q_col,T_col,QV_col,WC_col, &
783 & RimeF_col,QI_col,QR_col,QW_col, THICK_col, RHC_col, DPCOL, &
784 & pcond1d,pidep1d, &
785 & piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d,pimlt1d, &
786 & praut1d,pracw1d,prevp1d,pisub1d,pevap1d, DBZ_col,NR_col,NS_col, &
787 & vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d,INDEXR1d, & !jul28
788 & RFlag1d !jun01
789 REAL,DIMENSION(2) :: PRECtot,PRECmax
790 !-----------------------------------------------------------------------
791 !
792 DO J=JTS,JTE
793 DO I=ITS,ITE
794 LMH(I,J) = KTE-LOWLYR(I,J)+1
795 ENDDO
796 ENDDO
799 DO 98 J=JTS,JTE
800 DO 98 I=ITS,ITE
801 DO L=KTS,KTE
802 KFLIP=KTE+1-L
803 CWM(I,J,L)=CWM_PHY(I,KFLIP,J)
804 T(I,J,L)=T_PHY(I,KFLIP,J)
805 Q(I,J,L)=Q_PHY(I,KFLIP,J)
806 P(I,J,L)=P_PHY(I,KFLIP,J)
807 TLATGS(I,J,L)=TLATGS_PHY(I,KFLIP,J)
808 TRAIN(I,J,L)=TRAIN_PHY(I,KFLIP,J)
809 F_ice(L,I,J)=F_ice_PHY(I,KFLIP,J)
810 F_rain(L,I,J)=F_rain_PHY(I,KFLIP,J)
811 F_RimeF(L,I,J)=F_RimeF_PHY(I,KFLIP,J)
812 ENDDO
813 98 CONTINUE
815 DO 100 J=JTS,JTE
816 DO 100 I=ITS,ITE
817 LSFC=LMH(I,J) ! "L" of surface
818 !
819 DO K=KTS,KTE
820 KFLIP=KTE+1-K
821 DPCOL(K)=RHO_PHY(I,KFLIP,J)*GRAV*dz8w(I,KFLIP,J)
822 ENDDO
823 !
824 !
825 !--- Initialize column data (1D arrays)
826 !
827 IF (CWM(I,J,1) .LE. EPSQ) CWM(I,J,1)=EPSQ
828 F_ice(1,I,J)=1.
829 F_rain(1,I,J)=0.
830 F_RimeF(1,I,J)=1.
831 DO L=1,LSFC
832 !
833 !--- Pressure (Pa) = (Psfc-Ptop)*(ETA/ETA_sfc)+Ptop
834 !
835 P_col(L)=P(I,J,L)
836 !
837 !--- Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
838 !
839 THICK_col(L)=DPCOL(L)*RGRAV
840 T_col(L)=T(I,J,L)
841 TC=T_col(L)-T0C
842 QV_col(L)=max(EPSQ, Q(I,J,L))
843 IF (CWM(I,J,L) .LE. EPSQ1) THEN
844 WC_col(L)=0.
845 IF (TC .LT. T_ICE) THEN
846 F_ice(L,I,J)=1.
847 ELSE
848 F_ice(L,I,J)=0.
849 ENDIF
850 F_rain(L,I,J)=0.
851 F_RimeF(L,I,J)=1.
852 ELSE
853 WC_col(L)=CWM(I,J,L)
855 !-- Debug 20120111: TC==TC will fail if NaN
856 IF (WC_col(L)>QTwarn .AND. P_col(L)<Pwarn .AND. TC==TC) THEN
857 WRITE(0,*) 'WARN4: >1 g/kg condensate in stratosphere; I,J,L,TC,P,QT=', &
858 I,J,L,TC,.01*P_col(L),1000.*WC_col(L)
859 QTwarn=MAX(WC_col(L),10.*QTwarn)
860 Pwarn=MIN(P_col(L),0.5*Pwarn)
861 ENDIF
862 !-- TC/=TC will pass if TC is NaN
863 IF (WARN5 .AND. TC/=TC) THEN
864 WRITE(0,*) 'WARN5: NaN temperature; I,J,L,P=',I,J,L,.01*P_col(L)
865 WARN5=.FALSE.
866 ENDIF
868 ENDIF
869 IF (T_ICE<=-100.) F_ice(L,I,J)=0. !-- For no ice runs
870 !
871 !--- Determine composition of condensate in terms of
872 ! cloud water, ice, & rain
873 !
874 WC=WC_col(L)
875 QI=0.
876 QR=0.
877 QW=0.
878 Fice=F_ice(L,I,J)
879 Frain=F_rain(L,I,J)
880 IF (Fice .GE. 1.) THEN
881 QI=WC
882 ELSE IF (Fice .LE. 0.) THEN
883 QW=WC
884 ELSE
885 QI=Fice*WC
886 QW=WC-QI
887 ENDIF
888 IF (QW.GT.0. .AND. Frain.GT.0.) THEN
889 IF (Frain .GE. 1.) THEN
890 QR=QW
891 QW=0.
892 ELSE
893 QR=Frain*QW
894 QW=QW-QR
895 ENDIF
896 ENDIF
897 IF (QI .LE. 0.) F_RimeF(L,I,J)=1.
898 RimeF_col(L)=F_RimeF(L,I,J)
899 QI_col(L)=QI
900 QR_col(L)=QR
901 QW_col(L)=QW
902 !GFDL => New. Added RHC_col to allow for height- and grid-dependent values for
903 !GFDL the relative humidity threshold for condensation ("RHgrd")
904 !6/11/2010 mod - Use lower RHgrd_out threshold for < 850 hPa
905 !------------------------------------------------------------
906 IF(GID .EQ. 1 .AND. P_col(L)<P_RHgrd_out) THEN ! gopal's doing based on GFDL
907 RHC_col(L)=RHgrd_out
908 ELSE
909 RHC_col(L)=RHgrd_in
910 ENDIF
911 !------------------------------------------------------------
912 ENDDO
913 !
914 !#######################################################################
915 !
916 !--- Perform the microphysical calculations in this column
917 !
918 I_index=I
919 J_index=J
921 CALL EGCP01COLUMN ( ARAIN, ASNOW, DTPH, RHC_col, &
922 & I_index, J_index, LSFC, &
923 & P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
924 & THICK_col, WC_col, KTS,KTE, pcond1d,pidep1d, &
925 & piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d,pimlt1d, &
926 & praut1d,pracw1d,prevp1d,pisub1d,pevap1d, DBZ_col,NR_col,NS_col, &
927 & vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d,INDEXR1d, & !jul28
928 & RFlag1d) !jun01
930 !
931 !#######################################################################
932 !
933 !
934 !--- Update storage arrays
935 !
936 DO L=1,LSFC
937 TRAIN(I,J,L)=(T_col(L)-T(I,J,L))/DTPH
938 TLATGS(I,J,L)=T_col(L)-T(I,J,L)
939 T(I,J,L)=T_col(L)
940 Q(I,J,L)=QV_col(L)
941 CWM(I,J,L)=WC_col(L)
942 !
943 !--- REAL*4 array storage
944 !
945 IF (QI_col(L) .LE. EPSQ) THEN
946 F_ice(L,I,J)=0.
947 IF (T_col(L) .LT. T_ICEK) F_ice(L,I,J)=1.
948 F_RimeF(L,I,J)=1.
949 ELSE
950 F_ice(L,I,J)=MAX( 0., MIN(1., QI_col(L)/WC_col(L)) )
951 F_RimeF(L,I,J)=MAX(1., RimeF_col(L))
952 ENDIF
953 IF (QR_col(L) .LE. EPSQ) THEN
954 DUM=0
955 ELSE
956 DUM=QR_col(L)/(QR_col(L)+QW_col(L))
957 ENDIF
958 F_rain(L,I,J)=DUM
959 !
960 ENDDO
961 !
962 !--- Update accumulated precipitation statistics
963 !
964 !--- Surface precipitation statistics; SR is fraction of surface
965 ! precipitation (if >0) associated with snow
966 !
967 APREC(I,J)=(ARAIN+ASNOW)*RRHOL ! Accumulated surface precip (depth in m) !<--- Ying
968 PREC(I,J)=PREC(I,J)+APREC(I,J)
969 ACPREC(I,J)=ACPREC(I,J)+APREC(I,J)
970 IF(APREC(I,J) .LT. 1.E-8) THEN
971 SR(I,J)=0.
972 ELSE
973 SR(I,J)=RRHOL*ASNOW/APREC(I,J)
974 ENDIF
975 ! !
976 ! !--- Debug statistics
977 ! !
978 ! IF (PRINT_diag) THEN
979 ! PRECtot(1)=PRECtot(1)+ARAIN
980 ! PRECtot(2)=PRECtot(2)+ASNOW
981 ! PRECmax(1)=MAX(PRECmax(1), ARAIN)
982 ! PRECmax(2)=MAX(PRECmax(2), ASNOW)
983 ! ENDIF
986 !#######################################################################
987 !#######################################################################
988 !
989 100 CONTINUE ! End "I" & "J" loops
990 DO 101 J=JTS,JTE
991 DO 101 I=ITS,ITE
992 DO L=KTS,KTE
993 KFLIP=KTE+1-L
994 CWM_PHY(I,KFLIP,J)=CWM(I,J,L)
995 T_PHY(I,KFLIP,J)=T(I,J,L)
996 Q_PHY(I,KFLIP,J)=Q(I,J,L)
997 TLATGS_PHY(I,KFLIP,J)=TLATGS(I,J,L)
998 TRAIN_PHY(I,KFLIP,J)=TRAIN(I,J,L)
999 F_ice_PHY(I,KFLIP,J)=F_ice(L,I,J)
1000 F_rain_PHY(I,KFLIP,J)=F_rain(L,I,J)
1001 F_RimeF_PHY(I,KFLIP,J)=F_RimeF(L,I,J)
1002 ENDDO
1003 101 CONTINUE
1004 !
1005 END SUBROUTINE EGCP01DRV
1006 !
1007 !-----------------------------------------------------------------------
1008 !
1009 !###############################################################################
1010 ! ***** VERSION OF MICROPHYSICS DESIGNED FOR HIGHER RESOLUTION MESO ETA MODEL
1011 ! (1) Represents sedimentation by preserving a portion of the precipitation
1012 ! through top-down integration from cloud-top. Modified procedure to
1013 ! Zhao and Carr (1997).
1014 ! (2) Microphysical equations are modified to be less sensitive to time
1015 ! steps by use of Clausius-Clapeyron equation to account for changes in
1016 ! saturation mixing ratios in response to latent heating/cooling.
1017 ! (3) Prevent spurious temperature oscillations across 0C due to
1018 ! microphysics.
1019 ! (4) Uses lookup tables for: calculating two different ventilation
1020 ! coefficients in condensation and deposition processes; accretion of
1021 ! cloud water by precipitation; precipitation mass; precipitation rate
1022 ! (and mass-weighted precipitation fall speeds).
1023 ! (5) Assumes temperature-dependent variation in mean diameter of large ice
1024 ! (Houze et al., 1979; Ryan et al., 1996).
1025 ! -> 8/22/01: This relationship has been extended to colder temperatures
1026 ! to parameterize smaller large-ice particles down to mean sizes of MDImin,
1027 ! which is 50 microns reached at -55.9C.
1028 ! (6) Attempts to differentiate growth of large and small ice, mainly for
1029 ! improved transition from thin cirrus to thick, precipitating ice
1030 ! anvils.
1031 ! (7) Top-down integration also attempts to treat mixed-phase processes,
1032 ! allowing a mixture of ice and water. Based on numerous observational
1033 ! studies, ice growth is based on nucleation at cloud top &
1034 ! subsequent growth by vapor deposition and riming as the ice particles
1035 ! fall through the cloud. There are two modes of ice nucleation
1036 ! following Meyers et al. (JAM, 1992):
1037 ! a) Deposition & condensation freezing nucleation - eq. (2.4) when
1038 ! air is supersaturated w/r/t ice
1039 ! b) Contact freezing nucleation - eq. (2.6) in presence of cloud water
1040 ! (8) Depositional growth of newly nucleated ice is calculated for large time
1041 ! steps using Fig. 8 of Miller and Young (JAS, 1979), at 1 deg intervals
1042 ! using their ice crystal masses calculated after 600 s of growth in water
1043 ! saturated conditions. The growth rates are normalized by time step
1044 ! assuming 3D growth with time**1.5 following eq. (6.3) in Young (1993).
1045 ! (9) Ice precipitation rates can increase due to increase in response to
1046 ! cloud water riming due to (a) increased density & mass of the rimed
1047 ! ice, and (b) increased fall speeds of rimed ice.
1048 !###############################################################################
1049 !###############################################################################
1050 !
1051 SUBROUTINE EGCP01COLUMN ( ARAIN, ASNOW, DTPH, RHC_col, &
1052 & I_index, J_index, LSFC, &
1053 & P_col, QI_col, QR_col, Q_col, QW_col, RimeF_col, T_col, &
1054 & THICK_col, WC_col ,KTS,KTE,pcond1d,pidep1d, &
1055 & piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d,pimlt1d, &
1056 & praut1d,pracw1d,prevp1d,pisub1d,pevap1d, DBZ_col,NR_col,NS_col, &
1057 & vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d,INDEXR1d, & !jul28
1058 & RFlag1d) !jun01
1059 !
1060 !###############################################################################
1061 !###############################################################################
1062 !
1063 !-------------------------------------------------------------------------------
1064 !----- NOTE: Code is currently set up w/o threading!
1065 !-------------------------------------------------------------------------------
1066 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
1067 ! . . .
1068 ! SUBPROGRAM: Grid-scale microphysical processes - condensation & precipitation
1069 ! PRGRMMR: Ferrier ORG: W/NP22 DATE: 08-2001
1070 ! PRGRMMR: Jin (Modification for WRF structure)
1071 !-------------------------------------------------------------------------------
1072 ! ABSTRACT:
1073 ! * Merges original GSCOND & PRECPD subroutines.
1074 ! * Code has been substantially streamlined and restructured.
1075 ! * Exchange between water vapor & small cloud condensate is calculated using
1076 ! the original Asai (1965, J. Japan) algorithm. See also references to
1077 ! Yau and Austin (1979, JAS), Rutledge and Hobbs (1983, JAS), and Tao et al.
1078 ! (1989, MWR). This algorithm replaces the Sundqvist et al. (1989, MWR)
1079 ! parameterization.
1080 !-------------------------------------------------------------------------------
1081 !
1082 ! USAGE:
1083 ! * CALL EGCP01COLUMN_hr FROM SUBROUTINE EGCP01DRV
1084 !
1085 ! INPUT ARGUMENT LIST:
1086 ! DTPH - physics time step (s)
1087 ! I_index - I index
1088 ! J_index - J index
1089 ! LSFC - Eta level of level above surface, ground
1090 ! P_col - vertical column of model pressure (Pa)
1091 ! QI_col - vertical column of model ice mixing ratio (kg/kg)
1092 ! QR_col - vertical column of model rain ratio (kg/kg)
1093 ! Q_col - vertical column of model water vapor specific humidity (kg/kg)
1094 ! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
1095 ! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
1096 ! T_col - vertical column of model temperature (deg K)
1097 ! THICK_col - vertical column of model mass thickness (density*height increment)
1098 ! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
1099 ! RHC_col - vertical column of threshold relative humidity for onset of condensation (ratio) !GFDL
1100 !
1101 !
1102 ! OUTPUT ARGUMENT LIST:
1103 ! ARAIN - accumulated rainfall at the surface (kg)
1104 ! ASNOW - accumulated snowfall at the surface (kg)
1105 ! Q_col - vertical column of model water vapor specific humidity (kg/kg)
1106 ! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
1107 ! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
1108 ! QI_col - vertical column of model ice mixing ratio (kg/kg)
1109 ! QR_col - vertical column of model rain ratio (kg/kg)
1110 ! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
1111 ! T_col - vertical column of model temperature (deg K)
1112 ! DBZ_col - vertical column of radar reflectivity (dBZ)
1113 ! NR_col - vertical column of rain number concentration (m^-3)
1114 ! NS_col - vertical column of snow number concentration (m^-3)
1115 !
1116 ! OUTPUT FILES:
1117 ! NONE
1118 !
1119 ! Subprograms & Functions called:
1120 ! * Real Function CONDENSE - cloud water condensation
1121 ! * Real Function DEPOSIT - ice deposition (not sublimation)
1122 ! * Integer Function GET_INDEXR - estimate the mean size of raindrops (microns)
1123 !
1124 ! UNIQUE: NONE
1125 !
1126 ! LIBRARY: NONE
1127 !
1128 ! ATTRIBUTES:
1129 ! LANGUAGE: FORTRAN 90
1130 ! MACHINE : IBM SP
1131 !
1132 !-------------------------------------------------------------------------
1133 !--------------- Arrays & constants in argument list ---------------------
1134 !-------------------------------------------------------------------------
1135 !
1136 IMPLICIT NONE
1137 !
1138 INTEGER,INTENT(IN) :: KTS,KTE,I_index, J_index, LSFC
1139 REAL,INTENT(IN) :: DTPH
1140 REAL,INTENT(INOUT) :: ARAIN, ASNOW
1141 REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: P_col,QI_col,QR_col,RHC_col &
1142 & ,Q_col ,QW_col, RimeF_col, T_col, THICK_col,WC_col,pcond1d &
1143 & ,pidep1d,piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d &
1144 & ,pimlt1d,praut1d,pracw1d,prevp1d,pisub1d,pevap1d,DBZ_col,NR_col &
1145 & ,NS_col,vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d & !jun01
1146 & ,INDEXR1d,RFlag1d !jun01
1147 !
1148 !--------------------------------------------------------------------------------
1149 !--- The following arrays are integral calculations based on the mean
1150 ! snow/graupel diameters, which vary from 50 microns to 1000 microns
1151 ! (1 mm) at 1-micron intervals and assume exponential size distributions.
1152 ! The values are normalized and require being multipled by the number
1153 ! concentration of large ice (NLICE).
1154 !---------------------------------------
1155 ! - VENTI1 - integrated quantity associated w/ ventilation effects
1156 ! (capacitance only) for calculating vapor deposition onto ice
1157 ! - VENTI2 - integrated quantity associated w/ ventilation effects
1158 ! (with fall speed) for calculating vapor deposition onto ice
1159 ! - ACCRI - integrated quantity associated w/ cloud water collection by ice
1160 ! - MASSI - integrated quantity associated w/ ice mass
1161 ! - VSNOWI - mass-weighted fall speed of snow (large ice), used to calculate
1162 ! precipitation rates
1163 ! - VEL_RF - velocity increase of rimed particles as functions of crude
1164 ! particle size categories (at 0.1 mm intervals of mean ice particle
1165 ! sizes) and rime factor (different values of Rime Factor of 1.1**N,
1166 ! where N=0 to Nrime).
1167 !--------------------------------------------------------------------------------
1168 !--- The following arrays are integral calculations based on the mean
1169 ! rain diameters, which vary from 50 microns to 1000 microns
1170 ! (1 mm) at 1-micron intervals and assume exponential size distributions.
1171 ! The values are normalized and require being multiplied by the rain intercept
1172 ! (N0r).
1173 !---------------------------------------
1174 ! - VENTR1 - integrated quantity associated w/ ventilation effects
1175 ! (capacitance only) for calculating evaporation from rain
1176 ! - VENTR2 - integrated quantity associated w/ ventilation effects
1177 ! (with fall speed) for calculating evaporation from rain
1178 ! - ACCRR - integrated quantity associated w/ cloud water collection by rain
1179 ! - MASSR - integrated quantity associated w/ rain
1180 ! - VRAIN - mass-weighted fall speed of rain, used to calculate
1181 ! precipitation rates
1182 ! - RRATE - precipitation rates, which should also be equal to RHO*QR*VRAIN
1183 !
1184 !-------------------------------------------------------------------------
1185 !------- Key parameters, local variables, & important comments ---------
1186 !-----------------------------------------------------------------------
1187 !
1188 !--- TOLER => Tolerance or precision for accumulated precipitation
1189 !
1190 REAL, PARAMETER :: TOLER=5.E-7, C2=1./6., RHO0=1.194, &
1191 Xratio=.025, Zmin=0.01, DBZmin=-20.
1192 !
1193 !--- If BLEND=1:
1194 ! precipitation (large) ice amounts are estimated at each level as a
1195 ! blend of ice falling from the grid point above and the precip ice
1196 ! present at the start of the time step (see TOT_ICE below).
1197 !--- If BLEND=0:
1198 ! precipitation (large) ice amounts are estimated to be the precip
1199 ! ice present at the start of the time step.
1200 !
1201 !--- Extended to include sedimentation of rain on 2/5/01
1202 !
1203 REAL, PARAMETER :: BLEND=1.
1204 !
1205 !--- This variable is for debugging purposes (if .true.)
1206 !
1207 LOGICAL, PARAMETER :: PRINT_diag=.TRUE.
1208 !
1209 !-----------------------------------------------------------------------
1210 !--- Local variables
1211 !-----------------------------------------------------------------------
1212 !
1213 REAL :: EMAIRI, N0r, NLICE, NSmICE, NInuclei, Nrain, Nsnow, Nmix
1214 LOGICAL :: CLEAR, ICE_logical, DBG_logical, RAIN_logical, &
1215 STRAT, DRZL
1216 INTEGER :: INDEX_MY,INDEXR,INDEXR1,INDEXR2,INDEXS,IPASS,ITDX,IXRF,&
1217 & IXS,LBEF,L,INDEXRmin,INDEXS0C,IDR !mar03 !may20
1218 !
1219 !
1220 REAL :: ABI,ABW,AIEVP,ARAINnew,ASNOWnew,BLDTRH,BUDGET, &
1221 & CREVP,DELI,DELR,DELT,DELV,DELW,DENOMF, &
1222 & DENOMI,DENOMW,DENOMWI,DIDEP, &
1223 & DIEVP,DIFFUS,DLI,DTRHO,DUM,DUM1,DUM2,DUM3, &
1224 & DWV0,DWVI,DWVR,DYNVIS,ESI,ESW,FIR,FLIMASS, &
1225 & FWR,FWS,GAMMAR,GAMMAS, &
1226 & PCOND,PIACR,PIACW,PIACWI,PIACWR,PICND,PIDEP,PIDEP_max, &
1227 & PIEVP,PILOSS,PIMLT,PINIT,PP,PRACW,PRAUT,PREVP,PRLOSS, &
1228 & QI,QInew,QLICE,QR,QRnew,QSI,QSIgrd,QSInew,QSW,QSW0, &
1229 & QSWgrd,QSWnew,QT,QTICE,QTnew,QTRAIN,Q,QW,QWnew,Rcw, &
1230 & RFACTOR,RFmx,RFrime,RHO,RIMEF,RIMEF1,RQR,RR,RRHO,SFACTOR, &
1231 & TC,TCC,TFACTOR,THERM_COND,THICK,TK,TK2,TNEW, &
1232 & TOT_ICE,TOT_ICEnew,TOT_RAIN,TOT_RAINnew, &
1233 & VEL_INC,VENTR,VENTIL,VENTIS,VRAIN1,VRAIN2,VRIMEF,VSNOW, &
1234 & VSNOW1,WC,WCnew,WSgrd,WS,WSnew,WV,WVnew, &
1235 & XLI,XLIMASS,XRF, RHgrd, &
1236 & NSImax,QRdum,QSmICE,QLgIce,RQLICE,VCI,TIMLT, &
1237 & RQSnew,RQRnew,Zrain,Zsnow,Ztot,RHOX0C,RFnew,PSDEP,DELS !mar03 !apr22
1238 REAL, SAVE :: Revised_LICE=1.e-3 !-- kg/m**3
1239 !
1240 !#######################################################################
1241 !########################## Begin Execution ############################
1242 !#######################################################################
1243 !
1244 !
1245 ARAIN=0. ! Accumulated rainfall into grid box from above (kg/m**2)
1246 VRAIN1=0. ! Rain fall speeds into grib box from above (m/s)
1247 VSNOW1=0. ! Ice fall speeds into grib box from above (m/s)
1248 ASNOW=0. ! Accumulated snowfall into grid box from above (kg/m**2)
1249 INDEXS0C=MDImin ! Mean snow/graupel diameter just above (<0C) freezing level (height)
1250 RHOX0C=22.5 ! Estimated ice density at 0C (kg m^-3) !mar03
1251 TIMLT=0. ! Total ice melting in a layer (drizzle detection)
1252 STRAT=.FALSE. ! Stratiform rain DSD below melting level !may11
1253 DRZL=.FALSE. ! Drizzle DSD below melting level !may23
1254 !
1255 !-----------------------------------------------------------------------
1256 !------------ Loop from top (L=1) to surface (L=LSFC) ------------------
1257 !-----------------------------------------------------------------------
1258 !
1259 big_loop: DO L=1,LSFC
1260 pcond1d(L)=0.
1261 pidep1d(L)=0.
1262 piacw1d(L)=0.
1263 piacwi1d(L)=0.
1264 piacwr1d(L)=0.
1265 piacr1d(L)=0.
1266 picnd1d(L)=0.
1267 pievp1d(L)=0.
1268 pimlt1d(L)=0.
1269 praut1d(L)=0.
1270 pracw1d(L)=0.
1271 prevp1d(L)=0.
1272 pisub1d(L)=0.
1273 pevap1d(L)=0.
1274 vsnow1d(L)=0.
1275 vrain11d(L)=0.
1276 vrain21d(L)=0.
1277 vci1d(L)=0.
1278 NSmICE1d(L)=0.
1279 DBZ_col(L)=DBZmin
1280 NR_col(L)=0.
1281 NS_col(L)=0.
1282 INDEXR1d(L)=0.
1283 INDEXS1d(L)=0.
1284 RFlag1d(L)=0. !jun01
1285 !
1286 !--- Skip this level and go to the next lower level if no condensate
1287 ! and very low specific humidities
1288 !
1289 !--- Check if any rain is falling into layer from above
1290 !
1291 IF (ARAIN .GT. CLIMIT) THEN
1292 CLEAR=.FALSE.
1293 VRAIN1=0.
1294 ELSE
1295 CLEAR=.TRUE.
1296 ARAIN=0.
1297 ENDIF
1298 !
1299 !--- Check if any ice is falling into layer from above
1300 !
1301 !--- NOTE that "SNOW" in variable names is often synonomous with
1302 ! large, precipitation ice particles
1303 !
1304 IF (ASNOW .GT. CLIMIT) THEN
1305 CLEAR=.FALSE.
1306 VSNOW1=0.
1307 ELSE
1308 ASNOW=0.
1309 ENDIF
1310 !
1311 !-----------------------------------------------------------------------
1312 !------------ Proceed with cloud microphysics calculations -------------
1313 !-----------------------------------------------------------------------
1314 !
1315 TK=T_col(L) ! Temperature (deg K)
1316 TC=TK-T0C ! Temperature (deg C)
1317 PP=P_col(L) ! Pressure (Pa)
1318 Q=Q_col(L) ! Specific humidity of water vapor (kg/kg)
1319 WV=Q/(1.-Q) ! Water vapor mixing ratio (kg/kg)
1320 WC=WC_col(L) ! Grid-scale mixing ratio of total condensate (water or ice; kg/kg)
1321 RHgrd=RHC_col(L) ! Threshold relative humidity for the onset of condensation
1322 !
1323 !-----------------------------------------------------------------------
1324 !--- Moisture variables below are mixing ratios & not specifc humidities
1325 !-----------------------------------------------------------------------
1326 !
1327 !--- This check is to determine grid-scale saturation when no condensate is present
1328 !
1329 ESW=MIN(1000.*FPVS0(TK),0.99*PP) ! Saturation vapor pressure w/r/t water
1330 QSW=EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
1331 WS=QSW ! General saturation mixing ratio (water/ice)
1332 QSI=QSW ! Saturation mixing ratio w/r/t ice
1333 IF (TC .LT. 0.) THEN
1334 ESI=MIN(1000.*FPVS(TK),0.99*PP) ! Saturation vapor pressure w/r/t ice
1335 QSI=EPS*ESI/(PP-ESI) ! Saturation mixing ratio w/r/t water
1336 WS=QSI ! General saturation mixing ratio (water/ice)
1337 ENDIF
1338 !
1339 !--- Effective grid-scale Saturation mixing ratios
1340 !
1341 QSWgrd=RHgrd*QSW
1342 QSIgrd=RHgrd*QSI
1343 WSgrd=RHgrd*WS
1344 !
1345 !--- Check if air is subsaturated and w/o condensate
1346 !
1347 IF (WV.GT.WSgrd .OR. WC.GT.EPSQ) CLEAR=.FALSE.
1348 !
1349 !-----------------------------------------------------------------------
1350 !-- Loop to the end if in clear, subsaturated air free of condensate ---
1351 !-----------------------------------------------------------------------
1352 !
1353 IF (CLEAR) THEN
1354 STRAT=.FALSE. !- Reset stratiform rain flag
1355 DRZL=.FALSE. !- Reset drizzle flag
1356 INDEXRmin=MDRmin !- Reset INDEXRmin
1357 TIMLT=0. !- Reset accumulated ice melting
1358 CYCLE big_loop
1359 ENDIF
1360 !
1361 !-----------------------------------------------------------------------
1362 !--------- Initialize RHO, THICK & microphysical processes -------------
1363 !-----------------------------------------------------------------------
1364 !
1365 !
1366 !--- Virtual temperature, TV=T*(1./EPS-1)*Q, Q is specific humidity;
1367 ! (see pp. 63-65 in Fleagle & Businger, 1963)
1368 !
1369 RHO=PP/(RD*TK*(1.+EPS1*Q)) ! Air density (kg/m**3)
1370 RRHO=1./RHO ! Reciprocal of air density
1371 DTRHO=DTPH*RHO ! Time step * air density
1372 BLDTRH=BLEND*DTRHO ! Blend parameter * time step * air density
1373 THICK=THICK_col(L) ! Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
1374 !
1375 ARAINnew=0. ! Updated accumulated rainfall
1376 ASNOWnew=0. ! Updated accumulated snowfall
1377 QI=QI_col(L) ! Ice mixing ratio
1378 QInew=0. ! Updated ice mixing ratio
1379 QR=QR_col(L) ! Rain mixing ratio
1380 QRnew=0. ! Updated rain ratio
1381 QW=QW_col(L) ! Cloud water mixing ratio
1382 QWnew=0. ! Updated cloud water ratio
1383 !
1384 PCOND=0. ! Condensation (>0) or evaporation (<0) of cloud water (kg/kg)
1385 PIDEP=0. ! Deposition (>0) or sublimation (<0) of ice crystals (kg/kg)
1386 PINIT=0. ! Ice initiation (part of PIDEP calculation, kg/kg)
1387 PIACW=0. ! Cloud water collection (riming) by precipitation ice (kg/kg; >0)
1388 PIACWI=0. ! Growth of precip ice by riming (kg/kg; >0)
1389 PIACWR=0. ! Shedding of accreted cloud water to form rain (kg/kg; >0)
1390 PIACR=0. ! Freezing of rain onto large ice at supercooled temps (kg/kg; >0)
1391 PICND=0. ! Condensation (>0) onto wet, melting ice (kg/kg)
1392 PIEVP=0. ! Evaporation (<0) from wet, melting ice (kg/kg)
1393 PIMLT=0. ! Melting ice (kg/kg; >0)
1394 PRAUT=0. ! Cloud water autoconversion to rain (kg/kg; >0)
1395 PRACW=0. ! Cloud water collection (accretion) by rain (kg/kg; >0)
1396 PREVP=0. ! Rain evaporation (kg/kg; <0)
1397 NSmICE=0. ! Cloud ice number concentration (m^-3)
1398 Nrain=0. ! Rain number concentration (m^-3) !jul28 begin
1399 Nsnow=0. ! "Snow" number concentration (m^-3)
1400 RQRnew=0. ! Final rain content (kg/m**3)
1401 RQSnew=0. ! Final "snow" content (kg/m**3)
1402 Zrain=0. ! Radar reflectivity from rain (mm**6 m-3)
1403 Zsnow=0. ! Radar reflectivity from snow (mm**6 m-3)
1404 Ztot=0. ! Radar reflectivity from rain+snow (mm**6 m-3)
1405 INDEXR=MDRmin ! Mean diameter of rain (microns)
1406 INDEXR1=INDEXR ! 1st updated mean diameter of rain (microns)
1407 INDEXR2=INDEXR ! 2nd updated mean diameter of rain (microns)
1408 N0r=0. ! 1st estimate for rain intercept (m^-4)
1409 DUM1=MIN(0.,TC)
1410 DUM=XMImax*EXP(XMIexp*DUM1)
1411 INDEXS=MIN(MDImax, MAX(MDImin, INT(DUM) ) ) ! 1st estimate for mean diameter of snow (microns)
1412 VCI=0. ! Cloud ice fall speeds (m/s)
1413 VSNOW=0. ! "Snow" (snow/graupel/sleet/hail) fall speeds (m/s)
1414 VRAIN2=0. ! Rain fall speeds out of bottom of grid box (m/s)
1415 RimeF1=1. ! Rime Factor (ratio, >=1, defined below)
1416 !
1417 !--- Double check input hydrometeor mixing ratios
1418 !
1419 ! DUM=WC-(QI+QW+QR)
1420 ! DUM1=ABS(DUM)
1421 ! DUM2=TOLER*MIN(WC, QI+QW+QR)
1422 ! IF (DUM1 .GT. DUM2) THEN
1423 ! WRITE(0,"(/2(a,i4),a,i2)") '{@ i=',I_index,' j=',J_index,
1424 ! & ' L=',L
1425 ! WRITE(0,"(4(a12,g11.4,1x))")
1426 ! & '{@ TCold=',TC,'P=',.01*PP,'DIFF=',DUM,'WCold=',WC,
1427 ! & '{@ QIold=',QI,'QWold=',QW,'QRold=',QR
1428 ! ENDIF
1429 !
1430 !***********************************************************************
1431 !*********** MAIN MICROPHYSICS CALCULATIONS NOW FOLLOW! ****************
1432 !***********************************************************************
1433 !
1434 !--- Calculate a few variables, which are used more than once below
1435 !
1436 !--- Latent heat of vaporization as a function of temperature from
1437 ! Bolton (1980, JAS)
1438 !
1439 TK2=1./(TK*TK) ! 1./TK**2
1440 !
1441 !--- Basic thermodynamic quantities
1442 ! * DYNVIS - dynamic viscosity [ kg/(m*s) ]
1443 ! * THERM_COND - thermal conductivity [ J/(m*s*K) ]
1444 ! * DIFFUS - diffusivity of water vapor [ m**2/s ]
1445 !
1446 TFACTOR=SQRT(TK*TK*TK)/(TK+120.)
1447 DYNVIS=1.496E-6*TFACTOR
1448 THERM_COND=2.116E-3*TFACTOR
1449 DIFFUS=8.794E-5*TK**1.81/PP
1450 !
1451 !--- Air resistance term for the fall speed of ice following the
1452 ! basic research by Heymsfield, Kajikawa, others
1453 !
1454 GAMMAS=MIN(1.5, (1.E5/PP)**C1) !-- limited to 1.5x
1455 !
1456 !--- Air resistance for rain fall speed (Beard, 1985, JAS, p.470)
1457 !
1458 GAMMAR=(RHO0/RHO)**.4
1459 !
1460 !----------------------------------------------------------------------
1461 !------------- IMPORTANT MICROPHYSICS DECISION TREE -----------------
1462 !----------------------------------------------------------------------
1463 !
1464 !--- Determine if conditions supporting ice are present
1465 !
1466 IF (TC.LT.0. .OR. QI.GT. EPSQ .OR. ASNOW.GT.CLIMIT) THEN
1467 ICE_logical=.TRUE.
1468 ELSE
1469 ICE_logical=.FALSE.
1470 QLICE=0.
1471 QTICE=0.
1472 ENDIF
1473 IF (T_ICE <= -100.) THEN
1474 ICE_logical=.FALSE.
1475 QLICE=0.
1476 QTICE=0.
1477 ENDIF
1478 !
1479 !--- Determine if rain is present
1480 !
1481 RAIN_logical=.FALSE.
1482 IF (ARAIN.GT.CLIMIT .OR. QR.GT.EPSQ) RAIN_logical=.TRUE.
1483 !
1484 ice_test: IF (ICE_logical) THEN
1485 !
1486 !--- IMPORTANT: Estimate time-averaged properties.
1487 !
1488 !---
1489 ! -> Small ice particles are assumed to have a mean diameter of 50 microns.
1490 ! * QSmICE - estimated mixing ratio for small cloud ice
1491 !---
1492 ! * TOT_ICE - total mass (small & large) ice before microphysics,
1493 ! which is the sum of the total mass of large ice in the
1494 ! current layer and the input flux of ice from above
1495 ! * PILOSS - greatest loss (<0) of total (small & large) ice by
1496 ! sublimation, removing all of the ice falling from above
1497 ! and the ice within the layer
1498 ! * RimeF1 - Rime Factor, which is the mass ratio of total (unrimed & rimed)
1499 ! ice mass to the unrimed ice mass (>=1)
1500 ! * VrimeF - the velocity increase due to rime factor or melting (ratio, >=1)
1501 ! * VSNOW - Fall speed of rimed snow w/ air resistance correction
1502 ! * VCI - Fall speed of 50-micron ice crystals w/ air resistance correction
1503 ! * EMAIRI - equivalent mass of air associated layer and with fall of snow into layer
1504 ! * XLIMASS - used for debugging, associated with calculating large ice mixing ratio
1505 ! * FLIMASS - mass fraction of large ice
1506 ! * QTICE - time-averaged mixing ratio of total ice
1507 ! * QLICE - time-averaged mixing ratio of large ice
1508 ! * NLICE - time-averaged number concentration of large ice
1509 ! * NSmICE - number concentration of small ice crystals at current level
1510 ! * QSmICE - mixing ratio of small ice crystals at current level
1511 !---
1512 !--- Assumed number fraction of large ice particles to total (large & small)
1513 ! ice particles, which is based on a general impression of the literature.
1514 !
1515 NInuclei=0.
1516 NSmICE=0.
1517 QSmICE=0.
1518 Rcw=0.
1519 IF (TC<0.) THEN
1520 !
1521 !--- Max # conc of small ice crystals based on 10% of total ice content
1522 ! or the parameter NSI_max
1523 !
1524 NSImax=MAX(NSI_max, 0.1*RHO*QI/MASSI(MDImin) ) !aug27
1525 !
1526 !-- Specify Fletcher, Cooper, Meyers, etc. here for ice nuclei concentrations
1527 ! Cooper (1986): NInuclei=MIN(5.*EXP(-0.304*TC), NSImax)
1528 ! Fletcher (1962): NInuclei=MIN(0.01*EXP(-0.6*TC), NSImax)
1529 !
1530 !aug28: The formulas below mean that Fletcher is used for >-21C and Cooper at colder
1531 ! temperatures. In areas of high ice contents near the tops of deep convection,
1532 ! the number concentrations will be determined by the lower value of the "FQi"
1533 ! contribution to NSImax or Cooper.
1534 !
1535 NInuclei=MIN(0.01*EXP(-0.6*TC), NSImax) !aug28 - Fletcher (1962)
1536 NInuclei=MIN(5.*EXP(-0.304*TC), NInuclei) !aug28 - Cooper (1984)
1537 IF (QI>EPSQ) THEN
1538 DUM=RRHO*MASSI(MDImin)
1539 NSmICE=MIN(NInuclei, QI/DUM)
1540 QSmICE=NSmICE*DUM
1541 ENDIF ! End IF (QI>EPSQ)
1542 ENDIF ! End IF (TC<0.)
1543 init_ice: IF (QI<=EPSQ .AND. ASNOW<=CLIMIT) THEN
1544 TOT_ICE=0.
1545 PILOSS=0.
1546 RimeF1=1.
1547 VrimeF=1.
1548 VEL_INC=GAMMAS
1549 VSNOW=0.
1550 VSNOW1=0.
1551 VCI=0.
1552 EMAIRI=THICK
1553 XLIMASS=RimeF1*MASSI(INDEXS)
1554 FLIMASS=1.
1555 QLICE=0.
1556 RQLICE=0.
1557 QTICE=0.
1558 NLICE=0.
1559 ELSE init_ice
1560 !
1561 !--- For T<0C mean particle size follows Houze et al. (JAS, 1979, p. 160),
1562 ! converted from Fig. 5 plot of LAMDAs. Similar set of relationships
1563 ! also shown in Fig. 8 of Ryan (BAMS, 1996, p. 66).
1564 !
1565 !
1566 !sep10 - Start of changes described in (23) at top of code.
1567 !
1568 TOT_ICE=THICK*QI+BLEND*ASNOW
1569 PILOSS=-TOT_ICE/THICK
1570 QLgICE=MAX(0., QI-QSmICE) !-- 1st-guess estimate of large ice
1571 VCI=GAMMAS*VSNOWI(MDImin)
1572 !
1573 !-- Need to save this original value before two_pass iteration
1574 !
1575 LBEF=MAX(1,L-1)
1576 RimeF1=(RimeF_col(L)*THICK*QLgICE &
1577 & +RimeF_col(LBEF)*BLEND*ASNOW)/TOT_ICE
1578 !
1579 !mar08 see (32); !apr22a see (41) start - Estimate mean diameter (INDEXS in microns)
1580 IF (RimeF1>2.) THEN
1581 DUM3=RH_NgC*(RHO*QLgICE)**C1 !- convective mean diameter
1582 DUM2=RH_NgT*(RHO*QLgICE)**C1 !- transition mean diameter
1583 IF (RimeF1>=10.) THEN
1584 DUM=DUM3
1585 ELSE IF (RimeF1>=5.) THEN
1586 DUM=0.2*(RimeF1-5.) !- Blend at 5<=RF<10
1587 DUM=DUM3*DUM+DUM2*(1.-DUM)
1588 ELSE
1589 DUM1=REAL(INDEXS) !- stratiform mean diameter
1590 DUM=0.33333*(RimeF1-2.) !- Blend at 2<RF<5
1591 DUM=DUM2*DUM+DUM1*(1.-DUM)
1592 ENDIF
1593 INDEXS=MIN(MDImax, MAX(MDImin, INT(DUM) ) )
1594 ENDIF
1595 !may11 - begin
1596 EMAIRI=THICK+BLDTRH*VSNOW1
1597 QLICE=(THICK*QLgICE+BLEND*ASNOW)/EMAIRI !- Estimate of large ice
1598 !may17 - start; now calculated before the DO IPASS iteration
1599 RQLICE=RHO*QLICE
1600 QTICE=QLICE+QSmICE
1601 FLIMASS=QLICE/QTICE
1602 !may17 end
1603 !may11 - end
1604 !mar08 !apr22a end
1605 !
1606 !===========================================
1607 two_pass: DO IPASS=1,2
1608 !===========================================
1609 !
1610 !-- Prevent rime factor (RimeF1) from exceeding a maximum value, RFmx, in which
1611 ! the ice has an equivalent density near that of pure ice
1612 !
1613 DUM=1.E-6*REAL(INDEXS) !- Mean diameter in m
1614 RFmx=RFmx1*DUM*DUM*DUM/MASSI(INDEXS)
1615 RimeF1=MIN(RimeF1, RFmx)
1616 !
1617 vel_rime: IF (RimeF1<=1.) THEN
1618 RimeF1=1.
1619 VrimeF=1.
1620 ELSE vel_rime
1621 !--- Prevent rime factor (RimeF1) from exceeding a maximum value (RFmax)
1622 RimeF1=MIN(RimeF1, RFmax)
1623 IXS=MAX(2, MIN(INDEXS/100, 9))
1624 XRF=10.492*ALOG(RimeF1)
1625 IXRF=MAX(0, MIN(INT(XRF), Nrime))
1626 IF (IXRF .GE. Nrime) THEN
1627 VrimeF=VEL_RF(IXS,Nrime)
1628 ELSE
1629 VrimeF=VEL_RF(IXS,IXRF)+(XRF-FLOAT(IXRF))* &
1630 & (VEL_RF(IXS,IXRF+1)-VEL_RF(IXS,IXRF))
1631 ENDIF
1632 VrimeF=MAX(1., VrimeF)
1633 ENDIF vel_rime
1634 !
1635 !may17 - start
1636 VEL_INC=GAMMAS*VrimeF !- Normal rimed ice fall speeds
1637 VSNOW=VEL_INC*VSNOWI(INDEXS)
1638 IF (RimeF1>=5. .AND. INDEXS==MDImax .AND. RQLICE>RQhail) THEN
1639 !- Additional increase using Thompson graupel/hail fall speeds
1640 DUM=GAMMAS*AVhail*RQLICE**BVhail
1641 IF (DUM>VSNOW) THEN
1642 VSNOW=DUM
1643 VEL_INC=VSNOW/VSNOWI(INDEXS)
1644 ENDIF
1645 ENDIF
1646 XLIMASS=RimeF1*MASSI(INDEXS)
1647 NLICE=RQLICE/XLIMASS
1648 !
1649 !sep16 - End of change described in (26) at top of code.
1650 !
1651 !===========================================
1652 IF (IPASS>=2 .OR. &
1653 (NLICE>=NLImin .AND. NLICE<=NSI_max)) EXIT two_pass
1654 !may17 - end
1655 !===========================================
1656 !
1657 !--- Force NLICE to be between NLImin and NSI_max when IPASS=1
1658 !
1659 ! IF (PRINT_diag .AND. RQLICE>Revised_LICE) DUM2=NLICE !-- For debugging (see DUM2 below)
1660 NLICE=MAX(NLImin, MIN(NSI_max, NLICE) )
1661 !sep16 - End of changes
1662 !
1663 XLI=RQLICE/(NLICE*RimeF1) !- Mean mass of unrimed ice
1664 new_size: IF (XLI<=MASSI(MDImin) ) THEN
1665 INDEXS=MDImin
1666 ELSE IF (XLI<=MASSI(450) ) THEN new_size
1667 DLI=9.5885E5*XLI**.42066 ! DLI in microns
1668 INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
1669 ELSE IF (XLI<MASSI(MDImax) ) THEN new_size
1670 DLI=3.9751E6*XLI**.49870 ! DLI in microns
1671 INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
1672 ELSE new_size
1673 INDEXS=MDImax
1674 ! IF (PRINT_diag .AND. RQLICE>Revised_LICE) THEN
1675 ! WRITE(0,"(5(a12,g11.4,1x))") '{$ RimeF1=',RimeF1, &
1676 ! & ' RHO*QLICE=',RQLICE,' TC=',TC,' NLICE=',NLICE, &
1677 ! & ' NLICEold=',DUM2
1678 ! Revised_LICE=1.2*RQLICE
1679 ! ENDIF
1680 ENDIF new_size
1681 !===========================================
1682 ENDDO two_pass
1683 !===========================================
1684 ENDIF init_ice
1685 ENDIF ice_test
1686 !
1687 !mar03 !may11 !jun01 - start; see (34) above
1688 IF(TC<=0.) THEN
1689 STRAT=.FALSE.
1690 INDEXRmin=MDRmin
1691 TIMLT=0.
1692 INDEXS0C=INDEXS
1693 RHOX0C=22.5*MAX(1.,MIN(RimeF1,40.)) !- Estimated ice density at 0C (kg m^-3)
1694 ELSE ! TC>0.
1695 IF(.NOT.RAIN_logical) THEN
1696 STRAT=.FALSE. !- Reset STRAT
1697 INDEXRmin=MDRmin !- Reset INDEXRmin
1698 IF(.NOT.ICE_logical) TIMLT=0.
1699 ELSE
1700 !- STRAT=T for stratiform rain
1701 IF(TIMLT>EPSQ .AND. RHOX0C<=225.) THEN
1702 STRAT=.TRUE.
1703 ELSE
1704 STRAT=.FALSE. !- Reset STRAT
1705 INDEXRmin=MDRmin !- Reset INDEXRmin
1706 ENDIF
1707 IF(STRAT .AND. INDEXRmin<=MDRmin) THEN
1708 INDEXRmin=INDEXS0C*(0.001*RHOX0C)**C1
1709 INDEXRmin=MAX(MDRmin, MIN(INDEXRmin, INDEXRstrmax) )
1710 ENDIF
1711 ENDIF
1712 ENDIF
1713 !
1714 IF(STRAT .OR. TIMLT>EPSQ) THEN
1715 DRZL=.FALSE.
1716 ELSE
1717 !- DRZL=T for drizzle (no melted ice falling from above)
1718 DRZL=.TRUE. !mar30
1719 ENDIF
1720 !jun01 - end
1721 !
1722 !----------------------------------------------------------------------
1723 !--------------- Calculate individual processes -----------------------
1724 !----------------------------------------------------------------------
1725 !
1726 !--- Cloud water autoconversion to rain (PRAUT) and collection of cloud
1727 ! water by precipitation ice (PIACW)
1728 !
1729 IF (QW.GT.EPSQ .AND. TC.GE.T_ICE) THEN
1730 !-- The old autoconversion threshold returns
1731 DUM2=RHO*QW
1732 IF (DUM2>QAUT0) THEN
1733 !-- July 2010 version follows Liu & Daum (JAS, 2004) and Liu et al. (JAS, 2006)
1734 DUM2=DUM2*DUM2
1735 DUM=BRAUT*DUM2*QW
1736 DUM1=ARAUT*DUM2
1737 PRAUT=MIN(QW, DUM*(1.-EXP(-DUM1*DUM1)) )
1738 ENDIF
1739 IF (QLICE .GT. EPSQ) THEN
1740 !
1741 !--- Collection of cloud water by large ice particles ("snow")
1742 ! PIACWI=PIACW for riming, PIACWI=0 for shedding
1743 !
1744 FWS=MIN(1., CIACW*VEL_INC*NLICE*ACCRI(INDEXS) ) !jul28 (16)
1745 PIACW=FWS*QW
1746 IF (TC<0.) THEN
1747 PIACWI=PIACW !- Large ice riming
1748 Rcw=ARcw*DUM2**C1 !- Cloud droplet radius in microns
1749 ENDIF
1750 ENDIF ! End IF (QLICE .GT. EPSQ)
1751 ENDIF ! End IF (QW.GT.EPSQ .AND. TC.GE.T_ICE)
1752 !
1753 !----------------------------------------------------------------------
1754 !--- Calculate homogeneous freezing of cloud water (PIACW, PIACWI) and
1755 ! ice deposition (PIDEP), which also includes ice initiation (PINIT)
1756 !
1757 ice_only: IF (TC.LT.T_ICE .AND. (WV.GT.QSWgrd .OR. QW.GT.EPSQ)) THEN
1758 !
1759 !--- Adjust to ice saturation at T<T_ICE (-40C) if saturated w/r/t water
1760 ! or if cloud water is present (homogeneous glaciation).
1761 !
1762 PIACW=QW
1763 PIACWI=PIACW
1764 Rcw=0. ! Homogeneous freezing of cloud water adds to
1765 ! cloud ice, do not use to update RF for large ice
1766 DUM1=TK+XLF1*PIACW ! Updated (dummy) temperature (deg K)
1767 DUM2=WV ! Updated (dummy) water vapor mixing ratio
1768 DUM=MIN(1000.*FPVS(DUM1),0.99*PP) ! Updated (dummy) saturation vapor pressure w/r/t ice
1769 DUM=RHgrd*EPS*DUM/(PP-DUM) ! Updated (dummy) saturation mixing ratio w/r/t ice
1771 !-- Debug 20120111
1772 IF (WARN1 .AND. DUM1<XMIN) THEN
1773 WRITE(0,*) 'WARN1: Water saturation T<180K; I,J,L,TC,P=', &
1774 I_index,J_index,L,DUM1-T0C,.01*PP
1775 WARN1=.FALSE.
1776 ENDIF
1778 !-- Adjust to ice saturation if homogeneous freezing occurs
1779 IF (DUM2>DUM) &
1780 & PIDEP=DEPOSIT(PP,DUM1,DUM2,RHgrd,I_index,J_index,L)
1782 DWVi=0. ! Used only for debugging
1783 !
1784 ELSE IF (TC<0.) THEN ice_only
1785 !
1786 !--- These quantities are handy for ice deposition/sublimation
1787 ! PIDEP_max - max deposition or minimum sublimation to ice saturation
1788 !
1789 DENOMI=1.+XLS3*QSI*TK2
1790 ! DWVi=MIN(WV,QSW)-QSIgrd !may17
1791 DWVi=WV-QSIgrd !may17 !-- Speed up ice deposition
1792 PIDEP_max=MAX(PILOSS, DWVi/DENOMI)
1793 DIDEP=0. !-- Vapor deposition/sublimation onto existing ice
1794 IF (QTICE .GT. 0.) THEN
1795 !
1796 !--- Calculate ice deposition/sublimation
1797 ! * SFACTOR - [VEL_INC**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
1798 ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
1799 ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ;
1800 ! VENTIL, VENTIS - m**-2 ; VENTI1 - m ;
1801 ! VENTI2 - m**2/s**.5 ; DIDEP - unitless
1802 !
1803 ABI=1./(RHO*XLS2*QSI*TK2/THERM_COND+1./DIFFUS)
1804 !
1805 !--- VENTIL - Number concentration * ventilation factors for large ice
1806 !--- VENTIS - Number concentration * ventilation factors for small ice
1807 !
1808 !--- Variation in the number concentration of ice with time is not
1809 ! accounted for in these calculations (could be in the future).
1810 !
1811 DUM=(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
1812 SFACTOR=SQRT(GAMMAS)*DUM !-- GAMMAS (RF=1 for cloud ice)
1813 VENTIS=(VENTI1(MDImin)+SFACTOR*VENTI2(MDImin))*NSmICE
1814 SFACTOR=SQRT(VEL_INC)*DUM
1815 VENTIL=(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS))*NLICE
1816 DIDEP=ABI*(VENTIL+VENTIS)*DTPH
1817 ENDIF !-- IF (QTICE .GT. 0.) THEN
1818 !
1819 !--- WARNING: the Meyers et al. stuff is not used!
1820 !--- Two modes of ice nucleation from Meyers et al. (JAM, 1992):
1821 ! (1) Deposition & condensation freezing nucleation - eq. (2.4),
1822 ! requires ice supersaturation
1823 ! (2) Contact freezing nucleation - eq. (2.6), requires presence
1824 ! of cloud water
1825 !
1826 !--- Ice crystal growth during the physics time step is calculated using
1827 ! Miller & Young (1979, JAS) and is represented by MY_GROWTH(INDEX_MY),
1828 ! where INDEX_MY is tabulated air temperatures between -1 and -35 C.
1829 ! The original Miller & Young (MY) calculations only went down to -30C,
1830 ! so a fixed value is assumed at temperatures colder than -30C.
1831 !
1832 !--- Sensitivity test using:
1833 ! - Fletcher (1962) for ice initiation
1834 ! - Can occur only when there is water supersaturation and T<-12C.
1835 ! - Maximum cloud ice number concentration of 100 per liter
1836 !
1837 IF (WV>QSWgrd .AND. TC<T_ICE_init .AND. NSmICE<NInuclei) THEN
1838 INDEX_MY=MAX(MY_T1, MIN( INT(.5-TC), MY_T2 ) )
1839 !-- Only initiate small ice to get to NInuclei number concentrations
1840 DUM=MAX(0., NInuclei-NSmICE)
1841 !PINIT=MAX(0., DUM*MY_GROWTH_NMM(INDEX_MY)*RRHO)
1842 PINIT=MAX(0., DUM*MY_GROWTH(INDEX_MY)*RRHO)
1843 ENDIF
1844 !
1845 !--- Calculate PIDEP, but also account for limited water vapor supply
1846 !
1847 IF (DWVi>0.) THEN
1848 PIDEP=MIN(DWVI*DIDEP+PINIT, PIDEP_max)
1849 ELSE IF (DWVi<0.) THEN
1850 PIDEP=MAX(DWVi*DIDEP, PIDEP_max)
1851 ENDIF
1852 !
1853 ENDIF ice_only
1854 !
1855 !------------------------------------------------------------------------
1856 !--- Cloud water condensation
1857 !
1858 IF (TC>=T_ICE .AND. (QW>EPSQ .OR. WV>QSWgrd)) THEN
1859 !apr21 - start; see (39) above
1860 DUM1=QW-PIACWI
1861 DUM2=TK+XLS1*PIDEP+XLF1*PIACWI
1862 DUM3=WV-PIDEP
1863 PCOND=CONDENSE(PP,DUM1,DUM2,DUM3,RHgrd,I_index,J_index,L)
1864 !apr21 - end; see (39) above
1865 ENDIF
1866 !
1867 !--- Limit freezing of accreted rime to prevent temperature oscillations,
1868 ! a crude Schumann-Ludlam limit (p. 209 of Young, 1993).
1869 !
1870 TCC=TC+XLV1*PCOND+XLS1*PIDEP+XLF1*PIACWI
1871 IF (TCC>0.) THEN
1872 PIACWI=0.
1873 Rcw=0. !apr18 see (38)
1874 PIDEP=0. !apr18 see (38)
1875 TCC=TC+XLV1*PCOND !apr18 see (38)
1876 ENDIF
1877 !
1878 QSW0=0.
1879 IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical) THEN
1880 !
1881 !--- Calculate melting and evaporation/condensation
1882 ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ;
1883 ! VENTIL - m**-2 ; VENTI1 - m ;
1884 ! VENTI2 - m**2/s**.5 ; CIEVP - /s
1885 !
1886 SFACTOR=SQRT(VEL_INC)*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
1887 VENTIL=NLICE*(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS))
1888 AIEVP=VENTIL*DIFFUS*DTPH
1889 IF (AIEVP .LT. Xratio) THEN
1890 DIEVP=AIEVP
1891 ELSE
1892 DIEVP=1.-EXP(-AIEVP)
1893 ENDIF
1894 QSW0=EPS*ESW0/(PP-ESW0)
1895 DWV0=MIN(WV,QSW)-QSW0
1896 DUM=QW+PCOND
1897 IF (WV.LT.QSW .AND. DUM.LE.EPSQ) THEN
1898 !
1899 !--- Evaporation from melting snow (sink of snow) or shedding
1900 ! of water condensed onto melting snow (source of rain)
1901 !
1902 DUM=DWV0*DIEVP
1903 PIEVP=MAX( MIN(0., DUM), PILOSS)
1904 PICND=MAX(0., DUM)
1905 ENDIF ! End IF (WV.LT.QSW .AND. DUM.LE.EPSQ)
1906 PIMLT=THERM_COND*TCC*VENTIL*RRHO*DTPH/XLF
1907 !
1908 !--- Limit melting to prevent temperature oscillations across 0C
1909 !
1910 DUM1=MAX( 0., (TCC+XLV1*PIEVP)/XLF1 )
1911 PIMLT=MIN(PIMLT, DUM1)
1912 !
1913 !--- Limit loss of snow by melting (>0) and evaporation
1914 !
1915 DUM=PIEVP-PIMLT
1916 IF (DUM .LT. PILOSS) THEN
1917 DUM1=PILOSS/DUM
1918 PIMLT=PIMLT*DUM1
1919 PIEVP=PIEVP*DUM1
1920 ENDIF ! End IF (DUM .GT. QTICE)
1921 ENDIF ! End IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical)
1922 !
1923 !--- IMPORTANT: Estimate time-averaged properties.
1924 !
1925 ! * TOT_RAIN - total mass of rain before microphysics, which is the sum of
1926 ! the total mass of rain in the current layer and the input
1927 ! flux of rain from above
1928 ! * VRAIN1 - fall speed of rain into grid from above
1929 ! * VRAIN2 - fall speed of rain out of grid box to the level below
1930 ! * QTRAIN - time-averaged mixing ratio of rain (kg/kg)
1931 ! * PRLOSS - greatest loss (<0) of rain, removing all rain falling from
1932 ! above and the rain within the layer
1933 ! * RQR - rain content (kg/m**3)
1934 ! * INDEXR - mean size of rain drops to the nearest 1 micron in size
1935 ! * N0r - intercept of rain size distribution (typically 10**6 m**-4)
1936 !
1937 TOT_RAIN=0.
1938 QTRAIN=0.
1939 PRLOSS=0.
1940 RQR=0.
1941 do_rain: IF (RAIN_logical) THEN
1942 TOT_RAIN=THICK*QR+BLEND*ARAIN !aug27 begin mods
1943 PRLOSS=-TOT_RAIN/THICK
1944 !may11 - start
1945 QTRAIN=TOT_RAIN/(THICK+BLDTRH*VRAIN1) !- Rain mixing ratio
1946 RQR=RHO*QTRAIN !- Rain content
1947 IF(STRAT .AND. RQR>1.E-3) STRAT=.FALSE. !may31
1948 IF(DRZL .AND. PIMLT>EPSQ) DRZL=.FALSE. !jun01
1949 !
1950 DSD1: IF (RQR<=RQR_DRmin) THEN
1951 !-- Extremely light rain
1952 INDEXR=MDRmin
1953 N0r=RQR/MASSR(INDEXR)
1954 ELSE IF (DRZL) THEN DSD1
1955 !-- Drizzle - gradually increase N0r for rain contents <0.5 g m^-3 !may31
1956 DUM=MAX(1.0, 0.5E-3/RQR)
1957 N0r=MIN(1.E9, N0r0*DUM*SQRT(DUM) ) !- 8.e6 <= N0r <= 1.e9
1958 DUM1=RQR/(PI*RHOL*N0r)
1959 DUM=1.E6*SQRT(SQRT(DUM1))
1960 INDEXR=MAX(XMRmin, MIN(DUM, XMRmax) )
1961 N0r=RQR/MASSR(INDEXR)
1962 !may20 - start
1963 ELSE IF (RQR>=RQR_DRmax) THEN DSD1
1964 !-- Extremely heavy rain
1965 INDEXR=MDRmax
1966 N0r=RQR/MASSR(INDEXR)
1967 ELSE DSD1
1968 !-- Light to heavy rain
1969 N0r=N0r0
1970 DUM=CN0r0*SQRT(SQRT(RQR))
1971 INDEXR=MAX(XMRmin, MIN(DUM, XMRmax) )
1972 IF (STRAT .AND. INDEXR<INDEXRmin) THEN
1973 !-- Stratiform rain
1974 INDEXR=INDEXRmin
1975 N0r=RQR/MASSR(INDEXR)
1976 ENDIF
1977 ENDIF DSD1
1978 !may20 - end
1979 VRAIN2=GAMMAR*VRAIN(INDEXR) !-- VRAIN2 is used below
1980 !may11 - end
1981 !
1982 IF (TC .LT. T_ICE) THEN
1983 PIACR=-PRLOSS
1984 ELSE
1985 DWVr=WV-PCOND-QSW
1986 DUM=QW+PCOND
1987 IF (DWVr.LT.0. .AND. DUM.LE.EPSQ) THEN
1988 !
1989 !--- Rain evaporation
1990 !
1991 ! * RFACTOR - [GAMMAR**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
1992 ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
1993 !
1994 ! * Units: RFACTOR - s**.5/m ; ABW - m**2/s ; VENTR - m**-2 ;
1995 ! N0r - m**-4 ; VENTR1 - m**2 ; VENTR2 - m**3/s**.5 ;
1996 ! CREVP - unitless
1997 !
1998 RFACTOR=SQRT(GAMMAR)*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
1999 ABW=1./(RHO*XLV2*QSW*TK2/THERM_COND+1./DIFFUS)
2000 !
2001 !--- Note that VENTR1, VENTR2 lookup tables do not include the
2002 ! 1/Davg multiplier as in the ice tables
2003 !
2004 VENTR=N0r*(VENTR1(INDEXR)+RFACTOR*VENTR2(INDEXR))
2005 CREVP=ABW*VENTR*DTPH
2006 PREVP=MAX(DWVr*CREVP, PRLOSS)
2007 ELSE IF (QW .GT. EPSQ) THEN
2008 FWR=CRACW*GAMMAR*N0r*ACCRR(INDEXR)
2009 PRACW=MIN(1.,FWR)*QW
2010 ENDIF ! End IF (DWVr.LT.0. .AND. DUM.LE.EPSQ)
2011 !
2012 IF (ICE_logical .AND. TC<0. .AND. TCC<0.) THEN
2013 !
2014 !--- Biggs (1953) heteorogeneous freezing (e.g., Lin et al., 1983)
2015 ! - Rescaled mean drop diameter from microns (INDEXR) to mm (DUM) to prevent underflow
2016 !
2017 DUM=.001*FLOAT(INDEXR)
2018 DUM=(EXP(ABFR*TC)-1.)*DUM*DUM*DUM*DUM*DUM*DUM*DUM
2019 PIACR=MIN(CBFR*N0r*RRHO*DUM, QTRAIN)
2020 !- Start changes listed above in (33) - no collisional freezing of rain
2021 FIR=0.
2023 IF (QLICE .GT. EPSQ) THEN
2024 !
2025 !--- Freezing of rain by collisions w/ large ice
2026 !
2027 DUM=GAMMAR*VRAIN(INDEXR)
2028 DUM1=DUM-VSNOW
2029 !
2030 !--- DUM2 - Difference in spectral fall speeds of rain and
2031 ! large ice, parameterized following eq. (48) on p. 112 of
2032 ! Murakami (J. Meteor. Soc. Japan, 1990)
2033 !
2034 DUM2=SQRT(DUM1*DUM1+.04*DUM*VSNOW)
2035 DUM1=5.E-12*INDEXR*INDEXR+2.E-12*INDEXR*INDEXS &
2036 & +.5E-12*INDEXS*INDEXS
2037 FIR=MIN(1., CIACR*NLICE*DUM1*DUM2)
2038 !
2039 !--- Future? Should COLLECTION BY SMALL ICE BE INCLUDED???
2040 !
2041 PIACR=MIN(PIACR+FIR*QTRAIN, QTRAIN)
2042 ENDIF ! End IF (QLICE .GT. EPSQ)
2043 DUM=PREVP-PIACR
2044 If (DUM .LT. PRLOSS) THEN
2045 DUM1=PRLOSS/DUM
2046 PREVP=DUM1*PREVP
2047 PIACR=DUM1*PIACR
2048 ENDIF ! End If (DUM .LT. PRLOSS)
2049 ENDIF ! End IF (TC.LT.0. .AND. TCC.LT.0.)
2050 ENDIF ! End IF (TC .LT. T_ICE)
2051 ENDIF do_rain ! End IF (RAIN_logical)
2052 !
2053 INDEXR1=INDEXR
2054 !
2055 !----------------------------------------------------------------------
2056 !---------------------- Main Budget Equations -------------------------
2057 !----------------------------------------------------------------------
2058 !
2059 !
2060 !-----------------------------------------------------------------------
2061 !--- Update fields, determine characteristics for next lower layer ----
2062 !-----------------------------------------------------------------------
2063 !
2064 !--- Carefully limit sinks of cloud water
2065 !
2066 DUM1=PIACW+PRAUT+PRACW-MIN(0.,PCOND)
2067 IF (DUM1 .GT. QW) THEN
2068 DUM=QW/DUM1
2069 PIACW=DUM*PIACW
2070 PIACWI=DUM*PIACWI
2071 PRAUT=DUM*PRAUT
2072 PRACW=DUM*PRACW
2073 IF (PCOND .LT. 0.) PCOND=DUM*PCOND
2074 ENDIF
2075 PIACWR=PIACW-PIACWI ! TC >= 0C
2076 !
2077 !--- QWnew - updated cloud water mixing ratio
2078 !
2079 DELW=PCOND-PIACW-PRAUT-PRACW
2080 QWnew=QW+DELW
2081 IF (QWnew .LE. EPSQ) QWnew=0.
2082 IF (QW.GT.0. .AND. QWnew.NE.0.) THEN
2083 DUM=QWnew/QW
2084 IF (DUM .LT. TOLER) QWnew=0.
2085 ENDIF
2086 !
2087 !--- Update temperature and water vapor mixing ratios
2088 !
2089 DELT= XLV1*(PCOND+PIEVP+PICND+PREVP) &
2090 & +XLS1*PIDEP+XLF1*(PIACWI+PIACR-PIMLT)
2091 Tnew=TK+DELT
2092 !
2093 DELV=-PCOND-PIDEP-PIEVP-PICND-PREVP
2094 WVnew=WV+DELV
2095 !
2096 !--- Update ice mixing ratios
2097 !
2098 !---
2099 ! * TOT_ICEnew - total mass (small & large) ice after microphysics,
2100 ! which is the sum of the total mass of ice in the
2101 ! layer and the flux of ice out of the grid box below
2102 ! * RimeF - Rime Factor, which is the mass ratio of total (unrimed &
2103 ! rimed) ice mass to the unrimed ice mass (>=1)
2104 ! * QInew - updated mixing ratio of total (large & small) ice in layer
2105 ! * QLICEnew=FLIMASS*QInew, an estimate of the updated large ice mixing ratio
2106 !
2107 ! -> TOT_ICEnew=QInew*THICK+QLICEnew*BLDTRH*VSNOW+(QInew-QLICEnew)*BLDTRH*VCI
2108 ! =QInew*THICK+QInew*FLIMASS*BLDTRH*VSNOW+QInew*(1.-FLIMASS)*BLDTRH*VCI
2109 ! =QInew*THICK+QInew*BLDTRH*(FLIMASS*VSNOW+(1.-FLIMASS)*VCI)
2110 ! =QInew*(THICK+BLDTRH*(FLIMASS*VSNOW+(1.-FLIMASS)*VCI))
2111 ! -> Rearranging this equation gives:
2112 ! QInew=TOT_ICEnew/(THICK+BLDTRH*(FLIMASS*VSNOW+(1.-FLIMASS)*VCI))
2113 !
2114 ! * ASNOWnew - updated accumulation of snow and cloud ice at bottom of grid cell
2115 ! -> ASNOWnew=QInew*BLDTRH*(FLIMASS+VSNOW+(1.-FLIMASS)*VCI)
2116 !---
2117 !
2118 DELI=0.
2119 RimeF=1.
2120 IF (ICE_logical) THEN
2121 DELI=PIDEP+PIEVP+PIACWI+PIACR-PIMLT
2122 TOT_ICEnew=TOT_ICE+THICK*DELI
2123 IF (TOT_ICE.GT.0. .AND. TOT_ICEnew.NE.0.) THEN
2124 DUM=TOT_ICEnew/TOT_ICE
2125 IF (DUM .LT. TOLER) TOT_ICEnew=0.
2126 ENDIF
2127 IF (TOT_ICEnew .LE. CLIMIT) THEN
2128 TOT_ICEnew=0.
2129 RimeF=1.
2130 QInew=0.
2131 ASNOWnew=0.
2132 ELSE
2133 !
2134 !--- Update rime factor if appropriate
2135 !
2136 !apr22 - start; see (40) above
2137 IF (PINIT<=EPSQ) THEN
2138 PSDEP=MAX(0., PIDEP)
2139 ELSE
2140 !-- Assume vapor deposition is onto cloud ice if PINT>0
2141 PSDEP=0.
2142 ENDIF
2143 DELS=PIACWI+PIACR+PSDEP
2144 IF (DELS<=EPSQ .OR. Tnew>=T0C) THEN
2145 RimeF=RimeF1
2146 ELSE
2147 !
2148 !-----------------------------------------------------------------------
2149 ! RimeF is the new, updated rime factor; RimeF1 is the existing rime factor
2150 !
2151 ! From near line 1115:
2152 ! RimeF1=(RimeF_col(L)*THICK*QLgICE+RimeF_col(LBEF)*BLEND*ASNOW)/TOT_ICE
2153 ! RimeF1*TOT_ICE=RimeF_col(L)*THICK*QLgICE+RimeF_col(LBEF)*BLEND*ASNOW
2154 !
2155 ! TOT_ICEnew=TOT_ICE+THICK*DELI
2156 ! TOT_ICEnew=TOT_ICE+THICK*(PIDEP+PIEVP+PIACWI+PIACR-PIMLT)
2157 !
2158 ! But the following processes do not affect the rime factor (ice density):
2159 ! 1) PIEVP, evaporation from melting ice
2160 ! 2) PIDEP<0, sublimation of ice
2161 ! 3) PIMLT, melting of ice because it is shed to rain
2162 !
2163 ! So the final version is
2164 ! TOT_ICEnew=TOT_ICE+THICK*DELS,
2165 ! DELS=PSDEP+PIACWI+PIACR,
2166 ! PSDEP=MAX(0., PIDEP) only if PINIT<=EPSQ
2167 !
2168 !-----------------------------------------------------------------------
2169 ! Rime factor values associated with different processes.
2170 !-----------------------------------------------------------------------
2171 ! 1) RF=1 for growth by vapor deposition.
2172 ! 2) RFmx, the RF associated with an ice density of 900 kg m^-3, the assumed
2173 ! properties for the freezing of rain to ice (PIACR).
2174 ! 3) RFrime, the RF associated with the density of rimed ice, the assumed
2175 ! properties for the riming of cloud water onto ice (PIACWI).
2176 !
2177 ! The rime factor (RFnew) associated with the various microphysical processes is:
2178 ! RFnew=(1.*PSDEP+RFrime*PIACWI+RFmx*PIACR)/DELS,
2179 !
2180 ! The updated rime factor (RimeF) becomes:
2181 ! RimeF*TOT_ICEnew=RimeF1*TOT_ICE+THICK*DELS*RFnew,
2182 ! RimeF=(RimeF1*TOT_ICE+THICK*DELS*RFnew)/TOT_ICEnew
2183 !-----------------------------------------------------------------------
2184 !
2185 !-- Calculate density of rimed ice following Heymsfield and Pflaum (1985), where
2186 ! Rime density = 300.*(-Rcw*Vimpact/Tsfc)**0.44,
2187 ! in which Rcw is the mean diameter of the cloud droplets, Vimpact=VSNOW, and
2188 ! Tsfc (surface temperature of the particle, deg C) is approximated by TC.
2189 !
2190 ! Here the calculations are extended to temperatures as warm as -0.5C, whereas
2191 ! the original study only looked at ice particles whose surface temperature were
2192 ! colder than -5C. Rime densities will vary from 170 to 900 kg m^-3
2193 ! (Straka, 2009 textbook; p. 308).
2194 !
2195 DUM=1.E-6*REAL(INDEXS) !- Mean diameter in m
2196 DUM1=DUM*DUM*DUM/MASSI(INDEXS) !- Used to estimate ice densities
2197 RFmx=RFmx1*DUM1 !- Rime factor for density of pure ice (PIACR)
2198 IF (PIACWI>0. .AND. Rcw>0.) THEN
2199 DUM=-Rcw*VSNOW/MIN(-0.5,TC)
2200 DUM=MIN(12.14, MAX(0.275, DUM) )
2201 DUM=300.*DUM**0.44 !- Initial rime density
2202 DUM=MIN(900., MAX(170., DUM) ) !- Final rime density, kg m^-3
2203 RFrime=PI*DUM*DUM1 !- Rime factor for the density of rimed ice
2204 ELSE
2205 RFrime=1. !- Homogeneous freezing of cloud water does not
2206 !- modify RF, contributes to cloud ice
2207 ENDIF
2208 !
2209 RFnew=(PSDEP+RFrime*PIACWI+RFmx*PIACR)/DELS
2210 DUM2=TOT_ICE+THICK*DELS
2211 DUM3=RimeF1*TOT_ICE+RFnew*THICK*DELS
2212 IF (DUM2<=CLIMIT) THEN
2213 RimeF=RimeF1
2214 ELSE
2215 RimeF=MIN(RFmx, MAX(1., DUM3/DUM2) )
2216 ENDIF
2217 !apr22 - end; see (40) above
2218 ENDIF ! End IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ)
2219 !
2220 DUM1=BLDTRH*(FLIMASS*VSNOW+(1.-FLIMASS)*VCI)
2221 QInew=TOT_ICEnew/(THICK+DUM1)
2222 IF (QInew .LE. EPSQ) QInew=0.
2223 IF (QI.GT.0. .AND. QInew.NE.0.) THEN
2224 DUM=QInew/QI
2225 IF (DUM .LT. TOLER) QInew=0.
2226 ENDIF
2227 ASNOWnew=QInew*DUM1
2228 IF (ASNOW.GT.0. .AND. ASNOWnew.NE.0.) THEN
2229 DUM=ASNOWnew/ASNOW
2230 IF (DUM .LT. TOLER) ASNOWnew=0.
2231 ENDIF
2232 RQSnew=FLIMASS*RHO*QInew
2233 IF (RQSnew>0.) THEN !jul28 begin
2234 !mar03 - Use mean diameter for snow/graupel (INDEXS) from above
2235 Nsnow=RQSnew/(RimeF*MASSI(INDEXS))
2236 IF (RQSnew>0.0025 .AND. RimeF>5.) THEN
2237 !mar12 - Blend to Nsnow=1.e3 for 2.5 g m^-3 < RQSnew <= 5 g m^-3 when RF>5
2238 IF (RQSnew<=0.005) THEN
2239 DUM=MIN(1., 400.*RQSnew-1.) !- Blend at 2.5 < RQSnew (g m^-3) < 5
2240 Nsnow=1.E3*DUM+Nsnow*(1.-DUM) !- Final, blended Nsnow
2241 ELSE
2242 !-- When RF>5, this branch will produce 56.1, 62.2, & 69.1 dBZ reflectivities
2243 ! when RQSnew reaches 5, 7.5, and 10 g m^-3 (respectively) due to
2244 ! Nsnow being reduced to 1, 0.55, and 0.2 L^-1 (respectively). A 60 dBZ is reached
2245 ! when RQSnew=6.6 g m^-3 & Nsnow=0.712 L^-1 (solving a quadratic eqn).
2246 DUM=180.*(RQSnew-0.005) !- Steadily decrease Nsnow at > 5 g m^-3
2247 Nsnow=1.E3*(1.-MIN(DUM,0.8)) !- from 1 L^-1 down to 0.2 L^-1 at >=10 g m^-3
2248 ENDIF
2249 ENDIF
2250 Zsnow=Cdry*RQSnew*RQSnew/Nsnow
2251 ENDIF !jul28 end
2252 ENDIF ! End IF (TOT_ICEnew .LE. CLIMIT)
2253 ENDIF ! End IF (ICE_logical)
2254 !
2255 !--- Update rain mixing ratios
2256 !
2257 !---
2258 ! * TOT_RAINnew - total mass of rain after microphysics
2259 ! current layer and the input flux of ice from above
2260 ! * VRAIN2 - time-averaged fall speed of rain in grid and below
2261 ! * QRdum - first-guess estimate (dummy) rain mixing ratio in layer
2262 ! (uses rain fall speed from grid box above, VRAIN1)
2263 ! * QRnew - updated rain mixing ratio in layer
2264 ! -> TOT_RAINnew=QRnew*(THICK+BLDTRH*VRAIN2)
2265 ! * ARAINnew - updated accumulation of rain at bottom of grid cell
2266 !---
2267 !
2268 TIMLT=TIMLT+PIMLT*THICK !may31
2269 DELR=PRAUT+PRACW+PIACWR-PIACR+PIMLT+PREVP+PICND
2270 TOT_RAINnew=TOT_RAIN+THICK*DELR
2271 IF (TOT_RAIN.GT.0. .AND. TOT_RAINnew.NE.0.) THEN
2272 DUM=TOT_RAINnew/TOT_RAIN
2273 IF (DUM .LT. TOLER) TOT_RAINnew=0.
2274 ENDIF
2275 update_rn: IF (TOT_RAINnew .LE. CLIMIT) THEN
2276 TOT_RAINnew=0.
2277 VRAIN2=0.
2278 QRnew=0.
2279 ARAINnew=0.
2280 ELSE update_rn
2281 !may11 - start
2282 QRnew=TOT_RAINnew/(THICK+BLDTRH*VRAIN2)
2283 RQRnew=RHO*QRnew
2284 !may20 - start - made slight changes to speed up algorithm, plus remove
2285 ! block that increases N0r for heavy rain.
2286 IF(STRAT .AND. RQRnew>1.E-3) STRAT=.FALSE. !may31
2287 IF(DRZL .AND. TIMLT>EPSQ) DRZL=.FALSE. !jun01
2288 !
2289 !-- Iterate for consistent QRnew, RQRnew, N0r, INDEXR2, VRAIN2
2290 !
2291 INDEXR2=INDEXR1
2292 rain_pass: DO IPASS=1,3
2293 DSD2: IF (RQRnew<=RQR_DRmin) THEN
2294 !-- Extremely light rain
2295 INDEXR=MDRmin
2296 N0r=RQRnew/MASSR(INDEXR)
2297 ELSE IF (DRZL) THEN DSD2
2298 !-- Drizzle - gradually increase N0r for rain contents <0.5 g m^-3 !may31
2299 DUM=MAX(1.0, 0.5E-3/RQRnew)
2300 N0r=MIN(1.E9, N0r0*DUM*SQRT(DUM) ) !- 8.e6 <= N0r <= 1.e9
2301 DUM1=RQRnew/(PI*RHOL*N0r)
2302 DUM=1.E6*SQRT(SQRT(DUM1))
2303 INDEXR=MAX(XMRmin, MIN(DUM, XMRmax) )
2304 N0r=RQRnew/MASSR(INDEXR)
2305 ELSE IF (RQRnew>=RQR_DRmax) THEN
2306 !-- Extremely heavy rain
2307 INDEXR=MDRmax
2308 N0r=RQRnew/MASSR(INDEXR)
2309 ELSE DSD2
2310 !-- Light to heavy rain
2311 N0r=N0r0
2312 DUM=CN0r0*SQRT(SQRT(RQRnew))
2313 INDEXR=MAX(XMRmin, MIN(DUM, XMRmax) )
2314 IF (STRAT .AND. INDEXR<INDEXRmin) THEN
2315 !-- Stratiform rain
2316 INDEXR=INDEXRmin
2317 N0r=RQRnew/MASSR(INDEXR)
2318 ENDIF
2319 ENDIF DSD2
2320 VRAIN2=GAMMAR*VRAIN(INDEXR)
2321 QRnew=TOT_RAINnew/(THICK+BLDTRH*VRAIN2)
2322 RQRnew=RHO*QRnew
2323 IDR=ABS(INDEXR-INDEXR2)
2324 INDEXR2=INDEXR !jun13
2325 IF(IDR<=5) EXIT rain_pass !- Agreement within 5 microns
2326 !may20 - end
2327 ENDDO rain_pass
2328 !may11 - end
2329 IF (QRnew .LE. EPSQ) QRnew=0.
2330 IF (QR.GT.0. .AND. QRnew.NE.0.) THEN
2331 DUM=QRnew/QR
2332 IF (DUM .LT. TOLER) QRnew=0.
2333 ENDIF
2334 ARAINnew=BLDTRH*VRAIN2*QRnew
2335 IF (ARAIN.GT.0. .AND. ARAINnew.NE.0.) THEN
2336 DUM=ARAINnew/ARAIN
2337 IF (DUM .LT. TOLER) ARAINnew=0.
2338 ENDIF
2339 RQRnew=RHO*QRnew
2340 !may11 - start
2341 IF (RQRnew>EPSQ) THEN
2342 !may20 - remove may11 change in N0r calculation
2343 Nrain=N0r*1.E-6*REAL(INDEXR2)
2344 Zrain=Crain*RQRnew*RQRnew/Nrain !-- Rain reflectivity
2345 ENDIF
2346 !may11 - end
2347 ENDIF update_rn
2348 !
2349 WCnew=QWnew+QRnew+QInew
2350 !
2351 !----------------------------------------------------------------------
2352 !-------------- Begin debugging & verification ------------------------
2353 !----------------------------------------------------------------------
2354 !
2355 !--- QT, QTnew - total water (vapor & condensate) before & after microphysics, resp.
2356 !
2357 QT=THICK*(WV+WC)+ARAIN+ASNOW
2358 QTnew=THICK*(WVnew+WCnew)+ARAINnew+ASNOWnew
2359 BUDGET=QT-QTnew
2360 !
2361 !--- Additional check on budget preservation, accounting for truncation effects
2362 !
2363 DBG_logical=.FALSE.
2364 DUM=ABS(BUDGET)
2365 IF (DUM .GT. TOLER) THEN
2366 DUM=DUM/MIN(QT, QTnew)
2367 IF (DUM .GT. TOLER) DBG_logical=.TRUE.
2368 ENDIF
2369 !
2370 IF (PRINT_diag) THEN
2371 ESW=MIN(1000.*FPVS0(Tnew),0.99*PP)
2372 QSWnew=(RHgrd+0.001)*EPS*ESW/(PP-ESW)
2373 IF( (QWnew>EPSQ .OR. QRnew>EPSQ .OR. WVnew>QSWnew) &
2374 & .AND. TC<T_ICE) DBG_logical=.TRUE.
2375 ENDIF
2376 !
2377 IF ((WVnew.LT.EPSQ .OR. DBG_logical) .AND. PRINT_diag) THEN
2378 !
2379 WRITE(0,"(/2(a,i4),2(a,i2))") '{} i=',I_index,' j=',J_index,&
2380 & ' L=',L,' LSFC=',LSFC
2381 !
2382 ESW=MIN(1000.*FPVS0(Tnew),0.99*PP)
2383 QSWnew=EPS*ESW/(PP-ESW)
2384 IF (TC.LT.0. .OR. Tnew .LT. 0.) THEN
2385 ESI=MIN(1000.*FPVS(Tnew),0.99*PP)
2386 QSInew=EPS*ESI/(PP-ESI)
2387 ELSE
2388 QSI=QSW
2389 QSInew=QSWnew
2390 ENDIF
2391 WSnew=QSInew
2392 WRITE(0,"(4(a12,g11.4,1x))") &
2393 & '{} TCold=',TC,'TCnew=',Tnew-T0C,'P=',.01*PP,'RHO=',RHO, &
2394 & '{} THICK=',THICK,'RHold=',WV/WS,'RHnew=',WVnew/WSnew, &
2395 & 'RHgrd=',RHgrd, &
2396 & '{} RHWold=',WV/QSW,'RHWnew=',WVnew/QSWnew,'RHIold=',WV/QSI, &
2397 & 'RHInew=',WVnew/QSInew, &
2398 & '{} QSWold=',QSW,'QSWnew=',QSWnew,'QSIold=',QSI,'QSInew=',QSInew, &
2399 & '{} WSold=',WS,'WSnew=',WSnew,'WVold=',WV,'WVnew=',WVnew, &
2400 & '{} WCold=',WC,'WCnew=',WCnew,'QWold=',QW,'QWnew=',QWnew, &
2401 & '{} QIold=',QI,'QInew=',QInew,'QRold=',QR,'QRnew=',QRnew, &
2402 & '{} ARAINold=',ARAIN,'ARAINnew=',ARAINnew,'ASNOWold=',ASNOW, &
2403 & 'ASNOWnew=',ASNOWnew, &
2404 & '{} TOT_RAIN=',TOT_RAIN,'TOT_RAINnew=',TOT_RAINnew, &
2405 & 'TOT_ICE=',TOT_ICE,'TOT_ICEnew=',TOT_ICEnew, &
2406 & '{} BUDGET=',BUDGET,'QTold=',QT,'QTnew=',QTnew
2407 !
2408 WRITE(0,"(4(a12,g11.4,1x))") &
2409 & '{} DELT=',DELT,'DELV=',DELV,'DELW=',DELW,'DELI=',DELI, &
2410 & '{} DELR=',DELR,'PCOND=',PCOND,'PIDEP=',PIDEP,'PIEVP=',PIEVP, &
2411 & '{} PICND=',PICND,'PREVP=',PREVP,'PRAUT=',PRAUT,'PRACW=',PRACW, &
2412 & '{} PIACW=',PIACW,'PIACWI=',PIACWI,'PIACWR=',PIACWR,'PIMLT=', &
2413 & PIMLT, &
2414 & '{} PIACR=',PIACR
2415 !
2416 WRITE(0,"(4(a13,L2))") &
2417 & '{} ICE_logical=',ICE_logical,'RAIN_logical=',RAIN_logical, &
2418 & 'STRAT=',STRAT,'DRZL=',DRZL
2419 !
2420 IF (ICE_logical) WRITE(0,"(4(a12,g11.4,1x))") &
2421 & '{} RimeF1=',RimeF1,'GAMMAS=',GAMMAS,'VrimeF=',VrimeF, &
2422 & 'VSNOW=',VSNOW, &
2423 & '{} QSmICE=',QSmICE,'XLIMASS=',XLIMASS,'RQLICE=',RQLICE, &
2424 & 'QTICE=',QTICE, &
2425 & '{} NLICE=',NLICE,'NSmICE=',NSmICE,'PILOSS=',PILOSS, &
2426 & 'EMAIRI=',EMAIRI, &
2427 & '{} RimeF=',RimeF,'VCI=',VCI,'INDEXS=',FLOAT(INDEXS), &
2428 & 'FLIMASS=',FLIMASS, &
2429 & '{} Nsnow=',Nsnow,'Zsnow=',Zsnow,'TIMLT=',TIMLT,'RHOX0C=',RHOX0C, &
2430 & '{} INDEXS0C=',FLOAT(INDEXS0C)
2431 !
2432 IF (RAIN_logical) WRITE(0,"(4(a12,g11.4,1x))") &
2433 & '{} INDEXR1=',FLOAT(INDEXR1),'INDEXR=',FLOAT(INDEXR), &
2434 & 'GAMMAR=',GAMMAR,'N0r=',N0r, &
2435 & '{} VRAIN1=',VRAIN1,'VRAIN2=',VRAIN2,'QTRAIN=',QTRAIN,'RQR=',RQR, &
2436 & '{} PRLOSS=',PRLOSS,'VOLR1=',THICK+BLDTRH*VRAIN1, &
2437 & 'VOLR2=',THICK+BLDTRH*VRAIN2,'INDEXR2=',INDEXR2, &
2438 & '{} Nrain=',Nrain,'Zrain=',Zrain
2439 !
2440 IF (PRACW .GT. 0.) WRITE(0,"(a12,g11.4,1x)") '{} FWR=',FWR
2441 !
2442 IF (PIACR .GT. 0.) WRITE(0,"(a12,g11.4,1x)") '{} FIR=',FIR
2443 !
2444 DUM=PIMLT+PICND-PREVP-PIEVP
2445 IF (DUM.GT.0. .or. DWVi.NE.0.) &
2446 & WRITE(0,"(4(a12,g11.4,1x))") &
2447 & '{} TFACTOR=',TFACTOR,'DYNVIS=',DYNVIS, &
2448 & 'THERM_CON=',THERM_COND,'DIFFUS=',DIFFUS
2449 !
2450 IF (PREVP .LT. 0.) WRITE(0,"(4(a12,g11.4,1x))") &
2451 & '{} RFACTOR=',RFACTOR,'ABW=',ABW,'VENTR=',VENTR,'CREVP=',CREVP, &
2452 & '{} DWVr=',DWVr,'DENOMW=',DENOMW
2453 !
2454 IF (PIDEP.NE.0. .AND. DWVi.NE.0.) &
2455 & WRITE(0,"(4(a12,g11.4,1x))") &
2456 & '{} DWVi=',DWVi,'DENOMI=',DENOMI,'PIDEP_max=',PIDEP_max, &
2457 & 'SFACTOR=',SFACTOR, &
2458 & '{} ABI=',ABI,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
2459 & 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
2460 & '{} VENTIS=',VENTIS,'DIDEP=',DIDEP
2461 !
2462 IF (PIDEP.GT.0. .AND. PCOND.NE.0.) &
2463 & WRITE(0,"(4(a12,g11.4,1x))") &
2464 & '{} DENOMW=',DENOMW,'DENOMWI=',DENOMWI,'DENOMF=',DENOMF, &
2465 & 'DUM2=',PCOND-PIACW
2466 !
2467 IF (FWS .GT. 0.) WRITE(0,"(4(a12,g11.4,1x))") &
2468 & '{} FWS=',FWS
2469 !
2470 DUM=PIMLT+PICND-PIEVP
2471 IF (DUM.GT. 0.) WRITE(0,"(4(a12,g11.4,1x))") &
2472 & '{} SFACTOR=',SFACTOR,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
2473 & 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
2474 & '{} AIEVP=',AIEVP,'DIEVP=',DIEVP,'QSW0=',QSW0,'DWV0=',DWV0
2475 !
2476 ENDIF
2477 !
2478 !----------------------------------------------------------------------
2479 !------------------------- Update arrays ------------------------------
2480 !----------------------------------------------------------------------
2481 !
2482 T_col(L)=Tnew ! Update temperature
2483 Q_col(L)=max(EPSQ, WVnew/(1.+WVnew)) ! Update specific humidity
2484 WC_col(L)=max(EPSQ, WCnew) ! Update total condensate mixing ratio
2485 QI_col(L)=max(EPSQ, QInew) ! Update ice mixing ratio
2486 QR_col(L)=max(EPSQ, QRnew) ! Update rain mixing ratio
2487 QW_col(L)=max(EPSQ, QWnew) ! Update cloud water mixing ratio
2488 RimeF_col(L)=RimeF ! Update rime factor
2489 NR_col(L)=Nrain ! Update rain number concentration !jul28 begin
2490 NS_col(L)=Nsnow ! Update snow number concentration
2491 !aligo IF (RQRnew>RQRmix .AND. RQSnew>RQSmix .and. RimeF>10.) THEN
2492 ! Zsnow=Zsnow*Cwet !apr27
2493 !aligo ENDIF
2494 Ztot=MAX(Ztot, Zrain+Zsnow)
2495 Ztot=Zrain+Zsnow
2496 IF (Ztot>Zmin) DBZ_col(L)=10.*ALOG10(Ztot) ! Update radar reflectivity
2497 ASNOW=ASNOWnew ! Update accumulated snow
2498 VSNOW1=VSNOW ! Update total ice fall speed
2499 ARAIN=ARAINnew ! Update accumulated rain
2500 VRAIN1=VRAIN2 ! Update rain fall speed
2501 !
2502 !-- Assign microphysical processes and fall speeds to 1D array
2503 !
2504 if(pcond.gt.0)then
2505 pcond1d(L)=pcond
2506 elseif(pcond.lt.0)then
2507 pevap1d(L)=pcond
2508 endif
2509 if(pidep.gt.0)then
2510 pidep1d(L)=pidep
2511 elseif(pidep.lt.0)then
2512 pisub1d(L)=pidep
2513 endif
2514 piacw1d(L)=piacw
2515 piacwi1d(L)=piacwi
2516 piacwr1d(L)=piacwr
2517 piacr1d(L)=piacr
2518 picnd1d(L)=picnd
2519 pievp1d(L)=pievp
2520 pimlt1d(L)=pimlt
2521 praut1d(L)=praut
2522 pracw1d(L)=pracw
2523 prevp1d(L)=prevp
2524 if (qinew>EPSQ) then
2525 vsnow1d(L)=vsnow
2526 !sep22a - Start changes
2527 if (FLIMASS<1.) then
2528 vci1d(L)=vci
2529 NSmICE1d(L)=NSmICE
2530 endif
2531 !sep22a - End changes
2532 INDEXS1d(L)=FLOAT(INDEXS)
2533 endif
2534 if (qrnew>EPSQ) then
2535 vrain11d(L)=vrain1
2536 vrain21d(L)=vrain2
2537 INDEXR1d(L)=FLOAT(INDEXR2)
2538 !jun01 - start
2539 IDR=0
2540 IF(STRAT) IDR=1
2541 IF(DRZL) IDR=IDR+2
2542 RFlag1d(L)=IDR
2543 !jun01 - end
2544 endif
2545 !
2546 !#######################################################################
2547 !
2548 ENDDO big_loop
2549 !
2550 !#######################################################################
2551 !
2552 !-----------------------------------------------------------------------
2553 !--------------------------- Return to GSMDRIVE -----------------------
2554 !-----------------------------------------------------------------------
2555 !
2556 CONTAINS
2557 !#######################################################################
2558 !--------- Produces accurate calculation of cloud condensation ---------
2559 !#######################################################################
2560 !
2561 REAL FUNCTION CONDENSE (PP,QW,TK,WV,RHgrd,I,J,L)
2562 !
2563 !---------------------------------------------------------------------------------
2564 !------ The Asai (1965) algorithm takes into consideration the release of ------
2565 !------ latent heat in increasing the temperature & in increasing the ------
2566 !------ saturation mixing ratio (following the Clausius-Clapeyron eqn.). ------
2567 !---------------------------------------------------------------------------------
2568 !
2569 IMPLICIT NONE
2570 !
2571 INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
2572 REAL (KIND=HIGH_PRES), PARAMETER :: &
2573 & RHLIMIT=.001, RHLIMIT1=-RHLIMIT
2574 REAL (KIND=HIGH_PRES) :: COND, SSAT, WCdum
2575 !
2576 REAL,INTENT(IN) :: QW,PP,WV,TK,RHgrd
2577 REAL WVdum,Tdum,DWV,WS,ESW
2579 integer,INTENT(IN) :: I,J,L !-- Debug 20120111
2580 integer :: niter
2581 real :: DWVinp,SSATinp
2583 !
2584 !-----------------------------------------------------------------------
2585 !
2586 Tdum=TK
2587 WVdum=WV
2588 WCdum=QW
2589 ESW=MIN(1000.*FPVS0(Tdum),0.99*PP) ! Saturation vapor press w/r/t water
2590 WS=RHgrd*EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
2591 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
2592 SSAT=DWV/WS ! Supersaturation ratio
2593 CONDENSE=0.
2595 DWVinp=DWV !-- Debug 20120111
2596 SSATinp=SSAT
2598 DO NITER=1,MAX_ITERATIONS
2599 COND=DWV/(1.+XLV3*WS/(Tdum*Tdum)) ! Asai (1965, J. Japan)
2600 COND=MAX(COND, -WCdum) ! Limit cloud water evaporation
2601 Tdum=Tdum+XLV1*COND ! Updated temperature
2602 WVdum=WVdum-COND ! Updated water vapor mixing ratio
2603 WCdum=WCdum+COND ! Updated cloud water mixing ratio
2604 CONDENSE=CONDENSE+COND ! Total cloud water condensation
2605 ESW=MIN(1000.*FPVS0(Tdum),0.99*PP) ! Updated saturation vapor press w/r/t water
2606 WS=RHgrd*EPS*ESW/(PP-ESW) ! Updated saturation mixing ratio w/r/t water
2607 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
2608 SSAT=DWV/WS ! Grid-scale supersaturation ratio
2609 IF (SSAT>=RHLIMIT1 .AND. SSAT<=RHLIMIT) EXIT !-- Exit if SSAT is near 0
2610 IF (SSAT<RHLIMIT1 .AND. WCdum<=EPSQ) EXIT !-- Exit if SSAT<0 & no cloud water
2611 ENDDO
2613 !-- Debug 20120111
2614 IF (NITER>MAX_ITERATIONS) THEN
2615 !-- Too many iterations - indicates possible numerical instability
2616 IF (WARN3) THEN
2617 write(0,*) 'WARN3: Too many iterations in function CONDENSE; ', &
2618 ' I,J,L,TC,SSAT,QW,DWV=',I,J,L,TK-T0C,SSATinp,1000.*QW,DWVinp
2619 WARN3=.FALSE.
2620 ENDIF
2621 SSAT=0.
2622 CONDENSE=DWVinp
2623 ENDIF
2625 !
2626 END FUNCTION CONDENSE
2627 !
2628 !#######################################################################
2629 !---------------- Calculate ice deposition at T<T_ICE ------------------
2630 !#######################################################################
2631 !
2632 REAL FUNCTION DEPOSIT (PP,Tdum,WVdum,RHgrd,I,J,L) !-- Debug 20120111
2633 !
2634 !--- Also uses the Asai (1965) algorithm, but uses a different target
2635 ! vapor pressure for the adjustment
2636 !
2637 IMPLICIT NONE
2638 !
2639 INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
2640 REAL (KIND=HIGH_PRES), PARAMETER :: RHLIMIT=.001, &
2641 & RHLIMIT1=-RHLIMIT
2642 REAL (KIND=HIGH_PRES) :: DEP, SSAT
2643 !
2644 real,INTENT(IN) :: PP,RHgrd
2645 real,INTENT(INOUT) :: WVdum,Tdum
2646 real ESI,WS,DWV
2648 integer,INTENT(IN) :: I,J,L !-- Debug 20120111
2649 integer :: niter
2650 real :: Tinp,DWVinp,SSATinp
2652 !
2653 !-----------------------------------------------------------------------
2654 !
2655 ESI=MIN(1000.*FPVS(Tdum),0.99*PP) ! Saturation vapor press w/r/t ice
2656 WS=RHgrd*EPS*ESI/(PP-ESI) ! Saturation mixing ratio
2657 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
2658 SSAT=DWV/WS ! Supersaturation ratio
2659 DEPOSIT=0.
2661 Tinp=Tdum !-- Debug 20120111
2662 DWVinp=DWV
2663 SSATinp=SSAT
2665 DO NITER=1,MAX_ITERATIONS
2666 DEP=DWV/(1.+XLS3*WS/(Tdum*Tdum)) ! Asai (1965, J. Japan)
2667 Tdum=Tdum+XLS1*DEP ! Updated temperature
2668 WVdum=WVdum-DEP ! Updated ice mixing ratio
2669 DEPOSIT=DEPOSIT+DEP ! Total ice deposition
2670 ESI=MIN(1000.*FPVS(Tdum),0.99*PP) ! Updated saturation vapor press w/r/t ice
2671 WS=RHgrd*EPS*ESI/(PP-ESI) ! Updated saturation mixing ratio w/r/t ice
2672 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
2673 SSAT=DWV/WS ! Grid-scale supersaturation ratio
2674 IF (SSAT>=RHLIMIT1 .AND. SSAT<=RHLIMIT) EXIT !-- Exit if SSAT is near 0
2675 ENDDO
2677 !-- Debug 20120111
2678 IF (NITER>MAX_ITERATIONS) THEN
2679 !-- Too many iterations - indicates possible numerical instability
2680 IF (WARN2) THEN
2681 write(0,*) 'WARN2: Too many iterations in function DEPOSIT; ', &
2682 ' I,J,L,TC,SSAT,DWV=',I,J,L,Tinp-T0C,SSATinp,DWVinp
2683 WARN2=.FALSE.
2684 ENDIF
2685 SSAT=0.
2686 DEPOSIT=DWVinp
2687 ENDIF
2689 !
2690 END FUNCTION DEPOSIT
2691 !
2692 !#######################################################################
2693 !--- Used to calculate the mean size of rain drops (INDEXR) in microns
2694 !#######################################################################
2695 !
2696 !-- NOTE: this function is not used in this version.
2697 !
2698 INTEGER FUNCTION GET_INDEXR(RR)
2699 IMPLICIT NONE
2700 REAL, INTENT(IN) :: RR
2701 IF (RR .LE. RR_DRmin) THEN
2702 !
2703 !--- Assume fixed mean diameter of rain (0.05 mm) for low rain rates,
2704 ! instead vary N0r with rain rate
2705 !
2706 GET_INDEXR=MDRmin
2707 ELSE IF (RR .LE. RR_DR1) THEN
2708 !
2709 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2710 ! for mean diameters (Dr) between 0.05 and 0.10 mm:
2711 ! V(Dr)=5.6023e4*Dr**1.136, V in m/s and Dr in m
2712 ! RR = PI*1000.*N0r0*5.6023e4*Dr**(4+1.136) = 1.408e15*Dr**5.136,
2713 ! RR in kg/(m**2*s)
2714 ! Dr (m) = 1.123e-3*RR**.1947 -> Dr (microns) = 1.123e3*RR**.1947
2715 !
2716 GET_INDEXR=INT( 1.123E3*RR**.1947 + .5 )
2717 GET_INDEXR=MAX( MDRmin, MIN(GET_INDEXR, MDR1) )
2718 ELSE IF (RR .LE. RR_DR2) THEN
2719 !
2720 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2721 ! for mean diameters (Dr) between 0.10 and 0.20 mm:
2722 ! V(Dr)=1.0867e4*Dr**.958, V in m/s and Dr in m
2723 ! RR = PI*1000.*N0r0*1.0867e4*Dr**(4+.958) = 2.731e14*Dr**4.958,
2724 ! RR in kg/(m**2*s)
2725 ! Dr (m) = 1.225e-3*RR**.2017 -> Dr (microns) = 1.225e3*RR**.2017
2726 !
2727 GET_INDEXR=INT( 1.225E3*RR**.2017 + .5 )
2728 GET_INDEXR=MAX( MDR1, MIN(GET_INDEXR, MDR2) )
2729 ELSE IF (RR .LE. RR_DR3) THEN
2730 !
2731 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2732 ! for mean diameters (Dr) between 0.20 and 0.32 mm:
2733 ! V(Dr)=2831.*Dr**.80, V in m/s and Dr in m
2734 ! RR = PI*1000.*N0r0*2831.*Dr**(4+.80) = 7.115e13*Dr**4.80,
2735 ! RR in kg/(m**2*s)
2736 ! Dr (m) = 1.3006e-3*RR**.2083 -> Dr (microns) = 1.3006e3*RR**.2083
2737 !
2738 GET_INDEXR=INT( 1.3006E3*RR**.2083 + .5 )
2739 GET_INDEXR=MAX( MDR2, MIN(GET_INDEXR, MDR3) )
2740 ELSE IF (RR .LE. RR_DR4) THEN
2741 !
2742 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2743 ! for mean diameters (Dr) between 0.32 and 0.45 mm:
2744 ! V(Dr)=963.0*Dr**.666, V in m/s and Dr in m
2745 ! RR = PI*1000.*N0r0*963.0*Dr**(4+.666) = 2.4205e13*Dr**4.666,
2746 ! RR in kg/(m**2*s)
2747 ! Dr (m) = 1.354e-3*RR**.2143 -> Dr (microns) = 1.354e3*RR**.2143
2748 !
2749 GET_INDEXR=INT( 1.354E3*RR**.2143 + .5 )
2750 GET_INDEXR=MAX( MDR3, MIN(GET_INDEXR, MDR4) )
2751 ELSE IF (RR .LE. RR_DR5) THEN
2752 !
2753 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2754 ! for mean diameters (Dr) between 0.45 and 0.675 mm:
2755 ! V(Dr)=309.0*Dr**.5185, V in m/s and Dr in m
2756 ! RR = PI*1000.*N0r0*309.0*Dr**(4+.5185) = 7.766e12*Dr**4.5185,
2757 ! RR in kg/(m**2*s)
2758 ! Dr (m) = 1.404e-3*RR**.2213 -> Dr (microns) = 1.404e3*RR**.2213
2759 !
2760 GET_INDEXR=INT( 1.404E3*RR**.2213 + .5 )
2761 GET_INDEXR=MAX( MDR4, MIN(GET_INDEXR, MDR5) )
2762 ELSE IF (RR .LE. RR_DRmax) THEN
2763 !
2764 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2765 ! for mean diameters (Dr) between 0.675 and 1.0 mm:
2766 ! V(Dr)=85.81Dr**.343, V in m/s and Dr in m
2767 ! RR = PI*1000.*N0r0*85.81*Dr**(4+.343) = 2.1566e12*Dr**4.343,
2768 ! RR in kg/(m**2*s)
2769 ! Dr (m) = 1.4457e-3*RR**.2303 -> Dr (microns) = 1.4457e3*RR**.2303
2770 !
2771 GET_INDEXR=INT( 1.4457E3*RR**.2303 + .5 )
2772 GET_INDEXR=MAX( MDR5, MIN(GET_INDEXR, MDRmax) )
2773 ELSE
2774 !
2775 !--- Assume fixed mean diameter of rain (1.0 mm) for high rain rates,
2776 ! instead vary N0r with rain rate
2777 !
2778 GET_INDEXR=MDRmax
2779 ENDIF ! End IF (RR .LE. RR_DRmin) etc.
2780 !
2781 END FUNCTION GET_INDEXR
2782 !
2783 END SUBROUTINE EGCP01COLUMN
2785 !#######################################################################
2786 !------- Initialize constants & lookup tables for microphysics ---------
2787 !#######################################################################
2788 !
2790 ! SH 0211/2002
2792 !-----------------------------------------------------------------------
2793 SUBROUTINE fer_hires_init (GSMDT,DT,DELX,DELY,LOWLYR,restart, &
2794 !HWRF & MP_RESTART_STATE,TBPVS_STATE,TBPVS0_STATE, &
2795 & ALLOWED_TO_READ, &
2796 & IDS,IDE,JDS,JDE,KDS,KDE, &
2797 & IMS,IME,JMS,JME,KMS,KME, &
2798 & ITS,ITE,JTS,JTE,KTS,KTE, &
2799 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY)
2800 !-----------------------------------------------------------------------
2801 !-------------------------------------------------------------------------------
2802 !--- SUBPROGRAM DOCUMENTATION BLOCK
2803 ! PRGRMMR: Ferrier ORG: W/NP22 DATE: February 2001
2804 ! Jin ORG: W/NP22 DATE: January 2002
2805 ! (Modification for WRF structure)
2806 !
2807 !-------------------------------------------------------------------------------
2808 ! ABSTRACT:
2809 ! * Reads various microphysical lookup tables used in COLUMN_MICRO
2810 ! * Lookup tables were created "offline" and are read in during execution
2811 ! * Creates lookup tables for saturation vapor pressure w/r/t water & ice
2812 !-------------------------------------------------------------------------------
2813 !
2814 ! USAGE: CALL fer_hires_init FROM SUBROUTINE GSMDRIVE AT MODEL START TIME
2815 !
2816 ! INPUT ARGUMENT LIST:
2817 ! DTPH - physics time step (s)
2818 !
2819 ! OUTPUT ARGUMENT LIST:
2820 ! NONE
2821 !
2822 ! OUTPUT FILES:
2823 ! NONE
2824 !
2825 ! SUBROUTINES:
2826 ! MY_GROWTH_RATES - lookup table for growth of nucleated ice
2827 ! GPVS - lookup tables for saturation vapor pressure (water, ice)
2828 !
2829 ! UNIQUE: NONE
2830 !
2831 ! LIBRARY: NONE
2832 !
2833 ! COMMON BLOCKS:
2834 ! CMICRO_CONS - constants used in GSMCOLUMN
2835 ! CMY600 - lookup table for growth of ice crystals in
2836 ! water saturated conditions (Miller & Young, 1979)
2837 ! IVENT_TABLES - lookup tables for ventilation effects of ice
2838 ! IACCR_TABLES - lookup tables for accretion rates of ice
2839 ! IMASS_TABLES - lookup tables for mass content of ice
2840 ! IRATE_TABLES - lookup tables for precipitation rates of ice
2841 ! IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
2842 ! MAPOT - Need lat/lon grid resolution
2843 ! RVENT_TABLES - lookup tables for ventilation effects of rain
2844 ! RACCR_TABLES - lookup tables for accretion rates of rain
2845 ! RMASS_TABLES - lookup tables for mass content of rain
2846 ! RVELR_TABLES - lookup tables for fall speeds of rain
2847 ! RRATE_TABLES - lookup tables for precipitation rates of rain
2848 !
2849 ! ATTRIBUTES:
2850 ! LANGUAGE: FORTRAN 90
2851 ! MACHINE : IBM SP
2852 !
2853 !-----------------------------------------------------------------------
2854 !
2855 !
2856 !-----------------------------------------------------------------------
2857 IMPLICIT NONE
2858 !-----------------------------------------------------------------------
2859 !-------------------------------------------------------------------------
2860 !-------------- Parameters & arrays for lookup tables --------------------
2861 !-------------------------------------------------------------------------
2862 !
2863 !--- Common block of constants used in column microphysics
2864 !
2865 !WRF
2866 ! real DLMD,DPHD
2867 !WRF
2868 !
2869 !-----------------------------------------------------------------------
2870 !--- Parameters & data statement for local calculations
2871 !-----------------------------------------------------------------------
2872 !
2873 INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3
2874 !
2875 ! VARIABLES PASSED IN
2876 integer,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
2877 & ,IMS,IME,JMS,JME,KMS,KME &
2878 & ,ITS,ITE,JTS,JTE,KTS,KTE
2879 !WRF
2880 INTEGER, DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: LOWLYR
2881 !
2882 real, INTENT(IN) :: DELX,DELY
2883 !HWRF real,DIMENSION(*), INTENT(INOUT) :: MP_RESTART_STATE
2884 !HWRF real,DIMENSION(NX), INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
2885 real,DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(OUT),OPTIONAL :: &
2886 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
2887 real,INTENT(IN) :: DT,GSMDT
2888 LOGICAL,INTENT(IN) :: allowed_to_read,restart
2889 !
2890 !-----------------------------------------------------------------------
2891 ! LOCAL VARIABLES
2892 !-----------------------------------------------------------------------
2893 REAL :: BBFR,DTPH,DX,Thour_print,RDIS
2894 INTEGER :: I,IM,J,L,K,JTF,KTF,ITF
2895 INTEGER :: etampnew_unit1
2896 LOGICAL :: opened
2897 LOGICAL , EXTERNAL :: wrf_dm_on_monitor
2898 CHARACTER*80 errmess
2899 !
2900 !-----------------------------------------------------------------------
2901 !
2902 JTF=MIN0(JTE,JDE-1)
2903 KTF=MIN0(KTE,KDE-1)
2904 ITF=MIN0(ITE,IDE-1)
2905 !
2906 DO J=JTS,JTF
2907 DO I=ITS,ITF
2908 LOWLYR(I,J)=1
2909 ENDDO
2910 ENDDO
2911 !
2912 IF(.NOT.RESTART .AND. ALLOWED_TO_READ .AND. present(F_ICE_PHY)) THEN !HWRF
2913 CALL wrf_debug(1,'WARNING: F_ICE_PHY,F_RAIN_PHY AND F_RIMEF_PHY IS REINITIALIZED') !HWRF
2914 DO J = jts,jte
2915 DO K = kts,kte
2916 DO I= its,ite
2917 F_ICE_PHY(i,k,j)=0.
2918 F_RAIN_PHY(i,k,j)=0.
2919 F_RIMEF_PHY(i,k,j)=1.
2920 ENDDO
2921 ENDDO
2922 ENDDO
2923 ENDIF
2924 !
2925 !-----------------------------------------------------------------------
2926 IF(ALLOWED_TO_READ)THEN
2927 !-----------------------------------------------------------------------
2928 !
2929 DX=SQRT((DELX)**2+(DELY)**2)/1000. ! Model resolution at equator (km) !GFDL
2930 DX=MIN(100., MAX(5., DX) )
2931 !
2932 !-- Relative humidity threshold for the onset of grid-scale condensation
2933 !!-- 9/1/01: Assume the following functional dependence for 5 - 100 km resolution:
2934 !! RHgrd=0.90 for dx=100 km, 0.98 for dx=5 km, where
2935 ! RHgrd=0.90+.08*((100.-DX)/95.)**.5
2936 !
2937 DTPH=MAX(GSMDT*60.,DT)
2938 DTPH=NINT(DTPH/DT)*DT
2939 !
2940 !--- Create lookup tables for saturation vapor pressure w/r/t water & ice
2941 !
2942 CALL GPVS
2943 !
2944 !--- Read in various lookup tables
2945 !
2946 IF ( wrf_dm_on_monitor() ) THEN
2947 DO i = 31,99
2948 INQUIRE ( i , OPENED = opened )
2949 IF ( .NOT. opened ) THEN
2950 etampnew_unit1 = i
2951 GOTO 2061
2952 ENDIF
2953 ENDDO
2954 etampnew_unit1 = -1
2955 2061 CONTINUE
2956 ENDIF
2957 !
2958 CALL wrf_dm_bcast_bytes ( etampnew_unit1 , IWORDSIZE )
2959 !
2960 IF ( etampnew_unit1 < 0 ) THEN
2961 CALL wrf_error_fatal ( 'module_mp_fer_hires: fer_hires_init: Can not find unused fortran unit to read in lookup table.' )
2962 ENDIF
2963 !
2964 IF ( wrf_dm_on_monitor() ) THEN
2965 !!was OPEN (UNIT=1,FILE="eta_micro_lookup.dat",FORM="UNFORMATTED")
2966 !-- Use the new, expanded tables for rain (maximum mean diameter of 1 mm rather than 0.45 mm)
2967 !old OPEN(UNIT=etampnew_unit1,FILE="ETAMPNEW_DATA", &
2968 print*,'open ETAMPNEW_DATA.expanded_rain, in fer_hires'
2969 OPEN(UNIT=etampnew_unit1,FILE="ETAMPNEW_DATA.expanded_rain", &
2970 & FORM="UNFORMATTED",STATUS="OLD",ERR=9061)
2971 !
2972 READ(etampnew_unit1) VENTR1
2973 READ(etampnew_unit1) VENTR2
2974 READ(etampnew_unit1) ACCRR
2975 READ(etampnew_unit1) MASSR
2976 READ(etampnew_unit1) VRAIN
2977 READ(etampnew_unit1) RRATE
2978 READ(etampnew_unit1) VENTI1
2979 READ(etampnew_unit1) VENTI2
2980 READ(etampnew_unit1) ACCRI
2981 READ(etampnew_unit1) MASSI
2982 READ(etampnew_unit1) VSNOWI
2983 READ(etampnew_unit1) VEL_RF
2984 ! read(etampnew_unit1) my_growth ! Applicable only for DTPH=180 s
2985 CLOSE (etampnew_unit1)
2986 ENDIF
2987 !
2988 CALL wrf_dm_bcast_bytes ( VENTR1 , size ( VENTR1 ) * RWORDSIZE )
2989 CALL wrf_dm_bcast_bytes ( VENTR2 , size ( VENTR2 ) * RWORDSIZE )
2990 CALL wrf_dm_bcast_bytes ( ACCRR , size ( ACCRR ) * RWORDSIZE )
2991 CALL wrf_dm_bcast_bytes ( MASSR , size ( MASSR ) * RWORDSIZE )
2992 CALL wrf_dm_bcast_bytes ( VRAIN , size ( VRAIN ) * RWORDSIZE )
2993 CALL wrf_dm_bcast_bytes ( RRATE , size ( RRATE ) * RWORDSIZE )
2994 CALL wrf_dm_bcast_bytes ( VENTI1 , size ( VENTI1 ) * RWORDSIZE )
2995 CALL wrf_dm_bcast_bytes ( VENTI2 , size ( VENTI2 ) * RWORDSIZE )
2996 CALL wrf_dm_bcast_bytes ( ACCRI , size ( ACCRI ) * RWORDSIZE )
2997 CALL wrf_dm_bcast_bytes ( MASSI , size ( MASSI ) * RWORDSIZE )
2998 CALL wrf_dm_bcast_bytes ( VSNOWI , size ( VSNOWI ) * RWORDSIZE )
2999 CALL wrf_dm_bcast_bytes ( VEL_RF , size ( VEL_RF ) * RWORDSIZE )
3000 !
3001 !--- Calculates coefficients for growth rates of ice nucleated in water
3002 ! saturated conditions, scaled by physics time step (lookup table)
3003 !
3004 CALL MY_GROWTH_RATES (DTPH)
3005 !
3006 PI_E=ACOS(-1.)
3007 !
3008 !--- Constants associated with Biggs (1953) freezing of rain, as parameterized
3009 ! following Lin et al. (JCAM, 1983) & Reisner et al. (1998, QJRMS).
3010 !
3011 ABFR=-0.66
3012 BBFR=100.
3013 CBFR=20.*PI_E*PI_E*BBFR*RHOL*1.E-21*DTPH !mar03 - Bug fix in (33); include DTPH
3014 !
3015 !--- CIACW is used in calculating riming rates
3016 ! The assumed effective collection efficiency of cloud water rimed onto
3017 ! ice is =0.5 below:
3018 !
3019 CIACW=0.5*DTPH*0.25*PI_E !jul28 (16)
3020 !
3021 !--- CIACR is used in calculating freezing of rain colliding with large ice
3022 ! The assumed collection efficiency is 1.0
3023 !
3024 CIACR=PI_E*DTPH
3025 !
3026 !--- Based on rain lookup tables for mean diameters from 0.05 to 0.45 mm
3027 ! * Four different functional relationships of mean drop diameter as
3028 ! a function of rain rate (RR), derived based on simple fits to
3029 ! mass-weighted fall speeds of rain as functions of mean diameter
3030 ! from the lookup tables.
3031 !
3032 RR_DRmin=N0r0*RRATE(MDRmin) ! RR for mean drop diameter of .05 mm
3033 RR_DR1=N0r0*RRATE(MDR1) ! RR for mean drop diameter of .10 mm
3034 RR_DR2=N0r0*RRATE(MDR2) ! RR for mean drop diameter of .20 mm
3035 RR_DR3=N0r0*RRATE(MDR3) ! RR for mean drop diameter of .32 mm
3036 RR_DR4=N0r0*RRATE(MDR4) ! RR for mean drop diameter of .45 mm
3037 RR_DR5=N0r0*RRATE(MDR5) ! RR for mean drop diameter of .675 mm
3038 RR_DRmax=N0r0*RRATE(MDRmax) ! RR for mean drop diameter of 1.0 mm
3039 !
3040 RQR_DRmin=N0r0*MASSR(MDRmin) ! Rain content for mean drop diameter of .05 mm
3041 RQR_DRmax=N0r0*MASSR(MDRmax) ! Rain content for mean drop diameter of 1.0 mm
3042 C_NR=1./(PI*RHOL) !jul28
3043 Crain=720.E18*C_NR*C_NR !jul28
3044 C_N0r0=PI_E*RHOL*N0r0
3045 CN0r0=1.E6/SQRT(SQRT(C_N0r0))
3046 CN0r_DMRmin=1./(PI_E*RHOL*DMRmin*DMRmin*DMRmin*DMRmin)
3047 CN0r_DMRmax=1./(PI_E*RHOL*DMRmax*DMRmax*DMRmax*DMRmax)
3048 !
3049 !--- CRACW is used in calculating collection of cloud water by rain (an
3050 ! assumed collection efficiency of 1.0)
3051 !
3052 CRACW=DTPH*0.25*PI_E*1.0
3053 !
3054 ESW0=1000.*FPVS0(T0C) ! Saturation vapor pressure at 0C
3055 RFmax=1.1**Nrime ! Maximum rime factor allowed
3056 RFmx1=PI*900. ! Density near that of pure ice, 900 kg m^-3
3057 !
3058 !------------------------------------------------------------------------
3059 !--------------- Constants passed through argument list -----------------
3060 !------------------------------------------------------------------------
3061 !
3062 !--- Important parameters for self collection (autoconversion) of
3063 ! cloud water to rain.
3064 !
3065 !-- Relative dispersion == standard deviation of droplet spectrum / mean radius
3066 ! (see pp 1542-1543, Liu & Daum, JAS, 2004)
3067 RDIS=0.5 !-- relative dispersion of droplet spectrum
3068 BETA6=( (1.+3.*RDIS*RDIS)*(1.+4.*RDIS*RDIS)*(1.+5.*RDIS*RDIS)/ &
3069 & ((1.+RDIS*RDIS)*(1.+2.*RDIS*RDIS) ) )
3070 !-- Kappa=1.1e10 g^-2 cm^3 s^-1 after eq. (8b) on p.1105 of Liu et al. (JAS, 2006)
3071 ! => More extensive units conversion than can be described here to go from
3072 ! eq. (13) in Liu et al. (JAS, 2006) to what's programmed below. Note that
3073 ! the units used throughout the paper are in cgs units!
3074 !
3075 ARAUT=1.03e19/(NCW*SQRT(NCW))
3076 BRAUT=DTPH*1.1E10*BETA6/NCW
3077 !
3078 !--- QAUT0 is the *OLD* threshold cloud content for autoconversion to rain
3079 ! needed for droplets to reach a diameter of 20 microns (following
3080 ! Manton and Cotton, 1977; Banta and Hanson, 1987, JCAM). It's no longer
3081 ! used in this version, but the value is passed into radiation in case
3082 ! a ball park estimate is needed.
3083 !--- QAUT0=1.2567, 0.8378, or 0.4189 g/m**3 for droplet number concentrations
3084 ! of 300, 200, and 100 cm**-3, respectively
3085 !
3086 QAUT0=PI*RHOL*NCW*(20.E-6)**3/6. !-- legacy
3087 !
3088 !--- For calculating cloud droplet radius in microns, Rcw
3089 !
3090 ARcw=1.E6*(0.75/(PI*NCW*RHOL))**C1
3091 !
3092 !may20 - start
3093 !
3094 !- An explanation for the following settings assumed for "hail" or frozen drops (RF>10):
3095 ! RH_NgC=PI*500.*5.E3
3096 ! RHOg=500 kg m^-3, Ng=5.e3 m^-3 => "hail" content >7.85 g m^-3 for INDEXS=MDImax
3097 !- or -
3098 ! RH_NgC=PI*500.*1.E3
3099 ! RHOg=500 kg m^-3, Ng=1.e3 m^-3 => "hail" content >1.57 g m^-3 for INDEXS=MDImax
3100 !
3101 RH_NgC=PI*500.*1.E3 !- RHOg=500 kg m^-3, Ng=1.e3 m^-3
3102 RQhail=RH_NgC*(1.E-3)**3 !- Threshold "hail" content ! becomes 1.57 g m^-3
3103 Bvhail=0.82*C1 !- Exponent for Thompson graupel
3104 Avhail=1353.*(1./RH_NgC)**Bvhail !- 1353 ~ constant for Thompson graupel
3105 RH_NgC=1.E6*(1./RH_NgC)**C1 !mar08 (convection, convert to microns)
3106 !
3107 !- An explanation for the following settings assumed for graupel (RF>5):
3108 ! RH_NgT=PI*300.*10.E3
3109 ! RHOg=300 kg m^-3, Ng=10.e3 m^-3 => "graupel" content must exceed 9.43 g m^-3 for INDEXS=MDImax
3110 !- or -
3111 ! RH_NgT=PI*300.*5.E3
3112 ! RHOg=300 kg m^-3, Ng=5.e3 m^-3 => "graupel" content must exceed 4.71 g m^-3 for INDEXS=MDImax
3113 !
3114 RH_NgT=PI*300.*5.E3 !- RHOg=300 kg m^-3, Ng=5.e3 m^-3
3115 RH_NgT=1.E6*(1./RH_NgT)**C1 !mar08 (transition, convert to microns)
3116 !may20 - end
3117 !
3118 !--- For calculating snow optical depths by considering bulk density of
3119 ! snow based on emails from Q. Fu (6/27-28/01), where optical
3120 ! depth (T) = 1.5*SWP/(Reff*DENS), SWP is snow water path, Reff
3121 ! is effective radius, and DENS is the bulk density of snow.
3122 !
3123 ! SWP (kg/m**2)=(1.E-3 kg/g)*SWPrad, SWPrad in g/m**2 used in radiation
3124 ! T = 1.5*1.E3*SWPrad/(Reff*DENS)
3125 !
3126 ! See derivation for MASSI(INDEXS), note equal to RHO*QSNOW/NSNOW
3127 !
3128 ! SDENS=1.5e3/DENS, DENS=MASSI(INDEXS)/[PI_E*(1.E-6*INDEXS)**3]
3129 !
3130 DO I=MDImin,MDImax
3131 SDENS(I)=PI_E*1.5E-15*FLOAT(I*I*I)/MASSI(I)
3132 ENDDO
3133 !
3134 Thour_print=-DTPH/3600.
3137 ENDIF ! Allowed_to_read
3139 RETURN
3140 !
3141 !-----------------------------------------------------------------------
3142 !
3143 9061 CONTINUE
3144 WRITE( errmess , '(A,I4)' ) &
3145 'module_mp_hwrf: error opening ETAMPNEW_DATA on unit ' &
3146 &, etampnew_unit1
3147 CALL wrf_error_fatal(errmess)
3148 !
3149 !-----------------------------------------------------------------------
3150 END SUBROUTINE fer_hires_init
3151 !
3152 SUBROUTINE MY_GROWTH_RATES (DTPH)
3153 !
3154 !--- Below are tabulated values for the predicted mass of ice crystals
3155 ! after 600 s of growth in water saturated conditions, based on
3156 ! calculations from Miller and Young (JAS, 1979). These values are
3157 ! crudely estimated from tabulated curves at 600 s from Fig. 6.9 of
3158 ! Young (1993). Values at temperatures colder than -27C were
3159 ! assumed to be invariant with temperature.
3160 !
3161 !--- Used to normalize Miller & Young (1979) calculations of ice growth
3162 ! over large time steps using their tabulated values at 600 s.
3163 ! Assumes 3D growth with time**1.5 following eq. (6.3) in Young (1993).
3164 !
3165 IMPLICIT NONE
3166 !
3167 REAL,INTENT(IN) :: DTPH
3168 !
3169 REAL DT_ICE
3170 REAL,DIMENSION(35) :: MY_600
3171 !WRF
3172 !
3173 !-----------------------------------------------------------------------
3174 DATA MY_600 / &
3175 & 5.5e-8, 1.4E-7, 2.8E-7, 6.E-7, 3.3E-6, &
3176 & 2.E-6, 9.E-7, 8.8E-7, 8.2E-7, 9.4e-7, &
3177 & 1.2E-6, 1.85E-6, 5.5E-6, 1.5E-5, 1.7E-5, &
3178 & 1.5E-5, 1.E-5, 3.4E-6, 1.85E-6, 1.35E-6, &
3179 & 1.05E-6, 1.E-6, 9.5E-7, 9.0E-7, 9.5E-7, &
3180 & 9.5E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7, &
3181 & 9.E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7 / ! -31 to -35 deg C
3182 !
3183 !-----------------------------------------------------------------------
3184 !
3185 DT_ICE=(DTPH/600.)**1.5
3186 MY_GROWTH=DT_ICE*MY_600*1.E-3 !-- 20090714: Convert from g to kg
3187 !
3188 !-----------------------------------------------------------------------
3189 !
3190 END SUBROUTINE MY_GROWTH_RATES
3191 !
3192 !-----------------------------------------------------------------------
3193 !--------- Old GFS saturation vapor pressure lookup tables -----------
3194 !-----------------------------------------------------------------------
3195 !
3196 SUBROUTINE GPVS
3197 ! ******************************************************************
3198 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
3199 ! . . .
3200 ! SUBPROGRAM: GPVS COMPUTE SATURATION VAPOR PRESSURE TABLE
3201 ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
3202 !
3203 ! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE TABLE AS A FUNCTION OF
3204 ! TEMPERATURE FOR THE TABLE LOOKUP FUNCTION FPVS.
3205 ! EXACT SATURATION VAPOR PRESSURES ARE CALCULATED IN SUBPROGRAM FPVSX.
3206 ! THE CURRENT IMPLEMENTATION COMPUTES A TABLE WITH A LENGTH
3207 ! OF 7501 FOR TEMPERATURES RANGING FROM 180.0 TO 330.0 KELVIN.
3208 !
3209 ! PROGRAM HISTORY LOG:
3210 ! 91-05-07 IREDELL
3211 ! 94-12-30 IREDELL EXPAND TABLE
3212 ! 96-02-19 HONG ICE EFFECT
3213 ! 01-11-29 JIN MODIFIED FOR WRF
3214 !
3215 ! USAGE: CALL GPVS
3216 !
3217 ! SUBPROGRAMS CALLED:
3218 ! (FPVSX) - INLINABLE FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE
3219 !
3220 ! COMMON BLOCKS:
3221 ! COMPVS - SCALING PARAMETERS AND TABLE FOR FUNCTION FPVS.
3222 !
3223 ! ATTRIBUTES:
3224 ! LANGUAGE: FORTRAN 90
3225 !
3226 !$$$
3227 IMPLICIT NONE
3228 real :: X,XINC,T
3229 integer :: JX
3230 !----------------------------------------------------------------------
3231 XINC=(XMAX-XMIN)/(NX-1)
3232 C1XPVS=1.-XMIN/XINC
3233 C2XPVS=1./XINC
3234 C1XPVS0=1.-XMIN/XINC
3235 C2XPVS0=1./XINC
3236 !
3237 DO JX=1,NX
3238 X=XMIN+(JX-1)*XINC
3239 T=X
3240 TBPVS(JX)=FPVSX(T)
3241 TBPVS0(JX)=FPVSX0(T)
3242 ENDDO
3243 !
3244 END SUBROUTINE GPVS
3245 !-----------------------------------------------------------------------
3246 !***********************************************************************
3247 !-----------------------------------------------------------------------
3248 REAL FUNCTION FPVS(T)
3249 !-----------------------------------------------------------------------
3250 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
3251 ! . . .
3252 ! SUBPROGRAM: FPVS COMPUTE SATURATION VAPOR PRESSURE
3253 ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
3254 !
3255 ! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE FROM THE TEMPERATURE.
3256 ! A LINEAR INTERPOLATION IS DONE BETWEEN VALUES IN A LOOKUP TABLE
3257 ! COMPUTED IN GPVS. SEE DOCUMENTATION FOR FPVSX FOR DETAILS.
3258 ! INPUT VALUES OUTSIDE TABLE RANGE ARE RESET TO TABLE EXTREMA.
3259 ! THE INTERPOLATION ACCURACY IS ALMOST 6 DECIMAL PLACES.
3260 ! ON THE CRAY, FPVS IS ABOUT 4 TIMES FASTER THAN EXACT CALCULATION.
3261 ! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
3262 !
3263 ! PROGRAM HISTORY LOG:
3264 ! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
3265 ! 94-12-30 IREDELL EXPAND TABLE
3266 ! 96-02-19 HONG ICE EFFECT
3267 ! 01-11-29 JIN MODIFIED FOR WRF
3268 !
3269 ! USAGE: PVS=FPVS(T)
3270 !
3271 ! INPUT ARGUMENT LIST:
3272 ! T - REAL TEMPERATURE IN KELVIN
3273 !
3274 ! OUTPUT ARGUMENT LIST:
3275 ! FPVS - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
3276 !
3277 ! ATTRIBUTES:
3278 ! LANGUAGE: FORTRAN 90
3279 !$$$
3280 IMPLICIT NONE
3281 real,INTENT(IN) :: T
3282 real XJ
3283 integer :: JX
3284 !-----------------------------------------------------------------------
3285 IF (T>=XMIN .AND. T<=XMAX) THEN
3286 XJ=MIN(MAX(C1XPVS+C2XPVS*T,1.),FLOAT(NX))
3287 JX=MIN(XJ,NX-1.)
3288 FPVS=TBPVS(JX)+(XJ-JX)*(TBPVS(JX+1)-TBPVS(JX))
3289 ELSE IF (T>XMAX) THEN
3290 !-- Magnus Tetens formula for water saturation (Murray, 1967)
3291 ! (saturation vapor pressure in kPa)
3292 FPVS=0.61078*exp(17.2694*(T-273.16)/(T-35.86))
3293 ELSE
3294 !-- Magnus Tetens formula for ice saturation(Murray, 1967)
3295 ! (saturation vapor pressure in kPa)
3296 FPVS=0.61078*exp(21.8746*(T-273.16)/(T-7.66))
3297 ENDIF
3298 !
3299 END FUNCTION FPVS
3300 !-----------------------------------------------------------------------
3301 !-----------------------------------------------------------------------
3302 REAL FUNCTION FPVS0(T)
3303 !-----------------------------------------------------------------------
3304 IMPLICIT NONE
3305 real,INTENT(IN) :: T
3306 real :: XJ1
3307 integer :: JX1
3308 !-----------------------------------------------------------------------
3309 IF (T>=XMIN .AND. T<=XMAX) THEN
3310 XJ1=MIN(MAX(C1XPVS0+C2XPVS0*T,1.),FLOAT(NX))
3311 JX1=MIN(XJ1,NX-1.)
3312 FPVS0=TBPVS0(JX1)+(XJ1-JX1)*(TBPVS0(JX1+1)-TBPVS0(JX1))
3313 ELSE
3314 !-- Magnus Tetens formula for water saturation (Murray, 1967)
3315 ! (saturation vapor pressure in kPa)
3316 FPVS0=0.61078*exp(17.2694*(T-273.16)/(T-35.86))
3317 ENDIF
3318 !
3319 END FUNCTION FPVS0
3320 !-----------------------------------------------------------------------
3321 !***********************************************************************
3322 !-----------------------------------------------------------------------
3323 REAL FUNCTION FPVSX(T)
3324 !-----------------------------------------------------------------------
3325 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
3326 ! . . .
3327 ! SUBPROGRAM: FPVSX COMPUTE SATURATION VAPOR PRESSURE
3328 ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
3329 !
3330 ! ABSTRACT: EXACTLY COMPUTE SATURATION VAPOR PRESSURE FROM TEMPERATURE.
3331 ! THE WATER MODEL ASSUMES A PERFECT GAS, CONSTANT SPECIFIC HEATS
3332 ! FOR GAS AND LIQUID, AND NEGLECTS THE VOLUME OF THE LIQUID.
3333 ! THE MODEL DOES ACCOUNT FOR THE VARIATION OF THE LATENT HEAT
3334 ! OF CONDENSATION WITH TEMPERATURE. THE ICE OPTION IS NOT INCLUDED.
3335 ! THE CLAUSIUS-CLAPEYRON EQUATION IS INTEGRATED FROM THE TRIPLE POINT
3336 ! TO GET THE FORMULA
3337 ! PVS=PSATK*(TR**XA)*EXP(XB*(1.-TR))
3338 ! WHERE TR IS TTP/T AND OTHER VALUES ARE PHYSICAL CONSTANTS
3339 ! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
3340 !
3341 ! PROGRAM HISTORY LOG:
3342 ! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
3343 ! 94-12-30 IREDELL EXACT COMPUTATION
3344 ! 96-02-19 HONG ICE EFFECT
3345 ! 01-11-29 JIN MODIFIED FOR WRF
3346 !
3347 ! USAGE: PVS=FPVSX(T)
3348 ! REFERENCE: EMANUEL(1994),116-117
3349 !
3350 ! INPUT ARGUMENT LIST:
3351 ! T - REAL TEMPERATURE IN KELVIN
3352 !
3353 ! OUTPUT ARGUMENT LIST:
3354 ! FPVSX - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
3355 !
3356 ! ATTRIBUTES:
3357 ! LANGUAGE: FORTRAN 90
3358 !$$$
3359 IMPLICIT NONE
3360 !-----------------------------------------------------------------------
3361 real, parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
3362 , CLIQ=4.1855E+3,CVAP= 1.8460E+3 &
3363 , CICE=2.1060E+3,HSUB=2.8340E+6
3364 !
3365 real, parameter :: PSATK=PSAT*1.E-3
3366 real, parameter :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
3367 real, parameter :: DLDTI=CVAP-CICE &
3368 , XAI=-DLDTI/RV,XBI=XAI+HSUB/(RV*TTP)
3369 real T,TR
3370 !-----------------------------------------------------------------------
3371 TR=TTP/T
3372 !
3373 IF(T.GE.TTP)THEN
3374 FPVSX=PSATK*(TR**XA)*EXP(XB*(1.-TR))
3375 ELSE
3376 FPVSX=PSATK*(TR**XAI)*EXP(XBI*(1.-TR))
3377 ENDIF
3378 !
3379 END FUNCTION FPVSX
3380 !-----------------------------------------------------------------------
3381 !-----------------------------------------------------------------------
3382 REAL FUNCTION FPVSX0(T)
3383 !-----------------------------------------------------------------------
3384 IMPLICIT NONE
3385 real,parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
3386 , CLIQ=4.1855E+3,CVAP=1.8460E+3 &
3387 , CICE=2.1060E+3,HSUB=2.8340E+6
3388 real,PARAMETER :: PSATK=PSAT*1.E-3
3389 real,PARAMETER :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
3390 real,PARAMETER :: DLDTI=CVAP-CICE &
3391 , XAI=-DLDT/RV,XBI=XA+HSUB/(RV*TTP)
3392 real :: T,TR
3393 !-----------------------------------------------------------------------
3394 TR=TTP/T
3395 FPVSX0=PSATK*(TR**XA)*EXP(XB*(1.-TR))
3396 !
3397 END FUNCTION FPVSX0
3398 !
3399 END MODULE module_mp_fer_hires