pro quick_interp_tdm,year1,year2,out_prefix,dist,log=log,$
 outfac=outfac,gs=gs,area=area,check=check,angular=angular,$
 limit=limit,hires=hires,pts_prefix=pts_prefix
; 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,'  /hires'
 return
endif
;
close,/all
; 
; set up a few defaults
if n_elements(dist) eq 0 then xkm=500 else xkm=dist ; correlation decay distance
if n_elements(constrain) eq 0 then constrain=0 else constrain=1 ; obsolete 
if n_elements(log) eq 0 then log=0 else log=1 ; obsolete
if n_elements(outfac) eq 0 then outfac=1.0 ; multiplication factor for output files so that can be written
					   ;  as binary integer files (e.g. temperature would have an value
					   ;  of 10 for outfac)
if n_elements(pts_prefix) eq 0 then pts_prefix='../fort.' ; option to have an alternative input file prefix
if keyword_set(angular) eq 0 then angular=0 else angular=1 ; use angular distance interpolation rather 
							   ; than triangulation (slower!)
if angular eq 1 then print,'Angular distance weighting set'
limit=glimit(/all) ; sets limit to global field
if n_elements(gs) eq 1 then gs=[gs,gs] ; sets grid spacing for output grid, or if gs not set
if n_elements(gs) eq 0 then gs=[0.5,0.5] ; sets gridspacing to 0.5 deg lat/lon
if n_elements(out_prefix) eq 0 then begin ; set up output file prefix
 out_prefix=' '
 read,out_prefix,prompt='Enter prefix for ouput grid file: '
 out_prefix=strip(out_prefix)
endif
im1=1 & im2=12

;START OF ANNUAL LOOP
;
for iy=year1,year2 do begin 
 print,iy
 grid=fltarr(360/gs(0),180/gs(1),12)-9999  ; set up output grid
 zerogrid=fltarr(144,72) ; set up a "zero anomaly" grid for infilling spaces with missing data
			 ; missing data defined as areas of grid further than the decay distance from any
			 ; real station point

; START OF MONTHLY LOOP
 for im=im1,im2  do begin
  jm=im
  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 - that is with decay distance [IDL variable is decay]
  ;  this is done in the virtual Z device, and essentially finds all points on a 2.5 deg grid that 
  ;  fall outside station influence
  if keyword_set(test) eq 0 then begin
   set_plot,'Z' ; set plot window to "virtual"
   erase,255 ; clear current plot buffer, set backgroudn to white
   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; avoids a bug in IDL that throws out an occasional plot error in virtual window
   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()  ; read plot of missing data as grid and then (in lines below) find those areas that are
		;  white (i.e. value 255) and add these to an array of "missing" input points (pts2)
  set_plot,'X'
  pts2=fltarr(200000,5)
  pts2(0:nel,*)=pts1  ; add real input data to pts2
  nlat=n_elements(image(0,*)) & dlat=180.0/nlat
  nlon=n_elements(image(*,0)) & dlon=360.0/nlon
  nelold=nel
  ; add a random selection of missing points to pts2 (can have regular grid as it makes
  ; trianglulation unstable
  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
  ; write new pts2 to file - NB this is a file of monthly data - we are still within the year loop
  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  ; ANGULAR DISTANCE WEIGHTING
    bounds=[-180+gs(0),-90+gs(1),180-gs(0),90-gs(1)]
    r=area_grid(pts2(n,1),pts2(n,0),pts2(n,2),gs*2.0,bounds,dist,angular=angular)
    r=rebin(r,720,360)
   end
   else: begin
    if keyword_set(hires) eq 0 then begin  ; LOW RESOLUTION TRIANGULATION
     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)
    endif else begin                       ; HIGH RESOLUTION TRINAGULATION, THEN RESAMPLE TO LOW RES
     bounds=[-180+gs(0)/10.0,-90+gs(1)/10,180-gs(0)/10,90-gs(1)/10]
     triangulate,pts2(n,1),pts2(n,0),tr,b
     r=trigrid(pts2(n,1),pts2(n,0),pts2(n,2),tr,gs/5.0,bounds,$
      xgrid=xgrid,ygrid=ygrid)
     if keyword_set(check) then help,r
     r=rebin(r,720,360)
     if keyword_set(check) then help,r
    endelse
   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:
 if keyword_set(smooth) then dist_smooth,grid,dist
 flname=strip(out_prefix+string(iy,form='(i4)'))
 wrbin,grid*outfac,flname,/compress  ;WRITE BINARY GRID (INTEGER)
 if keyword_set(check) then goto,l2
endfor
l2:
end