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