!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE PCP_MXR                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE pcp_mxr (nx,ny,nz,t_3d,p_3d ,ref_3d                          & 1
           ,cldpcp_type_3d                                              &
           ,qr_3d,qs_3d,qh_3d,istatus )

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Perform 3D precipitation mixing ratio (in g/kg) analysis using
!  radar reflectivity data. For rain water, using Kessler (1969)
!  formula:
!            qr(g/kg) = a*(rho*arg)**b                  (1)
!
!  Here arg = Z (mm**6/m**3), and dBZ = 10log10 (arg).
!  Coeffcients a=17300.0, and b=7/4.
!  rho represents the air density.
!
!  For snow and hail, using Rogers and Yau (1989) formula:
!
!            qs(g/kg) = c*(rho*arg)**d                  (2)
!
!  where, c=38000.0,  d=2.2
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Jian Zhang)
!  06/13/96
!
!  MODIFICATION HISTORY:
!  07/30/97 (J. Zhang)
!           Added precipitation type in the argument list so that
!           mixing ratios of different precip. types can be computed.
!  09/04/97 (J. Zhang)
!            Changed the radar echo thresholds for inserting precip.
!            from radar reflectivities.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  INPUT:
  INTEGER :: nx,ny,nz              ! Model grid size
!
  REAL :: ref_3d(nx,ny,nz)  ! radar reflectivity (dBZ)
  REAL :: t_3d(nx,ny,nz)  ! Temperature (deg. Kelvin)
  REAL :: p_3d(nx,ny,nz)  ! Pressure (Pascal)

  INTEGER*2 cldpcp_type_3d(nx,ny,nz) ! cloud/precip type field
!
!  OUTPUT:
  INTEGER :: istatus
!
  REAL :: qr_3d(nx,ny,nz)     ! rain mixing ratio in (g/kg)
  REAL :: qs_3d(nx,ny,nz)     ! snow/sleet/frz-rain mixing ratio
                            ! in (g/kg)
  REAL :: qh_3d(nx,ny,nz)     ! hail mixing ratio in (g/kg)
!
!  LOCAL:
  REAL :: a,b,c,d                  ! Coef. for Z-qr relation.
  PARAMETER (a=17300.0, b=7.0/4.0)
  PARAMETER (c=38000.0, d=2.2)
  REAL :: rair                     ! Gas constant (J/deg/kg)
  PARAMETER (rair = 287.04)
  REAL :: thresh_ref
  PARAMETER (thresh_ref = 0.0)
  INTEGER :: pcptype
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k, iarg
  REAL :: arg,rhobar,br,dr
  PARAMETER (br=1.0/b, dr=1.0/d)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
  istatus=0
!
!-----------------------------------------------------------------------
!
!  Compute the precip mixing ratio in g/kg from radar reflectivity
!  factor following Kessler (1969) or Rogers and Yau (1989).
!
!-----------------------------------------------------------------------
!
  DO k = 1,nz-1
    DO j = 1,ny-1
      DO i = 1,nx-1
        IF (ref_3d(i,j,k) >= thresh_ref) THEN    ! valid radar refl.
          rhobar = p_3d(i,j,k)/rair/t_3d(i,j,k)
          arg = 10.0**(0.1*ref_3d(i,j,k))
          iarg = cldpcp_type_3d(i,j,k)
          pcptype = iarg/16              ! precip. type

          IF (pcptype == 0) THEN       ! no precip
            PRINT*,'+++ NOTE: radar echo though no precip. +++'
          ELSE IF (pcptype == 1.OR.pcptype == 3) THEN   ! rain or Z R
            qr_3d(i,j,k) = (arg/a)**br/rhobar
          ELSE IF (pcptype == 2) THEN                   ! snow
            qs_3d(i,j,k) = (arg/c)**dr/rhobar
          ELSE IF (pcptype == 4.OR.pcptype == 5) THEN   ! hail or sleet
            qh_3d(i,j,k) = (arg/c)**dr/rhobar
          ELSE                                          ! unknown
            PRINT*,'+++ NOTE: unknown precip type. +++'
          END IF
        ELSE
          qr_3d(i,j,k) = 0.
          qs_3d(i,j,k) = 0.
          qh_3d(i,j,k) = 0.
        END IF
      END DO ! k
    END DO ! i
  END DO ! j
!
!-----------------------------------------------------------------------
!
  istatus = 1
!
  RETURN
END SUBROUTINE pcp_mxr