interp_2d Subroutine

    pure subroutine interp_2d(me,x,y,f,istat)

    implicit none

    class(linear_interp_2d),intent(inout) :: me
    real(wp),intent(in)                   :: x
    real(wp),intent(in)                   :: y
    real(wp),intent(out)                  :: f     !! Interpolated \( f(x,y) \)
    integer,intent(out),optional          :: istat !! `0`  : no problems,
                                                   !! `-1` : class has not been initialized

    integer,dimension(2) :: ix, iy
    real(wp) :: p1, p2
    real(wp) :: q1, q2
    integer :: mflag
    real(wp) :: fx1, fx2

    if (me%initialized) then

        call dintrv(me%x,x,me%ilox,ix(1),ix(2),mflag)
        call dintrv(me%y,y,me%iloy,iy(1),iy(2),mflag)

        q1 = (x-me%x(ix(1)))/(me%x(ix(2))-me%x(ix(1)))
        q2 = (y-me%y(iy(1)))/(me%y(iy(2))-me%y(iy(1)))
        p1 = one-q1
        p2 = one-q2

        fx1 = p1*me%f(ix(1),iy(1)) + q1*me%f(ix(2),iy(1))
        fx2 = p1*me%f(ix(1),iy(2)) + q1*me%f(ix(2),iy(2))

        f = p2*fx1 + q2*fx2
        if (present(istat)) istat = 0

    else

        if (present(istat)) istat = -1
        f = zero

    end if

    end subroutine interp_2d