[825] | 1 | MODULE lpm_collision_kernels_mod |
---|
[790] | 2 | |
---|
| 3 | !------------------------------------------------------------------------------! |
---|
| 4 | ! Current revisions: |
---|
| 5 | ! ----------------- |
---|
[826] | 6 | ! |
---|
[790] | 7 | ! |
---|
| 8 | ! Former revisions: |
---|
| 9 | ! ----------------- |
---|
| 10 | ! $Id: lpm_collision_kernels.f90 826 2012-02-19 03:41:34Z raasch $ |
---|
| 11 | ! |
---|
[826] | 12 | ! 825 2012-02-19 03:03:44Z raasch |
---|
| 13 | ! routine renamed from wang_kernel to lpm_collision_kernels, |
---|
| 14 | ! turbulence_effects on collision replaced by wang_kernel |
---|
| 15 | ! |
---|
[800] | 16 | ! 799 2011-12-21 17:48:03Z franke |
---|
| 17 | ! speed optimizations and formatting |
---|
| 18 | ! Bugfix: iq=1 is not allowed (routine effic) |
---|
| 19 | ! Bugfix: replaced stop by ec=0.0 in case of very small ec (routine effic) |
---|
| 20 | ! |
---|
[791] | 21 | ! 790 2011-11-29 03:11:20Z raasch |
---|
| 22 | ! initial revision |
---|
[790] | 23 | ! |
---|
| 24 | ! Description: |
---|
| 25 | ! ------------ |
---|
| 26 | ! This routine calculates the effect of (SGS) turbulence on the collision |
---|
| 27 | ! efficiency of droplets. |
---|
| 28 | ! It is based on the original kernel developed by Wang (...) |
---|
| 29 | !------------------------------------------------------------------------------! |
---|
| 30 | |
---|
| 31 | USE arrays_3d |
---|
| 32 | USE cloud_parameters |
---|
| 33 | USE constants |
---|
| 34 | USE particle_attributes |
---|
| 35 | |
---|
| 36 | IMPLICIT NONE |
---|
| 37 | |
---|
| 38 | PRIVATE |
---|
| 39 | |
---|
[792] | 40 | PUBLIC colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi |
---|
[790] | 41 | |
---|
[792] | 42 | INTEGER, SAVE :: ip, jp, kp, pend, pstart, psum |
---|
[790] | 43 | |
---|
[792] | 44 | REAL, SAVE :: epsilon, eps2, urms, urms2 |
---|
| 45 | |
---|
| 46 | REAL, DIMENSION(:), ALLOCATABLE, SAVE :: winf |
---|
| 47 | REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: ec, ecf, gck |
---|
| 48 | |
---|
[790] | 49 | ! |
---|
| 50 | !-- Public interfaces |
---|
[792] | 51 | INTERFACE colker |
---|
| 52 | MODULE PROCEDURE colker |
---|
| 53 | END INTERFACE colker |
---|
[790] | 54 | |
---|
| 55 | INTERFACE effic |
---|
| 56 | MODULE PROCEDURE effic |
---|
| 57 | END INTERFACE effic |
---|
| 58 | |
---|
| 59 | INTERFACE fallg |
---|
| 60 | MODULE PROCEDURE fallg |
---|
| 61 | END INTERFACE fallg |
---|
| 62 | |
---|
[792] | 63 | INTERFACE phi |
---|
| 64 | MODULE PROCEDURE phi |
---|
| 65 | END INTERFACE phi |
---|
[790] | 66 | |
---|
[792] | 67 | INTERFACE turbsd |
---|
| 68 | MODULE PROCEDURE turbsd |
---|
| 69 | END INTERFACE turbsd |
---|
| 70 | |
---|
| 71 | INTERFACE turb_enhance_eff |
---|
| 72 | MODULE PROCEDURE turb_enhance_eff |
---|
| 73 | END INTERFACE turb_enhance_eff |
---|
| 74 | |
---|
| 75 | INTERFACE zhi |
---|
| 76 | MODULE PROCEDURE zhi |
---|
| 77 | END INTERFACE zhi |
---|
| 78 | |
---|
[790] | 79 | CONTAINS |
---|
| 80 | |
---|
| 81 | !------------------------------------------------------------------------------! |
---|
| 82 | ! SUBROUTINE for calculation of collision kernel |
---|
| 83 | !------------------------------------------------------------------------------! |
---|
[792] | 84 | SUBROUTINE colker( i1, j1, k1, kernel ) |
---|
[790] | 85 | |
---|
| 86 | USE arrays_3d |
---|
| 87 | USE cloud_parameters |
---|
| 88 | USE constants |
---|
[792] | 89 | USE cpulog |
---|
[790] | 90 | USE indices |
---|
[792] | 91 | USE interfaces |
---|
[790] | 92 | USE particle_attributes |
---|
| 93 | |
---|
| 94 | IMPLICIT NONE |
---|
| 95 | |
---|
[792] | 96 | INTEGER :: i, i1, j, j1, k1 |
---|
[790] | 97 | |
---|
| 98 | REAL, DIMENSION(prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1,prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1) :: kernel |
---|
| 99 | |
---|
[792] | 100 | ! CALL cpu_log( log_point_s(46), 'colker', 'start' ) |
---|
| 101 | |
---|
[790] | 102 | ip = i1 |
---|
| 103 | jp = j1 |
---|
| 104 | kp = k1 |
---|
| 105 | |
---|
| 106 | pstart = prt_start_index(kp,jp,ip) |
---|
| 107 | pend = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1 |
---|
| 108 | psum = prt_count(kp,jp,ip) |
---|
| 109 | |
---|
[792] | 110 | ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) ) |
---|
[790] | 111 | |
---|
[825] | 112 | IF ( wang_kernel ) THEN |
---|
[790] | 113 | |
---|
[792] | 114 | ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) ) |
---|
[790] | 115 | |
---|
[792] | 116 | epsilon = diss(kp,jp,ip) * 1.0E5 !dissipation rate in cm**2/s**-3 |
---|
| 117 | urms = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 ) |
---|
[790] | 118 | |
---|
[792] | 119 | IF ( epsilon <= 0.001 ) THEN |
---|
[790] | 120 | |
---|
| 121 | CALL fallg |
---|
| 122 | CALL effic |
---|
| 123 | |
---|
[792] | 124 | DO j = pstart, pend |
---|
| 125 | DO i = pstart, pend |
---|
| 126 | kernel(i,j) = pi * ( particles(j)%radius * 100.0 + & |
---|
| 127 | particles(i)%radius * 100.0 )**2 & |
---|
| 128 | * ec(i,j) * ABS( winf(j) - winf(i) ) |
---|
[790] | 129 | ENDDO |
---|
| 130 | ENDDO |
---|
| 131 | |
---|
| 132 | ELSE |
---|
| 133 | |
---|
[792] | 134 | CALL turbsd |
---|
[790] | 135 | CALL turb_enhance_eff |
---|
| 136 | CALL effic |
---|
| 137 | |
---|
[792] | 138 | DO j = pstart, pend, 1 |
---|
| 139 | DO i = pstart, pend, 1 |
---|
[790] | 140 | kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j) |
---|
| 141 | ENDDO |
---|
| 142 | ENDDO |
---|
[792] | 143 | |
---|
[790] | 144 | ENDIF |
---|
| 145 | |
---|
| 146 | DEALLOCATE(gck, ecf) |
---|
| 147 | |
---|
| 148 | ELSE |
---|
| 149 | |
---|
[792] | 150 | ! CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' ) |
---|
[790] | 151 | CALL fallg |
---|
[792] | 152 | ! CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' ) |
---|
| 153 | ! CALL cpu_log( log_point_s(51), 'colker_effic', 'start' ) |
---|
[790] | 154 | CALL effic |
---|
[792] | 155 | ! CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' ) |
---|
[790] | 156 | |
---|
[792] | 157 | DO j = pstart, pend |
---|
| 158 | DO i = pstart, pend |
---|
| 159 | kernel(i,j) = pi * ( particles(j)%radius * 100.0 + & |
---|
| 160 | particles(i)%radius * 100.0 )**2 & |
---|
| 161 | * ec(i,j) * ABS( winf(j) - winf(i) ) |
---|
[790] | 162 | ENDDO |
---|
| 163 | ENDDO |
---|
| 164 | |
---|
| 165 | ENDIF |
---|
| 166 | |
---|
[792] | 167 | DEALLOCATE( ec, winf ) |
---|
[790] | 168 | |
---|
[792] | 169 | ! CALL cpu_log( log_point_s(46), 'colker', 'stop' ) |
---|
[790] | 170 | |
---|
| 171 | END SUBROUTINE colker |
---|
| 172 | |
---|
| 173 | !------------------------------------------------------------------------------! |
---|
| 174 | ! SUBROUTINE for calculation of w, g and gck |
---|
| 175 | !------------------------------------------------------------------------------! |
---|
[792] | 176 | SUBROUTINE turbsd |
---|
| 177 | ! from Aayala 2008b, page 37ff, necessary input parameter water density, radii |
---|
| 178 | ! of droplets, air density, air viscosity, turbulent dissipation rate, |
---|
| 179 | ! taylor microscale reynolds number, gravitational acceleration |
---|
[799] | 180 | |
---|
[790] | 181 | USE constants |
---|
| 182 | USE cloud_parameters |
---|
| 183 | USE particle_attributes |
---|
| 184 | USE arrays_3d |
---|
| 185 | |
---|
| 186 | IMPLICIT NONE |
---|
| 187 | |
---|
[799] | 188 | REAL :: Relamda, & |
---|
[790] | 189 | Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be, & |
---|
| 190 | bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq, & |
---|
| 191 | vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, & |
---|
| 192 | wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp |
---|
| 193 | |
---|
[799] | 194 | REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens |
---|
[790] | 195 | |
---|
[799] | 196 | REAL, DIMENSION(pstart:pend) :: St, tau |
---|
| 197 | |
---|
| 198 | LOGICAL, SAVE :: first = .TRUE. |
---|
| 199 | |
---|
[790] | 200 | INTEGER :: i, j |
---|
| 201 | |
---|
[799] | 202 | ! |
---|
| 203 | !-- Initial assignment of constants |
---|
| 204 | IF ( first ) THEN |
---|
[790] | 205 | |
---|
[799] | 206 | first = .FALSE. |
---|
| 207 | airvisc = 0.1818 !dynamic viscosity in mg/cm*s |
---|
| 208 | airdens = 1.2250 !air density in mg/cm**3 |
---|
| 209 | waterdens = 1000.0 !water density in mg/cm**3 |
---|
| 210 | gravity = 980.6650 !in cm/s**2 |
---|
| 211 | anu = airvisc/airdens ! kinetic viscosity in cm**2/s |
---|
[790] | 212 | |
---|
[799] | 213 | ENDIF |
---|
[790] | 214 | |
---|
[799] | 215 | Relamda = urms**2*sqrt(15.0/epsilon/anu) |
---|
[790] | 216 | |
---|
[799] | 217 | Tl = urms**2/epsilon !in s |
---|
| 218 | Lf = 0.5 * (urms**3)/epsilon !in cm |
---|
| 219 | |
---|
[790] | 220 | tauk = (anu/epsilon)**0.5 !in s |
---|
[799] | 221 | eta = (anu**3/epsilon)**0.25 !in cm |
---|
[790] | 222 | vk = eta/tauk !in cm/s |
---|
| 223 | |
---|
| 224 | ao = (11.+7.*Relamda)/(205.+Relamda) |
---|
| 225 | |
---|
| 226 | lambda = urms * sqrt(15.*anu/epsilon) !in cm |
---|
| 227 | |
---|
| 228 | tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s |
---|
| 229 | |
---|
| 230 | CALL fallg !gives winf in cm/s |
---|
| 231 | |
---|
| 232 | DO i = pstart, pend |
---|
[799] | 233 | tau(i) = winf(i)/gravity !in s |
---|
[790] | 234 | St(i) = tau(i)/tauk |
---|
| 235 | ENDDO |
---|
| 236 | |
---|
| 237 | !***** TO CALCULATE wr ******************************** |
---|
| 238 | !from Aayala 2008b, page 38f |
---|
| 239 | |
---|
| 240 | z = tt/Tl |
---|
| 241 | be = sqrt(2.0)*lambda/Lf |
---|
| 242 | |
---|
[799] | 243 | bbb = sqrt(1.0-2.0*be**2) |
---|
[790] | 244 | d1 = (1.+bbb)/2.0/bbb |
---|
| 245 | e1 = Lf*(1.0+bbb)/2.0 !in cm |
---|
| 246 | d2 = (1.0-bbb)/2.0/bbb |
---|
| 247 | e2 = Lf*(1.0-bbb)/2.0 !in cm |
---|
| 248 | |
---|
[799] | 249 | ccc = sqrt(1.0-2.0*z**2) |
---|
[790] | 250 | b1 = (1.+ccc)/2./ccc |
---|
| 251 | c1 = Tl*(1.+ccc)/2. !in s |
---|
| 252 | b2 = (1.-ccc)/2./ccc |
---|
| 253 | c2 = Tl*(1.-ccc)/2. !in s |
---|
| 254 | |
---|
| 255 | DO i = pstart, pend |
---|
[799] | 256 | v1 = winf(i) !in cm/s |
---|
[790] | 257 | t1 = tau(i) !in s |
---|
| 258 | |
---|
| 259 | DO j = pstart,i |
---|
| 260 | rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm |
---|
[799] | 261 | v2 = winf(j) !in cm/s |
---|
[790] | 262 | t2 = tau(j) !in s |
---|
| 263 | |
---|
| 264 | v1xysq = b1*d1*PHI(c1,e1,v1,t1) & |
---|
| 265 | - b1*d2*PHI(c1,e2,v1,t1) & |
---|
| 266 | - b2*d1*PHI(c2,e1,v1,t1) & |
---|
| 267 | + b2*d2*PHI(c2,e2,v1,t1) |
---|
[799] | 268 | v1xysq = v1xysq * urms**2/t1 !in cm**2/s**2 |
---|
[790] | 269 | |
---|
| 270 | vrms1xy= sqrt(v1xysq) !in cm/s |
---|
| 271 | |
---|
| 272 | v2xysq = b1*d1*PHI(c1,e1,v2,t2) & |
---|
| 273 | - b1*d2*PHI(c1,e2,v2,t2) & |
---|
| 274 | - b2*d1*PHI(c2,e1,v2,t2) & |
---|
| 275 | + b2*d2*PHI(c2,e2,v2,t2) |
---|
[799] | 276 | v2xysq = v2xysq * urms**2/t2 !in cm**2/s**2 |
---|
[790] | 277 | |
---|
| 278 | vrms2xy= sqrt(v2xysq) !in cm/s |
---|
| 279 | |
---|
[799] | 280 | IF(winf(i).ge.winf(j)) THEN |
---|
| 281 | v1 = winf(i) |
---|
[790] | 282 | t1 = tau(i) |
---|
[799] | 283 | v2 = winf(j) |
---|
[790] | 284 | t2 = tau(j) |
---|
| 285 | ELSE |
---|
[799] | 286 | v1 = winf(j) |
---|
[790] | 287 | t1 = tau(j) |
---|
[799] | 288 | v2 = winf(i) |
---|
[790] | 289 | t2 = tau(i) |
---|
| 290 | ENDIF |
---|
| 291 | |
---|
| 292 | v1v2xy = b1*d1*ZHI(c1,e1,v1,t1,v2,t2) & |
---|
| 293 | - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) & |
---|
| 294 | - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) & |
---|
| 295 | + b2*d2*ZHI(c2,e2,v1,t1,v2,t2) |
---|
| 296 | fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2) |
---|
[799] | 297 | v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j) !in cm**2/s**2 |
---|
[790] | 298 | |
---|
[799] | 299 | wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy !in cm**2/s**2 |
---|
[790] | 300 | |
---|
| 301 | IF (wrtur2xy.le.0.0) wrtur2xy=0.0 |
---|
| 302 | |
---|
[799] | 303 | wrgrav2=pi/8.*(winf(j)-winf(i))**2 |
---|
[790] | 304 | |
---|
| 305 | wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2)) !in cm/s |
---|
| 306 | |
---|
| 307 | |
---|
| 308 | !***** TO CALCULATE gr ******************************** |
---|
| 309 | |
---|
| 310 | IF(St(j).gt.St(i)) THEN |
---|
| 311 | SSt = St(j) |
---|
| 312 | ELSE |
---|
| 313 | SSt = St(i) |
---|
| 314 | ENDIF |
---|
| 315 | |
---|
[799] | 316 | XX = -0.1988*SSt**4 + 1.5275*SSt**3 & |
---|
| 317 | -4.2942*SSt**2 + 5.3406*SSt |
---|
[790] | 318 | |
---|
| 319 | IF(XX.le.0.0) XX = 0.0 |
---|
| 320 | |
---|
| 321 | YY = 0.1886*exp(20.306/Relamda) |
---|
| 322 | |
---|
| 323 | c1_gr = XX/(gravity/(vk/tauk))**YY |
---|
| 324 | |
---|
[799] | 325 | ao_gr = ao + (pi/8.)*(gravity/(vk/tauk))**2 |
---|
[790] | 326 | fao_gr = 20.115 * (ao_gr/Relamda)**0.5 |
---|
| 327 | rc = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta !in cm |
---|
| 328 | |
---|
[799] | 329 | grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.) |
---|
[790] | 330 | IF (grFIN.lt.1.0) grFIN = 1.0 |
---|
| 331 | |
---|
[799] | 332 | gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN ! in cm**3/s |
---|
[790] | 333 | gck(j,i) = gck(i,j) |
---|
| 334 | |
---|
| 335 | ENDDO |
---|
| 336 | ENDDO |
---|
| 337 | |
---|
| 338 | END SUBROUTINE TurbSD |
---|
| 339 | |
---|
| 340 | !------------------------------------------------------------------------------! |
---|
| 341 | ! PHI as a function |
---|
| 342 | !------------------------------------------------------------------------------! |
---|
| 343 | REAL FUNCTION PHI(a,b,vsett,tau0) |
---|
| 344 | |
---|
| 345 | IMPLICIT NONE |
---|
| 346 | |
---|
| 347 | REAL :: a, aa1, b, vsett, tau0 |
---|
| 348 | |
---|
| 349 | aa1 = 1./tau0 + 1./a + vsett/b |
---|
| 350 | |
---|
[799] | 351 | PHI = 1./aa1 - vsett/2.0/b/aa1**2 !in s |
---|
[792] | 352 | |
---|
[790] | 353 | END FUNCTION PHI |
---|
| 354 | |
---|
| 355 | !------------------------------------------------------------------------------! |
---|
| 356 | ! ZETA as a function |
---|
| 357 | !------------------------------------------------------------------------------! |
---|
| 358 | REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2) |
---|
| 359 | |
---|
| 360 | IMPLICIT NONE |
---|
| 361 | |
---|
| 362 | REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2 |
---|
| 363 | |
---|
| 364 | aa1 = vsett2/b - 1./tau2 - 1./a |
---|
| 365 | aa2 = vsett1/b + 1./tau1 + 1./a |
---|
| 366 | aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2 |
---|
[799] | 367 | aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2 |
---|
[790] | 368 | aa5 = vsett2/b + 1./tau2 + 1./a |
---|
| 369 | aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2 |
---|
[799] | 370 | ZHI = (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2 & |
---|
| 371 | + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6 & |
---|
| 372 | + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2) & |
---|
[790] | 373 | * 1./2./b/aa3 ! in s**2 |
---|
[799] | 374 | |
---|
[790] | 375 | END FUNCTION ZHI |
---|
| 376 | |
---|
| 377 | !------------------------------------------------------------------------------! |
---|
| 378 | ! SUBROUTINE for calculation of terminal velocity winf |
---|
| 379 | !------------------------------------------------------------------------------! |
---|
| 380 | SUBROUTINE fallg |
---|
| 381 | |
---|
| 382 | USE constants |
---|
| 383 | USE cloud_parameters |
---|
| 384 | USE particle_attributes |
---|
| 385 | USE arrays_3d |
---|
| 386 | |
---|
| 387 | IMPLICIT NONE |
---|
| 388 | |
---|
| 389 | INTEGER :: i, j |
---|
| 390 | |
---|
[799] | 391 | LOGICAL, SAVE :: first = .TRUE. |
---|
[790] | 392 | |
---|
[799] | 393 | REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py |
---|
[790] | 394 | |
---|
[799] | 395 | REAL :: bond, x, xrey, y |
---|
| 396 | |
---|
| 397 | REAL, DIMENSION(1:7), SAVE :: b |
---|
| 398 | REAL, DIMENSION(1:6), SAVE :: c |
---|
| 399 | |
---|
| 400 | ! |
---|
| 401 | !-- Initial assignment of constants |
---|
| 402 | IF ( first ) THEN |
---|
| 403 | |
---|
| 404 | first = .FALSE. |
---|
| 405 | b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, & |
---|
[790] | 406 | -0.578878e-3,0.855176e-4,-0.327815e-5/) |
---|
[799] | 407 | c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0, & |
---|
[790] | 408 | -0.542819e-1, 0.238449e-2/) |
---|
| 409 | |
---|
[799] | 410 | eta = 1.818e-4 !in poise = g/(cm s) |
---|
| 411 | xlamb = 6.62e-6 !in cm |
---|
[790] | 412 | |
---|
[799] | 413 | rhow = 1.0 !in g/cm**3 |
---|
| 414 | rhoa = 1.225e-3 !in g/cm**3 |
---|
[790] | 415 | |
---|
[799] | 416 | grav = 980.665 !in cm/s**2 |
---|
| 417 | cunh = 1.257 * xlamb !in cm |
---|
| 418 | t0 = 273.15 !in K |
---|
| 419 | sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2 |
---|
| 420 | stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s) |
---|
| 421 | stb = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3 |
---|
| 422 | phy = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa)) |
---|
| 423 | py = phy**(1.0/6.0) |
---|
[790] | 424 | |
---|
[799] | 425 | ENDIF |
---|
| 426 | |
---|
[790] | 427 | !particle radius has to be in cm |
---|
| 428 | DO j = pstart, pend |
---|
| 429 | |
---|
| 430 | IF (particles(j)%radius*100.0 .le. 1.e-3) THEN |
---|
| 431 | |
---|
[799] | 432 | winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s |
---|
[790] | 433 | |
---|
| 434 | ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN |
---|
| 435 | |
---|
[799] | 436 | x = log(stb*(particles(j)%radius*100.0)**3) |
---|
[790] | 437 | y = 0.0 |
---|
| 438 | |
---|
| 439 | DO i = 1, 7 |
---|
| 440 | y = y + b(i) * (x**(i-1)) |
---|
| 441 | ENDDO |
---|
| 442 | |
---|
| 443 | xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y) |
---|
| 444 | winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0) !in cm/s |
---|
| 445 | |
---|
| 446 | ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN |
---|
| 447 | |
---|
| 448 | IF (particles(j)%radius*100.0.gt.0.35) THEN |
---|
| 449 | bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma |
---|
[799] | 450 | ELSE |
---|
| 451 | bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma |
---|
[790] | 452 | ENDIF |
---|
| 453 | |
---|
| 454 | x = log(16.0*bond*py/3.0) |
---|
| 455 | y = 0.0 |
---|
| 456 | |
---|
| 457 | DO i = 1, 6 |
---|
| 458 | y = y + c(i) * (x**(i-1)) |
---|
| 459 | ENDDO |
---|
| 460 | |
---|
| 461 | xrey = py*exp(y) |
---|
| 462 | |
---|
| 463 | IF (particles(j)%radius*100.0 .gt.0.35) THEN |
---|
| 464 | winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s |
---|
[799] | 465 | ELSE |
---|
| 466 | winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0) !in cm/s |
---|
[790] | 467 | ENDIF |
---|
| 468 | |
---|
| 469 | ENDIF |
---|
| 470 | ENDDO |
---|
[799] | 471 | RETURN |
---|
[790] | 472 | END SUBROUTINE fallg |
---|
| 473 | |
---|
| 474 | !------------------------------------------------------------------------------! |
---|
| 475 | ! SUBROUTINE for calculation of collision efficencies |
---|
| 476 | !------------------------------------------------------------------------------! |
---|
| 477 | |
---|
| 478 | SUBROUTINE effic |
---|
| 479 | |
---|
[792] | 480 | USE arrays_3d |
---|
| 481 | USE constants |
---|
| 482 | USE cloud_parameters |
---|
| 483 | USE particle_attributes |
---|
[790] | 484 | |
---|
[799] | 485 | !collision efficiencies of hall kernel |
---|
[790] | 486 | IMPLICIT NONE |
---|
| 487 | |
---|
| 488 | INTEGER :: i, ir, iq, j, k, kk |
---|
| 489 | |
---|
[792] | 490 | INTEGER, DIMENSION(:), ALLOCATABLE :: ira |
---|
[790] | 491 | |
---|
[792] | 492 | LOGICAL, SAVE :: first = .TRUE. |
---|
[790] | 493 | |
---|
[792] | 494 | REAL :: ek, particle_radius, pp, qq, rq |
---|
[790] | 495 | |
---|
[792] | 496 | REAL, DIMENSION(1:21), SAVE :: rat |
---|
| 497 | REAL, DIMENSION(1:15), SAVE :: r0 |
---|
| 498 | REAL, DIMENSION(1:15,1:21), SAVE :: ecoll |
---|
[790] | 499 | |
---|
[792] | 500 | ! |
---|
| 501 | !-- Initial assignment of constants |
---|
| 502 | IF ( first ) THEN |
---|
[790] | 503 | |
---|
[792] | 504 | first = .FALSE. |
---|
| 505 | r0 = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., & |
---|
| 506 | 300./) |
---|
| 507 | rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6, & |
---|
| 508 | 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/) |
---|
[790] | 509 | |
---|
[792] | 510 | ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, & |
---|
| 511 | 0.001,0.001,0.001,0.001,0.001,0.001/) |
---|
| 512 | ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, & |
---|
| 513 | 0.050,0.200,0.500,0.770,0.870,0.970/) |
---|
| 514 | ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, & |
---|
| 515 | 0.430,0.580,0.790,0.930,0.960,1.000/) |
---|
| 516 | ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, & |
---|
| 517 | 0.640,0.750,0.910,0.970,0.980,1.000/) |
---|
| 518 | ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, & |
---|
| 519 | 0.770,0.840,0.950,0.970,1.000,1.000/) |
---|
| 520 | ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, & |
---|
| 521 | 0.840,0.880,0.950,1.000,1.000,1.000/) |
---|
| 522 | ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, & |
---|
| 523 | 0.870,0.900,0.950,1.000,1.000,1.000/) |
---|
| 524 | ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, & |
---|
| 525 | 0.890,0.920,1.000,1.000,1.000,1.000/) |
---|
| 526 | ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, & |
---|
| 527 | 0.900,0.940,1.000,1.000,1.000,1.000/) |
---|
| 528 | ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, & |
---|
| 529 | 0.910,0.950,1.000,1.000,1.000,1.000/) |
---|
| 530 | ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, & |
---|
| 531 | 0.910,0.950,1.000,1.000,1.000,1.000/) |
---|
| 532 | ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, & |
---|
| 533 | 0.910,0.950,1.000,1.000,1.000,1.000/) |
---|
| 534 | ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, & |
---|
| 535 | 0.910,0.950,1.000,1.000,1.000,1.000/) |
---|
| 536 | ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, & |
---|
| 537 | 0.910,0.950,1.000,1.000,1.000,1.000/) |
---|
| 538 | ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, & |
---|
| 539 | 0.920,0.950,1.000,1.000,1.000,1.000/) |
---|
| 540 | ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, & |
---|
| 541 | 0.930,0.970,1.000,1.000,1.000,1.000/) |
---|
| 542 | ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, & |
---|
| 543 | 0.950,1.000,1.000,1.000,1.000,1.000/) |
---|
| 544 | ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, & |
---|
| 545 | 1.000,1.020,1.020,1.020,1.020,1.020/) |
---|
| 546 | ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, & |
---|
| 547 | 1.030,1.040,1.040,1.040,1.040,1.040/) |
---|
| 548 | ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, & |
---|
| 549 | 1.700,2.300,2.300,2.300,2.300,2.300/) |
---|
| 550 | ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, & |
---|
| 551 | 3.000,4.000,4.000,4.000,4.000,4.000/) |
---|
| 552 | ENDIF |
---|
[790] | 553 | |
---|
[792] | 554 | ! |
---|
| 555 | !-- Calculate the radius class index of particles with respect to array r |
---|
| 556 | ALLOCATE( ira(pstart:pend) ) |
---|
| 557 | DO j = pstart, pend |
---|
| 558 | particle_radius = particles(j)%radius * 1.0E6 |
---|
| 559 | DO k = 1, 15 |
---|
| 560 | IF ( particle_radius < r0(k) ) THEN |
---|
| 561 | ira(j) = k |
---|
| 562 | EXIT |
---|
| 563 | ENDIF |
---|
| 564 | ENDDO |
---|
| 565 | IF ( particle_radius >= r0(15) ) ira(j) = 16 |
---|
| 566 | ENDDO |
---|
[790] | 567 | |
---|
[792] | 568 | ! |
---|
| 569 | !-- Two-dimensional linear interpolation of the collision efficiency. |
---|
| 570 | !-- Radius has to be in µm |
---|
| 571 | DO j = pstart, pend |
---|
| 572 | DO i = pstart, j |
---|
| 573 | |
---|
| 574 | ir = ira(j) |
---|
| 575 | |
---|
| 576 | rq = particles(i)%radius / particles(j)%radius |
---|
| 577 | |
---|
| 578 | ! DO kk = 2, 21 |
---|
[799] | 579 | ! IF ( rq <= rat(kk) ) THEN |
---|
| 580 | ! iq = kk |
---|
| 581 | ! EXIT |
---|
| 582 | ! ENDIF |
---|
[792] | 583 | ! ENDDO |
---|
| 584 | |
---|
| 585 | iq = INT( rq * 20 ) + 1 |
---|
[799] | 586 | iq = MAX(iq , 2) |
---|
[792] | 587 | |
---|
| 588 | IF ( ir < 16 ) THEN |
---|
| 589 | |
---|
| 590 | IF ( ir >= 2 ) THEN |
---|
| 591 | pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / & |
---|
| 592 | ( r0(ir) - r0(ir-1) ) |
---|
| 593 | qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) |
---|
| 594 | ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1) & |
---|
| 595 | + pp * ( 1.0-qq ) * ecoll(ir,iq-1) & |
---|
| 596 | + qq * ( 1.0-pp ) * ecoll(ir-1,iq) & |
---|
| 597 | + pp * qq * ecoll(ir,iq) |
---|
[790] | 598 | ELSE |
---|
| 599 | qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) |
---|
| 600 | ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq) |
---|
| 601 | ENDIF |
---|
[792] | 602 | |
---|
[790] | 603 | ELSE |
---|
[792] | 604 | qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) |
---|
| 605 | ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq) |
---|
| 606 | ec(j,i) = MIN( ek, 1.0 ) |
---|
[790] | 607 | ENDIF |
---|
[792] | 608 | |
---|
[790] | 609 | ec(i,j) = ec(j,i) |
---|
[799] | 610 | IF ( ec(i,j) < 1.0E-20 ) ec(i,j) = 0.0 |
---|
[792] | 611 | |
---|
[790] | 612 | ENDDO |
---|
| 613 | ENDDO |
---|
[792] | 614 | |
---|
| 615 | DEALLOCATE( ira ) |
---|
| 616 | |
---|
[790] | 617 | END SUBROUTINE effic |
---|
| 618 | |
---|
| 619 | !------------------------------------------------------------------------------! |
---|
| 620 | ! SUBROUTINE for calculation of enhancement factor collision efficencies |
---|
| 621 | !------------------------------------------------------------------------------! |
---|
| 622 | SUBROUTINE turb_enhance_eff |
---|
| 623 | |
---|
| 624 | USE constants |
---|
| 625 | USE cloud_parameters |
---|
| 626 | USE particle_attributes |
---|
| 627 | USE arrays_3d |
---|
| 628 | |
---|
| 629 | IMPLICIT NONE |
---|
| 630 | |
---|
| 631 | INTEGER :: i, ik, ir, iq, j, k, kk |
---|
| 632 | |
---|
[799] | 633 | INTEGER, DIMENSION(:), ALLOCATABLE :: ira |
---|
[790] | 634 | |
---|
[799] | 635 | REAL :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3 |
---|
[790] | 636 | |
---|
[799] | 637 | LOGICAL, SAVE :: first = .TRUE. |
---|
| 638 | |
---|
| 639 | REAL, DIMENSION(1:11), SAVE :: rat |
---|
| 640 | REAL, DIMENSION(1:7), SAVE :: r0 |
---|
| 641 | REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400 |
---|
| 642 | |
---|
| 643 | ! |
---|
| 644 | !-- Initial assignment of constants |
---|
| 645 | IF ( first ) THEN |
---|
| 646 | |
---|
| 647 | first = .FALSE. |
---|
| 648 | |
---|
| 649 | r0 = (/10., 20., 30.,40., 50., 60.,100./) |
---|
| 650 | rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/) |
---|
| 651 | |
---|
[790] | 652 | ! 100 cm^2/s^3 |
---|
[799] | 653 | ecoll_100(:,1) = (/1.74, 1.74, 1.773, 1.49, 1.207, 1.207, 1.0 /) |
---|
| 654 | ecoll_100(:,2) = (/1.46, 1.46, 1.421, 1.245, 1.069, 1.069, 1.0 /) |
---|
| 655 | ecoll_100(:,3) = (/1.32, 1.32, 1.245, 1.123, 1.000, 1.000, 1.0 /) |
---|
| 656 | ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025, 1.025, 1.0 /) |
---|
| 657 | ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056, 1.056, 1.0 /) |
---|
| 658 | ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028, 1.028, 1.0 /) |
---|
| 659 | ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046, 1.046, 1.0 /) |
---|
| 660 | ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029, 1.029, 1.0 /) |
---|
| 661 | ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021, 1.021, 1.0 /) |
---|
| 662 | ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088, 1.088, 1.0 /) |
---|
| 663 | ecoll_100(:,11)= (/20.3, 20.3, 14.6 , 8.61, 2.60, 2.60 , 1.0 /) |
---|
[790] | 664 | |
---|
| 665 | ! 400 cm^2/s^3 |
---|
[799] | 666 | ecoll_400(:,1) = (/4.976, 4.976, 3.593, 2.519, 1.445, 1.445, 1.0 /) |
---|
| 667 | ecoll_400(:,2) = (/2.984, 2.984, 2.181, 1.691, 1.201, 1.201, 1.0 /) |
---|
| 668 | ecoll_400(:,3) = (/1.988, 1.988, 1.475, 1.313, 1.150, 1.150, 1.0 /) |
---|
| 669 | ecoll_400(:,4) = (/1.490, 1.490, 1.187, 1.156, 1.126, 1.126, 1.0 /) |
---|
| 670 | ecoll_400(:,5) = (/1.249, 1.249, 1.088, 1.090, 1.092, 1.092, 1.0 /) |
---|
| 671 | ecoll_400(:,6) = (/1.139, 1.139, 1.130, 1.091, 1.051, 1.051, 1.0 /) |
---|
| 672 | ecoll_400(:,7) = (/1.220, 1.220, 1.190, 1.138, 1.086, 1.086, 1.0 /) |
---|
| 673 | ecoll_400(:,8) = (/1.325, 1.325, 1.267, 1.165, 1.063, 1.063, 1.0 /) |
---|
| 674 | ecoll_400(:,9) = (/1.716, 1.716, 1.345, 1.223, 1.100, 1.100, 1.0 /) |
---|
| 675 | ecoll_400(:,10)= (/3.788, 3.788, 1.501, 1.311, 1.120, 1.120, 1.0 /) |
---|
| 676 | ecoll_400(:,11)= (/36.52, 36.52, 19.16, 22.80, 26.0, 26.0, 1.0 /) |
---|
[790] | 677 | |
---|
[799] | 678 | ENDIF |
---|
| 679 | |
---|
| 680 | ! |
---|
| 681 | !-- Calculate the radius class index of particles with respect to array r |
---|
| 682 | ALLOCATE( ira(pstart:pend) ) |
---|
| 683 | |
---|
| 684 | DO j = pstart, pend |
---|
| 685 | particle_radius = particles(j)%radius * 1.0E6 |
---|
| 686 | DO k = 1, 7 |
---|
| 687 | IF ( particle_radius < r0(k) ) THEN |
---|
| 688 | ira(j) = k |
---|
| 689 | EXIT |
---|
| 690 | ENDIF |
---|
| 691 | ENDDO |
---|
| 692 | IF ( particle_radius >= r0(7) ) ira(j) = 8 |
---|
| 693 | ENDDO |
---|
| 694 | |
---|
[790] | 695 | ! two-dimensional linear interpolation of the collision efficiency |
---|
| 696 | DO j = pstart, pend |
---|
| 697 | DO i = pstart, j |
---|
| 698 | |
---|
[799] | 699 | ir = ira(j) |
---|
[790] | 700 | |
---|
[799] | 701 | rq = particles(i)%radius/particles(j)%radius |
---|
[790] | 702 | |
---|
| 703 | DO kk = 2, 11 |
---|
[799] | 704 | IF ( rq <= rat(kk) ) THEN |
---|
| 705 | iq = kk |
---|
| 706 | EXIT |
---|
| 707 | ENDIF |
---|
[790] | 708 | ENDDO |
---|
| 709 | |
---|
| 710 | ! 0 cm2/s3 |
---|
| 711 | y1 = 1.0 |
---|
| 712 | ! 100 cm2/s3, 400 cm2/s3 |
---|
| 713 | IF (ir.lt.8) THEN |
---|
| 714 | IF (ir.ge.2) THEN |
---|
| 715 | pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1)) |
---|
| 716 | qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) |
---|
| 717 | y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+ & |
---|
| 718 | pp*(1.-qq)*ecoll_100(ir,iq-1)+ & |
---|
| 719 | qq*(1.-pp)*ecoll_100(ir-1,iq)+ & |
---|
| 720 | pp*qq*ecoll_100(ir,iq) |
---|
| 721 | y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+ & |
---|
| 722 | pp*(1.-qq)*ecoll_400(ir,iq-1)+ & |
---|
| 723 | qq*(1.-pp)*ecoll_400(ir-1,iq)+ & |
---|
| 724 | pp*qq*ecoll_400(ir,iq) |
---|
| 725 | ELSE |
---|
| 726 | qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) |
---|
| 727 | y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq) |
---|
| 728 | y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq) |
---|
| 729 | ENDIF |
---|
| 730 | ELSE |
---|
| 731 | qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) |
---|
| 732 | y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq) |
---|
| 733 | y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq) |
---|
| 734 | ENDIF |
---|
| 735 | ! linear interpolation |
---|
| 736 | ! dissipation rate in cm2/s3 |
---|
| 737 | x1 = 0.0 |
---|
| 738 | x2 = 100.0 |
---|
| 739 | x3 = 400.0 |
---|
| 740 | |
---|
| 741 | IF (epsilon.le.100.) THEN |
---|
| 742 | ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 & |
---|
| 743 | + (epsilon-0.)/(100.-0.) * y2 |
---|
| 744 | ELSE IF(epsilon.le.600.)THEN |
---|
| 745 | ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 & |
---|
| 746 | + (epsilon-100.)/(400.-100.) * y3 |
---|
| 747 | |
---|
| 748 | ELSE |
---|
| 749 | ecf(j,i) = (600.-400.)/(100.-400.) * y2 & |
---|
| 750 | + (600.-100.)/(400.-100.) * y3 |
---|
| 751 | ENDIF |
---|
| 752 | |
---|
| 753 | IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0 |
---|
| 754 | |
---|
| 755 | ecf(i,j)=ecf(j,i) |
---|
| 756 | ENDDO |
---|
| 757 | ENDDO |
---|
| 758 | |
---|
| 759 | RETURN |
---|
| 760 | END SUBROUTINE turb_enhance_eff |
---|
| 761 | |
---|
[825] | 762 | END MODULE lpm_collision_kernels_mod |
---|