s_compute_interpolation_points_non_uniform Subroutine

  subroutine s_compute_interpolation_points_non_uniform( self, tau )
    class(sll_t_spline_interpolator_1d), intent(in   ) :: self
    real(wp),               allocatable, intent(  out) :: tau(:)

   integer :: i, ntau
   real(wp), allocatable :: temp_knots(:)

   associate( nbasis   => self % bspl % nbasis, &
              ncells   => self % bspl % ncells, &
              degree   => self % bspl % degree, &
              xmin     => self % bspl % xmin  , &
              xmax     => self % bspl % xmax  , &
              dx       => self % dx           , &
              bc_xmin  => self %  bc_xmin     , &
              bc_xmax  => self %  bc_xmax     , &
              nbc_xmin => self % nbc_xmin     , &
              nbc_xmax => self % nbc_xmax )

   associate( breaks   => self % bspl % knots(1:ncells+1) )

      ! Determine size of tau and allocate tau
      ntau = nbasis - nbc_xmin - nbc_xmax
      allocate( tau(1:ntau) )

      ! Array of temporary knots needed to compute interpolation points
      ! using Greville-style averaging: tau(i) = average(temp_knots(i+1-degree:i))
      allocate( temp_knots( 2-degree : ntau ) )

      if ( bc_xmin == sll_p_periodic ) then

        associate( k => degree/2 )
          temp_knots(:) = self%bspl%knots(2-degree+k:ntau+k)
        end associate

      else

        associate( r => 2-degree, s => -nbc_xmin )
          select case (bc_xmin)
          case (sll_p_greville); temp_knots(r:s) = breaks(1)
          case (sll_p_hermite ); temp_knots(r:s) = 2.0_wp*breaks(1) - breaks(2+s-r:2:-1)
          end select
        end associate

        associate( r => -nbc_xmin+1, s => -nbc_xmin+1+ncells )
          temp_knots(r:s) = breaks(:)
        end associate

        associate( r => -nbc_xmin+1+ncells+1, s => ntau )
          select case (bc_xmax)
          case (sll_p_greville)
            temp_knots(r:s) = breaks(ncells+1)
          case (sll_p_hermite )
            temp_knots(r:s) = 2.0_wp*breaks(ncells+1) - breaks(ncells:ncells+r-s:-1)
          end select
        end associate

      end if

      ! Compute interpolation points using Greville-style averaging
      associate( inv_deg => 1.0_wp / real( degree, wp ) )
        do i = 1, ntau
          tau(i) = sum( temp_knots(i+1-degree:i) ) * inv_deg
        end do
      end associate

      ! Periodic case: apply periodic BCs to interpolation points
      if ( bc_xmin == sll_p_periodic ) then
        tau(:) = modulo( tau(:)-xmin, xmax-xmin ) + xmin

      ! Non-periodic case, odd degree: fix round-off issues
      else if ( self%odd == 1 ) then
        tau(1)    = xmin
        tau(ntau) = xmax
      end if

    end associate ! breaks
    end associate ! all other variables

  end subroutine s_compute_interpolation_points_non_uniform