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