function interpolate(locale, part1, llat, llon)


* learn current dimension setting

'q dims'
line2=sublin(result,2)
line3=sublin(result,3)
*say line2
origlon1=subwrd(line2,6)
origlon2=subwrd(line2,8)

*say line3
origlat1=subwrd(line3,6)
origlat2=subwrd(line3,8)

'q file'
line5=sublin(result,5)
xmax=subwrd(line5,3)
xmax=0.+xmax

* convert lat  lon into x y space

*say 'q w2gr 'llon' 'llat
    'q w2gr 'llon' 'llat
*say result
xdim=subwrd(result,3)
xdim=0.+xdim
if ( xdim <  1.) ; xdim=xmax+xdim ; endif
ydim=subwrd(result,6)

* nearest x dimensions :

x1= math_int(xdim) 
x2 = x1 + 1
if ( x2 > xmax ) ; x2  = x2 - xmax ; endif
*say "xdim = " xdim' 'x1' 'x2

* the weights are

xw1= xdim-x1 
xw2= 1.-xw1 

* nearest y dimensions :

y1= math_int(ydim) 
y2 = y1 + 1 

* the weights are

yw1= ydim-y1 
yw2= 1.-yw1 
   
_part2= '(x='x1',y='y1')*'xw2'*'yw2'+'
_part3= '(x='x2',y='y1')*'xw1'*'yw2'+'
_part5= '(x='x1',y='y2')*'xw2'*'yw1'+'
_part6= '(x='x2',y='y2')*'xw1'*'yw1

      part4=       part1%_part2
      part4= part4%part1%_part3
      part4= part4%part1%_part5
      part4= part4%part1%_part6

'set x 1'
'set y 1'
*say 'define 'locale' = 'part4
'define 'locale' = 'part4
*'d 'locale
*say subwrd(result,4)

* return to dimension setting as  before interpolation

'set lat 'origlat1' 'origlat2
'set lon 'origlon1' 'origlon2

return