eikcoefs Subroutine

public subroutine eikcoefs(ntheta_input, nperiod, outputs, job_id_in)

FIXME : Add documentation

Arguments

Type IntentOptional Attributes Name
integer, intent(inout) :: ntheta_input

Sets the number of theta grid points. If odd will be updated to be ntheta_input - 1

integer, intent(in) :: nperiod
type(eikcoefs_output_type), intent(out) :: outputs

A type containing all of the outputs of eikcoefs

integer, intent(in), optional :: job_id_in

job_id is the job id in a Trinity or list mode run. Only for debug


Contents

Source Code


Source Code

  subroutine eikcoefs (ntheta_input, nperiod, outputs, job_id_in)
    use constants, only: twopi
    use file_utils, only: open_output_file, close_output_file, finish_file_utils
    use file_utils, only: futils_initialized => initialized, init_file_utils
    use optionals, only: get_option_with_default
    implicit none
    !> Sets the number of theta grid points. If odd will be updated to be ntheta_input - 1
    integer, intent(in out) :: ntheta_input
    integer, intent(in) :: nperiod
    !> A type containing all of the outputs of eikcoefs
    type(eikcoefs_output_type), intent(out) :: outputs
    !> job_id is the job id in a Trinity or list mode run. Only for debug
    integer, intent(in), optional :: job_id_in
    logical :: start_file_utils, dummy_flag
    integer :: job_id, i, iBmax, fort11_unit
    real, allocatable, dimension(:) :: seik, rgrid, invR, trip, dsdthet, rmajor, dummy
    real, allocatable, dimension(:,:) :: thgrad, rpgrad, crpgrad, grads
    real, allocatable, dimension(:,:) :: gradrptot, gradstot, gradztot, gradthtot, gradbtot
    real :: shat, dbetadrho, dqdrp, qval, dum, bi
    real :: rp, pbar, drhodrp, drhodrhod, drhodpsi, dpsidrp, dpsidrho
    ! Used to translate definition of the poloidal angle such that the location
    ! of the maximium magnetic field is at pi (see section 3.2.3 of Ball MIT
    ! Masters thesis or section 2.4 of "GS2 analytic geometry specification")
    real :: thetaShift
    job_id = get_option_with_default(job_id_in, -1)
    start_file_utils = .not. futils_initialized

    if (start_file_utils) then
      ! Initialize file utils with a run name we provide
      ! and by setting trin_run prevent it from searching
      ! the command line arguments.
      write(*,*) 'WARNING: Automatically initializing file_utils&
      & in eikcoefs. This is a hack to &
      & support older programs like eiktest. You should&
      & initialize file utils before calling eikcoefs'
      call init_file_utils(dummy_flag, name='geometry_module', trin_run=.true., input=.false.)
    end if
    !Create module level output file and get unit
    call open_output_file(fort11_unit,".fort.11")

    !Compute the initial constants
    debug = (verb > 2)
    if (debug) writelots = .true.

    if (debug) write(6,*) "eikcoefs: local_eq=",local_eq, 'jid=',job_id

    if (local_eq) then
       if (irho /= 2) write(*,*) 'Forcing irho = 2 for local_eq'
       irho = 2
       if (bishop == 0) then
          write(*,*) 'Bishop = 0 not compatible with local_eq. Forcing bishop = 1'
          bishop = 1 !Not clear that this is a good choice for Miller
       end if
    end if

    if (debug) write(6,*) "eikcoefs: call check ",'jid=',job_id
    call check(gen_eq, efit_eq, ppl_eq, local_eq, dfit_eq, idfit_eq, chs_eq, gs2d_eq, transp_eq, bishop)

    call init_uniform_theta_grid(ntheta_input, outputs%theta, nperiod, ntgrid, nth)
    surf%nt = ntheta_input
    if (debug) write(6,*) 'eikcoefs: alloc_local_arrays'
    if (.not. allocated(seik)) call alloc_local_arrays(ntgrid)
    thetaShift = 0.0
    allocate(outputs%bmag(-ntgrid:ntgrid), outputs%bpol(-ntgrid:ntgrid))
    if (.not. allocated(geom)) call allocate_geom(geom)
    call initialise_geometry_and_get_bmag(geom, surf%r, thetaShift, outputs%theta, outputs%bmag, outputs%bpol)

    ! Start theta shift calc. This is used to modify the definition of
    ! theta such that the maxiumum magnetic field is at pi. This is
    ! necessary to correctly calculate bounce points. Can only shift with local_eq
    if (local_eq) then
       ! Locate the index of the maximum magnetic field in the range [-pi,pi)
       iBmax = maxloc(outputs%bmag(-nth:nth-1), dim = 1) - (nth + 1)

       ! If max is not at -pi and theta then shift theta.
       if (iBmax /= -nth) then
          thetaShift = outputs%theta(iBmax) - outputs%theta(-nth)
          call geom%finalise !Clean up geometry to allow reinitialisation
          write(*,*) "Bmax is not located at +/-pi, translating theta by",thetaShift," so that it is..."
          call initialise_geometry_and_get_bmag(geom, surf%r, thetaShift, outputs%theta, outputs%bmag, outputs%bpol)
       end if
    end if

    !     compute |grad rho|:
    allocate(outputs%grho(-ntgrid:ntgrid))
    outputs%grho(-nth:nth) = hypot(rpgrad(-nth:nth,1), rpgrad(-nth:nth,2)) * abs(drhodrp)
    call periodic_copy(outputs%grho, 0., nperiod)
    if (force_sym) call sym(outputs%grho, 0, ntgrid)

    if (writelots) write(fort11_unit,*) 'dpsidrp= ',dpsidrp
    drhodpsi = drhodrp / dpsidrp ; dpsidrho = 1.0 / drhodpsi

    !     compute derivatives of eikonal, grad B and grad Z
    if (bishop /= 0) then
       call bishop_grads(surf%r, surf%rmaj, outputs%bmag, outputs%bpol)
    else
       call bishop_large_aspect_ratio_approximation(surf%r, surf%dr, ntheta_input)
    end if

    if (writelots) then
       write(fort11_unit,*) 'i,theta,grads1,2=  '
       do i = -ntgrid, ntgrid
          write(fort11_unit,*) i, outputs%theta(i), grads(i, 1), grads(i, 2)
       end do
    end if
    !     compute total grad S including toroidal part
    call gradstottgrid(invR, grads, gradstot)

    gradrptot(-nth:nth, 1:2) = rpgrad(-nth:nth, :)
    call periodic_copy(gradrptot(:, 1), 0., nperiod)
    call periodic_copy(gradrptot(:, 2), 0., nperiod)
    gradrptot(:, 3) = 0.

    ! gradthtot is grad theta since ballooning, theta goes from
    ! -ntheta*nperiod -> ntheta*nperiod, but grad theta is periodic in theta
    gradthtot(-nth:nth, 1:2) = thgrad(-nth:nth, :)
    call periodic_copy(gradthtot(:, 1), 0., nperiod)
    call periodic_copy(gradthtot(:, 2), 0., nperiod)
    gradthtot(:, 3) = 0.

    !     compute the magnitude squared of grad S
    allocate(outputs%gds2(-ntgrid:ntgrid), outputs%gds21(-ntgrid:ntgrid), &
         outputs%gds22(-ntgrid:ntgrid), outputs%gds23(-ntgrid:ntgrid), &
         outputs%gds24(-ntgrid:ntgrid), outputs%gds24_noq(-ntgrid:ntgrid))
    call calculate_metric_terms(gradstot, gradrptot, gradthtot, dpsidrho, dqdrp, &
         qval, surf%r, outputs%bmag, force_sym, outputs%gds2, outputs%gds21, outputs%gds22, &
         outputs%gds23, outputs%gds24, outputs%gds24_noq)

    !     compute curvature, grad b, and coriolis drift terms
    allocate(outputs%gbdrift(-ntgrid:ntgrid), outputs%gbdrift0(-ntgrid:ntgrid))
    allocate(outputs%cvdrift(-ntgrid:ntgrid), outputs%cvdrift0(-ntgrid:ntgrid))
    allocate(outputs%cdrift(-ntgrid:ntgrid), outputs%cdrift0(-ntgrid:ntgrid))
    call drift(rgrid, rp, gradstot, gradrptot, gradztot, gradbtot, dqdrp, dpsidrho, dpsidrp,&
         outputs%bmag, invR, outputs%theta, dbetadrho, bi, force_sym, nperiod, &
         outputs%gbdrift, outputs%gbdrift0, outputs%cvdrift, outputs%cvdrift0, &
         outputs%cdrift, outputs%cdrift0)

    !Recompute the magnetic field along theta if use_large_aspect - should this come earlier?
    if (use_large_aspect) then
       outputs%bmag(-nth:nth)= 1 / (1 + rgrid(-nth:nth) * cos(outputs%theta(-nth:nth)) / surf%rmaj)
       call periodic_copy(outputs%bmag, 0., nperiod)
    end if
    if (force_sym) call sym(outputs%bmag, 0, ntgrid)

    !     compute the grad-parallel coefficient
    allocate(outputs%gradpar(-ntgrid:ntgrid))
    if (dfit_eq) then
       outputs%gradpar = trip / outputs%bmag
    else
       outputs%gradpar(-nth:nth) = -bi / (rmajor(-nth:nth)**2 * outputs%bmag(-nth:nth) * dsdthet(-nth:nth))
       call periodic_copy(outputs%gradpar, 0., nperiod)
    end if
    if (force_sym) call sym(outputs%gradpar, 0, ntgrid)

    ! Note : We should call this before we potentially redefine the theta grid
    ! due to equal_arc below as some specific equilibrium types rely on theta
    ! having a consistent definition with that used to calculate rgrid initially.
    ! In particular this is important for Efit and GS2D equilibria (in eeq). Other
    ! equilibrium types use interpolation so shouldn't depend on the specific theta grid
    ! (although it's possible there could still be some interaction, given that one of the
    ! outputs of the below 'aplot' is given by `seik` which is evaluated on the original
    ! grid rather than the "new" one defined below).
    ! @note Check if restriction on cartesian geo types still holds.
    allocate(outputs%rplot(-ntgrid:ntgrid), outputs%rprime(-ntgrid:ntgrid))
    allocate(outputs%zplot(-ntgrid:ntgrid), outputs%zprime(-ntgrid:ntgrid))
    allocate(outputs%aplot(-ntgrid:ntgrid), outputs%aprime(-ntgrid:ntgrid))
    call plotdata (rgrid, seik, grads, dpsidrho, surf%dr, outputs%theta, &
         outputs%rplot, outputs%zplot, outputs%aplot, &
         outputs%rprime, outputs%zprime, outputs%aprime)

    ! In the lines below we redefine theta to be theta', and gradpar
    ! to be gradpar'. This means that the values of theta in the gs2
    ! output file are not values of the poloidal angle but values of theta_prime
    if (equal_arc) then
       allocate(dummy(-ntgrid:ntgrid))
       call arclength (ntheta_input, nperiod, outputs%theta, outputs%gradpar, dummy)
       outputs%theta = dummy
    end if

    call calculate_surface_area(outputs%theta, outputs%gradpar, outputs%bmag, outputs%grho, &
         dpsidrho, surf%r, nth, outputs%surfarea, outputs%dvdrhon)
    if (writelots) write (fort11_unit,*) 'surfarea=', outputs%surfarea

    !Close output file
    call close_output_file(fort11_unit)

    outputs%kxfac = abs(qval * dpsidrho / surf%r)

    if (qval == 0. .and. (dfit_eq .or. idfit_eq)) then
       shat = max(1.e-7, shat)
       !RN> I think kxfac should include the sign. However, from GS2 convention,
       !RN> eventually kxfac = abs(kxfac). See get_grids in theta_grid. Even if the sign is
       !RN> effective, probably it has no effect due to symmetry. I just leave the sign here.
       outputs%kxfac = - dpsidrho
       outputs%gds22 = (outputs%grho * dpsidrho)**2
    end if

    if (start_file_utils) call finish_file_utils

    ! Pack data into output type
    outputs%ntheta = ntheta_input ; outputs%drhodpsi = drhodpsi ; outputs%rhoc = surf%r
    outputs%qsf = qval ; outputs%shat = shat ; outputs%dbetadrho = dbetadrho

  contains
    !> Initialises the specific geometry backend
    subroutine initialise_geometry_and_get_bmag(geom, rhoc, thetaShift, theta, bmag, bpol)
      use geo_utils, only: geo_input_type
      class(abstract_geo_type), intent(in out) :: geom
      real, intent(in) :: rhoc, thetaShift
      real, dimension(-ntgrid:), intent(in) :: theta
      real, dimension(-ntgrid:), intent(out) :: bmag, bpol
      type(geo_input_type) :: geom_inputs
      real :: avgrmid, rbar, B_T0
      integer :: fort78_unit, fort99_unit
      ! Initialise the specific geometry types (mostly reading data
      ! from file and calculating derived quantities).
      geom_inputs%surf = surf ; geom_inputs%thShift = thetaShift ; geom_inputs%big = big
      geom_inputs%eqfile = eqfile ; geom_inputs%eqnormfile = eqnormfile ; geom_inputs%theta = theta
      call geom%initialise(geom_inputs, rpmin, rpmax, surf%rmaj, B_T0, avgrmid)
      if (debug) write(6,*) "eikcoefs: rpmin,rpmax=", rpmin, rpmax
      if (.not. local_eq) surf%r_geo = geom%rcenter(rpmax) !Only used for writelots/checks
      if (surf%r_geo > surf%rmaj + 1.e-5) then
         write(6,*)
         write(6,*) 'Warning: the center of the LCFS is outboard of the '
         write(6,*) 'center of the reference surface.  This is probably wrong.'
         write(6,*)
      end if

      if (writelots) then
         write(fort11_unit,*) 'All lengths normalized to the avg midplane minor radius'
         if (irho /= 2) then
            write(fort11_unit,*) 'except rho, which always goes from 0 to 1.'
            write(fort11_unit,*) 'The rho normalization is determined by irho.'
         end if
         write(fort11_unit,*) 'avgrmid = ',avgrmid,' meters'
         if (.not. local_eq) write(fort11_unit,*) 'R_mag = ',surf%rmaj
         if (.not. local_eq) write(fort11_unit,*) 'B_T0 = ',B_T0
      end if

      if (debug) write(6,*) "eikcoefs: find rp ", 'jid=',job_id
      rp = rpofrho(rhoc) ; drhodrp = drho_drp(rp, surf%dr)
      if (debug) write(6,*) "eikcoefs: rp=",rp, 'jid=',job_id

      if (debug) write(6,*) "eikcoefs: find rgrid"
      call rtgrid(rgrid, rp, theta)
      call periodic_copy(rgrid, 0.0, nperiod)

      pbar = pbarfun(rgrid(0), theta(0))
      bi = geom%btori(pbar) ! Get the toroidal field function (i.e. fp = R B_T)

      if (writelots) then
         write(fort11_unit,*) 'rp of rho = ',rp
         write(fort11_unit,*) 'Rcenter(rho) = ',geom%rcenter(rp)
         write(fort11_unit,*) 'R_geo = ',surf%r_geo,' !or ',bi

         call open_output_file(fort78_unit,".fort.78")
         call open_output_file(fort99_unit,".fort.99")
         write(fort78_unit,*) '# theta_pol,  R,   Z'
         do i = -nth, nth
            write(fort78_unit,*) theta(i), geom%Rpos(rgrid(i), theta(i)), &
                 geom%Zpos(rgrid(i), theta(i))
            write(fort99_unit,*) 'i=', i, ' thet= ', theta(i), ' rgrid= ', rgrid(i)
         end do
         call close_output_file(fort78_unit) ; call close_output_file(fort99_unit)

         rbar = 0.
         do i = -nth, nth - 1
            rbar = rbar + rgrid(i) * (theta(i + 1) - theta(i))
         end do
         rbar = rbar / twopi
         write (*,*) 'r_bar = ', rbar
         write (*,*)

         write(fort11_unit,*) 'drhodrp=', drhodrp
         if (irho == 1) then
            drhodrhod = drho_drhod(rp, drhodrp, surf%dr)
            write(fort11_unit,*) 'drho_drhod=', drhodrhod
         end if
      end if

      if (debug) write(fort11_unit,*) 'eikcoefs: various logs', gen_eq, ppl_eq, chs_eq, efit_eq, dfit_eq, idfit_eq, gs2d_eq, transp_eq

      if (debug) write(6,*) "eikcoefs: call rmajortgrid"
      call rmajortgrid(rgrid, theta, rmajor) ; call periodic_copy(rmajor, 0., nperiod)
      if (use_large_aspect) then
         invR = 1.0 / surf%rmaj
      else
         call invRtgrid(rgrid, theta, invR) ; call periodic_copy(invR, 0., nperiod)
      end if

      !     compute gradient of rp in R, Z, phi (cartesian - crpgrad)
      ! Note that rpgrad is given in (R,Z,phi) coordinates if bishop == 0
      ! and (rho,l,phi) coordinates otherwise
      ! Also, rpgrad is grad psi if not local_eq and grad rho if local_eq -- MAB
      ! When bishop > 0
      ! rpgrad(:,1) is |grad rp| ; rpgrad(:,2) is 0
      ! crpgrad(:,1) is drp/dR ; crpgrad(:,2) is drp/dZ
      !
      ! bishop = 0 is high aspect ratio coefficients for debugging -- EGH
      call geom%gradient(rgrid, theta, crpgrad, 'P', dum, nth, ntgrid, .false.)
      call periodic_copy(crpgrad(:,1), 0., nperiod)
      call periodic_copy(crpgrad(:,2), 0., nperiod)

      !     compute coordinate gradients
      call geom%gradient(rgrid, theta, rpgrad, 'P', dum, nth, ntgrid, bishop /= 0)
      call geom%gradient(rgrid, theta, thgrad, 'T', dum, nth, ntgrid, bishop /= 0)
      call tripprod2dtgrid(rpgrad, thgrad, invR, trip)
      call periodic_copy(trip, 0., nperiod)

      !     compute eikonal S
      if (local_eq) qval = geom%qfun(pbar) !leq:qfun ignores input pbar
      if (debug) write(6,*) "eikcoefs: call eikonal"
      call eikonal(theta, invR, trip, qval, bi, seik, dsdthet, dpsidrp)
      if (debug) write(6,*) "eikcoefs: done eikonal"
      if (writelots) write(fort11_unit,*) 'q= ', qval
      if (local_eq) trip = trip * dpsidrp

      if (bishop > 0) then !Following bits _could_ be calculated for all bishop
         call B_pol(invR, crpgrad, dpsidrp, bpol, nth) !<Should this use rpgrad for local?
         call periodic_copy(bpol, 0., nperiod)
         call B_mod(invR, bpol, bmag, bi, nth)
         call periodic_copy(bmag, 0., nperiod)
      end if !Note bmag not set for bishop == 0
    end subroutine initialise_geometry_and_get_bmag

    subroutine bishop_grads(rhoc, rmaj, bmag, bpol)
      real, intent(in) :: rhoc, rmaj
      real, dimension(-ntgrid:), intent(in) :: bmag, bpol
      logical, parameter :: save_bishop = .true.
      integer :: bishop_unit
      real :: a_b, b_b, c_b, dp, dp_new, di, di_new, s_hat, s_hat_new, pressure, tmp
      real, allocatable, dimension(:) :: a_bish, b_bish, c_bish, da_bish, db_bish, dc_bish, &
           zoftheta, dzdl, dsdl, ltheta, th_bish, invRpol
      allocate(a_bish(-ntgrid:ntgrid), b_bish(-ntgrid:ntgrid), c_bish(-ntgrid:ntgrid), &
           da_bish(-ntgrid:ntgrid), db_bish(-ntgrid:ntgrid), dc_bish(-ntgrid:ntgrid), &
           zoftheta(-ntgrid:ntgrid), dzdl(-ntgrid:ntgrid), dsdl(-ntgrid:ntgrid), &
           ltheta(-ntgrid:ntgrid), th_bish(-ntgrid:ntgrid), invRpol(-ntgrid:ntgrid))

      call th_bishop(crpgrad, th_bish, nth)
      call periodic_copy(th_bish, 0., nperiod)
      call loftheta(rgrid, outputs%theta, ltheta)
      call periodic_copy(ltheta, ltheta(nth) - ltheta(-nth), nperiod)
      call inv_Rpol(outputs%theta, th_bish, ltheta, invRpol, nth)
      call periodic_copy(invRpol, 0., nperiod)
      if (debug) write(6,*) "eikcoefs:  gradients"
      if (debug) write(6,*) 'drhodpsi', drhodpsi
      if (debug) write(6,*) 'bi', bi
      if (verb > 5) write(6,*) 'rmajor', rmajor(:)
      if (verb > 5) write(6,*) 'bpol', bpol(:)
      if (verb > 5) write(6,*) rpgrad(-nth:nth,1)
      if (verb > 5) write(6,*) rpgrad(-nth:nth,2)
      if (verb > 5) write(6,*) thgrad(-nth:nth,1)
      if (verb > 5) write(6,*) thgrad(-nth:nth,2)
      if (verb > 5) write(6,*) 1.0 / trip(-nth:nth)
      if (verb > 5) write(6,*) -1.0 / (bpol(-nth:nth) * invRpol(-nth:nth))
      if (debug) write(6,*) "eikcoefs: end gradients"

      da_bish = (1 / rmajor**2) + (bi / (bpol * rmajor**2))**2
      db_bish = bi / (bpol * rmajor)**2
      dc_bish = bi * (sin(th_bish) + rmajor * invRpol) / (-bpol * rmajor**4)

      da_bish = da_bish / trip
      db_bish = db_bish / trip
      dc_bish = dc_bish / trip

      call integrate(da_bish, outputs%theta, a_bish, ntgrid)
      call integrate(db_bish, outputs%theta, b_bish, ntgrid)
      call integrate(dc_bish, outputs%theta, c_bish, ntgrid)

      a_b = a_bish(nth) - a_bish(-nth)
      b_b = b_bish(nth) - b_bish(-nth)
      c_b = c_bish(nth) - c_bish(-nth)

      a_bish = - a_bish * bpol * rmajor
      b_bish = - b_bish * bpol * rmajor
      c_bish = - c_bish * bpol * rmajor

      ! find shat, pressure' and report them
      dp = geom%dpfun(pbar) ; di = geom%dbtori(pbar)

      if (abs(qval) > epsilon(0.)) then
         tmp = rhoc / (qval * twopi * drhodpsi)
         s_hat = tmp * (a_b * di + b_b * dp + c_b * 2) !< Calculate the shear from equib.
         if (local_eq) s_hat = surf%shat !Can't determine this by calculation in local_eq
         if (debug) write(6,*) "a_b (f'), b_b (p'), c_b terms =", a_b*di, b_b*dp, c_b*2
         if (debug) write(6,*) 'qval, rho, drhodpsi, s_hat =', qval, rhoc, drhodpsi, s_hat
      else
         tmp = 1. ; s_hat = 0.
      end if

      write(fort11_unit,*)
      write(fort11_unit,*) 'Quantity, equilibrium value, value used'

      ! Decide how we want to set s_hat_new in the default case
      s_hat_new = s_hat
      select case (bishop)
      case (2, 3, 4, 5, 8, 9) !1, 6, 7 keep the initial value
         s_hat_new = surf%shat
      end select

      ! Decide how we want to set dp_new
      dp_new = dp
      select case (bishop)
      case (2)
         dp_new = p_prime_input
      case (3)
         pressure = geom%pfun(pbar)
         dp_new = -invLp_input * pressure * drhodpsi
      case (4, 6, 9)
         dp_new = 0.5 * beta_prime_input * drhodpsi
      case (5)
         dp_new = -alpha_input * drhodpsi / (qval**2 * rmaj)
      case (7)
         dp_new = dp * dp_mult
      case (8)
         dp_new = 0.5 * beta_prime_input * drhodpsi * dp_mult
      end select

      !Report on updated dp
      select case(bishop)
      case (2)
         write(fort11_unit,*) 'p_prime = ',dp,', ',dp_new
         if (writelots) write(*,*) 'p_prime = ',dp,', ',dp_new
      case (3)
         write(fort11_unit,*) '1/L_p = ',-dp / (drhodpsi * pressure), ', ', &
              -dp_new / (drhodpsi *pressure)
         if (writelots) write(*,*) '1/L_p = ',-dp / (drhodpsi * pressure), ', ', &
              -dp_new / (drhodpsi * pressure)
      case (5)
         write(fort11_unit,*) 'alpha = ', -qval**2 * rmaj * dp / drhodpsi,', ', &
              -qval**2 * rmaj * dp_new / drhodpsi
         if (writelots) write(*,*) 'alpha = ', -qval**2 * rmaj * dp / drhodpsi,', ', &
              -qval**2 * rmaj * dp_new / drhodpsi
      case default
         write(fort11_unit,*) 'd beta/d rho = ',2 * dp / drhodpsi,', ',2 * dp_new / drhodpsi
         if (writelots) write(*,*) 'd beta/d rho = ',2 * dp / drhodpsi,', ',2 * dp_new / drhodpsi
      end select

      write(fort11_unit,*) 's_hat ',s_hat,', ',s_hat_new
      if (writelots) write(*,*) 's_hat ',s_hat,', ',s_hat_new
      write(fort11_unit,*)

      if ((.not. local_eq) .and. (bishop == 1 .or. bishop > 9)) then
         di_new = di
      else
         di_new = (s_hat_new / tmp - 2 * c_b - dp_new * b_b) / a_b
      end if

      ! Calculate grad(S)
      ! here the coordinates are rho and l
      call gradl(ltheta, seik, dSdl, twopi * qval, nth)
      call periodic_copy(dSdl, 0., nperiod)
      grads(:,1) = a_bish * di_new + b_bish * dp_new + c_bish * 2
      grads(:,2) = dSdl

      ! Calculate grad(Z)
      ! obtain Z(theta) and hence dZ/dl
      call Ztgrid(rgrid, outputs%theta, Zoftheta)
      call periodic_copy(Zoftheta, 0., nperiod)
      call gradl(ltheta, Zoftheta, dZdl, 0., nth)
      call periodic_copy(dZdl, 0., nperiod)

      ! gradztot(:,1) is grad Z . rhohat
      gradztot(:,1) = crpgrad(:,2) / outputs%grho
      ! gradztot(:,2) is grad Z . lhat
      gradztot(:,2) = dZdl
      ! gradztot(:,3) is grad Z . phihat
      gradztot(:,3) = 0.0

      shat = s_hat_new ! Output
      dqdrp = s_hat_new * qval * drhodrp / rhoc ! Output
      dbetadrho = 2 * dp_new * dpsidrho ! Output

      call bishop_gradB(bmag, bpol, invRpol, th_bish, ltheta, rmajor, dp_new, bi, gradbtot)

      if (save_bishop) then
         call open_output_file(bishop_unit, '.bishop')
         write(bishop_unit, fmt="('#di, di_new = ',F12.7,' ',F12.7)") di, di_new
         write(bishop_unit, fmt="('#dp, dp_new = ',F12.7,' ',F12.7)") dp, dp_new
         write(bishop_unit, fmt="('#bi =', F12.7)") bi
         write(bishop_unit, fmt="('#drdpsi =', F12.7)") drhodpsi
         write(bishop_unit, fmt="('#a_b (fp)=', F20.10, ' ,b_b (pp)=', F20.10,' ,c_b  =', F20.10,' ,s_hat =', F20.10)") a_b*di, b_b*dp, c_b, s_hat
         write(bishop_unit, fmt ="('#  theta, a_bish,  b_bish,  c_bish, rmajor, Z, Bpol, invRpol, th_bish')")
         do i = -ntgrid, ntgrid
            write(bishop_unit, '( 10(2X, F25.17) )') outputs%theta(i), a_bish(i), b_bish(i), &
                 c_bish(i), rmajor(i), Zoftheta(i),  bpol(i), invRpol(i), th_bish(i)
         end do
         call close_output_file(bishop_unit)
      end if
    end subroutine bishop_grads

    !> @Note local_eq _must_ be false if we end in this routine. We make use of this.
    subroutine bishop_large_aspect_ratio_approximation(rhoc, dr, ntheta)
      real, intent(in) :: rhoc, dr
      integer, intent(in) :: ntheta
      real :: delrp, rp1, rp2, qval1, qval2, bi1, bi2, dqdrho, pbar_1, pbar_2
      real, dimension(:), allocatable :: seik1, seik2, dsdrp
      integer :: itot, k
      !     rp derivative
      delrp = dr / drhodrp
      rp1 = rp - delrp ; rp2 = rp + delrp
      pbar_1 = pbarfun(rp1, outputs%theta(0)) ; pbar_2 = pbarfun(rp2, outputs%theta(0))
      ! Calculate modified btori -- note we use the fact that rgrid == rp in
      ! seikon to avoid having to form rgrid1/rgrid2 here.
      bi1 = geom%btori(pbar_1) ; bi2 = geom%btori(pbar_2)

      allocate(seik1(-ntgrid:ntgrid), seik2(-ntgrid:ntgrid), dsdrp(-ntgrid:ntgrid))
      call seikon(rp1, outputs%theta, invR, qval1, bi1, nperiod, seik1)
      call seikon(rp2, outputs%theta, invR, qval2, bi2, nperiod, seik2)

      dqdrho = (qval2 - qval1) / (2 * dr)
      shat = dqdrho * rhoc / qval ! Output
      dqdrp = dqdrho * drhodrp ! Output
      dbetadrho = 2 * geom%dpfun(pbar) / drhodrp ! Output

      if (writelots) then
         write (fort11_unit,*) 'dqdrho=',dqdrho,' qval1,2= ',qval1, qval2
         write (fort11_unit,*) 's_hat= ',shat
         write (fort11_unit,*) 'i,rgrid,rgrid1,rgrid2:'
         do i = -nth, nth
            write (fort11_unit,*) i, rgrid(i), rp1, rp2
         end do
      end if

      do i = -nth, nth
         dsdrp(i) = (seik2(i) - seik1(i)) / (2 * delrp)
      end do
      call periodic_copy(dsdrp, -twopi * dqdrp, nperiod)

      !     compute gradient of S
      ! this is really gradient of alpha
      ! note that coordinates are R and Z (not rho and l, as in bishop case)
      do k = -nperiod+1, nperiod-1
         do i = -nth, nth
            itot = i + k * ntheta
            grads(itot, :) = rpgrad(i, :) * dsdrp(itot) + thgrad(i, :) * dsdthet(i)  ! Output
         end do
      end do
      ! gradztot(:,1) is grad Z . Rhat ! Output
      gradztot(:,1) = 0.
      ! gradztot(:,2) is grad Z . Zhat ! Output
      gradztot(:,2) = 1.
      ! gradztot(:,3) is grad Z . phihat ! Output
      gradztot(:,3) = 0.

      call geom%gradient(rgrid, outputs%theta, gradbtot(:, 1:2), 'B', rp, nth, ntgrid, bishop /= 0)
      gradbtot(:, 3) = 0.0
    end subroutine bishop_large_aspect_ratio_approximation

    !> FIXME : Add documentation
    subroutine alloc_local_arrays(n)
      implicit none
      integer, intent(in) :: n

      allocate(seik   (-n:n), &
           dsdthet    (-n:n), &
           invR       (-n:n), &
           trip       (-n:n), &
           rgrid      (-n:n), &
           rmajor     (-n:n))

      allocate(thgrad (-n:n,2), &
           rpgrad     (-n:n,2), &
           crpgrad    (-n:n,2), &
           grads      (-n:n,2))

      allocate(gradrptot (-n:n,3), &
           gradstot   (-n:n,3), &
           gradztot   (-n:n,3), &
           gradthtot  (-n:n,3), & ! MAB
           gradbtot   (-n:n,3))
    end subroutine alloc_local_arrays
  end subroutine eikcoefs