;========================================================================
; ncl-talk: 11 Dec 2019
; May I ask if anybody has a script to interpolate GFS (0.5 deg) 
; on the desired station points and also the desired vertical level? 
; I mean (both horizontal and vertical interpolation).
;========================================================================


;; Variable: tmp_gfs
;; Type: float
;; Dimensions and sizes:
;;    [initial_time0_hours | 1] x [lv_ISBL0 | 31] x [lat_0 | 181] x [lon_0 | 360]
;; Coordinates: 
;;    initial_time0_hours: [1915164..1915164]
;;    lv_ISBL0: [100..100000]                 <== Pa
;;    lat_0: [-90..90]                        <== S->N
;;    lon_0: [ 0..359]
;; Number Of Attributes: 11
;;    center :	US National Weather Service - NCEP (WMC)
;;    production_status :	Operational products
;;    long_name :	Temperature
;;    units :	K
;;    _FillValue :	1e+20
;;    grid_type :	Latitude/longitude
;;    parameter_discipline_and_category :	Meteorological products, Temperature
;;    parameter_template_discipline_category_number :	( 0, 0, 0, 0 )
;;    level_type :	Isobaric surface (Pa)
;;    forecast_time :	12
;;    forecast_time_units :	hours

;---Force time dimension [convenience]
  setfileoption("grb","SingleElementDimensions","Initial_time") ; initial_time0_hours

;---Specify several station locations
  lat_sta = (/ 30.37, 33.59, 36.4883/)
  lon_sta = (/ 48.23, 46.4 , 61.07  /)

;---Specify pressure levels to which data is to be interpolated 
  po   = (/850,700,600,500,400,350,300,275,250,225,200,150,100/)*100.0
  po@units = "Pa"

;---Open GFS file
  diri = "/Users/shea/Data/GRIB/"
  fili = "gfs.t12z.pgrb2.1p00.f012.grib2"
  pthi = diri+fili
  f    = addfile(pthi,"r")
  
;---Import variable; ::-1 means reorder upon input
  tmp_gfs  = f->TMP_P0_L100_GLL0(:,:,::-1,:) 
  printVarSummary(tmp_gfs)  ; see above 
  print("-----")

;---Print pressure levels [Pa]
; print(tmp_gfs&lv_ISBL0)   ; hPa
;;  100 200 300 500 700 1000 2000 3000 5000 7000 10000 15000 20000 25000 30000
;;35000 40000 45000 50000 55000 60000 65000 70000 75000 80000 85000 90000 
;;92500 95000 97500 100000

  lat_gfs = tmp_gfs&lat_0     ; clarity & convenience
  lon_gfs = tmp_gfs&lon_0
  prs_gfs = tmp_gfs&lv_ISBL0

;---Vertical interpolation at original grid resolution
  linlog  =  1                ; linear   ... -1 log [arbitrary]
  tmp_gfs_p0 = int2p_n_Wrap(prs_gfs,tmp_gfs,po,linlog,1)
  printVarSummary(tmp_gfs_p0)  ;  (1,13,181,360)
  print("-----")

;---Horizontal interpolation to station locations
  tmp_sta = linint2_points_Wrap(lon_gfs,lat_gfs,tmp_gfs_p0, True , lon_sta, lat_sta, 0)
  printVarSummary(tmp_sta) ; [initial_time0_hours | 1] x [LV_isbl0 | 13] x [pts | 3]
  print("-----")

  nt = 0
  print(sprintf("%7.0f", po)             +"  "+ \
        sprintf("%7.2f", tmp_sta(nt,:,0))+"  "+ \
        sprintf("%7.2f", tmp_sta(nt,:,1))+"  "+ \
        sprintf("%7.2f", tmp_sta(nt,:,2)))