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/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"


  vars = (/"PRES_221_SFC"/) 
 
   NARR_files = systemfunc("ls merged_AWIP32.2003072118.RS.sfc*")
   f=addfile(NARR_files(0)+".grb","r")
   vars_to_check = getfilevarnames(f)
;   print(vars_to_check)
   list_filevars(f)
;  print(lv_ISBL3)
   
do ifile = 0 ,  dimsizes(NARR_files)-1
   f=addfile(NARR_files(ifile)+".grb","r")

  res                       = True             ; plot mods desired
  res@cnFillMode       = "RasterFill" ;  = True
;  res@gsnPolar   = "NH"                          ; specify the hemisphere
;  res@mpMinLatF  = 60                           ; minimum lat to plot
;  res@mpCenterLonF = -90                          ; center longitude

  res@gsnFrame  = False    
  res@gsnMaximize = True
  res@gsnPaperOrientation   = "portrait"
  res@mpOutlineBoundarySets = "AllBoundaries"
  res@gsnSpreadColors       = True             ; use full range of colormap
  res@cnFillOn              = True             ; color plot desired
  res@cnLinesOn             = False            ; turn off contour lines
  res@cnLineLabelsOn        = False            ; turn off contour labels
;  res@cnFillDrawOrder  = "PostDraw"          ; areas before map gets set
  res@mpOutlineDrawOrder    = "PostDraw"      ; draw continental outline last
  res@mpLandFillColor = "Transparent"
  res@mpGeophysicalLineColor = "black"
;  res@mpGeophysicalLineThicknessF = 1.5          ; thickness of outlines
  res@mpOutlineOn = True 

  res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries
  res@tmXTOn = False            ; turn off top   labels
  res@tmYROn = False            ; turn off right labels
  res@lbTitleFontHeightF= .018 
  res@lbLabelFontHeightF = 0.018
  res@gsnRightStringFontHeightF = 0.018
  res@gsnLeftStringFontHeightF = 0.018
  res@cnInfoLabelOn       = False 
  res@cnInfoLabelOrthogonalPosF = 0.0
;  res@cnInfoLabelString   = "Min= $ZMN$ Max= $ZMX$"
  res@gsnAddCyclic       = True 
  
  res@mpLimitMode           = "LatLon"  
   res@vcRefMagnitudeF          = 8.             ; define vector ref mag
   res@vcRefLengthF             = 0.025          ; define length of vec ref
   res@vcRefAnnoOrthogonalPosF  = -1.            ; move ref vector
   res@vcRefAnnoParallelPosF  = 0.95            ; move ref vector
   res@vcMinDistanceF           = 0.025            ; larger means sparser
   res@vcLineArrowHeadMaxSizeF  = 0.0075          ; default: 0.05 (LineArrow), 0.012 (CurlyVector)
   res@vcGlyphStyle             = "CurlyVector"   ; default: "LineArrow"
   res@gsnScalarContour      = True               ; contours desired
      gsres               = True
      gsres@gsMarkerIndex = 16          ; circle at first
      gsres@gsMarkerThicknessF = 1
      gsres@gsMarkerSizeF = 0.01
      gsres@gsMarkerColor      = "black"
  tres       =  True
  tres@txFontHeightF = 0.018
;   res@mpCenterLonF  = 180.
;   res@mpMinLatF = 25
;   res@mpMaxLatF =45 
;   res@mpMinLonF = 360-130 
;   res@mpMaxLonF = 360-100 
;  res@mpLimitMode           = "LatLon"
  lat2d   = f->gridlat_221
  lon2d   = f->gridlon_221
  dimlc   = dimsizes(lat2d)
  nlat    = dimlc(0)
  mlon    = dimlc(1)
  res@mpLimitMode            = "Corners"            ; choose range of map
  res@mpLeftCornerLatF       = lat2d(0,0)
  res@mpLeftCornerLonF       = lon2d(0,0)
  res@mpRightCornerLatF      = lat2d(nlat-1,mlon-1)
  res@mpRightCornerLonF      = lon2d(nlat-1,mlon-1)
  res@tfDoNDCOverlay         = True  
  res@mpProjection           = "LambertConformal" 
  res@mpLambertParallel1F    = lat2d@mpLambertParallel1F
  res@mpLambertParallel2F    = lat2d@mpLambertParallel2F
  res@mpLambertMeridianF     = lat2d@mpLambertMeridianF

  res@cnLevelSelectionMode  = "AutomaticLevels"
;   gridlat_221 = f->gridlat_221
;   gridlon_221 = f->gridlon_221
;   printVarSummary(gridlat_221)
      var_value=f->PRES_221_SFC(:,:)
      printVarSummary(var_value)
;      var_value@lat2d = gridlat_221
;      var_value@lon2d = gridlon_221

    figurename = "wrfout_d01_PRES_221_SFC_"+ifile
      wks = gsn_open_wks("eps" ,figurename)           ; ps,pdf,x11,ncgm,eps
      gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map
;      res@gsnRightString    = var_value@long_name ;+ " at "+lv_ISBL3(ilevel)+"mb, "+ var_value@units
       res@cnLevelSelectionMode  = "ManualLevels"
       res@lbLabelAutoStride = True
       res@cnLevelSpacingF   =   1. 
      plot      = gsn_csm_contour_map(wks,var_value(::2,::2),res)
      frame(wks)
    delete(var_value)
   system("convert -density 300 -trim -resize 1000x800 "+figurename+".eps /nsftor/xhu/public_html/HurricaneImpactonO3/observations/"+figurename+".png")
   print("finish plotting "+figurename+".eps")
   system("rm "+figurename+".eps")

end do ; do ifile