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