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

;  a = addfile("/home/adam/Documents/wrf_270718_06z/rain2.nc","r")   ; Open a file
  a = addfile("./rain2.nc","r")   ; Open a file


res = True

res@cnFillOn                     = True
res@cnFillMode                   = "AreaFill"
res@cnFillPalette                = "prcp_1"

res@cnLevelSelectionMode         = "ExplicitLevels"
res@cnLevels                     = (/0.1,0.3,0.5,1,3,5,7,9,11,13,15,20,25,30,35/)

;res@cnMinLevelValF               = 0.1     ; This has no effect if cnLevelSelectionMode = "ExplicitLevels"
;res@cnRasterSmoothingOn          = True      ; this has no effect unless res@cnFillMode = "RasterFill"

res@cnConstFLabelBackgroundColor = False
res@cnOutOfRangePerimOn          = False
res@cnLinesOn                    = False

res@lbOrientation                = "Vertical"    ; default is "horizontal"

;
; This appears to be a WRF output file, so use wrf_map_resources to
; set the correct map projections. This function will set the two
; lambert parallels and meridian
;
res                              = wrf_map_resources(a,res)    
;print(res)    ; Uncomment these two lines if you want to see
;exit          ; what resources wrf_map_resources is setting.
res@tfDoNDCOverlay               = True        ; Tells NCL we are using a native map projection.
res@gsnAddCyclic                 = False       ; Don't try to add a longitude cyclic point
res@mpGeophysicalLineColor       = "Black"
res@mpGeophysicalLineThicknessF  = 1.5

;---Not needed since wrf_map_resources was used.
;res@mpProjection                 = "LambertConformal"
;res@mpLimitMode                  = "Corners"
;res@mpLambertParallel1F          = a@TRUELAT1
;res@mpLambertParallel2F          = a@TRUELAT2
;res@mpLambertParallel1F          = a@STAND_LON

res@mpFillOn                     = False
res@mpDataBaseVersion            = "HighRes"
res@mpOutlineDrawOrder           = "PostDraw"

res@tiMainString  = "Rain hourly"

type = "png"
wks = gsn_open_wks(type,"raintest")       ; Create a plot workstation

ntime = dimsizes(a->TOTRAIN(:,0,0))
total_rain = a->TOTRAIN
do i = 1, ntime-1
  ter  = total_rain(i,:,:) - total_rain(i-1,:,:)
;  res@gsnLeftString = sprinti("index %2d",i)

;
; The four corner settings are set by wrf_map_resources, so I commented
; them out here. You would need to set these inside the do loop only
; if the lat/lon domain changes for each iteration across time.
; Also, since we are plotting in the native map projection, the lat2d/lon2d
; special attributes are not needed.
;
;  lat = a->XLAT(i,:,:)
;  lon = a->XLONG(i,:,:)
;  dims = dimsizes(ter)
;  r = dimsizes(dims)
;  r1 = dims(r-1)
;  r2 = dims(r-2)
;
;  res@mpLeftCornerLatF  = lat(0,r1-1)
;  res@mpLeftCornerLonF  = lon(0,r1-1)
;  res@mpRightCornerLatF = lat(r2-1,0)
;  res@mpRightCornerLonF = lon(r2-1,0)
;  ter@lat2d = lat
;  ter@lon2d = lon

  plot = gsn_csm_contour_map(wks,ter,res)           ; Create plot
end do

end