subroutine da_cld_eff_radius(t,rho,qci,qrn,qsn,qgr,snow,xice,xland,method, & reff_water,reff_ice,reff_rain,reff_snow,reff_grau) !--------------------------------------------------------------------------- ! Purpose: compute effective radius of cloud liquid water and cloud ice ! ! METHOD: ! liquid water: adapted from WRFV2/phys/module_ra_cam.F, analytic formula following ! the formulation originally developed by J. T. Kiehl, and ! ice (method 1): adapted from WRFV2/phys/module_ra_cam.F, Kristjansson and Mitchell ! ice (method 2): WSM3, Hong et al., MWR 2004 ! rain/snow/graupel: assume exponential particle size distribution and ! spherical particles ! use Gauss-Laguerre Quadrature for integration ! HISTORY: 12/15/2008 effective radius unit is micron. Zhiquan Liu !--------------------------------------------------------------------------- integer, intent(in) :: method real*8, intent(in) :: t ! temperature real, intent(in) :: rho ! air density kg/m3 real*8, intent(in) :: qci ! cloud ice mixing ratio kg/kg real*8, intent(in) :: qrn ! cloud rain mixing ratio real*8, intent(in) :: qsn ! cloud snow mixing ratio real*8, intent(in) :: qgr ! cloud graupel mixing ratio real, intent(in) :: snow ! snow water equivalent real*8, intent(in) :: xice ! ice percentage real*8, intent(in) :: xland ! landsea percentage real*8, intent(out) :: reff_water ! effective radius liquid water real*8, intent(out) :: reff_ice ! effective radius ice real*8, intent(out) :: reff_rain ! effective radius rain real*8, intent(out) :: reff_snow ! effective radius snow real*8, intent(out) :: reff_grau ! effective radius graupel ! local variables integer :: index, nk real :: snowh, corr real, parameter :: rliqland = 8.0 ! liquid drop size if over land real, parameter :: rliqocean = 14.0 ! liquid drop size if over ocean real, parameter :: rliqice = 14.0 ! liquid drop size if over sea ice ! Tabulated values of re(T) in the temperature interval ! 180 K -- 274 K; hexagonal columns assumed: real, dimension(95), parameter :: retab = & (/ 5.92779, 6.26422, 6.61973, 6.99539, 7.39234, & 7.81177, 8.25496, 8.72323, 9.21800, 9.74075, 10.2930, & 10.8765, 11.4929, 12.1440, 12.8317, 13.5581, 14.2319, & 15.0351, 15.8799, 16.7674, 17.6986, 18.6744, 19.6955, & 20.7623, 21.8757, 23.0364, 24.2452, 25.5034, 26.8125, & 27.7895, 28.6450, 29.4167, 30.1088, 30.7306, 31.2943, & 31.8151, 32.3077, 32.7870, 33.2657, 33.7540, 34.2601, & 34.7892, 35.3442, 35.9255, 36.5316, 37.1602, 37.8078, & 38.4720, 39.1508, 39.8442, 40.5552, 41.2912, 42.0635, & 42.8876, 43.7863, 44.7853, 45.9170, 47.2165, 48.7221, & 50.4710, 52.4980, 54.8315, 57.4898, 60.4785, 63.7898, & 65.5604, 71.2885, 75.4113, 79.7368, 84.2351, 88.8833, & 93.6658, 98.5739, 103.603, 108.752, 114.025, 119.424, & 124.954, 130.630, 136.457, 142.446, 148.608, 154.956, & 161.503, 168.262, 175.248, 182.473, 189.952, 197.699, & 205.728, 214.055, 222.694, 231.661, 240.971, 250.639 /) ! ! constants from da_control.f: pi, ! real, parameter :: n0_rain = 0.08 ! cm(-4) real, parameter :: n0_snow = 0.04 ! cm(-4) real, parameter :: n0_grau = 0.04 ! cm(-4) real, parameter :: rho_w = 1000.0 ! kg m(-3) real, parameter :: rho_ice = 900.0 ! kg m(-3) real, parameter :: rho_snow = 100.0 ! kg m(-3) real, parameter :: rho_grau = 400.0 ! kg m(-3) ! Abscissas of Gauss-Laguerre Integration real, dimension(32) :: xk = (/ 0.0444893658333, 0.23452610952, & 0.576884629302, 1.07244875382, 1.72240877644, 2.52833670643, & 3.49221327285, 4.61645677223, 5.90395848335, 7.3581268086, & 8.98294126732, 10.783012089, 12.763745476, 14.9309117981, & 17.2932661372, 19.8536236493, 22.6357789624, 25.6201482024, & 28.8739336869, 32.3333294017, 36.1132042245, 40.1337377056, & 44.5224085362, 49.2086605665, 54.3501813324, 59.8791192845, & 65.9833617041, 72.6842683222, 80.1883747906, 88.735192639, & 98.8295523184, 111.751398227 /) ! total weights (weight*exp(xk)) of Gauss-Laguerre Integration real, dimension(32) :: totalw = (/ 0.114187105768, 0.266065216898, & 0.418793137325, 0.572532846497, 0.727648788453, 0.884536718946, & 1.04361887597, 1.20534920595, 1.37022171969, 1.53877595906, & 1.71164594592, 1.8895649683, 2.07318851235, 2.26590144444, & 2.46997418988, 2.64296709494, 2.76464437462, 3.22890542981, & 2.92019361963, 4.3928479809, 4.27908673189, 5.20480398519, & 5.11436212961, 4.15561492173, 6.19851060567, 5.34795780128, & 6.28339212457, 6.89198340969, 7.92091094244, 9.20440555803, & 11.1637432904, 15.3902417688 /) real, parameter :: limit = 1.0e-6 real :: piover6 ! pi/6 real :: sum1_rain, sum2_rain, sum1_snow, sum2_snow, sum1_grau, sum2_grau real :: lamda_rain, lamda_snow, lamda_grau real, dimension(32) :: psd_rain, psd_snow, psd_grau ! partical size distribution ! initialize reff_water = 0.0 reff_ice = 0.0 reff_rain = 0.0 reff_snow = 0.0 reff_grau = 0.0 lamda_rain = 0.0 lamda_snow = 0.0 lamda_grau = 0.0 ! cloud liquid effective radius snowh = 0.001 * snow ! here the snowh (in meter) is water-equivalent snow depth, ! which is different from the actually snow depth defined in the wrf output file ! Start with temperature-dependent value appropriate for continental air reff_water = rliqland + (rliqocean-rliqland) * min(1.0_8,max(0.0_8,(t_triple-t)*0.05_8)) ! Modify for snow depth over land reff_water = reff_water + (rliqocean-reff_water) * min(1.0,max(0.0,snowh*10.)) ! Ramp between polluted value over land to clean value over ocean. reff_water = reff_water + (rliqocean-reff_water) * min(1.0_8,max(0.0_8,1.0_8-xland)) ! Ramp between the resultant value and a sea ice value in the presence of ice. reff_water = reff_water + (rliqice-reff_water) * min(1.0_8,max(0.0_8,xice)) ! cloud ice effective radius if ( method == 1 ) then index = int(t-179.) index = min(max(index,1),94) corr = t - int(t) reff_ice = retab(index)*(1.-corr) + retab(index+1)*corr end if ! cloud ice effective radius ! rho*qci = 2.08*10**22 * Dice**8 if ( method == 2 ) then ! 0.5 for diameter - radii conversion ! 1.0e6 for meter - micron conversion ! 0.125 = 1/8 reff_ice = 1.0e6 * 0.5 * ( rho * qci / 2.08e22 ) ** 0.125 end if ! ! cloud rain/snow/graupel effective radius ! piover6 = pi/6. if ( qrn > limit ) then lamda_rain = (piover6*rho_w*n0_rain*rho/qrn)**0.25 end if if ( qsn > limit ) then lamda_snow = (piover6*rho_snow*n0_snow*rho/qsn)**0.25 end if if ( qgr > limit ) then lamda_grau = (piover6*rho_grau*n0_grau*rho/qgr)**0.25 end if sum1_rain = 0.0 sum2_rain = 0.0 sum1_snow = 0.0 sum2_snow = 0.0 sum1_grau = 0.0 sum2_grau = 0.0 if ( qrn > limit ) then do nk = 1, 32 psd_rain(nk) = n0_rain*exp(-2.0*lamda_rain*xk(nk)) sum1_rain = sum1_rain + totalw(nk) * (xk(nk)**3) * psd_rain(nk) sum2_rain = sum2_rain + totalw(nk) * (xk(nk)**2) * psd_rain(nk) end do reff_rain = 10000.0 * sum1_rain/sum2_rain ! micron end if if ( qsn > limit ) then do nk = 1, 32 psd_snow(nk) = n0_snow*exp(-2.0*lamda_snow*xk(nk)) sum1_snow = sum1_snow + totalw(nk) * (xk(nk)**3) * psd_snow(nk) sum2_snow = sum2_snow + totalw(nk) * (xk(nk)**2) * psd_snow(nk) end do reff_snow = 10000.0 * sum1_snow/sum2_snow ! micron end if if ( qgr > limit ) then do nk = 1, 32 psd_grau(nk) = n0_grau*exp(-2.0*lamda_grau*xk(nk)) sum1_grau = sum1_grau + totalw(nk) * (xk(nk)**3) * psd_grau(nk) sum2_grau = sum2_grau + totalw(nk) * (xk(nk)**2) * psd_grau(nk) end do reff_grau = 10000.0 * sum1_grau/sum2_grau ! micron end if end subroutine da_cld_eff_radius