function retrieve_hid,radar,sounding_file=sounding_file,$
                          ldr=ldr,qc=qc,POSTSCRIPT=postscript
    COMMON misc,TRUE,FALSE
    COMMON graphics,USE_Z_BUFFER,TC,XS,YS,cols,rows
    COMMON graphic_parms,scale,bottom,ncolors
    COMMON dirs,img_dir,uf_dir,hid_img_dir,hid_dir,dsd_dir
    !QUIET=1

;    print,'Working in npol_radial_hid!'
    if(NOT KEYWORD_SET(POSTSCRIPT)) then postscript=FALSE
    
;This is a program to run HID on polar-coordinate UF files.  The program
;makes use of the RSL in IDL TRMM GV library suite.   
;Written by Brenda Dolan
;10 July 2008
;Modifications by Timothy Lang, April 2011
;Modifications by Brenda Dolan, June 2012 to read NASA_QC NPOL files,
;as 
;well as the new HID beta functions (cdf_betas_sum.pro) from
;scattering 
;simulations, and the new fuzzy logic algorithm
;		cdf_fhc_all.pro
;
; INPUTS
;		  file: Specify a specific file.  If not included, it
;		  will 
;look for all *.uf files sounding_file: Specify the sounding file to
;be used in the temperature interpolation.  If nothing is specified, 
; it will use a default melting level and standard lapse rate.
;
; Sounding file is assumed to come from the University of Wyoming
; site, where the 2nd column is height in m (MSL) and the 3rd column
; is temperature in Cº.
;		   ldr: set this keyword if LDR is one of the
;		        available variables to be used in the HID.
;		    qc: set this keyword to use a flag field or the 
;                       sq, VR and KD fields to threshold bad data.
; OUTPUT
;	FH: Added to the UF file. Corresponds to the fuzzy logic HID field.
;
; *** Set to postscript
;
    USE_Z_BUFFER=TRUE
    bad=-32767.0
;
; *** Check to see if a structure was returned.
;
    a = size(radar,/type)
    if(a ne 8) then begin
        flag = 'No structure returned: ' + file
        return,flag
    endif
    date_time=rsl_datetime(radar,jdate=jdate)
    
;;; Change these varaible names to correspond to the names in the 
;;; UF file that is being read.
    
    dzname='CZ'
    drname='DR'
    kdname='KD'
    rhname='RH'
    vename='VR'
    rrname='RR'
    swname='SW'
    sqname='SQ'
    
;;; This will change the plot range.
    
    maxr=125
    minr=0
    
    n_elev=n_elements(radar.volume[0].sweep[0].ray.h.elev)
    n_bins=n_elements(radar.volume[0].sweep[0].ray[0].range)
    n_sweep=radar.volume[0].h.nsweeps
    gate_size=float(radar.volume[0].sweep[0].ray[0].h.gate_size)
    
    wh_zh=where(radar.volume.h.field_type eq dzname)
    DZ=radar.volume[wh_zh].sweep.ray.range
    bad=radar.volume[wh_zh].h.no_data_flag
    
;;;;;;;NEW;;;;;:Add SW and SQ to threshold based on these;;;;;;;;;
    
    wh_sw=where(radar.volume.h.field_type eq swname)
    IF wh_sw[0] ne -1 THEN BEGIN
        SW=radar.volume[wh_sw].sweep.ray.range
    ENDIF ELSE BEGIN
        SW=DZ*0.0+bad
    ENDELSE
    
    wh_sq=where(radar.volume.h.field_type eq sqname)
    IF wh_sq[0] ne -1 THEN BEGIN
        SQ=radar.volume[wh_sq].sweep.ray.range
    ENDIF ELSE BEGIN
        SQ=DZ*0.0+bad
    ENDELSE

;;;;;;;;;;;;;END NEW;;;;;;;;;;;;
    
    
    wh_zd=where(radar.volume.h.field_type eq drname)
    IF wh_zd[0] ne -1 THEN BEGIN
        DR=radar.volume[wh_zd].sweep.ray.range
    ENDIF ELSE BEGIN
        DR=DZ*0.0+bad
    ENDELSE
    
    wh_kd=where(radar.volume.h.field_type eq kdname)
    IF wh_kd[0] ne -1 THEN BEGIN
        KD=radar.volume[wh_kd].sweep.ray.range
    ENDIF ELSE BEGIN
        KD=DZ*0.0+bad
    ENDELSE
    
    wh_rh=where(radar.volume.h.field_type eq rhname)
    IF wh_rh[0] ne -1 THEN BEGIN
        RH=radar.volume[wh_rh].sweep.ray.range
    ENDIF ELSE BEGIN
        RH=DZ*0.0+bad
    ENDELSE
    
    wh_ve=where(radar.volume.h.field_type eq vename)
    IF wh_ve[0] ne -1 THEN BEGIN
        VE=radar.volume[wh_ve].sweep.ray.range
    ENDIF ELSE BEGIN
        VE=DZ*0.0+bad
    ENDELSE
    
    wh_rr=where(radar.volume.h.field_type eq rrname)
    IF wh_rr[0] ne -1 THEN BEGIN
        RR=radar.volume[wh_rr].sweep.ray.range
    ENDIF ELSE BEGIN
        RR=DZ*0.0+bad
    ENDELSE
    
    IF keyword_set(ldr) THEN BEGIN
        wh_lh=where(radar.volume.h.field_type eq 'LH')
        LD=radar.volume[wh_lh].sweep.ray.range
        fhc_vars={Zh:1,Zdr:1,Kdp:1,LDR:1,rho_hv:1,T:1}
        fhc_vars.LDR=1
    ENDIF ELSE BEGIN
        fhc_vars={Zh:1,Zdr:1,Kdp:1,LDR:0,rho_hv:1,T:1}
        fhc_vars.LDR=0        
    ENDELSE
    
;;;;NOW we need to clean up the data based on the Flag field from DROPS
    
    IF keyword_set(qc) THEN BEGIN
        wh_fl=where(radar.volume.h.field_type eq 'FL')
        IF wh_fl[0] ne -1 THEN BEGIN
            FL=radar.volume[wh_fl].sweep.ray.range
            thresh=where(fl eq 0)
        ENDIF ELSE BEGIN

;;;;;IF FL is not available, use a very rudimentary method of cleaning up the data
;;;;;;REPLACE: THRESHOLDING BASED ON KD, SQ, and SW/VE instead of Dz, ZDr;;;;;;
            
            thresh=where(SQ lt 0.2 OR kd eq bad or $
                         (sw lt 1.0 and ve lt 0.1 and ve gt -0.1 ))
;;;;;;ENDREPLACEMENT;;;;;;
            
;            print,'Using SQ / KD for thresholding!'
        ENDELSE
        
        dz[thresh]=bad
        dr[thresh]=bad
        kd[thresh]=bad
        rh[thresh]=bad
        ve[thresh]=bad
        rr[thresh]=bad
; write the masked dr and dz back into their structures...
        
        radar.volume[wh_zh].sweep.ray.range=DZ
        radar.volume[wh_zd].sweep.ray.range=DR
        radar.volume[wh_kd].sweep.ray.range=KD
        radar.volume[wh_rh].sweep.ray.range=RH
        radar.volume[wh_ve].sweep.ray.range=VE
        radar.volume[wh_rh].sweep.ray.range=RH
        IF wh_rr[0] ne -1 THEN BEGIN
            radar.volume[wh_rr].sweep.ray.range=RR
        ENDIF
    ENDIF
    nelev=n_elements(DZ(0,*,0))
    nbins=n_elements(DZ(*,0,0))
    nsweep=n_elements(DZ(0,0,*))
    
;;;;Calculate heights for each bin;;;;;;;;;;
    
    height=fltarr(nbins,nelev,nsweep)
    range=fltarr(nbins,nelev,nsweep)
    elev=fltarr(nbins,nelev,nsweep)
    relev=fltarr(nbins,nelev,nsweep)
    T_dat=fltarr(nbins,nelev,nsweep)
;range=findgen(n_bins)*150.
    
    FOR ii=0,nbins-1 DO BEGIN
        FOR jj=0,nelev-1 DO BEGIN
            FOR kk=0,nsweep-1 DO BEGIN
                range[ii,jj,kk]=float(ii)*gate_size/1000.
                elev1=float(radar.volume[0].sweep[kk].ray[jj].h.elev)
                IF elev1 ne bad THEN BEGIN
                    elev(ii,jj,kk)=elev1
                    relev(ii,jj,kk)=elev(ii,jj,kk)*3.14159/180.
                    IF elev(ii,jj,kk) lt 0.0 THEN BEGIN
                        elev(ii,jj,kk)=0.0
                        relev(ii,jj,kk)=0.0
                    ENDIF
        ;;;Assume a standard relationship betwene height, range and elevation angle.
                    height[ii,jj,kk]=(range[ii,jj,kk])^2/17000.+$
                      (range[ii,jj,kk])*SIN(relev[ii,jj,kk])
                    
                    IF height(ii,jj,kk) lt 0.0 THEN BEGIN
;                        print,'Uh-oh! sounding height is less than zero!'
                        break
                    ENDIF
                    
                ENDIF ELSE BEGIN
                    elev(ii,jj,kk)=0.0
                    relev(ii,jj,kk)=0.0
                    height[ii,jj,kk]=0.0
                ENDELSE
            ENDFOR
        ENDFOR
    ENDFOR
    
    
;aglheight=height;+0.5
    
    
;;;;This loop requires a univeristy of wyoming sounding, where height 
;;;;is in the 2nd column and temperature in the 3rd. Really,
;;;;temperature can come from any source for the HID it just
;;;;needs to end up the same size as the other HID input variables.
    
    line=''
    IF keyword_set(sounding_file) THEN BEGIN
;        print,'Using sounding!'
        openr,unit,sounding_file,/get_lun
        sdumh=fltarr(200)
        sdumt=fltarr(200)
        i=0
;
; *** Header lines
;
        readf,unit,line
        readf,unit,line
        WHILE ~EOF(unit) DO BEGIN
            readf,unit,line
            line_split=strsplit(line,'  ', /Extract)
            
            IF n_elements(line_split) gt 2 THEN BEGIN
                sdumh(i)=line_split[1]
                sdumt(i)=line_split[2]
            ENDIF ELSE BEGIN
                sdumh(i)=0.00
                sdumt(i)=0.00
            ENDELSE
            i=i+1
        ENDWHILE
        close, unit
        free_lun,unit
        
        use_temp=1
    ENDIF ELSE BEGIN
        
        use_temp=1
        print,'No sounding file, using typical lapse rate 6 deg / km'
        stop
;use_temp=0
        to=20.0
        zo=0.5
        sdumh=findgen(20)*1000+zo
        sdumt=findgen(20)*(-6)+to
    ENDELSE
    aglheight=height+0.3
    
;Now we need to interpolate the temperatures to our grid
    wh_dat=where(sdumh ne 0.00)
    sheight=sdumh(wh_dat)/1000.
    temp=sdumt(wh_dat)
;T_dat=interp_sounding(aglheight,sheight,temp)
    
;    print,'Interpolating sounding!'
    FOR ll=0,n_bins-1 do begin
        FOR j=0,n_elev-1 do begin
            FOR k=0,n_sweep-1 do begin
                diff=abs(aglheight(ll,j,k)-sheight)
                wh_m=where(diff eq min(diff)) ;
                T_dat(ll,j,k)=temp(wh_m[0])
            ENDFOR
        ENDFOR
    ENDFOR
;    print,'Done interpolating sounding!'
    use_temp=1
    
    NoMix=1            ;We no longer use mix types such as rain + hail
    
; Get which variables to use in the FHC

    fhc_vars.Zh=1
    fhc_vars.Zdr=1
    fhc_vars.Kdp=1
    fhc_vars.rho_hv=1
    fhc_vars.T=1
    var_names = ''
    
; Find which ones are going to be used.
    
    IF(fhc_vars.Zh)  THEN var_names += 'Zh'
    IF(fhc_vars.Zdr) THEN var_names += ', Zdr'
    IF(fhc_vars.Kdp) THEN var_names += ', Kdp'
    IF(fhc_vars.LDR) THEN var_names += ', LDR'
    IF(fhc_vars.rho_hv) THEN var_names += ', rho_hv'
    IF((fhc_vars.T)AND(use_temp)) THEN var_names += ', T'
    var_names = '('+var_names+')'
    
;        radar_data = {DZ:DZ,DR:DR,KD:KD,RH:RH,LD:LD}
    IF keyword_set(ldr) THEN BEGIN
        radar_data = {DZ:DZ,DR:DR,KD:KD,RH:RH,LD:LD}
    ENDIF ELSE BEGIN
        radar_data = {DZ:DZ,DR:DR,KD:KD,RH:RH}
    ENDELSE
    
;;; Get the membership beta functions.  The cdf_betas_sum is 
;;; the summer hid from simulations, as of June 2012. Specify 
;;; the /use_temp keyword to include temperature.
    
    cdf_betas_sum,mbf_sets,use_temp=use_temp
    
;;; Define the weights for each variable.  Although Zh and T are 
;;; given weights, a hybrid method is used in the actual algorithm 
;;; where the pol variables are weighted an the final score is
;;; multiplied by the Zh and T scores.
    
    weights= {W_Zh: 1.5,$       ; Same for all hydrotypes.
              W_Zdr: 0.8,$                      
              W_Kdp: 0.8,$
              W_rho: 0.4,$
              W_Ldr: 0.5,$
              W_T: 0.4}
    
;;; Call the HID subroutine, cdf_fhc_all. Pass the mbf_sets, 
;;; radar_data structure, variables to be used (fhc_vars), temperature 
;;; data, and weighting structure.
    
    FHTN = cdf_fhc_all(radar_data,mbf_sets,fhc_vars,T_dat, $
                       weights=weights,use_temp=use_temp)
    
;;Add the HID To the radar structure.
    
;    print,'DZNAME: ',dzname
    vol=rsl_get_volume(radar,dzname)

    rsl_add_volume,radar,vol,field='FH'
    wh_fht = where(radar.volume.h.field_type eq 'FH')
;    radar.volume[wh_fht].sweep.ray.range=*FHTN

    radar.volume[wh_fht].sweep.ray.range=FHTN 

;;; NEW: Specifiy hte maxrng to be used in the plot as the 
;;; maximum range.
    
    maxrng=max(range)
    maxhgt=15
    
;USE maxrng in all the rsl_plotsweep_from_radar calls in lines 412-458;
;;;;;;;;;;ENDNEW;;;;;;;;;;;;;;;
    
    dum='NPOL_'+date_time
    
;fields=[dzname,drname,kdname,rhname,'FH',vename,rrname]
;IF keyword_set(ldr) THEN BEGIN
;	fields=[dzname,drname,kdname,rhname,'FH',vename,rrname,'LD']
;ENDIF

    scan_type=radar.h.scan_mode
    
;Look for the azimuth or elevation angle.
    
    IF scan_type eq 'RHI' THEN BEGIN
;        FOR q=0,nsweep-1 DO BEGIN
        FOR q=0,0 DO BEGIN
;            plot_save_rhi,radar,drname,'ZD', q,maxhgt,maxrng,POSTSCRIPT=postscript
;            plot_save_rhi,radar,vename,'VR', q,maxhgt,maxrng,POSTSCRIPT=postscript
;            plot_save_rhi,radar,dzname,'CZ', q,maxhgt,maxrng,POSTSCRIPT=postscript
;            plot_save_rhi,radar,kdname,'KC',q,maxhgt,maxrng,POSTSCRIPT=postscript
;            plot_save_rhi,radar,'FH','FHSUM',q,maxhgt,maxrng,POSTSCRIPT=postscript
        ENDFOR
    ENDIF ELSE BEGIN
;
; *** Plot PPI sweeps
;            
        FOR ss=0,nsweep-1 DO BEGIN
;            plot_save_sweep,radar,drname,'ZD',ss,maxrng
;            plot_save_sweep,radar,vename,'VR',ss,maxrng
;            plot_save_sweep,radar,'RH'  ,'RH',ss,maxrng
;            plot_save_sweep,radar,rrname,'RR',ss,maxrng
;            plot_save_sweep,radar,kdname,'KD',ss,maxrng
;            plot_save_sweep,radar,dzname,'DZ',ss,maxrng
;            plot_save_sweep,radar,'FH'  ,'FHSUM',ss,maxrng
        ENDFOR
    ENDELSE

;;;;Write the data out to a new file with _hid.uf and a new variable 'FH'.
    
;    hid_uf_file= hid_dir + org_file_base + '.HID.uf.gz'
;    rsl_radar_to_uf_gzip,radar,hid_uf_file
;    print,' -> ' + hid_uf_file
;    print
    
;Do some cleanup so that if running a lot if files it doesn't run out of memory.
    
;    radar=0
;    *fhtn=0
    fhtn=0
    dz=0
    dr=0
    rh=0
    kd=0
    sheight=0
    temp=0
    sdumh=0
    sdumt=0
    elev=0
    relev=0
    height=0
    ve=0
    rr=0
    T_dat=0

    return,'OK'
    END