source: LMDZ5/branches/LMDZ5-dev032011/libf/cosp/zeff.F90 @ 2448

Last change on this file since 2448 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

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