interp_4d Subroutine

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

    implicit none

    class(linear_interp_4d),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(out)                  :: f     !! Interpolated \( f(x,y,z,q) \)
    integer,intent(out),optional          :: istat  !! `0`  : no problems,
                                                    !! `-1` : class has not been initialized

    integer,dimension(2) :: ix, iy, iz, iq
    real(wp) :: p1, p2, p3, p4
    real(wp) :: q1, q2, q3, q4
    integer :: mflag
    real(wp) :: fx111,fx211,fx121,fx221,fxy11,fxy21,fxyz1,&
                fx112,fx212,fx122,fx222,fxy12,fxy22,fxyz2

    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)

        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)))
        p1 = one-q1
        p2 = one-q2
        p3 = one-q3
        p4 = one-q4

        fx111 = p1*me%f(ix(1),iy(1),iz(1),iq(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(1))
        fx211 = p1*me%f(ix(1),iy(2),iz(1),iq(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(1))
        fx121 = p1*me%f(ix(1),iy(1),iz(2),iq(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(1))
        fx221 = p1*me%f(ix(1),iy(2),iz(2),iq(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(1))
        fx112 = p1*me%f(ix(1),iy(1),iz(1),iq(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(2))
        fx212 = p1*me%f(ix(1),iy(2),iz(1),iq(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(2))
        fx122 = p1*me%f(ix(1),iy(1),iz(2),iq(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(2))
        fx222 = p1*me%f(ix(1),iy(2),iz(2),iq(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(2))

        fxy11 = p2*fx111 + q2*fx211
        fxy21 = p2*fx121 + q2*fx221
        fxy12 = p2*fx112 + q2*fx212
        fxy22 = p2*fx122 + q2*fx222

        fxyz1 = p3*fxy11 + q3*fxy21
        fxyz2 = p3*fxy12 + q3*fxy22

        f = p4*fxyz1 + q4*fxyz2
        if (present(istat)) istat = 0

    else

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

    end if

    end subroutine interp_4d