function retrieve_dsd_parms_hid,radar,site,year,month,day,DPQC=dpqc
;+
;  This version utilizes the HID classification of each pixel to
;  determine whether or not to assign a rain rate.  
;
;- 
    COMMON misc,TRUE,FALSE
    COMMON graphics,USE_Z_BUFFER,TC,XS,YS
    COMMON graphic_parms,scale,bottom,ncolors
    COMMON colors,black,purple,blue,cyan,green,yellow,orange,red,white
    COMMON symbols,plot,asterisk,dot,diamond,triangle,square,exx
    COMMON linestyles,solid,dotted,dashed,dashdot,dashdotdot,longdash
    COMMON dirs,top_dir,img_dir
    if(NOT KEYWORD_SET(DPQC)) then DPQC=false

;
; *** Disdrometer-based Z-R coefficients
;
    a1 = 0.024285 & b1 = 0.6895 ; From Tokay DSD @ Kwaj
    a1 = 0.017    & b1 = 0.7143 ; FACE Z=300R^1.4
;    a1 = 0.036 & b1 = 0.625     ; Marshall-Palmer Z=200R^1.6

    a_dsd = (1.0D/a1)^(1./b1)
    b_dsd = (1.0D/b1)
;
; *** Get the necessary volumes to calculate parameters.
;
    dbz_field = 'DZ'
    a = where(radar.volume.h.field_type eq 'CZ',c)
    if(c eq 1) then dbz_field='CZ'
    dbz_vol = rsl_get_volume(radar,dbz_field)

    zdr_vol = rsl_get_volume(radar,'DR')
    kdp_vol = rsl_get_volume(radar,'KD')
    rhv_vol = rsl_get_volume(radar,'RH')
    rr1_vol = rsl_get_volume(radar,'RR')
    hid_vol = rsl_get_volume(radar,'FH')
;
; *** Get geometry info from radar structure.
;
    missing = radar.volume[0].h.no_data_flag
    nsweeps = radar.volume[0].h.nsweeps
    nrays = radar.volume[0].sweep[0].h.nrays
    nbins = radar.volume[0].sweep[0].ray[0].h.nbins
    method = intarr(nsweeps,nrays,nbins) + missing
;
; *** The DROPS based UF already has Chandra's RR volume, so create
; *** a new one... RP (for polzr).
;
    rsl_add_volume,radar,rr1_vol,field='RP'
    rrp_vol = rsl_get_volume(radar,'RR')
    rrp_vol.h.field_type = 'RP'
    for iswp=0,nsweeps-1 do rrp_vol.sweep[iswp].h.field_type = 'RP'
    
    for iswp=0,nsweeps-1 do begin
        range = rsl_get_range_from_sweep(dbz_vol.sweep[iswp])
;
; *** Compute the rain rate at each radar pixel.
;
        for iray=0,nrays-1 do begin
            for ibin=0L,nbins-1-1 do begin
;
; *** Extract the dbz, zdr and kdp values at this gate.
;
                dbz = dbz_vol.sweep[iswp].ray[iray].range[ibin]
                zdr = zdr_vol.sweep[iswp].ray[iray].range[ibin]
                kdp = kdp_vol.sweep[iswp].ray[iray].range[ibin]
                rhv = rhv_vol.sweep[iswp].ray[iray].range[ibin]
                rr1 = rr1_vol.sweep[iswp].ray[iray].range[ibin]
                hid = hid_vol.sweep[iswp].ray[iray].range[ibin]
;
; *** Check to see if the HID is ice species or liquid species.
; *** If ice, then do not calculate RR at this bin.
; ***
; *** HID Species
; ***     UC=0, DZ=1, RN=2, CR=3, DS=4, WS=5
; ***     VI=6, LDG=7, HDG=8, HA=9 and BD=10
;
                if(hid eq 3 or hid eq 4 or hid eq 6 or $
                   hid eq 7 or hid eq 8) then begin
                    rr = 0
                    goto,next_bin
                endif
skip_hid_check:
;
; *** 1) Zero out any negative Kdp
; *** 2) Set maximum dBZ
;
                if(kdp lt 0) then kdp = 0
                if(dbz ge 53) then dbz = 53
;
; *** Default Z-R
;
                if(dbz gt -20) then begin
                    x1 = 0.1D*(dbz - 10D*alog10(a_dsd))/b_dsd
                    rrp_vol.sweep[iswp].ray[iray].range[ibin] = 10D^(x1)
                    method[iswp,iray,ibin] = 1
                endif else begin
                    rrp_vol.sweep[iswp].ray[iray].range[ibin] = missing
                    method[iswp,iray,ibin] = 0
                    goto,next_bin
                endelse
;
; *** Check to see if DP variables can be used to calculate 
; *** the rain rate. We need D0, Nw and mu to calculate the 
; *** Polarimetric ZR
;
                zh = 10.^(0.1*dbz) ; Linear Z [mm^6 m^-3]
                xi_dr = 10.^(0.1*zdr) ; Linear Units
                d0 = 0
                dm = 0
                Nw = 0
                mu = 0
                beta = 0
                logNw = missing
;
; *** Light rain rates with noisy Zdr
;
                if(zdr ge -0.5 and zdr lt 0.2) then begin
                    d0 = 0.6096*(zh^0.0516)*(xi_dr^3.111)
                    nw = (21*zh)/d0^7.3529
                    mu = 3.0
                    dm = d0 * ((4+mu)/(3.67+mu))
                    method[iswp,iray,ibin] = 2
                    goto,calc_rain_rate
                endif
;
; *** Light rain rates with modest kdp
;
                if(kdp lt 0.3 and zdr ge 0.2 and zdr lt 0.5) then begin
                    x1 = ((zdr - 0.2)/0.3) * 1.81 * zdr^(0.486)
                    x2 = ((0.5-zdr)/0.3) * 0.6096*zh^(0.0516) * xi_dr^(3.111)
                    d0 = x1 + x2
                    Nw = (21*zh)/D0^7.3529
                    mu = 3.0
                    dm = d0 * ((4+mu)/(3.67+mu))
                    method[iswp,iray,ibin] = 3
                    goto,calc_rain_rate
                endif
;
; *** Moderate rain rates
;
                if(kdp lt 0.3 and zdr ge 0.5) then begin
                    d0 = 1.81*zdr^0.486
                    Nw = (21*zh)/d0^7.3529
                    mu = 3.0
                    dm = d0 * ((4+mu)/(3.67+mu))
                    method[iswp,iray,ibin] = 4
                    goto,calc_rain_rate
                endif
;
; *** Heavy rain rates
;
                if(dbz ge 35 and kdp ge 0.3 and zdr ge 0.2) then begin
                    
                    beta = 2.08*zh^(-0.365) * kdp^(0.38) * xi_dr^(0.965)
                    
                    a1 = 0.56 & b1 = 0.064 & c1 = 0.024*beta^(-1.42)
                    d0 = a1 * zh^b1 * xi_dr^c1
                    
                    a2 = 3.29 & b2 = 0.058 & c2 = -0.023 * beta^(-1.389)
                    logNw = a2 * zh^b2 * xi_dr^c2
                    nw = 10^(logNw)
                    
                    a3 = 203. * beta^(1.89)
                    b3 = 2.23 * beta^(0.0388)
                    c3 = 3.16 * beta^(-0.0463)
                    d3 = 0.374 * beta^(-0.355)
                    
                    x1 = a3 * (d0^b3)/(xi_dr-1)
                    x2 = c3 * (xi_dr^d3)
                    mu = x1 - x2
                    dm = d0 * ((4+mu)/(3.67+mu))

                    method[iswp,iray,ibin] = 5
                    goto,calc_rain_rate
                endif
calc_rain_rate:                            
;
; *** If we are unable to calculate the DSD parameters, then use
; *** the disdrometer and PMM-based estimates.
;
                if(nw eq missing) then begin
                    rrp_vol.sweep[iswp].ray[iray].range[ibin] = missing
                    goto,next_bin
                endif
;
; *** Set all pixels with KDP eq missing and RHV eq missing to missing
;
                rhv_thresh = 0.75
                kdp_thresh = -1.0
                if(rhv lt rhv_thresh and kdp lt kdp_thresh) then begin
                    rrp_vol.sweep[iswp].ray[iray].range[ibin] = missing
                    goto,next_bin
                endif
;
; *** Calculate the coefficient a' in Z = a' * R^1.5 using the DSD
; *** parameters. First, calculate f(mu)
;
                rr = get_polzr_rainrate(dbz,d0,Nw,mu,a_prime,missing=missing)
                rrp_vol.sweep[iswp].ray[iray].range[ibin] = rr
                if(rr lt 0) then goto,next_bin

                year   = radar.h.year
                month  = radar.h.month
                day    = radar.h.day
                hour   = radar.h.hour
                minute = radar.h.minute
                second = radar.h.sec
                
                format = '(i4.4,1x,5(i2.2,1x),2(I3,1x),6(F9.2,1x),i1)'
next_bin:                
            endfor              ; Range bins
        endfor                  ; Azimuths
    endfor                      ; Sweeps
;
; *** Now place rrp_vol back into radar structure
;
    a = where(radar.volume.h.field_type eq 'RP',c)
    radar.volume[a[0]] = rrp_vol

    flag = 'OK'
    return,flag
    end