Changeset 828 for palm/trunk/SOURCE/lpm_collision_kernels.f90
 Timestamp:
 Feb 21, 2012 12:00:36 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/lpm_collision_kernels.f90
r826 r828 1 MODULE lpm_collision_kernels_mod1 MODULE lpm_collision_kernels_mod 2 2 3 3 !! 4 4 ! Current revisions: 5 5 !  6 ! 6 ! code has been completely reformatted, routine colker renamed 7 ! recalculate_kernel, 8 ! routine init_kernels added, radius is now communicated to the collision 9 ! routines by array radclass 10 ! 11 ! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4 7 12 ! 8 13 ! Former revisions: … … 24 29 ! Description: 25 30 !  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 (...) 31 ! This module calculates collision efficiencies either due to pure gravitational 32 ! effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 24862507) or 33 ! including the effects of (SGS) turbulence (Wang kernel, see Wang and 34 ! Grabowski, 2009: Atmos. Sci. Lett., 10, 18). The original code has been 35 ! provided by L.P. Wang but is substantially reformatted and speed optimized 36 ! here. 37 ! 38 ! ATTENTION: 39 ! Physical quantities (like g, densities, etc.) used in this module still 40 ! have to be adjusted to those values used in the main PALM code. 41 ! Also, quantities in CGSunits should be converted to SIunits eventually. 29 42 !! 30 43 … … 33 46 USE constants 34 47 USE particle_attributes 48 USE pegrid 49 35 50 36 51 IMPLICIT NONE … … 38 53 PRIVATE 39 54 40 PUBLIC colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi 41 42 INTEGER, SAVE :: ip, jp, kp, pend, pstart, psum 43 44 REAL, SAVE :: epsilon, eps2, urms, urms2 45 46 REAL, DIMENSION(:), ALLOCATABLE, SAVE :: winf 47 REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: ec, ecf, gck 55 PUBLIC ckernel, init_kernels, rclass_lbound, rclass_ubound, & 56 recalculate_kernel 57 58 REAL :: epsilon, eps2, rclass_lbound, rclass_ubound, urms, urms2 59 60 REAL, DIMENSION(:), ALLOCATABLE :: epsclass, radclass, winf 61 REAL, DIMENSION(:,:), ALLOCATABLE :: ec, ecf, gck, hkernel, hwratio 62 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ckernel 63 64 SAVE 48 65 49 66 ! 50 67 ! Public interfaces 51 INTERFACE colker 52 MODULE PROCEDURE colker 53 END INTERFACE colker 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 63 INTERFACE phi 64 MODULE PROCEDURE phi 65 END INTERFACE phi 66 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 68 INTERFACE init_kernels 69 MODULE PROCEDURE init_kernels 70 END INTERFACE init_kernels 71 72 INTERFACE recalculate_kernel 73 MODULE PROCEDURE recalculate_kernel 74 END INTERFACE recalculate_kernel 75 78 76 79 77 CONTAINS 80 78 81 !! 82 ! SUBROUTINE for calculation of collision kernel 83 !! 84 SUBROUTINE colker( i1, j1, k1, kernel ) 79 80 SUBROUTINE init_kernels 81 !! 82 ! Initialization of the collision efficiency matrix with fixed radius and 83 ! dissipation classes, calculated at simulation start only. 84 !! 85 86 IMPLICIT NONE 87 88 INTEGER :: i, j, k 89 90 91 ! 92 ! Calculate collision efficiencies for fixed radius and dissipation 93 ! classes 94 IF ( collision_kernel(6:9) == 'fast' ) THEN 95 96 ALLOCATE( ckernel(1:radius_classes,1:radius_classes, & 97 0:dissipation_classes), epsclass(1:dissipation_classes), & 98 radclass(1:radius_classes) ) 99 100 ! 101 ! Calculate the radius class bounds with logarithmic distances 102 ! in the interval [1.0E6, 2.0E4] m 103 rclass_lbound = LOG( 1.0E6 ) 104 rclass_ubound = LOG( 2.0E4 ) 105 radclass(1) = 1.0E6 106 DO i = 2, radius_classes 107 radclass(i) = EXP( rclass_lbound + & 108 ( rclass_ubound  rclass_lbound ) * ( i1.0 ) /& 109 ( radius_classes  1.0 ) ) 110 ! IF ( myid == 0 ) THEN 111 ! PRINT*, 'i=', i, ' r = ', radclass(i)*1.0E6 112 ! ENDIF 113 ENDDO 114 ! 115 ! Collision routines expect radius to be in cm 116 radclass = radclass * 100.0 117 118 ! 119 ! Set the class bounds for dissipation in interval [0, 1000] cm**2/s**3 120 DO i = 1, dissipation_classes 121 epsclass(i) = 1000.0 * REAL( i ) / dissipation_classes 122 ! IF ( myid == 0 ) THEN 123 ! PRINT*, 'i=', i, ' eps = ', epsclass(i) 124 ! ENDIF 125 ENDDO 126 ! 127 ! Calculate collision efficiencies of the Wang/ayala kernel 128 ALLOCATE( ec(1:radius_classes,1:radius_classes), & 129 ecf(1:radius_classes,1:radius_classes), & 130 gck(1:radius_classes,1:radius_classes), & 131 winf(1:radius_classes) ) 132 133 DO k = 1, dissipation_classes 134 135 epsilon = epsclass(k) 136 urms = 202.0 * ( epsilon / 400.0 )**( 1.0 / 3.0 ) 137 138 CALL turbsd 139 CALL turb_enhance_eff 140 CALL effic 141 142 DO j = 1, radius_classes 143 DO i = 1, radius_classes 144 ckernel(i,j,k) = ec(i,j) * gck(i,j) * ecf(i,j) 145 ENDDO 146 ENDDO 147 148 ENDDO 149 150 ! 151 ! Calculate collision efficiencies of the Hall kernel 152 ALLOCATE( hkernel(1:radius_classes,1:radius_classes), & 153 hwratio(1:radius_classes,1:radius_classes) ) 154 155 CALL fallg 156 CALL effic 157 158 DO j = 1, radius_classes 159 DO i = 1, radius_classes 160 hkernel(i,j) = pi * ( radclass(j) + radclass(i) )**2 & 161 * ec(i,j) * ABS( winf(j)  winf(i) ) 162 ckernel(i,j,0) = hkernel(i,j) ! hall kernel stored on index 0 163 ENDDO 164 ENDDO 165 166 ! 167 ! Test output of efficiencies 168 IF ( j == 1 ) THEN 169 170 PRINT*, '*** Hall kernel' 171 WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes ) 172 DO j = 1, radius_classes 173 WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j), ( hkernel(i,j), i = 1,radius_classes ) 174 ENDDO 175 176 DO k = 1, dissipation_classes 177 DO i = 1, radius_classes 178 DO j = 1, radius_classes 179 IF ( hkernel(i,j) == 0.0 ) THEN 180 hwratio(i,j) = 9999999.9 181 ELSE 182 hwratio(i,j) = ckernel(i,j,k) / hkernel(i,j) 183 ENDIF 184 ENDDO 185 ENDDO 186 187 PRINT*, '*** epsilon = ', epsclass(k) 188 WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes ) 189 DO j = 1, radius_classes 190 ! WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( ckernel(i,j,k), i = 1,radius_classes ) 191 WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( hwratio(i,j), i = 1,radius_classes ) 192 ENDDO 193 ENDDO 194 195 ENDIF 196 197 DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf ) 198 199 ckernel = ckernel * 1.0E6 ! kernel is needed in m**3/s 200 201 ELSEIF( collision_kernel == 'hall' .OR. collision_kernel == 'wang' ) & 202 THEN 203 ! 204 ! Initial settings for Hall and WangKernel 205 ! To be done: move here parts from turbsd, fallg, ecoll, etc. 206 ENDIF 207 208 END SUBROUTINE init_kernels 209 210 211 !! 212 ! Calculation of collision kernels during each timestep and for each grid box 213 !! 214 SUBROUTINE recalculate_kernel( i1, j1, k1 ) 85 215 86 216 USE arrays_3d … … 94 224 IMPLICIT NONE 95 225 96 INTEGER :: i, i1, j, j1, k1 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 100 ! CALL cpu_log( log_point_s(46), 'colker', 'start' ) 101 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 110 ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) ) 111 112 IF ( wang_kernel ) THEN 113 114 ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) ) 115 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 ) 118 119 IF ( epsilon <= 0.001 ) THEN 120 121 CALL fallg 122 CALL effic 123 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) ) 129 ENDDO 226 INTEGER :: i, i1, j, j1, k1, pend, pstart 227 228 229 pstart = prt_start_index(k1,j1,i1) 230 pend = prt_start_index(k1,j1,i1) + prt_count(k1,j1,i1)  1 231 radius_classes = prt_count(k1,j1,i1) 232 233 ALLOCATE( ec(1:radius_classes,1:radius_classes), & 234 radclass(1:radius_classes), winf(1:radius_classes) ) 235 236 ! 237 ! Store particle radii on the radclass array. Collision routines 238 ! expect radii to be in cm. 239 radclass(1:radius_classes) = particles(pstart:pend)%radius * 100.0 240 241 epsilon = diss(k1,j1,i1) * 1.0E4 ! dissipation rate in cm**2/s**3 242 urms = 202.0 * ( epsilon / 400.0 )**( 0.33333333333 ) 243 244 IF ( wang_kernel .AND. epsilon > 0.001 ) THEN 245 ! 246 ! Call routines to calculate efficiencies for the Wang kernel 247 ALLOCATE( gck(1:radius_classes,1:radius_classes), & 248 ecf(1:radius_classes,1:radius_classes) ) 249 250 CALL turbsd 251 CALL turb_enhance_eff 252 CALL effic 253 254 DO j = 1, radius_classes 255 DO i = 1, radius_classes 256 ckernel(pstart+i1,pstart+j1,1) = ec(i,j) * gck(i,j) * ecf(i,j) 130 257 ENDDO 131 132 ELSE 133 134 CALL turbsd 135 CALL turb_enhance_eff 136 CALL effic 137 138 DO j = pstart, pend, 1 139 DO i = pstart, pend, 1 140 kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j) 141 ENDDO 258 ENDDO 259 260 DEALLOCATE( gck, ecf ) 261 262 ELSE 263 ! 264 ! Call routines to calculate efficiencies for the Hall kernel 265 CALL fallg 266 CALL effic 267 268 DO j = 1, radius_classes 269 DO i = 1, radius_classes 270 ckernel(pstart+i1,pstart+j1,1) = pi * & 271 ( radclass(j) + radclass(i) )**2 & 272 * ec(i,j) * ABS( winf(j)  winf(i) ) 142 273 ENDDO 143 144 ENDIF145 146 DEALLOCATE(gck, ecf)147 148 ELSE149 150 ! CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )151 CALL fallg152 ! CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )153 ! CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )154 CALL effic155 ! CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )156 157 DO j = pstart, pend158 DO i = pstart, pend159 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) )162 ENDDO163 274 ENDDO 164 275 165 276 ENDIF 166 277 167 DEALLOCATE( ec, winf ) 168 169 ! CALL cpu_log( log_point_s(46), 'colker', 'stop' ) 170 171 END SUBROUTINE colker 172 173 !! 174 ! SUBROUTINE for calculation of w, g and gck 278 ckernel = ckernel * 1.0E6 ! kernel is needed in m**3/s 279 280 DEALLOCATE( ec, radclass, winf ) 281 282 END SUBROUTINE recalculate_kernel 283 284 285 !! 286 ! Calculation of gck 287 ! This is from Aayala 2008b, page 37ff. 288 ! Necessary input parameters: water density, radii of droplets, air density, 289 ! air viscosity, turbulent dissipation rate, taylor microscale reynolds number, 290 ! gravitational acceleration > to be replaced by PALM parameters 175 291 !! 176 292 SUBROUTINE turbsd 177 ! from Aayala 2008b, page 37ff, necessary input parameter water density, radii178 ! of droplets, air density, air viscosity, turbulent dissipation rate,179 ! taylor microscale reynolds number, gravitational acceleration180 293 181 294 USE constants … … 186 299 IMPLICIT NONE 187 300 188 REAL :: Relamda, & 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 194 REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens 195 196 REAL, DIMENSION(pstart:pend) :: St, tau 301 INTEGER :: i, j 197 302 198 303 LOGICAL, SAVE :: first = .TRUE. 199 304 200 INTEGER :: i, j 201 202 ! 203 ! Initial assignment of constants 305 REAL :: ao, ao_gr, bbb, be, b1, b2, ccc, c1, c1_gr, c2, d1, d2, eta, & 306 e1, e2, fao_gr, fr, grfin, lambda, lambda_re, lf, rc, rrp, & 307 sst, tauk, tl, t2, tt, t1, vk, vrms1xy, vrms2xy, v1, v1v2xy, & 308 v1xysq, v2, v2xysq, wrfin, wrgrav2, wrtur2xy, xx, yy, z 309 310 REAL, SAVE :: airdens, airvisc, anu, gravity, waterdens 311 312 REAL, DIMENSION(1:radius_classes) :: st, tau 313 314 315 ! 316 ! Initial assignment of constants 204 317 IF ( first ) THEN 205 318 206 319 first = .FALSE. 207 airvisc = 0.1818 !dynamic viscosity in mg/cm*s208 airdens = 1.2250 !air density in mg/cm**3209 waterdens = 1000.0 !water density in mg/cm**3210 gravity = 980.6650 !in cm/s**2320 airvisc = 0.1818 ! dynamic viscosity in mg/cm*s 321 airdens = 1.2250 ! air density in mg/cm**3 322 waterdens = 1000.0 ! water density in mg/cm**3 323 gravity = 980.6650 ! in cm/s**2 211 324 anu = airvisc/airdens ! kinetic viscosity in cm**2/s 212 325 213 326 ENDIF 214 327 215 Relamda = urms**2*sqrt(15.0/epsilon/anu) 216 217 Tl = urms**2/epsilon !in s 218 Lf = 0.5 * (urms**3)/epsilon !in cm 219 220 tauk = (anu/epsilon)**0.5 !in s 221 eta = (anu**3/epsilon)**0.25 !in cm 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 233 tau(i) = winf(i)/gravity !in s 234 St(i) = tau(i)/tauk 328 lambda = urms * SQRT( 15.0 * anu / epsilon ) ! in cm 329 lambda_re = urms**2 * SQRT( 15.0 / epsilon / anu ) 330 tl = urms**2 / epsilon ! in s 331 lf = 0.5 * urms**3 / epsilon ! in cm 332 tauk = SQRT( anu / epsilon ) ! in s 333 eta = ( anu**3 / epsilon )**0.25 ! in cm 334 vk = eta / tauk ! in cm/s 335 336 ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re ) 337 tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk ! in s 338 339 CALL fallg ! gives winf in cm/s 340 341 DO i = 1, radius_classes 342 tau(i) = winf(i) / gravity ! in s 343 st(i) = tau(i) / tauk 235 344 ENDDO 236 345 237 !***** TO CALCULATE wr ******************************** 238 !from Aayala 2008b, page 38f 239 240 z = tt/Tl 241 be = sqrt(2.0)*lambda/Lf 242 243 bbb = sqrt(1.02.0*be**2) 244 d1 = (1.+bbb)/2.0/bbb 245 e1 = Lf*(1.0+bbb)/2.0 !in cm 246 d2 = (1.0bbb)/2.0/bbb 247 e2 = Lf*(1.0bbb)/2.0 !in cm 248 249 ccc = sqrt(1.02.0*z**2) 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 256 v1 = winf(i) !in cm/s 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 261 v2 = winf(j) !in cm/s 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) 268 v1xysq = v1xysq * urms**2/t1 !in cm**2/s**2 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) 276 v2xysq = v2xysq * urms**2/t2 !in cm**2/s**2 277 278 vrms2xy= sqrt(v2xysq) !in cm/s 279 280 IF(winf(i).ge.winf(j)) THEN 346 ! 347 ! Calculate wr (from Aayala 2008b, page 38f) 348 z = tt / tl 349 be = SQRT( 2.0 ) * lambda / lf 350 bbb = SQRT( 1.0  2.0 * be**2 ) 351 d1 = ( 1.0 + bbb ) / ( 2.0 * bbb ) 352 e1 = lf * ( 1.0 + bbb ) * 0.5 ! in cm 353 d2 = ( 1.0  bbb ) * 0.5 / bbb 354 e2 = lf * ( 1.0  bbb ) * 0.5 ! in cm 355 ccc = SQRT( 1.0  2.0 * z**2 ) 356 b1 = ( 1.0 + ccc ) * 0.5 / ccc 357 c1 = tl * ( 1.0 + ccc ) * 0.5 ! in s 358 b2 = ( 1.0  ccc ) * 0.5 / ccc 359 c2 = tl * ( 1.0  ccc ) * 0.5 ! in s 360 361 DO i = 1, radius_classes 362 363 v1 = winf(i) ! in cm/s 364 t1 = tau(i) ! in s 365 366 DO j = 1, i 367 rrp = radclass(i) + radclass(j) ! radius in cm 368 v2 = winf(j) ! in cm/s 369 t2 = tau(j) ! in s 370 371 v1xysq = b1 * d1 * phi(c1,e1,v1,t1)  b1 * d2 * phi(c1,e2,v1,t1) & 372  b2 * d1 * phi(c2,e1,v1,t1) + b2 * d2 * phi(c2,e2,v1,t1) 373 v1xysq = v1xysq * urms**2 / t1 ! in cm**2/s**2 374 vrms1xy = SQRT( v1xysq ) ! in cm/s 375 376 v2xysq = b1 * d1 * phi(c1,e1,v2,t2)  b1 * d2 * phi(c1,e2,v2,t2) & 377  b2 * d1 * phi(c2,e1,v2,t2) + b2 * d2 * phi(c2,e2,v2,t2) 378 v2xysq = v2xysq * urms**2 / t2 ! in cm**2/s**2 379 vrms2xy = SQRT( v2xysq ) ! in cm/s 380 381 IF ( winf(i) >= winf(j) ) THEN 281 382 v1 = winf(i) 282 383 t1 = tau(i) … … 290 391 ENDIF 291 392 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) 297 v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j) !in cm**2/s**2 298 299 wrtur2xy=vrms1xy**2 + vrms2xy**2  2.*v1v2xy !in cm**2/s**2 300 301 IF (wrtur2xy.le.0.0) wrtur2xy=0.0 302 303 wrgrav2=pi/8.*(winf(j)winf(i))**2 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) 393 v1v2xy = b1 * d1 * zhi(c1,e1,v1,t1,v2,t2)  & 394 b1 * d2 * zhi(c1,e2,v1,t1,v2,t2)  & 395 b2 * d1 * zhi(c2,e1,v1,t1,v2,t2) + & 396 b2 * d2* zhi(c2,e2,v1,t1,v2,t2) 397 fr = d1 * EXP( rrp / e1 )  d2 * EXP( rrp / e2 ) 398 v1v2xy = v1v2xy * fr * urms**2 / tau(i) / tau(j) ! in cm**2/s**2 399 wrtur2xy = vrms1xy**2 + vrms2xy**2  2.0 * v1v2xy ! in cm**2/s**2 400 IF ( wrtur2xy < 0.0 ) wrtur2xy = 0.0 401 wrgrav2 = pi / 8.0 * ( winf(j)  winf(i) )**2 402 wrfin = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) ) ! in cm/s 403 404 ! 405 ! Calculate gr 406 IF ( st(j) > st(i) ) THEN 407 sst = st(j) 312 408 ELSE 313 SSt = St(i)409 sst = st(i) 314 410 ENDIF 315 411 316 XX = 0.1988*SSt**4 + 1.5275*SSt**3 & 317 4.2942*SSt**2 + 5.3406*SSt 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 325 ao_gr = ao + (pi/8.)*(gravity/(vk/tauk))**2 326 fao_gr = 20.115 * (ao_gr/Relamda)**0.5 327 rc = sqrt( fao_gr * abs(St(j)St(i)) ) * eta !in cm 328 329 grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.) 330 IF (grFIN.lt.1.0) grFIN = 1.0 331 332 gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN ! in cm**3/s 412 xx = 0.1988 * sst**4 + 1.5275 * sst**3  4.2942 * sst**2 + & 413 5.3406 * sst 414 IF ( xx < 0.0 ) xx = 0.0 415 yy = 0.1886 * EXP( 20.306 / lambda_re ) 416 417 c1_gr = xx / ( gravity / vk * tauk )**yy 418 419 ao_gr = ao + ( pi / 8.0) * ( gravity / vk * tauk )**2 420 fao_gr = 20.115 * SQRT( ao_gr / lambda_re ) 421 rc = SQRT( fao_gr * ABS( st(j)  st(i) ) ) * eta ! in cm 422 423 grfin = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5 ) 424 IF ( grfin < 1.0 ) grfin = 1.0 425 426 gck(i,j) = 2.0 * pi * rrp**2 * wrfin * grfin ! in cm**3/s 333 427 gck(j,i) = gck(i,j) 334 428 … … 336 430 ENDDO 337 431 338 END SUBROUTINE TurbSD 339 340 !! 341 ! PHI as a function 342 !! 343 REAL FUNCTION PHI(a,b,vsett,tau0) 432 END SUBROUTINE turbsd 433 434 435 !! 436 ! phi as a function 437 !! 438 REAL FUNCTION phi( a, b, vsett, tau0 ) 344 439 345 440 IMPLICIT NONE 346 441 347 REAL :: a, aa1, b, vsett, tau0348 349 aa1 = 1. /tau0 + 1./a + vsett/b350 351 PHI = 1./aa1  vsett/2.0/b/aa1**2 !in s 352 353 END FUNCTION PHI 354 355 !! 356 ! ZETAas a function357 !! 358 REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2)442 REAL :: a, aa1, b, tau0, vsett 443 444 aa1 = 1.0 / tau0 + 1.0 / a + vsett / b 445 phi = 1.0 / aa1  0.5 * vsett / b / aa1**2 ! in s 446 447 END FUNCTION phi 448 449 450 !! 451 ! zeta as a function 452 !! 453 REAL FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 ) 359 454 360 455 IMPLICIT NONE 361 456 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 = (vsett1vsett2)/b + 1./tau1 + 1./tau2 367 aa4 = (vsett2/b)**2  (1./tau2 + 1./a)**2 368 aa5 = vsett2/b + 1./tau2 + 1./a 369 aa6 = 1./tau1  1./a + (1./tau2 + 1./a) * vsett1/vsett2 370 ZHI = (1./aa1  1./aa2) * (vsett1vsett2)/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) & 373 * 1./2./b/aa3 ! in s**2 374 375 END FUNCTION ZHI 376 377 !! 378 ! SUBROUTINE for calculation of terminal velocity winf 379 !! 380 SUBROUTINE fallg 457 REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, tau1, tau2, vsett1, vsett2 458 459 aa1 = vsett2 / b  1.0 / tau2  1.0 / a 460 aa2 = vsett1 / b + 1.0 / tau1 + 1.0 / a 461 aa3 = ( vsett1  vsett2 ) / b + 1.0 / tau1 + 1.0 / tau2 462 aa4 = ( vsett2 / b )**2  ( 1.0 / tau2 + 1.0 / a )**2 463 aa5 = vsett2 / b + 1.0 / tau2 + 1.0 / a 464 aa6 = 1.0 / tau1  1.0 / a + ( 1.0 / tau2 + 1.0 / a) * vsett1 / vsett2 465 zhi = (1.0 / aa1  1.0 / aa2 ) * ( vsett1  vsett2 ) * 0.5 / b / aa3**2 & 466 + (4.0 / aa4  1.0 / aa5**2  1.0 / aa1**2 ) * vsett2 * 0.5 / b /aa6& 467 + (2.0 * ( b / aa2  b / aa1 )  vsett1 / aa2**2 + vsett2 / aa1**2 )& 468 * 0.5 / b / aa3 ! in s**2 469 470 END FUNCTION zhi 471 472 473 !! 474 ! Calculation of terminal velocity winf 475 !! 476 SUBROUTINE fallg 381 477 382 478 USE constants … … 385 481 USE arrays_3d 386 482 387 IMPLICIT NONE 388 389 INTEGER :: i, j 390 391 LOGICAL, SAVE :: first = .TRUE. 392 393 REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py 394 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 483 IMPLICIT NONE 484 485 INTEGER :: i, j 486 487 LOGICAL, SAVE :: first = .TRUE. 488 489 REAL, SAVE :: cunh, eta, grav, phy, py, rhoa, rhow, sigma, stb, stok, & 490 t0, xlamb 491 492 REAL :: bond, x, xrey, y 493 494 REAL, DIMENSION(1:7), SAVE :: b 495 REAL, DIMENSION(1:6), SAVE :: c 496 497 ! 498 ! Initial assignment of constants 499 IF ( first ) THEN 500 501 first = .FALSE. 502 b = (/ 0.318657E1, 0.992696E0, 0.153193E2, 0.987059E3, & 503 0.578878E3, 0.855176E4, 0.327815E5 /) 504 c = (/ 0.500015E1, 0.523778E1, 0.204914E1, 0.475294E0, & 505 0.542819E1, 0.238449E2 /) 506 507 eta = 1.818E4 ! in poise = g/(cm s) 508 xlamb = 6.62E6 ! in cm 509 rhow = 1.0 ! in g/cm**3 510 rhoa = 1.225E3 ! in g/cm**3 511 grav = 980.665 ! in cm/s**2 512 cunh = 1.257 * xlamb ! in cm 513 t0 = 273.15 ! in K 514 sigma = 76.1  0.155 * ( 293.15  t0 ) ! in N/m = g/s**2 515 stok = 2.0 * grav * ( rhow  rhoa ) / ( 9.0 * eta ) ! in 1/(cm s) 516 stb = 32.0 * rhoa * ( rhow  rhoa) * grav / (3.0 * eta * eta) 517 ! in 1/cm**3 518 phy = sigma**3 * rhoa**2 / ( eta**4 * grav * ( rhow  rhoa ) ) 519 py = phy**( 1.0 / 6.0 ) 520 521 ENDIF 522 523 DO j = 1, radius_classes 524 525 IF ( radclass(j) <= 1.0E3 ) THEN 526 527 winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) ) ! in cm/s 528 529 ELSEIF ( radclass(j) > 1.0E3 .AND. radclass(j) <= 5.35E2 ) THEN 530 531 x = LOG( stb * radclass(j)**3 ) 532 y = 0.0 533 534 DO i = 1, 7 535 y = y + b(i) * x**(i1) 536 ENDDO 537 538 xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y ) 539 winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) ) ! in cm/s 540 541 ELSEIF ( radclass(j) > 5.35E2 ) THEN 542 543 IF ( radclass(j) > 0.35 ) THEN 544 bond = grav * ( rhow  rhoa ) * 0.35**2 / sigma 545 ELSE 546 bond = grav * ( rhow  rhoa ) * radclass(j)**2 / sigma 547 ENDIF 548 549 x = LOG( 16.0 * bond * py / 3.0 ) 550 y = 0.0 551 552 DO i = 1, 6 553 y = y + c(i) * x**(i1) 554 ENDDO 555 556 xrey = py * EXP( y ) 557 558 IF ( radclass(j) > 0.35 ) THEN 559 winf(j) = xrey * eta / ( 2.0 * rhoa * 0.35 ) ! in cm/s 560 ELSE 561 winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) ) ! in cm/s 562 ENDIF 563 564 ENDIF 565 566 ENDDO 567 568 END SUBROUTINE fallg 569 570 571 !! 572 ! Calculation of collision efficencies for the Hall kernel 573 !! 574 SUBROUTINE effic 575 576 USE arrays_3d 577 USE cloud_parameters 578 USE constants 579 USE particle_attributes 580 581 IMPLICIT NONE 582 583 INTEGER :: i, iq, ir, j, k, kk 584 585 INTEGER, DIMENSION(:), ALLOCATABLE :: ira 586 587 LOGICAL, SAVE :: first = .TRUE. 588 589 REAL :: ek, particle_radius, pp, qq, rq 590 591 REAL, DIMENSION(1:21), SAVE :: rat 592 REAL, DIMENSION(1:15), SAVE :: r0 593 REAL, DIMENSION(1:15,1:21), SAVE :: ecoll 594 595 ! 596 ! Initial assignment of constants 597 IF ( first ) THEN 403 598 404 599 first = .FALSE. 405 b = (/0.318657e1,0.992696e0,0.153193e2,0.987059e3, & 406 0.578878e3,0.855176e4,0.327815e5/) 407 c = (/0.500015e1,0.523778e1,0.204914e1,0.475294e0, & 408 0.542819e1, 0.238449e2/) 409 410 eta = 1.818e4 !in poise = g/(cm s) 411 xlamb = 6.62e6 !in cm 412 413 rhow = 1.0 !in g/cm**3 414 rhoa = 1.225e3 !in g/cm**3 415 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 * (rhowrhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3 422 phy = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhowrhoa)) 423 py = phy**(1.0/6.0) 424 425 ENDIF 426 427 !particle radius has to be in cm 428 DO j = pstart, pend 429 430 IF (particles(j)%radius*100.0 .le. 1.e3) THEN 431 432 winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s 433 434 ELSEIF (particles(j)%radius*100.0.gt.1.e3.and.particles(j)%radius*100.0.le.5.35e2) THEN 435 436 x = log(stb*(particles(j)%radius*100.0)**3) 437 y = 0.0 438 439 DO i = 1, 7 440 y = y + b(i) * (x**(i1)) 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.35e2) THEN 447 448 IF (particles(j)%radius*100.0.gt.0.35) THEN 449 bond = grav*(rhowrhoa) * 0.35 * 0.35/sigma 450 ELSE 451 bond = grav*(rhowrhoa)*(particles(j)%radius*100.0)**2/sigma 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**(i1)) 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 465 ELSE 466 winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0) !in cm/s 467 ENDIF 468 469 ENDIF 470 ENDDO 471 RETURN 472 END SUBROUTINE fallg 473 474 !! 475 ! SUBROUTINE for calculation of collision efficencies 476 !! 477 478 SUBROUTINE effic 479 480 USE arrays_3d 481 USE constants 482 USE cloud_parameters 483 USE particle_attributes 484 485 !collision efficiencies of hall kernel 486 IMPLICIT NONE 487 488 INTEGER :: i, ir, iq, j, k, kk 489 490 INTEGER, DIMENSION(:), ALLOCATABLE :: ira 491 492 LOGICAL, SAVE :: first = .TRUE. 493 494 REAL :: ek, particle_radius, pp, qq, rq 495 496 REAL, DIMENSION(1:21), SAVE :: rat 497 REAL, DIMENSION(1:15), SAVE :: r0 498 REAL, DIMENSION(1:15,1:21), SAVE :: ecoll 499 500 ! 501 ! Initial assignment of constants 502 IF ( first ) THEN 503 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/) 509 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 553 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 567 568 ! 569 ! Twodimensional 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 579 ! IF ( rq <= rat(kk) ) THEN 580 ! iq = kk 581 ! EXIT 582 ! ENDIF 583 ! ENDDO 584 585 iq = INT( rq * 20 ) + 1 586 iq = MAX(iq , 2) 587 588 IF ( ir < 16 ) THEN 589 590 IF ( ir >= 2 ) THEN 591 pp = ( ( particles(j)%radius * 1.0E06 )  r0(ir1) ) / & 592 ( r0(ir)  r0(ir1) ) 593 qq = ( rq rat(iq1) ) / ( rat(iq)  rat(iq1) ) 594 ec(j,i) = ( 1.0pp ) * ( 1.0qq ) * ecoll(ir1,iq1) & 595 + pp * ( 1.0qq ) * ecoll(ir,iq1) & 596 + qq * ( 1.0pp ) * ecoll(ir1,iq) & 597 + pp * qq * ecoll(ir,iq) 598 ELSE 599 qq = (rqrat(iq1))/(rat(iq)rat(iq1)) 600 ec(j,i) = (1.qq) * ecoll(1,iq1) + qq * ecoll(1,iq) 601 ENDIF 602 603 ELSE 604 qq = ( rq  rat(iq1) ) / ( rat(iq)  rat(iq1) ) 605 ek = ( 1.0  qq ) * ecoll(15,iq1) + qq * ecoll(15,iq) 606 ec(j,i) = MIN( ek, 1.0 ) 607 ENDIF 608 609 ec(i,j) = ec(j,i) 610 IF ( ec(i,j) < 1.0E20 ) ec(i,j) = 0.0 611 612 ENDDO 613 ENDDO 614 615 DEALLOCATE( ira ) 616 617 END SUBROUTINE effic 618 619 !! 620 ! SUBROUTINE for calculation of enhancement factor collision efficencies 621 !! 622 SUBROUTINE turb_enhance_eff 600 r0 = (/ 6.0, 8.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60., & 601 70.0, 100.0, 150.0, 200.0, 300.0 /) 602 rat = (/ 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, & 603 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, & 604 1.00 /) 605 606 ecoll(:,1) = (/0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, & 607 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001/) 608 ecoll(:,2) = (/0.003, 0.003, 0.003, 0.004, 0.005, 0.005, 0.005, & 609 0.010, 0.100, 0.050, 0.200, 0.500, 0.770, 0.870, 0.970/) 610 ecoll(:,3) = (/0.007, 0.007, 0.007, 0.008, 0.009, 0.010, 0.010, & 611 0.070, 0.400, 0.430, 0.580, 0.790, 0.930, 0.960, 1.000/) 612 ecoll(:,4) = (/0.009, 0.009, 0.009, 0.012, 0.015, 0.010, 0.020, & 613 0.280, 0.600, 0.640, 0.750, 0.910, 0.970, 0.980, 1.000/) 614 ecoll(:,5) = (/0.014, 0.014, 0.014, 0.015, 0.016, 0.030, 0.060, & 615 0.500, 0.700, 0.770, 0.840, 0.950, 0.970, 1.000, 1.000/) 616 ecoll(:,6) = (/0.017, 0.017, 0.017, 0.020, 0.022, 0.060, 0.100, & 617 0.620, 0.780, 0.840, 0.880, 0.950, 1.000, 1.000, 1.000/) 618 ecoll(:,7) = (/0.030, 0.030, 0.024, 0.022, 0.032, 0.062, 0.200, & 619 0.680, 0.830, 0.870, 0.900, 0.950, 1.000, 1.000, 1.000/) 620 ecoll(:,8) = (/0.025, 0.025, 0.025, 0.036, 0.043, 0.130, 0.270, & 621 0.740, 0.860, 0.890, 0.920, 1.000, 1.000, 1.000, 1.000/) 622 ecoll(:,9) = (/0.027, 0.027, 0.027, 0.040, 0.052, 0.200, 0.400, & 623 0.780, 0.880, 0.900, 0.940, 1.000, 1.000, 1.000, 1.000/) 624 ecoll(:,10)= (/0.030, 0.030, 0.030, 0.047, 0.064, 0.250, 0.500, & 625 0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/) 626 ecoll(:,11)= (/0.040, 0.040, 0.033, 0.037, 0.068, 0.240, 0.550, & 627 0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/) 628 ecoll(:,12)= (/0.035, 0.035, 0.035, 0.055, 0.079, 0.290, 0.580, & 629 0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/) 630 ecoll(:,13)= (/0.037, 0.037, 0.037, 0.062, 0.082, 0.290, 0.590, & 631 0.780, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/) 632 ecoll(:,14)= (/0.037, 0.037, 0.037, 0.060, 0.080, 0.290, 0.580, & 633 0.770, 0.890, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/) 634 ecoll(:,15)= (/0.037, 0.037, 0.037, 0.041, 0.075, 0.250, 0.540, & 635 0.760, 0.880, 0.920, 0.950, 1.000, 1.000, 1.000, 1.000/) 636 ecoll(:,16)= (/0.037, 0.037, 0.037, 0.052, 0.067, 0.250, 0.510, & 637 0.770, 0.880, 0.930, 0.970, 1.000, 1.000, 1.000, 1.000/) 638 ecoll(:,17)= (/0.037, 0.037, 0.037, 0.047, 0.057, 0.250, 0.490, & 639 0.770, 0.890, 0.950, 1.000, 1.000, 1.000, 1.000, 1.000/) 640 ecoll(:,18)= (/0.036, 0.036, 0.036, 0.042, 0.048, 0.230, 0.470, & 641 0.780, 0.920, 1.000, 1.020, 1.020, 1.020, 1.020, 1.020/) 642 ecoll(:,19)= (/0.040, 0.040, 0.035, 0.033, 0.040, 0.112, 0.450, & 643 0.790, 1.010, 1.030, 1.040, 1.040, 1.040, 1.040, 1.040/) 644 ecoll(:,20)= (/0.033, 0.033, 0.033, 0.033, 0.033, 0.119, 0.470, & 645 0.950, 1.300, 1.700, 2.300, 2.300, 2.300, 2.300, 2.300/) 646 ecoll(:,21)= (/0.027, 0.027, 0.027, 0.027, 0.027, 0.125, 0.520, & 647 1.400, 2.300, 3.000, 4.000, 4.000, 4.000, 4.000, 4.000/) 648 ENDIF 649 650 ! 651 ! Calculate the radius class index of particles with respect to array r 652 ALLOCATE( ira(1:radius_classes) ) 653 DO j = 1, radius_classes 654 particle_radius = radclass(j) * 1.0E4 ! in microm 655 DO k = 1, 15 656 IF ( particle_radius < r0(k) ) THEN 657 ira(j) = k 658 EXIT 659 ENDIF 660 ENDDO 661 IF ( particle_radius >= r0(15) ) ira(j) = 16 662 ENDDO 663 664 ! 665 ! Twodimensional linear interpolation of the collision efficiency. 666 ! Radius has to be in Âµm 667 DO j = 1, radius_classes 668 DO i = 1, j 669 670 ir = ira(j) 671 rq = radclass(i) / radclass(j) 672 iq = INT( rq * 20 ) + 1 673 iq = MAX( iq , 2) 674 675 IF ( ir < 16 ) THEN 676 IF ( ir >= 2 ) THEN 677 pp = ( ( radclass(j) * 1.0E04 )  r0(ir1) ) / & 678 ( r0(ir)  r0(ir1) ) 679 qq = ( rq rat(iq1) ) / ( rat(iq)  rat(iq1) ) 680 ec(j,i) = ( 1.0pp ) * ( 1.0qq ) * ecoll(ir1,iq1) & 681 + pp * ( 1.0qq ) * ecoll(ir,iq1) & 682 + qq * ( 1.0pp ) * ecoll(ir1,iq) & 683 + pp * qq * ecoll(ir,iq) 684 ELSE 685 qq = ( rq  rat(iq1) ) / ( rat(iq)  rat(iq1) ) 686 ec(j,i) = (1.0qq) * ecoll(1,iq1) + qq * ecoll(1,iq) 687 ENDIF 688 ELSE 689 qq = ( rq  rat(iq1) ) / ( rat(iq)  rat(iq1) ) 690 ek = ( 1.0  qq ) * ecoll(15,iq1) + qq * ecoll(15,iq) 691 ec(j,i) = MIN( ek, 1.0 ) 692 ENDIF 693 694 ec(i,j) = ec(j,i) 695 IF ( ec(i,j) < 1.0E20 ) ec(i,j) = 0.0 696 697 ENDDO 698 ENDDO 699 700 DEALLOCATE( ira ) 701 702 END SUBROUTINE effic 703 704 705 !! 706 ! Calculation of enhancement factor for collision efficencies due to turbulence 707 !! 708 SUBROUTINE turb_enhance_eff 623 709 624 710 USE constants … … 627 713 USE arrays_3d 628 714 629 IMPLICIT NONE 630 631 INTEGER :: i, ik, ir, iq, j, k, kk 632 633 INTEGER, DIMENSION(:), ALLOCATABLE :: ira 634 635 REAL :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3 636 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 652 ! 100 cm^2/s^3 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 /) 664 665 ! 400 cm^2/s^3 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 /) 677 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 695 ! twodimensional linear interpolation of the collision efficiency 696 DO j = pstart, pend 697 DO i = pstart, j 698 699 ir = ira(j) 700 701 rq = particles(i)%radius/particles(j)%radius 702 703 DO kk = 2, 11 704 IF ( rq <= rat(kk) ) THEN 705 iq = kk 706 EXIT 707 ENDIF 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(ir1))/(r0(ir)r0(ir1)) 716 qq = (rqrat(iq1))/(rat(iq)rat(iq1)) 717 y2= (1.pp)*(1.qq)*ecoll_100(ir1,iq1)+ & 718 pp*(1.qq)*ecoll_100(ir,iq1)+ & 719 qq*(1.pp)*ecoll_100(ir1,iq)+ & 720 pp*qq*ecoll_100(ir,iq) 721 y3= (1.pp)*(1.qq)*ecoll_400(ir1,iq1)+ & 722 pp*(1.qq)*ecoll_400(ir,iq1)+ & 723 qq*(1.pp)*ecoll_400(ir1,iq)+ & 724 pp*qq*ecoll_400(ir,iq) 725 ELSE 726 qq = (rqrat(iq1))/(rat(iq)rat(iq1)) 727 y2= (1.qq)*ecoll_100(1,iq1)+qq*ecoll_100(1,iq) 728 y3= (1.qq)*ecoll_400(1,iq1)+qq*ecoll_400(1,iq) 729 ENDIF 730 ELSE 731 qq = (rqrat(iq1))/(rat(iq)rat(iq1)) 732 y2 = (1.qq) * ecoll_100(7,iq1) + qq * ecoll_100(7,iq) 733 y3 = (1.qq) * ecoll_400(7,iq1) + 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) = (epsilon100.)/(0.100.) * y1 & 743 + (epsilon0.)/(100.0.) * y2 744 ELSE IF(epsilon.le.600.)THEN 745 ecf(j,i) = (epsilon400.)/(100.400.) * y2 & 746 + (epsilon100.)/(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 715 IMPLICIT NONE 716 717 INTEGER :: i, ik, iq, ir, j, k, kk 718 719 INTEGER, DIMENSION(:), ALLOCATABLE :: ira 720 721 REAL :: particle_radius, pp, qq, rq, x1, x2, x3, y1, y2, y3 722 723 LOGICAL, SAVE :: first = .TRUE. 724 725 REAL, DIMENSION(1:11), SAVE :: rat 726 REAL, DIMENSION(1:7), SAVE :: r0 727 REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400 728 729 ! 730 ! Initial assignment of constants 731 IF ( first ) THEN 732 733 first = .FALSE. 734 735 r0 = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 100.0 /) 736 rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /) 737 ! 738 ! In 100 cm^2/s^3 739 ecoll_100(:,1) = (/1.74, 1.74, 1.773, 1.49, 1.207, 1.207, 1.0 /) 740 ecoll_100(:,2) = (/1.46, 1.46, 1.421, 1.245, 1.069, 1.069, 1.0 /) 741 ecoll_100(:,3) = (/1.32, 1.32, 1.245, 1.123, 1.000, 1.000, 1.0 /) 742 ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025, 1.025, 1.0 /) 743 ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056, 1.056, 1.0 /) 744 ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028, 1.028, 1.0 /) 745 ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046, 1.046, 1.0 /) 746 ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029, 1.029, 1.0 /) 747 ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021, 1.021, 1.0 /) 748 ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088, 1.088, 1.0 /) 749 ecoll_100(:,11)= (/20.3, 20.3, 14.6 , 8.61, 2.60, 2.60 , 1.0 /) 750 ! 751 ! In 400 cm^2/s^3 752 ecoll_400(:,1) = (/4.976, 4.976, 3.593, 2.519, 1.445, 1.445, 1.0 /) 753 ecoll_400(:,2) = (/2.984, 2.984, 2.181, 1.691, 1.201, 1.201, 1.0 /) 754 ecoll_400(:,3) = (/1.988, 1.988, 1.475, 1.313, 1.150, 1.150, 1.0 /) 755 ecoll_400(:,4) = (/1.490, 1.490, 1.187, 1.156, 1.126, 1.126, 1.0 /) 756 ecoll_400(:,5) = (/1.249, 1.249, 1.088, 1.090, 1.092, 1.092, 1.0 /) 757 ecoll_400(:,6) = (/1.139, 1.139, 1.130, 1.091, 1.051, 1.051, 1.0 /) 758 ecoll_400(:,7) = (/1.220, 1.220, 1.190, 1.138, 1.086, 1.086, 1.0 /) 759 ecoll_400(:,8) = (/1.325, 1.325, 1.267, 1.165, 1.063, 1.063, 1.0 /) 760 ecoll_400(:,9) = (/1.716, 1.716, 1.345, 1.223, 1.100, 1.100, 1.0 /) 761 ecoll_400(:,10)= (/3.788, 3.788, 1.501, 1.311, 1.120, 1.120, 1.0 /) 762 ecoll_400(:,11)= (/36.52, 36.52, 19.16, 22.80, 26.0, 26.0, 1.0 /) 763 764 ENDIF 765 766 ! 767 ! Calculate the radius class index of particles with respect to array r0 768 ALLOCATE( ira(1:radius_classes) ) 769 770 DO j = 1, radius_classes 771 particle_radius = radclass(j) * 1.0E4 ! in microm 772 DO k = 1, 7 773 IF ( particle_radius < r0(k) ) THEN 774 ira(j) = k 775 EXIT 776 ENDIF 777 ENDDO 778 IF ( particle_radius >= r0(7) ) ira(j) = 8 779 ENDDO 780 781 ! 782 ! Twodimensional linear interpolation of the collision efficiencies 783 DO j = 1, radius_classes 784 DO i = 1, j 785 786 ir = ira(j) 787 rq = radclass(i) / radclass(j) 788 789 DO kk = 2, 11 790 IF ( rq <= rat(kk) ) THEN 791 iq = kk 792 EXIT 793 ENDIF 794 ENDDO 795 796 y1 = 1.0 ! 0 cm2/s3 797 ! 798 ! 100 cm2/s3, 400 cm2/s3 799 IF ( ir < 8 ) THEN 800 IF ( ir >= 2 ) THEN 801 pp = ( radclass(j)*1.0E4  r0(ir1) ) / ( r0(ir)  r0(ir1) ) 802 qq = ( rq  rat(iq1) ) / ( rat(iq)  rat(iq1) ) 803 y2 = ( 1.0pp ) * ( 1.0qq ) * ecoll_100(ir1,iq1) + & 804 pp * ( 1.0qq ) * ecoll_100(ir,iq1) + & 805 qq * ( 10.pp ) * ecoll_100(ir1,iq) + & 806 pp * qq * ecoll_100(ir,iq) 807 y3 = ( 1.0pp ) * ( 1.0qq ) * ecoll_400(ir1,iq1) + & 808 pp * ( 1.0qq ) * ecoll_400(ir,iq1) + & 809 qq * ( 1.0pp ) * ecoll_400(ir1,iq) + & 810 pp * qq * ecoll_400(ir,iq) 811 ELSE 812 qq = ( rq  rat(iq1) ) / ( rat(iq)  rat(iq1) ) 813 y2 = ( 1.0qq ) * ecoll_100(1,iq1) + qq * ecoll_100(1,iq) 814 y3 = ( 1.0qq ) * ecoll_400(1,iq1) + qq * ecoll_400(1,iq) 815 ENDIF 816 ELSE 817 qq = ( rq  rat(iq1) ) / ( rat(iq)  rat(iq1) ) 818 y2 = ( 1.0qq ) * ecoll_100(7,iq1) + qq * ecoll_100(7,iq) 819 y3 = ( 1.0qq ) * ecoll_400(7,iq1) + qq * ecoll_400(7,iq) 820 ENDIF 821 ! 822 ! Linear interpolation of dissipation rate in cm2/s3 823 IF ( epsilon <= 100.0 ) THEN 824 ecf(j,i) = ( epsilon  100.0 ) / ( 0.0  100.0 ) * y1 & 825 + ( epsilon  0.0 ) / ( 100.0  0.0 ) * y2 826 ELSEIF ( epsilon <= 600.0 ) THEN 827 ecf(j,i) = ( epsilon  400.0 ) / ( 100.0  400.0 ) * y2 & 828 + ( epsilon  100.0 ) / ( 400.0  100.0 ) * y3 829 ELSE 830 ecf(j,i) = ( 600.0  400.0 ) / ( 100.0  400.0 ) * y2 & 831 + ( 600.0  100.0 ) / ( 400.0  100.0 ) * y3 832 ENDIF 833 834 IF ( ecf(j,i) < 1.0 ) ecf(j,i) = 1.0 835 836 ecf(i,j) = ecf(j,i) 837 838 ENDDO 839 ENDDO 840 841 END SUBROUTINE turb_enhance_eff 761 842 762 843 END MODULE lpm_collision_kernels_mod
Note: See TracChangeset
for help on using the changeset viewer.