module sfincs_interface

   implicit none

   public :: get_neo_from_sfincs

   private

# ifdef USE_SFINCS
   integer :: nproc_sfincs
   integer :: irad_min, irad_max
   real :: Er_window
   logical :: includeXDotTerm
   logical :: includeElectricFieldTermInXiDot
!  logical :: includeRadialExBDrive
   integer :: magneticDriftScheme
   logical :: includePhi1
   logical :: includePhi1InKineticEquation
!  logical :: includePhi1InCollisionOperator
   integer :: geometryScheme
   integer :: VMECRadialOption
   integer :: coordinateSystem
   integer :: inputRadialCoordinate
   integer :: inputRadialCoordinateForGradients
   character(200) :: equilibriumFile
   real :: aHat, psiAHat, Delta
   real :: nu_n, dPhiHatdrN
   integer :: nxi, nx, ntheta, nzeta
   logical :: calculate_radial_electric_field
   logical :: read_sfincs_output_from_file
   logical :: sfincs_finished = .true.
   real, dimension(:), allocatable :: fprim_local, tprim_local
# endif

contains

   subroutine get_neo_from_sfincs(nradii, drho, f_neoclassical, phi_neoclassical, &
                                  dfneo_dalpha, dphineo_dalpha)

# ifdef USE_SFINCS
      use mp, only: proc0, iproc
      use mp, only: comm_split, comm_free
      use geometry, only: geo_surf
      use species, only: spec, nspec
# endif
      use mp, only: mp_abort
      use zgrid, only: nzgrid

      implicit none

      integer, intent(in) :: nradii
      real, intent(in) :: drho
      real, dimension(:, -nzgrid:, :, :, :, -nradii/2:), intent(out) :: f_neoclassical
      real, dimension(:, -nzgrid:, -nradii/2:), intent(out) :: phi_neoclassical
      real, dimension(:, -nzgrid:, :, :, :), intent(out) :: dfneo_dalpha
      real, dimension(:, -nzgrid:), intent(out) :: dphineo_dalpha

# ifdef USE_SFINCS
      integer :: sfincs_comm
      integer :: color, ierr
      integer :: irad
      real :: dPhiHatdrN_best_guess
      logical :: Er_converged
      integer :: nsfincs_calls

      if (.not. allocated(fprim_local)) allocate (fprim_local(nspec))
      if (.not. allocated(tprim_local)) allocate (tprim_local(nspec))

      if (proc0) call read_sfincs_parameters(nradii)
      call broadcast_sfincs_parameters
      if (iproc < nproc_sfincs) then
         color = 0
      else
         color = 1
      end if
      call comm_split(color, sfincs_comm, ierr)
      if (iproc < nproc_sfincs) then
         do irad = irad_min, irad_max
            ! get local values of -dlog(ns)/drho and -dlog(Ts)/drho
            ! using dlog(n)/drho = dlog(n0)/drho + delrho*d/drho(dlog(n)/drho)
            fprim_local = 1.0 / geo_surf%drhotordrho * (spec%fprim &
                                                        + irad * drho * (spec%fprim**2 - spec%d2ndr2) / geo_surf%drhotordrho)
            tprim_local = 1.0 / geo_surf%drhotordrho * (spec%tprim &
                                                        + irad * drho * (spec%tprim**2 - spec%d2Tdr2) / geo_surf%drhotordrho)
!          fprim_local = 1.0/geo_surf%drhotordrho*(spec%fprim - irad*drho*spec%d2ndr2)
!          tprim_local = 1.0/geo_surf%drhotordrho*(spec%tprim - irad*drho*spec%d2Tdr2)
            if (calculate_radial_electric_field) then

               ! get best guess at radial electric field
               ! using force balance with radial pressure gradient
               if (dPhiHatdrN > -9999.0) then
                  dPhiHatdrN_best_guess = dPhiHatdrN
               else
                  dPhiHatdrN_best_guess = fprim_local(1) + tprim_local(1)
               end if
               call iterate_sfincs_until_electric_field_converged(sfincs_comm, &
                                                                  irad, drho, irad_max, dPhiHatdrN_best_guess, &
                                                                  Er_converged, nsfincs_calls)

               if (proc0) then
                  write (*, *)
                  write (*, *) 'At irad= ', irad, 'Er_converged= ', Er_converged, &
                     'nsfincs_calls_required= ', nsfincs_calls, 'dPhiHatdrN= ', dPhiHatdrN
                  write (*, *)
               end if

               ! write_and_finish_sfincs manipulates sfincs output
               ! to get the neoclassical distribution function and potential
               ! on the stella (zed,alpha,vpa,mu) grid; it then
               ! deallocates sfincs arrays, etc. to make it ready
               ! for running again if need be
               call write_and_finish_sfincs(f_neoclassical(:, :, :, :, :, irad), &
                                            phi_neoclassical(:, :, irad), dfneo_dalpha, dphineo_dalpha, irad)
            else
               ! init_and_run_sfincs initializes sfincs,
               ! including passing geometry info if necessary;
               ! and runs sfincs (if requested)
               call init_and_run_sfincs(sfincs_comm, irad, drho, irad_max)
               ! write_and_finish_sfincs manipulates sfincs output
               ! to get the neoclassical distribution function and potential
               ! on the stella (zed,alpha,vpa,mu) grid; it then
               ! deallocates sfincs arrays, etc. to make it ready
               ! for running again if need be
               call write_and_finish_sfincs(f_neoclassical(:, :, :, :, :, irad), &
                                            phi_neoclassical(:, :, irad), dfneo_dalpha, dphineo_dalpha, irad)
            end if
         end do
      end if

      call comm_free(sfincs_comm, ierr)

      ! NB: NEED TO CHECK THIS BROADCAST OF SFINCS RESULTS
      do irad = irad_min, irad_max
         call broadcast_sfincs_output &
            (f_neoclassical(:, :, :, :, :, irad), phi_neoclassical(:, :, irad))
      end do
      call broadcast_sfincs_output(dfneo_dalpha, dphineo_dalpha)

      deallocate (fprim_local, tprim_local)

      if (proc0) then
         write (*, *)
         write (*, *) 'maxval(fneo): ', maxval(f_neoclassical(:, :, :, :, :, irad_min:irad_max)), &
            'maxval(phineo): ', maxval(phi_neoclassical(:, :, irad_min:irad_max))
         write (*, *)
      end if

      if ((irad_min /= -nradii / 2) .or. (irad_max /= nradii / 2)) &
           call mp_abort('WARNING: irad_min must equal -nradii/2 and irad_max must equal &
           & nradii/2 to proceed to stella calculation.  aborting.')

# else
      real :: dum
      ! this pointless dum assignment only here to avoid
      ! annoying warning messages during compilation
      ! about unused variable drho
      dum = drho
      f_neoclassical = 0.; phi_neoclassical = 0.
      dfneo_dalpha = 0.; dphineo_dalpha = 0.
      call mp_abort("to run with include_neoclassical_terms=.true., &
           & USE_SFINCS must be defined at compilation time.  Aborting.")
# endif

   end subroutine get_neo_from_sfincs

# ifdef USE_SFINCS
   subroutine iterate_sfincs_until_electric_field_converged( &
      sfincs_comm, irad, drho, nrad_max, dPhiHatdrN_best_guess, &
      dphiHatdrN_is_converged, number_of_sfincs_calls_for_convergence)

      use mp, only: proc0, iproc

      implicit none

      integer, intent(in) :: sfincs_comm, irad, nrad_max
      real, intent(in) :: drho, dPhiHatdrN_best_guess
      logical, intent(out) :: dPhiHatdrN_is_converged
      integer, intent(out) :: number_of_sfincs_calls_for_convergence

      integer :: itmax_bracket = 10
      integer :: itmax_root = 10
      real :: window = 0.3
      real :: tol = 0.1

      integer :: it
      real :: a, b, c, d, e, fa, fb, fc, p, q, r, s, tol1, xm, eps
      real :: converged_dPhiHatdrN

      dPhiHatdrN_is_converged = .false.

      a = dPhiHatdrN_best_guess * (1.0 - Er_window)
      b = dPhiHatdrN_best_guess * (1.0 + Er_window)
      ! initialize sfincs, run it, and return the total charge flux as fa
      call get_total_charge_flux(sfincs_comm, irad, drho, nrad_max, a, fa)
      call get_total_charge_flux(sfincs_comm, irad, drho, nrad_max, b, fb)
      number_of_sfincs_calls_for_convergence = 2
      do it = 1, itmax_bracket
         eps = epsilon(a)
         if ((fa > 0.0 .and. fb > 0.0) .or. (fa < 0.0 .and. fb < 0.0)) then
            if (proc0) then
               write (*, *)
               write (*, *) 'dPhiHatdrN values ', a, ' and ', b, ' do not bracket root.'
               write (*, *) 'flux at ', a, ' is ', fa, '.'
               write (*, *) 'flux at ', b, ' is ', fb, '.'
            end if
            a = a * (1.0 - Er_window)
            b = b * (1.0 + Er_window)
            if (proc0) then
               write (*, *) 'Trying again with values ', a, ' and ', b, ' .'
               write (*, *)
            end if
            call get_total_charge_flux(sfincs_comm, irad, drho, nrad_max, a, fa)
            call get_total_charge_flux(sfincs_comm, irad, drho, nrad_max, b, fb)
            number_of_sfincs_calls_for_convergence = number_of_sfincs_calls_for_convergence + 2

!           ! eliminate the endpoint corresonding to the flux that is furthest from zero in magnitude
!           if (abs(fa) > abs(fb)) then
!              ! keep b as an endpoint and eliminate a
!              a = b ; fa = fb
!              b = a*(1.0+window)
!              if (proc0) then
!                 write (*,*) 'Trying again with values ', a, ' and ', b, ' .'
!                 write (*,*)
!              end if
!              call get_total_charge_flux (sfincs_comm, irad, drho, nrad_max, b, fb)
!           else
!              ! keep a as an endpoint and eliminate b
!              b = a ; fb = fa
!              a = b*(1.0-window)
!              if (proc0) then
!                 write (*,*) 'Trying again with values ', a, ' and ', b, ' .'
!                 write (*,*)
!              end if
!              call get_total_charge_flux (sfincs_comm, irad, drho, nrad_max, a, fa)
!           end if
         else
            exit
         end if
      end do

      c = b
      fc = fb
      do it = 1, itmax_root
         if ((fb > 0.0 .and. fc > 0.0) .or. (fb < 0.0 .and. fc < 0.0)) then
            c = a
            fc = fa
            d = b - a
            e = d
         end if
         if (abs(fc) < abs(fb)) then
            a = b
            b = c
            c = a
            fa = fb
            fb = fc
            fc = fa
         end if
         tol1 = 2.0 * eps * abs(b) + 0.5 * tol
         xm = 0.5 * (c - b)
         if (abs(xm) <= tol1 .or. fb == 0.0) then
            converged_dPhiHatdrN = b
            dPhiHatdrN_is_converged = .true.
!          number_of_sfincs_calls_for_convergence = it+1
            exit
         end if
         if (abs(e) >= tol1 .and. abs(fa) > abs(fb)) then
            s = fb / fa
            if (a == c) then
               p = 2.0 * xm * s
               q = 1.0 - s
            else
               q = fa / fc
               r = fb / fc
               p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0))
               q = (q - 1.0) * (r - 1.0) * (s - 1.0)
            end if
            if (p > 0.0) q = -q
            p = abs(p)
            if (2.0 * p < min(3.0 * xm * q - abs(tol1 * q), abs(e * q))) then
               e = d
               d = p / q
            else
               d = xm
               e = d
            end if
         else
            d = xm
            e = d
         end if
         a = b
         fa = fb
         b = b + merge(d, sign(tol1, xm), abs(d) > tol1)
         call get_total_charge_flux(sfincs_comm, irad, drho, nrad_max, b, fb)
         number_of_sfincs_calls_for_convergence = number_of_sfincs_calls_for_convergence + 1
      end do

   end subroutine iterate_sfincs_until_electric_field_converged

   subroutine get_total_charge_flux(sfincs_comm, irad, drho, nrad_max, &
                                    dPhiHatdrN_in, total_charge_flux)

      use mp, only: iproc, broadcast_with_comm
      use sfincs_main, only: finish_sfincs
      use globalVariables, only: Zs, particleFlux_vd_psiHat
      use species, only: nspec

      implicit none

      integer, intent(in) :: sfincs_comm, irad, nrad_max
      real, intent(in) :: drho
      real, intent(in) :: dPhiHatdrN_in
      real, intent(out) :: total_charge_flux

      dPhiHatdrN = dPhiHatdrN_in
      call init_and_run_sfincs(sfincs_comm, irad, drho, nrad_max)

      call broadcast_with_comm(Zs, sfincs_comm)
      call broadcast_with_comm(particleFlux_vd_psiHat, sfincs_comm)

      total_charge_flux = sum(Zs(:nspec) * particleFlux_vd_psiHat(:nspec))

   end subroutine get_total_charge_flux

   subroutine init_and_run_sfincs(sfincs_comm, irad, drho, nrad_max)

      use mp, only: proc0, iproc
      use sfincs_main, only: init_sfincs, prepare_sfincs, run_sfincs, finish_sfincs
      use globalVariables, only: Zs, particleFlux_vd_psiHat
      use species, only: nspec

      implicit none

      integer, intent(in) :: sfincs_comm, irad, nrad_max
      real, intent(in) :: drho

      if (.not. sfincs_finished) then
         call finish_sfincs
         sfincs_finished = .true.
      end if
      call init_sfincs(sfincs_comm)
      call pass_inputoptions_to_sfincs(irad * drho)
      call pass_outputoptions_to_sfincs
      call prepare_sfincs
      ! if geometryScheme = 5, then sfincs will read in equilibrium
      ! parameters from vmec file separately
      ! otherwise, assume system is axisymmetric and pass geometry
      ! from stella (miller local equilibrium or similar)
      if (geometryScheme /= 5) call pass_geometry_to_sfincs(irad * drho)
      if (read_sfincs_output_from_file) then
         if (proc0) call read_sfincs_output(irad, nrad_max)
      else
         if (proc0) then
            write (*, *)
            write (*, *) 'Running sfincs at irad= ', irad, ', with dPhiHatdrN= ', dPhiHatdrN
         end if
         call run_sfincs
         if (proc0) then
            write (*, *) 'sfincs finished running.  total charge flux= ', sum(Zs(:nspec) * particleFlux_vd_psiHat(:nspec))
            write (*, *)
         end if
         ! write Phi1Hat and delta_f to file
         ! so we have the option of using it
         ! again without re-running sfincs
         if (proc0) call write_sfincs(irad, nrad_max)
      end if
      sfincs_finished = .false.

   end subroutine init_and_run_sfincs

   subroutine write_and_finish_sfincs(fneo, phineo, dfneo, dphineo, irad)

      use mp, only: proc0
      use sfincs_main, only: finish_sfincs
      use zgrid, only: nzgrid

      implicit none

      real, dimension(:, -nzgrid:, :, :, :), intent(out) :: fneo
      real, dimension(:, -nzgrid:), intent(out) :: phineo
      real, dimension(:, -nzgrid:, :, :, :), intent(in out) :: dfneo
      real, dimension(:, -nzgrid:), intent(in out) :: dphineo
      integer, intent(in) :: irad

      if (proc0) then
         ! only need to compute dfneo_dalpha and dphineo_dalpha
         ! for central radius and for stellarator calculation
         if (irad == 0 .and. geometryScheme == 5) then
            call get_sfincs_output(fneo, phineo, dfneo, dphineo)
         else
            call get_sfincs_output(fneo, phineo)
         end if
      end if

      call finish_sfincs

   end subroutine write_and_finish_sfincs

   subroutine read_sfincs_parameters(nradii)

      use constants, only: pi
      use mp, only: nproc
      use file_utils, only: input_unit_exist
      use species, only: nspec
      use physics_parameters, only: rhostar, vnew_ref
      use geometry, only: geo_surf, aref, bref

      implicit none

      integer, intent(in) :: nradii

      namelist /sfincs_input/ nproc_sfincs, &
         calculate_radial_electric_field, &
         includeXDotTerm, &
         includeElectricFieldTermInXiDot, &
         irad_min, irad_max, &
         magneticDriftScheme, &
         includePhi1, &
         includePhi1InKineticEquation, &
         !         includePhi1InCollisionOperator, &
         geometryScheme, &
         VMECRadialOption, &
         equilibriumFile, &
         coordinateSystem, &
         inputRadialCoordinate, &
         inputRadialCoordinateForGradients, &
         aHat, psiAHat, nu_N, nxi, nx, Delta, &
         dPhiHatdrN, &
         ntheta, nzeta, &
         read_sfincs_output_from_file, Er_window

      logical :: exist
      integer :: in_file

      ! if read_sfincs_output_from_file=.true.,
      ! will try to read in Phi1Hat and delta_f
      ! from pre-saved file named sfincs.output
      ! otherwise, run sfincs to compute these
      ! quantities on the fly
      read_sfincs_output_from_file = .false.
      ! number of processors to use for sfincs calculation
      nproc_sfincs = 1
      ! minimum and maximum radial index (irad=0 corresponds to central radius)
      irad_min = -nradii / 2; irad_max = nradii / 2
      ! if calculate_radial_electric_field, then
      ! will scan in radial electric field to find value
      ! for which ambipolarity is satisfied, and then
      ! use this value to obtain neoclassical fluxes,
      ! distribution function, and potential
      calculate_radial_electric_field = .true.
      ! do not include radial electric field term if set to .false.
      includeXDotTerm = .true.
      includeElectricFieldTermInXiDot = .true.
      ! include v_E . grad r term
!    includeRadialExBDrive = .true.
      ! no poloidal or toroidal magnetic drifts
      magneticDriftScheme = 0
      ! combo of next two variables means
      ! phi1 will be calculated via quasineutrality
      includePhi1 = .true.
      includePhi1InKineticEquation = .false.
!    includePhi1InCollisionOperator = .false.
      ! will be overridden by direct input of geometric quantities
      ! unless geometryScheme = 5 (vmec equilibrium)
      geometryScheme = 1
      ! only relevant if geometryScheme = 5
      ! radial option to use for vmec equilibrium
      ! 0 corresponds to using radial interpolation to get desired surface
      ! 1 corresponds to using nearest surface on VMEC HALF grid
      ! 2 corresponds to using nearest surface on VMEC FULL grid
      ! should not change this unless self-consistently change in the
      ! vmec input namelist
      VMECRadialOption = 0
      ! path of vmec equilibrium file
      equilibriumFile = 'wout_161s1.nc'
      ! seems to be a nonsensical option
      coordinateSystem = 3
      ! option 3 corresponds to using sqrt of toroidal flux
      ! normalized by toroidal flux enclosed by the LCFS
      inputRadialCoordinate = 3
      ! option 3 corresponds to same choice
      ! when calculating gradients of density, temperature, and potential
      inputRadialCoordinateForGradients = 3
      ! corresponds to r_LCFS as reference length in sfincs
      ! only used in sfincs when geometryScheme=1
      aHat = 1.0
      ! psitor_LCFS / (B_ref * a_ref^2)
      psiAHat = geo_surf%psitor_lcfs
      ! Delta is rho* = mref*vt_ref/(e*Bref*aref), with reference
      ! quantities given in SI units
      ! unless geometryScheme = 5, in which case Bref=1T
      ! and aref = 1m (these are hardwired in sfincs)
      ! set negative to allow check later to see if any value given in input file
      Delta = -1.0
      ! nu_n = nu_ref * aref/vt_ref
      ! nu_ref = 4*sqrt(2*pi)*nref*e**4*loglam/(3*sqrt(mref)*Tref**3/2)
      ! (with nref, Tref, and mref in Gaussian units)
      ! set negative to allow check later to see if any value given in input file
      nu_n = -1.0
      ! radial derivative of normalized phi
      dPhiHatdrN = -9999.9
      Er_window = 0.3
      ! number of spectral coefficients in pitch angle
      nxi = 48
      ! number of speeds
      nx = 12
      ! number of poloidal angles
      Ntheta = 65
      ! number of toroidal angles, 1 is appropriate for tokamak
      Nzeta = 1

      in_file = input_unit_exist("sfincs_input", exist)
      if (exist) read (unit=in_file, nml=sfincs_input)

      if (nproc_sfincs > nproc) then
         write (*, *) 'requested number of processors for sfincs is greater &
              & than total processor count.'
         write (*, *) 'allocating ', nproc, ' processors for sfincs.'
      end if

      if (Delta < 0.0) then
         Delta = rhostar
         ! if geometryScheme=5, Bref=1T and aref=1m are hard-wired in sfincs
         ! but these are not the values used in stella to define rhostar
         if (geometryScheme == 5) Delta = rhostar * bref * aref
      end if

      if (nu_n < 0.0) then
         nu_n = vnew_ref * (4./(3.*sqrt(pi)))
         ! if geometryScheme=5, aref=1m is hard-wired in sfincs
         ! but this is not the value used in stella
         if (geometryScheme == 5) nu_n = nu_n / aref
      end if

! FLAG -- NOT YET SURE IF THIS SHOULD BE HERE
!    if (nspec == 1 .and. includePhi1) then
!       write (*,*) 'includePhi1 = .true. is incompatible with a single-species run.'
!       write (*,*) 'forcing includePhi1 = .false.'
!       includePhi1 = .false.
!    end if

      ! ensure that ntheta and nzeta are odd for SFINCS
      ntheta = 2 * (ntheta / 2) + 1
      nzeta = 2 * (nzeta / 2) + 1

   end subroutine read_sfincs_parameters

   subroutine broadcast_sfincs_parameters

      use mp, only: broadcast
      use physics_parameters, only: rhostar

      implicit none

      call broadcast(read_sfincs_output_from_file)
      call broadcast(nproc_sfincs)
      call broadcast(irad_min)
      call broadcast(irad_max)
      call broadcast(calculate_radial_electric_field)
      call broadcast(includeXDotTerm)
      call broadcast(includeElectricFieldTermInXiDot)
!    call broadcast (includeRadialExBDrive)
      call broadcast(magneticDriftScheme)
      call broadcast(includePhi1)
      call broadcast(includePhi1InKineticEquation)
!    call broadcast (includePhi1InCollisionOperator)
      call broadcast(geometryScheme)
      call broadcast(VMECRadialOption)
      call broadcast(equilibriumFile)
      call broadcast(coordinateSystem)
      call broadcast(inputRadialCoordinate)
      call broadcast(inputRadialCoordinateForGradients)
      call broadcast(aHat)
      call broadcast(psiAHat)
      call broadcast(Delta)
      call broadcast(nu_N)
      call broadcast(dPhiHatdrN)
      call broadcast(Er_window)
      call broadcast(nxi)
      call broadcast(nx)
      call broadcast(ntheta)
      call broadcast(nzeta)

      write (*, *) 'Delta stella', Delta, rhostar

   end subroutine broadcast_sfincs_parameters

   subroutine pass_inputoptions_to_sfincs(delrho)

      use mp, only: mp_abort
      use geometry, only: geo_surf
      use species, only: spec, nspec
      use zgrid, only: nzed
      use physics_parameters, only: nine, tite
      use globalVariables, only: includeXDotTerm_sfincs => includeXDotTerm
      use globalVariables, only: includeElectricFieldTermInXiDot_sfincs => includeElectricFieldTermInXiDot
!    use globalVariables, only: includeRadialExBDrive_sfincs => includeRadialExBDrive
      use globalVariables, only: magneticDriftScheme_sfincs => magneticDriftScheme
      use globalVariables, only: includePhi1_sfincs => includePhi1
      use globalVariables, only: includePhi1InKineticEquation_sfincs => includePhi1InKineticEquation
!    use globalVariables, only: includePhi1InCollisionOperator_sfincs => includePhi1InCollisionOperator
      use globalVariables, only: geometryScheme_sfincs => geometryScheme
      use globalVariables, only: equilibriumFile_sfincs => equilibriumFile
      use globalVariables, only: VMECRadialOption_sfincs => VMECRadialOption
      use globalVariables, only: coordinateSystem_sfincs => coordinateSystem
      use globalVariables, only: RadialCoordinate => inputRadialCoordinate
      use globalVariables, only: RadialCoordinateForGradients => inputRadialCoordinateForGradients
      use globalVariables, only: rN_wish
      use globalVariables, only: Nspecies, nHats, THats, MHats, Zs
      use globalVariables, only: adiabaticNHat, adiabaticTHat, adiabaticZ
      use globalVariables, only: nxi_sfincs => Nxi
      use globalVariables, only: nx_sfincs => Nx
      use globalVariables, only: ntheta_sfincs => Ntheta
      use globalVariables, only: nzeta_sfincs => Nzeta
      use globalVariables, only: dnHatdrNs, dTHatdrNs
      use globalVariables, only: aHat_sfincs => aHat
      use globalVariables, only: psiAHat_sfincs => psiAHat
      use globalVariables, only: Delta_sfincs => Delta
      use globalVariables, only: nu_n_sfincs => nu_n
      use globalVariables, only: dPhiHatdrN_sfincs => dPhiHatdrN
      use globalVariables, only: withAdiabatic

      implicit none

      real, intent(in) :: delrho

      includeXDotTerm_sfincs = includeXDotTerm
      includeElectricFieldTermInXiDot_sfincs = includeElectricFieldTermInXiDot
!    includeRadialExBDrive_sfincs = includeRadialExBDrive
      magneticDriftScheme_sfincs = magneticDriftScheme
      includePhi1_sfincs = includePhi1
      includePhi1InKineticEquation_sfincs = includePhi1InKineticEquation
!    includePhi1InCollisionOperator_sfincs = includePhi1InCollisionOperator
      geometryScheme_sfincs = geometryScheme
      VMECRadialOption_sfincs = VMECRadialOption
      equilibriumFile_sfincs = trim(equilibriumFile)
      coordinateSystem_sfincs = coordinateSystem
      RadialCoordinate = inputRadialCoordinate
      RadialCoordinateForGradients = inputRadialCoordinateForGradients
      Nspecies = nspec
      nHats(:nspec) = spec%dens * (1.0 - delrho * spec%fprim)
      THats(:nspec) = spec%temp * (1.0 - delrho * spec%tprim)
      mHats(:nspec) = spec%mass
      Zs(:nspec) = spec%z
      nzeta_sfincs = nzeta
      ntheta_sfincs = ntheta
      nx_sfincs = nx
      nxi_sfincs = nxi
      aHat_sfincs = aHat
      psiAHat_sfincs = psiAHat
      Delta_sfincs = Delta
      nu_n_sfincs = nu_n
      dPhiHatdrN_sfincs = dPhiHatdrN
      if (nspec == 1) then
         withAdiabatic = .true.
         adiabaticNHat = nHats(1) / nine
         adiabaticTHat = THats(1) / tite
         adiabaticZ = -1
      end if

      if (inputRadialCoordinate == 3) then
         rN_wish = geo_surf%rhotor + delrho * geo_surf%drhotordrho
      else
         call mp_abort('only inputRadialCoordinate=3 currently supported. aborting.')
      end if
      if (inputRadialCoordinateForGradients == 3) then
         ! radial density gradient with respect to rhotor = sqrt(psitor/psitor_LCFS)
         ! normalized by reference density (not species density)
         dnHatdrNs(:nspec) = -spec%dens * fprim_local
         ! radial temperature gradient with respect to rhotor = sqrt(psitor/psitor_LCFS)
         ! normalized by reference tmperatures (not species temperature)
         dTHatdrNs(:nspec) = -spec%temp * tprim_local
      else
         call mp_abort('only inputRadialCoordinateForGradients=3 currently supported. aborting.')
      end if

   end subroutine pass_inputoptions_to_sfincs

   subroutine pass_outputoptions_to_sfincs
      use export_f, only: export_f_theta_option
      use export_f, only: export_f_zeta_option
      use export_f, only: export_f_xi_option
      use export_f, only: export_f_x_option
      use export_f, only: export_delta_f, export_full_f
      use export_f, only: export_f_stella
      implicit none
      export_f_theta_option = 0
      export_f_zeta_option = 0
      export_f_xi_option = 0
      export_f_x_option = 0
      export_delta_f = .true.
      export_full_f = .false.
      export_f_stella = .true.
   end subroutine pass_outputoptions_to_sfincs

   ! if this subroutine is being called, then
   ! using sfincs in tokamak geometry
   ! so zed in stella is theta
   subroutine pass_geometry_to_sfincs(delrho)

      use constants, only: pi
      use splines, only: linear_interp_periodic
      use zgrid, only: nz2pi, zed
      use geometry, only: bmag, dbdzed, gradpar
      use geometry, only: dBdrho, d2Bdrdth, dgradpardrho, dIdrho
      use geometry, only: geo_surf
      use globalVariables, only: BHat
      use globalVariables, only: dBHatdtheta
      use globalVariables, only: iota
      use globalVariables, only: DHat
      use globalVariables, only: BHat_sup_theta
      use globalVariables, only: BHat_sub_zeta
      use export_f, only: export_f_theta

      implicit none

      real, intent(in) :: delrho

      integer :: nzeta = 1
      integer :: nzpi
      real :: q_local
      real, dimension(:), allocatable :: B_local, dBdz_local, gradpar_local
      real, dimension(:), allocatable :: zed_stella
      real, dimension(:), allocatable :: theta_sfincs

      nzpi = nz2pi / 2
      allocate (B_local(-nzpi:nzpi))
      allocate (dBdz_local(-nzpi:nzpi))
      allocate (gradpar_local(-nzpi:nzpi))
      allocate (theta_sfincs(ntheta))
      allocate (zed_stella(-nzpi:nzpi))

      call init_zero_arrays

      ! first get some geometric quantities at this radius
      ! for theta from -pi to pi
      q_local = geo_surf%qinp * (1.0 + delrho * geo_surf%shat / geo_surf%rhoc)
      B_local = bmag(1, -nzpi:nzpi) + delrho * dBdrho(-nzpi:nzpi)
      dBdz_local = dbdzed(1, -nzpi:nzpi) + delrho * d2Bdrdth(-nzpi:nzpi)
      gradpar_local = gradpar(-nzpi:nzpi) + delrho * dgradpardrho(-nzpi:nzpi)

      zed_stella = zed(-nzpi:nzpi) + pi
      theta_sfincs = export_f_theta(:ntheta)

      iota = 1./q_local

      ! interpolate from stella zed-grid to sfincs theta grid
      ! point at -pi (stella) is same as point at 0 (sfincs)
      BHat(1, 1) = B_local(-nzpi)
      call linear_interp_periodic(zed_stella, B_local, theta_sfincs(2:), BHat(2:, 1))
      ! FLAG -- needs to be changed for stellarator runs
      BHat = spread(BHat(:, 1), 2, nzeta)

      dBHatdtheta(1, 1) = dBdz_local(-nzpi)
      call linear_interp_periodic(zed_stella, dBdz_local, theta_sfincs(2:), dBHatdtheta(2:, 1))
      dBHatdtheta = spread(dBHatdtheta(:, 1), 2, nzeta)

      ! this is bhat . grad theta
      BHat_sup_theta(1, 1) = B_local(-nzpi) * gradpar_local(-nzpi)
      call linear_interp_periodic(zed_stella, B_local * gradpar_local, theta_sfincs(2:), BHat_sup_theta(2:, 1))
      BHat_sup_theta = spread(BHat_sup_theta(:, 1), 2, nzeta)
      ! this is I(psi) / (aref*Bref)
      BHat_sub_zeta = geo_surf%rgeo + delrho * dIdrho
      ! this is grad psitor . (grad theta x grad zeta)
      ! note that + sign below relies on B = I grad zeta + grad zeta x grad psi
      DHat = q_local * BHat_sup_theta

      deallocate (B_local, dBdz_local, gradpar_local)
      deallocate (theta_sfincs, zed_stella)

   end subroutine pass_geometry_to_sfincs

   subroutine init_zero_arrays
      use globalVariables, only: dBHatdzeta
      use globalVariables, only: dBHatdpsiHat
      use globalVariables, only: BHat_sup_zeta
      use globalVariables, only: BHat_sub_psi
      use globalVariables, only: BHat_sub_theta
      use globalVariables, only: dBHat_sub_psi_dtheta
      use globalVariables, only: dBHat_sub_psi_dzeta
      use globalVariables, only: dBHat_sub_theta_dpsiHat
      use globalVariables, only: dBHat_sub_theta_dzeta
      use globalVariables, only: dBHat_sub_zeta_dpsiHat
      use globalVariables, only: dBHat_sub_zeta_dtheta
      use globalVariables, only: dBHat_sup_theta_dpsiHat
      use globalVariables, only: dBHat_sup_theta_dzeta
      use globalVariables, only: dBHat_sup_zeta_dpsiHat
      use globalVariables, only: dBHat_sup_zeta_dtheta
      implicit none
      dBHatdzeta = 0.
      dBHatdpsiHat = 0.
      BHat_sup_zeta = 0.
      BHat_sub_psi = 0.
      BHat_sub_theta = 0.
      dBHat_sub_psi_dtheta = 0.
      dBHat_sub_psi_dzeta = 0.
      dBHat_sub_theta_dpsiHat = 0.
      dBHat_sub_theta_dzeta = 0.
      dBHat_sub_zeta_dpsiHat = 0.
      dBHat_sub_zeta_dtheta = 0.
      dBHat_sup_theta_dpsiHat = 0.
      dBHat_sup_theta_dzeta = 0.
      dBHat_sup_zeta_dpsiHat = 0.
      dBHat_sup_zeta_dtheta = 0.
   end subroutine init_zero_arrays

   subroutine get_sfincs_output(f_neoclassical, phi_neoclassical, &
                                dfneo_dalpha, dphineo_dalpha)

      use constants, only: pi
      use sort, only: sort_array_ascending, unsort_array_ascending
      use species, only: nspec
      use zgrid, only: nzgrid, nz2pi
      use export_f, only: h_sfincs => delta_f
      use globalVariables, only: Phi1Hat
      use parameters_kxky_grids, only: nalpha

      implicit none

      real, dimension(:, -nzgrid:, :, :, :), intent(out) :: f_neoclassical
      real, dimension(:, -nzgrid:), intent(out) :: phi_neoclassical
      real, dimension(:, -nzgrid:, :, :, :), intent(out), optional :: dfneo_dalpha
      real, dimension(:, -nzgrid:), intent(out), optional :: dphineo_dalpha

      integer :: i, j
      integer :: ialpha, iz, is

      real, dimension(:), allocatable :: zed_stella
      integer :: nfp, nfp_stella
      integer, dimension(:), allocatable :: nzed_per_field_period
      real, dimension(:, :), allocatable :: zed_stella_by_field_period
      real, dimension(:, :), allocatable :: alpha_like_stella
      integer, dimension(:, :), allocatable :: alpha_sort_map

      integer :: nzed_sfincs, nalpha_sfincs
      real, dimension(:), allocatable :: zed_sfincs, alpha_like_sfincs
      real, dimension(:, :), allocatable :: phi_sfincs

      real, dimension(:, :), allocatable :: phi_stella_zgrid, phi_stella
      real, dimension(:, :), allocatable :: dphi_dalpha_stella_zgrid
      real, dimension(:, :), allocatable :: tmp_sfincs, tmp_stella_zgrid, tmp_stella
      real, dimension(:, :), allocatable :: dh_stella_zgrid
      real, dimension(:, :, :, :), allocatable :: h_stella, dh_stella

      ! zed coordinate in stella is zeta when simulating stellarators (using vmec)
      ! and theta otherwise.  this leads to some complications, treated below

      allocate (zed_stella(nz2pi))
      allocate (alpha_like_stella(nalpha, nz2pi))
      allocate (alpha_sort_map(nalpha, nz2pi))
      ! obtain theta and zeta grids used in stella, and assign them
      ! to the zed and alpha-like coordinates.
      ! for stellarator, zed=zeta and alpha_like=theta.
      ! for tokamak, zed=theta and alpha_like = zeta.
      call get_stella_theta_zeta_grids(alpha_like_stella, zed_stella)
!    do iz = 1, nz2pi
!       do ialpha = 1, size(alpha_like_stella,1)
!          write (*,*) 'unsorted_alpha_like_stella', alpha_like_stella(ialpha,iz)
!       end do
!       write (*,*)
!    end do
      ! rearrange alpha_like_stella to be in ascending order and store map
      ! so that sorting can be undone later
      do iz = 1, nz2pi
         call sort_array_ascending(alpha_like_stella(:, iz), alpha_sort_map(:, iz))
      end do
!    do iz = 1, nz2pi
!       do ialpha = 1, size(alpha_like_stella,1)
!          write (*,*) 'sorted_alpha_like_stella', alpha_like_stella(ialpha,iz)
!       end do
!       write (*,*)
!    end do
      ! obtain the number of alpha-like and zed grid points to use in sfincs theta-zeta grid
      ! also obtain the number of field periods per 2*pi segment in zed
      call get_sfincs_theta_zeta_grid_sizes(nalpha_sfincs, nzed_sfincs, nfp)
      ! obtain the alpha-like and zed coordinate grids
      ! note that additional points are added at periodic points
      ! that are not sampled in sfincs
      allocate (zed_sfincs(nzed_sfincs))
      allocate (alpha_like_sfincs(nalpha_sfincs))
      call get_sfincs_theta_zeta_grids(alpha_like_sfincs, zed_sfincs)

      allocate (phi_sfincs(nalpha_sfincs, nzed_sfincs))
      call get_sfincs_field_theta_zeta(Phi1Hat, phi_sfincs)

      ! this is the number of field periods included in stella
      ! simulation domain
      nfp_stella = int((zed_stella(nz2pi) * nfp - 100.*epsilon(0.)) / (2.*pi)) + 1

      ! obtain the number of zed grid points in each field period within the zed domain
      allocate (nzed_per_field_period(nfp_stella))
      call get_nzed_per_field_period(zed_stella, nfp, nzed_per_field_period)
      ! obtain the zed grid within each field period
      allocate (zed_stella_by_field_period(maxval(nzed_per_field_period), nfp_stella))
      call sort_zed_by_field_period(zed_stella, nzed_per_field_period, nfp, zed_stella_by_field_period)

      ! interpolate phi from sfincs zed grid to stella zed grid
      allocate (phi_stella_zgrid(nalpha_sfincs, nz2pi))
      call get_field_on_stella_zed_grid(phi_sfincs, nfp_stella, nfp, nzed_per_field_period, &
                                        nalpha_sfincs, zed_sfincs, zed_stella_by_field_period, phi_stella_zgrid)

      allocate (phi_stella(nalpha, nz2pi))
      ! interpolate onto stella (sorted) alpha grid
      call get_field_stella(phi_stella_zgrid, alpha_like_sfincs, alpha_like_stella, phi_stella)
      ! need to remap from ascending (sorted) alpha grid to original ordering
      do iz = 1, nz2pi
         call unsort_array_ascending(phi_stella(:, iz), alpha_sort_map(:, iz))
      end do
      call get_field_on_extended_zed(phi_stella, phi_neoclassical)

      if (present(dphineo_dalpha)) then
         allocate (dphi_dalpha_stella_zgrid(nalpha_sfincs, nz2pi))
         call get_dfield_dalpha(phi_stella_zgrid, alpha_like_sfincs, dphi_dalpha_stella_zgrid)
         call get_field_stella(dphi_dalpha_stella_zgrid, alpha_like_sfincs, alpha_like_stella, phi_stella)
         do iz = 1, nz2pi
            call unsort_array_ascending(phi_stella(:, iz), alpha_sort_map(:, iz))
         end do
         call get_field_on_extended_zed(phi_stella, dphineo_dalpha)
         deallocate (dphi_dalpha_stella_zgrid)
      end if

      deallocate (phi_stella, phi_stella_zgrid, phi_sfincs)

      allocate (tmp_sfincs(nalpha_sfincs, nzed_sfincs))
      allocate (tmp_stella_zgrid(nalpha_sfincs, nz2pi))
      allocate (tmp_stella(nalpha, nz2pi))
      allocate (h_stella(nalpha, nz2pi, size(h_sfincs, 4), size(h_sfincs, 5)))
      if (present(dfneo_dalpha)) then
         allocate (dh_stella_zgrid(nalpha_sfincs, nz2pi))
         allocate (dh_stella(nalpha, nz2pi, size(h_sfincs, 4), size(h_sfincs, 5)))
      end if

      do is = 1, nspec
         do i = 1, size(h_sfincs, 5)
            do j = 1, size(h_sfincs, 4)
               ! re-order theta and zeta indices for sfincs h to ensure alpha-like coordinate
               ! appears before zed coordinate
               call get_sfincs_field_theta_zeta(h_sfincs(is, :, :, j, i), tmp_sfincs)
               ! interpolate onto stella zed grid
               call get_field_on_stella_zed_grid(tmp_sfincs, nfp_stella, nfp, nzed_per_field_period, &
                                                 nalpha_sfincs, zed_sfincs, zed_stella_by_field_period, tmp_stella_zgrid)
               ! interpolate onto (sorted) stella alpha-like coordinate
               call get_field_stella(tmp_stella_zgrid, alpha_like_sfincs, alpha_like_stella, tmp_stella)
               do iz = 1, nz2pi
                  call unsort_array_ascending(tmp_stella(:, iz), alpha_sort_map(:, iz))
               end do
               ! use periodicity to copy onto extended zed grid if nperiod > 1
               call get_field_on_extended_zed(tmp_stella, h_stella(:, :, j, i))

               if (present(dfneo_dalpha)) then
                  call get_dfield_dalpha(tmp_stella_zgrid, alpha_like_sfincs, dh_stella_zgrid)
                  call get_field_stella(dh_stella_zgrid, alpha_like_sfincs, &
                                        alpha_like_stella, tmp_stella)
                  do iz = 1, nz2pi
                     call unsort_array_ascending(tmp_stella(:, iz), alpha_sort_map(:, iz))
                  end do
                  call get_field_on_extended_zed(tmp_stella, dh_stella(:, :, j, i))
               end if
            end do
         end do

         do ialpha = 1, nalpha
            do iz = -nzgrid, nzgrid
               call sfincs_vspace_to_stella_vspace(ialpha, iz, is, h_stella(ialpha, iz + nzgrid + 1, :, :), &
                                                   phi_neoclassical(ialpha, iz), f_neoclassical(ialpha, iz, :, :, is))
               if (present(dfneo_dalpha)) &
                  call sfincs_vspace_to_stella_vspace(ialpha, iz, is, &
                                                      dh_stella(ialpha, iz + nzgrid + 1, :, :), dphineo_dalpha(ialpha, iz), &
                                                      dfneo_dalpha(ialpha, iz, :, :, is))
            end do
         end do

      end do

      deallocate (tmp_sfincs, tmp_stella_zgrid, tmp_stella)
      deallocate (zed_stella, alpha_like_stella)
      deallocate (alpha_sort_map)
      deallocate (zed_sfincs, alpha_like_sfincs)
      deallocate (nzed_per_field_period, zed_stella_by_field_period)
      deallocate (h_stella)
      if (present(dfneo_dalpha)) deallocate (dh_stella, dh_stella_zgrid)

   end subroutine get_sfincs_output

   subroutine get_stella_theta_zeta_grids(alpha_like_stella, zed_stella)

      use constants, only: pi
      use zgrid, only: nz2pi, zed
      use geometry, only: dzetadz
      use geometry, only: alpha
      use parameters_kxky_grids, only: nalpha
      use globalVariables, only: iota

      implicit none

      real, dimension(:, :), intent(out) :: alpha_like_stella
      real, dimension(:), intent(out) :: zed_stella

      integer :: nzpi

      nzpi = nz2pi / 2

      ! convert from scaled zed grid on [-pi,pi]
      ! to un-scaled grid with lower bound of zero
      ! note that dzetadz=1 unless geometryScheme=5 (VMEC)
      ! for geometryScheme=5, this will get extended zeta domain
      ! with lower bound of 0
      zed_stella = (zed(-nzpi:nzpi) + pi) / dzetadz

      ! if geometryScheme is 5, then using vmec geo
      ! and thus zed in stella is scaled zeta
      ! otherwise zed in stella is theta
      if (geometryScheme == 5) then
         ! alpha_like coordinate is theta = alpha + iota*zeta
         alpha_like_stella = spread(alpha, 2, nz2pi) + spread(iota * zed_stella, 1, nalpha)
      else
         ! alpha_like coordinate is zeta = (theta-alpha)/iota
         alpha_like_stella = (spread(zed_stella, 1, nalpha) - spread(alpha, 2, nz2pi)) / iota
      end if

      ! restrict alpha to [0,2*pi]
      alpha_like_stella = modulo(alpha_like_stella, 2.*pi)

   end subroutine get_stella_theta_zeta_grids

   subroutine get_nzed_per_field_period(zed_stella, nfp, nzed_per_field_period)

      use constants, only: pi
      use zgrid, only: nz2pi

      implicit none

      real, dimension(:), intent(in) :: zed_stella
      integer, intent(in) :: nfp
      integer, dimension(:), intent(out)  :: nzed_per_field_period

      integer :: ifp, iz

      nzed_per_field_period = 0

      ifp = 1; iz = 1
      do while (iz <= nz2pi)
         if (zed_stella(iz) <= ifp * 2.*pi / nfp) then
            nzed_per_field_period(ifp) = nzed_per_field_period(ifp) + 1
            iz = iz + 1
         else
            ifp = ifp + 1
         end if
      end do

   end subroutine get_nzed_per_field_period

   subroutine sort_zed_by_field_period(zed_ext, nzed_per_fp, nfp, zed_by_fp)

      use constants, only: pi

      implicit none

      real, dimension(:), intent(in) :: zed_ext
      integer, dimension(:), intent(in) :: nzed_per_fp
      integer, intent(in) :: nfp
      real, dimension(:, :), intent(out) :: zed_by_fp

      integer :: ifp, llim, ulim

      ulim = 0
      do ifp = 1, size(zed_by_fp, 2)
         llim = ulim + 1
         ulim = llim + nzed_per_fp(ifp) - 1
         zed_by_fp(:nzed_per_fp(ifp), ifp) = modulo(zed_ext(llim:ulim) - 100.*epsilon(0.), 2.*pi / nfp)
      end do

      ! avoid special case of setting zed = 0 to zed=2*pi/nfp
      zed_by_fp(1, 1) = 0.0

   end subroutine sort_zed_by_field_period

   subroutine get_sfincs_theta_zeta_grid_sizes(nalpha_sfincs, nzed_sfincs, nfp)

      use constants, only: pi
      use export_f, only: export_f_zeta

      implicit none

      integer, intent(out) :: nalpha_sfincs, nzed_sfincs, nfp

      ! if geometryScheme is 5, then using vmec geo
      ! and thus zed in stella is scaled zeta
      ! otherwise zed in stella is theta
      if (geometryScheme == 5) then
         ! note that zeta grid in sfincs only covers one field period of the stellarator
         ! get the number of field periods from the sfincs zeta grid
         nfp = nint(2.*pi / (export_f_zeta(nzeta) + export_f_zeta(2)))
         ! zed coordinate is zeta
         nzed_sfincs = nzeta + 1
         ! alpha_like coordinate is theta = alpha + iota*zeta
         nalpha_sfincs = ntheta + 1
      else
         nfp = 1
         ! zed coordinate is theta
         nzed_sfincs = ntheta + 1
         ! alpha_like coordinate is zeta = (theta-alpha)*q
         nalpha_sfincs = nzeta
      end if

   end subroutine get_sfincs_theta_zeta_grid_sizes

   subroutine get_sfincs_theta_zeta_grids(alpha_like_sfincs, zed_sfincs)

      use export_f, only: export_f_theta, export_f_zeta

      implicit none

      real, dimension(:), intent(out) :: zed_sfincs, alpha_like_sfincs

      if (geometryScheme == 5) then
         ! zed is zeta.  it goes from 0 to 2*pi/nfp - dzeta in sfincs,
         ! where nfp is the number of field periods
         zed_sfincs(:nzeta) = export_f_zeta(:nzeta)
         ! add in point at 2*pi/nfp using periodicity
         zed_sfincs(nzeta + 1) = zed_sfincs(nzeta) + zed_sfincs(2)

         ! alpha-like coordinate is theta.  it goes from 0 to 2*pi-dtheta in sfincs
         alpha_like_sfincs(:ntheta) = export_f_theta(:ntheta)
         ! add in point at 2*pi using periodicity
         alpha_like_sfincs(ntheta + 1) = export_f_theta(ntheta) + export_f_theta(2)
      else
         ! zed is theta.  it goes from 0 to 2*pi-dtheta in sfincs
         zed_sfincs(:ntheta) = export_f_theta(:ntheta)
         ! add in point at 2*pi using periodicity
         zed_sfincs(ntheta + 1) = export_f_theta(ntheta) + export_f_theta(2)

         ! alpha-like coordinate is zeta.  there should only be 1 zeta for tokamak calculation
         alpha_like_sfincs(:nzeta) = export_f_zeta(:nzeta)
      end if

   end subroutine get_sfincs_theta_zeta_grids

   subroutine get_sfincs_field_theta_zeta(field_in, field_out)

      implicit none

      real, dimension(:, :), intent(in) :: field_in
      real, dimension(:, :), intent(out) :: field_out

      integer :: itheta, izeta

      ! want phi from sfincs such that alpha-like coordinate
      ! appears in first index and zed coordinate in second index
      if (geometryScheme == 5) then
         ! get phi on sfincs (theta,zeta) grid
         field_out(:ntheta, :nzeta) = field_in
         ! use periodicity in theta to add point at theta=2*pi
         field_out(ntheta + 1, :nzeta) = field_out(1, :nzeta)
         ! use periodicity in zeta to add point at zeta = 2*pi/nfp
         field_out(:, nzeta + 1) = field_out(:, 1)
      else
         ! get phi on sfincs (zeta,theta) grid
         do izeta = 1, nzeta
            do itheta = 1, ntheta
               field_out(izeta, itheta) = field_in(itheta, izeta)
            end do
         end do
         ! use periodicity in theta to add point at 2*pi
         field_out(:, ntheta + 1) = field_out(:, 1)
      end if

   end subroutine get_sfincs_field_theta_zeta

   subroutine get_field_on_stella_zed_grid(field_sfincs, nfp_stella, nfp, nzed_per_field_period, &
                                           nalpha_sfincs, zed_sfincs, zed_stella_by_field_period, field_stella_zgrid)

      use constants, only: pi
      use splines, only: linear_interp_periodic

      implicit none

      real, dimension(:, :), intent(in) :: field_sfincs
      integer, intent(in) :: nfp_stella, nfp, nalpha_sfincs
      integer, dimension(:), intent(in) :: nzed_per_field_period
      real, dimension(:), intent(in) ::  zed_sfincs
      real, dimension(:, :), intent(in) :: zed_stella_by_field_period
      real, dimension(:, :), intent(out) :: field_stella_zgrid

      integer :: ialpha, ifp
      integer :: llim, ulim

      do ialpha = 1, nalpha_sfincs
         ulim = 0
         do ifp = 1, nfp_stella
            llim = ulim + 1
            ulim = llim + nzed_per_field_period(ifp) - 1
            call linear_interp_periodic(zed_sfincs, field_sfincs(ialpha, :), &
                                        zed_stella_by_field_period(:nzed_per_field_period(ifp), ifp), &
                                        field_stella_zgrid(ialpha, llim:ulim), 2.*pi / nfp)
         end do
      end do

   end subroutine get_field_on_stella_zed_grid

   subroutine get_field_stella(field_stella_zgrid, alpha_like_sfincs, alpha_like_stella, field_stella)

      use splines, only: linear_interp_periodic
      use zgrid, only: nz2pi

      implicit none

      real, dimension(:, :), intent(in) :: field_stella_zgrid
      real, dimension(:), intent(in) :: alpha_like_sfincs
      real, dimension(:, :), intent(in) :: alpha_like_stella
      real, dimension(:, :), intent(out) :: field_stella

      integer :: iz

      do iz = 1, nz2pi
         call linear_interp_periodic(alpha_like_sfincs, field_stella_zgrid(:, iz), &
                                     alpha_like_stella(:, iz), field_stella(:, iz))
      end do

   end subroutine get_field_stella

   subroutine get_field_on_extended_zed(field_stella, field_neoclassical)

      use zgrid, only: nzgrid, nz2pi, nperiod
      use parameters_kxky_grids, only: nalpha

      implicit none

      real, dimension(:, :), intent(in) :: field_stella
      real, dimension(:, -nzgrid:), intent(out) :: field_neoclassical

      integer :: ialpha
      integer :: iz_low, iz_up
      integer :: ip

      ! need to account for cases with nperiod > 1
      do ialpha = 1, nalpha
         iz_low = -nzgrid
         iz_up = -nzgrid + nz2pi - 1
         field_neoclassical(ialpha, iz_low:iz_up) = field_stella(ialpha, :)
         ! if nperiod > 1 need to make copies of
         ! neoclassical potential for other 2pi segments
         if (nperiod > 1) then
            do ip = 2, 2 * nperiod - 1
               iz_low = iz_up + 1
               iz_up = iz_low + nz2pi - 2
               field_neoclassical(ialpha, iz_low:iz_up) = field_stella(ialpha, 2:)
            end do
         end if
      end do

   end subroutine get_field_on_extended_zed

   subroutine sfincs_vspace_to_stella_vspace(ialpha, iz, is, h_stella, phi_neoclassical, f_neoclassical)

      use constants, only: pi
      use species, only: spec
      use vpamu_grids, only: nvpa, nvgrid, nmu
      use vpamu_grids, only: vpa, vperp2
      use vpamu_grids, only: maxwell_mu, maxwell_vpa
      use globalVariables, only: nxi_sfincs => nxi
      use globalVariables, only: nx_sfincs => nx
      use globalVariables, only: x_sfincs => x
      use xGrid, only: xGrid_k

      implicit none

      integer, intent(in) :: ialpha, iz, is
      real, dimension(:, :), intent(in) :: h_stella
      real, intent(in) :: phi_neoclassical
      real, dimension(:, :), intent(out) :: f_neoclassical

      integer :: iv, imu, ixi
      real, dimension(1) :: x_stella
      integer, dimension(2) :: sgnvpa
      integer :: nxi_stella
      real, dimension(:), allocatable :: xi_stella, hstella
      real, dimension(:, :), allocatable :: xsfincs_to_xstella, legpoly
      real, dimension(:), allocatable :: htmp

      ! each (vpa,mu) pair in stella specifies a speed
      ! on each speed arc, there are two (vpa,mu) pairs:
      ! one each corresponding to a given +/- vpa
      nxi_stella = 2
      allocate (xi_stella(nxi_stella))
      allocate (hstella(nxi_stella))
      allocate (legpoly(nxi_stella, 0:nxi_sfincs - 1))

      allocate (htmp(nxi_sfincs))
      allocate (xsfincs_to_xstella(1, nx_sfincs))

      sgnvpa(1) = 1; sgnvpa(2) = -1

      ! h_stella is on the sfincs energy grid
      ! but is spectral in pitch-angle
      do imu = 1, nmu
         ! loop over positive vpa values
         ! negative vpa values will be treated inside loop using symmetry
         do iv = nvgrid + 1, nvpa
            ! x_stella is the speed
            ! corresponding to this (vpa,mu) grid point
            x_stella = sqrt(vpa(iv)**2 + vperp2(ialpha, iz, imu))
            ! xi_stella contains the two pitch angles (+/-)vpa/v
            ! corresponding to this (vpa,mu) grid point
            xi_stella = sgnvpa * vpa(iv) / x_stella(1)

            ! set up matrix that interpolates from sfincs speed grid
            ! to the speed corresponding to this (vpa,mu) grid point
            call polynomialInterpolationMatrix(nx_sfincs, 1, &
                                               x_sfincs, x_stella, exp(-x_sfincs * x_sfincs) * (x_sfincs**xGrid_k), &
                                               exp(-x_stella * x_stella) * (x_stella**xGrid_k), xsfincs_to_xstella)

            ! do the interpolation
            do ixi = 1, nxi_sfincs
               htmp(ixi) = sum(h_stella(ixi, :) * xsfincs_to_xstella(1, :))
            end do

            ! next need to Legendre transform in pitch-angle
            ! first evaluate Legendre polynomials at requested pitch angles
            call legendre(xi_stella, legpoly)

            ! then do the transforms
            call legendre_transform(legpoly, htmp, hstella)

            f_neoclassical(nvpa - iv + 1, imu) = hstella(2)
            f_neoclassical(iv, imu) = hstella(1)

         end do
         ! h_sfincs is H_nc / (nref/vt_ref^3), with H_nc the non-Boltzmann part of F_nc
         ! NB: n_ref, etc. is fixed in stella to be the reference density
         ! at the central sfincs simulation; i.e., it does not vary with radius
         ! to be consistent with stella distribution functions,
         ! want H_nc / (n_s / vt_s^3 * pi^(3/2)) / exp(-v^2/vts^2)
         f_neoclassical(:, imu) = f_neoclassical(:, imu) &
                                  * pi**1.5 * spec(is)%stm**3 / spec(is)%dens

         ! phi_sfincs is e phi / Tref as long as alpha=1 (default)
         ! need to multiply by Z_s * Tref/T_s
         f_neoclassical(:, imu) = f_neoclassical(:, imu) &
                                  - phi_neoclassical * spec(is)%zt * maxwell_vpa * maxwell_mu(ialpha, iz, imu)
      end do

      deallocate (xi_stella, hstella, legpoly)
      deallocate (xsfincs_to_xstella)
      deallocate (htmp)

   end subroutine sfincs_vspace_to_stella_vspace

   subroutine get_dfield_dalpha(field, alpha_like_sfincs, dfield_dalpha)

      use zgrid, only: nz2pi

      implicit none

      real, dimension(:, :), intent(in) :: field
      real, dimension(:), intent(in) :: alpha_like_sfincs
      real, dimension(:, :), intent(out) :: dfield_dalpha

      integer :: iz

      do iz = 1, nz2pi
         call get_periodic_derivative(field(:, iz), alpha_like_sfincs, dfield_dalpha(:, iz))
      end do

   end subroutine get_dfield_dalpha

   ! 4th order accurate, centered differences, assumes
   ! first and last elements of f are equal (periodic)
   subroutine get_periodic_derivative(f, x, dfdx)

      implicit none

      real, dimension(:), intent(in) :: f, x
      real, dimension(:), intent(out) :: dfdx

      integer :: i, n
      real, dimension(:), allocatable :: fp

      n = size(x)

      allocate (fp(-1:n + 2))
      ! extend f using periodicity, as these additional points
      ! required near boundaries for finite differences
      fp(1:n) = f
      fp(-1:0) = f(n - 1:n)
      fp(n + 1:n + 2) = f(1:2)

      do i = 1, n
         dfdx(i) = (0.25 * (fp(i - 2) - fp(i + 2)) + 2.0 * (fp(i + 1) - fp(i - 1))) / 3.0
      end do

      deallocate (fp)

   end subroutine get_periodic_derivative

!   subroutine bilinear_interpolation (alpha_in, zed_in, phi_in, alpha_out, zed_out, phi_out)

!     use zgrid, only: nz2pi

!     implicit none

!     real, dimension (:,:), intent (in) :: alpha_sfincs, phi_sfincs
!     real, dimension (:), intent (in) :: zed_sfincs
!     real, dimension (:), intent (in) :: alpha_stella, zed_stella
!     real, dimension (:,:), intent (out) :: phi_stella

!     integer :: ialpha, iz

!     do ialpha = 1, nalpha
!        do iz = 1, nz2pi
!           call find_sfincs_cell (alpha_stella(ialpha), zed_stella(iz), alpha_sfincs, zed_sfincs, &
!                alpha_grids, zed_grids)
!        end do
!     end do

!   contains

!     subroutine find_sfincs_cell (alpha_target, zed_target, alpha, zed, ialpha_out, ized_out)

!       use zgrid, only: nz2pi

!       implicit none

!       real, intent (in) :: alpha_target, zed_target
!       real, dimension (:,:), intent (in) :: alpha
!       real, dimension (:), intent (in) :: zed
!       integer, dimension (2), intent (out) :: ialpha_out, ized_out

!       integer :: iz, ia

!       do iz = 1, nz2pi
!          do ia = 1, size(alpha,1)
!             if (
!          end do
!       end do

!     end subroutine find_sfincs_cell

!   end subroutine bilinear_interpolation

   ! returns the Legendre polynomials (legp)
   ! on requested grid (x)
   subroutine legendre(x, legp)

      implicit none

      real, dimension(:), intent(in) :: x
      real, dimension(:, 0:), intent(out) :: legp

      integer :: n, idx

      n = size(legp, 2) - 1

      legp(:, 0) = 1.0
      legp(:, 1) = x

      do idx = 2, n
         legp(:, idx) = ((2.*idx - 1.) * x * legp(:, idx - 1) + (1.-idx) * legp(:, idx - 2)) / idx
      end do

   end subroutine legendre

   subroutine legendre_transform(legp, coefs, func)

      implicit none

      real, dimension(:, 0:), intent(in) :: legp
      real, dimension(:), intent(in) :: coefs
      real, dimension(:), intent(out) :: func

      integer :: i

      func = 0.
      do i = 1, size(coefs)
         func = func + legp(:, i - 1) * coefs(i)
      end do

   end subroutine legendre_transform

   subroutine broadcast_sfincs_output(fneo, phineo)
      use mp, only: broadcast
      use zgrid, only: nzgrid
      implicit none
      real, dimension(-nzgrid:, :, :, :, :), intent(in out) :: fneo
      real, dimension(-nzgrid:, :), intent(in out) :: phineo
      call broadcast(fneo)
      call broadcast(phineo)
   end subroutine broadcast_sfincs_output

   subroutine write_sfincs(irad, nrad_max)

      use species, only: nspec
      use globalVariables, only: Phi1Hat
      use export_f, only: export_f_zeta, export_f_theta
      use export_f, only: delta_f

      implicit none

      integer, intent(in) :: irad, nrad_max

      integer :: unit = 999
      integer :: izeta, itheta, is, i, j
      character(1) :: irad_str

      write (irad_str, '(i0)') irad + min(nrad_max, 1)
      open (unit=unit, file='sfincs.output.rad'//irad_str, status='replace', action='write')

      do izeta = 1, nzeta
         do itheta = 1, ntheta
            write (unit, '(a8,3e13.5,i3)') 'Phi1Hat', export_f_theta(itheta), export_f_zeta(izeta), Phi1Hat(itheta, izeta), irad
         end do
         write (unit, *)
      end do
      write (unit, *)

      do is = 1, nspec
         do i = 1, size(delta_f, 5)
            do j = 1, size(delta_f, 4)
               do izeta = 1, nzeta
                  do itheta = 1, ntheta
                     write (unit, '(a8,3e13.5,4i5)') 'delta_f', export_f_theta(itheta), export_f_zeta(izeta), &
                        delta_f(is, itheta, izeta, j, i), i, j, is, irad
                  end do
                  write (unit, *)
               end do
            end do
         end do
      end do
      write (unit, *)

      close (unit)

   end subroutine write_sfincs

   subroutine read_sfincs_output(irad, nrad_max)

      use species, only: nspec
      use globalVariables, only: Phi1Hat
      use export_f, only: export_f_zeta, export_f_theta
      use export_f, only: delta_f

      implicit none

      integer, intent(in) :: irad, nrad_max

      integer :: unit = 999
      integer :: izeta, itheta, is, i, j
      character(8) :: dum
      character(1) :: irad_str

      write (irad_str, '(i0)') irad + nrad_max
      open (unit=unit, file='sfincs.output.rad'//irad_str, status='old', action='read')

      do izeta = 1, nzeta
         do itheta = 1, ntheta
            read (unit, *) dum, export_f_theta(itheta), export_f_zeta(izeta), Phi1Hat(itheta, izeta), dum
         end do
         read (unit, *)
      end do
      read (unit, *)

      do is = 1, nspec
         do i = 1, size(delta_f, 5)
            do j = 1, size(delta_f, 4)
               do izeta = 1, nzeta
                  do itheta = 1, ntheta
                     read (unit, *) dum, export_f_theta(itheta), export_f_zeta(izeta), &
                        delta_f(is, itheta, izeta, j, i), dum, dum, dum, dum
                  end do
                  read (unit, *)
               end do
            end do
         end do
      end do
      read (unit, *)

      close (unit)

   end subroutine read_sfincs_output

# endif

end module sfincs_interface
