;;; NCL script for creating graph of 2D time-varying data from RCM
;;; Plots each .nc file in datadir and datadir/clim

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"

  ;;;;;
  model    =   "WRFP"
  datadir  = "~/WRFP/ncep/table3"
  outdir   = "~/WRFP/ncep/output-3"
  ;;;;;

  datafiles = systemfunc("ls "+datadir+"/*"+model+"*.nc")
  climfiles = systemfunc("ls "+outdir+"/clim/*.nc")
  allfiles = array_append_record(datafiles, climfiles, 0)
  nf       = dimsizes(allfiles)


  
  do f = 0, nf-1
  
    filename = allfiles(f)
    basename = systemfunc("basename "+filename)
   ;; get variable name from filename.
    ugh = stringtochar(basename)
    blah = indStrSubset(basename, "_")
    varname  = chartostring(ugh(0:blah(0)-1))
 

    ;; open file
    fin = addfile(filename, "r")
    
    ;; read in coordinate variables first
    
    lat = fin->lat
    lon = fin->lon
    time = fin->time
    
    nx = dimsizes(fin->xc)
    ny = dimsizes(fin->yc)
    nt = dimsizes(time)
    last = nt-1

    
    ;; read in data.
    
    data = fin->$varname$(last,:,:)
  
    
    ;; hook up coordinate variables
    
    data@lat2d = lat
    data@lon2d = lon
  

    ;; read in map projection info

    pname = data@grid_mapping
    proj = fin->$pname$


    ;; set plotting parameters
    
    res = True
    
    
    ;; map projection stuff.
    
    if (lower_case(proj@grid_mapping_name) .eq. "lambert_conformal_conic") then
      res@mpProjection	= "LambertConformal"
      res@mpLambertMeridianF	= proj@longitude_of_central_meridian
      res@mpLambertParallel1F	= proj@standard_parallel(0)
      res@mpLambertParallel2F	= proj@standard_parallel(1)
      res@tfDoNDCOverlay = True   ;; always for Lambert conformal
    end if
    
    if (lower_case(proj@grid_mapping_name) .eq. "transverse_mercator") then
      res@mpProjection	= "Mercator"
      res@mpCenterLatF	= proj@latitude_of_projection_origin
      res@mpCenterLonF	= proj@longitude_of_central_meridian
    end if

    if (lower_case(pname) .eq. "polar_stereographic") then
      res@mpProjection= "Stereographic"
      res@mpCenterLonF= proj@straight_vertical_longitude_from_pole - 360
      res@mpCenterLatF  = proj@latitude_of_projection_origin
    end if
  
    ;; /map projection stuff
      

    ;; map boundaries
    ;; expand slighltly to be sure we plot things on the borders
    res@mpLimitMode             = "Corners"
    res@mpLeftCornerLatF        = lat(0,0) -1 ;; SW corner
    res@mpLeftCornerLonF        = lon(0,0) -1 ;; SW corner
    res@mpRightCornerLatF       = lat(yc|ny-1,xc|nx-1) + 1 ;; NE corner
    res@mpRightCornerLonF       = lon(yc|ny-1,xc|nx-1) + 1 ;; NE corner

    ;; general plotting stuff
    
    res@mpFillOn                = False
    
    res@gsnMaximize             = True
    
    res@cnFillOn                = True
    res@cnLinesOn               = False
    res@cnFillMode              = "RasterFill"
    
    res@lbLabelAutoStride       = True
    res@lbBoxLinesOn            = False
    
    res@gsnSpreadColors         = True
      
    date = ut_calendar(time(last),-2)
    res@tiMainString            = basename + " - "+ date
    
      
    ;; create the plot
    
    wks = gsn_open_wks("ps", outdir+"/"+basename)
    gsn_define_colormap(wks,"nrl_sirkes") ;; 21 colors + fg/bg
      
    plot = gsn_csm_contour_map(wks, data, res)
  
    print(systemfunc("date")+" - " +basename+" done")
    
      
    ;; cleanup 
  
    delete(data)
    delete(fin)
    delete(lat)
    delete(lon )
    delete(plot)
    delete(pname)
    delete(proj)
    delete(res)
    delete(time)
    delete(wks)
    
    delete(filename)
    delete(basename)
    delete(ugh)
    delete(blah)
    delete(varname)
   
end do

exit