[2428] | 1 | ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ |
---|
| 2 | ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/quickbeam/zeff.f90 $ |
---|
[1262] | 3 | subroutine zeff(freq,D,N,nsizes,k2,tt,ice,xr,z_eff,z_ray,kr,qe,qs,rho_e) |
---|
| 4 | use math_lib |
---|
| 5 | use optics_lib |
---|
| 6 | implicit none |
---|
| 7 | |
---|
| 8 | ! Purpose: |
---|
| 9 | ! Simulates radar return of a volume given DSD of spheres |
---|
| 10 | ! Part of QuickBeam v1.03 by John Haynes |
---|
| 11 | ! http://reef.atmos.colostate.edu/haynes/radarsim |
---|
| 12 | ! |
---|
| 13 | ! Inputs: |
---|
| 14 | ! [freq] radar frequency (GHz) |
---|
| 15 | ! [D] discrete drop sizes (um) |
---|
| 16 | ! [N] discrete concentrations (cm^-3 um^-1) |
---|
| 17 | ! [nsizes] number of discrete drop sizes |
---|
| 18 | ! [k2] |K|^2, -1=use frequency dependent default |
---|
[2428] | 19 | ! [tt] hydrometeor temperature (K) |
---|
[1262] | 20 | ! [ice] indicates volume consists of ice |
---|
| 21 | ! [xr] perform Rayleigh calculations? |
---|
| 22 | ! [qe] if using a mie table, these contain ext/sca ... |
---|
| 23 | ! [qs] ... efficiencies; otherwise set to -1 |
---|
| 24 | ! [rho_e] medium effective density (kg m^-3) (-1 = pure) |
---|
| 25 | ! |
---|
| 26 | ! Outputs: |
---|
| 27 | ! [z_eff] unattenuated effective reflectivity factor (mm^6/m^3) |
---|
| 28 | ! [z_ray] reflectivity factor, Rayleigh only (mm^6/m^3) |
---|
| 29 | ! [kr] attenuation coefficient (db km^-1) |
---|
| 30 | ! |
---|
| 31 | ! Created: |
---|
| 32 | ! 11/28/05 John Haynes (haynes@atmos.colostate.edu) |
---|
| 33 | |
---|
| 34 | ! ----- INPUTS ----- |
---|
| 35 | integer, intent(in) :: ice, xr |
---|
| 36 | integer, intent(in) :: nsizes |
---|
| 37 | real*8, intent(in) :: freq,D(nsizes),N(nsizes),tt,qe(nsizes), & |
---|
| 38 | qs(nsizes), rho_e(nsizes) |
---|
| 39 | real*8, intent(inout) :: k2 |
---|
| 40 | |
---|
| 41 | ! ----- OUTPUTS ----- |
---|
| 42 | real*8, intent(out) :: z_eff,z_ray,kr |
---|
| 43 | |
---|
| 44 | ! ----- INTERNAL ----- |
---|
| 45 | integer :: & |
---|
[2428] | 46 | correct_for_rho ! correct for density flag |
---|
[1262] | 47 | real*8, dimension(nsizes) :: & |
---|
[2428] | 48 | D0, & ! D in (m) |
---|
| 49 | N0, & ! N in m^-3 m^-1 |
---|
| 50 | sizep, & ! size parameter |
---|
| 51 | qext, & ! extinction efficiency |
---|
| 52 | qbsca, & ! backscatter efficiency |
---|
| 53 | rho_ice, & ! bulk density ice (kg m^-3) |
---|
| 54 | f ! ice fraction |
---|
| 55 | real*8, dimension(nsizes) :: xtemp |
---|
[1262] | 56 | real*8 :: & |
---|
[2428] | 57 | wl, & ! wavelength (m) |
---|
[1262] | 58 | cr ! kr(dB/km) = cr * kr(1/km) |
---|
| 59 | complex*16 :: & |
---|
[2428] | 60 | m ! complex index of refraction of bulk form |
---|
[1262] | 61 | complex*16, dimension(nsizes) :: & |
---|
[2428] | 62 | m0 ! complex index of refraction |
---|
[1262] | 63 | |
---|
| 64 | integer*4 :: i,one |
---|
| 65 | real*8 :: pi |
---|
| 66 | real*8 :: eta_sum, eta_mie, const, z0_eff, z0_ray, k_sum, & |
---|
[2428] | 67 | n_r, n_i, dqv(1), dqsc, dg, dph(1) |
---|
[1262] | 68 | integer*4 :: err |
---|
| 69 | complex*16 :: Xs1(1), Xs2(1) |
---|
| 70 | |
---|
| 71 | one=1 |
---|
| 72 | pi = acos(-1.0) |
---|
| 73 | rho_ice(:) = 917 |
---|
| 74 | z0_ray = 0.0 |
---|
| 75 | |
---|
| 76 | ! // conversions |
---|
[2428] | 77 | D0 = d*1E-6 ! m |
---|
| 78 | N0 = n*1E12 ! 1/(m^3 m) |
---|
| 79 | wl = 2.99792458/(freq*10) ! m |
---|
[1262] | 80 | |
---|
| 81 | ! // dielectric constant |k^2| defaults |
---|
| 82 | if (k2 < 0) then |
---|
| 83 | k2 = 0.933 |
---|
| 84 | if (abs(94.-freq) < 3.) k2=0.75 |
---|
| 85 | if (abs(35.-freq) < 3.) k2=0.88 |
---|
| 86 | if (abs(13.8-freq) < 3.) k2=0.925 |
---|
| 87 | endif |
---|
| 88 | |
---|
| 89 | if (qe(1) < -9) then |
---|
| 90 | |
---|
| 91 | ! // get the refractive index of the bulk hydrometeors |
---|
| 92 | if (ice == 0) then |
---|
| 93 | call m_wat(freq,tt,n_r,n_i) |
---|
| 94 | else |
---|
| 95 | call m_ice(freq,tt,n_r,n_i) |
---|
| 96 | endif |
---|
| 97 | m = cmplx(n_r,-n_i) |
---|
| 98 | m0(:) = m |
---|
| 99 | |
---|
| 100 | correct_for_rho = 0 |
---|
| 101 | if ((ice == 1) .and. (minval(rho_e) >= 0)) correct_for_rho = 1 |
---|
| 102 | |
---|
| 103 | ! :: correct refractive index for ice density if needed |
---|
| 104 | if (correct_for_rho == 1) then |
---|
| 105 | f = rho_e/rho_ice |
---|
| 106 | m0 = ((2+m0**2+2*f*(m0**2-1))/(2+m0**2+f*(1-m0**2)))**(0.5) |
---|
| 107 | endif |
---|
| 108 | |
---|
| 109 | ! :: Mie calculations |
---|
| 110 | sizep = (pi*D0)/wl |
---|
| 111 | dqv(1) = 0. |
---|
| 112 | do i=1,nsizes |
---|
| 113 | call mieint(sizep(i), m0(i), one, dqv, qext(i), dqsc, qbsca(i), & |
---|
| 114 | dg, xs1, xs2, dph, err) |
---|
| 115 | end do |
---|
| 116 | |
---|
| 117 | else |
---|
| 118 | ! // Mie table used |
---|
| 119 | |
---|
| 120 | qext = qe |
---|
| 121 | qbsca = qs |
---|
| 122 | |
---|
| 123 | endif |
---|
| 124 | |
---|
| 125 | ! // eta_mie = 0.25*sum[qbsca*pi*D^2*N(D)*deltaD] |
---|
| 126 | ! <--------- eta_sum ---------> |
---|
| 127 | ! // z0_eff = (wl^4/!pi^5)*(1./k2)*eta_mie |
---|
| 128 | eta_sum = 0. |
---|
| 129 | if (size(D0) == 1) then |
---|
| 130 | eta_sum = qbsca(1)*(n(1)*1E6)*D0(1)**2 |
---|
| 131 | else |
---|
[2428] | 132 | xtemp = qbsca*N0*D0**2 |
---|
| 133 | call avint(xtemp,D0,nsizes,D0(1),D0(size(D0,1)),eta_sum) |
---|
[1262] | 134 | endif |
---|
| 135 | |
---|
| 136 | eta_mie = eta_sum*0.25*pi |
---|
| 137 | const = (wl**4/pi**5)*(1./k2) |
---|
| 138 | z0_eff = const*eta_mie |
---|
| 139 | |
---|
| 140 | ! // kr = 0.25*cr*sum[qext*pi*D^2*N(D)*deltaD] |
---|
| 141 | ! <---------- k_sum ---------> |
---|
| 142 | k_sum = 0. |
---|
| 143 | if (size(D0) == 1) then |
---|
| 144 | k_sum = qext(1)*(n(1)*1E6)*D0(1)**2 |
---|
| 145 | else |
---|
[2428] | 146 | xtemp = qext*N0*D0**2 |
---|
| 147 | call avint(xtemp,D0,nsizes,D0(1),D0(size(D0,1)),k_sum) |
---|
[1262] | 148 | endif |
---|
| 149 | cr = 10./log(10.) |
---|
| 150 | kr = k_sum*0.25*pi*(1000.*cr) |
---|
[2428] | 151 | |
---|
[1262] | 152 | ! // z_ray = sum[D^6*N(D)*deltaD] |
---|
| 153 | if (xr == 1) then |
---|
| 154 | z0_ray = 0. |
---|
| 155 | if (size(D0) == 1) then |
---|
| 156 | z0_ray = (n(1)*1E6)*D0(1)**6 |
---|
| 157 | else |
---|
[2428] | 158 | xtemp = N0*D0**6 |
---|
| 159 | call avint(xtemp,D0,nsizes,D0(1),D0(size(D0)),z0_ray) |
---|
[1262] | 160 | endif |
---|
| 161 | endif |
---|
| 162 | |
---|
| 163 | ! // convert to mm^6/m^3 |
---|
| 164 | z_eff = z0_eff*1E18 ! 10.*alog10(z0_eff*1E18) |
---|
| 165 | z_ray = z0_ray*1E18 ! 10.*alog10(z0_ray*1E18) |
---|
| 166 | |
---|
| 167 | end subroutine zeff |
---|