Linear interpolation

A more accurate interpolation method is linear interpolation. To interpolate linearly, we locate each data coordinate between two model points on the grid. The data ordinate is defined by drawing a straight line through the model values or, in other words, by weighting them proportionally to the distance from the coordinate. The adjoint uses the same weights to distribute the data ordinate between the two grid points.

The interfaces in module lint1 are exactly the same as those in bin1, but the actual implementation details are slightly different. Using the similar interfaces makes these modules replaceable; in any application, one of them can be replaced by the other with a simple change of the use statement. The nearest neighbor interpolation is faster, while the linear interpolation provides a more accurate result.

module lint1	
  use adj_mod

integer,                             private :: nd
  integer, dimension (:), allocatable, private :: x
  real,    dimension (:), allocatable, private :: w1, w2
  logical, dimension (:), allocatable, private :: m

contains

subroutine lint1_init (coord, o1, d1, n1, weight)
    real, dimension (:), intent (in) :: coord, weight
    real,                intent (in) :: o1, d1
    integer,             intent (in) :: n1
    optional                         :: weight

nd = size (coord) 
    if (.not. allocated (x)) allocate (x (nd), m (nd), w1 (nd), w2 (nd))

w1 = (coord - o1) / d1
    x = w1     ; w2  = w1 - x
    x = x  + 1 ; w1 = 1. -  w2

m = (x < 1 .or. x >= n1)
    if (present (weight)) then
        m = m .or. (weight == 0.)
       where (.not. m)
          w1 = w1 * weight ; w2 = w2 * weight
       end where
    end if
  end subroutine lint1_init

function lint1_op (adj, add, mod, ord) result (stat)
    integer              :: stat
    logical, intent (in) :: adj, add
    real, dimension (:)  :: mod, ord

integer              :: i, i1, i2

call adjnull (adj, add, mod, ord)

do i = 1, nd ; if (m (i)) cycle
       i1 = x (i) ; i2 = i1 + 1
       if (adj) then 
          mod (i1) = mod (i1) + w1 (i) * ord (i )
          mod (i2) = mod (i2) + w2 (i) * ord (i )
       else
          ord (i ) = ord (i ) + w1 (i) * mod (i1) + &
                                w2 (i) * mod (i2)
       end if
    end do

stat = 0
  end function lint1_op

subroutine lint1_close ()
    deallocate (x, m, w1, w2)
  end subroutine lint1_close

end module lint1