; filepath = "/home/debian/WRF/CASES/kinugawa/output/ \
;              3.33km_1day_FNL_mp10_origin_SST_10output/\
;              wrfout_d01_2015-09-09_00:00:00"
              
begin

;---Set bounding box
; bounding box for cover kinu watershed
  lats		= (/ 30.00, 40.00 /)
  lons		= (/ 132.00, 142.00 /)

  
;---We generate plots, but what kind do we prefer?
  type				= "x11"
  output_name 			= "./" 
  wks 				= gsn_open_wks(type,output_name)
  ;gsn_define_colormap(wks,"rainbow")  ; overwrite the .hluresfile color map


;---Set some basic resources
  res 				= True
  res@NoHeaderFooter		= True            ; Switch headers and footers off

  pltres 			= True
  pltres@PanelPlot 		= True      ; Indicate these plots are to be paneled.

  mpres 			= True
  mpres@gsnMaximize          = True
  mpres@mpMinLatF            = 20.0
  mpres@mpMaxLatF            = 50.0
  mpres@mpMinLonF            = 120
  mpres@mpMaxLonF            = 150
  mpres@mpOutlineOn          = True
  mpres@mpFillOn       		= False
  mpres@mpDataBaseVersion    	= "MediumRes"


;---Get xlat,xlon from model for interpolation
  filepath	= "./wrfout_d01_2015-09-08_12:00:00"
  a		= addfile(filepath,"r")
  time		= 0
  xlat2d_m	= wrf_user_getvar(a,"XLAT",time)
  xlon2d_m	= wrf_user_getvar(a,"XLONG",time)
  xland		= wrf_user_getvar(a,"XLAND",time)
  dims2d	= dimsizes(xlat2d_m)
  
  ;printVarSummary(xlat2d_m)
  ;printVarSummary(xlon2d_m)
  ;printVarSummary(xland)
  ;print(dims2d)
  ;print (xlat2d_m)

;---Set zoomin indices
; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y)
; subtract one since we want to use it as an index in NCL
  loc		= wrf_user_ll_to_ij(a, lons, lats, True)
  loc		= loc - 1
  x_start	= loc(0,0)
  x_end		= loc(0,1)
  y_start	= loc(1,0)
  y_end		= loc(1,1)
  mpres@ZoomIn	= True
  mpres@Xstart	= x_start
  mpres@Ystart	= y_start
  mpres@Xend	= x_end
  mpres@Yend	= y_end
  
  print ("---Bounding Box Information---")
  print ("x_start = " + x_start + "	x_end = " + x_end)
  print ("y_start = " + y_start + "	y_end = " + y_end)
  print ("------------------------------")
  
  
;---Loop of time   
  times		= wrf_user_getvar(a,"times",-1)  ; get all times in the file
  plots		= new ( 4, graphic )
  
  do itime = 0,0
    if (itime.eq.0)then
      asc_file = "./crain_composite_2015.09.09.00.00_UTC"
    end if

    stime = get_cpu_time()
    data	= asciiread(asc_file, -1, "float")
    header	= readAsciiHead(asc_file,2)
    date	= str_get_field(header(1),2," ")
      
; read data
    radar	= data(7::3)
    lat_obs	= data(8::3)
    lon_obs	= data(9::3)
      
; These long_name and units attributes are not necessary. I
; did this so I got better output from print statements below.

    radar@long_name	= "cband radar observations"
    radar@units	= "mm/h"
    lat_obs@long_name	= "observed latitudes"
    lat_obs@units	= "degrees_north"
    lon_obs@long_name	= "observed longitudes"
    lon_obs@units	= "degrees_east"

    radar@_FillValue	= -9999
    radar@_FillValue	= -999

    etime = get_cpu_time()
    print("Reading the data took " + (etime-stime) + " seconds.")

    printVarSummary(radar)
    printVarSummary(lat_obs)
    printVarSummary(lon_obs)
    printMinMax(radar,0)
    printMinMax(lat_obs,0)
    printMinMax(lon_obs,0)
      
;---Get slice of lat/lon data at halfway points.
    nlat   = dims2d(0)
    nlon   = dims2d(1)
    xlat_m = xlat2d_m(:,nlon/2)
    xlon_m = xlon2d_m(nlat/2,:)

    ;print (nlat+ " "+ nlon)
    printVarSummary(xlat_m)
    printVarSummary(xlon_m)
    printMinMax(xlat_m,0)
    printMinMax(xlon_m,0)
    
    rscan = (/1,0.5,0.25/)
;    rscan = (/0.8,0.4,0.20/)
;    rscan = (/0.25,0.125,0.05/)
;    rscan   = (/1.0/)
    grid    = obj_anal_ic(lon_obs,lat_obs,radar,xlon_m,xlat_m,rscan,False)

    printVarSummary(grid)
    printMinMax(grid,0)
    ;print (grid)

    
;---Generate contours.
    res1			= res
    res1@cnFillOn		= True          
    res1@cnLevelSelectionMode	= "ExplicitLevels"
    res1@cnFillOn		= True
    res1@cnFillMode		= "RasterFill"     ; handles memory better, but can produce blocky contours

    res1@cnLinesOn		= False
    res@cnLevels		= (/0.1, 1, 5, 10, 20, 50, 100, 210/)
    res@cnFillColors		= (/"white","cyan", "blue","green","yellow","darkorange", \
				    "coral","red","blueviolet"/)

;    var_zoom	= grid(15:40,15:40)
    var_zoom	= grid(y_start:y_end,x_start:x_end)
    printVarSummary(var_zoom)
    printMinMax(var_zoom,0)
    contour	= wrf_contour(a,wks,var_zoom,res1)
  
;---Overlay contours on a map
    ;pltres@CommonTitle	= True
    ;pltres@PlotTitle	= date
    ;plot	= wrf_map_overlays(a,wks,(/contour/),pltres,mpres)

      
;---End loop of time
  end do
  
  
end