interp_6d Subroutine

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

    implicit none

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

    integer,dimension(2) :: ix, iy, iz, iq, ir, is
    real(wp) :: p1, p2, p3, p4, p5, p6
    real(wp) :: q1, q2, q3, q4, q5, q6
    integer :: mflag
    real(wp) :: fx11111, fx21111, fx12111, fx22111, fx11211, fx21211, fx12211, &
                fx22211, fxy1111, fxy2111, fxy1211, fxy2211, fxyz111, fxyz211, &
                fxyzq11, fx11121, fx21121, fx12121, fx22121, fx11221, fx21221, &
                fx12221, fx22221, fxy1121, fxy2121, fxy1221, fxy2221, fxyz121, &
                fxyz221, fxyzq21, fx11112, fx21112, fx12112, fx22112, fx11212, &
                fx21212, fx12212, fx22212, fxy1112, fxy2112, fxy1212, fxy2212, &
                fxyz112, fxyz212, fxyzq12, fx11122, fx21122, fx12122, fx22122, &
                fx11222, fx21222, fx12222, fx22222, fxy1122, fxy2122, fxy1222, &
                fxy2222, fxyz122, fxyz222, fxyzq22, fxyzqr1, fxyzqr2

    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)
        call dintrv(me%s,s,me%ilos,is(1),is(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)))
        q6 = (s-me%s(is(1)))/(me%s(is(2))-me%s(is(1)))
        p1 = one-q1
        p2 = one-q2
        p3 = one-q3
        p4 = one-q4
        p5 = one-q5
        p6 = one-q6

        fx11111 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(1),is(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(1),is(1))
        fx21111 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(1),is(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(1),is(1))
        fx12111 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(1),is(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(1),is(1))
        fx22111 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(1),is(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(1),is(1))
        fx11211 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(1),is(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(1),is(1))
        fx21211 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(1),is(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(1),is(1))
        fx12211 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(1),is(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(1),is(1))
        fx22211 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(1),is(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(1),is(1))
        fx11121 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(2),is(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(2),is(1))
        fx21121 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(2),is(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(2),is(1))
        fx12121 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(2),is(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(2),is(1))
        fx22121 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(2),is(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(2),is(1))
        fx11221 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(2),is(1)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(2),is(1))
        fx21221 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(2),is(1)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(2),is(1))
        fx12221 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(2),is(1)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(2),is(1))
        fx22221 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(2),is(1)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(2),is(1))
        fx11112 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(1),is(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(1),is(2))
        fx21112 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(1),is(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(1),is(2))
        fx12112 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(1),is(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(1),is(2))
        fx22112 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(1),is(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(1),is(2))
        fx11212 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(1),is(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(1),is(2))
        fx21212 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(1),is(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(1),is(2))
        fx12212 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(1),is(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(1),is(2))
        fx22212 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(1),is(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(1),is(2))
        fx11122 = p1*me%f(ix(1),iy(1),iz(1),iq(1),ir(2),is(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(1),ir(2),is(2))
        fx21122 = p1*me%f(ix(1),iy(2),iz(1),iq(1),ir(2),is(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(1),ir(2),is(2))
        fx12122 = p1*me%f(ix(1),iy(1),iz(2),iq(1),ir(2),is(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(1),ir(2),is(2))
        fx22122 = p1*me%f(ix(1),iy(2),iz(2),iq(1),ir(2),is(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(1),ir(2),is(2))
        fx11222 = p1*me%f(ix(1),iy(1),iz(1),iq(2),ir(2),is(2)) + q1*me%f(ix(2),iy(1),iz(1),iq(2),ir(2),is(2))
        fx21222 = p1*me%f(ix(1),iy(2),iz(1),iq(2),ir(2),is(2)) + q1*me%f(ix(2),iy(2),iz(1),iq(2),ir(2),is(2))
        fx12222 = p1*me%f(ix(1),iy(1),iz(2),iq(2),ir(2),is(2)) + q1*me%f(ix(2),iy(1),iz(2),iq(2),ir(2),is(2))
        fx22222 = p1*me%f(ix(1),iy(2),iz(2),iq(2),ir(2),is(2)) + q1*me%f(ix(2),iy(2),iz(2),iq(2),ir(2),is(2))

        fxy1111 = p2*fx11111 + q2*fx21111
        fxy2111 = p2*fx12111 + q2*fx22111
        fxy1211 = p2*fx11211 + q2*fx21211
        fxy2211 = p2*fx12211 + q2*fx22211
        fxy1121 = p2*fx11121 + q2*fx21121
        fxy2121 = p2*fx12121 + q2*fx22121
        fxy1221 = p2*fx11221 + q2*fx21221
        fxy2221 = p2*fx12221 + q2*fx22221
        fxy1112 = p2*fx11112 + q2*fx21112
        fxy2112 = p2*fx12112 + q2*fx22112
        fxy1212 = p2*fx11212 + q2*fx21212
        fxy2212 = p2*fx12212 + q2*fx22212
        fxy1122 = p2*fx11122 + q2*fx21122
        fxy2122 = p2*fx12122 + q2*fx22122
        fxy1222 = p2*fx11222 + q2*fx21222
        fxy2222 = p2*fx12222 + q2*fx22222

        fxyz111 = p3*fxy1111 + q3*fxy2111
        fxyz211 = p3*fxy1211 + q3*fxy2211
        fxyz121 = p3*fxy1121 + q3*fxy2121
        fxyz221 = p3*fxy1221 + q3*fxy2221
        fxyz112 = p3*fxy1112 + q3*fxy2112
        fxyz212 = p3*fxy1212 + q3*fxy2212
        fxyz122 = p3*fxy1122 + q3*fxy2122
        fxyz222 = p3*fxy1222 + q3*fxy2222

        fxyzq11 = p4*fxyz111 + q4*fxyz211
        fxyzq21 = p4*fxyz121 + q4*fxyz221
        fxyzq12 = p4*fxyz112 + q4*fxyz212
        fxyzq22 = p4*fxyz122 + q4*fxyz222

        fxyzqr1 = p5*fxyzq11 + q5*fxyzq21
        fxyzqr2 = p5*fxyzq12 + q5*fxyzq22

        f = p6*fxyzqr1 + q6*fxyzqr2
        if (present(istat)) istat = 0

    else

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

    end if

    end subroutine interp_6d