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