;http://www.ncl.ucar.edu/Document/Functions/Contributed/obj_anal_ic_Wrap.shtml
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

; 203005
;Assume zlon, zlat and zVal are one-dimensional (1D) arrays; let lon and lat be 1D arrays specifying grid locations. 

;rscan
;A one-dimensional array of length K specifying the successive radii of influence. The maximum size of rscan is 10. Typically, rscan contains one-to-four elements. Must be expressed in degrees of latitude and must be monotonically decreasing. eg: rscan = (/10, 5, 3/)

;option
;If option=False, this function will operate under default mode. If option=True, then this variable may have associated with it the attribute blend_wgt. This specifies how successive new estimates are 'blended'(local smoother) with previously derived estimates. These values must be 0.0 < blend_wgt < 1.0. A value of 1.0 means that the current interpolated value will be used. Otherwise:

;      new_value(n) = blend_wgt(n)*current_estimate + (1-blend_wgt(n))*value(n-1)
;The size of blend_wgt must be the same size as rscan.

;  grid = obj_anal_ic(zlon,zlat,zVal, lon,lat, rscan, False) will use default behavior and return a 2D array of size N x M M.

begin

; ascii file
a_asc = (/"jan26-28_2015.csv"/)
;---Read the data in as a giant 1D array
    values = asciiread(a_asc,-1,"float")
    nvals  = dimsizes(values)
    print("nvals = " + nvals)
;   nx=dims(1)
;   ny=dims(0)
;---Calculate ntime based on the number of values we read in.
;   ntxy  = nblk+1
; create variables to hold the ascii file variables.
;   num_obs = 203004
    num_obs = 2904
    num_rows =  3
;
      lat_o    = new((/num_obs/),float)
      lon_o    = new((/num_obs/),float)
      obs_o    = new((/num_obs/),float)
      z1 = asciiread("jan26-28_2015.txt",(/num_obs,num_rows/),"float")
      lat_o(:)=z1(:,0)
      lon_o(:)=z1(:,1)
      obs_o(:)=z1(:,2)
;     print(z1)
;     print("obs")
      print(obs_o)
      filename = "xlat_xlon_01_26_2015.nc"
      fin = addfile(filename,"r")
      xlat2d_m = fin->xlat  
      xlon2d_m = fin->xlong
      dims2d=dimsizes(xlat2d_m)
      ny = dims2d(0)
      nx = dims2d(1)
      xlat_m  = new((/ny/),float)      
      xlon_m  = new((/nx/),float)      
      xlat_m = xlat2d_m(:,nx/2)
      xlon_m = xlon2d_m(nx/2,:)
      print("xlat_m = " + xlat_m)
      print("xlon_m = " + xlon_m)
      

;  rscan = (/1, 0.5, 0.2/)
   rscan = (/2, 1.0, 0.2/)
   rscan = (/1/)
   print("lon_o = " + lon_o)
   print("lat_o = " + lat_o)
   print("obs_o = " + obs_o)
;  grid = obj_anal_ic(zlon,zlat,zVal, lon,lat, rscan, False)
   grid = obj_anal_ic(lon_o,lat_o,obs_o, xlon_m,xlat_m, rscan, False)
   dims = dimsizes(grid)
   print(dims)
   printMinMax(grid,False)
   do j=1,ny-1
   do i=1,nx-1
   if (grid(j,i).gt.0)then 
    print("grid = " + grid(j,i))
   end if
   end do
   end do
;  print("grid  = " + grid)



end