;***************************************************************************
; read RHB surface met 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 RHB Sfc Met data into arrays and plot a variable.
; See http://www.esrl.noaa.gov/psd/enso/rapid_response/data_pub/
;*************************************************************************** 
; written by CAS 8/31/2016. V1.0
;*************************************************************************** 

begin

; --- Read Data ----------------------------------------;

  filein="psd-sfcmet_ronaldhbrown_ship_20160216.na"
  titledate = filein
  cnlvl=systemfunc("/usr/local/bin/wc "+ filein)
  nlvl=stringtointeger(cnlvl)
  head1 = readAsciiHead(filein,1)     ; there are two header lines
  strs = str_split(head1, " ")
  nheads = stringtointeger(strs(0))
  print(nheads+" ")

  print(nlvl+" obs+"+nheads)
  skiprows=nheads+1
  ncol = 19
  nlvl=nlvl-skiprows
  print(nlvl+" lines of obs")

  TestData = readAsciiTable (filein , (/ncol/), "double",skiprows) 
  TestData@_FillValue=9999
  TestData@missing_value=9999
  
  time   = TestData (:,0)     ; time [UTCday)
  year   = TestData (:,1)     ; years
  day    = TestData (:,2)     ; day
  hour   = TestData (:,3)     ; hour
  minute = TestData (:,4)     ; minutes
  seconds= TestData (:,5)     ; seconds
  p      = TestData (:,7)     ; air pressure     [mb] or [hPa] 
  tc     = TestData (:,8)     ; temperature  [C]    
  rh     = TestData (:,9)     ; relative humidity [%]
  wspd   = TestData (:,10)    ; wind speed   [m/s]    
  wdir   = TestData (:,11)    ; meteorological wind dir [from degrees]
  sst    = TestData (:,12)    ; SST [C]
  flux_d = TestData (:,13)    ; SW flux down [W m-2]
  lat    = TestData (:,14)    ; latitude [degree_north]
  lon    = TestData (:,14)    ; longitude [degree_north]
  plat_c = TestData (:,15)    ; platform_course [degree]
  plat_a = TestData (:,16)    ; platform_azimuth_angle [degree]
  plat_s = TestData (:,17)    ; platform_speed_wrt_ground  [m/s]
  
  
  i=0

;*************************
; create plot
;*************************

  res=True
  wks  = gsn_open_wks ("png", "rhb_sfcmet")
  res@tiMainString     = "Time(UTC day) vs tc(C) RHB"          ; add title
  plot  = gsn_csm_xy (wks,time,tc,res) ; create plot
  max_p = max(p)
  min_p = min(p)   
  print(max_p+" max p")
  print(min_p+" min p")  

 end