initialize_4d Subroutine

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

    implicit none

    class(linear_interp_4d),intent(inout)  :: me
    real(wp),dimension(:),intent(in)       :: x
    real(wp),dimension(:),intent(in)       :: y
    real(wp),dimension(:),intent(in)       :: z
    real(wp),dimension(:),intent(in)       :: q
    real(wp),dimension(:,:,:,:),intent(in) :: f
    integer,intent(out)                    :: istat !! `0`   : no problems,
                                                    !! `1`   : `x` is not strictly increasing,
                                                    !! `2`   : `y` is not strictly increasing,
                                                    !! `3`   : `z` is not strictly increasing,
                                                    !! `4`   : `q` is not strictly increasing,
                                                    !! `10`  : `x` is not equal to size(f,1),
                                                    !! `20`  : `y` is not equal to size(f,2),
                                                    !! `30`  : `z` is not equal to size(f,3),
                                                    !! `40`  : `q` is not equal to size(f,4),
                                                    !! `100` : cannot use linear interpolation for only one point.

    call me%destroy()

    istat = 0

    if (istat==0 .and. size(x)/=size(f,1)) istat = 10
    if (istat==0 .and. size(y)/=size(f,2)) istat = 20
    if (istat==0 .and. size(z)/=size(f,3)) istat = 30
    if (istat==0 .and. size(q)/=size(f,4)) istat = 40

    if (istat==0) then
        call me%check_inputs(x=x,y=y,z=z,q=q,ierr=istat)
        if (istat==0) then
            allocate(me%f(size(x),size(y),size(z),size(q))); me%f = f
            allocate(me%x(size(x))); me%x = x
            allocate(me%y(size(y))); me%y = y
            allocate(me%z(size(z))); me%z = z
            allocate(me%q(size(q))); me%q = q
            me%initialized = .true.
        end if
    end if

    end subroutine initialize_4d