PartMC  2.6.1
coag_kernel_sedi.F90
Go to the documentation of this file.
1 ! Copyright (C) 2005-2011 Nicole Riemer and Matthew West
2 ! Copyright (C) Andreas Bott
3 ! Licensed under the GNU General Public License version 2 or (at your
4 ! option) any later version. See the file COPYING for details.
5 
6 !> \file
7 !> The pmc_coag_kernel_sedi module.
8 !!
9 !! Contains code based on \c coad1d.f by Andreas Bott
10 !! - http://www.meteo.uni-bonn.de/mitarbeiter/ABott/
11 !! - Released under the GPL to Nicole Riemer (personal communication)
12 !! - A. Bott, A flux method for the numerical solution of the
13 !! stochastic collection equation, J. Atmos. Sci. 55, 2284-2293,
14 !! 1998.
15 
16 !> Gravitational sedimentation coagulation kernel.
17 module pmc_coag_kernel_sedi
18 
19  use pmc_env_state
20  use pmc_constants
21  use pmc_aero_data
22  use pmc_aero_particle
23 
24 contains
25 
26 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
27 
28  !> Sedimentation coagulation kernel.
29  subroutine kernel_sedi(aero_particle_1, aero_particle_2, &
30  aero_data, env_state, k)
31 
32  !> First particle.
33  type(aero_particle_t), intent(in) :: aero_particle_1
34  !> Second particle.
35  type(aero_particle_t), intent(in) :: aero_particle_2
36  !> Aerosol data.
37  type(aero_data_t), intent(in) :: aero_data
38  !> Environment state.
39  type(env_state_t), intent(in) :: env_state
40  !> Kernel \c k(a,b) (m^3/s).
41  real(kind=dp), intent(out) :: k
42 
43  call kernel_sedi_helper(aero_particle_volume(aero_particle_1), &
44  aero_particle_volume(aero_particle_2), aero_data, env_state%temp, &
45  env_state%pressure, k)
46 
47  end subroutine kernel_sedi
48 
49 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
50 
51  !> Minimum and maximum values of the sedimentation coagulation.
52  subroutine kernel_sedi_minmax(v1, v2, aero_data, env_state, k_min, k_max)
53 
54  !> Volume of first particle (m^3).
55  real(kind=dp), intent(in) :: v1
56  !> Volume of second particle (m^3).
57  real(kind=dp), intent(in) :: v2
58  !> Aerosol data.
59  type(aero_data_t), intent(in) :: aero_data
60  !> Environment state.
61  type(env_state_t), intent(in) :: env_state
62  !> Minimum kernel \c k(a,b) (m^3/s).
63  real(kind=dp), intent(out) :: k_min
64  !> Maximum kernel \c k(a,b) (m^3/s).
65  real(kind=dp), intent(out) :: k_max
66 
67  call kernel_sedi_helper(v1, v2, aero_data, env_state%temp, &
68  env_state%pressure, k_min)
69  k_max = k_min
70 
71  end subroutine kernel_sedi_minmax
72 
73 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
74 
75  !> Helper function that does the actual sedimentation kernel computation.
76  !!
77  !! Helper function. Do not call directly. Instead use kernel_sedi().
78  subroutine kernel_sedi_helper(v1, v2, aero_data, temp, pressure, k)
79 
80  !> Volume of first particle (m^3).
81  real(kind=dp), intent(in) :: v1
82  !> Volume of second particle (m^3).
83  real(kind=dp), intent(in) :: v2
84  !> Aerosol data.
85  type(aero_data_t), intent(in) :: aero_data
86  !> Temperature (K).
87  real(kind=dp), intent(in) :: temp
88  !> Pressure (Pa).
89  real(kind=dp), intent(in) :: pressure
90  !> Kernel k(a,b) (m^3/s).
91  real(kind=dp), intent(out) :: k
92 
93  real(kind=dp) r1, r2, winf1, winf2, ec
94 
95  r1 = aero_data_vol2rad(aero_data, v1) ! m
96  r2 = aero_data_vol2rad(aero_data, v2) ! m
97  call fall_g(aero_data_vol_to_mobility_rad(aero_data, v1, temp, pressure), &
98  winf1) ! winf1 in m/s
99  call fall_g(aero_data_vol_to_mobility_rad(aero_data, v2, temp, pressure), &
100  winf2) ! winf2 in m/s
101  call effic(r1, r2, ec) ! ec is dimensionless
102  k = ec * const%pi * (r1 + r2)**2 * abs(winf1 - winf2)
103 
104  end subroutine kernel_sedi_helper
105 
106 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
107 
108  !> Finds the terminal velocity of a particle based on its size.
109  subroutine fall_g(r, w_inf)
110 
111  !> Particle mobility radius (m).
112  real(kind=dp), intent(in) :: r
113  !> Terminal velocity (m/s).
114  real(kind=dp), intent(out) :: w_inf
115 
116  ! terminal velocity of falling drops
117  real(kind=dp) eta, xlamb, rhow, rhoa, grav, cunh, t0, sigma
118  real(kind=dp) stok, stb, phy, py, rr, x, y, xrey, bond
119  integer i
120  real(kind=dp) b(7),c(6)
121  data b /-0.318657d1,0.992696d0,-0.153193d-2,-0.987059d-3, &
122  -0.578878d-3,0.855176d-4,-0.327815d-5/
123  data c /-0.500015d1,0.523778d1,-0.204914d1,0.475294d0, &
124  -0.542819d-1,0.238449d-2/
125 
126  eta = 1.818d-4
127  xlamb = 6.62d-6
128  rhow = 1d0
129  rhoa = 1.225d-3
130  grav = 980.665d0
131  cunh = 1.257d0 * xlamb
132  t0 = 273.15d0
133  sigma = 76.1d0 - 0.155d0 * (293.15d0 - t0)
134  stok = 2d0 * grav * (rhow - rhoa) / (9d0 * eta)
135  stb = 32d0 * rhoa * (rhow - rhoa) * grav / (3d0 * eta * eta)
136  phy = sigma * sigma * sigma * rhoa * rhoa &
137  / (eta**4 * grav * (rhow - rhoa))
138  py = phy**(1d0/6d0)
139 
140  ! rr: radius in cm-units
141  rr = r * 1d2
142 
143  if (rr .le. 1d-3) then
144  w_inf = stok * (rr * rr + cunh * rr)
145  elseif (rr .gt. 1d-3 .and. rr .le. 5.35d-2) then
146  x = log(stb * rr * rr * rr)
147  y = 0d0
148  do i = 1,7
149  y = y + b(i) * (x**(i - 1))
150  end do
151  xrey = (1d0 + cunh/rr) * exp(y)
152  w_inf = xrey * eta / (2d0 * rhoa * rr)
153  elseif (rr .gt. 5.35d-2) then
154  bond = grav * (rhow - rhoa) * rr**2 / sigma
155  if (rr .gt. 0.35d0) then
156  bond = grav * (rhow - rhoa) * 0.35d0**2 / sigma
157  end if
158  x = log(16d0 * bond * py / 3d0)
159  y = 0d0
160  do i = 1,6
161  y = y + c(i) * (x**(i - 1))
162  end do
163  xrey = py * exp(y)
164  w_inf = xrey * eta / (2d0 * rhoa * rr)
165  if (rr .gt. 0.35d0) then
166  w_inf = xrey * eta / (2d0 * rhoa * 0.35d0)
167  end if
168  end if
169  w_inf = w_inf / 100d0
170 
171  end subroutine fall_g
172 
173 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
174 
175  !> Coagulation efficiency.
176  !!
177  !! Determines the chance that two particles will actually coagulate,
178  !! given that they approach close enough to do so.
179  subroutine effic(r1, r2, ec)
180 
181  !> Geometric radius of first particle (m).
182  real(kind=dp), intent(in) :: r1
183  !> Geometric radius of second particle (m).
184  real(kind=dp), intent(in) :: r2
185  !> Collision efficiency (dimensionless).
186  real(kind=dp), intent(out) :: ec
187 
188  real(kind=dp) :: r_small, r_big, rq, p, q, ek
189  integer :: k, ir, kk, iq
190  ! collision efficiencies of hall kernel
191  real(kind=dp) :: rat(21),r0(15),ecoll(15,21)
192 
193  data r0 /6.0d0,8.0d0,10.0d0,15.0d0,20.0d0,25.0d0,30.0d0,40.0d0 &
194  ,50.0d0,60.0d0,70.0d0,100.0d0,150.0d0,200.0d0,300.0d0/
195  data rat /0.0d0,0.05d0,0.1d0,0.15d0,0.2d0,0.25d0,0.3d0,0.35d0 &
196  ,0.4d0,0.45d0,0.5d0,0.55d0,0.6d0,0.65d0,0.7d0,0.75d0,0.8d0 &
197  ,0.85d0,0.9d0,0.95d0,1.0d0/
198  ! two-dimensional linear interpolation of the collision efficiency
199  data ecoll /0.001d0,0.001d0,0.001d0,0.001d0,0.001d0,0.001d0 &
200  ,0.001d0,0.001d0,0.001d0,0.001d0 ,0.001d0,0.001d0,0.001d0 &
201  ,0.001d0,0.001d0,0.003d0,0.003d0,0.003d0,0.004d0,0.005d0 &
202  ,0.005d0,0.005d0,0.010d0,0.100d0,0.050d0,0.200d0,0.500d0 &
203  ,0.770d0,0.870d0,0.970d0 ,0.007d0,0.007d0,0.007d0,0.008d0 &
204  ,0.009d0,0.010d0,0.010d0,0.070d0,0.400d0,0.430d0 ,0.580d0 &
205  ,0.790d0,0.930d0,0.960d0,1.000d0,0.009d0,0.009d0,0.009d0 &
206  ,0.012d0,0.015d0 ,0.010d0,0.020d0,0.280d0,0.600d0,0.640d0 &
207  ,0.750d0,0.910d0,0.970d0,0.980d0,1.000d0 ,0.014d0,0.014d0 &
208  ,0.014d0,0.015d0,0.016d0,0.030d0,0.060d0,0.500d0,0.700d0 &
209  ,0.770d0 ,0.840d0,0.950d0,0.970d0,1.000d0,1.000d0,0.017d0 &
210  ,0.017d0,0.017d0,0.020d0,0.022d0 ,0.060d0,0.100d0,0.620d0 &
211  ,0.780d0,0.840d0,0.880d0,0.950d0,1.000d0,1.000d0,1.000d0 &
212  ,0.030d0,0.030d0,0.024d0,0.022d0,0.032d0,0.062d0,0.200d0 &
213  ,0.680d0,0.830d0,0.870d0 ,0.900d0,0.950d0,1.000d0,1.000d0 &
214  ,1.000d0,0.025d0,0.025d0,0.025d0,0.036d0,0.043d0 ,0.130d0 &
215  ,0.270d0,0.740d0,0.860d0,0.890d0,0.920d0,1.000d0,1.000d0 &
216  ,1.000d0,1.000d0 ,0.027d0,0.027d0,0.027d0,0.040d0,0.052d0 &
217  ,0.200d0,0.400d0,0.780d0,0.880d0,0.900d0 ,0.940d0,1.000d0 &
218  ,1.000d0,1.000d0,1.000d0,0.030d0,0.030d0,0.030d0,0.047d0 &
219  ,0.064d0 ,0.250d0,0.500d0,0.800d0,0.900d0,0.910d0,0.950d0 &
220  ,1.000d0,1.000d0,1.000d0,1.000d0 ,0.040d0,0.040d0,0.033d0 &
221  ,0.037d0,0.068d0,0.240d0,0.550d0,0.800d0,0.900d0,0.910d0 &
222  ,0.950d0,1.000d0,1.000d0,1.000d0,1.000d0,0.035d0,0.035d0 &
223  ,0.035d0,0.055d0,0.079d0 ,0.290d0,0.580d0,0.800d0,0.900d0 &
224  ,0.910d0,0.950d0,1.000d0,1.000d0,1.000d0,1.000d0 ,0.037d0 &
225  ,0.037d0,0.037d0,0.062d0,0.082d0,0.290d0,0.590d0,0.780d0 &
226  ,0.900d0,0.910d0 ,0.950d0,1.000d0,1.000d0,1.000d0,1.000d0 &
227  ,0.037d0,0.037d0,0.037d0,0.060d0,0.080d0 ,0.290d0,0.580d0 &
228  ,0.770d0,0.890d0,0.910d0,0.950d0,1.000d0,1.000d0,1.000d0 &
229  ,1.000d0 ,0.037d0,0.037d0,0.037d0,0.041d0,0.075d0,0.250d0 &
230  ,0.540d0,0.760d0,0.880d0,0.920d0 ,0.950d0,1.000d0,1.000d0 &
231  ,1.000d0,1.000d0,0.037d0,0.037d0,0.037d0,0.052d0,0.067d0 &
232  ,0.250d0,0.510d0,0.770d0,0.880d0,0.930d0,0.970d0,1.000d0 &
233  ,1.000d0,1.000d0,1.000d0 ,0.037d0,0.037d0,0.037d0,0.047d0 &
234  ,0.057d0,0.250d0,0.490d0,0.770d0,0.890d0,0.950d0 ,1.000d0 &
235  ,1.000d0,1.000d0,1.000d0,1.000d0,0.036d0,0.036d0,0.036d0 &
236  ,0.042d0,0.048d0 ,0.230d0,0.470d0,0.780d0,0.920d0,1.000d0 &
237  ,1.020d0,1.020d0,1.020d0,1.020d0,1.020d0 ,0.040d0,0.040d0 &
238  ,0.035d0,0.033d0,0.040d0,0.112d0,0.450d0,0.790d0,1.010d0 &
239  ,1.030d0 ,1.040d0,1.040d0,1.040d0,1.040d0,1.040d0,0.033d0 &
240  ,0.033d0,0.033d0,0.033d0,0.033d0 ,0.119d0,0.470d0,0.950d0 &
241  ,1.300d0,1.700d0,2.300d0,2.300d0,2.300d0,2.300d0,2.300d0 &
242  ,0.027d0,0.027d0,0.027d0,0.027d0,0.027d0,0.125d0,0.520d0 &
243  ,1.400d0,2.300d0,3.000d0 ,4.000d0,4.000d0,4.000d0,4.000d0 &
244  ,4.000d0/
245 
246  r_small = min(r1 * 1d6, r2 * 1d6) ! um
247  r_big = max(r1 * 1d6, r2 * 1d6) ! um
248  rq = r_small / r_big
249 
250  ir = 1
251  do k = 1, 15
252  if (r_big .gt. r0(k)) then
253  ir = k + 1
254  end if
255  end do
256 
257  iq = 1
258  do kk = 1,21
259  if (rq .gt. rat(kk)) then
260  iq = kk + 1
261  end if
262  end do
263 
264  if (ir .lt. 16) then
265  if (ir .ge. 2) then
266  p = (r_big - r0(ir - 1)) / (r0(ir) - r0(ir - 1))
267  q = (rq - rat(iq - 1)) / (rat(iq) - rat(iq - 1))
268  ec = (1d0 - p) * (1d0 - q) * ecoll(ir - 1, iq - 1) &
269  + p * (1d0 - q) * ecoll(ir, iq - 1) &
270  + q * (1d0 - p) * ecoll(ir - 1, iq) &
271  + p * q * ecoll(ir, iq)
272  else
273  q = (rq - rat(iq - 1)) / (rat(iq) - rat(iq - 1))
274  ec = (1d0 - q) * ecoll(1, iq - 1) + q * ecoll(1, iq)
275  end if
276  else
277  q = (rq - rat(iq - 1)) / (rat(iq) - rat(iq - 1))
278  ek = (1d0 - q) * ecoll(15, iq - 1) + q * ecoll(15, iq)
279  ec = min(ek, 1d0)
280  end if
281 
282  if (ec .lt. 1d-20) stop 99
283 
284  end subroutine effic
285 
286 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
287 
288 end module pmc_coag_kernel_sedi
pmc_coag_kernel_sedi::effic
subroutine effic(r1, r2, ec)
Coagulation efficiency.
Definition: coag_kernel_sedi.F90:180
pmc_aero_particle::aero_particle_t
Single aerosol particle data structure.
Definition: aero_particle.F90:26
pmc_aero_data::aero_data_vol_to_mobility_rad
real(kind=dp) function aero_data_vol_to_mobility_rad(aero_data, v, temp, pressure)
Convert mass-equivalent volume (m^3) to mobility equivalent radius (m).
Definition: aero_data.F90:174
pmc_aero_particle
The aero_particle_t structure and associated subroutines.
Definition: aero_particle.F90:9
pmc_constants::dp
integer, parameter dp
Kind of a double precision real number.
Definition: constants.F90:12
pmc_env_state::env_state_t
Current environment state.
Definition: env_state.F90:29
pmc_coag_kernel_sedi::fall_g
subroutine fall_g(r, w_inf)
Finds the terminal velocity of a particle based on its size.
Definition: coag_kernel_sedi.F90:110
pmc_constants::const
type(const_t), save const
Fixed variable for accessing the constant's values.
Definition: constants.F90:73
pmc_env_state
The env_state_t structure and associated subroutines.
Definition: env_state.F90:9
pmc_aero_data::aero_data_vol2rad
real(kind=dp) elemental function aero_data_vol2rad(aero_data, v)
Convert mass-equivalent volume (m^3) to geometric radius (m).
Definition: aero_data.F90:92
pmc_coag_kernel_sedi
Gravitational sedimentation coagulation kernel.
Definition: coag_kernel_sedi.F90:17
pmc_aero_data::aero_data_t
Aerosol material properties and associated data.
Definition: aero_data.F90:49
pmc_constants
Physical constants.
Definition: constants.F90:9
pmc_coag_kernel_sedi::kernel_sedi_helper
subroutine kernel_sedi_helper(v1, v2, aero_data, temp, pressure, k)
Helper function that does the actual sedimentation kernel computation.
Definition: coag_kernel_sedi.F90:79
pmc_coag_kernel_sedi::kernel_sedi
subroutine kernel_sedi(aero_particle_1, aero_particle_2, aero_data, env_state, k)
Sedimentation coagulation kernel.
Definition: coag_kernel_sedi.F90:31
pmc_aero_data
The aero_data_t structure and associated subroutines.
Definition: aero_data.F90:9
pmc_coag_kernel_sedi::kernel_sedi_minmax
subroutine kernel_sedi_minmax(v1, v2, aero_data, env_state, k_min, k_max)
Minimum and maximum values of the sedimentation coagulation.
Definition: coag_kernel_sedi.F90:53
pmc_aero_particle::aero_particle_volume
elemental real(kind=dp) function aero_particle_volume(aero_particle)
Total volume of the particle (m^3).
Definition: aero_particle.F90:266
PartMC 2.6.1 documentation generated by doxygen 1.8.20