run_eikcoefs_and_resample Subroutine

  subroutine run_eikcoefs_and_resample(ntheta_input, ntheta_target, nperiod, results, job_id)
    implicit none
    !> This is the number of theta grid points per two pi period used in the
    !> geometry calculation (actually one less than this).
    integer, intent(in out) :: ntheta_input
    !> This is the desired number of theta grid points per two pi period
    !> in the resulting geometrical coefficients used elsewhere.
    integer, intent(in) :: ntheta_target
    integer, intent(in) :: nperiod
    type(eikcoefs_output_type), intent(out) :: results
    integer, intent(in), optional :: job_id
    type(eikcoefs_output_type) :: full_results
    real, dimension(:), allocatable :: input_index_space, target_index_space
    integer :: ig, ntgrid_target, ntgrid_input
    real :: spacing

    ! First run eikcoefs using requested resolution
    call eikcoefs(ntheta_input, nperiod, full_results, job_id)

    ntgrid_input = (size(full_results%theta) - 1) / 2
    ntgrid_target =(ntheta_target*(2*nperiod-1)) / 2

    ! Set scalar outputs
    results%ntheta = ntheta_target
    results%drhodpsi = full_results%drhodpsi
    results%kxfac = full_results%kxfac
    results%qsf = full_results%qsf
    results%surfarea = full_results%surfarea
    results%dvdrhon = full_results%dvdrhon
    results%shat = full_results%shat
    results%dbetadrho = full_results%dbetadrho
    results%rhoc = full_results%rhoc

    ! Make an index space array for the original theta
    ! grid. Spans [0,1] in size(theta)=2*ntgrid_input points
    allocate(input_index_space(-ntgrid_input:ntgrid_input))
    spacing = 1.0 / (2 * ntgrid_input)
    do ig = -ntgrid_input, ntgrid_input
       input_index_space(ig) = (ig + ntgrid_input) * spacing
    end do

    ! Make new theta grid
    allocate(target_index_space(-ntgrid_target:ntgrid_target))

    ! First make target index space grid
    spacing = 1.0 / (2 * ntgrid_target)
    do ig = -ntgrid_target, ntgrid_target
       target_index_space(ig) = (ig + ntgrid_target) * spacing
    end do

    ! Now interpolate theta grid onto target index space grid
    call resample( input_index_space, full_results%theta, target_index_space, results%theta)

    ! Now interpolate other output quantities onto target theta grid
    call resample(full_results%theta, full_results%bmag, results%theta, results%bmag)
    call resample(full_results%theta, full_results%gradpar, results%theta, results%gradpar)
    call resample(full_results%theta, full_results%grho, results%theta, results%grho)
    call resample(full_results%theta, full_results%cdrift, results%theta, results%cdrift)
    call resample(full_results%theta, full_results%cvdrift, results%theta, results%cvdrift)
    call resample(full_results%theta, full_results%gbdrift, results%theta, results%gbdrift)
    call resample(full_results%theta, full_results%cdrift0, results%theta, results%cdrift0)
    call resample(full_results%theta, full_results%cvdrift0, results%theta, results%cvdrift0)
    call resample(full_results%theta, full_results%gbdrift0, results%theta, results%gbdrift0)
    call resample(full_results%theta, full_results%gds2, results%theta, results%gds2)
    call resample(full_results%theta, full_results%gds21, results%theta, results%gds21)
    call resample(full_results%theta, full_results%gds22, results%theta, results%gds22)
    call resample(full_results%theta, full_results%gds23, results%theta, results%gds23)
    call resample(full_results%theta, full_results%gds24, results%theta, results%gds24)
    call resample(full_results%theta, full_results%gds24_noq, results%theta, results%gds24_noq)
    call resample(full_results%theta, full_results%rplot, results%theta, results%rplot)
    call resample(full_results%theta, full_results%zplot, results%theta, results%zplot)
    call resample(full_results%theta, full_results%aplot, results%theta, results%aplot)
    call resample(full_results%theta, full_results%rprime, results%theta, results%rprime)
    call resample(full_results%theta, full_results%zprime, results%theta, results%zprime)
    call resample(full_results%theta, full_results%aprime, results%theta, results%aprime)
    call resample(full_results%theta, full_results%bpol, results%theta, results%bpol)

  contains
    subroutine resample(x_orig, y_orig, x_new, y_new)
      use splines, only: new_spline, delete_spline, splint, spline
      real, dimension(:), intent(in) :: x_orig, y_orig
      real, dimension(-ntgrid_target:), intent(in) :: x_new
      real, dimension(:), allocatable, intent(out) ::  y_new
      type(spline) :: the_spline
      integer :: i, lower, upper

      ! Spline the old data
      the_spline = new_spline(x_orig, y_orig)
      lower = lbound(x_new, 1) ; upper = ubound(x_new, 1)

      allocate(y_new(lower:upper))

      do i = lower, upper
         y_new(i) = splint(x_new(i), the_spline)
      end do

      ! Get rid of the spline
      call delete_spline(the_spline)
    end subroutine resample
  end subroutine run_eikcoefs_and_resample