pro quick_interp_tdm2,year1,year2,out_prefix,dist,$
 log=log,gs=gs,area=area,check=check,angular=angular,nostn=nostn,$
 synth_prefix=synth_prefix,limit=limit,hires=hires,pts_prefix=pts_prefix,$
 anomfac=anomfac,synthfac=synthfac,binfac=binfac,actfac=actfac,$
 dumpbin=dumpbin,dumpmon=dumpmon,dumpglo=dumpglo,info=info,ratio=ratio

; 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 ; multi factor for output bin files to write them
					   ;  as binary integer files (e.g. T would have an value
					   ;  of 10 for binfac)
if n_elements(actfac) eq 0 then actfac=1.0 ; multi factor to obtain actual (i.e. =*1) values,
					   ; required for dumping to .mon and .glo
if n_elements(synthfac) eq 0 then synthfac=1.0 ; multi factor to obtain synth file actual values
if n_elements(anomfac) eq 0 then anomfac=1.0 ; multi factor to obtain anom file actual values
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 output 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']

get_lun,LogLun
openw,LogLun,'/cru/tyn1/f014/scratch/log-qit2.dat'

if keyword_set(info) then printf,-1,"      MEAN    AV MAG   STD DEV",format="(18x,a30)"

relaxval=0.0 & if keyword_set(ratio) then relaxval=1.0

;START OF ANNUAL LOOP
;
for iy=year1,year2 do begin 
 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)
   rdbin,dummygrid,synthfile,gridsize=2.5,/quiet
   synthmiss=where(dummygrid eq -9999,nsynthmiss)
   if (nsynthmiss gt 0) then dummygrid(synthmiss)=relaxval
   dummygrid=dummygrid/synthfac
 endif else begin
   dummygrid=fltarr(144,72,12) ; 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
   dummygrid(*,*,*)=relaxval
 endelse
 
 for im=im1,im2  do begin
  jm=im

  if n_elements(nostn) eq 0 then begin
    ptsfile=pts_prefix+string(iy,".",im,".txt",form='(i4.4,a1,i2.2,a4)')
  
    UnixCommand = "cat " + ptsfile + " | wc -l"
    spawn, UnixCommand, NCount  
    NLine=fix(strtrim(NCount))
    NStn=NLine(0)
    pts1=fltarr(NStn,5)
  
    openr, lunMon, ptsfile, /get_lun
    for XLine = 0, NStn-1 do begin
      readf, lunMon, lat,lon,elv,dat,wmo, format="(2f8.2,f8.1,f13.5,i7)"
      pts1(XLine,0)=lat & pts1(XLine,1)=lon & pts1(XLine,2)=(dat*anomfac) & pts1(XLine,3)=elv
    endfor
    free_lun, lunMon
    if keyword_set(info) then begin
      mean=mean(pts1(*,2)) & avmag=mean(abs(pts1(*,2))) & stddev=stddev(pts1(*,2))
      printf,-1,"data ",iy," month ", im, mean, avmag, stddev,format="(a5,i4,a7,i2,3f10.4)"
    endif
  endif else begin
    pts1=fltarr(2,5)
    pts1(*,0)=90.0 & pts1(*,1)=0.0 & pts1(*,2)=relaxval & pts1(*,3)=0.0 
  endelse
;---------------------------------------------------
  ; 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
  
  dummymax=max(dummygrid(*,*,(im-1)))
  
  while dummymax gt -9999 do begin
    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
        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
  
    for lon=0,nlon-1 do begin		; add a selection of synth/missing points to pts2 
      for lat=0,nlat-1 do begin		;   (cannot have regular grid as it makes tri unstable)
        if (image(lon,lat) ne 160 $	; if new synth, use genuines only, else use all
        	and (dummymax eq 0.0 or dummygrid(lon,lat,(im-1)) ne 0.0)) then begin
          nel=nel+1
          pts2(nel,0:2)=[((lat*dlat)+(dlat/2.0)- 91+(2*randomu(seed))), $
      			 ((lon*dlon)+(dlon/2.0)-181+(2*randomu(seed))),dummygrid(lon,lat,(im-1))]
        endif 
      endfor
    endfor
    
    pts2=pts2(0:nel,*)
    
    if (dummymax gt relaxval) then begin
      dummygrid(*,*,(im-1))=relaxval	; on next iteration, fill with relaxation
      pts1=pts2 & pts2=0
    endif else begin
      dummygrid(*,*,(im-1))=-9999	; no more iterations
      pts1=0
    endelse
    
    dummymax=max(dummygrid(*,*,(im-1)))
  endwhile
  
  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)
  dumpvalid=where(r NE 0,dumpvalidn)
  
  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 (gs[0] eq 5.0) then SaveGrid=23
    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!
    
    if keyword_set(info) then begin
      mean=mean(r(dumpvalid)) & avmag=mean(abs(r(dumpvalid))) & stddev=stddev(r(dumpvalid))
      printf,-1,"grid "," non-zero",mean,avmag,stddev," cells=",dumpvalidn,format="(a5,4x,a9,3f10.4,a,i9)"
    endif
  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		; @@@@@@@@@ problem with save/load mon @@@@@@@@@
  dumpvalid=where(grid NE 0,dumpvalidn)
  mean=mean(grid(dumpvalid)) & avmag=mean(abs(grid(dumpvalid))) & stddev=stddev(grid(dumpvalid))
  printf,-1,"grid ",iy," non-zero",mean,avmag,stddev," cells=",dumpvalidn,format="(a5,i4,a,3f10.4,a,i9)"

  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
  dumpvalid=where(fix(grid*binfac) NE 0,dumpvalidn)
  mean=mean(grid(dumpvalid)) & avmag=mean(abs(grid(dumpvalid))) & stddev=stddev(grid(dumpvalid))
  printf,-1,"grid ",iy," non-zero",mean,avmag,stddev," cells=",dumpvalidn,format="(a5,i4,a,3f10.4,a,i9)"

  flname=strip(out_prefix+string(iy,form='(i4)'))
  wrbin,grid*binfac,flname 	; WRITE BINARY GRID (INTEGER), no compression
 endif
 
 if keyword_set(check) then goto,l2
endfor
l2:

close,LogLun
free_lun,LogLun

end