function retrieve_dp_rainrates,radar,site,year,month,day
;+
;  This version utilizes the HID classification of each pixel to
;  determine whether or not to assign a rain rate.  
; *** HID
; *** Check to see if the HID is ice or liquid species.
; *** If ice, then do not calculate RR at this bin.
; *** -1: Bad Data
; ***  0: Unclassified  [UC]
; ***  1: Drizzle       [DZ]
; ***  2: Rain          [RN]
; ***  3: Ice Crystals  [CR]
; ***  4: Dry Snow      [DS]
; ***  5: Wet Snow      [WS]
; ***  6: Vertical Ice  [VI]
; ***  7: Low Den Graup [LDG]
; ***  8: Hi Den Graup  [HDG]
; ***  9: Hail          [HA]
; *** 10: Big Drops     [BD]
;- 
    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
;
; *** Get the necessary volumes to calculate parameters.
;
    zfield = 'DZ'
    a = where(radar.volume.h.field_type eq 'CZ',c)
    if(c eq 1) then zfield='CZ'
    dbz_vol = rsl_get_volume(radar,zfield)

    zdr_vol = rsl_get_volume(radar,'DR')
    kdp_vol = rsl_get_volume(radar,'KD')
    rhv_vol = rsl_get_volume(radar,'RH')
    rr_vol  = rsl_get_volume(radar,'CZ')
    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
;
; *** RT: Rain estimator type
;
;    rain_field = 'RT'
;    rsl_add_volume,radar,rrt_vol,field=rain_field
;    rrt_vol = rsl_get_volume(radar,rain_field)
;    rrt_vol.h.field_type = rain_field
;    rrt_vol.sweep.h = rr_vol.sweep.h
;    rrt_vol.sweep.h.field_type = rain_field
;    rrt_vol.sweep.ray.h = rr_vol.sweep.ray.h
;    for iswp=0,nsweeps-1 do begin
;        rrt_vol.sweep[iswp].ray.range = rrt_vol.sweep[iswp].ray.range*0 + missing
;    endfor
;
; *** RG: Gatlin IFloodS Rain estimator 
;
;    rain_field = 'RG'
;    rsl_add_volume,radar,rrg_vol,field=rain_field
;    rrg_vol = rsl_get_volume(radar,rain_field)
;    rrg_vol.h.field_type = rain_field
;    rrg_vol.sweep.h = rr_vol.sweep.h
;    rrg_vol.sweep.h.field_type = rain_field
;    rrg_vol.sweep.ray.h = rr_vol.sweep.ray.h
;    for iswp=0,nsweeps-1 do begin
;        rrg_vol.sweep[iswp].ray.range = rrg_vol.sweep[iswp].ray.range*0 + missing
;    endfor
;
; *** RC: Cifelli et al. 2002
;
    rain_field = 'RC'
    rsl_add_volume,radar,rrc_vol,field=rain_field
    rrc_vol = rsl_get_volume(radar,rain_field)
    rrc_vol.h.field_type = rain_field
    rrc_vol.sweep.h = rr_vol.sweep.h
    rrc_vol.sweep.h.field_type = rain_field
    rrc_vol.sweep.ray.h = rr_vol.sweep.ray.h
    for iswp=0,nsweeps-1 do begin
        rrc_vol.sweep[iswp].ray.range = rrc_vol.sweep[iswp].ray.range*0 + missing
    endfor
    rain_field = 'MC'
    rsl_add_volume,radar,mc_vol,field=rain_field
    mc_vol = rsl_get_volume(radar,rain_field)
    mc_vol.h.field_type = rain_field
    mc_vol.sweep.h = rr_vol.sweep.h
    mc_vol.sweep.h.field_type = rain_field
    mc_vol.sweep.ray.h = rr_vol.sweep.ray.h
    for iswp=0,nsweeps-1 do begin
        mc_vol.sweep[iswp].ray.range = mc_vol.sweep[iswp].ray.range*0 + missing
    endfor

;
; *** RP: PolZR
;
    rain_field = 'RP'
    rsl_add_volume,radar,rrp_vol,field=rain_field
    rrp_vol = rsl_get_volume(radar,rain_field)
    rrp_vol.h.field_type = rain_field
    rrp_vol.sweep.h = rr_vol.sweep.h
    rrp_vol.sweep.h.field_type = rain_field
    rrp_vol.sweep.ray.h = rr_vol.sweep.ray.h
    for iswp=0,nsweeps-1 do begin
        rrp_vol.sweep[iswp].ray.range = rrp_vol.sweep[iswp].ray.range*0 + missing
    endfor
    rain_field = 'MP'
    rsl_add_volume,radar,mp_vol,field=rain_field
    mp_vol = rsl_get_volume(radar,rain_field)
    mp_vol.h.field_type = rain_field
    mp_vol.sweep.h = rr_vol.sweep.h
    mp_vol.sweep.h.field_type = rain_field
    mp_vol.sweep.ray.h = rr_vol.sweep.ray.h
    for iswp=0,nsweeps-1 do begin
        mp_vol.sweep[iswp].ray.range = mp_vol.sweep[iswp].ray.range*0 + missing
    endfor

;
; *** Process each sweep
;
    for iswp=0,nsweeps-1 do begin
        range = rsl_get_range_from_sweep(rr_vol.sweep[iswp])
        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]
                hid = hid_vol.sweep[iswp].ray[iray].range[ibin]
;
; *** 1) Zero out any negative Kdp
; *** 2) Set maximum dBZ
;
                if(kdp lt 0) then kdp = 0   ; Cifelli 2011
                if(dbz ge 53) then dbz = 53 ; Giagrande et al.
;
; *** If HID is: IC/3, DS/4, WS/5, VI/6, LDG/7, HDG/8, or  HA/9, do not
; *** assign a rainrate. 
; *** NOTE: meth1 and meth1 provide the rain estimator type which
; ***       we will place in a sweep structure as well.
;
                if(hid eq -1 or hid eq 3 or hid eq 4 or hid eq 5 or $
                   hid eq 6 or hid eq 7 or hid eq 8 or hid eq 9) then begin
                    rrc = [missing, missing]
;                    rrg = missing
                    rrp = [missing, missing]
;                    rrt = missing
                endif else begin
                    rrc = get_cif2002_rainrate(site,dbz,zdr,kdp,rhv,hid,meth1,MISSING=missing)
                    rrp = get_bringi_rainrate(dbz,zdr,kdp,rhv,hid,meth2,MISSING=missing)
;                    rrg = get_gatlin_rainrate(dbz,zdr,kdp,rhv,hid,meth3,MISSING=missing)
;                    rrt = meth2 ; PolZR only
                 endelse

;                rrt_vol.sweep[iswp].ray[iray].range[ibin] = rrt
;                print,rrc
                if(rrc[0] le 0) then rrc[0] = missing
                rrc_vol.sweep[iswp].ray[iray].range[ibin] = rrc[0]
		mc_vol.sweep[iswp].ray[iray].range[ibin] = rrc[1]
;		print,mc_vol.sweep[iswp].ray[iray].range[ibin]

;                if(rrg le 0) then rrg = missing
;                rrg_vol.sweep[iswp].ray[iray].range[ibin] = rrg

                if(rrp[0] le 0) then rrp[0] = missing
                rrp_vol.sweep[iswp].ray[iray].range[ibin] = rrp[0]
		mp_vol.sweep[iswp].ray[iray].range[ibin] = rrp[1]
;               print,'rrp =  ',rrp_vol.sweep[iswp].ray[iray].range[ibin]
;               print,'mp =  ',mp_vol.sweep[iswp].ray[iray].range[ibin]

            endfor              ; Range bins
        endfor                  ; Azimuths
    endfor                      ; Sweeps
;
; *** Now place rrN_vol back into radar structure
;
    a = where(radar.volume.h.field_type eq 'RC',c)
    radar.volume[a[0]] = rrc_vol

    a = where(radar.volume.h.field_type eq 'MC',c)
        radar.volume[a[0]] = mc_vol

;    a = where(radar.volume.h.field_type eq 'RG',c)
;    radar.volume[a[0]] = rrg_vol

;    a = where(radar.volume.h.field_type eq 'RT',c)
;    radar.volume[a[0]] = rrt_vol

    a = where(radar.volume.h.field_type eq 'RP',c)
    radar.volume[a[0]] = rrp_vol

    a = where(radar.volume.h.field_type eq 'MP',c)
        radar.volume[a[0]] = mp_vol

    flag = 'OK'
    return,flag
    end