da_cld_eff_radius.inc
References to this file elsewhere.
1 subroutine da_cld_eff_radius(t,rho,qci,qrn,qsn,qgr,snow,xice,xland,method, &
2 reff_water,reff_ice,reff_rain,reff_snow,reff_grau)
3
4 !---------------------------------------------------------------------------
5 ! Purpose: compute effective radius of cloud liquid water and cloud ice
6 !
7 ! METHOD:
8 ! liquid water: adapted from WRFV2/phys/module_ra_cam.F, analytic formula following
9 ! the formulation originally developed by J. T. Kiehl, and
10 ! ice (method 1): adapted from WRFV2/phys/module_ra_cam.F, Kristjansson and Mitchell
11 ! ice (method 2): WSM3, Hong et al., MWR 2004
12 ! rain/snow/graupel: assume exponential particle size distribution and
13 ! spherical particles
14 ! use Gauss-Laguerre Quadrature for integration
15 !---------------------------------------------------------------------------
16
17 integer, intent(in) :: method
18 real, intent(in) :: t ! temperature
19 real, intent(in) :: rho ! air density kg/m3
20 real, intent(in) :: qci ! cloud ice mixing ratio kg/kg
21 real, intent(in) :: qrn ! cloud rain mixing ratio
22 real, intent(in) :: qsn ! cloud snow mixing ratio
23 real, intent(in) :: qgr ! cloud graupel mixing ratio
24 real, intent(in) :: snow ! snow water equivalent
25 real, intent(in) :: xice ! ice percentage
26 real, intent(in) :: xland ! landsea percentage
27 real, intent(out) :: reff_water ! effective radius liquid water
28 real, intent(out) :: reff_ice ! effective radius ice
29 real, intent(out) :: reff_rain ! effective radius rain
30 real, intent(out) :: reff_snow ! effective radius snow
31 real, intent(out) :: reff_grau ! effective radius graupel
32
33 ! local variables
34 integer :: index, nk
35 real :: snowh, corr
36 real, parameter :: rliqland = 8.0 ! liquid drop size if over land
37 real, parameter :: rliqocean = 14.0 ! liquid drop size if over ocean
38 real, parameter :: rliqice = 14.0 ! liquid drop size if over sea ice
39 ! Tabulated values of re(T) in the temperature interval
40 ! 180 K -- 274 K; hexagonal columns assumed:
41 real, dimension(95), parameter :: retab = &
42 (/ 5.92779, 6.26422, 6.61973, 6.99539, 7.39234, &
43 7.81177, 8.25496, 8.72323, 9.21800, 9.74075, 10.2930, &
44 10.8765, 11.4929, 12.1440, 12.8317, 13.5581, 14.2319, &
45 15.0351, 15.8799, 16.7674, 17.6986, 18.6744, 19.6955, &
46 20.7623, 21.8757, 23.0364, 24.2452, 25.5034, 26.8125, &
47 27.7895, 28.6450, 29.4167, 30.1088, 30.7306, 31.2943, &
48 31.8151, 32.3077, 32.7870, 33.2657, 33.7540, 34.2601, &
49 34.7892, 35.3442, 35.9255, 36.5316, 37.1602, 37.8078, &
50 38.4720, 39.1508, 39.8442, 40.5552, 41.2912, 42.0635, &
51 42.8876, 43.7863, 44.7853, 45.9170, 47.2165, 48.7221, &
52 50.4710, 52.4980, 54.8315, 57.4898, 60.4785, 63.7898, &
53 65.5604, 71.2885, 75.4113, 79.7368, 84.2351, 88.8833, &
54 93.6658, 98.5739, 103.603, 108.752, 114.025, 119.424, &
55 124.954, 130.630, 136.457, 142.446, 148.608, 154.956, &
56 161.503, 168.262, 175.248, 182.473, 189.952, 197.699, &
57 205.728, 214.055, 222.694, 231.661, 240.971, 250.639 /)
58 !
59 ! constants from da_control.f: pi,
60 !
61 real, parameter :: n0_rain = 0.08 ! cm(-4)
62 real, parameter :: n0_snow = 0.04 ! cm(-4)
63 real, parameter :: n0_grau = 0.04 ! cm(-4)
64 real, parameter :: rho_w = 1000.0 ! kg m(-3)
65 real, parameter :: rho_ice = 900.0 ! kg m(-3)
66 real, parameter :: rho_snow = 100.0 ! kg m(-3)
67 real, parameter :: rho_grau = 400.0 ! kg m(-3)
68
69 ! Abscissas of Gauss-Laguerre Integration
70 real, dimension(32) :: xk = (/ 0.0444893658333, 0.23452610952, &
71 0.576884629302, 1.07244875382, 1.72240877644, 2.52833670643, &
72 3.49221327285, 4.61645677223, 5.90395848335, 7.3581268086, &
73 8.98294126732, 10.783012089, 12.763745476, 14.9309117981, &
74 17.2932661372, 19.8536236493, 22.6357789624, 25.6201482024, &
75 28.8739336869, 32.3333294017, 36.1132042245, 40.1337377056, &
76 44.5224085362, 49.2086605665, 54.3501813324, 59.8791192845, &
77 65.9833617041, 72.6842683222, 80.1883747906, 88.735192639, &
78 98.8295523184, 111.751398227 /)
79
80 ! total weights (weight*exp(xk)) of Gauss-Laguerre Integration
81 real, dimension(32) :: totalw = (/ 0.114187105768, 0.266065216898, &
82 0.418793137325, 0.572532846497, 0.727648788453, 0.884536718946, &
83 1.04361887597, 1.20534920595, 1.37022171969, 1.53877595906, &
84 1.71164594592, 1.8895649683, 2.07318851235, 2.26590144444, &
85 2.46997418988, 2.64296709494, 2.76464437462, 3.22890542981, &
86 2.92019361963, 4.3928479809, 4.27908673189, 5.20480398519, &
87 5.11436212961, 4.15561492173, 6.19851060567, 5.34795780128, &
88 6.28339212457, 6.89198340969, 7.92091094244, 9.20440555803, &
89 11.1637432904, 15.3902417688 /)
90
91 real, parameter :: limit = 1.0e-6
92 real :: piover6 ! pi/6
93 real :: sum1_rain, sum2_rain, sum1_snow, sum2_snow, sum1_grau, sum2_grau
94 real :: lamda_rain, lamda_snow, lamda_grau
95 real, dimension(32) :: psd_rain, psd_snow, psd_grau ! partical size distribution
96
97 ! initialize
98 reff_water = 0.0
99 reff_ice = 0.0
100 reff_rain = 0.0
101 reff_snow = 0.0
102 reff_grau = 0.0
103 lamda_rain = 0.0
104 lamda_snow = 0.0
105 lamda_grau = 0.0
106
107 ! cloud liquid effective radius
108
109 snowh = 0.001 * snow ! here the snowh (in meter) is water-equivalent snow depth,
110 ! which is different from the actually snow depth defined in the wrf output file
111
112 ! Start with temperature-dependent value appropriate for continental air
113 reff_water = rliqland + (rliqocean-rliqland) * min(1.0,max(0.0,(t_triple-t)*0.05))
114 ! Modify for snow depth over land
115 reff_water = reff_water + (rliqocean-reff_water) * min(1.0,max(0.0,snowh*10.))
116 ! Ramp between polluted value over land to clean value over ocean.
117 reff_water = reff_water + (rliqocean-reff_water) * min(1.0,max(0.0,1.0-xland))
118 ! Ramp between the resultant value and a sea ice value in the presence of ice.
119 reff_water = reff_water + (rliqice-reff_water) * min(1.0,max(0.0,xice))
120
121 reff_water = reff_water * 0.001 ! micron to mm
122 write(555,*) 'reff_water=',reff_water
123
124 ! cloud ice effective radius
125
126 if ( method == 1 ) then
127 index = int(t-179.)
128 index = min(max(index,1),94)
129 corr = t - int(t)
130 reff_ice = retab(index)*(1.-corr) + retab(index+1)*corr
131 end if
132
133 ! cloud ice effective radius
134 ! rho*qci = 2.08*10**22 * Dice**8
135
136 if ( method == 2 ) then
137 ! 0.5 for diameter - radii conversion
138 ! 1.0e6 for meter - micron conversion
139 ! 0.125 = 1/8
140 reff_ice = 1.0e6 * 0.5 * ( rho * qci / 2.08e22 ) ** 0.125
141 end if
142
143 reff_ice = reff_ice * 0.001 ! micron to mm
144 write(555,*) 'reff_ice=',reff_ice
145 !
146 ! cloud rain/snow/graupel effective radius
147 !
148 piover6 = pi/6.
149 if ( qrn > limit ) then
150 lamda_rain = (piover6*rho_w*n0_rain*rho/qrn)**0.25
151 end if
152 if ( qsn > limit ) then
153 lamda_snow = (piover6*rho_snow*n0_snow*rho/qsn)**0.25
154 end if
155 if ( qgr > limit ) then
156 lamda_grau = (piover6*rho_grau*n0_grau*rho/qgr)**0.25
157 end if
158 sum1_rain = 0.0
159 sum2_rain = 0.0
160 sum1_snow = 0.0
161 sum2_snow = 0.0
162 sum1_grau = 0.0
163 sum2_grau = 0.0
164 if ( qrn > limit ) then
165 do nk = 1, 32
166 psd_rain(nk) = n0_rain*exp(-2.0*lamda_rain*xk(nk))
167 sum1_rain = sum1_rain + totalw(nk) * (xk(nk)**3) * psd_rain(nk)
168 sum2_rain = sum2_rain + totalw(nk) * (xk(nk)**2) * psd_rain(nk)
169 end do
170 reff_rain = 10. * sum1_rain/sum2_rain ! mm
171 write(555,*) 'reff_rain=',reff_rain
172 end if
173 if ( qsn > limit ) then
174 do nk = 1, 32
175 psd_snow(nk) = n0_snow*exp(-2.0*lamda_snow*xk(nk))
176 sum1_snow = sum1_snow + totalw(nk) * (xk(nk)**3) * psd_snow(nk)
177 sum2_snow = sum2_snow + totalw(nk) * (xk(nk)**2) * psd_snow(nk)
178 end do
179 reff_snow = 10. * sum1_snow/sum2_snow ! mm
180 write(555,*) 'reff_snow=',reff_snow
181 end if
182 if ( qgr > limit ) then
183 do nk = 1, 32
184 psd_grau(nk) = n0_grau*exp(-2.0*lamda_grau*xk(nk))
185 sum1_grau = sum1_grau + totalw(nk) * (xk(nk)**3) * psd_grau(nk)
186 sum2_grau = sum2_grau + totalw(nk) * (xk(nk)**2) * psd_grau(nk)
187 end do
188 reff_grau = 10. * sum1_grau/sum2_grau ! mm
189 write(555,*) 'reff_grau=',reff_grau
190 end if
191
192 end subroutine da_cld_eff_radius