;Ku-band precipitation radar (KuPR). - sensitive to strong rain
;13.6 GHz
;Swath Width=245 kilometers (km)
;Range Resolution=250 meters (m)
;NS=normal scan swath=49=245km

;  Range resolution 5km horizontal, 125m vertical

;****************************************
;  INPUTS
;  filename='/uufs/chpc.utah.edu/common/home/mace-group5/gpm/GPM_L2/GPM_2AKu.06/2A.GPM.Ku.V8-20180723.20190125-S072422-E085656.027888.V06A.SUB.HDF5'
;  sub_pt=24

pro read_gpm_2aku,filename,sub_pt,forbit,$
  lat_sub,lon_sub,colon_sub,julian_day,dbz,lat,lon,colon,nbin,$
  dbz_measured

;  GPM orbit number from filename
parts=strsplit(file_basename(filename),'.',/extract)
forbit=parts[5]

;  First part of the filename string
part1=file_dirname(filename,/mark_directory)
filesd=filename

;  Pick up another one for this orbit, needed more points
if forbit eq '028102' or forbit eq '027964' then begin 
  forbit0=forbit-1
  file_string=part1+'*.'+string(forbit0,format='(I06)')+'*.HDF5'
  filesd0=file_search(file_string,count=num_files0)
  if num_files0 gt 0 then begin
    filesd=[filesd0,filesd]
  endif
endif 

if forbit eq '028162' or forbit eq '028193' or $
   forbit eq '022353' then begin
  forbit2=forbit+1
  file_string=part1+'*.'+string(forbit2,format='(I06)')+'*.HDF5'
  filesd2=file_search(file_string,count=num_files2)
  if num_files2 gt 0 then begin
    filesd=[filesd,filesd2]
  endif  
endif

num_files=n_elements(filesd)

;******************
;  Read data from hdf5 file
;******************
colon=!null
;presult=h5_parse(filename)
for i=0,num_files-1 do begin
  ;  Open an existing hdf5 file
  print,filesd[i]
  file_id=h5f_open(filesd[i])

  ;  Group NS
  ns_id=h5g_open(file_id,'NS')
  dset_id=h5d_open(ns_id,'Latitude')
  lat1=h5d_read(dset_id)
  h5d_close,dset_id
  dset_id=h5d_open(ns_id,'Longitude')
  lon1=h5d_read(dset_id)
  h5d_close,dset_id

  ;  Group NS/SLV
  slv_id=h5g_open(ns_id,'SLV')
  dset_id=h5d_open(slv_id,'zFactorCorrected')
  dbz1=h5d_read(dset_id)
  h5d_close,dset_id
  h5g_close,slv_id

  ;  Group NS/ScanTime
  time_id=h5g_open(ns_id,'ScanTime')
  dset_id=h5d_open(time_id,'DayOfYear')
  doy=h5d_read(dset_id)
  h5d_close,dset_id
  dset_id=h5d_open(time_id,'Year')
  year=h5d_read(dset_id)
  h5d_close,dset_id
  dset_id=h5d_open(time_id,'Month')
  month=h5d_read(dset_id)
  h5d_close,dset_id
  dset_id=h5d_open(time_id,'DayOfMonth')
  day=h5d_read(dset_id)
  h5d_close,dset_id
  dset_id=h5d_open(time_id,'Hour')
  hour=h5d_read(dset_id)
  h5d_close,dset_id
  dset_id=h5d_open(time_id,'Minute')
  minute=h5d_read(dset_id)
  h5d_close,dset_id
  dset_id=h5d_open(time_id,'Second')
  second=h5d_read(dset_id)
  h5d_close,dset_id
  dset_id=h5d_open(time_id,'MilliSecond')
  msecond=h5d_read(dset_id)
  h5d_close,dset_id
  h5g_close,time_id

  ;  Group NS/PRE
  pre_id=h5g_open(ns_id,'PRE')
  dset_id=h5d_open(pre_id,'elevation')
  elevation1=h5d_read(dset_id)  ;meters
  dset_id=h5d_open(pre_id,'zFactorMeasured')
  dbz_measured1=h5d_read(dset_id)
  h5d_close,dset_id
  h5g_close,pre_id

  h5g_close,ns_id
  h5f_close,file_id
 
  ;****************
  ;  Transpose data to be nscan, npix, nheight
  ;****************
  
  ;  transpose for concatanating
  lat1=transpose(lat1)
  lon1=transpose(lon1)
  dbz1=transpose(dbz1)
  dbz_measured1=transpose(dbz_measured1)
  s=size(dbz1,/dimensions)
  ;print,s
  nbin=s[2]  ;176
  nray=s[1]  ;49
  ;nscan=s[0] ;7937

  ;  Non transpose way
  ;s=size(dbz1,/dimensions)
  ;print,s
  ;nbin=s[0]  ;176
  ;nray=s[1]  ;49
  ;nscan=s[2] ;7937

  ;*****************
  ;  Caculate colongitude because the pacific crosses the 180 line
  ;*****************
  colon1=lon1
  result=where(lon1 lt 0 and lon1 ne -9999.90,count)
  if count gt 0 then colon1[result]=colon1[result]+360.0

  ;*****************
  ;  Pull out subsatellite point
  ;*****************

  ;  transpose
  lat_sub1=reform(lat1[*,sub_pt])
  lon_sub1=reform(lon1[*,sub_pt])
  colon_sub1=reform(colon1[*,sub_pt])
  
  ;  non transpose
  ;lat_sub1=reform(lat1[sub_pt,*])
  ;lon_sub1=reform(lon1[sub_pt,*])
  ;colon_sub1=reform(colon1[sub_pt,*])

  ;*****************
  ;  Calculate time
  ;*****************
  msecond=msecond*1e-3
  second=second+msecond
  julian_day1=julday(month,day,year,hour,minute,second)

  if colon eq !null then begin
    colon=colon1
    lon=lon1
    lat=lat1
    dbz=dbz1
    dbz_measured=dbz_measured1
    julian_day=julian_day1
    lat_sub=lat_sub1
    lon_sub=lon_sub1
    colon_sub=colon_sub1   
  endif else begin
    print,'appending second file'
    colon=[colon,colon1]
    lon=[lon,lon1]
    lat=[lat,lat1]
    dbz=[dbz,dbz1]
    dbz_measured=[dbz_measured,dbz_measured1]
    julian_day=[julian_day,julian_day1]
    lat_sub=[lat_sub,lat_sub1]
    lon_sub=[lon_sub,lon_sub1]
    colon_sub=[colon_sub,colon_sub1]   
  endelse
  
endfor

lat=transpose(lat)
lon=transpose(lon)
colon=transpose(colon)
dbz=transpose(dbz)
dbz_measured=transpose(dbz_measured)

end