;************************************
; CSM_Graphics: transect_1.ncl
;************************************
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("h_avg_Y0191_D000.00.nc","r")
  t  = in->T(0,:,:,:)
;************************************
; calculate great circle along transect
;************************************
  leftlat  =  -60.0
  rightlat =  -30.0

  leftlon  =  -60.0
  rightlon =   20.0

  npts     =   100                    ; 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_t,t&lat_t,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&z_t           ; coordinate variable for 0th dimension only
;********************************
; create plot
;********************************
  wks = gsn_open_wks("ps","trans")       ; open a ps file
  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          = (/leftlat +", "+leftlon,rightlat+", "+rightlon/)

  res@cnFillOn            = True         ; turn on color
  res@lbLabelAutoStride   = True         ; nice label bar label stride
  res@gsnSpreadColors     = True         ; use full range of colormap
  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
;********************************
  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@gsnDraw     = False               ; don't draw yet
  mres@tiMainString= "Transect Location" ; title
  map = gsn_csm_map_ce(wks,mres)         ; create map

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