;************************************
; trans_1.ncl
;************************************
;
; Concepts illustrated:
;   - Calculating and plotting a transect
;   - Using gc_latlon to calculate a great circle path
;   - Using linint2_points to interpolate rectilinear grid values to set of lat/lon points
;   - Attaching polylines to a map plot
;   - Explicitly setting tickmarks and labels on the bottom X axis
;************************************
; These files are loaded by default in NCL V6.2.0 and newer
; 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
  in = addfile("gfs.t06z.pgrb2b.0p25.f012.grib2","r")
  t  = in->RH_P0_L100_GLL0(:,::-1,:)

  print(t&lv_ISBL0)

;************************************
; calculate great circle along transect
;************************************
  leftlat  =  14.0
  rightlat =  18.2

  leftlon  =  120.9
  rightlon =  121.47

  npts     =   1000                    ; number of points in resulting transect

  dist     = gc_latlon(leftlat,leftlon,rightlat,rightlon,npts,2)
  points   = ispan(0,npts-1,1)*1.0 
;********************************
; interpolate data to great circle
;********************************
  trans   = linint2_points(t&lon_0,t&lat_0,t,True,dist@gclon,dist@gclat,2)
  copy_VarAtts(t,trans)          ; copy attributes

  trans!0      = "z_t"           ; create named dimension and assign
  trans&z_t    = t&lv_ISBL0           ; coordinate variable for 0th dimension only
;********************************
; create plot
;********************************
  wks = gsn_open_wks("png","trans")       ; send graphics to PNG file

  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          = (/leftlat +", "+leftlon,rightlat+", "+rightlon/)

 res@tmYBMode            = "Explicit" 
 res@tmYBLabels  = (/875,825,775,725,675,625,575,525,475,425,375,325,275,225,175,125,70,50,30,20,10/)



  res@cnFillOn            = True         ; turn on color
  res@cnFillPalette       = "BlAqGrYeOrReVi200"  ; set color map
  res@cnLinesOn           = False        ; turn off countour lines
  res@lbOrientation       = "vertical"   ; vertical label bar
  res@pmLabelBarOrthogonalPosF = -0.05        ; move label bar closer to plot
 
  res@tiMainString        = "Transect"   ; add title
  res@tiXAxisString       = "lat/lon along transect"
  res@trYReverse          = True         ; reverse y axis
;  res@trXReverse          = True         ; reverse x axis (neg longitudes)
  res@cnLevelSpacingF     = 1.0          ; set contour spacing
  
  plot = gsn_csm_contour(wks,trans,res)  ; create plot
;********************************
; show transect on a map
;********************************
  mres              = True                ; plot mods desired
  mres@gsnFrame     = False               ; don't turn page yet
  mres@gsnDraw      = False               ; don't draw yet
  mres@tiMainString = "Transect Location" ; title
  map = gsn_csm_map(wks,mres)         ; create map

; add polyline to map
  pres                  = True           ; polyline mods desired
  pres@gsLineColor      = "red"          ; color of lines
  pres@gsLineThicknessF = 2.0            ; line thickness   
  id = gsn_add_polyline(wks,map,(/leftlon,rightlon/),(/leftlat,rightlat/),pres)
  draw(map)    ; draws map and polyline
  frame(wks)
end