source: LMDZ5/trunk/libf/phylmd/cosp/zeff.F90 @ 3132

Last change on this file since 3132 was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

  • mise a jour des sources pour ISCCP, CALIPSO et PARASOL
  • prise en compte des changements de phases pour les nuages (Calipso)
  • rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 4.7 KB
Line 
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 $
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
19!   [tt]        hydrometeor temperature (K)
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 :: &
46  correct_for_rho        ! correct for density flag
47  real*8, dimension(nsizes) :: &
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
56  real*8 :: &
57  wl, &                  ! wavelength (m)
58  cr                            ! kr(dB/km) = cr * kr(1/km)
59  complex*16 :: &
60  m                 ! complex index of refraction of bulk form
61  complex*16, dimension(nsizes) :: &
62  m0                ! complex index of refraction
63 
64  integer*4 :: i,one
65  real*8 :: pi
66  real*8 :: eta_sum, eta_mie, const, z0_eff, z0_ray, k_sum, &
67            n_r, n_i, dqv(1), dqsc, dg, dph(1)
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
77  D0 = d*1E-6            ! m
78  N0 = n*1E12            ! 1/(m^3 m)
79  wl = 2.99792458/(freq*10)   ! m
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
132    xtemp = qbsca*N0*D0**2
133    call avint(xtemp,D0,nsizes,D0(1),D0(size(D0,1)),eta_sum)
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
146    xtemp = qext*N0*D0**2
147    call avint(xtemp,D0,nsizes,D0(1),D0(size(D0,1)),k_sum)
148  endif
149  cr = 10./log(10.)
150  kr = k_sum*0.25*pi*(1000.*cr)
151     
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
158      xtemp = N0*D0**6
159      call avint(xtemp,D0,nsizes,D0(1),D0(size(D0)),z0_ray)
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
Note: See TracBrowser for help on using the repository browser.