!wdt-caps sra
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE PCP_MXR_FERRIER             ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE pcp_mxr_ferrier (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 Ferrier et al (1995)
!  formulation:
!
!
!     For rain water:
!
!          18
!        10   * 720                              1.75
!  Zer = --------------------------- * (rho * qr)
!          1.75      0.75       1.75
!        pi     * N0r     * rhor
!
!
!     For dry snow (t <= 0 C):
!
!
!          18           2                     0.25
!        10  * 720 * |K|                * rhos
!                       ice                                    1.75
!  Zes = ----------------------------------------- * (rho * qs)     t <= 0 C
!          1.75         2          0.75       2
!        pi        * |K|      * N0s     * rhoi
!                     water
!
!
!     For wet snow (t >= 0 C):
!
!
!              18
!            10   * 720                                 1.75
!  Zes =     ----------------------------   * (rho * qs)            t >  0 C
!              1.75      0.75       1.75
!            pi     * N0s     * rhos
!
!
!     For hail water:
!
!
!          /   18                       \ 0.95
!         /  10   * 720                  \              1.6625
!  Zeh = |   ---------------------------- | * (rho * qh)
!         \    1.75      0.75       1.75 /
!          \ pi     * N0h     * rhoh    /
!
!  Here Zx (mm**6/m**3, x=r,s,h), and dBZ = 10log10 (Zx).
!  rho represents the air density, rhor,rhos,rhoh are the density of
!  rain, snow and hail respectively. Other variables are all constants
!  for this scheme, see below.
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: (Donghai Wang and Eric Kemp)
!  07/20/2000
!
!  MODIFICATION HISTORY:
!
!  11/09/2000 Keith Brewster
!  Moved some parameters with real-valued exponentiation to be
!  computed at runtime due to compiler complaint.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  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)
  REAL :: rair                     ! Gas constant (J/deg/kg)
  PARAMETER (rair = 287.04)
  REAL :: thresh_ref
  PARAMETER (thresh_ref = 0.0)
  INTEGER :: pcptype
  REAL :: zerf,zesnegf,zesposf,zehf


  REAL :: a4over7
  PARAMETER (a4over7 = 4.0/7.0)
  REAL :: zehpower
  PARAMETER (zehpower = 1.0/1.6625)
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k, iarg
  REAL :: ze,rho,tc
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Intiailize constant factors in the Ze terms for rain, snow and hail,
!  respectively, in Ferrier.
!
!-----------------------------------------------------------------------
!
  istatus=0
  zerf = (1.0E+18 * 720 )                                    &
             / (3.1415926**1.75 * 8.0E+06 ** 0.75 * 1000 ** 1.75 )      &
             / ( 1000. ** 1.75)

  zesnegf = ((1.0E+18 * 720 * 0.176 * 100 ** 0.25 )          &
             /(3.1415926**1.75 * 0.93 * 3.0E+06 ** 0.75 * 917 ** 2))    &
             / ( 1000. ** 1.75)

  zesposf = ((1.0E+18 * 720 )                                &
             / (3.1415926**1.75 * 3.0E+06 ** 0.75 * 100 ** 1.75))       &
             / ( 1000. ** 1.75)

  zehf = (((1.0E+18 * 720 )                                  &
      / (3.1415926**1.75 * 4.0E+04 ** 0.75 * 913 ** 1.75 )) ** 0.95)    &
      / ( 1000. ** 1.6625)

!-----------------------------------------------------------------------
!
!  Compute the precip mixing ratio in g/kg from radar reflectivity
!  factor following Ferrier et al (1995).
!
!-----------------------------------------------------------------------
!
  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.
          rho = p_3d(i,j,k)/rair/t_3d(i,j,k)
          ze = 10.0**(0.1*ref_3d(i,j,k))
          iarg = cldpcp_type_3d(i,j,k)
          pcptype = iarg/16              ! precip. type
          tc = t_3d(i,j,k) - 273.15

          IF (pcptype == 0 .OR. rho <= 0.0) THEN       ! no precip
            PRINT*,'+++ NOTE: radar echo though no precip. +++'
          ELSE IF (pcptype == 1.OR.pcptype == 3) THEN   ! rain or Z R
!        print*,'calculating qr...'
            qr_3d(i,j,k) = (ze / zerf ) ** a4over7 / rho
          ELSE IF (pcptype == 2) THEN                   ! snow
            IF (tc <= 0.0) THEN     !dry snow
!          print*,'calculating dry qs...'
              qs_3d(i,j,k) = (ze / zesnegf) ** a4over7 / rho
            ELSE IF (tc > 0.0) THEN     !wet snow
!          print*,'calculating wet qs...'
              qs_3d(i,j,k) = (ze / zesposf) ** a4over7 / rho
              qr_3d(i,j,k) = (ze / zerf ) ** a4over7 / rho
              qr_3d(i,j,k) = qr_3d(i,j,k) - qs_3d(i,j,k)
            END IF
          ELSE IF (pcptype == 4) THEN   ! sleet
!          print*,'calculating wet sleet qs...'
            qs_3d(i,j,k) = (ze / zesposf) ** a4over7 / rho
            qr_3d(i,j,k) = (ze / zerf ) ** a4over7 / rho
            qr_3d(i,j,k) = qr_3d(i,j,k) - qs_3d(i,j,k)
          ELSE IF (pcptype == 5) THEN   ! hail
!      else if (pcptype.eq.4.or.pcptype.eq.5) then   ! hail or sleet
!        print*,'calculating qh...'
            qh_3d(i,j,k) = (ze / zehf) ** zehpower / rho
          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
  PRINT*,'Finish Ferrier ...'
!
!-----------------------------------------------------------------------
!
  istatus = 1
!
  RETURN
END SUBROUTINE pcp_mxr_ferrier