bspline_oo_module.f90 Source File


This file depends on

sourcefile~~bspline_oo_module.f90~~EfferentGraph sourcefile~bspline_oo_module.f90 bspline_oo_module.f90 sourcefile~bspline_kinds_module.f90 bspline_kinds_module.F90 sourcefile~bspline_oo_module.f90->sourcefile~bspline_kinds_module.f90 sourcefile~bspline_sub_module.f90 bspline_sub_module.f90 sourcefile~bspline_oo_module.f90->sourcefile~bspline_sub_module.f90 sourcefile~bspline_sub_module.f90->sourcefile~bspline_kinds_module.f90

Files dependent on this one

sourcefile~~bspline_oo_module.f90~~AfferentGraph sourcefile~bspline_oo_module.f90 bspline_oo_module.f90 sourcefile~bspline_module.f90 bspline_module.f90 sourcefile~bspline_module.f90->sourcefile~bspline_oo_module.f90

Source Code

!*****************************************************************************************
!> author: Jacob Williams
!  license: BSD
!  date: 12/6/2015
!
!  Object-oriented style wrappers to [[bspline_sub_module]].
!  This module provides classes ([[bspline_1d(type)]], [[bspline_2d(type)]],
!  [[bspline_3d(type)]], [[bspline_4d(type)]], [[bspline_5d(type)]], and [[bspline_6d(type)]])
!  which can be used instead of the main subroutine interface.

    module bspline_oo_module

    use bspline_kinds_module, only: wp, ip
    use,intrinsic :: iso_fortran_env, only: error_unit
    use bspline_sub_module

    implicit none

    private

    integer(ip),parameter :: int_size     = storage_size(1_ip,kind=ip)   !! size of a default integer [bits]
    integer(ip),parameter :: logical_size = storage_size(.true.,kind=ip) !! size of a default logical [bits]
    integer(ip),parameter :: real_size    = storage_size(1.0_wp,kind=ip) !! size of a `real(wp)` [bits]

    type,public,abstract :: bspline_class
        !! Base class for the b-spline types
        private
        integer(ip) :: inbvx = 1_ip  !! internal variable used by [[dbvalu]] for efficient processing
        integer(ip) :: iflag = 1_ip  !! saved `iflag` from the list routine call.
        logical :: initialized = .false. !! true if the class is initialized and ready to use
        logical :: extrap = .false. !! if true, then extrapolation is allowed during evaluation
    contains
        private
        procedure,non_overridable :: destroy_base  !! destructor for the abstract type
        procedure,non_overridable :: set_extrap_flag !! internal routine to set the `extrap` flag
        procedure(destroy_func),deferred,public :: destroy  !! destructor
        procedure(size_func),deferred,public :: size_of !! size of the structure in bits
        procedure,public,non_overridable :: status_ok  !! returns true if the last `iflag` status code was `=0`.
        procedure,public,non_overridable :: status_message => get_bspline_status_message  !! retrieve the last
                                                                                          !! status message
        procedure,public,non_overridable :: clear_flag => clear_bspline_flag  !! to reset the `iflag` saved in the class.
    end type bspline_class

    abstract interface

        pure subroutine destroy_func(me)
        !! interface for bspline destructor routines
        import :: bspline_class
        implicit none
        class(bspline_class),intent(inout) :: me
        end subroutine destroy_func

        pure function size_func(me) result(s)
        !! interface for size routines
        import :: bspline_class,ip
        implicit none
        class(bspline_class),intent(in) :: me
        integer(ip) :: s !! size of the structure in bits
        end function size_func

    end interface

    type,extends(bspline_class),public :: bspline_1d
        !! Class for 1d b-spline interpolation.
        !!
        !!@note The 1D class also contains two methods
        !!      for computing definite integrals.
        private
        integer(ip) :: nx = 0_ip  !! Number of \(x\) abcissae
        integer(ip) :: kx = 0_ip  !! The order of spline pieces in \(x\)
        real(wp),dimension(:),allocatable :: bcoef  !! array of coefficients of the b-spline interpolant
        real(wp),dimension(:),allocatable :: tx  !! The knots in the \(x\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: work_val_1   !! [[db1val] work array of dimension `3*kx`
        contains
        private
        generic,public :: initialize => initialize_1d_auto_knots,initialize_1d_specify_knots
        procedure :: initialize_1d_auto_knots
        procedure :: initialize_1d_specify_knots
        procedure,public :: evaluate => evaluate_1d
        procedure,public :: destroy => destroy_1d
        procedure,public :: size_of => size_1d
        procedure,public :: integral => integral_1d
        procedure,public :: fintegral => fintegral_1d
        final :: finalize_1d
    end type bspline_1d

    type,extends(bspline_class),public :: bspline_2d
        !! Class for 2d b-spline interpolation.
        private
        integer(ip) :: nx = 0_ip  !! Number of \(x\) abcissae
        integer(ip) :: ny = 0_ip  !! Number of \(y\) abcissae
        integer(ip) :: kx = 0_ip  !! The order of spline pieces in \(x\)
        integer(ip) :: ky = 0_ip  !! The order of spline pieces in \(y\)
        real(wp),dimension(:,:),allocatable :: bcoef  !! array of coefficients of the b-spline interpolant
        real(wp),dimension(:),allocatable :: tx  !! The knots in the \(x\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: ty  !! The knots in the \(y\) direction for the spline interpolant
        integer(ip) :: inbvy = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloy = 1_ip  !! internal variable used for efficient processing
        real(wp),dimension(:),allocatable :: work_val_1  !! [[db2val] work array of dimension `ky`
        real(wp),dimension(:),allocatable :: work_val_2  !! [[db2val] work array of dimension `3_ip*max(kx,ky)`
        contains
        private
        generic,public :: initialize => initialize_2d_auto_knots,initialize_2d_specify_knots
        procedure :: initialize_2d_auto_knots
        procedure :: initialize_2d_specify_knots
        procedure,public :: evaluate => evaluate_2d
        procedure,public :: destroy => destroy_2d
        procedure,public :: size_of => size_2d
        final :: finalize_2d
    end type bspline_2d

    type,extends(bspline_class),public :: bspline_3d
        !! Class for 3d b-spline interpolation.
        private
        integer(ip) :: nx = 0_ip  !! Number of \(x\) abcissae
        integer(ip) :: ny = 0_ip  !! Number of \(y\) abcissae
        integer(ip) :: nz = 0_ip  !! Number of \(z\) abcissae
        integer(ip) :: kx = 0_ip  !! The order of spline pieces in \(x\)
        integer(ip) :: ky = 0_ip  !! The order of spline pieces in \(y\)
        integer(ip) :: kz = 0_ip  !! The order of spline pieces in \(z\)
        real(wp),dimension(:,:,:),allocatable :: bcoef  !! array of coefficients of the b-spline interpolant
        real(wp),dimension(:),allocatable :: tx  !! The knots in the \(x\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: ty  !! The knots in the \(y\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: tz  !! The knots in the \(z\) direction for the spline interpolant
        integer(ip) :: inbvy = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: inbvz = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloy = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloz = 1_ip  !! internal variable used for efficient processing
        real(wp),dimension(:,:),allocatable :: work_val_1  !! [[db3val] work array of dimension `ky,kz`
        real(wp),dimension(:),allocatable   :: work_val_2  !! [[db3val] work array of dimension `kz`
        real(wp),dimension(:),allocatable   :: work_val_3  !! [[db3val] work array of dimension `3_ip*max(kx,ky,kz)`
        contains
        private
        generic,public :: initialize => initialize_3d_auto_knots,initialize_3d_specify_knots
        procedure :: initialize_3d_auto_knots
        procedure :: initialize_3d_specify_knots
        procedure,public :: evaluate => evaluate_3d
        procedure,public :: destroy => destroy_3d
        procedure,public :: size_of => size_3d
        final :: finalize_3d
    end type bspline_3d

    type,extends(bspline_class),public :: bspline_4d
        !! Class for 4d b-spline interpolation.
        private
        integer(ip) :: nx = 0_ip  !! Number of \(x\) abcissae
        integer(ip) :: ny = 0_ip  !! Number of \(y\) abcissae
        integer(ip) :: nz = 0_ip  !! Number of \(z\) abcissae
        integer(ip) :: nq = 0_ip  !! Number of \(q\) abcissae
        integer(ip) :: kx = 0_ip  !! The order of spline pieces in \(x\)
        integer(ip) :: ky = 0_ip  !! The order of spline pieces in \(y\)
        integer(ip) :: kz = 0_ip  !! The order of spline pieces in \(z\)
        integer(ip) :: kq = 0_ip  !! The order of spline pieces in \(q\)
        real(wp),dimension(:,:,:,:),allocatable :: bcoef  !! array of coefficients of the b-spline interpolant
        real(wp),dimension(:),allocatable :: tx  !! The knots in the \(x\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: ty  !! The knots in the \(y\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: tz  !! The knots in the \(z\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: tq  !! The knots in the \(q\) direction for the spline interpolant
        integer(ip) :: inbvy = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: inbvz = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: inbvq = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloy  = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloz  = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloq  = 1_ip  !! internal variable used for efficient processing
        real(wp),dimension(:,:,:),allocatable :: work_val_1  !! [[db4val]] work array of dimension `ky,kz,kq`
        real(wp),dimension(:,:),allocatable   :: work_val_2  !! [[db4val]] work array of dimension `kz,kq`
        real(wp),dimension(:),allocatable     :: work_val_3  !! [[db4val]] work array of dimension `kq`
        real(wp),dimension(:),allocatable     :: work_val_4  !! [[db4val]] work array of dimension `3_ip*max(kx,ky,kz,kq)`
        contains
        private
        generic,public :: initialize => initialize_4d_auto_knots,initialize_4d_specify_knots
        procedure :: initialize_4d_auto_knots
        procedure :: initialize_4d_specify_knots
        procedure,public :: evaluate => evaluate_4d
        procedure,public :: destroy => destroy_4d
        procedure,public :: size_of => size_4d
        final :: finalize_4d
    end type bspline_4d

    type,extends(bspline_class),public :: bspline_5d
        !! Class for 5d b-spline interpolation.
        private
        integer(ip) :: nx = 0_ip  !! Number of \(x\) abcissae
        integer(ip) :: ny = 0_ip  !! Number of \(y\) abcissae
        integer(ip) :: nz = 0_ip  !! Number of \(z\) abcissae
        integer(ip) :: nq = 0_ip  !! Number of \(q\) abcissae
        integer(ip) :: nr = 0_ip  !! Number of \(r\) abcissae
        integer(ip) :: kx = 0_ip  !! The order of spline pieces in \(x\)
        integer(ip) :: ky = 0_ip  !! The order of spline pieces in \(y\)
        integer(ip) :: kz = 0_ip  !! The order of spline pieces in \(z\)
        integer(ip) :: kq = 0_ip  !! The order of spline pieces in \(q\)
        integer(ip) :: kr = 0_ip  !! The order of spline pieces in \(r\)
        real(wp),dimension(:,:,:,:,:),allocatable :: bcoef  !! array of coefficients of the b-spline interpolant
        real(wp),dimension(:),allocatable :: tx  !! The knots in the \(x\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: ty  !! The knots in the \(y\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: tz  !! The knots in the \(z\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: tq  !! The knots in the \(q\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: tr  !! The knots in the \(r\) direction for the spline interpolant
        integer(ip) :: inbvy = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: inbvz = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: inbvq = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: inbvr = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloy  = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloz  = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloq  = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: ilor  = 1_ip  !! internal variable used for efficient processing
        real(wp),dimension(:,:,:,:),allocatable :: work_val_1  !! [[db5val]] work array of dimension `ky,kz,kq,kr`
        real(wp),dimension(:,:,:),allocatable   :: work_val_2  !! [[db5val]] work array of dimension `kz,kq,kr`
        real(wp),dimension(:,:),allocatable     :: work_val_3  !! [[db5val]] work array of dimension `kq,kr`
        real(wp),dimension(:),allocatable       :: work_val_4  !! [[db5val]] work array of dimension `kr`
        real(wp),dimension(:),allocatable       :: work_val_5  !! [[db5val]] work array of dimension `3_ip*max(kx,ky,kz,kq,kr)`
        contains
        private
        generic,public :: initialize => initialize_5d_auto_knots,initialize_5d_specify_knots
        procedure :: initialize_5d_auto_knots
        procedure :: initialize_5d_specify_knots
        procedure,public :: evaluate => evaluate_5d
        procedure,public :: destroy => destroy_5d
        procedure,public :: size_of => size_5d
        final :: finalize_5d
    end type bspline_5d

    type,extends(bspline_class),public :: bspline_6d
        !! Class for 6d b-spline interpolation.
        private
        integer(ip) :: nx = 0_ip  !! Number of \(x\) abcissae
        integer(ip) :: ny = 0_ip  !! Number of \(y\) abcissae
        integer(ip) :: nz = 0_ip  !! Number of \(z\) abcissae
        integer(ip) :: nq = 0_ip  !! Number of \(q\) abcissae
        integer(ip) :: nr = 0_ip  !! Number of \(r\) abcissae
        integer(ip) :: ns = 0_ip  !! Number of \(s\) abcissae
        integer(ip) :: kx = 0_ip  !! The order of spline pieces in \(x\)
        integer(ip) :: ky = 0_ip  !! The order of spline pieces in \(y\)
        integer(ip) :: kz = 0_ip  !! The order of spline pieces in \(z\)
        integer(ip) :: kq = 0_ip  !! The order of spline pieces in \(q\)
        integer(ip) :: kr = 0_ip  !! The order of spline pieces in \(r\)
        integer(ip) :: ks = 0_ip  !! The order of spline pieces in \(s\)
        real(wp),dimension(:,:,:,:,:,:),allocatable :: bcoef  !! array of coefficients of the b-spline interpolant
        real(wp),dimension(:),allocatable :: tx  !! The knots in the \(x\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: ty  !! The knots in the \(y\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: tz  !! The knots in the \(z\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: tq  !! The knots in the \(q\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: tr  !! The knots in the \(r\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: ts  !! The knots in the \(s\) direction for the spline interpolant
        integer(ip) :: inbvy = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: inbvz = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: inbvq = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: inbvr = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: inbvs = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloy  = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloz  = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloq  = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: ilor  = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: ilos  = 1_ip  !! internal variable used for efficient processing
        real(wp),dimension(:,:,:,:,:),allocatable :: work_val_1  !! [[db6val]] work array of dimension `ky,kz,kq,kr,ks`
        real(wp),dimension(:,:,:,:),allocatable   :: work_val_2  !! [[db6val]] work array of dimension `kz,kq,kr,ks`
        real(wp),dimension(:,:,:),allocatable     :: work_val_3  !! [[db6val]] work array of dimension `kq,kr,ks`
        real(wp),dimension(:,:),allocatable       :: work_val_4  !! [[db6val]] work array of dimension `kr,ks`
        real(wp),dimension(:),allocatable         :: work_val_5  !! [[db6val]] work array of dimension `ks`
        real(wp),dimension(:),allocatable         :: work_val_6  !! [[db6val]] work array of dimension `3_ip*max(kx,ky,kz,kq,kr,ks)`
        contains
        private
        generic,public :: initialize => initialize_6d_auto_knots,initialize_6d_specify_knots
        procedure :: initialize_6d_auto_knots
        procedure :: initialize_6d_specify_knots
        procedure,public :: evaluate => evaluate_6d
        procedure,public :: destroy => destroy_6d
        procedure,public :: size_of => size_6d
        final :: finalize_6d
    end type bspline_6d

    interface bspline_1d
        !! Constructor for [[bspline_1d(type)]]
        procedure :: bspline_1d_constructor_empty,&
                     bspline_1d_constructor_auto_knots,&
                     bspline_1d_constructor_specify_knots
    end interface
    interface bspline_2d
        !! Constructor for [[bspline_2d(type)]]
        procedure :: bspline_2d_constructor_empty,&
                     bspline_2d_constructor_auto_knots,&
                     bspline_2d_constructor_specify_knots
    end interface
    interface bspline_3d
        !! Constructor for [[bspline_3d(type)]]
        procedure :: bspline_3d_constructor_empty,&
                     bspline_3d_constructor_auto_knots,&
                     bspline_3d_constructor_specify_knots
    end interface
    interface bspline_4d
        !! Constructor for [[bspline_4d(type)]]
        procedure :: bspline_4d_constructor_empty,&
                     bspline_4d_constructor_auto_knots,&
                     bspline_4d_constructor_specify_knots
    end interface
    interface bspline_5d
        !! Constructor for [[bspline_5d(type)]]
        procedure :: bspline_5d_constructor_empty,&
                     bspline_5d_constructor_auto_knots,&
                     bspline_5d_constructor_specify_knots
    end interface
    interface bspline_6d
        !! Constructor for [[bspline_6d(type)]]
        procedure :: bspline_6d_constructor_empty,&
                     bspline_6d_constructor_auto_knots,&
                     bspline_6d_constructor_specify_knots
    end interface

    contains
!*****************************************************************************************

!*****************************************************************************************
!>
!  This routines returns true if the `iflag` code from the last
!  routine called was `=0`. Maybe of the routines have output `iflag`
!  variables, so they can be checked explicitly, or this routine
!  can be used.
!
!  If the class is initialized using a function constructor, then
!  this is the only way to know if it was properly initialized,
!  since those are pure functions with not output `iflag` arguments.
!
!  If `status_ok=.false.`, then the error message can be
!  obtained from the [[get_bspline_status_message]] routine.
!
!  Note: after an error condition, the [[clear_bspline_flag]] routine
!  can be called to reset the `iflag` to 0.

    elemental function status_ok(me) result(ok)

    implicit none

    class(bspline_class),intent(in) :: me
    logical                         :: ok

    ok = ( me%iflag == 0_ip )

    end function status_ok
!*****************************************************************************************

!*****************************************************************************************
!>
!  This sets the `iflag` variable in the class to `0`
!  (which indicates that everything is OK). It can be used
!  after an error is encountered.

    elemental subroutine clear_bspline_flag(me)

    implicit none

    class(bspline_class),intent(inout) :: me

    me%iflag = 0_ip

    end subroutine clear_bspline_flag
!*****************************************************************************************

!*****************************************************************************************
!>
!  Get the status message from a [[bspline_class]] routine call.
!
!  If `iflag` is not included, then the one in the class is used (which
!  corresponds to the last routine called.)
!  Otherwise, it will convert the
!  input `iflag` argument into the appropriate message.
!
!  This is a wrapper for [[get_status_message]].

    pure function get_bspline_status_message(me,iflag) result(msg)

    implicit none

    class(bspline_class),intent(in) :: me
    character(len=:),allocatable    :: msg    !! status message associated with the flag
    integer(ip),intent(in),optional :: iflag  !! the corresponding status code

    if (present(iflag)) then
        msg = get_status_message(iflag)
    else
        msg = get_status_message(me%iflag)
    end if

    end function get_bspline_status_message
!*****************************************************************************************

!*****************************************************************************************
!>
!  Actual size of a [[bspline_1d]] structure in bits.

    pure function size_1d(me) result(s)

    implicit none

    class(bspline_1d),intent(in) :: me
    integer(ip) :: s !! size of the structure in bits

    s = 2_ip*int_size + logical_size + 2_ip*int_size

    if (allocated(me%bcoef))      s = s + real_size*size(me%bcoef,kind=ip)
    if (allocated(me%tx))         s = s + real_size*size(me%tx,kind=ip)
    if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,kind=ip)

    end function size_1d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Actual size of a [[bspline_2d]] structure in bits.

    pure function size_2d(me) result(s)

    implicit none

    class(bspline_2d),intent(in) :: me
    integer(ip) :: s !! size of the structure in bits

    s = 2_ip*int_size + logical_size + 6_ip*int_size

    if (allocated(me%bcoef))      s = s + real_size*size(me%bcoef,1_ip,kind=ip)*&
                                                    size(me%bcoef,2_ip,kind=ip)
    if (allocated(me%tx))         s = s + real_size*size(me%tx,kind=ip)
    if (allocated(me%ty))         s = s + real_size*size(me%ty,kind=ip)
    if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,kind=ip)
    if (allocated(me%work_val_2)) s = s + real_size*size(me%work_val_2,kind=ip)

    end function size_2d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Actual size of a [[bspline_3d]] structure in bits.

    pure function size_3d(me) result(s)

    implicit none

    class(bspline_3d),intent(in) :: me
    integer(ip) :: s !! size of the structure in bits

    s = 2_ip*int_size + logical_size + 10_ip*int_size

    if (allocated(me%bcoef))      s = s + real_size*size(me%bcoef,1_ip,kind=ip)*&
                                                    size(me%bcoef,2_ip,kind=ip)*&
                                                    size(me%bcoef,3_ip,kind=ip)
    if (allocated(me%tx))         s = s + real_size*size(me%tx,kind=ip)
    if (allocated(me%ty))         s = s + real_size*size(me%ty,kind=ip)
    if (allocated(me%tz))         s = s + real_size*size(me%tz,kind=ip)
    if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,1_ip,kind=ip)*&
                                                    size(me%work_val_1,2_ip,kind=ip)
    if (allocated(me%work_val_2)) s = s + real_size*size(me%work_val_2,kind=ip)
    if (allocated(me%work_val_3)) s = s + real_size*size(me%work_val_3,kind=ip)

    end function size_3d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Actual size of a [[bspline_4d]] structure in bits.

    pure function size_4d(me) result(s)

    implicit none

    class(bspline_4d),intent(in) :: me
    integer(ip) :: s !! size of the structure in bits

    s = 2_ip*int_size + logical_size + 14_ip*int_size

    if (allocated(me%bcoef))      s = s + real_size*size(me%bcoef,1_ip,kind=ip)*&
                                                    size(me%bcoef,2_ip,kind=ip)*&
                                                    size(me%bcoef,3_ip,kind=ip)*&
                                                    size(me%bcoef,4_ip,kind=ip)
    if (allocated(me%tx))         s = s + real_size*size(me%tx,kind=ip)
    if (allocated(me%ty))         s = s + real_size*size(me%ty,kind=ip)
    if (allocated(me%tz))         s = s + real_size*size(me%tz,kind=ip)
    if (allocated(me%tq))         s = s + real_size*size(me%tq,kind=ip)
    if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,1_ip,kind=ip)*&
                                                    size(me%work_val_1,2_ip,kind=ip)*&
                                                    size(me%work_val_1,3_ip,kind=ip)
    if (allocated(me%work_val_2)) s = s + real_size*size(me%work_val_2,1_ip,kind=ip)*&
                                                    size(me%work_val_2,2_ip,kind=ip)
    if (allocated(me%work_val_3)) s = s + real_size*size(me%work_val_3,kind=ip)
    if (allocated(me%work_val_4)) s = s + real_size*size(me%work_val_4,kind=ip)

    end function size_4d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Actual size of a [[bspline_5d]] structure in bits.

    pure function size_5d(me) result(s)

    implicit none

    class(bspline_5d),intent(in) :: me
    integer(ip) :: s !! size of the structure in bits

    s = 2_ip*int_size + logical_size + 18_ip*int_size

    if (allocated(me%bcoef))      s = s + real_size*size(me%bcoef,1_ip,kind=ip)*&
                                                    size(me%bcoef,2_ip,kind=ip)*&
                                                    size(me%bcoef,3_ip,kind=ip)*&
                                                    size(me%bcoef,4_ip,kind=ip)*&
                                                    size(me%bcoef,5_ip,kind=ip)
    if (allocated(me%tx))         s = s + real_size*size(me%tx,kind=ip)
    if (allocated(me%ty))         s = s + real_size*size(me%ty,kind=ip)
    if (allocated(me%tz))         s = s + real_size*size(me%tz,kind=ip)
    if (allocated(me%tq))         s = s + real_size*size(me%tq,kind=ip)
    if (allocated(me%tr))         s = s + real_size*size(me%tr,kind=ip)
    if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,1_ip,kind=ip)*&
                                                    size(me%work_val_1,2_ip,kind=ip)*&
                                                    size(me%work_val_1,3_ip,kind=ip)*&
                                                    size(me%work_val_1,4_ip,kind=ip)
    if (allocated(me%work_val_2)) s = s + real_size*size(me%work_val_2,1_ip,kind=ip)*&
                                                    size(me%work_val_2,2_ip,kind=ip)*&
                                                    size(me%work_val_2,3_ip,kind=ip)
    if (allocated(me%work_val_3)) s = s + real_size*size(me%work_val_3,1_ip,kind=ip)*&
                                                    size(me%work_val_3,2_ip,kind=ip)
    if (allocated(me%work_val_4)) s = s + real_size*size(me%work_val_4,kind=ip)
    if (allocated(me%work_val_5)) s = s + real_size*size(me%work_val_5,kind=ip)

    end function size_5d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Actual size of a [[bspline_6d]] structure in bits.

    pure function size_6d(me) result(s)

    implicit none

    class(bspline_6d),intent(in) :: me
    integer(ip) :: s !! size of the structure in bits

    s = 2_ip*int_size + logical_size + 22_ip*int_size

    if (allocated(me%bcoef))      s = s + real_size*size(me%bcoef,1_ip,kind=ip)*&
                                                    size(me%bcoef,2_ip,kind=ip)*&
                                                    size(me%bcoef,3_ip,kind=ip)*&
                                                    size(me%bcoef,4_ip,kind=ip)*&
                                                    size(me%bcoef,5_ip,kind=ip)*&
                                                    size(me%bcoef,6,kind=ip)
    if (allocated(me%tx))         s = s + real_size*size(me%tx,kind=ip)
    if (allocated(me%ty))         s = s + real_size*size(me%ty,kind=ip)
    if (allocated(me%tz))         s = s + real_size*size(me%tz,kind=ip)
    if (allocated(me%tq))         s = s + real_size*size(me%tq,kind=ip)
    if (allocated(me%tr))         s = s + real_size*size(me%tr,kind=ip)
    if (allocated(me%ts))         s = s + real_size*size(me%ts,kind=ip)
    if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,1_ip,kind=ip)*&
                                                    size(me%work_val_1,2_ip,kind=ip)*&
                                                    size(me%work_val_1,3_ip,kind=ip)*&
                                                    size(me%work_val_1,4_ip,kind=ip)*&
                                                    size(me%work_val_1,5_ip,kind=ip)
    if (allocated(me%work_val_2)) s = s + real_size*size(me%work_val_2,1_ip,kind=ip)*&
                                                    size(me%work_val_2,2_ip,kind=ip)*&
                                                    size(me%work_val_2,3_ip,kind=ip)*&
                                                    size(me%work_val_2,4_ip,kind=ip)
    if (allocated(me%work_val_3)) s = s + real_size*size(me%work_val_3,1_ip,kind=ip)*&
                                                    size(me%work_val_3,2_ip,kind=ip)*&
                                                    size(me%work_val_3,3_ip,kind=ip)
    if (allocated(me%work_val_4)) s = s + real_size*size(me%work_val_4,1_ip,kind=ip)*&
                                                    size(me%work_val_4,2_ip,kind=ip)
    if (allocated(me%work_val_5)) s = s + real_size*size(me%work_val_5,kind=ip)
    if (allocated(me%work_val_6)) s = s + real_size*size(me%work_val_6,kind=ip)

    end function size_6d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Destructor for contents of the base [[bspline_class]] class.
!  (this routine is called by the extended classes).

    pure subroutine destroy_base(me)

    implicit none

    class(bspline_class),intent(inout) :: me

    me%inbvx = 1_ip
    me%iflag = 1_ip
    me%initialized = .false.
    me%extrap = .false.

    end subroutine destroy_base
!*****************************************************************************************

!*****************************************************************************************
!>
!  Destructor for [[bspline_1d]] class.

    pure subroutine destroy_1d(me)

    implicit none

    class(bspline_1d),intent(inout) :: me

    call me%destroy_base()

    me%nx = 0_ip
    me%kx = 0_ip
    if (allocated(me%bcoef))      deallocate(me%bcoef)
    if (allocated(me%tx))         deallocate(me%tx)
    if (allocated(me%work_val_1)) deallocate(me%work_val_1)

    end subroutine destroy_1d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Destructor for [[bspline_2d]] class.

    pure subroutine destroy_2d(me)

    implicit none

    class(bspline_2d),intent(inout) :: me

    call me%destroy_base()

    me%nx    = 0_ip
    me%ny    = 0_ip
    me%kx    = 0_ip
    me%ky    = 0_ip
    me%inbvy = 1_ip
    me%iloy  = 1_ip
    if (allocated(me%bcoef))      deallocate(me%bcoef)
    if (allocated(me%tx))         deallocate(me%tx)
    if (allocated(me%ty))         deallocate(me%ty)
    if (allocated(me%work_val_1)) deallocate(me%work_val_1)
    if (allocated(me%work_val_2)) deallocate(me%work_val_2)

    end subroutine destroy_2d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Destructor for [[bspline_3d]] class.

    pure subroutine destroy_3d(me)

    implicit none

    class(bspline_3d),intent(inout) :: me

    call me%destroy_base()

    me%nx    = 0_ip
    me%ny    = 0_ip
    me%nz    = 0_ip
    me%kx    = 0_ip
    me%ky    = 0_ip
    me%kz    = 0_ip
    me%inbvy = 1_ip
    me%inbvz = 1_ip
    me%iloy  = 1_ip
    me%iloz  = 1_ip
    if (allocated(me%bcoef))      deallocate(me%bcoef)
    if (allocated(me%tx))         deallocate(me%tx)
    if (allocated(me%ty))         deallocate(me%ty)
    if (allocated(me%tz))         deallocate(me%tz)
    if (allocated(me%work_val_1)) deallocate(me%work_val_1)
    if (allocated(me%work_val_2)) deallocate(me%work_val_2)
    if (allocated(me%work_val_3)) deallocate(me%work_val_3)

    end subroutine destroy_3d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Destructor for [[bspline_4d]] class.

    pure subroutine destroy_4d(me)

    implicit none

    class(bspline_4d),intent(inout) :: me

    me%nx    = 0_ip
    me%ny    = 0_ip
    me%nz    = 0_ip
    me%nq    = 0_ip
    me%kx    = 0_ip
    me%ky    = 0_ip
    me%kz    = 0_ip
    me%kq    = 0_ip
    me%inbvy = 1_ip
    me%inbvz = 1_ip
    me%inbvq = 1_ip
    me%iloy  = 1_ip
    me%iloz  = 1_ip
    me%iloq  = 1_ip
    if (allocated(me%bcoef))      deallocate(me%bcoef)
    if (allocated(me%tx))         deallocate(me%tx)
    if (allocated(me%ty))         deallocate(me%ty)
    if (allocated(me%tz))         deallocate(me%tz)
    if (allocated(me%tq))         deallocate(me%tq)
    if (allocated(me%work_val_1)) deallocate(me%work_val_1)
    if (allocated(me%work_val_2)) deallocate(me%work_val_2)
    if (allocated(me%work_val_3)) deallocate(me%work_val_3)
    if (allocated(me%work_val_4)) deallocate(me%work_val_4)

    end subroutine destroy_4d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Destructor for [[bspline_5d]] class.

    pure subroutine destroy_5d(me)

    implicit none

    class(bspline_5d),intent(inout) :: me

    me%nx    = 0_ip
    me%ny    = 0_ip
    me%nz    = 0_ip
    me%nq    = 0_ip
    me%nr    = 0_ip
    me%kx    = 0_ip
    me%ky    = 0_ip
    me%kz    = 0_ip
    me%kq    = 0_ip
    me%kr    = 0_ip
    me%inbvy = 1_ip
    me%inbvz = 1_ip
    me%inbvq = 1_ip
    me%inbvr = 1_ip
    me%iloy  = 1_ip
    me%iloz  = 1_ip
    me%iloq  = 1_ip
    me%ilor  = 1_ip
    if (allocated(me%bcoef))      deallocate(me%bcoef)
    if (allocated(me%tx))         deallocate(me%tx)
    if (allocated(me%ty))         deallocate(me%ty)
    if (allocated(me%tz))         deallocate(me%tz)
    if (allocated(me%tq))         deallocate(me%tq)
    if (allocated(me%tr))         deallocate(me%tr)
    if (allocated(me%work_val_1)) deallocate(me%work_val_1)
    if (allocated(me%work_val_2)) deallocate(me%work_val_2)
    if (allocated(me%work_val_3)) deallocate(me%work_val_3)
    if (allocated(me%work_val_4)) deallocate(me%work_val_4)
    if (allocated(me%work_val_5)) deallocate(me%work_val_5)

    end subroutine destroy_5d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Destructor for [[bspline_6d]] class.

    pure subroutine destroy_6d(me)

    implicit none

    class(bspline_6d),intent(inout) :: me

    me%nx    = 0_ip
    me%ny    = 0_ip
    me%nz    = 0_ip
    me%nq    = 0_ip
    me%nr    = 0_ip
    me%ns    = 0_ip
    me%kx    = 0_ip
    me%ky    = 0_ip
    me%kz    = 0_ip
    me%kq    = 0_ip
    me%kr    = 0_ip
    me%ks    = 0_ip
    me%inbvy = 1_ip
    me%inbvz = 1_ip
    me%inbvq = 1_ip
    me%inbvr = 1_ip
    me%inbvs = 1_ip
    me%iloy  = 1_ip
    me%iloz  = 1_ip
    me%iloq  = 1_ip
    me%ilor  = 1_ip
    me%ilos  = 1_ip
    if (allocated(me%bcoef))      deallocate(me%bcoef)
    if (allocated(me%tx))         deallocate(me%tx)
    if (allocated(me%ty))         deallocate(me%ty)
    if (allocated(me%tz))         deallocate(me%tz)
    if (allocated(me%tq))         deallocate(me%tq)
    if (allocated(me%tr))         deallocate(me%tr)
    if (allocated(me%ts))         deallocate(me%ts)
    if (allocated(me%work_val_1)) deallocate(me%work_val_1)
    if (allocated(me%work_val_2)) deallocate(me%work_val_2)
    if (allocated(me%work_val_3)) deallocate(me%work_val_3)
    if (allocated(me%work_val_4)) deallocate(me%work_val_4)
    if (allocated(me%work_val_5)) deallocate(me%work_val_5)
    if (allocated(me%work_val_6)) deallocate(me%work_val_6)

    end subroutine destroy_6d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Finalizer for [[bspline_1d]] class. Just a wrapper for [[destroy_1d]].
    pure elemental subroutine finalize_1d(me)
        type(bspline_1d),intent(inout) :: me; call me%destroy()
    end subroutine finalize_1d
!*****************************************************************************************
!*****************************************************************************************
!>
!  Finalizer for [[bspline_2d]] class. Just a wrapper for [[destroy_2d]].
    pure elemental subroutine finalize_2d(me)
        type(bspline_2d),intent(inout) :: me; call me%destroy()
    end subroutine finalize_2d
!*****************************************************************************************
!*****************************************************************************************
!>
!  Finalizer for [[bspline_3d]] class. Just a wrapper for [[destroy_3d]].
    pure elemental subroutine finalize_3d(me)
        type(bspline_3d),intent(inout) :: me; call me%destroy()
    end subroutine finalize_3d
!*****************************************************************************************
!*****************************************************************************************
!>
!  Finalizer for [[bspline_4d]] class. Just a wrapper for [[destroy_4d]].
    pure elemental subroutine finalize_4d(me)
        type(bspline_4d),intent(inout) :: me; call me%destroy()
    end subroutine finalize_4d
!*****************************************************************************************
!*****************************************************************************************
!>
!  Finalizer for [[bspline_5d]] class. Just a wrapper for [[destroy_5d]].
    pure elemental subroutine finalize_5d(me)
        type(bspline_5d),intent(inout) :: me; call me%destroy()
    end subroutine finalize_5d
!*****************************************************************************************
!*****************************************************************************************
!>
!  Finalizer for [[bspline_6d]] class. Just a wrapper for [[destroy_6d]].
    pure elemental subroutine finalize_6d(me)
        type(bspline_6d),intent(inout) :: me; call me%destroy()
    end subroutine finalize_6d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Sets the `extrap` flag in the class.

    pure subroutine set_extrap_flag(me,extrap)

    implicit none

    class(bspline_class),intent(inout) :: me
    logical,intent(in),optional :: extrap  !! if not present, then False is used

    if (present(extrap)) then
        me%extrap = extrap
    else
        me%extrap = .false.
    end if

    end subroutine set_extrap_flag
!*****************************************************************************************

!*****************************************************************************************
!>
!  It returns an empty [[bspline_1d]] type. Note that INITIALIZE still
!  needs to be called before it can be used.
!  Not really that useful except perhaps in some OpenMP applications.

    pure elemental function bspline_1d_constructor_empty() result(me)

    implicit none

    type(bspline_1d) :: me

    end function bspline_1d_constructor_empty
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_1d]] type (auto knots).
!  This is a wrapper for [[initialize_1d_auto_knots]].

    pure function bspline_1d_constructor_auto_knots(x,fcn,kx,extrap) result(me)

    implicit none

    type(bspline_1d)                 :: me
    real(wp),dimension(:),intent(in) :: x      !! `(nx)` array of \(x\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in) :: fcn    !! `(nx)` array of function values to interpolate. `fcn(i)` should
                                               !! contain the function value at the point `x(i)`
    integer(ip),intent(in)           :: kx     !! The order of spline pieces in \(x\)
                                               !! ( \( 2 \le k_x < n_x \) )
                                               !! (order = polynomial degree + 1)
    logical,intent(in),optional      :: extrap !! if true, then extrapolation is allowed
                                               !! (default is false)

    call initialize_1d_auto_knots(me,x,fcn,kx,me%iflag,extrap)

    end function bspline_1d_constructor_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_1d]] type (user-specified knots).
!  This is a wrapper for [[initialize_1d_specify_knots]].

    pure function bspline_1d_constructor_specify_knots(x,fcn,kx,tx,extrap) result(me)

    implicit none

    type(bspline_1d)                 :: me
    real(wp),dimension(:),intent(in) :: x     !! `(nx)` array of \(x\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in) :: fcn   !! `(nx)` array of function values to interpolate. `fcn(i)` should
                                              !! contain the function value at the point `x(i)`
    integer(ip),intent(in)           :: kx    !! The order of spline pieces in \(x\)
                                              !! ( \( 2 \le k_x < n_x \) )
                                              !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in) :: tx    !! The `(nx+kx)` knots in the \(x\) direction
                                              !! for the spline interpolant.
                                              !! Must be non-decreasing.
    logical,intent(in),optional      :: extrap !! if true, then extrapolation is allowed
                                               !! (default is false)

    call initialize_1d_specify_knots(me,x,fcn,kx,tx,me%iflag,extrap)

    end function bspline_1d_constructor_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_1d]] type (with automatically-computed knots).
!  This is a wrapper for [[db1ink]].

    pure subroutine initialize_1d_auto_knots(me,x,fcn,kx,iflag,extrap)

    implicit none

    class(bspline_1d),intent(inout)  :: me
    real(wp),dimension(:),intent(in) :: x     !! `(nx)` array of \(x\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in) :: fcn   !! `(nx)` array of function values to interpolate. `fcn(i)` should
                                              !! contain the function value at the point `x(i)`
    integer(ip),intent(in)           :: kx    !! The order of spline pieces in \(x\)
                                              !! ( \( 2 \le k_x < n_x \) )
                                              !! (order = polynomial degree + 1)
    integer(ip),intent(out)          :: iflag !! status flag (see [[db1ink]])
    logical,intent(in),optional      :: extrap !! if true, then extrapolation is allowed
                                               !! (default is false)

    integer(ip) :: iknot
    integer(ip) :: nx

    call me%destroy()

    nx = size(x,kind=ip)

    me%nx = nx
    me%kx = kx

    allocate(me%tx(nx+kx))
    allocate(me%bcoef(nx))
    allocate(me%work_val_1(3_ip*kx))

    iknot = 0_ip         !knot sequence chosen by db1ink

    call db1ink(x,nx,fcn,kx,iknot,me%tx,me%bcoef,iflag)

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

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

    end subroutine initialize_1d_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_1d]] type (with user-specified knots).
!  This is a wrapper for [[db1ink]].

    pure subroutine initialize_1d_specify_knots(me,x,fcn,kx,tx,iflag,extrap)

    implicit none

    class(bspline_1d),intent(inout)  :: me
    real(wp),dimension(:),intent(in) :: x     !! `(nx)` array of \(x\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in) :: fcn   !! `(nx)` array of function values to interpolate. `fcn(i)` should
                                              !! contain the function value at the point `x(i)`
    integer(ip),intent(in)           :: kx    !! The order of spline pieces in \(x\)
                                              !! ( \( 2 \le k_x < n_x \) )
                                              !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in) :: tx    !! The `(nx+kx)` knots in the \(x\) direction
                                              !! for the spline interpolant.
                                              !! Must be non-decreasing.
    integer(ip),intent(out)          :: iflag !! status flag (see [[db1ink]])
    logical,intent(in),optional      :: extrap !! if true, then extrapolation is allowed
                                               !! (default is false)

    integer(ip) :: nx

    call me%destroy()

    nx = size(x,kind=ip)

    call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,iflag=iflag)

    if (iflag == 0_ip) then

        me%nx = nx
        me%kx = kx

        allocate(me%tx(nx+kx))
        allocate(me%bcoef(nx))
        allocate(me%work_val_1(3_ip*kx))

        me%tx = tx

        call db1ink(x,nx,fcn,kx,1_ip,me%tx,me%bcoef,iflag)

        call me%set_extrap_flag(extrap)

    end if

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

    end subroutine initialize_1d_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Evaluate a [[bspline_1d]] interpolate.  This is a wrapper for [[db1val]].

    pure subroutine evaluate_1d(me,xval,idx,f,iflag)

    implicit none

    class(bspline_1d),intent(inout) :: me
    real(wp),intent(in)             :: xval  !! \(x\) coordinate of evaluation point.
    integer(ip),intent(in)          :: idx   !! \(x\) derivative of piecewise polynomial to evaluate.
    real(wp),intent(out)            :: f     !! interpolated value
    integer(ip),intent(out)         :: iflag !! status flag (see [[db1val]])

    if (me%initialized) then
        call db1val(xval,idx,me%tx,me%nx,me%kx,me%bcoef,f,iflag,&
                    me%inbvx,me%work_val_1,extrap=me%extrap)
    else
        iflag = 1_ip
    end if
    me%iflag = iflag

    end subroutine evaluate_1d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Evaluate a [[bspline_1d]] definite integral.  This is a wrapper for [[db1sqad]].

    pure subroutine integral_1d(me,x1,x2,f,iflag)

    implicit none

    class(bspline_1d),intent(inout) :: me
    real(wp),intent(in)             :: x1    !! left point of interval
    real(wp),intent(in)             :: x2    !! right point of interval
    real(wp),intent(out)            :: f     !! integral of the b-spline over \( [x_1, x_2] \)
    integer(ip),intent(out)         :: iflag !! status flag (see [[db1sqad]])

    if (me%initialized) then
        call db1sqad(me%tx,me%bcoef,me%nx,me%kx,x1,x2,f,iflag,me%work_val_1)
    else
        iflag = 1_ip
    end if
    me%iflag = iflag

    end subroutine integral_1d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Evaluate a [[bspline_1d]] definite integral.  This is a wrapper for [[db1fqad]].

    subroutine fintegral_1d(me,fun,idx,x1,x2,tol,f,iflag)

    implicit none

    class(bspline_1d),intent(inout) :: me
    procedure(b1fqad_func)          :: fun   !! external function of one argument for the
                                             !! integrand `bf(x)=fun(x)*dbvalu(tx,bcoef,nx,kx,idx,x,inbv)`
    integer(ip),intent(in)          :: idx   !! order of the spline derivative, `0 <= idx <= k-1`
                                             !! `idx=0` gives the spline function
    real(wp),intent(in)             :: x1    !! left point of interval
    real(wp),intent(in)             :: x2    !! right point of interval
    real(wp),intent(in)             :: tol   !! desired accuracy for the quadrature
    real(wp),intent(out)            :: f     !! integral of `bf(x)` over \( [x_1, x_2] \)
    integer(ip),intent(out)         :: iflag !! status flag (see [[db1sqad]])

    if (me%initialized) then
        call db1fqad(fun,me%tx,me%bcoef,me%nx,me%kx,idx,x1,x2,tol,f,iflag,me%work_val_1)
    else
        iflag = 1_ip
    end if
    me%iflag = iflag

    end subroutine fintegral_1d
!*****************************************************************************************

!*****************************************************************************************
!>
!  It returns an empty [[bspline_2d]] type. Note that INITIALIZE still
!  needs to be called before it can be used.
!  Not really that useful except perhaps in some OpenMP applications.

    elemental function bspline_2d_constructor_empty() result(me)

    implicit none

    type(bspline_2d) :: me

    end function bspline_2d_constructor_empty
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_2d]] type (auto knots).
!  This is a wrapper for [[initialize_2d_auto_knots]].

    pure function bspline_2d_constructor_auto_knots(x,y,fcn,kx,ky,extrap) result(me)

    implicit none

    type(bspline_2d)                   :: me
    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) :: fcn   !! `(nx,ny)` matrix of function values to interpolate.
                                                !! `fcn(i,j)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`)
    integer(ip),intent(in)             :: kx    !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)             :: ky    !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    logical,intent(in),optional      :: extrap  !! if true, then extrapolation is allowed
                                                !! (default is false)

    call initialize_2d_auto_knots(me,x,y,fcn,kx,ky,me%iflag,extrap)

    end function bspline_2d_constructor_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_2d]] type (user-specified knots).
!  This is a wrapper for [[initialize_2d_specify_knots]].

    pure function bspline_2d_constructor_specify_knots(x,y,fcn,kx,ky,tx,ty,extrap) result(me)

    implicit none

    type(bspline_2d)                   :: me
    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) :: fcn   !! `(nx,ny)` matrix of function values to interpolate.
                                                !! `fcn(i,j)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`)
    integer(ip),intent(in)             :: kx    !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)             :: ky    !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)   :: tx    !! The `(nx+kx)` knots in the \(x\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)   :: ty    !! The `(ny+ky)` knots in the \(y\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    logical,intent(in),optional      :: extrap  !! if true, then extrapolation is allowed
                                                !! (default is false)

    call initialize_2d_specify_knots(me,x,y,fcn,kx,ky,tx,ty,me%iflag,extrap)

    end function bspline_2d_constructor_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_2d]] type (with automatically-computed knots).
!  This is a wrapper for [[db2ink]].

    pure subroutine initialize_2d_auto_knots(me,x,y,fcn,kx,ky,iflag,extrap)

    implicit none

    class(bspline_2d),intent(inout)    :: me
    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) :: fcn   !! `(nx,ny)` matrix of function values to interpolate.
                                                !! `fcn(i,j)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`)
    integer(ip),intent(in)             :: kx    !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)             :: ky    !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(out)            :: iflag !! status flag (see [[db2ink]])
    logical,intent(in),optional        :: extrap !! if true, then extrapolation is allowed
                                                 !! (default is false)

    integer(ip) :: iknot
    integer(ip) :: nx,ny

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)

    me%nx = nx
    me%ny = ny

    me%kx = kx
    me%ky = ky

    allocate(me%tx(nx+kx))
    allocate(me%ty(ny+ky))
    allocate(me%bcoef(nx,ny))
    allocate(me%work_val_1(ky))
    allocate(me%work_val_2(3_ip*max(kx,ky)))

    iknot = 0_ip         !knot sequence chosen by db2ink

    call db2ink(x,nx,y,ny,fcn,kx,ky,iknot,me%tx,me%ty,me%bcoef,iflag)

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

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

    end subroutine initialize_2d_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_2d]] type (with user-specified knots).
!  This is a wrapper for [[db2ink]].

    pure subroutine initialize_2d_specify_knots(me,x,y,fcn,kx,ky,tx,ty,iflag,extrap)

    implicit none

    class(bspline_2d),intent(inout)    :: me
    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) :: fcn   !! `(nx,ny)` matrix of function values to interpolate.
                                                !! `fcn(i,j)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`)
    integer(ip),intent(in)             :: kx    !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)             :: ky    !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)   :: tx    !! The `(nx+kx)` knots in the \(x\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)   :: ty    !! The `(ny+ky)` knots in the \(y\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    integer(ip),intent(out)            :: iflag !! status flag (see [[db2ink]])
    logical,intent(in),optional      :: extrap  !! if true, then extrapolation is allowed
                                                !! (default is false)

    integer(ip) :: nx,ny

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)

    call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,&
                                  ny=ny,ky=ky,ty=ty,&
                                  iflag=iflag)

    if (iflag == 0_ip) then

        me%nx = nx
        me%ny = ny

        me%kx = kx
        me%ky = ky

        allocate(me%tx(nx+kx))
        allocate(me%ty(ny+ky))
        allocate(me%bcoef(nx,ny))
        allocate(me%work_val_1(ky))
        allocate(me%work_val_2(3_ip*max(kx,ky)))

        me%tx = tx
        me%ty = ty

        call db2ink(x,nx,y,ny,fcn,kx,ky,1_ip,me%tx,me%ty,me%bcoef,iflag)

        call me%set_extrap_flag(extrap)

    end if

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

    end subroutine initialize_2d_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Evaluate a [[bspline_2d]] interpolate.  This is a wrapper for [[db2val]].

    pure subroutine evaluate_2d(me,xval,yval,idx,idy,f,iflag)

    implicit none

    class(bspline_2d),intent(inout) :: me
    real(wp),intent(in)             :: xval  !! \(x\) coordinate of evaluation point.
    real(wp),intent(in)             :: yval  !! \(y\) coordinate of evaluation point.
    integer(ip),intent(in)           :: idx   !! \(x\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)           :: idy   !! \(y\) derivative of piecewise polynomial to evaluate.
    real(wp),intent(out)            :: f     !! interpolated value
    integer(ip),intent(out)         :: iflag !! status flag (see [[db2val]])

    if (me%initialized) then
        call db2val(xval,yval,&
                    idx,idy,&
                    me%tx,me%ty,&
                    me%nx,me%ny,&
                    me%kx,me%ky,&
                    me%bcoef,f,iflag,&
                    me%inbvx,me%inbvy,me%iloy,&
                    me%work_val_1,me%work_val_2,&
                    extrap=me%extrap)
    else
        iflag = 1_ip
    end if

    me%iflag = iflag

    end subroutine evaluate_2d
!*****************************************************************************************

!*****************************************************************************************
!>
!  It returns an empty [[bspline_3d]] type. Note that INITIALIZE still
!  needs to be called before it can be used.
!  Not really that useful except perhaps in some OpenMP applications.

    elemental function bspline_3d_constructor_empty() result(me)

    implicit none

    type(bspline_3d) :: me

    end function bspline_3d_constructor_empty
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_3d]] type (auto knots).
!  This is a wrapper for [[initialize_3d_auto_knots]].

    pure function bspline_3d_constructor_auto_knots(x,y,z,fcn,kx,ky,kz,extrap) result(me)

    implicit none

    type(bspline_3d)                     :: me
    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) :: fcn !! `(nx,ny,nz)` matrix of function values to interpolate.
                                                !! `fcn(i,j,k)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`,`z(k)`)
    integer(ip),intent(in)               :: kx  !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)               :: ky  !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)               :: kz  !! The order of spline pieces in \(z\)
                                                !! ( \( 2 \le k_z < n_z \) )
                                                !! (order = polynomial degree + 1)
    logical,intent(in),optional       :: extrap !! if true, then extrapolation is allowed
                                                !! (default is false)

    call initialize_3d_auto_knots(me,x,y,z,fcn,kx,ky,kz,me%iflag,extrap)

    end function bspline_3d_constructor_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_3d]] type (user-specified knots).
!  This is a wrapper for [[initialize_3d_specify_knots]].

    pure function bspline_3d_constructor_specify_knots(x,y,z,fcn,kx,ky,kz,tx,ty,tz,extrap) result(me)

    implicit none

    type(bspline_3d)                     :: me
    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) :: fcn !! `(nx,ny,nz)` matrix of function values to interpolate.
                                                !! `fcn(i,j,k)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`,`z(k)`)
    integer(ip),intent(in)               :: kx  !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)               :: ky  !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)               :: kz  !! The order of spline pieces in \(z\)
                                                !! ( \( 2 \le k_z < n_z \) )
                                                !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)     :: tx  !! The `(nx+kx)` knots in the \(x\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)     :: ty  !! The `(ny+ky)` knots in the \(y\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)     :: tz  !! The `(nz+kz)` knots in the \(z\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    logical,intent(in),optional :: extrap       !! if true, then extrapolation is allowed
                                                !! (default is false)

    call initialize_3d_specify_knots(me,x,y,z,fcn,kx,ky,kz,tx,ty,tz,me%iflag,extrap)

    end function bspline_3d_constructor_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_3d]] type (with automatically-computed knots).
!  This is a wrapper for [[db3ink]].

    pure subroutine initialize_3d_auto_knots(me,x,y,z,fcn,kx,ky,kz,iflag,extrap)

    implicit none

    class(bspline_3d),intent(inout)      :: me
    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) :: fcn !! `(nx,ny,nz)` matrix of function values to interpolate.
                                                !! `fcn(i,j,k)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`,`z(k)`)
    integer(ip),intent(in)               :: kx  !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)               :: ky  !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)               :: kz  !! The order of spline pieces in \(z\)
                                                !! ( \( 2 \le k_z < n_z \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(out) :: iflag            !! status flag (see [[db3ink]])
    logical,intent(in),optional :: extrap       !! if true, then extrapolation is allowed
                                                !! (default is false)

    integer(ip) :: iknot
    integer(ip) :: nx,ny,nz

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)
    nz = size(z,kind=ip)

    me%nx = nx
    me%ny = ny
    me%nz = nz

    me%kx = kx
    me%ky = ky
    me%kz = kz

    allocate(me%tx(nx+kx))
    allocate(me%ty(ny+ky))
    allocate(me%tz(nz+kz))
    allocate(me%bcoef(nx,ny,nz))
    allocate(me%work_val_1(ky,kz))
    allocate(me%work_val_2(kz))
    allocate(me%work_val_3(3_ip*max(kx,ky,kz)))

    iknot = 0_ip         !knot sequence chosen by db3ink

    call db3ink(x,nx,y,ny,z,nz,&
                fcn,&
                kx,ky,kz,&
                iknot,&
                me%tx,me%ty,me%tz,&
                me%bcoef,iflag)

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

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

    end subroutine initialize_3d_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_3d]] type (with user-specified knots).
!  This is a wrapper for [[db3ink]].

    pure subroutine initialize_3d_specify_knots(me,x,y,z,fcn,kx,ky,kz,tx,ty,tz,iflag,extrap)

    implicit none

    class(bspline_3d),intent(inout)      :: me
    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) :: fcn !! `(nx,ny,nz)` matrix of function values to interpolate.
                                                !! `fcn(i,j,k)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`,`z(k)`)
    integer(ip),intent(in)               :: kx  !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)               :: ky  !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)               :: kz  !! The order of spline pieces in \(z\)
                                                !! ( \( 2 \le k_z < n_z \) )
                                                !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)     :: tx  !! The `(nx+kx)` knots in the \(x\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)     :: ty  !! The `(ny+ky)` knots in the \(y\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)     :: tz  !! The `(nz+kz)` knots in the \(z\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    integer(ip),intent(out)  :: iflag               !! status flag (see [[db3ink]])
    logical,intent(in),optional  :: extrap      !! if true, then extrapolation is allowed
                                                !! (default is false)

    integer(ip) :: nx,ny,nz

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)
    nz = size(z,kind=ip)

    call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,&
                                  ny=ny,ky=ky,ty=ty,&
                                  nz=nz,kz=kz,tz=tz,&
                                  iflag=iflag)

    if (iflag == 0_ip) then

        me%nx = nx
        me%ny = ny
        me%nz = nz

        me%kx = kx
        me%ky = ky
        me%kz = kz

        allocate(me%tx(nx+kx))
        allocate(me%ty(ny+ky))
        allocate(me%tz(nz+kz))
        allocate(me%bcoef(nx,ny,nz))
        allocate(me%work_val_1(ky,kz))
        allocate(me%work_val_2(kz))
        allocate(me%work_val_3(3_ip*max(kx,ky,kz)))

        me%tx = tx
        me%ty = ty
        me%tz = tz

        call db3ink(x,nx,y,ny,z,nz,&
                    fcn,&
                    kx,ky,kz,&
                    1_ip,&
                    me%tx,me%ty,me%tz,&
                    me%bcoef,iflag)

        call me%set_extrap_flag(extrap)

    end if

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

    end subroutine initialize_3d_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Evaluate a [[bspline_3d]] interpolate.  This is a wrapper for [[db3val]].

    pure subroutine evaluate_3d(me,xval,yval,zval,idx,idy,idz,f,iflag)

    implicit none

    class(bspline_3d),intent(inout) :: me
    real(wp),intent(in)             :: xval  !! \(x\) coordinate of evaluation point.
    real(wp),intent(in)             :: yval  !! \(y\) coordinate of evaluation point.
    real(wp),intent(in)             :: zval  !! \(z\) coordinate of evaluation point.
    integer(ip),intent(in)          :: idx   !! \(x\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idy   !! \(y\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idz   !! \(z\) derivative of piecewise polynomial to evaluate.
    real(wp),intent(out)            :: f     !! interpolated value
    integer(ip),intent(out)         :: iflag !! status flag (see [[db3val]])

    if (me%initialized) then
        call db3val(xval,yval,zval,&
                    idx,idy,idz,&
                    me%tx,me%ty,me%tz,&
                    me%nx,me%ny,me%nz,&
                    me%kx,me%ky,me%kz,&
                    me%bcoef,f,iflag,&
                    me%inbvx,me%inbvy,me%inbvz,&
                    me%iloy,me%iloz,&
                    me%work_val_1,me%work_val_2,me%work_val_3,&
                    extrap=me%extrap)
    else
        iflag = 1_ip
    end if

    me%iflag = iflag

    end subroutine evaluate_3d
!*****************************************************************************************

!*****************************************************************************************
!>
!  It returns an empty [[bspline_4d]] type. Note that INITIALIZE still
!  needs to be called before it can be used.
!  Not really that useful except perhaps in some OpenMP applications.

    elemental function bspline_4d_constructor_empty() result(me)

    implicit none

    type(bspline_4d) :: me

    end function bspline_4d_constructor_empty
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_4d]] type (auto knots).
!  This is a wrapper for [[initialize_4d_auto_knots]].

    pure function bspline_4d_constructor_auto_knots(x,y,z,q,fcn,kx,ky,kz,kq,extrap) result(me)

    implicit none

    type(bspline_4d)                       :: me
    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) :: fcn !! `(nx,ny,nz,nq)` matrix of function values to interpolate.
                                                  !! `fcn(i,j,k,l)` should contain the function value at the
                                                  !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`)
    integer(ip),intent(in)                 :: kx  !! The order of spline pieces in \(x\)
                                                  !! ( \( 2 \le k_x < n_x \) )
                                                  !! (order = polynomial degree + 1)
    integer(ip),intent(in)                 :: ky  !! The order of spline pieces in \(y\)
                                                  !! ( \( 2 \le k_y < n_y \) )
                                                  !! (order = polynomial degree + 1)
    integer(ip),intent(in)                 :: kz  !! The order of spline pieces in \(z\)
                                                  !! ( \( 2 \le k_z < n_z \) )
                                                  !! (order = polynomial degree + 1)
    integer(ip),intent(in)                 :: kq  !! The order of spline pieces in \(q\)
                                                  !! ( \( 2 \le k_q < n_q \) )
                                                  !! (order = polynomial degree + 1)
    logical,intent(in),optional      :: extrap    !! if true, then extrapolation is allowed
                                                  !! (default is false)

    call initialize_4d_auto_knots(me,x,y,z,q,fcn,kx,ky,kz,kq,me%iflag,extrap)

    end function bspline_4d_constructor_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_4d]] type (user-specified knots).
!  This is a wrapper for [[initialize_4d_specify_knots]].

    pure function bspline_4d_constructor_specify_knots(x,y,z,q,fcn,kx,ky,kz,kq,&
                                                        tx,ty,tz,tq,extrap) result(me)

    implicit none

    type(bspline_4d)                           :: me
    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)     :: fcn !! `(nx,ny,nz,nq)` matrix of function values to interpolate.
                                                      !! `fcn(i,j,k,l)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`)
    integer(ip),intent(in)                     :: kx  !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ky  !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kz  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kq  !! The order of spline pieces in \(q\)
                                                      !! ( \( 2 \le k_q < n_q \) )
                                                      !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)           :: tx  !! The `(nx+kx)` knots in the \(x\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: ty  !! The `(ny+ky)` knots in the \(y\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tz  !! The `(nz+kz)` knots in the \(z\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tq  !! The `(nq+kq)` knots in the \(q\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    logical,intent(in),optional      :: extrap        !! if true, then extrapolation is allowed
                                                      !! (default is false)

    call initialize_4d_specify_knots(me,x,y,z,q,fcn,kx,ky,kz,kq,tx,ty,tz,tq,me%iflag,extrap)

    end function bspline_4d_constructor_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_4d]] type (with automatically-computed knots).
!  This is a wrapper for [[db4ink]].

    pure subroutine initialize_4d_auto_knots(me,x,y,z,q,fcn,kx,ky,kz,kq,iflag,extrap)

    implicit none

    class(bspline_4d),intent(inout)            :: me
    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)     :: fcn !! `(nx,ny,nz,nq)` matrix of function values to interpolate.
                                                      !! `fcn(i,j,k,l)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`)
    integer(ip),intent(in)                     :: kx  !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ky  !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kz  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kq  !! The order of spline pieces in \(q\)
                                                      !! ( \( 2 \le k_q < n_q \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(out)  :: iflag                     !! status flag (see [[db4ink]])
    logical,intent(in),optional  :: extrap            !! if true, then extrapolation is allowed
                                                      !! (default is false)

    integer(ip) :: iknot
    integer(ip) :: nx,ny,nz,nq

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)
    nz = size(z,kind=ip)
    nq = size(q,kind=ip)

    me%nx = nx
    me%ny = ny
    me%nz = nz
    me%nq = nq

    me%kx = kx
    me%ky = ky
    me%kz = kz
    me%kq = kq

    allocate(me%tx(nx+kx))
    allocate(me%ty(ny+ky))
    allocate(me%tz(nz+kz))
    allocate(me%tq(nq+kq))
    allocate(me%bcoef(nx,ny,nz,nq))
    allocate(me%work_val_1(ky,kz,kq))
    allocate(me%work_val_2(kz,kq))
    allocate(me%work_val_3(kq))
    allocate(me%work_val_4(3_ip*max(kx,ky,kz,kq)))

    iknot = 0_ip         !knot sequence chosen by db4ink

    call db4ink(x,nx,y,ny,z,nz,q,nq,&
                fcn,&
                kx,ky,kz,kq,&
                iknot,&
                me%tx,me%ty,me%tz,me%tq,&
                me%bcoef,iflag)

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

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

    end subroutine initialize_4d_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_4d]] type (with user-specified knots).
!  This is a wrapper for [[db4ink]].

    pure subroutine initialize_4d_specify_knots(me,x,y,z,q,fcn,&
                                                kx,ky,kz,kq,tx,ty,tz,tq,iflag,extrap)

    implicit none

    class(bspline_4d),intent(inout)            :: me
    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)     :: fcn !! `(nx,ny,nz,nq)` matrix of function values to interpolate.
                                                      !! `fcn(i,j,k,l)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`)
    integer(ip),intent(in)                     :: kx  !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ky  !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kz  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kq  !! The order of spline pieces in \(q\)
                                                      !! ( \( 2 \le k_q < n_q \) )
                                                      !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)           :: tx  !! The `(nx+kx)` knots in the \(x\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: ty  !! The `(ny+ky)` knots in the \(y\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tz  !! The `(nz+kz)` knots in the \(z\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tq  !! The `(nq+kq)` knots in the \(q\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    integer(ip),intent(out)  :: iflag                     !! status flag (see [[db4ink]])
    logical,intent(in),optional  :: extrap            !! if true, then extrapolation is allowed
                                                      !! (default is false)

    integer(ip) :: nx,ny,nz,nq

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)
    nz = size(z,kind=ip)
    nq = size(q,kind=ip)

    call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,&
                                  ny=ny,ky=ky,ty=ty,&
                                  nz=nz,kz=kz,tz=tz,&
                                  nq=nq,kq=kq,tq=tq,&
                                  iflag=iflag)

    if (iflag == 0_ip) then

        me%nx = nx
        me%ny = ny
        me%nz = nz
        me%nq = nq

        me%kx = kx
        me%ky = ky
        me%kz = kz
        me%kq = kq

        allocate(me%tx(nx+kx))
        allocate(me%ty(ny+ky))
        allocate(me%tz(nz+kz))
        allocate(me%tq(nq+kq))
        allocate(me%bcoef(nx,ny,nz,nq))
        allocate(me%work_val_1(ky,kz,kq))
        allocate(me%work_val_2(kz,kq))
        allocate(me%work_val_3(kq))
        allocate(me%work_val_4(3_ip*max(kx,ky,kz,kq)))

        me%tx = tx
        me%ty = ty
        me%tz = tz
        me%tq = tq

        call db4ink(x,nx,y,ny,z,nz,q,nq,&
                    fcn,&
                    kx,ky,kz,kq,&
                    1_ip,&
                    me%tx,me%ty,me%tz,me%tq,&
                    me%bcoef,iflag)

        call me%set_extrap_flag(extrap)

    end if

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

    end subroutine initialize_4d_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Evaluate a [[bspline_4d]] interpolate.  This is a wrapper for [[db4val]].

    pure subroutine evaluate_4d(me,xval,yval,zval,qval,idx,idy,idz,idq,f,iflag)

    implicit none

    class(bspline_4d),intent(inout) :: me
    real(wp),intent(in)             :: xval  !! \(x\) coordinate of evaluation point.
    real(wp),intent(in)             :: yval  !! \(y\) coordinate of evaluation point.
    real(wp),intent(in)             :: zval  !! \(z\) coordinate of evaluation point.
    real(wp),intent(in)             :: qval  !! \(q\) coordinate of evaluation point.
    integer(ip),intent(in)          :: idx   !! \(x\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idy   !! \(y\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idz   !! \(z\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idq   !! \(q\) derivative of piecewise polynomial to evaluate.
    real(wp),intent(out)            :: f     !! interpolated value
    integer(ip),intent(out)         :: iflag !! status flag (see [[db4val]])

    if (me%initialized) then
        call db4val(xval,yval,zval,qval,&
                    idx,idy,idz,idq,&
                    me%tx,me%ty,me%tz,me%tq,&
                    me%nx,me%ny,me%nz,me%nq,&
                    me%kx,me%ky,me%kz,me%kq,&
                    me%bcoef,f,iflag,&
                    me%inbvx,me%inbvy,me%inbvz,me%inbvq,&
                    me%iloy,me%iloz,me%iloq,&
                    me%work_val_1,me%work_val_2,me%work_val_3,me%work_val_4,&
                    extrap=me%extrap)
    else
        iflag = 1_ip
    end if

    me%iflag = iflag

    end subroutine evaluate_4d
!*****************************************************************************************

!*****************************************************************************************
!>
!  It returns an empty [[bspline_5d]] type. Note that INITIALIZE still
!  needs to be called before it can be used.
!  Not really that useful except perhaps in some OpenMP applications.

    elemental function bspline_5d_constructor_empty() result(me)

    implicit none

    type(bspline_5d) :: me

    end function bspline_5d_constructor_empty
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_5d]] type (auto knots).
!  This is a wrapper for [[initialize_5d_auto_knots]].

    pure function bspline_5d_constructor_auto_knots(x,y,z,q,r,fcn,kx,ky,kz,kq,kr,extrap) result(me)

    implicit none

    type(bspline_5d)                           :: me
    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)` matrix of function values to interpolate.
                                                      !! `fcn(i,j,k,l,m)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`)
    integer(ip),intent(in)                     :: kx  !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ky  !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kz  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kq  !! The order of spline pieces in \(q\)
                                                      !! ( \( 2 \le k_q < n_q \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kr  !! The order of spline pieces in \(r\)
                                                      !! ( \( 2 \le k_r < n_r \) )
                                                      !! (order = polynomial degree + 1)
    logical,intent(in),optional      :: extrap        !! if true, then extrapolation is allowed
                                                      !! (default is false)

    call initialize_5d_auto_knots(me,x,y,z,q,r,fcn,kx,ky,kz,kq,kr,me%iflag,extrap)

    end function bspline_5d_constructor_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_5d]] type (user-specified knots).
!  This is a wrapper for [[initialize_5d_specify_knots]].

    pure function bspline_5d_constructor_specify_knots(x,y,z,q,r,fcn,&
                                                        kx,ky,kz,kq,kr,&
                                                        tx,ty,tz,tq,tr,extrap) result(me)

    implicit none

    type(bspline_5d)                           :: me
    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)` matrix of function values to interpolate.
                                                      !! `fcn(i,j,k,l,m)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`)
    integer(ip),intent(in)                     :: kx  !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ky  !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kz  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kq  !! The order of spline pieces in \(q\)
                                                      !! ( \( 2 \le k_q < n_q \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kr  !! The order of spline pieces in \(r\)
                                                      !! ( \( 2 \le k_r < n_r \) )
                                                      !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)           :: tx  !! The `(nx+kx)` knots in the \(x\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: ty  !! The `(ny+ky)` knots in the \(y\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tz  !! The `(nz+kz)` knots in the \(z\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tq  !! The `(nq+kq)` knots in the \(q\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tr  !! The `(nr+kr)` knots in the \(r\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    logical,intent(in),optional      :: extrap        !! if true, then extrapolation is allowed
                                                      !! (default is false)

    call initialize_5d_specify_knots(me,x,y,z,q,r,fcn,kx,ky,kz,kq,kr,tx,ty,tz,tq,tr,me%iflag,extrap)

    end function bspline_5d_constructor_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_5d]] type (with automatically-computed knots).
!  This is a wrapper for [[db5ink]].

    pure subroutine initialize_5d_auto_knots(me,x,y,z,q,r,fcn,kx,ky,kz,kq,kr,iflag,extrap)

    implicit none

    class(bspline_5d),intent(inout)            :: me
    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)` matrix of function values to interpolate.
                                                      !! `fcn(i,j,k,l,m)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`)
    integer(ip),intent(in)                     :: kx  !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ky  !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kz  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kq  !! The order of spline pieces in \(q\)
                                                      !! ( \( 2 \le k_q < n_q \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kr  !! The order of spline pieces in \(r\)
                                                      !! ( \( 2 \le k_r < n_r \) )
                                                      !! (order = polynomial degree + 1)
    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) :: iknot
    integer(ip) :: nx,ny,nz,nq,nr

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)
    nz = size(z,kind=ip)
    nq = size(q,kind=ip)
    nr = size(r,kind=ip)

    me%nx = nx
    me%ny = ny
    me%nz = nz
    me%nq = nq
    me%nr = nr

    me%kx = kx
    me%ky = ky
    me%kz = kz
    me%kq = kq
    me%kr = kr

    allocate(me%tx(nx+kx))
    allocate(me%ty(ny+ky))
    allocate(me%tz(nz+kz))
    allocate(me%tq(nq+kq))
    allocate(me%tr(nr+kr))
    allocate(me%bcoef(nx,ny,nz,nq,nr))
    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)))

    iknot = 0_ip         !knot sequence chosen by db5ink

    call db5ink(x,nx,y,ny,z,nz,q,nq,r,nr,&
                fcn,&
                kx,ky,kz,kq,kr,&
                iknot,&
                me%tx,me%ty,me%tz,me%tq,me%tr,&
                me%bcoef,iflag)

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

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

    end subroutine initialize_5d_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_5d]] type (with user-specified knots).
!  This is a wrapper for [[db5ink]].

    pure subroutine initialize_5d_specify_knots(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
    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)` matrix of function values to interpolate.
                                                      !! `fcn(i,j,k,l,m)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`)
    integer(ip),intent(in)                     :: kx  !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ky  !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kz  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kq  !! The order of spline pieces in \(q\)
                                                      !! ( \( 2 \le k_q < n_q \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kr  !! The order of spline pieces in \(r\)
                                                      !! ( \( 2 \le k_r < n_r \) )
                                                      !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)           :: tx  !! The `(nx+kx)` knots in the \(x\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: ty  !! The `(ny+ky)` knots in the \(y\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tz  !! The `(nz+kz)` knots in the \(z\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tq  !! The `(nq+kq)` knots in the \(q\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tr  !! The `(nr+kr)` knots in the \(r\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    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) :: nx,ny,nz,nq,nr

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)
    nz = size(z,kind=ip)
    nq = size(q,kind=ip)
    nr = size(r,kind=ip)

    call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,&
                                  ny=ny,ky=ky,ty=ty,&
                                  nz=nz,kz=kz,tz=tz,&
                                  nq=nq,kq=kq,tq=tq,&
                                  nr=nr,kr=kr,tr=tr,&
                                  iflag=iflag)

    if (iflag == 0_ip) then

        me%nx = nx
        me%ny = ny
        me%nz = nz
        me%nq = nq
        me%nr = nr

        me%kx = kx
        me%ky = ky
        me%kz = kz
        me%kq = kq
        me%kr = kr

        allocate(me%tx(nx+kx))
        allocate(me%ty(ny+ky))
        allocate(me%tz(nz+kz))
        allocate(me%tq(nq+kq))
        allocate(me%tr(nr+kr))
        allocate(me%bcoef(nx,ny,nz,nq,nr))
        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,nx,y,ny,z,nz,q,nq,r,nr,&
                    fcn,&
                    kx,ky,kz,kq,kr,&
                    1_ip,&
                    me%tx,me%ty,me%tz,me%tq,me%tr,&
                    me%bcoef,iflag)

        call me%set_extrap_flag(extrap)

    end if

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

    end subroutine initialize_5d_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Evaluate a [[bspline_5d]] interpolate.  This is a wrapper for [[db5val]].

    pure subroutine evaluate_5d(me,xval,yval,zval,qval,rval,idx,idy,idz,idq,idr,f,iflag)

    implicit none

    class(bspline_5d),intent(inout) :: me
    real(wp),intent(in)             :: xval  !! \(x\) coordinate of evaluation point.
    real(wp),intent(in)             :: yval  !! \(y\) coordinate of evaluation point.
    real(wp),intent(in)             :: zval  !! \(z\) coordinate of evaluation point.
    real(wp),intent(in)             :: qval  !! \(q\) coordinate of evaluation point.
    real(wp),intent(in)             :: rval  !! \(r\) coordinate of evaluation point.
    integer(ip),intent(in)          :: idx   !! \(x\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idy   !! \(y\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idz   !! \(z\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idq   !! \(q\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idr   !! \(r\) derivative of piecewise polynomial to evaluate.
    real(wp),intent(out)            :: f     !! interpolated value
    integer(ip),intent(out)         :: iflag !! status flag (see [[db5val]])

    if (me%initialized) then
        call db5val(xval,yval,zval,qval,rval,&
                    idx,idy,idz,idq,idr,&
                    me%tx,me%ty,me%tz,me%tq,me%tr,&
                    me%nx,me%ny,me%nz,me%nq,me%nr,&
                    me%kx,me%ky,me%kz,me%kq,me%kr,&
                    me%bcoef,f,iflag,&
                    me%inbvx,me%inbvy,me%inbvz,me%inbvq,me%inbvr,&
                    me%iloy,me%iloz,me%iloq,me%ilor,&
                    me%work_val_1,me%work_val_2,me%work_val_3,me%work_val_4,me%work_val_5,&
                    extrap=me%extrap)
    else
        iflag = 1_ip
    end if

    me%iflag = iflag

    end subroutine evaluate_5d
!*****************************************************************************************

!*****************************************************************************************
!>
!  It returns an empty [[bspline_6d]] type. Note that INITIALIZE still
!  needs to be called before it can be used.
!  Not really that useful except perhaps in some OpenMP applications.

    elemental function bspline_6d_constructor_empty() result(me)

    implicit none

    type(bspline_6d) :: me

    end function bspline_6d_constructor_empty
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_6d]] type (auto knots).
!  This is a wrapper for [[initialize_6d_auto_knots]].

    pure function bspline_6d_constructor_auto_knots(x,y,z,q,r,s,fcn,&
                                                    kx,ky,kz,kq,kr,ks,extrap) result(me)

    implicit none

    type(bspline_6d)                           :: me
    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)           :: s   !! `(ns)` array of \(s\) abcissae. Must be strictly increasing.
    real(wp),dimension(:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,ns)` matrix of function values to interpolate.
                                                      !! `fcn(i,j,k,l,m,n)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`,`s(n)`)
    integer(ip),intent(in)                      :: kx  !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                      :: ky  !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                      :: kz  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                      :: kq  !! The order of spline pieces in \(q\)
                                                      !! ( \( 2 \le k_q < n_q \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                      :: kr  !! The order of spline pieces in \(r\)
                                                      !! ( \( 2 \le k_r < n_r \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                      :: ks  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    logical,intent(in),optional      :: extrap        !! if true, then extrapolation is allowed
                                                      !! (default is false)

    call initialize_6d_auto_knots(me,x,y,z,q,r,s,fcn,kx,ky,kz,kq,kr,ks,me%iflag,extrap)

    end function bspline_6d_constructor_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_6d]] type (user-specified knots).
!  This is a wrapper for [[initialize_6d_specify_knots]].

    pure function bspline_6d_constructor_specify_knots(x,y,z,q,r,s,fcn,&
                                                        kx,ky,kz,kq,kr,ks,&
                                                        tx,ty,tz,tq,tr,ts,extrap) result(me)

    implicit none

    type(bspline_6d)                           :: me
    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)           :: s   !! `(ns)` array of \(s\) abcissae. Must be strictly increasing.
    real(wp),dimension(:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,ns)` matrix of function values to interpolate.
                                                      !! `fcn(i,j,k,l,m,n)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`,`s(n)`)
    integer(ip),intent(in)                     :: kx  !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ky  !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kz  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kq  !! The order of spline pieces in \(q\)
                                                      !! ( \( 2 \le k_q < n_q \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kr  !! The order of spline pieces in \(r\)
                                                      !! ( \( 2 \le k_r < n_r \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ks  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)           :: tx  !! The `(nx+kx)` knots in the \(x\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: ty  !! The `(ny+ky)` knots in the \(y\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tz  !! The `(nz+kz)` knots in the \(z\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tq  !! The `(nq+kq)` knots in the \(q\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tr  !! The `(nr+kr)` knots in the \(r\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: ts  !! The `(ns+ks)` knots in the \(s\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    logical,intent(in),optional      :: extrap        !! if true, then extrapolation is allowed
                                                      !! (default is false)

    call initialize_6d_specify_knots(me,x,y,z,q,r,s,fcn,&
                                        kx,ky,kz,kq,kr,ks,&
                                        tx,ty,tz,tq,tr,ts,me%iflag,extrap)

    end function bspline_6d_constructor_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_6d]] type (with automatically-computed knots).
!  This is a wrapper for [[db6ink]].

    pure subroutine initialize_6d_auto_knots(me,x,y,z,q,r,s,fcn,&
                                                kx,ky,kz,kq,kr,ks,iflag,extrap)

    implicit none

    class(bspline_6d),intent(inout)            :: me
    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)           :: s   !! `(ns)` array of \(s\) abcissae. Must be strictly increasing.
    real(wp),dimension(:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,ns)` matrix of function values to interpolate.
                                                      !! `fcn(i,j,k,l,m,n)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`,`s(n)`)
    integer(ip),intent(in)                     :: kx  !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ky  !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kz  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kq  !! The order of spline pieces in \(q\)
                                                      !! ( \( 2 \le k_q < n_q \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kr  !! The order of spline pieces in \(r\)
                                                      !! ( \( 2 \le k_r < n_r \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ks  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(out)  :: iflag                     !! status flag (see [[db6ink]])
    logical,intent(in),optional  :: extrap            !! if true, then extrapolation is allowed
                                                      !! (default is false)

    integer(ip) :: iknot
    integer(ip) :: nx,ny,nz,nq,nr,ns

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)
    nz = size(z,kind=ip)
    nq = size(q,kind=ip)
    nr = size(r,kind=ip)
    ns = size(s,kind=ip)

    me%nx = nx
    me%ny = ny
    me%nz = nz
    me%nq = nq
    me%nr = nr
    me%ns = ns

    me%kx = kx
    me%ky = ky
    me%kz = kz
    me%kq = kq
    me%kr = kr
    me%ks = ks

    allocate(me%tx(nx+kx))
    allocate(me%ty(ny+ky))
    allocate(me%tz(nz+kz))
    allocate(me%tq(nq+kq))
    allocate(me%tr(nr+kr))
    allocate(me%ts(ns+ks))
    allocate(me%bcoef(nx,ny,nz,nq,nr,ns))
    allocate(me%work_val_1(ky,kz,kq,kr,ks))
    allocate(me%work_val_2(kz,kq,kr,ks))
    allocate(me%work_val_3(kq,kr,ks))
    allocate(me%work_val_4(kr,ks))
    allocate(me%work_val_5(ks))
    allocate(me%work_val_6(3_ip*max(kx,ky,kz,kq,kr,ks)))

    iknot = 0_ip         !knot sequence chosen by db6ink

    call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,&
                fcn,&
                kx,ky,kz,kq,kr,ks,&
                iknot,&
                me%tx,me%ty,me%tz,me%tq,me%tr,me%ts,&
                me%bcoef,iflag)

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

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

    end subroutine initialize_6d_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_6d]] type (with user-specified knots).
!  This is a wrapper for [[db6ink]].

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

    implicit none

    class(bspline_6d),intent(inout)            :: me
    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)           :: s   !! `(ns)` array of \(s\) abcissae. Must be strictly increasing.
    real(wp),dimension(:,:,:,:,:,:),intent(in) :: fcn !! `(nx,ny,nz,nq,nr,ns)` matrix of function values to interpolate.
                                                      !! `fcn(i,j,k,l,m,n)` should contain the function value at the
                                                      !! point (`x(i)`,`y(j)`,`z(k)`,`q(l)`,`r(m)`,`s(n)`)
    integer(ip),intent(in)                     :: kx  !! The order of spline pieces in \(x\)
                                                      !! ( \( 2 \le k_x < n_x \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ky  !! The order of spline pieces in \(y\)
                                                      !! ( \( 2 \le k_y < n_y \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kz  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kq  !! The order of spline pieces in \(q\)
                                                      !! ( \( 2 \le k_q < n_q \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: kr  !! The order of spline pieces in \(r\)
                                                      !! ( \( 2 \le k_r < n_r \) )
                                                      !! (order = polynomial degree + 1)
    integer(ip),intent(in)                     :: ks  !! The order of spline pieces in \(z\)
                                                      !! ( \( 2 \le k_z < n_z \) )
                                                      !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)           :: tx  !! The `(nx+kx)` knots in the \(x\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: ty  !! The `(ny+ky)` knots in the \(y\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tz  !! The `(nz+kz)` knots in the \(z\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tq  !! The `(nq+kq)` knots in the \(q\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: tr  !! The `(nr+kr)` knots in the \(r\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)           :: ts  !! The `(ns+ks)` knots in the \(s\) direction
                                                      !! for the spline interpolant.
                                                      !! Must be non-decreasing.
    integer(ip),intent(out)  :: iflag                     !! status flag (see [[db6ink]])
    logical,intent(in),optional  :: extrap            !! if true, then extrapolation is allowed
                                                      !! (default is false)

    integer(ip) :: nx,ny,nz,nq,nr,ns

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)
    nz = size(z,kind=ip)
    nq = size(q,kind=ip)
    nr = size(r,kind=ip)
    ns = size(s,kind=ip)

    call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,&
                                  ny=ny,ky=ky,ty=ty,&
                                  nz=nz,kz=kz,tz=tz,&
                                  nq=nq,kq=kq,tq=tq,&
                                  nr=nr,kr=kr,tr=tr,&
                                  ns=ns,ks=ks,ts=ts,&
                                  iflag=iflag)

    if (iflag == 0_ip) then

        me%nx = nx
        me%ny = ny
        me%nz = nz
        me%nq = nq
        me%nr = nr
        me%ns = ns

        me%kx = kx
        me%ky = ky
        me%kz = kz
        me%kq = kq
        me%kr = kr
        me%ks = ks

        allocate(me%tx(nx+kx))
        allocate(me%ty(ny+ky))
        allocate(me%tz(nz+kz))
        allocate(me%tq(nq+kq))
        allocate(me%tr(nr+kr))
        allocate(me%ts(ns+ks))
        allocate(me%bcoef(nx,ny,nz,nq,nr,ns))
        allocate(me%work_val_1(ky,kz,kq,kr,ks))
        allocate(me%work_val_2(kz,kq,kr,ks))
        allocate(me%work_val_3(kq,kr,ks))
        allocate(me%work_val_4(kr,ks))
        allocate(me%work_val_5(ks))
        allocate(me%work_val_6(3_ip*max(kx,ky,kz,kq,kr,ks)))

        me%tx = tx
        me%ty = ty
        me%tz = tz
        me%tq = tq
        me%tr = tr
        me%ts = ts

        call db6ink(x,nx,y,ny,z,nz,q,nq,r,nr,s,ns,&
                    fcn,&
                    kx,ky,kz,kq,kr,ks,&
                    1_ip,&
                    me%tx,me%ty,me%tz,me%tq,me%tr,me%ts,&
                    me%bcoef,iflag)

        call me%set_extrap_flag(extrap)

    end if

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

    end subroutine initialize_6d_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Evaluate a [[bspline_6d]] interpolate.  This is a wrapper for [[db6val]].

    pure subroutine evaluate_6d(me,xval,yval,zval,qval,rval,sval,idx,idy,idz,idq,idr,ids,f,iflag)

    implicit none

    class(bspline_6d),intent(inout) :: me
    real(wp),intent(in)             :: xval  !! \(x\) coordinate of evaluation point.
    real(wp),intent(in)             :: yval  !! \(y\) coordinate of evaluation point.
    real(wp),intent(in)             :: zval  !! \(z\) coordinate of evaluation point.
    real(wp),intent(in)             :: qval  !! \(q\) coordinate of evaluation point.
    real(wp),intent(in)             :: rval  !! \(r\) coordinate of evaluation point.
    real(wp),intent(in)             :: sval  !! \(s\) coordinate of evaluation point.
    integer(ip),intent(in)          :: idx   !! \(x\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idy   !! \(y\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idz   !! \(z\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idq   !! \(q\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: idr   !! \(r\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)          :: ids   !! \(s\) derivative of piecewise polynomial to evaluate.
    real(wp),intent(out)            :: f     !! interpolated value
    integer(ip),intent(out)         :: iflag !! status flag (see [[db6val]])

    if (me%initialized) then
        call db6val(xval,yval,zval,qval,rval,sval,&
                    idx,idy,idz,idq,idr,ids,&
                    me%tx,me%ty,me%tz,me%tq,me%tr,me%ts,&
                    me%nx,me%ny,me%nz,me%nq,me%nr,me%ns,&
                    me%kx,me%ky,me%kz,me%kq,me%kr,me%ks,&
                    me%bcoef,f,iflag,&
                    me%inbvx,me%inbvy,me%inbvz,me%inbvq,me%inbvr,me%inbvs,&
                    me%iloy,me%iloz,me%iloq,me%ilor,me%ilos,&
                    me%work_val_1,me%work_val_2,me%work_val_3,me%work_val_4,me%work_val_5,me%work_val_6,&
                    extrap=me%extrap)
    else
        iflag = 1_ip
    end if

    me%iflag = iflag

    end subroutine evaluate_6d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Error checks for the user-specified knot vector sizes.
!
!@note If more than one is the wrong size, then the `iflag` error code will
!      correspond to the one with the highest rank.

    pure subroutine check_knot_vectors_sizes(nx,ny,nz,nq,nr,ns,&
                                             kx,ky,kz,kq,kr,ks,&
                                             tx,ty,tz,tq,tr,ts,iflag)

    implicit none

    integer(ip),intent(in),optional           :: nx
    integer(ip),intent(in),optional           :: ny
    integer(ip),intent(in),optional           :: nz
    integer(ip),intent(in),optional           :: nq
    integer(ip),intent(in),optional           :: nr
    integer(ip),intent(in),optional           :: ns
    integer(ip),intent(in),optional           :: kx
    integer(ip),intent(in),optional           :: ky
    integer(ip),intent(in),optional           :: kz
    integer(ip),intent(in),optional           :: kq
    integer(ip),intent(in),optional           :: kr
    integer(ip),intent(in),optional           :: ks
    real(wp),dimension(:),intent(in),optional :: tx
    real(wp),dimension(:),intent(in),optional :: ty
    real(wp),dimension(:),intent(in),optional :: tz
    real(wp),dimension(:),intent(in),optional :: tq
    real(wp),dimension(:),intent(in),optional :: tr
    real(wp),dimension(:),intent(in),optional :: ts
    integer(ip),intent(out)                   :: iflag  !! 0 if everything is OK

    iflag = 0_ip

    if (present(nx) .and. present(kx) .and. present(tx)) then
        if (size(tx,kind=ip)/=(nx+kx)) then
            iflag = 501_ip  ! tx is not the correct size (nx+kx)
        end if
    end if

    if (present(ny) .and. present(ky) .and. present(ty)) then
        if (size(ty,kind=ip)/=(ny+ky)) then
            iflag = 502_ip  ! ty is not the correct size (ny+ky)
        end if
    end if

    if (present(nz) .and. present(kz) .and. present(tz)) then
        if (size(tz,kind=ip)/=(nz+kz)) then
            iflag = 503_ip  ! tz is not the correct size (nz+kz)
        end if
    end if

    if (present(nq) .and. present(kq) .and. present(tq)) then
        if (size(tq,kind=ip)/=(nq+kq)) then
            iflag = 504_ip  ! tq is not the correct size (nq+kq)
        end if
    end if

    if (present(nr) .and. present(kr) .and. present(tr)) then
        if (size(tr,kind=ip)/=(nr+kr)) then
            iflag = 505_ip  ! tr is not the correct size (nr+kr)
        end if
    end if

    if (present(ns) .and. present(ks) .and. present(ts)) then
        if (size(ts,kind=ip)/=(ns+ks)) then
            iflag = 506_ip  ! ts is not the correct size (ns+ks)
        end if
    end if

    end subroutine check_knot_vectors_sizes
!*****************************************************************************************

!*****************************************************************************************
    end module bspline_oo_module
!*****************************************************************************************