c-----------------------------------------------------------------------c
c
c Radar reflectivity is obtained from the following equation.
c
c   Z = integral[0 to inf] N(D) D^6 dD  (from Rogers and Yau, Cloud Physics book)
c
c If N(D) = No exp(-slp*D), then the solution for Z is
c
c   Z = No * 6!/slp^7  =  720. No/slp^7
c
c So if your model uses the exponential decay size distribution of 
c  precipitation, then you can use this equation.  Our subroutine for 
c  this follows (there are some details regarding snow and hail versus rain).
c
c If your model uses another approach for N(D), then I suggest you follow
c  the first equation and perform the integration accordingly.  Be sure
c  to account for different refraction and Mie scattering for snow and/or
c  hail.  This is used in the following routine.
c
c-----------------------------------------------------------------------c
      subroutine compute_reflec( density, temperature, ref,
     &                        qr_in, qs_in, qh_in, nx, ny, nz )
      implicit none
c Input
      integer nx,ny,nz
      real density(nx,ny,nz),            ! air density (kg/m3)
     _     temperature(nx,ny,nz)         ! temperature (deg C)
      real qr_in(nx,ny,nz),              ! rain mixing ratio (kg/kg)
           qs_in(nx,ny,nz),              ! snow mixing ratio (kg/kg)
           qh_in(nx,ny,nz)               ! hail mixing ratio (kg/kg)
c Output
      real ref(nx,ny,nz)                 ! reflectivity (dBZ)
C NCLEND

c From Lou Wicker, who stole it from Jerry Straka.
c 17 December 2001
c Computes radar reflectivity from model water variables.
c
c Algorithm is derived from Jerry Strakas code
c NOTE:  An important trick to get these calculations to work is to
c use double precision variables for all the intermediate calculations.
c If you dont, you get overflow and underflow error when computing
c the logrithms.  Ugh. 
c
c Units are MKS, and for most accurate results, make sure that the
c number concentrations and densities of rain, snow and hail are the
c same as the model used in producing the fields.
c
c Modified by Wei Wang, April 22, 2003
c Make the code able to take correct account of snow in ncep3 and ncep5 schemes
c
c-----------------------------------------------------------------------

      integer i,j,k

      real*8 qr, qs, qh, den, temp, reflect
      real*8 xcnoh,xcnos,zsdryc,zswetc,zhdryc,zhwetc,dadr,dads,dadh
      real*8 cnoh,cnos,cnor,cr1,hwdn,swdn,rwdn,hwdnsq,swdnsq,rwdnsq
      real*8 dhmin,dielf,pie,qrmin,qsmin,qhmin,tfr,tfrh,zrc
      real*8 tem,gmh,gms,gmr,zrain,zswet,zhwet,zsdry,zhdry

      real*8 reflectmin
      real   dbz

        parameter ( reflectmin = 0.001 )

c
c constants
c
c This code is based on the size distributions of the precipitation fields
c  (rain, snow, hail or graupel) to be of the form:
c  N = No exp(-slp*D) dD     where D is diameter
c
      cnoh = 4.0e+04              ! No for hail (/m^4)
      cnor = 8.0e+06              ! No for rain (/m^4)
      cnos = 3.0e+06              ! No for snow (/m^4)
      cr1  = 7.2e+20              ! 6!=720 and unit conversion=1e18
      hwdn = 660.0                ! density of hail -- change this to your model configuration
      swdn = 100.0                ! density of snow
      rwdn = 1000.0               ! density of rain
      hwdnsq = hwdn**2
      swdnsq = swdn**2
      rwdnsq = rwdn**2
      dhmin = 0.005
      dielf = 0.21/0.93           ! snow or hail to rain refractivity ratio
      pie   = 4.0*atan(1.0)
      qrmin = 1.0e-05
      qsmin = 1.0e-06
      qhmin = 1.0e-05
      tfr   = 273.16
      tfrh  = tfr - 8.0
      zrc   = cr1*cnor

c define double precision mixing ratios

      do k=1,nz
      do j=1,ny
      do i=1,nx

      den = density(i,j,k)
      temp = temperature(i,j,k)+tfr

      qr = amax1(0.,qr_in(i,j,k))
      qs = amax1(0.,qs_in(i,j,k))
      qh = amax1(0.,qh_in(i,j,k))


c adjust No for snow and hail  (not sure what the origin of this is)
      if ( temp .lt. tfr ) then
       xcnoh = cnoh*exp(-0.025*(temp-tfr))
       xcnos = cnos*exp(-0.038*(temp-tfr))
      else
       xcnoh = cnoh*exp(-0.075*(temp-tfr))
       xcnos = cnos*exp(-0.088*(temp-tfr))
      end if
c
c slope of size distribution except for mixing ratio (i.e. the constant part
c                                                      of slp parameter)
      dadh = ( den /(pie*hwdn*xcnoh) )**.25
      dads = ( den /(pie*swdn*xcnos) )**.25
      dadr = ( den /(pie*rwdn*cnor)  )**.25
c
c constants in the equation
      zhdryc = dielf*cr1*(hwdnsq/rwdnsq)*xcnoh
      zhwetc = cr1*(xcnoh)**.95
      zsdryc = dielf*cr1*swdnsq/rwdnsq*xcnos
      zswetc = cr1*xcnos

      zrain = 0.
      zswet = 0.
      zsdry = 0.
      zhwet = 0.
      zhdry = 0.
      gmr = 0.0
      gms = 0.0
      gmh = 0.0

c-----------------------------------------------------------------------
c Compute Z for each hydrometeor
c-----------------------------------------------------------------------

      reflect = reflectmin

c Rain
c..slope for rain
       if ( qr .gt. qrmin ) then
        gmr = dadr*(qr)**.25
       end if

       zrain = zrc*gmr**7

c Snow
c..slope for snow
       if ( qs .gt. qsmin ) then
        gms = dads*(qs)**.25
       endif

       if ( temp .lt. tfr ) then
         zsdry = zsdryc*gms**7
       else
         zswet = zswetc*gms**7
       end if

c Hail/graupel
c..slope for hail/graupel
       if ( qh .gt. qhmin ) then
        gmh = dadh*(qh)**.25
       endif

c..computation for dry and wet hail reflectivity including mie scattering
c  parameterization (Smith, 1975) 
       if ( temp .gt. tfr ) then
        zhwet = (zhwetc*(gmh**7))**.95
       else
        zhdry = zhdryc*gmh**7
       endif

       
c-----------------------------------------------------------------------
c Compute final reflectivity value
c-----------------------------------------------------------------------
      dbz = zrain + zswet + zsdry + zhwet + zhdry

      if( dbz .gt. 0. ) then
       reflect = 10.0 * alog10(dbz)
       ref(i,j,k) = max(reflect, reflectmin)
      else
       ref(i,j,k) = 0.
      endif

      enddo
      enddo
      enddo


      return
      end