;***************************************************************************
; read CXI sonde data. format: NASA-Ames
; 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 CXISLAND sonde data into arrays.
; See http://www.esrl.noaa.gov/psd/enso/rapid_response/data_pub/
;*************************************************************************** 
; written by CAS 8/31/2016. V1.0
;*************************************************************************** 

 begin
  
  filein="psd-radiosonde_kiritimati_20160327112057_level2_cora.na"
  cnlvl=systemfunc("/usr/local/bin/wc "+ filein)
  nlvl=stringtointeger(cnlvl) 
  print(nlvl+" obs+header")

  head1 = readAsciiHead(filein,1)     ; there are two header lines
  strs = str_split(head1, " ")
  nheads = stringtointeger(strs(0))
  print(nlvl+" obs+"+nheads)
  skiprows=nheads+1

  nlvl=nlvl-skiprows
  print(nlvl+" lines of obs")

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

  TestData = readAsciiTable (filein , (/ncol/), "float",skiprows) 

; order: Surface is 'bottom'  eg: 1000,950,935,897,...

  TestData@_FillValue=-9999.
  TestData@missing_value=-9999.

  tt   = TestData (:,0)    ; time since launch (sec)
  z    = TestData (:,1)    ; geopotential [gpm]   
  p    = TestData (:,2)    ; pressure     [mb] or [hPa]
  tc   = TestData (:,3)    ; temperature  [C]    
  tdc  = TestData (:,4)    ; dew pt temp  [C]
  rh   = TestData (:,5)    ; rh [%]
  wspd = TestData (:,6)    ; wind speed   [knots or m/s]    
  wdir = TestData (:,7)    ; meteorological wind dir
  u    = TestData (:,8)    ; u (m/s)
  v    = TestData (:,9)    ; v (m/s)
  xlat = TestData (:,10)   ; Lat
  xlon = TestData (:,11)   ; Lon
  gps  = TestData (:,12)   ; GPS altitude (m)
  upsp = TestData (:,13)   ; upward platform speed w/respect to groumd
  fz   = TestData (:,14)   ; z status_flag 
  fp   = TestData (:,15)   ; air_pressure status_flag
  ftc  = TestData (:,16)   ; air_temperature status_flag
  ftdc = TestData (:,17)   ; dewpoint status_flag
  fw   = TestData (:,18)   ; wind status_flag
  fgps = TestData (:,19)   ; gps status_flag  
  
  p@missing_value=9999.
  p@_FillValue=9999.

  tc@missing_value=999.
  tc@_FillValue=999.

  tdc@missing_value=999.
  tdc@_FillValue=999.

  wdir@missing_value=999.
  wdir@_FillValue=999.

  wspd@missing_value=999.
  wspd@_FillValue=999.

  z@missing_value=99999.
  z@_FillValue=99999.
  
  wspd=wspd*1.9438    

   print("read this data:")
   print("Level 0: "+p(0)+"mb "+tc(0)+"C "+tdc(0)+"C "+z(0)+"m "+wspd(0)+"kts "+wdir(0)+"deg")

 end