# Nearest-neighbor interpolation would do this:  data = model( 1.5 + (t-t0)/dt)
#                             This is likewise but with _linear_ interpolation.
module lint {
use tridiag
integer :: n1
logical :: inv
real    :: o1,d1
real, parameter, private :: eps = 0.01              # regularization parameter
real, dimension (:), pointer :: coordinate, weight
real, dimension (n1), allocatable :: diag           # diagonal 
real, dimension (n1), allocatable :: offd           # off-diagonal
#%  _init (inv, o1,d1,n1, coordinate, weight)
#%  _lop  ( mm, dd)
integer i, im,  id
real    f2, fx,gx, wmax
wmax = eps*maxval(weight*weight)
diag = 2.*wmax; offd = -wmax   # regularization
do id= 1, size(dd) {
        f2 = (coordinate(id)-o1)/d1;     i=f2  ;   im= 1+i
        if( 1<=im .and. im< size(mm)) { fx=(f2-i)*weight(id); gx= (1.-fx)*weight(id)        
                if( adj) {
                        mm(im  ) +=  gx * dd(id)
                        mm(im+1) +=  fx * dd(id)
			if (inv) {
				diag (im  ) += gx * gx
				diag (im+1) += fx * fx
				offd (im  ) += gx * fx
			}
		} else {
                        dd(id)   +=  gx * mm(im)  +  fx * mm(im+1)
		}
	}
}
if (adj .and. inv) call tris (diag, offd, mm) # tridiagonal solver
}