pro quick_interp_tdm,year1,year2,out_prefix,dist,$
 log=log,binfac=binfac,actfac=actfac,gs=gs,area=area,check=check,angular=angular,synth_prefix=synth_prefix,$
 limit=limit,hires=hires,pts_prefix=pts_prefix,dumpbin=dumpbin,dumpmon=dumpmon,dumpglo=dumpglo
; runs Idl trigrid interpolation of anomaly files plus synthetic anomaly
; grids 
; first reads in anomaly file data
; the adds dummy gridpoints that are 
;  further than distance (dist)
;  from any of the observed data
; TDM: the dummy grid points default to zero, but if the synth_prefix files are present in call,
;  the synthetic data from these grids are read in and used instead
if n_params() lt 1 then begin
 print,'quick_interp,year1,year2,out_prefix,dist,'
 print,'  binfac=binfac,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(binfac) eq 0 then binfac=1.0 ; multiplication factor for output bin files to write them
					   ;  as binary integer files (e.g. temperature would have an value
					   ;  of 10 for binfac)
if n_elements(actfac) eq 0 then actfac=1.0 ; multiplication factor to obtain actual (i.e. =*1) values,
					   ; required for dumping to .mon and .glo
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

monthname=strarr(12)
monthname[*]=['01','02','03','04','05','06','07','08','09','10','11','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
 if keyword_set(synth_prefix) then begin
   synthfile=strtrim(synth_prefix,2)+strtrim(string(iy),2)
   print,synthfile
   rdbin,dummygrid,synthfile,gridsize=2.5
 endif else begin
   dummygrid=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
 endelse
 
; START OF MONTHLY LOOP
 for im=im1,im2  do begin
  jm=im
  plotoff
  print,'YEAR: ',iy,'     MONTH: ',im
;read in observed anomalies
  ptsfile=pts_prefix+string(iy,im,form='(i4.4,i2.2)')
  print,strtrim(ptsfile,2)
  rptsf,pts1,5,file=ptsfile,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,dummygrid(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(dumpglo) then begin
    r=r*actfac
    SaveGrid=-1
    if (gs[0] eq 0.5) then SaveGrid=12
    if (gs[0] eq 2.5) then SaveGrid=22
    if (SaveGrid eq -1) then print, "@@@@@ grid not found for save to .glo @@@@@"
    SaveFile=strip(out_prefix+monthname(im-1)+'.'+string(iy,form='(i4)')+'.glo')
    SaveTitle=strip(monthname(im-1)+string(iy,form='(i5)'))
    SaveGlo,SaveGrid,r,CallFile=Savefile,CallTitle=SaveTitle		; plot routine - change at will!
  endif
  
  if keyword_set(check) then goto,l1
 endfor
 l1:
 if keyword_set(smooth) then dist_smooth,grid,dist
 
 if keyword_set(dumpmon) then begin
  savefile=strip(out_prefix+string(iy,form='(i4)')+'.mon')
  print,savefile
  openw, lunSave, savefile, /get_lun
  
  for xmon = 0, 11 do begin
   r=grid(*,*,xmon)
   r=r*actfac
  
   for xlat = (180/gs(0))-1, 0, -1 do begin
     for xlon = (360/(gs(0)*2)), (360/gs(0))-1, 6 do begin
       xlon0=xlon+0 & xlon1=xlon+1 & xlon2=xlon+2 & xlon3=xlon+3 & xlon4=xlon+4 & xlon5=xlon+5
       printf,lunSave,r(xlon0,xlat),r(xlon1,xlat),r(xlon2,xlat), $
       		r(xlon3,xlat),r(xlon4,xlat),r(xlon5,xlat),format='(6e12.4)'
     endfor
     for xlon = 0, (360/(gs(0)*2))-1, 6 do begin
       xlon0=xlon+0 & xlon1=xlon+1 & xlon2=xlon+2 & xlon3=xlon+3 & xlon4=xlon+4 & xlon5=xlon+5
       printf,lunSave,r(xlon0,xlat),r(xlon1,xlat),r(xlon2,xlat), $
       		r(xlon3,xlat),r(xlon4,xlat),r(xlon5,xlat),format='(6e12.4)'
     endfor
   endfor
  endfor
  free_lun, lunSave
 endif
 
 if (keyword_set(dumpbin)) then begin
  flname=strip(out_prefix+string(iy,form='(i4)'))
  print,flname
  wrbin,grid*binfac,flname ;WRITE BINARY GRID (INTEGER), no compression
;  wrbin,grid*binfac,flname,/compress  ;WRITE BINARY GRID (INTEGER)
 endif
 
 if keyword_set(check) then goto,l2
endfor
l2:
end