initialize_5d_specify_knots_multi Subroutine

    pure subroutine initialize_5d_specify_knots_multi(me,x,y,z,q,r,fcn,kx,ky,kz,kq,kr,tx,ty,tz,tq,tr,iflag,extrap)

    implicit none

    class(bspline_5d),intent(inout) :: me  !! the spline
    real(wp),dimension(:),intent(in) :: x  !! `(nx)` array of \(x\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in) :: y  !! `(ny)` array of \(y\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in) :: z  !! `(nz)` array of \(z\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in) :: q  !! `(nq)` array of \(q\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in) :: r  !! `(nr)` array of \(r\) abcissae. Must be strictly increasing.
    real(wp),dimension(:,:,:,:,:,:),intent(in) :: fcn  !! `(nx,ny,nz,nq,nr,nf)` array of function values to interpolate
                                                  !! (the last dimension is the output index)
    integer(ip),intent(in) :: kx  !! The order of spline pieces in \(x\)
    integer(ip),intent(in) :: ky  !! The order of spline pieces in \(y\)
    integer(ip),intent(in) :: kz  !! The order of spline pieces in \(z\)
    integer(ip),intent(in) :: kq  !! The order of spline pieces in \(q\)
    integer(ip),intent(in) :: kr  !! The order of spline pieces in \(r\)
    real(wp),dimension(:),intent(in) :: tx  !! The `(nx+kx)` knots in the \(x\) direction.
    real(wp),dimension(:),intent(in) :: ty  !! The `(ny+ky)` knots in the \(y\) direction.
    real(wp),dimension(:),intent(in) :: tz  !! The `(nz+kz)` knots in the \(z\) direction.
    real(wp),dimension(:),intent(in) :: tq  !! The `(nq+kq)` knots in the \(q\) direction.
    real(wp),dimension(:),intent(in) :: tr  !! The `(nr+kr)` knots in the \(r\) direction.
    integer(ip),intent(out) :: iflag  !! status flag (see [[db5ink]])
    logical,intent(in),optional :: extrap  !! if true, then extrapolation is allowed
                                           !! (default is false)

    integer(ip),dimension(5) :: n  !! number of points in each direction
    integer(ip) :: nf  !! number of outputs

    call me%destroy()

    n(1) = size(x,kind=ip)
    n(2) = size(y,kind=ip)
    n(3) = size(z,kind=ip)
    n(4) = size(q,kind=ip)
    n(5) = size(r,kind=ip)
    nf = size(fcn,6,kind=ip)

    call check_knot_vectors_sizes(nx=n(1),kx=kx,tx=tx,ny=n(2),ky=ky,ty=ty,nz=n(3),kz=kz,tz=tz,nq=n(4),kq=kq,tq=tq,nr=n(5),kr=kr,tr=tr,iflag=iflag)
    if (iflag/=0_ip) then
        me%iflag = iflag
        return
    end if

    me%nx = n(1)
    me%ny = n(2)
    me%nz = n(3)
    me%nq = n(4)
    me%nr = n(5)
    me%kx = kx
    me%ky = ky
    me%kz = kz
    me%kq = kq
    me%kr = kr
    allocate(me%tx(n(1)+kx))
    allocate(me%ty(n(2)+ky))
    allocate(me%tz(n(3)+kz))
    allocate(me%tq(n(4)+kq))
    allocate(me%tr(n(5)+kr))
    allocate(me%bcoef(n(1),n(2),n(3),n(4),n(5),nf))
    allocate(me%work_val_1(ky,kz,kq,kr))
    allocate(me%work_val_2(kz,kq,kr))
    allocate(me%work_val_3(kq,kr))
    allocate(me%work_val_4(kr))
    allocate(me%work_val_5(3_ip*max(kx,ky,kz,kq,kr)))

    me%tx = tx
    me%ty = ty
    me%tz = tz
    me%tq = tq
    me%tr = tr
    call db5ink(x,n(1),y,n(2),z,n(3),q,n(4),r,n(5),fcn,kx,ky,kz,kq,kr,1_ip,me%tx,me%ty,me%tz,me%tq,me%tr,me%bcoef,iflag)

    if (iflag==0_ip) then
        me%nf = nf
        call me%set_extrap_flag(extrap)
    end if

    me%initialized = iflag==0_ip
    me%iflag = iflag

    end subroutine initialize_5d_specify_knots_multi