module int1 {                                     # generic 1-d interpolation
  integer                                 :: nd, nf
  integer, dimension (nd),    allocatable :: x
  logical, dimension (nd),    allocatable :: m
  real,    dimension (nf,nd), allocatable :: w
#% _init (coord, o1, d1, n1, interp, nd, nf, weight, ww)
    real, dimension (:),   intent (in)  :: coord, weight
    real, dimension (:,:), intent (out) :: ww
    real,                  intent (in)  :: o1, d1
    integer,               intent (in)  :: n1
    optional                            :: weight, ww
    interface {
       integer function interp (x, w) {
         real, intent (in)   :: x
         real, dimension (:) :: w
       }
    }
    integer                             :: id, ix, i1, i2, stat, j
    real                                :: rx
    if (present (ww)) ww = 0.
    do id = 1, nd {
       rx = (coord (id) - o1)/d1 ; ix = rx 
       rx -= ix ; x (id) = ix - 0.5*nf + 1.
       m (id) = .true. ; w (:, id) = 0. 
       i1 = x (id) + 1 ; i2 = x (id) + nf
       if (i1 < 1 .or. i2 > n1) cycle
       m (id) = .false. ; stat = interp (rx, w (:,id))
       if (present (weight)) w (:,id) = w (:,id) * weight (id)
       if (present (ww)) 
          do j = 1, nf
             ww (i1:i2+1-j,j) += w (1:1+nf-j,id) * w (j:nf,id)
    }
#% _lop (mod, ord)
    integer :: i, i1, i2, j
    do i = 1, nd { if (m (i)) cycle 
       i1 = x (i) + 1 
       i2 = x (i) + nf
       if (adj) 
          mod (i1:i2) +=  w (:,i)  * ord (i)
       else
          ord (i) += sum (mod (i1:i2) * w (:,i))
    }
}