get_weights Subroutine

  subroutine get_weights (maxpts_in, llim, ulim, nodes, wgts, ndiv, divmax, err_flag)
    implicit none
    integer, intent (in) :: maxpts_in
    real, intent (in) :: llim, ulim
    real, dimension (:), intent (in) :: nodes
    real, dimension (:), intent (out) :: wgts
    logical, intent (out) :: err_flag
    integer, intent (out) :: ndiv, divmax
    integer :: npts, rmndr, basepts, divrmndr, base_idx, idiv, epts, maxpts
    integer, dimension (:), allocatable :: divpts
    real :: wgt_max

    err_flag = .false.

    ! npts is the number of grid points in the integration interval
    npts = size(nodes)

    ! maxpts is the max number of pts in an integration subinterval. It
    ! is the lower of npts and maxpts_in (which in practice is always
    ! the input nmax).
    maxpts = min(maxpts_in, npts)

    do
       wgts = 0.0; divmax = npts
       wgt_max = abs(ulim - llim) * wgt_fac / npts

       ! only need to subdivide integration interval if maxpts < npts
       ! This will be the typical case unless nmax < ntheta.
       if (maxpts >= npts) then
          call get_intrvl_weights (llim, ulim, nodes, wgts)
       else
          rmndr = mod(npts - maxpts, maxpts - 1)

          ! if rmndr is 0, then each subinterval contains maxpts pts
          if (rmndr == 0) then
             ! ndiv is the number of subintervals
             ndiv = (npts - maxpts) / (maxpts - 1) + 1
             allocate (divpts(ndiv))
             ! divpts is an array containing the number of pts for each subinterval
             divpts = maxpts
          else
             ndiv = (npts - maxpts) / (maxpts - 1) + 2
             allocate (divpts(ndiv))
             ! epts is the effective number of pts after taking into account double
             ! counting of some grid points (those that are boundaries of subintervals
             ! are used twice)
             epts = npts + ndiv - 1
             basepts = epts / ndiv
             divrmndr = mod(epts, ndiv)

             ! determines if all intervals have same number of pts
             if (divrmndr == 0) then
                divpts = basepts
             else
                ! First divrmndr intervals get one extra point each
                divpts(:divrmndr) = basepts + 1
                divpts(divrmndr+1:) = basepts
             end if
          end if

          base_idx = 0

          ! loop calls subroutine to get weights for each subinterval
          do idiv = 1, ndiv
             if (idiv == 1) then
                call get_intrvl_weights (llim, nodes(base_idx+divpts(idiv)), &
                     nodes(base_idx+1:base_idx+divpts(idiv)), &
                     wgts(base_idx+1:base_idx+divpts(idiv)))
             else if (idiv == ndiv) then
                call get_intrvl_weights (nodes(base_idx+1), ulim, &
                     nodes(base_idx+1:base_idx+divpts(idiv)), &
                     wgts(base_idx+1:base_idx+divpts(idiv)))
             else
                call get_intrvl_weights (nodes(base_idx+1), nodes(base_idx+divpts(idiv)), &
                     nodes(base_idx+1:base_idx+divpts(idiv)), &
                     wgts(base_idx+1:base_idx+divpts(idiv)))
             end if
             base_idx = base_idx + divpts(idiv) - 1
          end do

          divmax = maxval(divpts)
          deallocate (divpts)
       end if

       ! check to make sure the weights do not get too large @Warning:
       ! This looks a bit unusual -- we're taking the magnitude of the
       ! maximum value, whilst usually we would find the maximum value
       ! of the magnitude.
       if (abs(maxval(wgts)) > wgt_max) then
          if (maxpts < 3) then
             ! If we can't split any more then set the error flag and leave.
             err_flag = .true.
             exit
          end if
          ! Otherwise we can reduce the allowed number of points per interval
          ! by one and try again.
          maxpts = divmax - 1
       else
          exit
       end if
    end do
  end subroutine get_weights