;***************************************************************************
; read CXI sonde 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 G-IV sonde data into arrays.
; See http://www.esrl.noaa.gov/psd/enso/rapid_response/data_pub/
;*************************************************************************** 
; written by CAS 9/21/2016. V1.0
;*************************************************************************** 

 begin
  
  filein="D20160121_231857_PQC.frd"
  cnlvl=systemfunc("/usr/local/bin/wc "+ filein)
  nlvl=stringtointeger(cnlvl) 
  skiprows=21
  nlvl=nlvl-skiprows
  print(nlvl+" lines of observations from "+filein)

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

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

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

  TestData@_FillValue=-999.
  TestData@missing_value=-999.

  obs  = TestData (:,0)    ; sensor observation
  tt   = TestData (:,1)    ; time since launch (sec)  
  p    = TestData (:,2)    ; pressure     [mb / hPa]
  tc   = TestData (:,3)    ; temperature  [C]    
  rh  = TestData  (:,4)    ; dew pt temp  [C]   
  z    = TestData (:,5)    ; geopotential [gpm]
  wspd = TestData (:,7)    ; wind speed   [knots or m/s]    
  wdir = TestData (:,6)    ; meteorological wind dir
  u    = TestData (:,8)    ; u [m/s]
  v    = TestData (:,9)    ; v [m/s]
  ns   = TestData (:,10)   ; ??
  wz   = TestData (:,11)   ; w [m/s]
  zw   = TestData (:,12)   ; zw [m]  
  fp   = TestData (:,13)   ; air_pressure status_flag
  ftc  = TestData (:,14)   ; air_temperature status_flag
  ftz  = TestData (:,15)   ; z status_flag
  fw   = TestData (:,16)   ; w status_flag
  xlat = TestData (:,17)   ; Lat [N]
  xlon = TestData (:,18)   ; Lon [E]
  tdc  = tc

; calculated td. note NCL's Td calculate sets the missing incorrectly.

  tdctmp = dewtemp_trh(tc+273.15,rh)
  tdctmp = tdctmp-273.15
  tdc    = where(tdctmp.lt.-9999,tdctmp,TestData@_FillValue)

  ; convert wspd to knots
  
  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