;**************************************************************
; [1] Read Grib file
; [2] Using Nearest neighbor function "triple2grid",
;     interpolate to a 2x2 grid
; [3] Specify required lat/lon start/end points for a cross section
; [4] Plot cross section

;**************************************************************
; Only one variable is done here. However, this could readily
; be expanded to multiple variables.
;**************************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
  diri   = "./"                                      ; input directory
  fili   = "merged_AWIP32.1979010100.3D.NARR.grb"    ; input GRIB

;*******************************************
; open file and read in data
;*******************************************
  f       = addfile (diri+fili, "r")
  lat2d   = f->gridlat_221                 ; (277,349)
  lon2d   = f->gridlon_221
  lev     = f->lv_ISBL0                    ; (29)

  print("lat2d: min="+min(lat2d)+"   max="+max(lat2d))
  print("lev:   min="+min(lev)  +"   max="+max(lev))
                                          
  x       = f->U_GRD_221_ISBL              ; (29,277,349)    
                                           ; (lv_ISBL0, gridx_221, gridy_221)
  printVarSummary(x)
  printMinMax(x, True)     

  dimx    = dimsizes( x )
  klev    = dimx(0)
  nlat    = dimx(1)
  mlon    = dimx(2)

;********************************************
; Before interpolating make all lon >= 0.0
; The reason is to get around the dateline longitude sign switch
; going from positive [eg 150E] to negative [eg -170]
;********************************************
; use the "where" function when version a035 is available
;********************************************

;;lon2d   = where(lon2d.ge.0, lon2d, lon2d+360.)   ; a035

  LON1D   = ndtooned(lon2d)
  ii      = ind(LON1D.lt.0.)
  LON1D(ii) = LON1D(ii) + 360.
  lon2d   = onedtond( LON1D, dimsizes(lon2d) )

  print(" min="+min(lat2d)+"   max="+max(lat2d))   ; min=0.89728   max=85.3335
  print(" min="+min(lon2d)+"   max="+max(lon2d))   ; min=148.64    max=357.433

;********************************************
; specify new variable grid properties
; that original data will be regridded to
;********************************************
  NLAT          = 44                                 ; 2x2
  MLON          = 105

  dlat          = 2.
  lat           = ispan ( 0,NLAT-1,1 )*dlat          ; 0 to 86
  lat!0         = "lat"
  lat@units     = "degrees_north"
  lat&lat       = lat

  dlon          = 2.
  lon           = ispan ( 0,MLON-1,1 )*dlon + 148.   ; 148 to 358
  lon!0         = "lon"
  lon@units     = "degrees_east"
  lon&lon       = lon

  lev!0         = "lev"
;                           0    1    2
  xNEW           = new ((/klev,NLAT,MLON/), typeof(x),getFillValue(x))
  xNEW!0         = "lev"
  xNEW&lev       =  lev  
  xNEW!1         = "lat"
  xNEW&lat       =  lat
  xNEW!2         = "lon"
  xNEW&lon       =  lon
  xNEW@units     = x@units
  xNEW@long_name = x@long_name
  xNEW@time      = x@initial_time

;********************************************
; interpolate: triple2grid was designed for random
; data to a grid. It does not take advantage of any
; grid structure for efficiency. It can be slow.
;********************************************

  do kl=0,klev-1  
     print("kl="+kl+"  lev="+lev(kl))
     xNEW(kl,:,:) = triple2grid(ndtooned(lon2d),ndtooned(lat2d),ndtooned( x(kl,:,:) ) \
                               ,lon,lat, False)
  end do
  copy_VarAtts(x, xNEW)
  printVarSummary(xNEW)
  printMinMax(xNEW, True)     

;********************************************
; Create a Cross Section using a Great Circle path
; between 2 user specified points
;********************************************
; Note: user can specify any sequence of (lat,lon) pairs
; http://www.ncl.ucar.edu/Document/Functions/Built-in/linint2_points.shtml
; A great circle path is done for convenience.
;********************************************
  lat_a  =   50.0
  lon_a  =  260.0
  lat_b  =   30.0
  lon_b  =  280.0

  npts   =    50                 ; number of points in resulting transect
  dist   = gc_latlon(lat_a,lon_a,lat_b,lon_b,npts,2)
  points = ispan(0,npts-1,1)*1.0 
  print(dist@gclat+"  "+dist@gclon)

;********************************
; interpolate data to great circle
;********************************

  trans   = linint2_points(lon,lat,xNEW,False,dist@gclon,dist@gclat,2)
  copy_VarAtts(xNEW,trans)          ; copy attributes
  trans!0      = "lev"           ; create named dimension and assign
  trans&lev    =  lev            ; coordinate variable for 0th dimension only
  printVarSummary(trans)

;********************************************
; create plots
;********************************************
  wks = gsn_open_wks("ps","narr_4")    
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; choose colormap

  res                     = True         ; plot mods desired
  res@tmXBMode            = "Explicit"   ; explicitly label x-axis
  res@tmXBValues          = (/points(0),points(npts-1)/) ; points to label
                                         ; label values
  res@tmXBLabels          = (/lat_a +","+lon_a,lat_b+", "+lon_b/)

  res@cnFillOn            = True         ; turn on color
  res@lbLabelAutoStride   = True         ; nice label bar label stride
 ;res@lbBoxLinesOn        = False
  res@gsnSpreadColors     = True         ; use full range of colormap
  res@cnLinesOn           = False        ; turn off countour lines
  res@pmLabelBarHeightF   = 0.1 
 
  res@tiMainString        = "Transect"   ; add title
  plot = gsn_csm_pres_hgt(wks,trans,res) ; just for fun

;********************************
; show transect on a map
;********************************
  i                = NhlNewColor(wks,0.8,0.8,0.8) ;add gray to continents
  mres             = True                ; plot mods desired
  mres@gsnFrame    = False               ; don't turn page yet
  mres@mpMaxLatF   =  60.        
  mres@mpMaxLonF   =  300
  mres@mpMinLatF   =  20.
  mres@mpMinLonF   =  230.
  mres@tiMainString= "Transect Location" ; title
  map = gsn_csm_map_ce(wks,mres)         ; draw map
                                         ; add polyline  
  pres                  = True           ; polyline mods desired
  pres@gsLineColor      = "red"          ; color of lines
  pres@gsLineThicknessF =  2.0           ; line thickness   
  gsn_polyline(wks,map,(/lon_a,lon_b/),(/lat_a,lat_b/),pres)
  frame(wks)

end