;---------------------------------------------------------------------------
;
; 2012-dec-27	Demo program to plot original Daymet tiles,
;		  and show the surrounding no-data areas.
;		By Dave Allured, NOAA/PSD/CIRES Climate Analysis Branch.
;
; Usage:
;
; 1.  Go to a directory containing a number of Daymet tile
; directories.
;
; 2.  Make a text file "list1", containing the tile directory
; names.  E.g.:  ls -1d 1*_1980 > list1
;
; 3.  Run this script in the same directory.  The program is
; currently set to output frames in an X11 window.  Click on
; each frame to advance to the next frame.
;
;---------------------------------------------------------------------------

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"    

begin
  dirs   = asciiread ("list1", -1, "string")
  nplots = dimsizes (dirs)
  var    = "prcp"

  wks = gsn_open_wks ("x11", "daymet-orig")
  gsn_define_colormap (wks, "WhViBlGrYeOrRe")	; set color map
  
  do pi = 0, nplots-1
    infile = dirs(pi) + "/" + var + ".nc"
    f = addfile (infile, "r")

    x = f->$var$(24,:,:)			; read sample data grid
    date = "January 25"				; note the corresponding date
    
    x@lat2d = f->lat				; read 2-D coordinates
    x@lon2d = f->lon				; and attach to data array

    print ("Tile = " + dirs(pi) + sprintf (", data range = %0.9g", min(x)) \
      + sprintf (" to %0.9g", max(x)) )
    
    x = where (ismissing (x), -900, x)		; make no-data areas visible

    res				= True
    res@mpMinLatF		=   39.0	; Wyoming and vicinity
    res@mpMaxLatF		=   47.0
    res@mpMinLonF		= -113.5
    res@mpMaxLonF		= -103.5

    res@cnFillOn		= True
    res@cnLinesOn		= False
    res@gsnSpreadColors		= True
    res@gsnSpreadColorStart	= 8
  
    res@cnLevelSelectionMode	= "ManualLevels"
    res@cnLevelSpacingF		=  1		; set fixed contour range
    res@cnMinLevelValF		= -1		; for precip
    res@cnMaxLevelValF		= 15

    res@mpFillOn		= False		; turn off map fill
    res@mpOutlineOn		= True		; show state outlines
    res@mpOutlineBoundarySets	= "geophysicalandusstates"

    res@tiMainString		= "Daymet Tile " + dirs(pi) + ", " + date
    res@tiMainFontHeightF	= 0.020
  
    plot = gsn_csm_contour_map_ce (wks, x, res)
    delete (x)
  end do
end