!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! GNU General Public License !!
!! !!
!! This file is part of the Flexible Modeling System (FMS). !!
!! !!
!! FMS is free software; you can redistribute it and/or modify !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation. !!
!! !!
!! FMS is distributed in the hope that it will be useful, !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !!
!! GNU General Public License for more details. !!
!! !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to: !!
!! Free Software Foundation, Inc. !!
!! 59 Temple Place, Suite 330 !!
!! Boston, MA 02111-1307 USA !!
!! or see: !!
!! http://www.gnu.org/licenses/gpl.txt !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------
!------------ FMS version number and tagname for this file -----------
! $Id: zeff.f90,v 1.1.2.1.2.1 2009/08/10 10:48:14 rsh Exp $
! $Name: mom4p1_pubrel_dec2009_nnz $
subroutine zeff(freq,D,N,nsizes,k2,tt,ice,xr,z_eff,z_ray,kr,qe,qs,rho_e)
use math_lib
use optics_lib
implicit none
! Purpose:
! Simulates radar return of a volume given DSD of spheres
! Part of QuickBeam v1.03 by John Haynes
! http://reef.atmos.colostate.edu/haynes/radarsim
!
! Inputs:
! [freq] radar frequency (GHz)
! [D] discrete drop sizes (um)
! [N] discrete concentrations (cm^-3 um^-1)
! [nsizes] number of discrete drop sizes
! [k2] |K|^2, -1=use frequency dependent default
! [tt] hydrometeor temperature (C)
! [ice] indicates volume consists of ice
! [xr] perform Rayleigh calculations?
! [qe] if using a mie table, these contain ext/sca ...
! [qs] ... efficiencies; otherwise set to -1
! [rho_e] medium effective density (kg m^-3) (-1 = pure)
!
! Outputs:
! [z_eff] unattenuated effective reflectivity factor (mm^6/m^3)
! [z_ray] reflectivity factor, Rayleigh only (mm^6/m^3)
! [kr] attenuation coefficient (db km^-1)
!
! Created:
! 11/28/05 John Haynes (haynes@atmos.colostate.edu)
! ----- INPUTS -----
integer, intent(in) :: ice, xr
integer, intent(in) :: nsizes
real*8, intent(in) :: freq,D(nsizes),N(nsizes),tt,qe(nsizes), &
qs(nsizes), rho_e(nsizes)
real*8, intent(inout) :: k2
! ----- OUTPUTS -----
real*8, intent(out) :: z_eff,z_ray,kr
! ----- INTERNAL -----
integer :: &
correct_for_rho ! correct for density flag
real*8, dimension(nsizes) :: &
D0, & ! D in (m)
N0, & ! N in m^-3 m^-1
sizep, & ! size parameter
qext, & ! extinction efficiency
qbsca, & ! backscatter efficiency
rho_ice, & ! bulk density ice (kg m^-3)
f ! ice fraction
real*8 :: &
wl, & ! wavelength (m)
cr ! kr(dB/km) = cr * kr(1/km)
complex*16 :: &
m ! complex index of refraction of bulk form
complex*16, dimension(nsizes) :: &
m0 ! complex index of refraction
integer*4 :: i,one
real*8 :: pi
real*8 :: eta_sum, eta_mie, const, z0_eff, z0_ray, k_sum, &
n_r, n_i, dqv(1), dqxt, dqsc, dbsc, dg, dph(1)
integer*4 :: err
complex*16 :: Xs1(1), Xs2(1)
one=1
pi = acos(-1.0)
rho_ice(:) = 917
z0_ray = 0.0
! // conversions
D0 = d*1E-6 ! m
N0 = n*1E12 ! 1/(m^3 m)
wl = 2.99792458/(freq*10) ! m
! // dielectric constant |k^2| defaults
if (k2 < 0) then
k2 = 0.933
if (abs(94.-freq) < 3.) k2=0.75
if (abs(35.-freq) < 3.) k2=0.88
if (abs(13.8-freq) < 3.) k2=0.925
endif
if (qe(1) < -9) then
! // get the refractive index of the bulk hydrometeors
if (ice == 0) then
call m_wat(freq,tt,n_r,n_i)
else
call m_ice(freq,tt,n_r,n_i)
endif
m = cmplx(n_r,-n_i)
m0(:) = m
correct_for_rho = 0
if ((ice == 1) .and. (minval(rho_e) >= 0)) correct_for_rho = 1
! :: correct refractive index for ice density if needed
if (correct_for_rho == 1) then
f = rho_e/rho_ice
m0 = ((2+m0**2+2*f*(m0**2-1))/(2+m0**2+f*(1-m0**2)))**(0.5)
endif
! :: Mie calculations
sizep = (pi*D0)/wl
dqv(1) = 0.
do i=1,nsizes
call mieint(sizep(i), m0(i), one, dqv, qext(i), dqsc, qbsca(i), &
dg, xs1, xs2, dph, err)
end do
else
! // Mie table used
qext = qe
qbsca = qs
endif
! // eta_mie = 0.25*sum[qbsca*pi*D^2*N(D)*deltaD]
! <--------- eta_sum --------->
! // z0_eff = (wl^4/!pi^5)*(1./k2)*eta_mie
eta_sum = 0.
if (size(D0) == 1) then
eta_sum = qbsca(1)*(n(1)*1E6)*D0(1)**2
else
call avint(qbsca*N0*D0**2,D0,nsizes,D0(1),D0(size(D0,1)),eta_sum)
endif
eta_mie = eta_sum*0.25*pi
const = (wl**4/pi**5)*(1./k2)
z0_eff = const*eta_mie
! // kr = 0.25*cr*sum[qext*pi*D^2*N(D)*deltaD]
! <---------- k_sum --------->
k_sum = 0.
if (size(D0) == 1) then
k_sum = qext(1)*(n(1)*1E6)*D0(1)**2
else
call avint(qext*N0*D0**2,D0,nsizes,D0(1),D0(size(D0,1)),k_sum)
endif
cr = 10./log(10.)
kr = k_sum*0.25*pi*(1000.*cr)
! // z_ray = sum[D^6*N(D)*deltaD]
if (xr == 1) then
z0_ray = 0.
if (size(D0) == 1) then
z0_ray = (n(1)*1E6)*D0(1)**6
else
call avint(N0*D0**6,D0,nsizes,D0(1),D0(size(D0)),z0_ray)
endif
endif
! // convert to mm^6/m^3
z_eff = z0_eff*1E18 ! 10.*alog10(z0_eff*1E18)
z_ray = z0_ray*1E18 ! 10.*alog10(z0_ray*1E18)
end subroutine zeff