MODULE wang_kernel_mod !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: wang_kernel.f90 800 2011-12-21 18:13:51Z suehring $ ! ! 799 2011-12-21 17:48:03Z franke ! speed optimizations and formatting ! Bugfix: iq=1 is not allowed (routine effic) ! Bugfix: replaced stop by ec=0.0 in case of very small ec (routine effic) ! ! 790 2011-11-29 03:11:20Z raasch ! initial revision ! ! Description: ! ------------ ! This routine calculates the effect of (SGS) turbulence on the collision ! efficiency of droplets. ! It is based on the original kernel developed by Wang (...) !------------------------------------------------------------------------------! USE arrays_3d USE cloud_parameters USE constants USE particle_attributes IMPLICIT NONE PRIVATE PUBLIC colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi INTEGER, SAVE :: ip, jp, kp, pend, pstart, psum REAL, SAVE :: epsilon, eps2, urms, urms2 REAL, DIMENSION(:), ALLOCATABLE, SAVE :: winf REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: ec, ecf, gck ! !-- Public interfaces INTERFACE colker MODULE PROCEDURE colker END INTERFACE colker INTERFACE effic MODULE PROCEDURE effic END INTERFACE effic INTERFACE fallg MODULE PROCEDURE fallg END INTERFACE fallg INTERFACE phi MODULE PROCEDURE phi END INTERFACE phi INTERFACE turbsd MODULE PROCEDURE turbsd END INTERFACE turbsd INTERFACE turb_enhance_eff MODULE PROCEDURE turb_enhance_eff END INTERFACE turb_enhance_eff INTERFACE zhi MODULE PROCEDURE zhi END INTERFACE zhi CONTAINS !------------------------------------------------------------------------------! ! SUBROUTINE for calculation of collision kernel !------------------------------------------------------------------------------! SUBROUTINE colker( i1, j1, k1, kernel ) USE arrays_3d USE cloud_parameters USE constants USE cpulog USE indices USE interfaces USE particle_attributes IMPLICIT NONE INTEGER :: i, i1, j, j1, k1 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 ! CALL cpu_log( log_point_s(46), 'colker', 'start' ) ip = i1 jp = j1 kp = k1 pstart = prt_start_index(kp,jp,ip) pend = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1 psum = prt_count(kp,jp,ip) ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) ) IF ( turbulence_effects_on_collision ) THEN ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) ) epsilon = diss(kp,jp,ip) * 1.0E5 !dissipation rate in cm**2/s**-3 urms = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 ) IF ( epsilon <= 0.001 ) THEN CALL fallg CALL effic DO j = pstart, pend DO i = pstart, pend kernel(i,j) = pi * ( particles(j)%radius * 100.0 + & particles(i)%radius * 100.0 )**2 & * ec(i,j) * ABS( winf(j) - winf(i) ) ENDDO ENDDO ELSE CALL turbsd CALL turb_enhance_eff CALL effic DO j = pstart, pend, 1 DO i = pstart, pend, 1 kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j) ENDDO ENDDO ENDIF DEALLOCATE(gck, ecf) ELSE ! CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' ) CALL fallg ! CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' ) ! CALL cpu_log( log_point_s(51), 'colker_effic', 'start' ) CALL effic ! CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' ) DO j = pstart, pend DO i = pstart, pend kernel(i,j) = pi * ( particles(j)%radius * 100.0 + & particles(i)%radius * 100.0 )**2 & * ec(i,j) * ABS( winf(j) - winf(i) ) ENDDO ENDDO ENDIF DEALLOCATE( ec, winf ) ! CALL cpu_log( log_point_s(46), 'colker', 'stop' ) END SUBROUTINE colker !------------------------------------------------------------------------------! ! SUBROUTINE for calculation of w, g and gck !------------------------------------------------------------------------------! SUBROUTINE turbsd ! 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 USE constants USE cloud_parameters USE particle_attributes USE arrays_3d IMPLICIT NONE REAL :: Relamda, & Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be, & bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq, & vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, & wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens REAL, DIMENSION(pstart:pend) :: St, tau LOGICAL, SAVE :: first = .TRUE. INTEGER :: i, j ! !-- Initial assignment of constants IF ( first ) THEN first = .FALSE. airvisc = 0.1818 !dynamic viscosity in mg/cm*s airdens = 1.2250 !air density in mg/cm**3 waterdens = 1000.0 !water density in mg/cm**3 gravity = 980.6650 !in cm/s**2 anu = airvisc/airdens ! kinetic viscosity in cm**2/s ENDIF Relamda = urms**2*sqrt(15.0/epsilon/anu) Tl = urms**2/epsilon !in s Lf = 0.5 * (urms**3)/epsilon !in cm tauk = (anu/epsilon)**0.5 !in s eta = (anu**3/epsilon)**0.25 !in cm vk = eta/tauk !in cm/s ao = (11.+7.*Relamda)/(205.+Relamda) lambda = urms * sqrt(15.*anu/epsilon) !in cm tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s CALL fallg !gives winf in cm/s DO i = pstart, pend tau(i) = winf(i)/gravity !in s St(i) = tau(i)/tauk ENDDO !***** TO CALCULATE wr ******************************** !from Aayala 2008b, page 38f z = tt/Tl be = sqrt(2.0)*lambda/Lf bbb = sqrt(1.0-2.0*be**2) d1 = (1.+bbb)/2.0/bbb e1 = Lf*(1.0+bbb)/2.0 !in cm d2 = (1.0-bbb)/2.0/bbb e2 = Lf*(1.0-bbb)/2.0 !in cm ccc = sqrt(1.0-2.0*z**2) b1 = (1.+ccc)/2./ccc c1 = Tl*(1.+ccc)/2. !in s b2 = (1.-ccc)/2./ccc c2 = Tl*(1.-ccc)/2. !in s DO i = pstart, pend v1 = winf(i) !in cm/s t1 = tau(i) !in s DO j = pstart,i rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm v2 = winf(j) !in cm/s t2 = tau(j) !in s v1xysq = b1*d1*PHI(c1,e1,v1,t1) & - b1*d2*PHI(c1,e2,v1,t1) & - b2*d1*PHI(c2,e1,v1,t1) & + b2*d2*PHI(c2,e2,v1,t1) v1xysq = v1xysq * urms**2/t1 !in cm**2/s**2 vrms1xy= sqrt(v1xysq) !in cm/s v2xysq = b1*d1*PHI(c1,e1,v2,t2) & - b1*d2*PHI(c1,e2,v2,t2) & - b2*d1*PHI(c2,e1,v2,t2) & + b2*d2*PHI(c2,e2,v2,t2) v2xysq = v2xysq * urms**2/t2 !in cm**2/s**2 vrms2xy= sqrt(v2xysq) !in cm/s IF(winf(i).ge.winf(j)) THEN v1 = winf(i) t1 = tau(i) v2 = winf(j) t2 = tau(j) ELSE v1 = winf(j) t1 = tau(j) v2 = winf(i) t2 = tau(i) ENDIF v1v2xy = b1*d1*ZHI(c1,e1,v1,t1,v2,t2) & - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) & - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) & + b2*d2*ZHI(c2,e2,v1,t1,v2,t2) fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2) v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j) !in cm**2/s**2 wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy !in cm**2/s**2 IF (wrtur2xy.le.0.0) wrtur2xy=0.0 wrgrav2=pi/8.*(winf(j)-winf(i))**2 wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2)) !in cm/s !***** TO CALCULATE gr ******************************** IF(St(j).gt.St(i)) THEN SSt = St(j) ELSE SSt = St(i) ENDIF XX = -0.1988*SSt**4 + 1.5275*SSt**3 & -4.2942*SSt**2 + 5.3406*SSt IF(XX.le.0.0) XX = 0.0 YY = 0.1886*exp(20.306/Relamda) c1_gr = XX/(gravity/(vk/tauk))**YY ao_gr = ao + (pi/8.)*(gravity/(vk/tauk))**2 fao_gr = 20.115 * (ao_gr/Relamda)**0.5 rc = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta !in cm grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.) IF (grFIN.lt.1.0) grFIN = 1.0 gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN ! in cm**3/s gck(j,i) = gck(i,j) ENDDO ENDDO END SUBROUTINE TurbSD !------------------------------------------------------------------------------! ! PHI as a function !------------------------------------------------------------------------------! REAL FUNCTION PHI(a,b,vsett,tau0) IMPLICIT NONE REAL :: a, aa1, b, vsett, tau0 aa1 = 1./tau0 + 1./a + vsett/b PHI = 1./aa1 - vsett/2.0/b/aa1**2 !in s END FUNCTION PHI !------------------------------------------------------------------------------! ! ZETA as a function !------------------------------------------------------------------------------! REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2) IMPLICIT NONE REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2 aa1 = vsett2/b - 1./tau2 - 1./a aa2 = vsett1/b + 1./tau1 + 1./a aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2 aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2 aa5 = vsett2/b + 1./tau2 + 1./a aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2 ZHI = (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2 & + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6 & + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2) & * 1./2./b/aa3 ! in s**2 END FUNCTION ZHI !------------------------------------------------------------------------------! ! SUBROUTINE for calculation of terminal velocity winf !------------------------------------------------------------------------------! SUBROUTINE fallg USE constants USE cloud_parameters USE particle_attributes USE arrays_3d IMPLICIT NONE INTEGER :: i, j LOGICAL, SAVE :: first = .TRUE. REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py REAL :: bond, x, xrey, y REAL, DIMENSION(1:7), SAVE :: b REAL, DIMENSION(1:6), SAVE :: c ! !-- Initial assignment of constants IF ( first ) THEN first = .FALSE. b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, & -0.578878e-3,0.855176e-4,-0.327815e-5/) c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0, & -0.542819e-1, 0.238449e-2/) eta = 1.818e-4 !in poise = g/(cm s) xlamb = 6.62e-6 !in cm rhow = 1.0 !in g/cm**3 rhoa = 1.225e-3 !in g/cm**3 grav = 980.665 !in cm/s**2 cunh = 1.257 * xlamb !in cm t0 = 273.15 !in K sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2 stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s) stb = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3 phy = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa)) py = phy**(1.0/6.0) ENDIF !particle radius has to be in cm DO j = pstart, pend IF (particles(j)%radius*100.0 .le. 1.e-3) THEN winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN x = log(stb*(particles(j)%radius*100.0)**3) y = 0.0 DO i = 1, 7 y = y + b(i) * (x**(i-1)) ENDDO xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y) winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0) !in cm/s ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN IF (particles(j)%radius*100.0.gt.0.35) THEN bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma ELSE bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma ENDIF x = log(16.0*bond*py/3.0) y = 0.0 DO i = 1, 6 y = y + c(i) * (x**(i-1)) ENDDO xrey = py*exp(y) IF (particles(j)%radius*100.0 .gt.0.35) THEN winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s ELSE winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0) !in cm/s ENDIF ENDIF ENDDO RETURN END SUBROUTINE fallg !------------------------------------------------------------------------------! ! SUBROUTINE for calculation of collision efficencies !------------------------------------------------------------------------------! SUBROUTINE effic USE arrays_3d USE constants USE cloud_parameters USE particle_attributes !collision efficiencies of hall kernel IMPLICIT NONE INTEGER :: i, ir, iq, j, k, kk INTEGER, DIMENSION(:), ALLOCATABLE :: ira LOGICAL, SAVE :: first = .TRUE. REAL :: ek, particle_radius, pp, qq, rq REAL, DIMENSION(1:21), SAVE :: rat REAL, DIMENSION(1:15), SAVE :: r0 REAL, DIMENSION(1:15,1:21), SAVE :: ecoll ! !-- Initial assignment of constants IF ( first ) THEN first = .FALSE. r0 = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., & 300./) 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, & 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) 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/) ENDIF ! !-- Calculate the radius class index of particles with respect to array r ALLOCATE( ira(pstart:pend) ) DO j = pstart, pend particle_radius = particles(j)%radius * 1.0E6 DO k = 1, 15 IF ( particle_radius < r0(k) ) THEN ira(j) = k EXIT ENDIF ENDDO IF ( particle_radius >= r0(15) ) ira(j) = 16 ENDDO ! !-- Two-dimensional linear interpolation of the collision efficiency. !-- Radius has to be in µm DO j = pstart, pend DO i = pstart, j ir = ira(j) rq = particles(i)%radius / particles(j)%radius ! DO kk = 2, 21 ! IF ( rq <= rat(kk) ) THEN ! iq = kk ! EXIT ! ENDIF ! ENDDO iq = INT( rq * 20 ) + 1 iq = MAX(iq , 2) IF ( ir < 16 ) THEN IF ( ir >= 2 ) THEN pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / & ( r0(ir) - r0(ir-1) ) qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1) & + pp * ( 1.0-qq ) * ecoll(ir,iq-1) & + qq * ( 1.0-pp ) * ecoll(ir-1,iq) & + pp * qq * ecoll(ir,iq) ELSE qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq) ENDIF ELSE qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) ) ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq) ec(j,i) = MIN( ek, 1.0 ) ENDIF ec(i,j) = ec(j,i) IF ( ec(i,j) < 1.0E-20 ) ec(i,j) = 0.0 ENDDO ENDDO DEALLOCATE( ira ) END SUBROUTINE effic !------------------------------------------------------------------------------! ! SUBROUTINE for calculation of enhancement factor collision efficencies !------------------------------------------------------------------------------! SUBROUTINE turb_enhance_eff USE constants USE cloud_parameters USE particle_attributes USE arrays_3d IMPLICIT NONE INTEGER :: i, ik, ir, iq, j, k, kk INTEGER, DIMENSION(:), ALLOCATABLE :: ira REAL :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3 LOGICAL, SAVE :: first = .TRUE. REAL, DIMENSION(1:11), SAVE :: rat REAL, DIMENSION(1:7), SAVE :: r0 REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400 ! !-- Initial assignment of constants IF ( first ) THEN first = .FALSE. r0 = (/10., 20., 30.,40., 50., 60.,100./) rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/) ! 100 cm^2/s^3 ecoll_100(:,1) = (/1.74, 1.74, 1.773, 1.49, 1.207, 1.207, 1.0 /) ecoll_100(:,2) = (/1.46, 1.46, 1.421, 1.245, 1.069, 1.069, 1.0 /) ecoll_100(:,3) = (/1.32, 1.32, 1.245, 1.123, 1.000, 1.000, 1.0 /) ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025, 1.025, 1.0 /) ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056, 1.056, 1.0 /) ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028, 1.028, 1.0 /) ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046, 1.046, 1.0 /) ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029, 1.029, 1.0 /) ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021, 1.021, 1.0 /) ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088, 1.088, 1.0 /) ecoll_100(:,11)= (/20.3, 20.3, 14.6 , 8.61, 2.60, 2.60 , 1.0 /) ! 400 cm^2/s^3 ecoll_400(:,1) = (/4.976, 4.976, 3.593, 2.519, 1.445, 1.445, 1.0 /) ecoll_400(:,2) = (/2.984, 2.984, 2.181, 1.691, 1.201, 1.201, 1.0 /) ecoll_400(:,3) = (/1.988, 1.988, 1.475, 1.313, 1.150, 1.150, 1.0 /) ecoll_400(:,4) = (/1.490, 1.490, 1.187, 1.156, 1.126, 1.126, 1.0 /) ecoll_400(:,5) = (/1.249, 1.249, 1.088, 1.090, 1.092, 1.092, 1.0 /) ecoll_400(:,6) = (/1.139, 1.139, 1.130, 1.091, 1.051, 1.051, 1.0 /) ecoll_400(:,7) = (/1.220, 1.220, 1.190, 1.138, 1.086, 1.086, 1.0 /) ecoll_400(:,8) = (/1.325, 1.325, 1.267, 1.165, 1.063, 1.063, 1.0 /) ecoll_400(:,9) = (/1.716, 1.716, 1.345, 1.223, 1.100, 1.100, 1.0 /) ecoll_400(:,10)= (/3.788, 3.788, 1.501, 1.311, 1.120, 1.120, 1.0 /) ecoll_400(:,11)= (/36.52, 36.52, 19.16, 22.80, 26.0, 26.0, 1.0 /) ENDIF ! !-- Calculate the radius class index of particles with respect to array r ALLOCATE( ira(pstart:pend) ) DO j = pstart, pend particle_radius = particles(j)%radius * 1.0E6 DO k = 1, 7 IF ( particle_radius < r0(k) ) THEN ira(j) = k EXIT ENDIF ENDDO IF ( particle_radius >= r0(7) ) ira(j) = 8 ENDDO ! two-dimensional linear interpolation of the collision efficiency DO j = pstart, pend DO i = pstart, j ir = ira(j) rq = particles(i)%radius/particles(j)%radius DO kk = 2, 11 IF ( rq <= rat(kk) ) THEN iq = kk EXIT ENDIF ENDDO ! 0 cm2/s3 y1 = 1.0 ! 100 cm2/s3, 400 cm2/s3 IF (ir.lt.8) THEN IF (ir.ge.2) THEN pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1)) qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+ & pp*(1.-qq)*ecoll_100(ir,iq-1)+ & qq*(1.-pp)*ecoll_100(ir-1,iq)+ & pp*qq*ecoll_100(ir,iq) y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+ & pp*(1.-qq)*ecoll_400(ir,iq-1)+ & qq*(1.-pp)*ecoll_400(ir-1,iq)+ & pp*qq*ecoll_400(ir,iq) ELSE qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq) y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq) ENDIF ELSE qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1)) y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq) y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq) ENDIF ! linear interpolation ! dissipation rate in cm2/s3 x1 = 0.0 x2 = 100.0 x3 = 400.0 IF (epsilon.le.100.) THEN ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 & + (epsilon-0.)/(100.-0.) * y2 ELSE IF(epsilon.le.600.)THEN ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 & + (epsilon-100.)/(400.-100.) * y3 ELSE ecf(j,i) = (600.-400.)/(100.-400.) * y2 & + (600.-100.)/(400.-100.) * y3 ENDIF IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0 ecf(i,j)=ecf(j,i) ENDDO ENDDO RETURN END SUBROUTINE turb_enhance_eff END MODULE wang_kernel_mod