load "/usr/local/ncl/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/contributed.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin

  diri_500 = "/thumper/shomtm62/daily_netCDF_zg500/Composite/"   ; Obtain 500hpa input
  diri_u = "/thumper/shomtm62/daily_netCDF_Uwind/"   ; Obtain wind  input
  diri_v = "/thumper/shomtm62/daily_netCDF_Vwind/"
  prfx       = "ECPC"
  prfx1      = "DJF"

  u_wnd = (diri_u + prfx + "_" + prfx1 + ".nc")
  v_wnd = (diri_v + prfx + "_" + prfx1 + ".nc")
  hpa = (diri_500 + prfx + "_" + prfx1 + ".nc")
  
  fin_500 = addfile (hpa, "r")
  fin_u = addfile(u_wnd, "r")
  fin_v = addfile(v_wnd, "r")

  ut = fin_u->uas                ; Define input file
  u = short2flt(ut)            ; Convert from "short" to "float" if needed
  vt = fin_v->vas
  v = short2flt(vt)
  hgt0 = fin_500->zg500
  hgt1 = short2flt(hgt0) 

  lat = fin_u->lat              ; and predefine dimensions
  lon = fin_u->lon

  u_Comp = dim_avg_n_Wrap(u,0)             ; Create composite
  v_Comp = dim_avg_n_Wrap(v,0)
  hgt_Comp = dim_avg_n_Wrap(hgt1,0)

;************************************************
; create plot
;************************************************
  wks_type = "ps"
  wks = gsn_open_wks(wks_type, prfx + "_" + prfx1 + "_Overlay_Composite")      ; open a ps file
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map

  winds = True             ; plot mods desired 
  winds@gsnDraw = False
  winds@gsnFrame = False
  hgt = winds
  res = winds

  lat2d = lat
  lon2d = lon
  u_Comp@lat2d = lat
  u_Comp@lon2d = lon

;  winds@vcRefAnnoOn             = False           ; Turn off legend
  winds@vcRefMagnitudeF         = 5.0             ; define vector ref mag
  winds@vcRefLengthF            = 0.035            ; define length of vec ref
  winds@vcGlyphStyle            = "FillArrow"    ; turn on filled vectors
  winds@vcMinDistanceF          = 0.017
  winds@gsnAddCyclic            = False         ; Regional data

  winds@vcVectorDrawOrder 	= "Predraw"
  winds@tiMainString         = ""
  winds@gsnLeftString       = ""
  winds@gsnRightString      = ""

;***********************************************
; Manual Levels for 500hpa Geopotential Heights
;***********************************************
  hgt@cnLevelSelectionMode = "ManualLevels"   ; set manual contour levels
  hgt@cnMinLevelValF              = 5000.
  hgt@cnMaxLevelValF              = 6000.
  hgt@cnLevelSpacingF             = 50.       ; Contour interval.
  hgt@lbLabelStride              = 2         ; Coutour Stride.

  hgt@gsnAddCyclic          = False
  hgt@gsnSpreadColors       = True             ; use full range of colormap
  hgt@cnFillOn              = True             ; color plot desired
  hgt@cnLinesOn             = False            ; turn off contour lines
  hgt@cnLineLabelsOn        = False            ; turn off contour labels

  hgt@cnFillDrawOrder 	    = "Predraw"
  hgt@tiMainString         = ""
  hgt@gsnLeftString       = ""
  hgt@gsnRightString      = ""
  
  if (prfx.eq."MM5I").OR.(prfx.eq."UW").OR.(prfx.eq."NARR2").OR.(prfx.eq."WRFP").OR.(prfx.eq."WRFG").OR.(prfx.eq."NARR").OR.(prfx.eq."CCSM").OR.(prfx.eq."GFDL") then

   print("Projection used == Lambert Conformal")

   res@mpLimitMode           = "Corners"      ; Corners method of zooming
   res@mpLeftCornerLatF      = 22.5             ; Need to be this method for
   res@mpLeftCornerLonF      = -120.5           ; ALL map projections
   res@mpRightCornerLatF     = 48.5
   res@mpRightCornerLonF     = -61
   res@mpProjection        = "LambertConformal"
   res@mpLambertParallel1F = 30.
   res@mpLambertParallel2F = 60.
   res@mpLambertMeridianF  = -97.0

  end if

  if (prfx.eq."CRCM").OR.(prfx.eq."ECPC").OR.(prfx.eq."ECP2") then

   print("Projection used == Stereographic")

   res@mpLimitMode           = "Corners"
   res@mpLeftCornerLatF      = 22      
   res@mpLeftCornerLonF      = -118   
   res@mpRightCornerLatF     = 48
   res@mpRightCornerLonF     = -60
   res@mpProjection         = "Stereographic"
   res@mpRelativeCenterLon = True
   res@mpCenterLonF = -97
   res@mpRelativeCenterLat = True
   res@mpCenterLatF = 90.

  end if

  if (prfx.eq."HRM3").OR.(prfx.eq."RCM3") then

   print("Projection used == Mercator")

   res@mpLimitMode           = "Corners" 
   res@mpLeftCornerLatF      = 22       
   res@mpLeftCornerLonF      = -120.5  
   res@mpRightCornerLatF     = 48.5
   res@mpRightCornerLonF     = -61
   res@mpProjection       = "Mercator"
   res@mpCenterLonF        =  -97
   res@mpCenterLatF       =  47.5

  end if

  res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
  res@mpOutlineDrawOrder    = "PostDraw"      ; draw continental outline last
  res@tiMainString 	    = "500hpa heights and Winds"
  res@mpFillOn 		    = False

  wind_plot = gsn_csm_vector(wks,u_Comp(:,:),v_Comp(:,:),winds) ; plot
  hgt_plot = gsn_csm_contour(wks,hgt_Comp(:,:),hgt)
  map = gsn_csm_map(wks,res)

  overlay(map,hgt_plot)
  overlay(map,wind_plot)
  maximize_output(wks,True)

  print("")
  print("Input data = " + prfx + " and " + prfx1)

 end