interp_5d Subroutine

    pure subroutine interp_5d(me,x,y,z,q,r,f,istat)

    implicit none

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

    integer,dimension(2) :: ix, iy, iz, iq, ir
    real(wp) :: p1, p2, p3, p4, p5
    real(wp) :: q1, q2, q3, q4, q5
    integer :: mflag
    real(wp) :: fx1111, fx2111, fx1211, fx2211, fx1121, fx2121, fx1221, fx2221, &
                fxy111, fxy211, fxy121, fxy221, fxyz11, fxyz21, fxyzq1, fx1112, &
                fx2112, fx1212, fx2212, fx1122, fx2122, fx1222, fx2222, fxy112, &
                fxy212, fxy122, fxy222, fxyz12, fxyz22, fxyzq2

    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)
        call dintrv(me%z,z,me%iloz,iz(1),iz(2),mflag)
        call dintrv(me%q,q,me%iloq,iq(1),iq(2),mflag)
        call dintrv(me%r,r,me%ilor,ir(1),ir(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)))
        q3 = (z-me%z(iz(1)))/(me%z(iz(2))-me%z(iz(1)))
        q4 = (q-me%q(iq(1)))/(me%q(iq(2))-me%q(iq(1)))
        q5 = (r-me%r(ir(1)))/(me%r(ir(2))-me%r(ir(1)))
        p1 = one-q1
        p2 = one-q2
        p3 = one-q3
        p4 = one-q4
        p5 = one-q5

        fx1111 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(1))
        fx2111 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(1))
        fx1211 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(1))
        fx2211 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(1))
        fx1121 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(1))
        fx2121 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(1))
        fx1221 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(1))
        fx2221 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(1))
        fx1112 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(2))
        fx2112 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(2))
        fx1212 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(2))
        fx2212 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(2))
        fx1122 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(2))
        fx2122 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(2))
        fx1222 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(2))
        fx2222 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(2))

        fxy111 = p2*fx1111 + q2*fx2111
        fxy211 = p2*fx1211 + q2*fx2211
        fxy121 = p2*fx1121 + q2*fx2121
        fxy221 = p2*fx1221 + q2*fx2221
        fxy112 = p2*fx1112 + q2*fx2112
        fxy212 = p2*fx1212 + q2*fx2212
        fxy122 = p2*fx1122 + q2*fx2122
        fxy222 = p2*fx1222 + q2*fx2222

        fxyz11 = p3*fxy111 + q3*fxy211
        fxyz21 = p3*fxy121 + q3*fxy221
        fxyz12 = p3*fxy112 + q3*fxy212
        fxyz22 = p3*fxy122 + q3*fxy222

        fxyzq1 = p4*fxyz11 + q4*fxyz21
        fxyzq2 = p4*fxyz12 + q4*fxyz22

        f = p5*fxyzq1 + q5*fxyzq2
        if (present(istat)) istat = 0

    else

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

    end if

    end subroutine interp_5d