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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

;======================================================================
;                         The main code
;======================================================================

begin

;---A few constants
;   dir                = "/net/home/eos/ean2/SNODAS_processed/"
    dir                = "./"
    
;---Zoom in on smaller area so we can see the grid structure.
    minlat = 45    ; min(lat)
    maxlat = 45.5  ; max(lat)
    minlon = -70   ; min(lon)
    maxlon = -69.5 ; max(lon)
    bndadd = 0.25

   ;-- load avg snow depth file 
    a                   = addfile(dir+"NE_SnowDays_Nov-Apr_SNWZGT200_mm_2004-2017.nc","r")
    avgSnowDays		= a->SnowDays({minlat:maxlat},{minlon:maxlon})

    printVarSummary(avgSnowDays)

    lat = a->lat
    lon	= a->lon

;----------------------------------------------------------------------
; Plotting options section
;----------------------------------------------------------------------

       pltType = "png"                                  ; plot type
       pltDir = "./"                                   ; plot directory
       pltName = "plot_with_latlon_dots"
       pltPath = pltDir+pltName                        ; plot path

   wks = gsn_open_wks(pltType,pltPath)              ; create workstation for ps file

       gsn_define_colormap(wks,"MPL_cool")      ; define color table 

        res = True
	
	res@gsnMaximize		 = True
        res@gsnDraw              = False
        res@gsnFrame             = False

        res@cnFillOn             = True              ; color plot desired
        res@cnLinesOn            = True             ; turn off contour lines
        res@cnLineLabelsOn       = True             ; turn off contour labels
        res@cnInfoLabelOn        = False             ; turn off info label (top labels of indvid. plots)
        res@cnFillMode           = "RasterFill"      ; turn raster on
        res@cnLevelSelectionMode = "ManualLevels"    ; Set contour levels manually
        res@cnMinLevelValF       = 0                 ; minimum contour, mm
        res@cnMaxLevelValF       = 180               ; maximum contour, mm 
	res@cnLevelSpacingF	 = 20		     ; countour interval (100 mm)
	res@lbLabelBarOn	 = True
        res@pmLabelBarOrthogonalPosF = 0.2
	res@lbOrientation	 = "vertical"
	res@mpOutlineOn		 = True
	res@mpDataBaseVersion	 = "MediumRes"
	res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
	res@mpProjection	 = "CylindricalEquidistant"
       
        res@mpLimitMode           = "LatLon"           ; required
        res@mpMinLatF             = minlat ; -bndadd
        res@mpMaxLatF             = maxlat ; +bndadd
        res@mpMinLonF             = minlon ; -bndadd
        res@mpMaxLonF             = maxlon ; +bndadd
        res@mpCenterLonF          = (minlon + maxlon)*0.5
        res@mpCenterLatF          = (minlat + maxlat)*0.5
	res@pmTickMarkDisplayMode = "Always"
	res@tmXTOn		  = False
	res@tmYLOn		  = False
	
	res@gsnLeftString         = ""                ; Turn off left subtitle
        res@gsnRightString        = ""                ; Turn off right subtitle
        res@gsnMajorLatSpacing    = 1
        res@gsnMajorLonSpacing    = 1
        res@gsnMinorLonSpacing    = 1
	
        res@gsnAddCyclic          = False             ; regional grid (changes central meridian)/xwo

        msg_col   = "red"
        valid_col = "black"

        res@tiMainString = valid_col + " = valid, " + msg_col + "= missing"
        plot = gsn_csm_contour_map(wks,avgSnowDays,res)
;
; Before we draw the plot, add filled dots at the lat/lon 
; locations of the grid, using different colors for valid 
; values and missing values.  The gsn_coordinates procedure
; will draw the plot and advance the frame.
;
        mkres                          = True
        mkres@gsMarkerIndex            = 16    ; filled dot
        mkres@gsMarkerSizeF            = 5.0
        mkres@gsnCoordsMissingColor    = msg_col
        mkres@gsnCoordsNonMissingColor = valid_col
        gsn_coordinates(wks,plot,avgSnowDays,mkres)

end