MODULE wang_kernel_mod !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: wang_kernel.f90 791 2011-11-29 03:33:42Z raasch $ ! ! 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 TurbSD, turb_enhance_eff, effic, fallg, PHI, ZHI, colker INTEGER, SAVE :: ip, jp, kp, pstart, pend, psum REAL, SAVE :: epsilon, urms REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: gck, ec, ecf REAL, DIMENSION(:), ALLOCATABLE, SAVE :: winf REAL, SAVE :: eps2, urms2 ! !-- Public interfaces INTERFACE TurbSD MODULE PROCEDURE TurbSD END INTERFACE TurbSD INTERFACE PHI MODULE PROCEDURE PHI END INTERFACE PHI INTERFACE ZHI MODULE PROCEDURE ZHI END INTERFACE ZHI INTERFACE turb_enhance_eff MODULE PROCEDURE turb_enhance_eff END INTERFACE turb_enhance_eff INTERFACE effic MODULE PROCEDURE effic END INTERFACE effic INTERFACE fallg MODULE PROCEDURE fallg END INTERFACE fallg INTERFACE colker MODULE PROCEDURE colker END INTERFACE colker CONTAINS !------------------------------------------------------------------------------! ! SUBROUTINE for calculation of collision kernel !------------------------------------------------------------------------------! SUBROUTINE colker(i1,j1,k1,kernel) USE arrays_3d USE cloud_parameters USE constants USE indices 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 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(winf(pstart:pend), ec(pstart:pend,pstart:pend)) IF ( turbulence_effects_on_collision ) THEN ALLOCATE(gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend)) epsilon = diss(kp,jp,ip) * 10000.0 !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, 1 DO i = pstart, pend, 1 kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * 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 fallg CALL effic DO j = pstart, pend, 1 DO i = pstart, pend, 1 kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * ec(i,j) * abs(winf(j)-winf(i)) ENDDO ENDDO ENDIF DEALLOCATE(winf, ec) RETURN 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 :: airvisc, airdens, waterdens, gravity, anu, 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, DIMENSION(pstart:pend) :: vST, tau, St INTEGER :: i, j 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 Relamda = urms**2.0*sqrt(15.0/epsilon/anu) Tl = urms**2.0/epsilon !in s Lf = 0.5 * (urms**3.0)/epsilon !in cm tauk = (anu/epsilon)**0.5 !in s eta = (anu**3.0/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 vST(i) = winf(i) !in cm/s tau(i) = vST(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.0) 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.0) 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 = vST(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 = vST(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.0/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.0/t2 !in cm**2/s**2 vrms2xy= sqrt(v2xysq) !in cm/s IF(vST(i).ge.vST(j)) THEN v1 = vST(i) t1 = tau(i) v2 = vST(j) t2 = tau(j) ELSE v1 = vST(j) t1 = tau(j) v2 = vST(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.0/tau(i)/tau(j) !in cm**2/s**2 wrtur2xy=vrms1xy**2.0 + vrms2xy**2.0 - 2.*v1v2xy !in cm**2/s**2 IF (wrtur2xy.le.0.0) wrtur2xy=0.0 wrgrav2=pi/8.*(vST(j)-vST(i))**2.0 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.0 + 1.5275*SSt**3.0 & -4.2942*SSt**2.0 + 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.0 fao_gr = 20.115 * (ao_gr/Relamda)**0.5 rc = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta !in cm grFIN = ((eta**2.0+rc**2.0)/(rrp**2.0+rc**2.0))**(c1_gr/2.) IF (grFIN.lt.1.0) grFIN = 1.0 gck(i,j) = 2. * pi * rrp**2.0 * wrFIN * grFIN ! in cm**3/s gck(j,i) = gck(i,j) ENDDO ENDDO RETURN 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.0 !in s RETURN 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.0 - (1./tau2 + 1./a)**2.0 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.0 & + (4./aa4 - 1./aa5**2.0 - 1./aa1**2.0) * vsett2/2./b/aa6 & + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2.0 + vsett2/aa1**2.0) & * 1./2./b/aa3 ! in s**2 RETURN 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 REAL :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py, & x, y, xrey, bond REAL, DIMENSION(1:7) :: b REAL, DIMENSION(1:6) :: c REAL, DIMENSION(1:20) :: rat REAL, DIMENSION(1:15) :: r0 REAL, DIMENSION(1:15,1:20) :: ecoll 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.0) * (rhoa**2.0)/((eta**4.0) * grav * (rhow-rhoa)) py = phy**(1.0/6.0) !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.0) 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 bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2.0/sigma IF (particles(j)%radius*100.0.gt.0.35) THEN bond = grav*(rhow-rhoa) * 0.35 * 0.35/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) winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0) !in cm/s IF (particles(j)%radius*100.0 .gt.0.35) THEN winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s ENDIF ENDIF ENDDO RETURN END SUBROUTINE fallg !------------------------------------------------------------------------------! ! SUBROUTINE for calculation of collision efficencies !------------------------------------------------------------------------------! SUBROUTINE effic USE constants USE cloud_parameters USE particle_attributes USE arrays_3d !collision efficiencies of hall kernel IMPLICIT NONE INTEGER :: i, ir, iq, j, k, kk REAL :: ek, rq, pp, qq REAL, DIMENSION(1:21) :: rat REAL, DIMENSION(1:15) :: r0 REAL, DIMENSION(1:15,1:21) :: ecoll 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/) ! two-dimensional linear interpolation of the collision efficiency ! radius has to be in µm DO j = pstart, pend DO i = pstart, j DO k = 2, 15 IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN ir = k ELSEIF ((particles(j)%radius*1.0E06).gt.r0(15)) THEN ir = 16 ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN ir = 1 ENDIF ENDDO rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06) DO kk = 2, 21 IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk ENDDO IF (ir.lt.16) 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)) ec(j,i) = (1.-pp) * (1.-qq) * ecoll(ir-1,iq-1) + & pp * (1.-qq) * ecoll(ir,iq-1) + & qq * (1.-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.-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).lt.1.e-20) stop 99 ENDDO ENDDO RETURN 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 !collision efficiencies of hall kernel IMPLICIT NONE INTEGER :: i, ik, ir, iq, j, k, kk REAL :: rq, y1, pp, qq, y2, y3, x1, x2, x3 REAL, DIMENSION(1:11) :: rat REAL, DIMENSION(1:7) :: r0 REAL, DIMENSION(1:7,1:11) :: ecoll_100, ecoll_400 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 /) ! two-dimensional linear interpolation of the collision efficiency DO j = pstart, pend DO i = pstart, j DO k = 2, 7 IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN ir = k ELSEIF ((particles(j)%radius*1.0E06).gt.r0(7)) THEN ir = 8 ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN ir = 1 ENDIF ENDDO rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06) DO kk = 2, 11 IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk 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