geometry_miller.f90 Source File


Source Code

module geometry_miller

   use common_types, only: flux_surface_type

   implicit none

   public :: init_local_defaults
   public :: read_local_parameters
   public :: communicate_parameters_multibox
   public :: get_local_geo
   public :: finish_local_geo
   public :: local

   private
   
   ! These routines are always called on the first processor only
   logical :: debug = .false.

   integer :: nzed_local
   real :: rhoc, rmaj, shift
   real :: kappa, kapprim
   real :: tri, triprim
   real :: betaprim, betadbprim
   real :: qinp, shat, d2qdr2
   real :: rgeo
   real :: dpsipdrho, d2psidr2, dpsipdrho_psi0
   real :: psitor_lcfs
   real :: rhotor, drhotordrho, dIdrho, dI
   real :: rhoc0
   logical :: write_profile_variation, read_profile_variation
   logical :: load_psi0_variables

   integer :: nz, nz2pi

   real :: bi, dqdr, d2Idr2
   real, dimension(:), allocatable :: grho, bmag, grho_psi0, bmag_psi0, gradpar
   real, dimension(:), allocatable :: gradpararc, arc
   real, dimension(:), allocatable :: gds2, gds21, gds22
   real, dimension(:), allocatable :: gds23, gds24
   real, dimension(:), allocatable :: gbdrift0, gbdrift
   real, dimension(:), allocatable :: cvdrift0, cvdrift
   real, dimension(:), allocatable :: d2Rdth2, d2Zdth2, d2Rdrdth, d2Zdrdth
   real, dimension(:), allocatable :: gpsi, dBdrho, d2Bdrdth
   real, dimension(:), allocatable :: dgradpardrho, dgradparBdrho, dBdth, gradparb
   real, dimension(:), allocatable :: dcvdrift0drho, dgbdrift0drho, theta
   real, dimension(:), allocatable :: varthet, dvarthdr, gradrho_gradthet, cross, d2varthdr2
   real, dimension(:), allocatable :: gradthet2, gradalph_gradthet, gradrho_gradalph, gradalph2
   real, dimension(:), allocatable :: d2Bdr2, d2Rdr2, d2Zdr2, drz, drzdth
   real, dimension(:), allocatable :: d2Rdr2dth, d2Zdr2dth, d2gpsidr2, dcrossdr
   real, dimension(:), allocatable :: dcvdriftdrho, dgbdriftdrho
   real, dimension(:), allocatable :: dgds2dr, dgds21dr, dgds22dr
   real, dimension(:), allocatable :: dgr2dr, dgpsi2dr
   real, dimension(:), allocatable :: dgrgt, dgt2, dgagr, dgagt, dga2
   real, dimension(:, :), allocatable :: Rr, Zr

   real, dimension(:), allocatable :: jacrho, delthet, djacdrho, djacrdrho
   real, dimension(:), allocatable :: d2jacdr2, dRdrho, dZdrho, dRdth, dZdth

   real, dimension(:), allocatable :: d2R, d2Z

   type(flux_surface_type) :: local

   logical :: defaults_initialized = .false.

contains

   !============================================================================
   !=========================== INIT LOCAL DEFAULTS ============================
   !============================================================================ 
   subroutine init_local_defaults

      implicit none

      if (defaults_initialized) return
      defaults_initialized = .true.

      nzed_local = 128
      rhoc = 0.5
      rhoc0 = 0.5
      rmaj = 3.0
      rgeo = 3.0
      qinp = 1.4
      shat = 0.8
      shift = 0.0
      kappa = 0.0
      kapprim = 0.0
      tri = 0.0
      triprim = 0.0
      ! betaprim = -(4pi/Bref^2)*d(ptot)/drho
      betaprim = 0.0
      ! betadbprim = -(4pi/Bref^2)*d^2ptot/drho^2
      betadbprim = 0.0
      d2qdr2 = 0.0
      d2psidr2 = 0.0
      read_profile_variation = .false.
      write_profile_variation = .false.
      load_psi0_variables = .true.

      ! only needed for sfincs when not using
      ! geo info from file
      rhotor = rhoc
      psitor_lcfs = 1.0
      drhotordrho = 1.0

   end subroutine init_local_defaults

   !============================================================================
   !=========================== READ LOCAL PARAMETERS ==========================
   !============================================================================ 
   subroutine read_local_parameters(nzed, nzgrid, local_out)

      use file_utils, only: input_unit_exist
      use common_types, only: flux_surface_type

      implicit none

      type(flux_surface_type), intent(out) :: local_out
      integer, intent(in) :: nzed, nzgrid

      real :: dum
      integer :: in_file, np, j
      logical :: exist
      
      namelist /millergeo_parameters/ rhoc, rmaj, shift, qinp, shat, &
         kappa, kapprim, tri, triprim, rgeo, betaprim, &
         betadbprim, d2qdr2, d2psidr2, &
         nzed_local, read_profile_variation, write_profile_variation


      if (debug) write (*, *) 'geometry_miller::read_local_parameters'
      call init_local_defaults

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

      local%rhoc = rhoc
      local%rmaj = rmaj
      local%rgeo = rgeo
      local%shift = shift
      local%kappa = kappa
      local%kapprim = kapprim
      local%qinp = qinp
      local%shat = shat
      local%tri = tri
      local%triprim = triprim
      local%betaprim = betaprim
      local%betadbprim = betadbprim
      local%d2qdr2 = d2qdr2
      local%d2psidr2 = d2psidr2
      local%zed0_fac = 1.0

      ! following two variables are not inputs
      local%dr = 1.e-3 * (rhoc / rmaj)
      local%rhotor = rhotor
      local%psitor_lcfs = psitor_lcfs
      local%drhotordrho = drhotordrho
      local%dpsitordrho = 0.0
      local%d2psitordrho2 = 0.0

      ! the next three variablaes are for multibox simulations
      ! with radial variation
      local%rhoc_psi0 = rhoc
      local%qinp_psi0 = qinp
      local%shat_psi0 = shat

      ! first get nperiod corresponding to input number of grid points
      nz2pi = nzed / 2
      np = (nzgrid - nz2pi) / nzed + 1

      ! now switch to using (possible higher resolution) local grid
      nz2pi = nzed_local / 2
      ! this is the equivalent of nzgrid on the local grid
      nz = nz2pi + nzed_local * (np - 1)

      ! initialize to zero
      ! will be overwritten if reading in from file
      ! only relevant for profile variation tests
      ! these needs to be deallocated somewhere
      allocate (d2R(-nz:nz))
      allocate (d2Z(-nz:nz))
      allocate (bmag_psi0(-nz:nz))
      allocate (grho_psi0(-nz:nz))
      d2R = 0.; d2Z = 0.; dI = 0.

      if (read_profile_variation) then
         open (1002, file='RZ.in', status='old')
         read (1002, '(12e13.5)') rhoc0, dI, qinp, shat, d2qdr2, kappa, kapprim, tri, triprim, &
            betaprim, betadbprim, dpsipdrho_psi0
         do j = -nz, nz
            read (1002, '(5e13.5)') dum, d2R(j), d2Z(j), bmag_psi0(j), grho_psi0(j)
         end do
         close (1002)
         local%qinp = qinp + shat * qinp / rhoc0 * (local%rhoc - rhoc0) &
                      + 0.5 * (local%rhoc - rhoc0)**2 * d2qdr2
         local%shat = (local%rhoc / local%qinp) &
                      * (shat * qinp / rhoc0 + (local%rhoc - rhoc0) * d2qdr2)
         local%kappa = kappa + kapprim * (local%rhoc - rhoc0)
         local%tri = tri + triprim * (local%rhoc - rhoc0)
         local%betaprim = betaprim + betadbprim * (local%rhoc - rhoc0)

         local%rhoc_psi0 = rhoc0
         local%qinp_psi0 = qinp
         local%shat_psi0 = shat

         load_psi0_variables = .false.
      end if

      local_out = local

   end subroutine read_local_parameters

   !============================================================================
   !=========================== COMMUNICATE PARAMETERS =========================
   !============================================================================ 
   subroutine communicate_parameters_multibox(surf, drl, drr)
      use mp, only: job, scope, mp_abort, &
                    crossdomprocs, subprocs, &
                    send, receive
      use job_manage, only: njobs
      use common_types, only: flux_surface_type

      implicit none

      real, optional, intent(in) :: drl, drr
      type(flux_surface_type), intent(inout) :: surf

      real :: lrhoc, lqinp, lshat, lkappa, ltri, lbetaprim
      real :: rrhoc, rqinp, rshat, rkappa, rtri, rbetaprim
      real :: dqdr
      real :: rhoc_psi0, qinp_psi0, shat_psi0

      !FLAG DSO -  I think d2psidrho2 needs to be communicated, but
      !            I'm unsure what quantity needs to be updated

      if (debug) write (*, *) 'geometry_miller::communicate_parameters_multibox'
      if (job == 1) then
         dqdr = local%shat * local%qinp / local%rhoc

         lrhoc = local%rhoc + drl
         lqinp = local%qinp + drl * dqdr + 0.5 * drl**2 * local%d2qdr2
         lshat = (lrhoc / lqinp) * (dqdr + drl * local%d2qdr2)
         lkappa = kappa + drl * kapprim
         ltri = tri + drl * triprim
         lbetaprim = betaprim + drl * betadbprim

         rrhoc = local%rhoc + drr
         rqinp = local%qinp + drr * dqdr + 0.5 * drr**2 * local%d2qdr2
         rshat = (rrhoc / rqinp) * (dqdr + drr * local%d2qdr2)
         rkappa = kappa + drr * kapprim
         rtri = tri + drr * triprim
         rbetaprim = betaprim + drr * betadbprim
      end if

      call scope(crossdomprocs)

      if (job == 1) then
         call send(lrhoc, 0, 120)
         call send(lqinp, 0, 121)
         call send(lshat, 0, 122)
         call send(lkappa, 0, 123)
         call send(ltri, 0, 124)
         call send(lbetaprim, 0, 125)
         call send(local%rhoc, 0, 126)
         call send(d2R, 0, 127)
         call send(d2Z, 0, 128)
         call send(dIdrho, 0, 129)
         call send(rhoc, 0, 130)
         call send(qinp, 0, 131)
         call send(shat, 0, 132)
         call send(dpsipdrho, 0, 133)
         call send(bmag, 0, 134)
         call send(grho, 0, 135)

         call send(rrhoc, njobs - 1, 220)
         call send(rqinp, njobs - 1, 221)
         call send(rshat, njobs - 1, 222)
         call send(rkappa, njobs - 1, 223)
         call send(rtri, njobs - 1, 224)
         call send(rbetaprim, njobs - 1, 225)
         call send(local%rhoc, njobs - 1, 226)
         call send(d2R, njobs - 1, 227)
         call send(d2Z, njobs - 1, 228)
         call send(dIdrho, njobs - 1, 229)
         call send(rhoc, njobs - 1, 230)
         call send(qinp, njobs - 1, 231)
         call send(shat, njobs - 1, 232)
         call send(dpsipdrho, njobs - 1, 233)
         call send(bmag, njobs - 1, 234)
         call send(grho, njobs - 1, 235)
         rhoc_psi0 = rhoc
         qinp_psi0 = qinp
         shat_psi0 = shat
         local%rhoc_psi0 = rhoc_psi0
         local%qinp_psi0 = qinp_psi0
         local%shat_psi0 = shat_psi0
      elseif (job == 0) then
         call receive(rhoc, 1, 120)
         call receive(qinp, 1, 121)
         call receive(shat, 1, 122)
         call receive(kappa, 1, 123)
         call receive(tri, 1, 124)
         call receive(betaprim, 1, 125)
         call receive(rhoc0, 1, 126)
         call receive(d2R, 1, 127)
         call receive(d2Z, 1, 128)
         call receive(dI, 1, 129)
         call receive(rhoc_psi0, 1, 130)
         call receive(qinp_psi0, 1, 131)
         call receive(shat_psi0, 1, 132)
         call receive(dpsipdrho_psi0, 1, 133)
         call receive(bmag_psi0, 1, 134)
         call receive(grho_psi0, 1, 135)
         local%rhoc = rhoc
         local%qinp = qinp
         local%shat = shat
         local%kappa = kappa
         local%tri = tri
         local%betaprim = betaprim
         local%rhoc_psi0 = rhoc_psi0
         local%qinp_psi0 = qinp_psi0
         local%shat_psi0 = shat_psi0

         load_psi0_variables = .false.
      elseif (job == njobs - 1) then
         call receive(rhoc, 1, 220)
         call receive(qinp, 1, 221)
         call receive(shat, 1, 222)
         call receive(kappa, 1, 223)
         call receive(tri, 1, 224)
         call receive(betaprim, 1, 225)
         call receive(rhoc0, 1, 226)
         call receive(d2R, 1, 227)
         call receive(d2Z, 1, 228)
         call receive(dI, 1, 229)
         call receive(rhoc_psi0, 1, 230)
         call receive(qinp_psi0, 1, 231)
         call receive(shat_psi0, 1, 232)
         call receive(dpsipdrho_psi0, 1, 233)
         call receive(bmag_psi0, 1, 234)
         call receive(grho_psi0, 1, 235)
         local%rhoc = rhoc
         local%qinp = qinp
         local%shat = shat
         local%kappa = kappa
         local%tri = tri
         local%betaprim = betaprim
         local%rhoc_psi0 = rhoc_psi0
         local%qinp_psi0 = qinp_psi0
         local%shat_psi0 = shat_psi0

         load_psi0_variables = .false.
      end if

      surf%rhoc = local%rhoc
      surf%qinp = local%qinp
      surf%shat = local%shat
      surf%kappa = local%kappa
      surf%tri = local%tri
      surf%betaprim = local%betaprim
      surf%rhoc_psi0 = rhoc_psi0
      surf%qinp_psi0 = qinp_psi0
      surf%shat_psi0 = shat_psi0

      call scope(subprocs)

   end subroutine communicate_parameters_multibox

   !============================================================================
   !========================= CALCULATE MILLER GEOMETRY ========================
   !============================================================================ 
	! gradpar_out is called b_dot_grad_zeta in geometry
   subroutine get_local_geo(nzed, nzgrid, zed_in, zed_equal_arc, &
                            dpsipdrho_out, dpsipdrho_psi0_out, dIdrho_out, grho_out, &
                            bmag_out, bmag_psi0_out, &
                            gds2_out, gds21_out, gds22_out, gds23_out, gds24_out, gradpar_out, &
                            gbdrift0_out, gbdrift_out, cvdrift0_out, cvdrift_out, &
                            dBdrho_out, d2Bdrdth_out, dgradpardrho_out, &
                            btor_out, rmajor_out, &
                            dcvdrift0drho_out, dcvdriftdrho_out, &
                            dgbdrift0drho_out, dgbdriftdrho_out, &
                            dgds2dr_out, dgds21dr_out, &
                            dgds22dr_out, djacdrho_out)

      use constants, only: pi
      use splines, only: geo_spline
      use file_utils, only: run_name

      implicit none

      integer, intent(in) :: nzed, nzgrid
      real, dimension(-nzgrid:), intent(in) :: zed_in
      logical, intent(in) :: zed_equal_arc
      real, intent(out) :: dpsipdrho_out, dpsipdrho_psi0_out, dIdrho_out
      real, dimension(-nzgrid:), intent(out) :: grho_out, &
                                                bmag_out, bmag_psi0_out, &
                                                gds2_out, gds21_out, gds22_out, gds23_out, gds24_out, &
                                                gradpar_out, gbdrift0_out, &
                                                gbdrift_out, cvdrift0_out, cvdrift_out, &
                                                dBdrho_out, d2Bdrdth_out, dgradpardrho_out, &
                                                btor_out, rmajor_out, &
                                                dcvdrift0drho_out, dcvdriftdrho_out, &
                                                dgbdrift0drho_out, dgbdriftdrho_out, &
                                                dgds2dr_out, dgds21dr_out, &
                                                dgds22dr_out, &
                                                djacdrho_out

      integer :: nr, np
      integer :: i, j
      real :: rmin, dum
      real, dimension(3) :: dr
      real, allocatable, dimension(:) :: zed_arc
      character(len=512) :: filename
      integer :: n
      
      ! Track code
      if (debug) write (*, *) 'geometry_miller::get_local_geo'
      
      ! number of grid points used for radial derivatives
      nr = 3

      ! first get nperiod corresponding to input number of grid points
      nz2pi = nzed / 2
      np = (nzgrid - nz2pi) / nzed + 1

      ! now switch to using (possible higher resolution) local grid
      nz2pi = nzed_local / 2
      ! this is the equivalent of nzgrid on the local grid
      nz = nz2pi + nzed_local * (np - 1)

      ! Allocate arrays
      if (debug) write (*, *) 'geometry_miller::get_local_geo'
      call allocate_arrays(nr, nz)

      dqdr = local%shat * local%qinp / local%rhoc

      dr(1) = -local%dr
      dr(2) = 0.
      dr(3) = local%dr

      do j = -nz, nz
         theta(j) = j * (2 * np - 1) * pi / real(nz)
         do i = 1, 3
            rmin = local%rhoc + dr(i)
            Rr(i, j) = Rpos(rmin, theta(j), j)
            Zr(i, j) = Zpos(rmin, theta(j), j)
         end do
      end do

      if (.not. allocated(delthet)) allocate (delthet(-nz:nz - 1))
      ! get delta theta as a function of theta
      delthet = theta(-nz + 1:) - theta(:nz - 1)
      
      
      if (debug) write (*, *) 'geometry_miller::radial_derivatives'

      ! get dR/drho and dZ/drho
      call get_drho(Rr, dRdrho)
      call get_drho(Zr, dZdrho)

      ! get dR/dtheta and dZ/dtheta
      call get_dthet(Rr(2, :), dRdth)
      call get_dthet(Zr(2, :), dZdth)

      ! get second derivatives of R and Z with respect to theta
      call get_d2dthet2(Rr(2, :), d2Rdth2)
      call get_d2dthet2(Zr(2, :), d2Zdth2)
      ! get mixed theta and rho derivatives of R and Z
      call get_dthet(dRdrho, d2Rdrdth)
      call get_dthet(dZdrho, d2Zdrdth)

      ! get the Jacobian of the transformation from (rho,theta,zeta) to (R,Z,zeta)
      ! this is what I call jacr or jacrho in following comments
      ! as opposed to jacobian, which is for tranformation from (psi,theta,zeta) to (R,Z,zeta)
      call get_jacrho

      ! theta_integrate returns integral from 0 -> 2*pi
      ! note that dpsipdrho here is an intermediary
      ! that requires manipulation to get final dpsipdrho
      call theta_integrate(jacrho(-nz2pi:nz2pi) / Rr(2, -nz2pi:nz2pi)**2, dpsipdrho)
      dpsipdrho = dpsipdrho / (2.*pi)

      ! get dpsinorm/drho = (I/2*pi*q)*int_0^{2*pi} dthet jacrho/R**2

      ! if using input.profiles, we are given
      ! dpsitordrho and must use it to compute rgeo
      if (abs(local%dpsitordrho) > epsilon(0.)) then
         local%rgeo = local%dpsitordrho / dpsipdrho
         dpsipdrho = local%dpsitordrho / local%qinp
         local%d2psidr2 = (local%d2psitordrho2 - local%dpsitordrho * local%shat / local%rhoc) &
                          / local%qinp
         ! I=Btor*R is a flux function
         ! bi = I/(Btor(psi,theta of Rgeo)*a) = Rgeo/a
         bi = local%rgeo
      else
         ! otherwise, we are given rgeo
         ! and must use it to compute dpsipdrho

         ! I=Btor*R is a flux function
         ! bi = I/(Btor(psi,theta of Rgeo)*a) = Rgeo/a
         bi = local%rgeo + dI * (rhoc - rhoc0)
         dpsipdrho = dpsipdrho * bi / local%qinp
      end if

!    ! get dpsinorm/drho
!    call get_dpsipdrho (dpsipdrho)

      ! get |grad rho| and |grad psi|
      call get_gradrho(dpsipdrho, grho)

      ! quantity needed in calculation of dI/drho and djacrho/drho
      drz = (dRdrho * dRdth + dZdrho * dZdth) / jacrho
      call get_dthet(drz, drzdth)

      ! get dI/drho
      call get_dIdrho(dpsipdrho, grho, dIdrho)
      dIdrho_out = dIdrho

      ! get djacobian/drho*dpsi/drho and djacr/drho
      call get_djacdrho(dpsipdrho, dIdrho, grho)

      ! get d2R/drho2 and d2Z/drho2
      call get_d2RZdr2

      d2R = d2Rdr2
      d2Z = d2Zdr2

      ! get theta derivative of d2R/drho2 and d2Z/drho2
      call get_dthet(d2Rdr2, d2Rdr2dth)
      call get_dthet(d2Zdr2, d2Zdr2dth)

      ! calculate the magnitude of B (normalized by B(psi,theta corresponding to Rgeo))
      ! B/B0 = sqrt(I**2 + |grad psi|**2)/R
      bmag = sqrt(bi**2 + gpsi**2) / Rr(2, :)

      ! the next line is for multibox runs
      if (load_psi0_variables) then
         dpsipdrho_psi0 = dpsipdrho
         bmag_psi0 = bmag
         grho_psi0 = grho
      end if

      if (write_profile_variation) then
         open (1002, file='RZ.out', status='unknown')
         write (1002, '(12e13.5)') local%rhoc, dIdrho, local%qinp, local%shat, local%d2qdr2, &
            local%kappa, local%kapprim, &
            local%tri, local%triprim, &
            local%betaprim, local%betadbprim, dpsipdrho
         do j = -nz, nz
            write (1002, '(5e13.5)') theta(j), d2Rdr2(j), d2Zdr2(j), bmag(j), grho(j)
         end do
         close (1002)
      end if

      ! get dB/dtheta
      call get_dthet(bmag, dbdth)

      ! calculate b . grad theta
      gradpar = dpsipdrho / (bmag * jacrho)
      ! b . grad B
      gradparb = gradpar * dBdth

      ! get d|grad rho|^2/drho and d|grad psi|^2/drho
      call get_dgr2dr(dpsipdrho, grho)

      ! get dB/drho and d2B/drho2
      call get_dBdrho(bmag, dIdrho)

      ! d (b . grad theta) / drho
      dgradpardrho = -gradpar * (dBdrho / bmag + djacdrho / jacrho)

      ! get d/dtheta (dB/drho)
      call get_dthet(dBdrho, d2Bdrdth)

      ! d(b . grad B)/drho
      dgradparBdrho = dgradpardrho * dBdth + gradpar * d2Bdrdth

      ! obtain varthet = (I/(q*(dpsi/dr)) * int_0^theta dtheta' jacrho/R^2
      call get_varthet(dpsipdrho)

      ! obtain dvarthet/drho
      call get_dvarthdr(dpsipdrho, dIdrho)

      ! get |grad theta|^2, grad r . grad theta, grad alpha . grad theta, etc.
      call get_graddotgrad(dpsipdrho, grho)
      
      if (debug) write (*, *) 'geometry_miller::get_gds'
      call get_gds(gds2, gds21, gds22, gds23, gds24)

      ! this is (grad alpha x B) . grad theta
      cross = dpsipdrho * (gradrho_gradalph * gradalph_gradthet - gradalph2 * gradrho_gradthet)

      ! note that the definitions of gbdrift, gbdrift0, dgbdriftdr and dgbdrift0dr
      ! are such that it gets multiplied by vperp2, not mu.  This is in contrast to
      ! Michael's GS3 notes

      ! this is bhat/B x (grad B/B) . grad alpha * 2 * dpsiN/drho
      gbdrift = 2.0 * (-dBdrho + cross * dBdth * dpsipdrho / bmag**2) / bmag
      ! this is bhat/B x (bhat . grad bhat) . grad alpha * 2 * dpsiN/drho
      ! this is assuming betaprim = 4*pi*ptot/B0^2 * (-d ln ptot / drho)
      cvdrift = (gbdrift + 2.0 * local%betaprim / bmag**2)

      ! this is 2 *(bhat/B x grad B / B) . (grad q) * dpsiN/drho / (bhat . grad B)
      ! same as usual GS2 definition once bhat . grad B is added in below
      cvdrift0 = -2.*bi * dqdr / bmag**2

      ! this is 2*dpsiN/drho times the rho derivative (bhat/B x grad B / B) . (grad q)
      dcvdrift0drho = cvdrift0 * (dgradparbdrho + gradparb * (dIdrho / bi - 2.*dBdrho / bmag - local%d2psidr2 / dpsipdrho)) &
                      - 2.*bi * gradparb * local%d2qdr2 / bmag**2
      ! this is 2*dpsiN/drho/B times the rho derivative of (bhat x gradB/B) . (grad q)
      ! note that there's an extra factor of 1/B that's not expanded due to v_perp -> mu
      dgbdrift0drho = cvdrift0 * (dgradparbdrho + gradparb * (dIdrho / bi - dBdrho / bmag - local%d2psidr2 / dpsipdrho)) &
                      - 2.*bi * gradparb * local%d2qdr2 / bmag**2

      cvdrift0 = cvdrift0 * gradparb
      ! this is 2 * dpsiN/drho * (bhat/B x gradB/B) . (grad q)
      gbdrift0 = cvdrift0

      ! get d^2I/drho^2 and d^2 Jac / dr^2
      if (debug) write (*, *) 'geometry_miller::get_d2Idr2_d2jacdr2'
      call get_d2Idr2_d2jacdr2(grho, dIdrho)

      ! get d^2varhteta/drho^2
      call get_d2varthdr2(dpsipdrho, dIdrho)

      ! get d2B/drho^2
      call get_d2Bdr2(bmag, dIdrho)

      ! get d/dr [(grad alpha x B) . grad theta]
      call get_dcrossdr(dpsipdrho, dIdrho, grho)

      ! dgbdriftdrho is d/drho [(bhat/B x (grad B) . grad alpha) * 2 * dpsiN/drho] / B
      ! note that there's an extra factor of 1/B that's not expanded due to v_perp -> mu
      dgbdriftdrho = 2.0 * (local%d2psidr2 * dBdrho / dpsipdrho - d2Bdr2 &
                            + dpsipdrho * (dcrossdr * dBdth + cross * (d2Bdrdth - 2.*dBdth * dBdrho / bmag)) / bmag**2) / bmag
      ! dcvdriftdrho is d/drho (bhat/B x [bhat . grad bhat] . grad alpha) * 2 * dpsiN/drho
      dcvdriftdrho = dgbdriftdrho - gbdrift * dBdrho / bmag &
                     + 2.0 * local%betadbprim / bmag**2 - 4.0 * local%betaprim * dBdrho / bmag**3 &
                     - 2.0 * local%betaprim * local%d2psidr2 / dpsipdrho

      !the next two sets of lines are corrections needed for the side boxes in a multibox simulation
      !gbdrift  = gbdrift *(dpsipdrho_psi0/dpsipdrho)*(bmag/bmag_psi0)
      !gbdrift0 = gbdrift0*(dpsipdrho_psi0/dpsipdrho)*(bmag/bmag_psi0)
      gbdrift = gbdrift * (dpsipdrho_psi0 / dpsipdrho)
      gbdrift0 = gbdrift0 * (dpsipdrho_psi0 / dpsipdrho)
      cvdrift = cvdrift * (dpsipdrho_psi0 / dpsipdrho)
      cvdrift0 = cvdrift0 * (dpsipdrho_psi0 / dpsipdrho)

      !dgbdriftdrho  = dgbdriftdrho *(dpsipdrho_psi0/dpsipdrho)*(bmag/bmag_psi0)
      !dgbdrift0drho = dgbdrift0drho*(dpsipdrho_psi0/dpsipdrho)*(bmag/bmag_psi0)
      dgbdriftdrho = dgbdriftdrho * (dpsipdrho_psi0 / dpsipdrho)
      dgbdrift0drho = dgbdrift0drho * (dpsipdrho_psi0 / dpsipdrho)
      dcvdriftdrho = dcvdriftdrho * (dpsipdrho_psi0 / dpsipdrho)
      dcvdrift0drho = dcvdrift0drho * (dpsipdrho_psi0 / dpsipdrho)

      ! interpolate here
      if (debug) write (*, *) 'geometry_miller::zed_equal_arc'
      if (zed_equal_arc) then
         call theta_integrate(1./gradpar, dum)
         gradpararc = (theta(nz) - theta(-nz)) / ((2 * np - 1) * dum)
         call theta_integrate_indef(gradpararc / gradpar, arc)

         allocate (zed_arc(-nzgrid:nzgrid))

         call geo_spline(arc, theta, zed_in, zed_arc)

         call geo_spline(theta, grho_psi0, zed_arc, grho_out) !grho is used to normalize fluxes
         call geo_spline(theta, bmag, zed_arc, bmag_out)
         call geo_spline(theta, bmag_psi0, zed_arc, bmag_psi0_out)
         call geo_spline(theta, gds2, zed_arc, gds2_out)
         call geo_spline(theta, gds21, zed_arc, gds21_out)
         call geo_spline(theta, gds22, zed_arc, gds22_out)
         call geo_spline(theta, gds21, zed_arc, gds23_out)
         call geo_spline(theta, gds21, zed_arc, gds24_out)
         call geo_spline(theta, gradpararc, zed_arc, gradpar_out)
         call geo_spline(theta, gbdrift, zed_arc, gbdrift_out)
         call geo_spline(theta, gbdrift0, zed_arc, gbdrift0_out)
         call geo_spline(theta, cvdrift, zed_arc, cvdrift_out)
         call geo_spline(theta, cvdrift0, zed_arc, cvdrift0_out)
         call geo_spline(theta, dBdrho, zed_arc, dBdrho_out)
         call geo_spline(theta, d2Bdrdth, zed_arc, d2Bdrdth_out)
         call geo_spline(theta, dgradpardrho, zed_arc, dgradpardrho_out)
         call geo_spline(theta, Rr(2, :), zed_arc, rmajor_out)
         call geo_spline(theta, dcvdriftdrho, zed_arc, dcvdriftdrho_out)
         call geo_spline(theta, dgbdriftdrho, zed_arc, dgbdriftdrho_out)
         call geo_spline(theta, dcvdrift0drho, zed_arc, dcvdrift0drho_out)
         call geo_spline(theta, dgbdrift0drho, zed_arc, dgbdrift0drho_out)
         call geo_spline(theta, dgds2dr, zed_arc, dgds2dr_out)
         call geo_spline(theta, dgds21dr, zed_arc, dgds21dr_out)
         call geo_spline(theta, dgds22dr, zed_arc, dgds22dr_out)
         call geo_spline(theta, djacdrho / dpsipdrho, zed_arc, djacdrho_out)

         deallocate (zed_arc)
      else
         call geo_spline(theta, grho_psi0, zed_in, grho_out) !grho is used to normalize fluxes
         call geo_spline(theta, bmag, zed_in, bmag_out)
         call geo_spline(theta, bmag_psi0, zed_in, bmag_psi0_out)
         call geo_spline(theta, gds2, zed_in, gds2_out)
         call geo_spline(theta, gds21, zed_in, gds21_out)
         call geo_spline(theta, gds22, zed_in, gds22_out)
         call geo_spline(theta, gds21, zed_in, gds23_out)
         call geo_spline(theta, gds21, zed_in, gds24_out)
         call geo_spline(theta, gradpar, zed_in, gradpar_out)
         call geo_spline(theta, gbdrift, zed_in, gbdrift_out)
         call geo_spline(theta, gbdrift0, zed_in, gbdrift0_out)
         call geo_spline(theta, cvdrift, zed_in, cvdrift_out)
         call geo_spline(theta, cvdrift0, zed_in, cvdrift0_out)
         call geo_spline(theta, dBdrho, zed_in, dBdrho_out)
         call geo_spline(theta, d2Bdrdth, zed_in, d2Bdrdth_out)
         call geo_spline(theta, dgradpardrho, zed_in, dgradpardrho_out)
         call geo_spline(theta, Rr(2, :), zed_in, rmajor_out)
         call geo_spline(theta, dcvdriftdrho, zed_in, dcvdriftdrho_out)
         call geo_spline(theta, dgbdriftdrho, zed_in, dgbdriftdrho_out)
         call geo_spline(theta, dcvdrift0drho, zed_in, dcvdrift0drho_out)
         call geo_spline(theta, dgbdrift0drho, zed_in, dgbdrift0drho_out)
         call geo_spline(theta, dgds2dr, zed_in, dgds2dr_out)
         call geo_spline(theta, dgds21dr, zed_in, dgds21dr_out)
         call geo_spline(theta, dgds22dr, zed_in, dgds22dr_out)
         call geo_spline(theta, djacdrho / dpsipdrho, zed_in, djacdrho_out)
      end if

      ! get the toroidal component of the magnetic field
      ! btor = B_toroidal/Bref = I/R Bref = rgeo * a/R
      btor_out = bi / rmajor_out

      dpsipdrho_out = dpsipdrho
      dpsipdrho_psi0_out = dpsipdrho_psi0

      n = 12
      do i = -nz, nz
         theta(i) = round(theta(i), n)
         Rr(2, i) = round(Rr(2, i), n)
         dRdrho(i) = round(dRdrho(i), n)
         d2Rdr2(i) = round(d2Rdr2(i), n) 
         dRdth(i) = round(dRdth(i), n) 
         d2Rdrdth(i) = round(d2Rdrdth(i), n) 
         dZdrho(i) = round(dZdrho(i), n)
         d2Zdr2(i) = round(d2Zdr2(i), n)
         dZdth(i) = round(dZdth(i), n)
         d2Zdrdth(i) = round(d2Zdrdth(i), n) 
         bmag(i) = round(bmag(i), n) 
         dBdrho(i) = round(dBdrho(i), n)
         d2Bdr2(i) = round(d2Bdr2(i), n)
         dBdth(i) = round(dBdth(i), n) 
         d2Bdrdth(i) = round(d2Bdrdth(i), n)
         varthet(i) = round(varthet(i), n) 
         dvarthdr(i) = round(dvarthdr(i), n)
         d2varthdr2(i) = round(d2varthdr2(i), n)
         jacrho(i) = round(jacrho(i), n)
         djacrdrho(i) = round(djacrdrho(i), n) 
         djacdrho(i) = round(djacdrho(i), n) 
         d2jacdr2(i) = round(d2jacdr2(i), n)
         grho(i) = round(grho(i), n)
         dgr2dr(i) = round(dgr2dr(i), n)
         gradthet2(i) = round(gradthet2(i), n)
         dgt2(i) = round(dgt2(i), n)
         gradrho_gradthet(i) = round(gradrho_gradthet(i), n)
         dgrgt(i) = round(dgrgt(i), n)
         gradalph_gradthet(i) = round(gradalph_gradthet(i), n)
         dgagt(i) = round(dgagt(i), n) 
         gradrho_gradalph(i) = round(gradrho_gradalph(i), n) 
         dgagr(i) = round(dgagr(i), n) 
         gradalph2(i) = round(gradalph2(i), n) 
         dga2(i) = round(dga2(i), n)
         cross(i) = round(cross(i), n) 
         dcrossdr(i) = round(dcrossdr(i), n) 
         gbdrift0(i) = round(gbdrift0(i), n) 
         dgbdrift0drho(i) = round(dgbdrift0drho(i), n) 
         cvdrift0(i) = round(cvdrift0(i), n) 
         dcvdrift0drho(i) = round(dcvdrift0drho(i), n)
         gbdrift(i) = round(gbdrift(i), n)
         dgbdriftdrho(i) = round(dgbdriftdrho(i), n)
         cvdrift(i) = round(cvdrift(i), n)
         dcvdriftdrho(i) = round(dcvdriftdrho(i), n)
         drzdth(i) = round(drzdth(i), n)
         gradpar(i) = round(gradpar(i), n)
         dgradpardrho(i) = round(dgradpardrho(i), n)
         gradparB(i) = round(gradparB(i), n)
         dgradparBdrho(i) = round(dgradparBdrho(i), n)
         gds2(i) = round(gds2(i), n)
         dgds2dr(i) = round(dgds2dr(i), n)
         gds21(i) = round(gds21(i), n)
         dgds21dr(i) = round(dgds21dr(i), n)
         gds22(i) = round(gds22(i), n)
         dgds22dr(i) = round(dgds22dr(i), n)
         gds23(i) = round(gds23(i), n)
         gds24(i) = round(gds24(i), n) 
         Zr(2, i) = round(Zr(2,i), n)
      end do

      
      if (debug) write (*, *) 'geometry_miller::write_geometry_miller_txt_files'
      filename = "geometry_miller."//trim(run_name)//".input"
      open (1002, file=trim(filename), status='unknown')
      write (1002, '(5a16)') '#1.rhoc', '2.rmaj', '3.rgeo', '4.shift', '5.qinp'
      write (1002, '(5e16.8)') local%rhoc, local%rmaj, local%rgeo, local%shift, local%qinp
      write (1002, *)
      write (1002, '(5a16)') '#6.shat', '7.kappa', '8.kapprim', '9.tri', '10.triprim'
      write (1002, '(5e16.8)') local%shat, local%kappa, local%kapprim, local%tri, local%triprim
      write (1002, *)
      write (1002, '(5a16)') '11.betaprim', '12.dpsitordrho', '13.rhotor', &
         '14.drhotordrho', '15.d2qdr2'
      write (1002, '(5e16.8)') local%betaprim, local%dpsitordrho, local%rhotor, &
         local%drhotordrho, local%d2qdr2
      write (1002, *)
      write (1002, '(3a16)') '16.d2psidr2', '17.betadbprim', '18.psitor_lcfs'
      write (1002, '(3e16.8)') local%d2psidr2, local%betadbprim, local%psitor_lcfs
      close (1002)
      filename = "geometry_miller."//trim(run_name)//".output"
      open (1001, file=trim(filename), status='unknown')
      write (1001, '(a9,e18.9,a11,e18.9,a11,e18.9)') '#dI/dr: ', dIdrho, 'd2I/dr2: ', d2Idr2, 'dpsi/dr: ', dpsipdrho
      write (1001, '(58a15)') '#1.theta', '2.R', '3.dR/dr', '4.d2Rdr2', '5.dR/dth', &
         '6.d2Rdrdth', '7.dZ/dr', '8.d2Zdr2', '9.dZ/dth', '10.d2Zdrdth', &
         '11.bmag', '12.dBdr', '13.d2Bdr2', '14.dB/dth', '15.d2Bdrdth', &
         '16.varthet', '17.dvarthdr', '18.d2varthdr2', '19.jacr', '20.djacrdr', &
         '21.djacdrho', '22.d2jacdr2', '23.grho2', '24.dgr2dr', '25.gthet2', &
         '26.dgt2', '27.grgthet', '28.dgrgt', '29.galphgth', '30.dgagt', &
         '31.grgalph', '32.dgagr', '33.galph2', '34.dga2', '35.cross', &
         '36.dcrossdr', '37.gbdrift0', '38.dgbdrift0', '39.cvdrift0', '40.dcvdrift0', &
         '41.gbdrift', '42.dgbdrift', '43.cvdrift', '44.dcvdrift', '45.drzdth', &
         '46.gradpar', '47.dgpardr', '48.gradparB', '49.dgparBdr', '50.gds2', &
         '51.dgds2dr', '52.gds21', '53.dgds21dr', '54.gds22', '55.dgds22dr', &
         '56.gds23', '57.gds24', '58.Zr'

      if (debug) write (*, *) 'geometry_miller::write_geometry_miller_txt_files::start_loop'
      if (debug) then 
         i = nz
         write(*,*) 'a', theta(i), Rr(2, i), dRdrho(i), d2Rdr2(i), dRdth(i)
         write(*,*) 'b', d2Rdrdth(i), dZdrho(i), d2Zdr2(i), dZdth(i), d2Zdrdth(i)
         write(*,*) 'c', bmag(i), dBdrho(i), d2Bdr2(i), dBdth(i), d2Bdrdth(i)
         write(*,*) 'd', varthet(i), dvarthdr(i), d2varthdr2(i), jacrho(i), djacrdrho(i)
         write(*,*) 'e', djacdrho(i), d2jacdr2(i), grho(i)**2, dgr2dr(i), gradthet2(i)
         write(*,*) 'f', dgt2(i), gradrho_gradthet(i), dgrgt(i), gradalph_gradthet(i), dgagt(i)
         write(*,*) 'g', gradrho_gradalph(i), dgagr(i), gradalph2(i), dga2(i), cross(i)
         write(*,*) 'h', dcrossdr(i), gbdrift0(i), dgbdrift0drho(i), cvdrift0(i), dcvdrift0drho(i)
         write(*,*) 'i', gbdrift(i), dgbdriftdrho(i), cvdrift(i), dcvdriftdrho(i), drzdth(i)
         write(*,*) 'j', gradpar(i), dgradpardrho(i), gradparB(i), dgradparBdrho(i), gds2(i)
         write(*,*) 'k', dgds2dr(i), gds21(i), dgds21dr(i), gds22(i), dgds22dr(i), gds23(i), gds24(i)
         write(*,*) 'l', Zr(2, i)
      end if
      do i = -nz, nz
         write (1001, '(59e18.9)') theta(i), Rr(2, i), dRdrho(i), d2Rdr2(i), dRdth(i), &
            d2Rdrdth(i), dZdrho(i), d2Zdr2(i), dZdth(i), d2Zdrdth(i), &
            bmag(i), dBdrho(i), d2Bdr2(i), dBdth(i), d2Bdrdth(i), &
            varthet(i), dvarthdr(i), d2varthdr2(i), jacrho(i), djacrdrho(i), &
            djacdrho(i), d2jacdr2(i), grho(i)**2, dgr2dr(i), gradthet2(i), &
            dgt2(i), gradrho_gradthet(i), dgrgt(i), gradalph_gradthet(i), dgagt(i), &
            gradrho_gradalph(i), dgagr(i), gradalph2(i), dga2(i), cross(i), &
            dcrossdr(i), gbdrift0(i), dgbdrift0drho(i), cvdrift0(i), dcvdrift0drho(i), &
            gbdrift(i), dgbdriftdrho(i), cvdrift(i), dcvdriftdrho(i), drzdth(i), &
            gradpar(i), dgradpardrho(i), gradparB(i), dgradparBdrho(i), gds2(i), &
            dgds2dr(i), gds21(i), dgds21dr(i), gds22(i), dgds22dr(i), gds23(i), gds24(i), &
            Zr(2, i)
      end do
      close (1001)
      if (debug) write (*, *) 'geometry_miller::write_geometry_miller_txt_files_finished'

      defaults_initialized = .false.

   end subroutine get_local_geo

   !============================================================================
   !============================== ALLOCATE ARRAYS =============================
   !============================================================================ 
   subroutine allocate_arrays(nr, nz)

      implicit none

      integer, intent(in) :: nr, nz

      ! periodic quantities can be computed on 2*pi grid and replicated
      allocate (grho(-nz:nz), bmag(-nz:nz), gradpar(-nz:nz)); grho = 0.0; bmag = 0.0; gradpar = 0.0
      allocate (gds2(-nz:nz), gds21(-nz:nz), gds22(-nz:nz), gds23(-nz:nz), gds24(-nz:nz))
      gds2 = 0.0; gds21 = 0.0; gds22 = 0.0; gds23 = 0.0; gds24 = 0.0
      allocate (gbdrift0(-nz:nz), gbdrift(-nz:nz)); gbdrift0 = 0.0; gbdrift = 0.0
      allocate (cvdrift0(-nz:nz), cvdrift(-nz:nz)); cvdrift0 = 0.0; cvdrift = 0.0
      allocate (Rr(nr, -nz:nz), Zr(nr, -nz:nz)); Rr = 0.0; Zr = 0.0
      allocate (jacrho(-nz:nz), djacdrho(-nz:nz), djacrdrho(-nz:nz), d2jacdr2(-nz:nz))
      jacrho = 0.0; djacdrho = 0.0; djacrdrho = 0.0; d2jacdr2 = 0.0
      allocate (d2Rdrdth(-nz:nz), d2Zdrdth(-nz:nz), gpsi(-nz:nz)); d2Rdrdth = 0.0; d2Zdrdth = 0.0; gpsi = 0.0
      allocate (dBdrho(-nz:nz), dgradpardrho(-nz:nz)); dBdrho = 0.0; dgradpardrho = 0.0
      allocate (d2Bdrdth(-nz:nz), dgradparBdrho(-nz:nz), dBdth(-nz:nz), gradparb(-nz:nz))
      d2Bdrdth = 0.0; dgradparBdrho = 0.0; dBdth = 0.0; gradparb = 0.0
      allocate (dcvdrift0drho(-nz:nz), dgbdrift0drho(-nz:nz)); dcvdrift0drho = 0.0; dgbdrift0drho = 0.0
      allocate (theta(-nz:nz)); theta = 0.0
      allocate (gradpararc(-nz:nz)); gradpararc = 0.0
      allocate (arc(-nz:nz)); arc = 0.0
      allocate (dRdrho(-nz:nz), dZdrho(-nz:nz), dRdth(-nz:nz), dZdth(-nz:nz))
      dRdrho = 0.0; dZdrho = 0.0; dRdth  = 0.0; dZdth= 0.0
      allocate (gradrho_gradthet(-nz:nz), gradthet2(-nz:nz), dgr2dr(-nz:nz), dgpsi2dr(-nz:nz))
      gradrho_gradthet = 0.0; gradthet2 = 0.0; dgr2dr = 0.0; dgpsi2dr = 0.0
      allocate (dgrgt(-nz:nz), dgt2(-nz:nz), dgagr(-nz:nz), dgagt(-nz:nz), dga2(-nz:nz))
      dgrgt = 0.0; dgt2 = 0.0; dgagr = 0.0; dgagt = 0.0; dga2 = 0.0
      allocate (d2Rdr2(-nz:nz), d2Zdr2(-nz:nz), d2Bdr2(-nz:nz)); d2Rdr2 = 0.0; d2Zdr2 = 0.0; d2Bdr2 = 0.0
      allocate (drz(-nz:nz), drzdth(-nz:nz), d2Rdr2dth(-nz:nz), d2Zdr2dth(-nz:nz))
      drz = 0.0; drzdth = 0.0; d2Rdr2dth = 0.0; d2Zdr2dth = 0.0
      allocate (d2Rdth2(-nz:nz), d2Zdth2(-nz:nz)); d2Rdth2 = 0.0; d2Zdth2 = 0.0
      allocate (d2gpsidr2(-nz:nz)); d2gpsidr2 = 0.0
      allocate (gradalph_gradthet(-nz:nz), gradalph2(-nz:nz), gradrho_gradalph(-nz:nz))
      gradalph_gradthet = 0.0; gradalph2 = 0.0; gradrho_gradalph = 0.0
      allocate (dgds2dr(-nz:nz), dgds21dr(-nz:nz)); dgds2dr = 0.0; dgds21dr = 0.0
      allocate (dgds22dr(-nz:nz)); dgds22dr = 0.0
      allocate (dcvdriftdrho(-nz:nz), dgbdriftdrho(-nz:nz)); dcvdriftdrho = 0.0; dgbdriftdrho = 0.0
      allocate (varthet(-nz:nz), dvarthdr(-nz:nz), d2varthdr2(-nz:nz)); varthet = 0.0; dvarthdr = 0.0; d2varthdr2 = 0.0
      allocate (cross(-nz:nz)); cross = 0.0
      allocate (dcrossdr(-nz:nz)); dcrossdr = 0.0

   end subroutine allocate_arrays

   !============================================================================
   !============================= DEALLOCATE ARRAYS ============================
   !============================================================================ 
   subroutine deallocate_arrays

      implicit none

      deallocate (grho)
      deallocate (bmag)
      deallocate (gradpar)

      deallocate (gds2)
      deallocate (gds21)
      deallocate (gds22)
      deallocate (gds23)
      deallocate (gds24)
      deallocate (gbdrift0)
      deallocate (gbdrift)
      deallocate (cvdrift0)
      deallocate (cvdrift)
      deallocate (Rr, Zr)
      deallocate (jacrho, djacdrho, djacrdrho, d2jacdr2)
      deallocate (d2Rdrdth, d2Zdrdth, gpsi)
      deallocate (dBdrho, dgradpardrho)
      deallocate (d2Bdrdth, dgradparBdrho, dBdth, gradparb)
      deallocate (dcvdrift0drho, dgbdrift0drho)
      deallocate (theta)
      deallocate (gradpararc)
      deallocate (arc)
      deallocate (dRdrho, dZdrho, dRdth, dZdth)
      deallocate (gradrho_gradthet, gradthet2, dgr2dr, dgpsi2dr)
      deallocate (dgrgt, dgt2, dgagr, dgagt, dga2)
      deallocate (d2Rdr2, d2Zdr2, d2Bdr2)
      deallocate (drz, drzdth, d2Rdr2dth, d2Zdr2dth)
      deallocate (d2Rdth2, d2Zdth2)
      deallocate (d2gpsidr2)
      deallocate (gradalph_gradthet, gradalph2, gradrho_gradalph)
      deallocate (dgds2dr, dgds21dr)
      deallocate (dgds22dr)
      deallocate (dcvdriftdrho, dgbdriftdrho)
      deallocate (varthet, dvarthdr, d2varthdr2)
      deallocate (cross)
      deallocate (dcrossdr)
      deallocate (d2R, d2Z)
      if (allocated(delthet)) deallocate (delthet)
      if (allocated(bmag_psi0)) deallocate (bmag_psi0)
      if (allocated(grho_psi0)) deallocate (grho_psi0)

   end subroutine deallocate_arrays

   !============================================================================
   !================================== FINISH =================================
   !============================================================================ 
   subroutine finish_local_geo

      implicit none

      call deallocate_arrays

   end subroutine finish_local_geo
   
   
   !============================================================================
   !=============================== CALCULATIONS ===============================
   !============================================================================ 

   ! takes in f(r), with r given at three radial locations
   ! and returns df = df/dr at the middle radius
   subroutine get_drho(f, df)

      implicit none

      real, dimension(:, -nz:), intent(in) :: f
      real, dimension(-nz:), intent(out) :: df
      
      df = 0.5 * (f(3, :) - f(1, :)) / local%dr
      
   end subroutine get_drho

   ! given function f(theta), calculate second derivative
   ! of f with respect to theta
   ! second order accurate, with equal grid spacing assumed
   subroutine get_d2dthet2(f, d2f)

      implicit none

      real, dimension(-nz:), intent(in) :: f
      real, dimension(-nz:), intent(out) :: d2f
      
      ! assuming equal grid spacing in theta here
      d2f(-nz + 1:nz - 1) = (f(:nz - 2) - 2.*f(-nz + 1:nz - 1) + f(-nz + 2:)) / delthet(-nz + 1:nz - 1)**2

      ! use periodicity at boundary
      d2f(-nz) = (f(nz - 1) - 2.*f(-nz) + f(-nz + 1)) / delthet(-nz + 1)**2
      d2f(nz) = d2f(-nz)
      
   end subroutine get_d2dthet2

   ! given function f(theta:-pi->pi), calculate theta derivative
   ! second order accurate, with equal grid spacing assumed
   ! assumes periodic in theta -- may need to change this in future
   subroutine get_dthet(f, df)

      implicit none

      real, dimension(-nz:), intent(in) :: f
      real, dimension(-nz:), intent(out) :: df
      
      ! assuming equal grid spacing in theta here
      df(-nz + 1:nz - 1) = (f(-nz + 2:) - f(:nz - 2)) / (delthet(:nz - 2) + delthet(-nz + 1:))

      ! use periodicity at boundary
      df(-nz) = (f(-nz + 1) - f(nz - 1)) / (delthet(-nz) + delthet(nz - 1))
      df(nz) = df(-nz)

   end subroutine get_dthet

   subroutine get_jacrho

      implicit none

      ! jacrho = R*(dR/drho * dZ/dtheta - dR/dtheta * dZ/drho)
      jacrho = Rr(2, :) * (dRdrho * dZdth - dRdth * dZdrho)

   end subroutine get_jacrho

!   ! get dpsinorm/drho = (I/2*pi*q)*int_0^{2*pi} dthet jacrho/R**2
!   subroutine get_dpsipdrho (dpsipdrho)

!     use constants, only: pi

!     implicit none

!     real, intent (out) :: dpsipdrho

!     ! theta_integrate returns integral from 0 -> 2*pi
!     call theta_integrate (jacrho(-nz2pi:nz2pi)/Rr(2,-nz2pi:nz2pi)**2, dpsipdrho)

!     ! integration done using trapezoidal rule
!     dpsipdrho = dpsipdrho*bi/(2.*pi*local%qinp)

!   end subroutine get_dpsipdrho

   subroutine get_gradrho(dpsipdrho, grho)

      implicit none

      real, intent(in) :: dpsipdrho
      real, dimension(-nz:), intent(out) :: grho

      grho = Rr(2, :) * sqrt(dRdth**2 + dZdth**2) / jacrho
      gpsi = grho * dpsipdrho

   end subroutine get_gradrho

   subroutine get_dIdrho(dpsipdrho, grho, dIdrho)

      use constants, only: pi

      implicit none

      real, intent(in) :: dpsipdrho
      real, dimension(-nz:), intent(in) :: grho
      real, intent(out) :: dIdrho

      real :: num1, num2, denom
      real, dimension(:), allocatable :: dum

      allocate (dum(-nz:nz)); dum = 0.

      dum = jacrho * (1.0 + (bi / gpsi)**2) / Rr(2, :)**2
      call theta_integrate(dum(-nz2pi:nz2pi), denom)

      dum = jacrho * (2.*dRdrho / Rr(2, :) + dqdr / local%qinp) / Rr(2, :)**2
      call theta_integrate(dum(-nz2pi:nz2pi), num1)

      ! betaprim below is (4*pi*ptot/B0^2)*(-d ln ptot / drho)
      dum = (-2.*(dRdth * d2Rdrdth + dZdth * d2Zdrdth) / jacrho &
             + drzdth + local%betaprim * jacrho / dpsipdrho**2) / grho**2
      call theta_integrate(dum(-nz2pi:nz2pi), num2)

      dIdrho = bi * (num1 + num2) / denom

      deallocate (dum)

   end subroutine get_dIdrho

   subroutine get_djacdrho(dpsipdrho, dIdrho, grho)

      implicit none

      real, intent(in) :: dpsipdrho, dIdrho
      real, dimension(-nz:), intent(in) :: grho

      ! this is dpsi/dr * d/dr (jacobian)
      ! betaprim below is (4*pi*ptot/B0^2)*(-d ln ptot / drho)
      djacdrho = (Rr(2, :) / grho)**2 * (2.*(dRdth * d2Rdrdth + dZdth * d2Zdrdth) / jacrho &
                                         - drzdth + jacrho * (bi * dIdrho / Rr(2, :)**2 - local%betaprim) / dpsipdrho**2)

      ! this is d/dr (jacobian_r)
      djacrdrho = djacdrho + jacrho * local%d2psidr2 / dpsipdrho

   end subroutine get_djacdrho

   subroutine get_d2RZdr2

      implicit none

      ! get factor common to both d2R/drho2 and d2Z/drho2
      d2Rdr2 = ((djacrdrho - jacrho * dRdrho / Rr(2, :)) / Rr(2, :) &
                - dRdrho * d2Zdrdth + dZdrho * d2Rdrdth) / (dRdth**2 + dZdth**2)

      d2Zdr2 = -d2Rdr2 * dRdth
      d2Rdr2 = d2Rdr2 * dZdth

   end subroutine get_d2RZdr2

   subroutine get_dgr2dr(dpsipdrho, grho)

      implicit none

      real, intent(in) :: dpsipdrho
      real, dimension(-nz:), intent(in) :: grho

      dgr2dr = 2.*(grho**2 * (dRdrho / Rr(2, :) - djacrdrho / jacrho) &
                   + (Rr(2, :) / jacrho)**2 * (dRdth * d2Rdrdth + d2Zdrdth * dZdth))

      dgpsi2dr = 2.*(gpsi**2 * (dRdrho / Rr(2, :) - djacdrho / jacrho) &
                     + (Rr(2, :) / jacrho)**2 * (dRdth * d2Rdrdth + d2Zdrdth * dZdth) * dpsipdrho**2)

   end subroutine get_dgr2dr

   subroutine get_graddotgrad(dpsipdrho, grho)

      implicit none

      real, intent(in) :: dpsipdrho
      real, dimension(-nz:), intent(in) :: grho

      ! grad theta . grad theta
      gradthet2 = (Rr(2, :) / jacrho)**2 * (dRdrho**2 + dZdrho**2)
      ! grad rho . grad theta
      gradrho_gradthet = -(Rr(2, :) / jacrho)**2 * (dRdrho * dRdth + dZdrho * dZdth)

      ! grad alpha . grad theta
      gradalph_gradthet = -(varthet * dqdr + local%qinp * dvarthdr) * gradrho_gradthet &
                          - bi * jacrho / (dpsipdrho * Rr(2, :)**2) * gradthet2
      ! grad rho . grad alpha
      gradrho_gradalph = -(varthet * dqdr + local%qinp * dvarthdr) * grho**2 &
                         - bi * jacrho / (dpsipdrho * Rr(2, :)**2) * gradrho_gradthet
      ! grad alpha . grad alpha
      gradalph2 = (1./Rr(2, :)**2) + ((varthet * dqdr + local%qinp * dvarthdr) * grho)**2 &
                  + 2.*bi * jacrho * (varthet * dqdr + local%qinp * dvarthdr) * gradrho_gradthet / (dpsipdrho * Rr(2, :)**2) &
                  + (bi * jacrho / (dpsipdrho * Rr(2, :)**2))**2 * gradthet2

   end subroutine get_graddotgrad

   subroutine get_gds(gds2, gds21, gds22, gds23, gds24)

      implicit none

      real, dimension(-nz:), intent(out) :: gds2, gds21, gds22, gds23, gds24

      ! |grad alpha|^2 * (dpsiN/drho)^2 (dpsiN/drho factor accounts for ky normalization)
      gds2 = gradalph2 * dpsipdrho_psi0**2
      ! (grad q . grad alpha) * (dpsiN/drho)^2
      gds21 = gradrho_gradalph * dqdr * dpsipdrho_psi0**2
      ! |grad q|^2 * (dpsiN/drho)^2
      gds22 = (grho * dpsipdrho_psi0 * dqdr)**2
      ! (grad rho . grad theta * |grad alpha|^2 - grad alpha . grad theta * grad rho . grad alpha) * (dpsiN/drho)^2 / B^2
      gds23 = (gradrho_gradthet * gradalph2 - gradalph_gradthet * gradrho_gradalph) * (dpsipdrho_psi0 / bmag)**2
      ! (grad rho . grad theta * grad rho . grad alpha - grad alpha . grad theta * |grad rho|^2) * (dpsiN/drho)^2 / B^2 * q/rho
      gds24 = (gradrho_gradthet * gradrho_gradalph - gradalph_gradthet * grho**2) &
              * (dpsipdrho_psi0 / bmag)**2 * (local%qinp_psi0 / local%rhoc_psi0)

      ! note that kperp2 = (n0/a)^2*(drho/dpsiN)^2*(gds2 + 2*theta0*gds21 + theta0^2*gds22)
      ! theta0 = kx/(ky*shat)

   end subroutine get_gds

   subroutine get_dBdrho(bmag, dIdrho)

      implicit none

      real, dimension(-nz:), intent(in) :: bmag
      real, intent(in) :: dIdrho

      ! dB/drho
      dBdrho = (bi * dIdrho + 0.5 * dgpsi2dr) / (bmag * Rr(2, :)**2) &
               - bmag * dRdrho / Rr(2, :)

   end subroutine get_dBdrho

   subroutine get_varthet(dpsipdrho)

      implicit none

      real, intent(in) :: dpsipdrho

      call theta_integrate_indef(jacrho / Rr(2, :)**2, varthet)
      varthet = bi * varthet / (dpsipdrho * local%qinp)

   end subroutine get_varthet

   subroutine get_dvarthdr(dpsipdrho, dIdrho)

      implicit none

      real, intent(in) :: dpsipdrho, dIdrho

      real, dimension(-nz:nz) :: dum

      dum = bi * jacrho * (dIdrho / bi - dqdr / local%qinp + djacdrho / jacrho &
                           - 2.*dRdrho / Rr(2, :)) / Rr(2, :)**2
      call theta_integrate_indef(dum, dvarthdr)
      dvarthdr = dvarthdr / (dpsipdrho * local%qinp)

   end subroutine get_dvarthdr

   subroutine get_d2Idr2_d2jacdr2(grho, dIdrho)

      use constants, only: pi

      implicit none

      real, dimension(-nz:), intent(in) :: grho
      real, intent(in) :: dIdrho

      real :: denom, num1, num2, num3, num4
      real, dimension(-nz:nz) :: tmp, tmp2

      ! denom is the denominator in the expression for d^2 I / dr^2
      tmp = jacrho / Rr(2, :)**2 * (1.0 + (bi / gpsi)**2)
      call theta_integrate(tmp(-nz2pi:nz2pi), denom)
      denom = denom / bi

      d2jacdr2 = dIdrho * bi * jacrho / gpsi**2 &
                 * (dIdrho / bi + djacrdrho / jacrho - dgpsi2dr / gpsi**2 &
                    - 2.*dRdrho / Rr(2, :))

      tmp = -d2jacdr2 / Rr(2, :)**2 - dIdrho * jacrho / (bi * Rr(2, :)**2) &
            * (djacrdrho / jacrho - dIdrho / bi - 2.*dRdrho / Rr(2, :))
      call theta_integrate(tmp(-nz2pi:nz2pi), num1)

      ! tmp = -jacrho/(dpsipdrho*Rr(2,:)**2)*(djacdrho/jacrho - 2.*dRdrho/Rr(2,:))
      ! call theta_integrate (tmp(-nz2pi:nz2pi), num2)
      ! d2jacdr2 = d2jacdr2 - tmp*Rr(2,:)**2*local%d2psidr2
      ! num2 = local%d2psidr2 * (2*pi*local%qinp/bi*(dqdr/local%qinp - dIdrho/bi) + num2)

      tmp = (d2Rdr2 * dRdth + dRdrho * d2Rdrdth + d2Zdr2 * dZdth + dZdrho * d2Zdrdth) / jacrho &
            - djacrdrho * (dRdrho * dRdth + dZdrho * dZdth) / jacrho**2
      call get_dthet(tmp, tmp2)
      tmp = (tmp2 - 2./jacrho * (-djacrdrho / jacrho * (dRdth * d2Rdrdth + dZdth * d2Zdrdth) &
                                 + d2Rdrdth**2 + dRdth * d2Rdr2dth + d2Zdrdth**2 + dZdth * d2Zdr2dth)) / grho**2 &
            - dgr2dr * (drzdth - 2./jacrho * (dRdth * d2Rdrdth + dZdth * d2Zdrdth)) / grho**4
      call theta_integrate(tmp(-nz2pi:nz2pi), num2)
      d2jacdr2 = d2jacdr2 - tmp * Rr(2, :)**2

      tmp = jacrho * (local%betadbprim + local%betaprim * (djacrdrho / jacrho - dgpsi2dr / gpsi**2)) / gpsi**2
      call theta_integrate(tmp(-nz2pi:nz2pi), num3)
      !FLAG - next negative sign?
      d2jacdr2 = d2jacdr2 - tmp * Rr(2, :)**2

      tmp = jacrho / Rr(2, :)**2 * (2.*d2Rdr2 / Rr(2, :) - 2.*(dRdrho / Rr(2, :))**2 &
                                    + local%d2qdr2 / local%qinp - (dqdr / local%qinp)**2 + (2 * dRdrho / Rr(2, :) + dqdr / local%qinp) &
                                    * (djacrdrho / jacrho - 2.*dRdrho / Rr(2, :)))
      call theta_integrate(tmp(-nz2pi:nz2pi), num4)

      d2Idr2 = (num1 + num2 + num3 + num4) / denom
!    d2jacdr2 = d2jacdr2 + bi*jacrho/(gpsi*Rr(2,:))**2*d2Idr2 + 2.*djacdrho*dRdrho/Rr(2,:)**3
      d2jacdr2 = d2jacdr2 + bi * jacrho / gpsi**2 * d2Idr2 + 2.*djacdrho * dRdrho / Rr(2, :)

   end subroutine get_d2Idr2_d2jacdr2

   subroutine get_d2varthdr2(dpsipdrho, dIdrho)

      implicit none

      real, intent(in) :: dpsipdrho, dIdrho

      real, dimension(-nz:nz) :: dum

      dum = bi * jacrho / (local%qinp * dpsipdrho * Rr(2, :)**2) * ((dIdrho / bi - dqdr / local%qinp &
                                                                    !    dum = bi*jacrho/(local%qinp*Rr(2,:)**2)*( (dIdrho/bi - dqdr/local%qinp &
                                                                    + djacdrho / jacrho - 2.*dRdrho / Rr(2, :))**2 &
                                                                   + d2Idr2 / bi - (dIdrho / bi)**2 - local%d2qdr2 / local%qinp &
                                                                   + (dqdr / local%qinp)**2 + d2jacdr2 / jacrho - (djacdrho / jacrho)**2 &
                                                                   - djacdrho * local%d2psidr2 / (dpsipdrho * jacrho) &
                                                                   - 2.*d2Rdr2 / Rr(2, :) + 2.*(dRdrho / Rr(2, :))**2)
      call theta_integrate_indef(dum, d2varthdr2)

   end subroutine get_d2varthdr2

   subroutine get_d2Bdr2(bmag, dIdrho)

      implicit none

      real, dimension(-nz:), intent(in) :: bmag
      real, intent(in) :: dIdrho

      ! d2gpsidr2 = 2.*( dgr2dr*(dRdrho/Rr(2,:) - djacdrho/jacrho) &
      !      + grho**2*(d2Rdr2/Rr(2,:) - (dRdrho/Rr(2,:))**2 - d2jacdr2/jacrho &
      !      + djacdrho*djacrdrho/jacrho**2) + (Rr(2,:)/jacrho)**2 &
      !      * (dRdth**2 + dRdth*d2Rdr2dth + dZdth**2 + dZdth*d2Zdr2dth &
      !      + 2.*(dRdrho/Rr(2,:) - djacrdrho/jacrho)*(dRdth*d2Rdrdth+dZdth*d2Zdrdth)) )
      d2gpsidr2 = 2.*(dRdrho / Rr(2, :) - djacdrho / jacrho) * dgpsi2dr &
                  + 2.*gpsi**2 * (d2Rdr2 / Rr(2, :) - (dRdrho / Rr(2, :))**2 - d2jacdr2 / jacrho + djacdrho * djacrdrho / jacrho**2) &
                  + 2.*(Rr(2, :) * gpsi / jacrho)**2 * (d2Rdrdth**2 + dRdth * d2Rdr2dth + d2Zdrdth**2 + dZdth * d2Zdr2dth &
                                                        + 2.*(dRdth * d2Rdrdth + dZdth * d2Zdrdth) * (dRdrho / Rr(2, :) - djacdrho / jacrho))

      ! d2gpsidr2 = 2.*(dpsipdrho*Rr(2,:)/jacrho)**2 &
      !      * (2.*(dRdrho/Rr(2,:)-djacdrho/jacrho) &
      !      * ((dRdrho/Rr(2,:)-djacdrho/jacrho)*(dRdth**2+dZdth**2) &
      !      + 2.*(dRdth*d2Rdrdth+dZdth*d2Zdrdth)) &
      !      + (dRdth**2+dZdth**2)*(d2rdr2/Rr(2,:) - (dRdrho/Rr(2,:))**2 &
      !      - d2jacdr2/jacrho + (djacdrho/jacrho)**2) &
      !      + d2Rdrdth**2 + dRdth*d2Rdr2dth + d2Zdrdth**2 + dZdth*d2Zdr2dth) &
      !      + 4.*dpsipdrho*local%d2psidr2*dgr2dr &
      !      + 2.*grho**2*(local%d2psidr2**2 + dpsipdrho*local%d3psidr3)

      ! get d/drho (dB/drho)
      d2Bdr2 = -dBdrho * dRdrho / Rr(2, :) + bmag * (dRdrho / Rr(2, :))**2 &
               - bmag * d2Rdr2 / Rr(2, :) + 0.5 * (2.*(dIdrho**2 + bi * d2Idr2) &
                                                   + d2gpsidr2) / (bmag * Rr(2, :)**2) &
               - (dBdrho + bmag * dRdrho / Rr(2, :)) * (2.*dRdrho / Rr(2, :) + dBdrho / bmag)

   end subroutine get_d2Bdr2

   subroutine get_dcrossdr(dpsipdrho, dIdrho, grho)

      implicit none

      real, intent(in) :: dpsipdrho, dIdrho
      real, dimension(-nz:), intent(in) :: grho

      ! dgr2 = d/drho (|grad rho|^2)
      ! dgr2 = 2.*(Rr(2,:)/jacrho)**2*((dRdrho/Rr(2,:)-djacdrho/jacrho)*(dRdth**2+dZdth**2) &
      !      + dRdth*d2Rdrdth + dZdth*d2Zdrdth)
      ! dgrgt = d/drho (grad rho . grad theta)
!    dgrgt = -(Rr(2,:)/jacrho)**2*(2.*(dRdrho/Rr(2,:)-djacdrho/jacrho)*(dRdrho*dRdth+dZdrho*dZdth) &
!         + d2Rdr2*dRdth+dRdrho*d2Rdrdth+d2Zdr2*dZdth+dZdrho*d2Zdrdth)
      dgrgt = 2.*gradrho_gradthet * (dRdrho / Rr(2, :) - djacrdrho / jacrho) &
              - (Rr(2, :) / jacrho)**2 * (d2Rdr2 * dRdth + dRdrho * d2Rdrdth + d2Zdr2 * dZdth + dZdrho * d2Zdrdth)
      ! dgt2 = d/drho (|grad theta|^2)
      dgt2 = 2.*(Rr(2, :) / jacrho)**2 * ((dRdrho / Rr(2, :) - djacrdrho / jacrho) * (dRdrho**2 + dZdrho**2) &
                                          + dRdrho * d2Rdr2 + dZdrho * d2Zdr2)
      ! this is d/drho (|grad alph|^2)
      ! will later multiply it by 0.5*dpsipdrho**2
      dga2 = -2 * dRdrho / Rr(2, :)**3 + dgr2dr * (varthet * dqdr + local%qinp * dvarthdr)**2 &
             + (2.0 * grho**2 * (varthet * dqdr + local%qinp * dvarthdr) &
                + 2.*bi * jacrho * gradrho_gradthet / (dpsipdrho * Rr(2, :)**2)) &
             * (local%d2qdr2 * varthet + 2.*dqdr * dvarthdr + local%qinp * d2varthdr2) &
             + 2.*(varthet * dqdr + local%qinp * dvarthdr) * bi * jacrho / (dpsipdrho * Rr(2, :)**2) &
             * (dgrgt + gradrho_gradthet * (dIdrho / bi + djacdrho / jacrho - 2.*dRdrho / Rr(2, :))) &
             + (bi * jacrho / (dpsipdrho * Rr(2, :)**2))**2 * (dgt2 + 2.*gradthet2 * (dIdrho / bi + djacdrho / jacrho &
                                                                                     - 2.*dRdrho / Rr(2, :)))

      ! dgagr = d/drho (grad alpha . grad rho)
      dgagr = -grho**2 * (2.*dvarthdr * dqdr + varthet * local%d2qdr2 + local%qinp * d2varthdr2) &
              - dgr2dr * (varthet * dqdr + local%qinp * dvarthdr) - bi * jacrho / (dpsipdrho * Rr(2, :)**2) &
              * (dgrgt + gradrho_gradthet * (dIdrho / bi + djacdrho / jacrho - 2.*dRdrho / Rr(2, :)))

      ! dgagt = d/drho (grad alpha . grad theta)
      dgagt = -gradrho_gradthet * (2.*dvarthdr * dqdr + varthet * local%d2qdr2 + local%qinp * d2varthdr2) &
              - dgrgt * (varthet * dqdr + local%qinp * dvarthdr) - bi * jacrho / (dpsipdrho * Rr(2, :)**2) &
              * (dgt2 + gradthet2 * (dIdrho / bi + djacdrho / jacrho - 2.*dRdrho / Rr(2, :)))

      ! dcrossdr = d/drho [(grad alpha x B) . grad theta)]
      dcrossdr = dpsipdrho * (dgagr * gradalph_gradthet + gradrho_gradalph * dgagt &
                             - dga2 * gradrho_gradthet - gradalph2 * dgrgt) + local%d2psidr2 * cross / dpsipdrho

      ! this is (dpsi/drho)^2*d|grad alpha|^2/dr
      dgds2dr = dga2 * dpsipdrho_psi0**2
      ! this is (dpsi/drho)^2*d(grad alpha . grad q)/dr
      ! note that there will be multiplication by 2 in dist_fn.fpp
      dgds21dr = (dgagr * dqdr + local%d2qdr2 * gradrho_gradalph) * dpsipdrho_psi0**2
      ! this is (dpsi/drho)^2*d(|grad q|^2)/dr
      dgds22dr = (dqdr**2 * dgr2dr + 2.*grho**2 * dqdr * local%d2qdr2) * dpsipdrho_psi0**2

      ! note that dkperp2/dr = (n0/a)^2*(drho/dpsiN)^2*(dgds2dr + 2*theta0*dgds21dr + theta0^2*dgds22dr)

   end subroutine get_dcrossdr

   subroutine theta_integrate(integrand, integral)

      implicit none

      real, dimension(-nz2pi:), intent(in) :: integrand
      real, intent(out) :: integral

      ! use trapezoidal rule to integrate in theta
      integral = 0.5 * sum(delthet(-nz2pi:nz2pi - 1) * (integrand(-nz2pi:nz2pi - 1) + integrand(-nz2pi + 1:nz2pi)))

   end subroutine theta_integrate

   ! get indefinite integral of integrand
   subroutine theta_integrate_indef(integrand, integral)

      implicit none

      real, dimension(-nz:), intent(in) :: integrand
      real, dimension(-nz:), intent(out) :: integral

      integer :: i
      
      ! use trapezoidal rule to integrate in theta
      integral(0) = 0.0
      do i = 1, nz
         integral(i) = integral(i - 1) + 0.5 * delthet(i - 1) * (integrand(i - 1) + integrand(i))
      end do
      do i = -1, -nz, -1
         integral(i) = integral(i + 1) - 0.5 * delthet(i) * (integrand(i + 1) + integrand(i))
      end do
      
   end subroutine theta_integrate_indef

   function Rpos(r, theta, j)

      use constants, only: pi

      integer, intent(in) :: j
      real, intent(in) :: r, theta
      real :: Rpos
      real :: g, gp, dr
      integer :: i

      dr = r - local%rhoc

! For Y Xiao:
!    g = local%delp/local%rhoc + local%d * sin(theta)**2
!    Rpos = local%rmaj*(1.+r*(cos(theta)-g)-g*dr)

      g = cos(theta + local%tri * sin(theta))
      gp = -sin(theta + local%tri * sin(theta)) &
           * local%triprim * sin(theta)

      ! allow for strange specification of R_psi
      if (j == nz + 1) then
         i = -nz
      else
         i = j
      end if

      ! second line here is (1/2)*(r-r0)**2*d2R/dr|_r0
      ! note that d2R=0 unless read_profile_variation = T in input file
      Rpos = local%rmaj + local%shift * dr + g * local%rhoc + (g + local%rhoc * gp) * dr &
             + 0.5 * (r - rhoc0)**2 * d2R(i)

   end function Rpos

   function Zpos(r, theta, j)

      integer, intent(in) :: j
      real, intent(in) :: r, theta
      real :: Zpos, dr
      integer :: i

      ! allow for strange specification of Z_psi
      if (j == nz + 1) then
         i = -nz
      else
         i = j
      end if

      dr = r - local%rhoc
      ! note that d2Z=0 unless read_profile_variation=T in input file
      Zpos = local%kappa * sin(theta) * local%rhoc + (local%rhoc * local%kapprim + local%kappa) * sin(theta) * dr &
             + 0.5 * (r - rhoc0)**2 * d2Z(i)

   end function Zpos

   function mod2pi(theta)

      real, intent(in) :: theta
      real :: pi, th, mod2pi
      real, parameter :: theta_tol = 1.e-6
      logical :: out

      pi = 2.*acos(0.)

      if (theta <= pi .and. theta >= -pi) then
         mod2pi = theta
         return
      end if

      if (theta - theta_tol <= pi .and. theta >= -pi) then
         mod2pi = pi
         return
      end if

      if (theta <= pi .and. theta + theta_tol >= -pi) then
         mod2pi = -pi
         return
      end if

      th = theta
      out = .true.
      do while (out)
         if (th > pi) th = th - 2.*pi
         if (th < -pi) th = th + 2.*pi
         if (th <= pi .and. th >= -pi) out = .false.
      end do
      mod2pi = th

   end function mod2pi

   function round(val, n)
     implicit none
     real :: val, round
     real :: scaled, remainder
     integer :: n
     integer :: sgn

     scaled = val*10.0**n
     sgn = sign(1.0, scaled)
     remainder = modulo(abs(scaled), 10.0)
     round = (scaled - sgn * remainder )/ 10.0**n
     !aint(val*10.0**n)/10.0**n
   end function round

end module geometry_miller