;**************************************************
; skewt_4.ncl
;
; Concepts illustrated:
;   - Drawing Skew-T plots
;   - Drawing two raobs on the same Skew-T plot
;   - Customizing the background of a Skew-T plot
;**************************************************
;
; This file is loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
;
; This file still has to be loaded maually


 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"

 begin

; --- Read Data ----------------------------------------;
  diri = "./"
  fili = "sounding_201603312345Z.txt"
  pthi = diri+fili
  NLVL = numAsciiRow(pthi)
  ncol = 6
  TestData = asciiread (diri+fili , (/NLVL,ncol/), "float")
  p_obs    = TestData (:,0)
  tc_obs   = TestData (:,1)
  tdc_obs  = TestData (:,2)
  z_obs    = TestData (:,3)
  ;wspd     = TestData (:,4)
  ;wdir     = TestData (:,5)
  u_obs = TestData (:,4)
  v_obs = TestData (:,5)

  Fill = -32768.0              ; _FillValue must be the same type as variable
  p_obs@_FillValue = Fill
  tc_obs@_FillValue = Fill
  tdc_obs@_FillValue = Fill
  z_obs@_FillValue = Fill
  ;wdir@_FillValue = Fill
  ;wspd@_FillValue = Fill
  u_obs@_FillValue = Fill
  v_obs@_FillValue = Fill
                               ; Eliminate all "p_obs" where 'p' is missing
                               ; use overwrite syntax  [  :=   ]
  ip = ind(.not.ismissing(p_obs))
  p_obs    := p_obs(ip)
  tc_obs   := tc_obs(ip)
  tdc_obs  := tdc_obs(ip)
  z_obs    := z_obs(ip)
  u_obs    := u_obs(ip)
  v_obs    := v_obs(ip)

;*************************
; create plot
;*************************
  wks  = gsn_open_wks ("png", "skewt_sound1")  ; send graphics to PNG file

; --- Create background skew-T and plot sounding----------------

  skewtOpts                 = True
  skewtOpts@DrawColAreaFill = True    ; default is False
  skewtOpts@DrawStandardAtm   = True    ; draw standard atm on plot
  skewtOpts@DrawHeightScale  = True    ; default is False
  skewtOpts@DrawHeightScaleFt = False ;plot height scale in km
  skewtOpts@DrawFahrenheit   = False   ; default is True
  skewtOpts@PrintOpts         = False   ;
  skewtOpts@tiMainString    = "CAPRICORN 2016 sounding vs. WRF sounding 31 March 2016"

  skewt_bkgd = skewT_BackGround (wks, skewtOpts)
  draw (skewt_bkgd)

  dataOpts                  = True     ; options describing data and ploting
  dataOpts@PrintZ           = True    ; do not print Z
  dataOpts@DrawFahrenheit   = False    ; default is True
  dataOpts@ThermoInfo       = False   ; print thermodynamic info

;*****************************************************************
; First sounding
;*****************************************************************
  dataOpts@colTemperature   = "blue"
  dataOpts@colDewPt         = "blue"
  dataOpts@Parcel   = 1
  dataOpts@WspdWdir  = False  ; wind speed and dir [else: u,v]
  dataOpts@HspdHdir  = True   ; wind speed and dir [else: u,v]
  dataOpts@PlotWindH = False  ; plot wind barbs at h lvls [pibal; special]
 ;dataOpts@xpWind           = 40     ; move to left [default is 45]
  dataOpts@DrawWindBarbThk = 1.0     ; wind barb thickness
  dataOpts@Wthin           = 50
  skewt_data = skewT_PlotData(wks, skewt_bkgd, p_obs, (tc_obs)-273.16, (tdc_obs)-273.16, z_obs/1000., u_obs, v_obs, dataOpts)
  draw (skewt_data)

frame(wks)

 end