; *****************************COPYRIGHT******************************
;
; (c) CROWN COPYRIGHT 2000, METEOROLOGICAL OFFICE, All Rights Reserved.
;
; ******************************COPYRIGHT******************************
;+ 
; ----------------------------------------------------------------------------+
; FUNCTION:  READ_CHILRADAR
;
; PURPOSE:   Reads Chilboton radar data in ASCII or NetCDF format
;
; ARGUMENTS: 
;  IN: filenames - an array of filenames
; OUT: radardata - a structure holding the radar data
;
; OPTIONAL ARGUMENTS:
;  maxnumlevs=maxnumlevs:   maximum number of levels in the data
;  maxnumtimes=maxnumtimes: maximum number of times in the data
; 
; OUPUT DATA STRUCTURE
; Data structure for output:
;  year:        year
;  month:       month
;  day:         day 
;  numlevs:     number of range gates (i.e. no. of data values in a profile)
;  numtimes:    number of rays (i.e. number of profiles)
;  nodatamdi    missing data indicator for no data (e.g. radar turned off)
;  nosigmdi     missing data indicator for no signal (i.e. beyond sensitivity)
;  radarid:     radar frequency (35GHz or 94GHz)
;  latitude:    latitude of radar
;  longitude:   longitude of radar
;  altitude:    altitude of radar
;  time:        time (decimal hours since midnight GMT)
;  height:      height above mean sea level (km)
;  refval:      calibrated radar reflectivity factor (dBZ)
;  (Note that cloud-free pixels are indicated by the value -99.00) 
;
; History:
; 
; 12/04/00    Version 1.0 first released. Author Richard Forbes
; 01/11/00    Included option to read NetCDF format files. RMF
;
;-----------------------------------------------------------------------------+

FUNCTION read_chilradar,filelist,fileformat=fileformat,$
         maxnumlevs=maxnumlevs,maxnumtimes=maxnumtimes

IF n_elements(maxnumlevs)  EQ 0 THEN maxnumlevs=218
IF n_elements(maxnumtimes) EQ 0 THEN maxnumtimes=2880 ;(1 day, 30 sec interval)

num_files = N_ELEMENTS(filelist)
nodatamdi = -1.0E+30
nosigmdi  = -999.0
year      = ' '
month     = ' '
day       = ' '
numlevs   = 0
numtimes  = 0
radarid   = ' '
latitude  = 0.0
longitude = 0.0
altitude  = 0.0
hour      = STRARR(maxnumtimes)
minute    = STRARR(maxnumtimes)
; Set up time array for one day
; this assumes that the time interval between observations is 30 seconds
time      = FINDGEN(maxnumtimes)*30.0/(24*3600.0)
height    = FINDGEN(maxnumlevs)
refval    = FLTARR(maxnumtimes,maxnumlevs)
refval(*,*) = nodatamdi
  
  tmp = {chilrad $
      ,year:year $
      ,month:month $
      ,day:day $
      ,numlevs:numlevs $
      ,numtimes:numtimes $
      ,nosigmdi:nosigmdi $
      ,nodatamdi:nodatamdi $  
      ,radarid:radarid $
      ,latitude:latitude $
      ,longitude:longitude $
      ,altitude:altitude $
      ,hour:hour $
      ,minute:minute $
      ,time:time $
      ,height:height $
      ,refval:refval }

radardata = REPLICATE(tmp,num_files)

;=======================================================================
; Loop over files 
;=======================================================================
FOR filenum = 0, num_files-1 DO BEGIN

  file = filelist(filenum)

  ; if fileformat is not set, make an intelligent guess based on file suffix
  IF n_elements(fileformat) EQ 0 THEN BEGIN
    IF STRPOS(file,'.nc') NE -1 OR STRPOS(file,'.ncdf') NE -1 THEN $
       fileformat=1 ELSE fileformat=0
  ENDIF

;-----------------------------------------------------------------------
; Read in ASCII format files
;-----------------------------------------------------------------------

  IF fileformat EQ 0 THEN BEGIN

    PRINT,'Reading ASCII format data from file '+STRCOMPRESS(file,/remove_all)

    ; Read header
    status = DC_READ_FIXED(file, $
         year,month,day,numlevs,numtimes,radarid,$
         Format='(A4,A2,A2,I6,I7,A3)')

    ; Check on array sizes 
    IF numlevs GT maxnumlevs THEN $
      PRINT,'WARNING numlevs ('+STRTRIM(numlevs,2)+') GT maxnumlevs (',$
             STRTRIM(maxnumlevs,2)+')'
    IF numtimes GT maxnumtimes THEN $
      PRINT,'WARNING numtimes ('+STRTRIM(numtimes,2)+') GT maxnumtimes (',$
             STRTRIM(maxnumtimes,2)+')'
   
    ; Set up arrays
    numvals = LONG(numtimes)*LONG(numlevs)
    time1   = FLTARR(numvals)
    height1 = FLTARR(numvals)
    refval1 = FLTARR(numvals)

    ; Read all data in a file
    status = DC_READ_FIXED(file, $
             time1,height1,refval1,$
             Format='(F7.5,F6.3,F7.3)',/Column,NSkip=1)

    IF status lt 0 THEN PRINT,'Error reading file: status=',status
 
    ; Set up time variables
    time2 = REFORM(time1,numlevs,numtimes)

    hour(0:numtimes-1) = STRTRIM(FIX(time2(0,0:numtimes-1)),2)
    index = WHERE(STRLEN(hour) EQ 1)
    IF index(0) NE -1 THEN hour(index) = '0'+hour(index)

    minute(0:numtimes-1) = STRTRIM(FIX((time2(0,0:numtimes-1) MOD 1.0)*60.0),2)
    index = WHERE(STRLEN(minute) EQ 1)
    IF index(0) NE -1 THEN minute(index) = '0'+minute(index)

    ; set "time" and "height" to be 1D arrays, set refval to be a 2D array
    ; with x=time and y=height
    time(0:numtimes-1)  = time2(0,0:numtimes-1)
    height(0:numlevs-1) = height1(0:numlevs-1)
    refval(0:numtimes-1,0:numlevs-1) = $
      TRANSPOSE(REFORM(refval1,numlevs,numtimes))

  ENDIF ; on fileformt = 0 i.e. ASCII

;-----------------------------------------------------------------------
; Read NetCDF format file
;-----------------------------------------------------------------------

  IF fileformat EQ 1 THEN BEGIN

    ; start up HDF module to enable access to the NetCDF routines if 1st file
    IF filenum EQ 0 THEN BEGIN 
      @hdf_startup
    ENDIF

    PRINT,'Reading NetCDF format data from file '+STRCOMPRESS(file,/remove_all)

    ; open the file, and get an id for the file (like a unit number)
    cdfid = NCOPEN(file, NC_NOWRITE)

    ; get id's for time and range *dimensions*
    time_id = NCDIMID(cdfid, "time")
    range_id = NCDIMID(cdfid, "range")

    ; Use enquiry function to get array dimensions
    st = NCDIMINQ(cdfid,time_id, "time", time_size)
    st = NCDIMINQ(cdfid,range_id, "range", range_size)

    ; Need to predefine types (and size if arrays) of variables
    frequency = 0.0
    latitude  = 0.0
    longitude = 0.0
    altitude  = 0.0
    time1     = fltarr(time_size)
    height1   = fltarr(range_size)
    nez       = fltarr(time_size)
    ;Note!  Dimensions in PV-WAVE are in reverse order of those in CDL file!!
    refval1   = fltarr(range_size, time_size) 

    ; Get date information
    year = [0]
    st   = ncattget(cdfid, NC_GLOBAL, "year", year)
    year = STRING(year(0))
    month =[0]
    st    = ncattget(cdfid, NC_GLOBAL, "month", month)
    month = STRING(month(0))
    day  =[0]
    st   = ncattget(cdfid, NC_GLOBAL, "day", day)
    day  = STRING(day(0))

    ; Get radar frequency information
    freq_varid = ncvarid(cdfid, 'frequency')
    st = NCVARGET1(cdfid, freq_varid, 1, frequency)

    ; Get latitude and longitude information
    lat_varid = ncvarid(cdfid, 'latitude')
    st = NCVARGET1(cdfid, lat_varid, 1, latitude)
    lon_varid = ncvarid(cdfid, 'longitude')
    st = NCVARGET1(cdfid, lon_varid, 1, longitude)

    ; Get radar altitude information
    alt_varid = ncvarid(cdfid, 'altitude')
    st = NCVARGET1(cdfid, alt_varid, 1, altitude)            
      
    ; Get time data vector
    ;----------------------------
    time_varid = ncvarid(cdfid, 'time')
    st = NCVARGET(cdfid, time_varid, 0, time_size, time1)

    ; Get range data vector
    ;----------------------------
    range_varid = NCVARID(cdfid,"range")
    st = NCVARGET(cdfid, range_varid, 0, range_size, height1)

    ; Get reflectivity data array
    ;----------------------------
    zh_varid = ncvarid(cdfid, 'Zh')
    st = NCVARGET(cdfid, zh_varid, [0,0], [range_size, time_size], refval1)

    ; Close netcdf file
    st = NCCLOSE(cdfid)
    
    numlevs  = range_size
    numtimes = time_size
    radarid  = STRTRIM(FIX(frequency),2)

    ; set up time arrays    
    hour(0:numtimes-1) = STRTRIM(FIX(time1(0:numtimes-1)),2)
    index = WHERE(STRLEN(hour) EQ 1)
    IF index(0) NE -1 THEN hour(index) = '0'+hour(index)

    minute(0:numtimes-1) = STRTRIM(FIX((time1(0:numtimes-1) MOD 1.0)*60.0),2)
    index = WHERE(STRLEN(minute) EQ 1)
    IF index(0) NE -1 THEN minute(index) = '0'+minute(index)

    ; work out index for times in the file (not all times may be present) 
    timeindex = FIX((time1(0:numtimes-1)-time1(0))*3600/30.0+0.1)

    height = height1(0:numlevs-1)
    refval(timeindex,0:numlevs-1) = TRANSPOSE(refval1(0:numlevs-1,0:numtimes-1))

  ENDIF ; on fileformat = 1 i.e. NetCDF file

;-----------------------------------------------------------------------
; Put data into a structure
;-----------------------------------------------------------------------

    tmp = {chilrad $
      ,year:year $
      ,month:month $
      ,day:day $
      ,numlevs:numlevs $
      ,numtimes:numtimes $
      ,nosigmdi:nosigmdi $
      ,nodatamdi:nodatamdi $
      ,radarid:radarid $
      ,latitude:latitude $
      ,longitude:longitude $
      ,altitude:altitude $
      ,hour:hour $
      ,minute:minute $
      ,time:time $
      ,height:height $
      ,refval:refval }

    radardata(filenum) = tmp      

ENDFOR ; on filenum

RETURN, radardata

END