SUBROUTINE pcp_mxr (nx,ny,nz,t_3d,p_3d ,ref_3d &
,cldpcp_type_3d &
,qr_3d,qs_3d,qg_3d,istatus )
!
!$$$ subprogram documentation block
! . . . .
! subprogram: pcp_mxr calculates hydrometeor mixing ratios based on Kessler radar reflectivity equations
!
! PRGMMR: ORG: DATE:
!
! ABSTRACT:
! This subroutine calculate precipitation based on Kessler radar reflectivity equations
!
! PROGRAM HISTORY LOG:
! 2009-01-20 Hu Add NCO document block
!
!
! input argument list:
! nx - no. of lons on subdomain (buffer points on ends)
! ny - no. of lats on subdomain (buffer points on ends)
! nz - no. of levels
! t_3d - 3D background temperature (K)
! p_3d - 3D background pressure (hPa)
! ref_3d - 3D reflectivity in analysis grid (dBZ)
! cldpcp_type_3d - 3D precipitation type
!
! output argument list:
! qr_3d - rain mixing ratio (g/kg)
! qs_3d - snow mixing ratio (g/kg)
! qg_3d - graupel/hail mixing ratio (g/kg)
! istatus -
!
! USAGE:
! INPUT FILES:
!
! OUTPUT FILES:
!
! REMARKS:
! Old documents from CAPS
!
!-----------------------------------------------------------------------
!
! 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 graupel/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.
!
!-----------------------------------------------------------------------
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90
! MACHINE: Linux cluster (WJET)
!
!$$$
!
!_____________________________________________________________________
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
use kinds, only: r_single,i_kind, r_kind
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
integer(i_kind),intent(in) :: nx,ny,nz ! Model grid size
!
REAL(r_kind), intent(in) :: ref_3d(nx,ny,nz)! radar reflectivity (dBZ)
real(r_single),intent(in) :: t_3d(nx,ny,nz) ! Temperature (deg. Kelvin)
real(r_single),intent(in) :: p_3d(nx,ny,nz) ! Pressure (Pascal)
integer(i_kind),intent(in):: cldpcp_type_3d(nx,ny,nz) ! cloud/precip type field
!
! OUTPUT:
INTEGER(i_kind),intent(out) :: istatus
!
REAL(r_single),intent(out) :: qr_3d(nx,ny,nz)! rain mixing ratio in (g/kg)
REAL(r_single),intent(out) :: qs_3d(nx,ny,nz)! snow/sleet/frz-rain mixing ratio
! in (g/kg)
REAL(r_single),intent(out) :: qg_3d(nx,ny,nz)! graupel/hail mixing ratio in (g/kg)
!
! LOCAL:
REAL(r_kind) :: a,b,c,d ! Coef. for Z-qr relation.
PARAMETER (a=17300.0_r_kind, b=7.0/4.0_r_kind)
PARAMETER (c=38000.0_r_kind, d=2.2_r_kind)
REAL(r_kind) :: rair ! Gas constant (J/deg/kg)
PARAMETER (rair = 287.04_r_kind)
REAL(r_kind) :: thresh_ref
PARAMETER (thresh_ref = 0.0_r_kind)
INTEGER(i_kind) :: pcptype
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER(i_kind) :: i,j,k, iarg
REAL(r_kind) :: arg,rhobar,br,dr
PARAMETER (br=1.0_r_kind/b, dr=1.0_r_kind/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 = 2,ny-1
DO i = 2,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_r_kind**(0.1_r_kind*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
qg_3d(i,j,k) = (arg/c)**dr/rhobar
ELSE ! unknown
PRINT*,'+++ NOTE: unknown precip type. +++'
END IF
ELSE
qr_3d(i,j,k) = 0._r_kind
qs_3d(i,j,k) = 0._r_kind
qg_3d(i,j,k) = 0._r_kind
END IF
END DO ! k
END DO ! i
END DO ! j
!
!-----------------------------------------------------------------------
!
istatus = 1
!
RETURN
END SUBROUTINE pcp_mxr
!
SUBROUTINE pcp_mxr_ferrier (nx,ny,nz,t_3d,p_3d ,ref_3d &
,cldpcp_type_3d &
,qr_3d,qs_3d,qg_3d,istatus,mype )
!
!$$$ subprogram documentation block
! . . . .
! subprogram: pcp_mxr calculate hydrometeor type based on ferrier radar reflectivity equations
!
! PRGMMR: ORG: DATE:
!
! ABSTRACT:
! This subroutine calculate precipitation based on ferrier radar reflectivity equations
!
! PROGRAM HISTORY LOG:
! 2009-01-20 Hu Add NCO document block
!
!
! input argument list:
! nx - no. of lons on subdomain (buffer points on ends)
! ny - no. of lats on subdomain (buffer points on ends)
! nz - no. of levels
! t_3d - 3D background temperature (K)
! p_3d - 3D background pressure (hPa)
! ref_3d - 3D reflectivity in analysis grid (dBZ)
! cldpcp_type_3d - 3D precipitation type
!
! output argument list:
! qr_3d - rain mixing ratio (g/kg)
! qs_3d - snow mixing ratio (g/kg)
! qg_3d - graupel/hail mixing ratio (g/kg)
! istatus -
!
! USAGE:
! INPUT FILES:
!
! OUTPUT FILES:
!
! REMARKS:
! Old document from CAPS
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! 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 * qg)
! \ 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.
!
! 04/07/2003 Keith Brewster
! Restructured code to make more tractable.and consistent with
! the reflec_ferrier subroutine.
!
!-----------------------------------------------------------------------
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90
! MACHINE: Linux cluster (WJET)
!
!$$$
!
!_____________________________________________________________________
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
use kinds, only: r_single,i_kind, r_kind
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
INTEGER(i_kind),intent(in) :: nx,ny,nz ! Model grid size
!
REAL(r_kind), intent(inout) :: ref_3d(nx,ny,nz)! radar reflectivity (dBZ)
REAL(r_single), intent(in) :: t_3d(nx,ny,nz) ! Temperature (deg. Kelvin)
REAL(r_single), intent(in) :: p_3d(nx,ny,nz) ! Pressure (Pascal)
INTEGER(i_kind),intent(in) :: cldpcp_type_3d(nx,ny,nz) ! cloud/precip type field
INTEGER(i_kind),intent(in) :: mype
!
! OUTPUT:
INTEGER(i_kind),intent(out):: istatus
!
REAL(r_single),intent(out) :: qr_3d(nx,ny,nz) ! rain mixing ratio in (g/kg)
REAL(r_single),intent(out) :: qs_3d(nx,ny,nz) ! snow/sleet/frz-rain mixing ratio
! in (g/kg)
REAL(r_single),intent(out) :: qg_3d(nx,ny,nz) ! graupel/hail mixing ratio
! in (g/kg)
!
REAL(r_kind),PARAMETER :: ki2 = 0.176_r_kind ! Dielectric factor for ice if other
! than melted drop diameters are used.
REAL(r_kind),PARAMETER :: kw2=0.93_r_kind ! Dielectric factor for water.
REAL(r_kind),PARAMETER :: m3todBZ=1.0E+18_r_kind ! Conversion factor from m**3 to
! mm**6 m**-3.
REAL(r_kind),PARAMETER :: Zefact=720.0_r_kind ! Multiplier for Ze components.
REAL(r_kind),PARAMETER :: lg10div=0.10_r_kind ! Log10 multiplier (1/10)
REAL(r_kind),PARAMETER :: pi=3.1415926_r_kind! Pi.
REAL(r_kind),PARAMETER :: N0r=8.0E+06_r_kind ! Intercept parameter in 1/(m^4) for rain.
REAL(r_kind),PARAMETER :: N0s=3.0E+06_r_kind ! Intercept parameter in 1/(m^4) for snow.
REAL(r_kind),PARAMETER :: N0h=4.0E+04_r_kind ! Intercept parameter in 1/(m^4) for graupel/hail.
REAL(r_kind),PARAMETER :: N0xpowf=3.0/7.0_r_kind ! Power to which N0r,N0s & N0h are
! raised.
REAL(r_kind),PARAMETER :: K2powf=4.0/7.0_r_kind ! Power to which K-squared
! of ice, water are raised
REAL(r_kind),PARAMETER :: zkpowf=4.0/7.0_r_kind ! Power to which Zk is raised
REAL(r_kind),PARAMETER :: zepowf=4.0/7.0_r_kind ! Power to which Ze is raised
REAL(r_kind),PARAMETER :: zehpowf=(4.0/7.0)*1.0526_r_kind ! Power to which Zeh is raised
REAL(r_kind),PARAMETER :: rhoi=917._r_kind ! Density of ice (kg m**-3)
REAL(r_kind),PARAMETER :: rhor=1000._r_kind ! Density of rain (kg m**-3)
REAL(r_kind),PARAMETER :: rhos=100._r_kind ! Density of snow (kg m**-3)
REAL(r_kind),PARAMETER :: rhoh=913._r_kind ! Density of graupel/hail (kg m**-3)
REAL(r_kind),PARAMETER :: rhoipowf=8.0/7.0_r_kind ! Power to which rhoi is raised.
REAL(r_kind),PARAMETER :: rhospowf=1.0/7.0_r_kind ! Power to which rhos is raised.
REAL(r_kind), PARAMETER :: rd=287.0_r_kind ! Gas constant for dry air (m**2/(s**2*K))
REAL(r_kind), PARAMETER :: thresh_ref = 0.0_r_kind
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER(i_kind) :: i,j,k, iarg
INTEGER(i_kind) :: pcptype
REAL(r_kind) :: zkconst,zerf,zesnegf,zesposf,zehf,rfract
REAL(r_kind) :: ze,zer,zeh,zes,rho,tc
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Intiailize constant factors in the Ze terms for rain, snow and graupel/hail,
! respectively, in Ferrier.
!
! These are the inverse of those presented in the reflec_ferrier function.
!
!-----------------------------------------------------------------------
!
istatus=0
zkconst = (Zefact*m3todBZ) ** zkpowf
zerf=1000._r_kind*(pi * (N0r**N0xpowf) * rhor )/zkconst
zesnegf=1000._r_kind*(pi*(kw2**k2powf)*(N0s**N0xpowf)*(rhoi**rhoipowf)) / &
( zkconst * (ki2**k2powf) * (rhos**rhospowf) )
zesposf=1000._r_kind*( pi * (N0s**N0xpowf) * rhos) / zkconst
zehf=1000._r_kind*( pi * (N0h**N0xpowf) * rhoh) / zkconst
!-----------------------------------------------------------------------
!
! Compute the precip mixing ratio in g/kg from radar reflectivity
! factor following Ferrier et al (1995).
!
!-----------------------------------------------------------------------
!
DO k = 2,nz-1
DO j = 2,ny-1
DO i = 2,nx-1
IF (ref_3d(i,j,k) >= thresh_ref) THEN ! valid radar refl.
rho = p_3d(i,j,k)/(rd*t_3d(i,j,k))
ze = 10.0_r_kind**(0.1_r_kind*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_r_kind
!mhu temporal fix
IF (tc <= 0.0_r_kind) THEN
qs_3d(i,j,k) = zesnegf * (ze**zepowf) / rho
qr_3d(i,j,k) = 0.0_r_kind
ELSE IF (tc < 5.0_r_kind) THEN !wet snow
rfract=0.20_r_kind*tc
zer=rfract*ze
zes=(1.-rfract)*ze
! qs_3d(i,j,k) = zesposf * (zes**zepowf) / rho
! qr_3d(i,j,k) = zerf * (zer**zepowf) / rho
qs_3d(i,j,k) = zesnegf * (zes**zepowf) / rho
qr_3d(i,j,k) = zerf * (zer**zepowf) / rho
else
qr_3d(i,j,k) = zerf * (ze**zepowf) / rho
qs_3d(i,j,k) = 0.0_r_kind
ENDIF
cycle
!mhu
IF (pcptype == 1) THEN ! rain
qr_3d(i,j,k) = zerf * (ze**zepowf) / rho
ELSE IF (pcptype == 2) THEN ! snow
IF (tc <= 0.0_r_kind) THEN !dry snow
qs_3d(i,j,k) = zesnegf * (ze**zepowf) / rho
ELSE IF (tc < 5.0_r_kind) THEN !wet snow
rfract=0.20_r_kind*tc
zer=rfract*ze
zes=(1.-rfract)*ze
qs_3d(i,j,k) = zesposf * (zes**zepowf) / rho
qr_3d(i,j,k) = zerf * (zer**zepowf) / rho
ELSE
qr_3d(i,j,k) = zerf * (ze**zepowf) / rho
END IF
ELSE IF (pcptype == 3) THEN ! ZR
qr_3d(i,j,k) = zerf * (ze**zepowf) / rho
ELSE IF (pcptype == 4) THEN ! sleet
IF (tc <= 0.0_r_kind) THEN ! graupel/hail category
qg_3d(i,j,k) = zehf * (ze**zehpowf) / rho
ELSE IF( tc < 10._r_kind ) THEN
rfract=0.10_r_kind*tc
zer=rfract*ze
zeh=(1.-rfract)*ze
qr_3d(i,j,k) = zerf * (zer**zepowf) / rho
qg_3d(i,j,k) = zehf * (zeh**zehpowf) / rho
ELSE
qr_3d(i,j,k) = zerf * (ze**zepowf) / rho
END IF
ELSE IF (pcptype == 5) THEN ! graupel/hail
qg_3d(i,j,k) = zehf * (ze**zehpowf) / rho
ELSE ! unknown
IF (tc <= 0.0_r_kind) THEN !dry snow
qs_3d(i,j,k) = zesnegf * (ze**zepowf) / rho
ELSE IF ( tc < 5.0_r_kind ) THEN !wet snow
rfract=0.20_r_kind*tc
zer=rfract*ze
zes=(1.-rfract)*ze
qs_3d(i,j,k) = zesposf * (zes**zepowf) / rho
qr_3d(i,j,k) = zerf * (zer**zepowf) / rho
ELSE ! rain
qr_3d(i,j,k) = zerf * (ze**zepowf) / rho
END IF
END IF
ELSE
qr_3d(i,j,k) = -999._r_kind
qs_3d(i,j,k) = -999._r_kind
qg_3d(i,j,k) = -999._r_kind
END IF
END DO ! k
END DO ! i
END DO ! j
! PRINT*,'Finish Ferrier ...'
!
!-----------------------------------------------------------------------
!
istatus = 1
!
RETURN
END SUBROUTINE pcp_mxr_ferrier