pro quick_interp,year1,year2,out_prefix,dist,log=log,$
 outfac=outfac,gs=gs,area=area,check=check,angular=angular,$
 limit=limit,update=update
; runs Idl trigrid interpolation of anomaly files plus synthetic anomaly
; grids 
; first reads in anomaly file data
; the adds dummy zero value gridpoints that are 
;  further than distance (dist)
;  from any of the observed data
if n_params() lt 1 then begin
 print,'quick_interp,year1,year2,out_prefix,dist,'
 print,'  outfac=outfac,gs=gs,/area,/log,/check'
 print,'  /update,limit=limit'
 return
endif
;
close,/all
if n_elements(dist) eq 0 then xkm=500 else xkm=dist
if n_elements(constrain) eq 0 then constrain=0 else constrain=1
if n_elements(log) eq 0 then log=0 else log=1
if n_elements(outfac) eq 0 then outfac=1.0
if n_elements(pts_prefix) eq 0 then pts_prefix='../fort.'
if keyword_set(angular) eq 0 then angular=0 else angular=1
if angular eq 1 then print,'Angular distance weighting set'
if n_elements(limit) lt 3 then limit=glimit(/all)
if n_elements(gs) eq 1 then gs=[gs,gs]
if n_elements(gs) eq 0 then gs=[0.5,0.5]
if n_elements(out_prefix) eq 0 then begin
 out_prefix=' '
 read,out_prefix,prompt='Enter prefix for ouput grid file: '
 out_prefix=strip(out_prefix)
endif
im1=1 & im2=12
for iy=year1,year2 do begin
 print,iy
 grid=fltarr(360/gs(0),180/gs(1),12)-9999
 zerogrid=fltarr(144,72)
 for im=im1,im2  do begin
  jm=im
  if keyword_set(hutch) then spawn,'rm  spl'+$
   string(iy,im,form='(i4.4,i2.2)')+'.res '
  plotoff
  print,'YEAR: ',iy,'     MONTH: ',im
;read in observed anomalies
  rptsf,pts1,5,file=pts_prefix+string(iy,im,form='(i4.4,i2.2)'),format='(3f8.2,f8.4,i8)'   
  print,rnge(pts1(*,0)),rnge(pts1(*,1))
;---------------------------------------------------
  ; map all points with influence radii
  if keyword_set(test) eq 0 then begin
   set_plot,'Z'
   erase,255
   device,set_res=[144,72]
  endif else begin
   initx
   set_plot,'x'
   window,0,xsize=144,ysize=72
  endelse
  lim=glimit(/all)
  nel=n_elements(pts1(*,0))-1
  map_set,limit=lim,/noborder,/isotro,/advance,xmargin=[0,0],ymargin=[0,0],pos=[0,0,1,1]
  a=findgen(33)*!pi*2/32
  for i=0.0,nel do begin
   x=cos(a)*(xkm/110.0)*(1.0/cos(!pi*pts1(i,0)/180.0))+pts1(i,1)
   x=x<179.9 & x=x>(-179.9)
   y=sin(a)*(xkm/110.0)+pts1(i,0)
   y=y>(-89.9) & y=y<89.9
   catch,error_value
   if error_value ne 0 then begin
    print,!err_string
    print,'polyfill error',pts1(i,1),rnge(x),pts1(i,0),rnge(y),form='(a15,6f8.2)'
    error_value=0
    i=i+1
   goto,skip_poly
   endif
   polyfill,x,y,color=160
   skip_poly:
  endfor
  catch,/cancel
  image=tvrd()
  set_plot,'X'
  pts2=fltarr(200000,5)
  pts2(0:nel,*)=pts1
  nlat=n_elements(image(0,*)) & dlat=180.0/nlat
  nlon=n_elements(image(*,0)) & dlon=360.0/nlon
  nelold=nel
  for lon=0,nlon-1 do begin
   for lat=0,nlat-1 do begin
    if image(lon,lat) ne 160 and randomu(seed) lt 0.200 then begin
      nel=nel+1
      pts2(nel,0:2)=[lat*dlat+dlat/2.0-90,lon*dlon+dlon/2.0-180,zerogrid(lon,lat)]
    endif 
   endfor
  endfor
  pts2=pts2(0:nel,*)
  help,pts2
  openw,1,pts_prefix+string(iy,im,form='(i4.4,i2.2)')+'a'
  for ii=0,nel do printf,1,pts2(ii,*),format='(3f8.2,f8.4,i8)'   
  close,1
;----------------------------
  n=where(pts2(*,2) gt -9999)
  case 1 of
   keyword_set(area): begin
    bounds=[-180+gs(0)/2.0,-90+gs(1)/2.0,180-gs(0)/2.0,90-gs(1)/2.0]
    r=area_grid(pts2(n,1),pts2(n,0),pts2(n,2),gs,bounds,angular=angular)
   end
   else: begin
    bounds=[-180+gs(0)/2.0,-90+gs(1)/2.0,180-gs(0)/2.0,90-gs(1)/2.0]
    triangulate,pts2(n,1),pts2(n,0),tr,b
    r=trigrid(pts2(n,1),pts2(n,0),pts2(n,2),tr,gs,bounds,$
     xgrid=xgrid,ygrid=ygrid)
   end
  endcase
  r=r<9999 & r=r>(-9998)
  grid(*,*,im-1)=r
  print,rnge(pts2(*,2)),rnge(r)
  if keyword_set(check) then goto,l1
 endfor
 l1:
 flname=strip(out_prefix+string(iy,form='(i4)'))
 wrbin,grid*outfac,flname,/compress
 if keyword_set(check) then goto,l2
endfor
l2:
end