; Vector with black and white contour plot for RegCM3 output
;
; MA Snyder and NS Diffenbaugh - May 15,2004

; pathnames must be changed based on local installation of NCL
load "/usr/local/ncarg/scripts/gsn_code.ncl"
load "/usr/local/ncarg/scripts/gsn_csm.ncl"
load "/usr/local/ncarg/scripts/contributed.ncl"
load "/usr/local/ncarg/scripts/shea_util.ncl"

begin

; Set path to input file
diri = "/vol1/msnyder/RegCM3/oct6/RegCM/Cal_Ncep_opt/"
fili = "ATM1960AVG.nc"

; Set path to HEAD_OUT.nc
diri1 = "/vol1/msnyder/RegCM3/oct6/RegCM/Cal_Ncep_opt/"
fili1 = "HEAD_OUT.nc"

; Pointers to input files
f  = addfile(diri+fili,"r")
f1 = addfile(diri1+fili1,"r")


; Read in variables
var    = f->U(0,0,:,:)      ; (time,lev,lat,lon)
var2   = f->V(0,0,:,:)      ; (time,lev,lat,lon)
var3   = f->TK(0,0,:,:)     ; (time,lev,lat,lon)
lat2d  = f1->XLAT(0,:,:)    ; (time,lat,lon)
lon2d  = f1->XLON(0,:,:)    ; (time,lat,lon)

; Get dimensions from variable
dimvar = dimsizes(var)
jlat   = dimvar(0)
ilon   = dimvar(1)

 wks  = gsn_open_wks ("ps", "vec_cn_bw")            ; open workstation
 gsn_define_colormap(wks,"gui_default")
 
  res                        = True               ; plot mods desired
 
  res@cnFillOn               = False
  res@cnLinesOn              = True
  res@cnLineLabelsOn         = True

  res@gsnScalarContour       = True

  res@mpGridAndLimbOn        = True               ; lat/lon lines on
  res@mpGridSpacingF         = 10.0               ; lat/lon lines every 10 deg
  res@mpGridLineDashPattern  = 2                  ; lat/lon lines as dashed
  res@pmTickMarkDisplayMode  = "Always"           ; turn on tickmarks
  
  res@tiMainString           = "RegCM3: vec_cn_bw"
  res@tiMainFontHeightF      = 0.020              ; smaller title
  
  res@gsnAddCyclic           = False              ; regional data

; !!!!! any plot of data that is on a native grid, must use the "corners"
; method of zooming in on map.
  
  res@mpLimitMode       = "Corners"            ; choose range of map
  res@mpLeftCornerLatF  = lat2d(0,0)
  res@mpLeftCornerLonF  = lon2d(0,0)
  res@mpRightCornerLatF = lat2d(jlat-1,ilon-1)
  res@mpRightCornerLonF = lon2d(jlat-1,ilon-1)
  
; The following 4 pieces of information are REQUIRED to properly display
; data on a native lambert conformal grid.  This data should be specified
; somewhere in the model itself.
  
  res@mpProjection        = "LambertConformal"
  res@mpLambertParallel1F = 30.
  res@mpLambertParallel2F = 60.
  res@mpLambertMeridianF  = -121.5
  
; usually, when data is placed onto a map, it is TRANSFORMED to the specified
; projection. Since this model is already on a native lambert conformal grid,
; we want to turn OFF the tranformation.
  
  res@tfDoNDCOverlay = True

; Turn on State borders

  res@mpPerimOn             = True
  res@mpFillOn              = False              ; turn off map fill
  res@mpOutlineBoundarySets = "GeophysicalAndUSStates"

; Vector plot specific resources

  res@vcRefMagnitudeF         = 10.0             ; define vector ref mag
  res@vcRefLengthF            = 0.045            ; define length of vec ref
  res@vcGlyphStyle            = "CurlyVector"    ; turn on curly vectors
  res@vcMinDistanceF          = 0.017            
  res@vcVectorDrawOrder       = "PostDraw"       ; draw vectors last
  
  plot = gsn_csm_vector_scalar_map(wks,var,var2,var3,res)     ; Draw contours over a map.

end