;***************************************************************************
; read GlobalHawk dropsonde data. format: FRD format
; PSD provides this code to make it easier to read the data. We do
; not support NCL and advise those with questions to read NCL webpages
; and/or contact their support staff.
;
; This code reads the Global Hawk sonde data into arrays.
; See http://www.esrl.noaa.gov/psd/enso/rapid_response/data_pub/
;*************************************************************************** 
; written by CAS 9/29/2016. V1.0
;*************************************************************************** 

 begin
  
  filein="D20160215_214000_P.QC.eol"
  cnlvl=systemfunc("/usr/local/bin/wc "+ filein)
  nlvl=stringtointeger(cnlvl) 
  skiprows=14
  nlvl=nlvl-skiprows
  print(nlvl+" lines of observations from "+filein)

  ncol = 17
  x = ispan (1,nlvl,1)

  TestData = readAsciiTable (filein , (/ncol/), "float",skiprows) 
  TestData@_FillValue=-999.
  TestData@missing_value=-999.

  tt   = TestData (:,0)    ; sensor observation
  ohr  = TestData (:,1)    ; hr obs
  omin = TestData (:,2)    ; min obs
  osec = TestData (:,3)    ; sec obs
  p    = TestData (:,4)    ; pressure     [mb / hPa]
  tc   = TestData (:,5)    ; temperature  [C]
  tdc  = TestData (:,6)    ; dew point temperature  [C]      
  rh   = TestData (:,7)    ; dew pt temp  [C]
  u    = TestData (:,8)    ; u [m/s]
  v    = TestData (:,9)    ; v [m/s]
  wspd = TestData (:,10)    ; wind speed   [knots or m/s]    
  wdir = TestData (:,11)    ; meteorological wind dir  
  dz   = TestData (:,12)    ; delta geopotential [m/s]
  z    = TestData (:,13)    ; geopotential [m]
  xlon = TestData (:,14)   ; Lon [E]
  xlat = TestData (:,15)   ; Lat [N]
  gps  = TestData (:,16)   ; GPSAlt [m]


  wspd=wspd*1.9438    

   print("Read the data:")
   print("First line of observations: "+p(0)+"mb "+tc(0)+"C "+tdc(0)+"C "+z(0)+"m "+wspd(0)+"kts "+wdir(0)+"deg")
   
iz = ind(.not.ismissing(p))  ; iz will have indices of non-missing values
                               ; iz(0)=0, iz(1)=2, iz(2)=4

kuse=iz(0)
   print("First non-missing pressure line: "+p(kuse)+"mb "+tc(kuse)+"C "+tdc(kuse)+"C "+z(kuse)+"m "+wspd(kuse)+"kts "+wdir(kuse)+"deg")
 end