grid_utilities.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------!
!> Numerical utilities related to the grids and different coordinate systems.
!------------------------------------------------------------------------------!
module grid_utilities
#include <PB3D_macros.h>
#include <wrappers.h>
    use str_utilities
    use output_ops
    use messages
    use num_vars, only: dp, pi, max_str_ln, iu
    use grid_vars, only: grid_type
 
    implicit none
    private
    public coord_f2e, coord_e2f, calc_xyz_grid, calc_eqd_grid, extend_grid_f, &
        &calc_int_vol, copy_grid, trim_grid, untrim_grid, calc_n_par_x_rich, &
        &calc_vec_comp, nufft, find_compr_range, calc_arc_angle, calc_tor_diff
#if ldebug
    public debug_calc_int_vol, debug_calc_vec_comp
#endif
    
    ! global variables
#if ldebug
    !> \ldebug
    logical :: debug_calc_int_vol = .false.                                     !< plot debug information for calc_int_vol()
    !> \ldebug
    logical :: debug_calc_vec_comp = .false.                                    !< plot debug information for calc_vec_comp()
#endif
    
    ! interfaces
    
    !> \public Converts Flux coordinates
    !!      \f$\left(r,\theta,\zeta\right)_\text{F}\f$
    !!      to Equilibrium coordinates
    !!      \f$\left(r,\theta,\zeta\right)_\text{E}\f$.
    !!
    !! Optionally, two  arrays \c  r_F_array and \c  r_E_array can  be provided,
    !! that define the mapping between the both coordinate system.
    !! 
    !! Standard, for E, the poloidal or toroidal normalized flux is used and for
    !! F, the poloidal or toroidal flux in F coordinates, divided by \f$2\pi\f$.
    !!
    !! \note For VMEC,  it can be slow,  as the zero of  a non-linear expression
    !! must be sought. This is done currently using num_ops.calc_zero_hh().
    !!
    !! \return ierr
    interface coord_f2e
        !> \public
        module procedure coord_f2e_r
        !> \public
        module procedure coord_f2e_rtz
    end interface
    
    !> \public Converts Equilibrium coordinates
    !!      \f$\left(r,\theta,\zeta\right)_\text{E}\f$.
    !!      to Flux coordinates
    !!      \f$\left(r,\theta,\zeta\right)_\text{F}\f$
    !!
    !! Optionally, two  arrays \c  r_E_array and \c  r_F_array can  be provided,
    !! that define the mapping between the both coordinate system.
    !!
    !! Standard, for E, the poloidal or toroidal normalized flux is used and for
    !! F, the poloidal or toroidal flux in F coordinates, divided by \f$2\pi\f$.
    !!
    !! \return ierr
    interface coord_e2f
        !> \public
        module procedure coord_e2f_r
        !> \public
        module procedure coord_e2f_rtz
    end interface
    
    !> \public Calculate grid  of equidistant points, where  optionally the last
    !! point can be excluded.
    interface calc_eqd_grid
        !> \public
        module procedure calc_eqd_grid_1d
        !> \public
        module procedure calc_eqd_grid_3d
    end interface
    
    !> \public Calculates the toroidal difference  for a magnitude calculated on
    !! three toroidal points: two extremities and one in the middle.
    !!
    !! This is done using the formula
    !!  \f[ \frac{b-a}{b+a} , \f]
    !! for the relative  difference between \f$a\f$ and \f$b\f$.  This is useful
    !! when it is used to calculate a  toroidal ripple and these are the extreme
    !! points.
    !! 
    !! The procedure also needs the map  between the flux poloidal angle and the
    !! geometrical poloidal angle.
    !!
    !! In a  first step the quantity  is interpolated on an  equidistant grid in
    !! the geometrical poloidal angle.
    !!
    !! The  difference  is then  calculated  for  values of  constant  geometric
    !! poloidal angle.
    !!
    !! Finally, this result is transformed back to the Flux coordinates.
    !!
    !! \note
    !!  -# The theta map should have first and last point equal.
    !!  -# This  routine should  be used  only for  periodic quantities  (so not
    !!  for  some  of the  Flux  quantities  that  are defined  along  generally
    !!  non-rational magnetic field lines).
    !!
    !! \return ierr
    interface calc_tor_diff
        !> \public
        module procedure calc_tor_diff_0d
        !> \public
        module procedure calc_tor_diff_2d
    end interface
    
contains
    !> \private version withi r, theta and zeta
    integer function coord_f2e_rtz(grid_eq,r_F,theta_F,zeta_F,r_E,&
        &theta_E,zeta_E,r_F_array,r_E_array,ord) result(ierr)
        use num_vars, only: tol_zero, eq_style
        use vmec_utilities, only: fourier2real
        use vmec_vars, only: l_v_c, l_v_s, is_asym_v
        
        character(*), parameter :: rout_name = 'coord_F2E_rtz'
        
        ! input / output
        type(grid_type), intent(in) :: grid_eq                                  !< equilibrium grid (for normal local limits and possibly loc_r)
        real(dp), intent(in) :: r_f(:)                                          !< (local) \f$r_\text{F}\f$
        real(dp), intent(in) :: theta_f(:,:,:)                                  !< \f$\theta_\text{F}\f$
        real(dp), intent(in) :: zeta_f(:,:,:)                                   !< \f$\zeta_\text{F}\f$
        real(dp), intent(inout) :: r_e(:)                                       !< (local) \f$r_\text{E}\f$
        real(dp), intent(inout) :: theta_e(:,:,:)                               !< \f$\theta_\text{E}\f$
        real(dp), intent(inout) :: zeta_e(:,:,:)                                !< \f$\zeta_\text{E}\f$
        real(dp), intent(in), optional, target :: r_f_array(:)                  !< optional array that defines mapping between two coord. systems
        real(dp), intent(in), optional, target :: r_e_array(:)                  !< optional array that defines mapping between two coord. systems
        integer, intent(in), optional :: ord                                    !< order (only used for eq_style 1)
        
        ! local variables (also used in child functions)
        character(len=max_str_ln) :: err_msg                                    ! error message
        integer :: dims(3)                                                      ! dimensions of the grid
        real(dp), allocatable :: l_v_c_loc(:,:)                                 ! local version of L_V_c
        real(dp), allocatable :: l_v_s_loc(:,:)                                 ! local version of L_V_s
        integer :: ord_loc                                                      ! local ord
        
        ! initialize ierr
        ierr = 0
        
        ! set up array sizes
        dims = shape(theta_e)
        
        ! set up local ord
        ord_loc = 1
        if (present(ord)) ord_loc = ord
        
        ! tests
        if (dims(3).ne.size(r_f) .or. dims(3).ne.size(r_e)) then
            ierr = 1
            err_msg = 'r_F and r_E need to have the correct dimensions'
            chckerr(err_msg)
        end if
        if (present(r_f_array).neqv.present(r_e_array)) then
            ierr = 1
            err_msg = 'both r_F_array and r_E_array have to be provided'
            chckerr(err_msg)
        end if
        
        ! convert normal position
        ierr = coord_f2e_r(grid_eq,r_f,r_e,r_f_array,r_e_array)
        chckerr('')
        
        ! choose which equilibrium style is being used:
        !   1:  VMEC
        !   2:  HELENA
        select case (eq_style)
            case (1)                                                            ! VMEC
                ierr = coord_f2e_vmec()
                chckerr('')
            case (2)                                                            ! HELENA
                ! trivial HELENA uses flux coordinates
                theta_e = theta_f
                zeta_e = zeta_f
        end select
    contains
        ! VMEC version
        !> \private
        integer function coord_f2e_vmec() result(ierr)
            use num_vars, only: norm_disc_prec_eq, rank
            use num_ops, only: calc_zero_hh
            use vmec_vars, only: mnmax_v
            use num_utilities, only: spline
            
            character(*), parameter :: rout_name = 'coord_F2E_VMEC'
            
            ! local variables
            integer :: id                                                       ! counter
            real(dp), pointer :: loc_r_f(:) => null()                           ! loc_r in F coords.
            real(dp), allocatable :: theta_e_guess_3d(:,:,:)                    ! guess for theta_E
            real(dp), allocatable :: lam(:,:,:,:)                               ! lambda and derivative
            
            ! initialize ierr
            ierr = 0
            
            ! set up loc_r_F
            if (present(r_f_array)) then
                loc_r_f => r_f_array
            else
                loc_r_f => grid_eq%loc_r_F
            end if
            
            ! the toroidal angle is trivial
            zeta_e = - zeta_f                                                   ! conversion VMEC LH -> RH coord. system
            
            ! allocate local copies of L_V_c and L_V_s
            allocate(l_v_c_loc(mnmax_v,dims(3)))
            allocate(l_v_s_loc(mnmax_v,dims(3)))
            
            ! interpolate L_V_c and L_V_s at requested normal points r_F
            do id = 1,mnmax_v
                ierr = spline(loc_r_f,l_v_c(id,grid_eq%i_min:grid_eq%i_max,0),&
                    &r_f,l_v_c_loc(id,:),ord=norm_disc_prec_eq)
                chckerr('')
                ierr = spline(loc_r_f,l_v_s(id,grid_eq%i_min:grid_eq%i_max,0),&
                    &r_f,l_v_s_loc(id,:),ord=norm_disc_prec_eq)
                chckerr('')
            end do
            
            ! set up guess
            allocate(theta_e_guess_3d(dims(1),dims(2),dims(3)))
            allocate(lam(dims(1),dims(2),dims(3),0:1))
            
            ! set up guess: try theta_F - lambda(theta_F)/(1+Dlambda(theta_F)),
            ! which is what you get when you approximate lambda(theta_V) as
            ! lambda(theta_F) + Dlambda(theta_F) (theta_V-theta_F)
            ! (lambda(theta_V) would be the exact solution)
            do id = 0,1
                ierr = fourier2real(l_v_c_loc,l_v_s_loc,theta_f,zeta_e,&
                    &lam(:,:,:,id),sym=[is_asym_v,.true.],deriv=[id,0])
                chckerr('')
            end do
            theta_e_guess_3d = theta_f - lam(:,:,:,0)/(1+lam(:,:,:,1))
            
            ! the poloidal angle has to be found as the zero of
            !   f = theta_F - theta_E - lambda
            err_msg = calc_zero_hh(dims,theta_e,fun_pol,ord_loc,&
                &theta_e_guess_3d,output=.true.)
            if (err_msg.ne.'') then
                ierr = 1
                chckerr(err_msg)
            end if
            
            ! do a check whether the result is indeed zero
            if (maxval(abs(fun_pol(dims,theta_e,0))).gt.tol_zero*100) then
                call plot_hdf5('f','ERROR_coord_F2E_VMEC_'//trim(i2str(rank)),&
                    &fun_pol(dims,theta_e,0))
                err_msg = 'In coord_F2E_VMEC, calculating whether f=0 as a &
                    &check, yields a deviation from 0 of '//trim(r2strt(&
                    &maxval(abs(fun_pol(dims,theta_e,0)))))//'. See output.'
                ierr = 1
                chckerr(err_msg)
            end if
            
            ! clean up
            nullify(loc_r_f)
        end function coord_f2e_vmec
        
        ! function  that  returns  f  =  theta_F  -  theta_V  -  lambda  or  its
        ! derivatives in the poloidal direction.  It uses zeta_E (= zeta_V) from
        ! the parent function.
        !> \private
        function fun_pol(dims,theta_E_in,dpol)
            character(*), parameter :: rout_name = 'fun_pol'
            
            ! input / output
            integer, intent(in) :: dims(3)
            real(dp), intent(in) :: theta_e_in(dims(1),dims(2),dims(3))
            integer, intent(in) :: dpol
            real(dp) :: fun_pol(dims(1),dims(2),dims(3))
            
            ! local variables
            real(dp) :: lam(dims(1),dims(2),dims(3))
            
            ! initialize fun_pol
            fun_pol = 0.0_dp
            
            ! calculate lambda
            ierr = fourier2real(l_v_c_loc,l_v_s_loc,theta_e_in,zeta_e,lam,&
                &sym=[is_asym_v,.true.],deriv=[dpol,0])
            chckerr('')
            
            ! calculate the output function
            if (dpol.eq.0) then
                fun_pol = theta_f - theta_e_in - lam
            else if (dpol.eq.1) then
                fun_pol = -1 - lam
            else if (dpol.gt.1) then
                fun_pol = - lam
            end if
        end function fun_pol
    end function coord_f2e_rtz
    !> \private version with only r
    integer function coord_f2e_r(grid_eq,r_F,r_E,r_F_array,r_E_array) &
        &result(ierr)                                                           ! version with only r
        use num_vars, only: norm_disc_prec_eq
        use num_utilities, only: spline
        
        character(*), parameter :: rout_name = 'coord_F2E_r'
        
        ! input / output
        type(grid_type), intent(in) :: grid_eq                                  !< equilibrium grid (for possibly loc_r)
        real(dp), intent(in) :: r_f(:)                                          !< (local) \f$r_\text{F}\f$
        real(dp), intent(inout) :: r_e(:)                                       !< (local) \f$r_\text{E}\f$
        real(dp), intent(in), optional, target :: r_f_array(:)                  !< optional array that defines mapping between two coord. systems
        real(dp), intent(in), optional, target :: r_e_array(:)                  !< optional array that defines mapping between two coord. systems
        
        ! local variables
        character(len=max_str_ln) :: err_msg                                    ! error message
        integer :: n_r                                                          ! dimension of the grid
        real(dp), pointer :: loc_r_e(:) => null()                               ! loc_r in E coords.
        real(dp), pointer :: loc_r_f(:) => null()                               ! loc_r in F coords.
        
        ! initialize ierr
        ierr = 0
        
        ! set up array sizes
        n_r = size(r_f)
        
        ! tests
        if (n_r.ne.size(r_e)) then
            ierr = 1
            err_msg = 'r_F and r_E need to have the correct dimensions'
            chckerr(err_msg)
        end if
        if (present(r_f_array).neqv.present(r_e_array)) then
            ierr = 1
            err_msg = 'both r_F_array and r_E_array have to be provided'
            chckerr(err_msg)
        end if
        
        ! set up loc_r_E and loc_r_F
        if (present(r_f_array)) then
            loc_r_e => r_e_array
            loc_r_f => r_f_array
        else
            loc_r_e => grid_eq%loc_r_E
            loc_r_f => grid_eq%loc_r_F
        end if
        
        ! convert normal position
        ierr = spline(loc_r_f,loc_r_e,r_f,r_e,ord=norm_disc_prec_eq)
        chckerr('')
        
        ! clean up
        nullify(loc_r_e,loc_r_f)
    end function coord_f2e_r
    
    !> version with r, theta and zeta
    integer function coord_e2f_rtz(grid_eq,r_E,theta_E,zeta_E,r_F,&
        &theta_F,zeta_F,r_E_array,r_F_array) result(ierr)                       ! version with r, theta and zeta
        use num_vars, only: eq_style
        
        character(*), parameter :: rout_name = 'coord_E2F_rtz'
        
        ! input / output
        type(grid_type), intent(in) :: grid_eq                                  !< equilibrium grid (for normal local limits)
        real(dp), intent(in) :: r_e(:)                                          !< (local) \f$r_\text{E}\f$
        real(dp), intent(in) :: theta_e(:,:,:)                                  !< \f$\theta_\text{E}\f$
        real(dp), intent(in) :: zeta_e(:,:,:)                                   !< \f$\zeta_\text{E}\f$
        real(dp), intent(inout) :: r_f(:)                                       !< (local) \f$r_\text{F}\f$
        real(dp), intent(inout) :: theta_f(:,:,:)                               !< \f$\theta_\text{F}\f$
        real(dp), intent(inout) :: zeta_f(:,:,:)                                !< \f$\zeta_\text{F}\f$
        real(dp), intent(in), optional, target :: r_e_array(:)                  !< optional array that defines mapping between two coord. systems
        real(dp), intent(in), optional, target :: r_f_array(:)                  !< optional array that defines mapping between two coord. systems
        
        ! local variables (also used in child functions)
        character(len=max_str_ln) :: err_msg                                    ! error message
        integer :: dims(3)                                                      ! dimensions of the grid
        
        ! initialize ierr
        ierr = 0
        
        ! set up array sizes
        dims = shape(theta_f)
        
        ! tests
        if (dims(3).ne.size(r_f) .or. dims(3).ne.size(r_e)) then
            ierr = 1
            err_msg = 'r_F and r_E need to have the correct dimensions'
            chckerr(err_msg)
        end if
        if (present(r_f_array).neqv.present(r_e_array)) then
            ierr = 1
            err_msg = 'both r_F_array and r_E_array have to be provided'
            chckerr(err_msg)
        end if
        
        ! convert normal position
        ierr = coord_e2f_r(grid_eq,r_e,r_f,r_e_array,r_f_array)
        chckerr('')
        
        ! choose which equilibrium style is being used:
        !   1:  VMEC
        !   2:  HELENA
        select case (eq_style)
            case (1)                                                            ! VMEC
                ierr = coord_e2f_vmec()
                chckerr('')
            case (2)                                                            ! HELENA
                ! trivial HELENA uses flux coordinates
                theta_f = theta_e
                zeta_f = zeta_e
        end select
    contains
        ! VMEC version
        !> \private
        integer function coord_e2f_vmec() result(ierr)
            use num_vars, only: norm_disc_prec_eq
            use vmec_utilities, only: fourier2real
            use vmec_vars, only: mnmax_v, l_v_c, l_v_s, is_asym_v
            use num_utilities, only: spline
            
            character(*), parameter :: rout_name = 'coord_E2F_VMEC'
            
            ! local variables
            integer :: id                                                       ! counter
            real(dp), allocatable :: l_v_c_loc(:,:), l_v_s_loc(:,:)             ! local version of L_V_c and L_V_s
            real(dp), allocatable :: lam(:,:,:)                                 ! lambda
            real(dp), pointer :: loc_r_e(:) => null()                           ! loc_r in E coords.
            
            ! initialize ierr
            ierr = 0
            
            ! set up loc_r_E
            if (present(r_e_array)) then
                loc_r_e => r_e_array
            else
                loc_r_e => grid_eq%loc_r_E
            end if
            
            ! the toroidal angle is trivial
            zeta_f = - zeta_e                                                   ! conversion VMEC LH -> RH coord. system
            
            ! allocate local copies of L_V_c and L_V_s and lambda
            allocate(l_v_c_loc(mnmax_v,dims(3)))
            allocate(l_v_s_loc(mnmax_v,dims(3)))
            allocate(lam(dims(1),dims(2),dims(3)))
            
            ! interpolate L_V_c and L_V_s at requested normal points r_E
            do id = 1,mnmax_v
                ierr = spline(loc_r_e,l_v_c(id,grid_eq%i_min:grid_eq%i_max,0),&
                    &r_e,l_v_c_loc(id,:),ord=norm_disc_prec_eq)
                chckerr('')
                ierr = spline(loc_r_e,l_v_s(id,grid_eq%i_min:grid_eq%i_max,0),&
                    &r_e,l_v_s_loc(id,:),ord=norm_disc_prec_eq)
                chckerr('')
            end do
            
            ! calculate lambda
            ierr = fourier2real(l_v_c_loc,l_v_s_loc,theta_e,zeta_e,lam,&
                &sym=[is_asym_v,.true.])
            chckerr('')
            
            ! the poloidal angle has to be found as
            !   theta_F = theta_E + lambda
            theta_f = theta_e + lam
            
            ! clean up
            nullify(loc_r_e)
        end function coord_e2f_vmec
    end function coord_e2f_rtz
    !> version with only r
    integer function coord_e2f_r(grid_eq,r_E,r_F,r_E_array,r_F_array) &
        &result(ierr)                                                           ! version with only r
        use num_vars, only: norm_disc_prec_eq
        use num_utilities, only: spline
        
        character(*), parameter :: rout_name = 'coord_E2F_r'
        
        ! input / output
        type(grid_type), intent(in) :: grid_eq                                  !< equilibrium grid (for normal local limits)
        real(dp), intent(in) :: r_e(:)                                          !< (local) \f$r_\text{E}\f$
        real(dp), intent(inout) :: r_f(:)                                       !< (local) \f$r_\text{F}\f$
        real(dp), intent(in), optional, target :: r_e_array(:)                  !< optional array that defines mapping between two coord. systems
        real(dp), intent(in), optional, target :: r_f_array(:)                  !< optional array that defines mapping between two coord. systems
        
        ! local variables
        character(len=max_str_ln) :: err_msg                                    ! error message
        integer :: n_r                                                          ! dimension of the grid
        real(dp), pointer :: loc_r_e(:) => null()                               ! loc_r in E coords.
        real(dp), pointer :: loc_r_f(:) => null()                               ! loc_r in F coords.
        
        ! initialize ierr
        ierr = 0
        
        ! set up array sizes
        n_r = size(r_e)
        
        ! tests
        if (n_r.ne.size(r_f)) then
            ierr = 1
            err_msg = 'r_E and r_F need to have the correct dimensions'
            chckerr(err_msg)
        end if
        if (present(r_e_array).neqv.present(r_f_array)) then
            ierr = 1
            err_msg = 'both r_E_array and r_F_array have to be provided'
            chckerr(err_msg)
        end if
        
        ! set up loc_r_E and loc_r_F
        if (present(r_f_array)) then
            loc_r_e => r_e_array
            loc_r_f => r_f_array
        else
            loc_r_e => grid_eq%loc_r_E
            loc_r_f => grid_eq%loc_r_F
        end if
        
        ! convert normal position
        ierr = spline(loc_r_e,loc_r_f,r_e,r_f,ord=norm_disc_prec_eq)
        chckerr('')
        
        ! clean up
        nullify(loc_r_e,loc_r_f)
    end function coord_e2f_r
 
    !> \private 3-D version
    integer function calc_eqd_grid_3d(var,min_grid,max_grid,grid_dim,&
        &excl_last) result(ierr)
        character(*), parameter :: rout_name = 'calc_eqd_grid_3D'
        
        ! input and output
        real(dp), intent(inout) :: var(:,:,:)                                   !< output
        real(dp), intent(in) :: min_grid                                        !< min. of angles [\f$\pi\f$]
        real(dp), intent(in) :: max_grid                                        !< max. of angles [\f$\pi\f$]
        integer, intent(in) :: grid_dim                                         !< in which dimension to create the grid
        logical, intent(in), optional :: excl_last                              !< .true. if last point excluded
        
        ! local variables
        integer :: id                                                           ! counter
        character(len=max_str_ln) :: err_msg                                    ! error message
        logical :: excl_last_loc                                                ! local copy of excl_last
        integer :: grid_size                                                    ! nr. of points
        integer :: grid_size_mod                                                ! grid_size or grid_size + 1
        
        ! initialize ierr
        ierr = 0
        
        ! tests
        if (grid_dim.lt.1 .or. grid_dim.gt.3) then
            ierr = 1
            err_msg = 'grid_dim has to point to a dimension going from 1 to 3'
            chckerr(err_msg)
        end if
        
        ! set up grid_size
        grid_size = size(var,grid_dim)
        
        ! test some values
        if (grid_size.lt.1) then
            err_msg = 'The angular array has to have a length of at &
                &least 1'
            ierr = 1
            chckerr(err_msg)
        end if
        
        ! set up local excl_last
        excl_last_loc = .false.
        if (present(excl_last)) excl_last_loc = excl_last
        
        ! add 1 to modified grid size if last point is to be excluded
        if (excl_last_loc) then
            grid_size_mod = grid_size + 1
        else
            grid_size_mod = grid_size
        end if
        
        ! initialize output vector
        var = 0.0_dp
        
        ! calculate grid points
        if (grid_size.eq.1) then
            var = min_grid                                                      ! = max_grid
        else
            if (grid_dim.eq.1) then
                do id = 1,grid_size
                    var(id,:,:) = min_grid + &
                        &(max_grid-min_grid)*(id-1)/(grid_size_mod-1)
                end do
            else if (grid_dim.eq.2) then
                do id = 1,grid_size
                    var(:,id,:) = min_grid + &
                        &(max_grid-min_grid)*(id-1)/(grid_size_mod-1)
                end do
            else
                do id = 1,grid_size
                    var(:,:,id) = min_grid + &
                        &(max_grid-min_grid)*(id-1)/(grid_size_mod-1)
                end do
            end if
        end if
    end function calc_eqd_grid_3d
    !> \private 1-D version
    integer function calc_eqd_grid_1d(var,min_grid,max_grid,excl_last) &
        &result(ierr)
        character(*), parameter :: rout_name = 'calc_eqd_grid_1D'
        
        ! input and output
        real(dp), intent(inout) :: var(:)                                       !< output
        real(dp), intent(in) :: min_grid                                        !< min. of angles [\f$\pi\f$]
        real(dp), intent(in) :: max_grid                                        !< max. of angles [\f$\pi\f$]
        logical, intent(in), optional :: excl_last                              !< .true. if last point excluded
        
        ! local variables
        real(dp), allocatable :: var_3d(:,:,:)                                  ! 3D version of var
        
        ! initialize ierr
        ierr = 0
        
        ! set up var_3D
        allocate(var_3d(size(var),1,1))
        
        ! call 3D version
        ierr = calc_eqd_grid_3d(var_3d,min_grid,max_grid,1,excl_last)
        chckerr('')
        
        ! update var
        var = var_3d(:,1,1)
    end function calc_eqd_grid_1d
 
    !> \private 2-D version
    integer function calc_tor_diff_2d(v_com,theta,norm_disc_prec,absolute,r) &
        &result(ierr)
        use num_vars, only: min_theta_plot, max_theta_plot, tol_zero
        use num_utilities, only: spline, order_per_fun
        
        character(*), parameter :: rout_name = 'calc_tor_diff_2D'
        
        ! input / output
        real(dp), intent(inout) :: v_com(:,:,:,:,:)                             !< covariant and contravariant components <tt>(dim1,dim2,dim3,3,2)</tt>
        real(dp), intent(in) :: theta(:,:,:)                                    !< geometric poloidal angle
        integer, intent(in) :: norm_disc_prec                                   !< precision for normal derivatives
        logical, intent(in), optional :: absolute                               !< calculate absolute, not relative, difference
        real(dp), intent(in), optional :: r(:)                                  !< normal positions for theta
        
        ! local variables
        integer :: id, jd, kd, ld                                               ! counters
        integer :: n_pol                                                        ! number of poloidal points to be used (1 less than total)
        integer :: istat                                                        ! status
        real(dp), allocatable :: theta_eqd(:)                                   ! equidistant grid
        real(dp), allocatable :: v_com_interp(:,:)                              ! interpolated local v_com for outer points
        real(dp), allocatable :: theta_ord(:)                                   ! ordered theta
        real(dp), allocatable :: v_com_ord(:)                                   ! ordered v_com
        logical :: absolute_loc                                                 ! local absolute
        
        ! initialize ierr
        ierr = 0
        
        ! set up variables
        n_pol = size(theta,1)-1
        allocate(theta_eqd(n_pol))
        allocate(v_com_interp(n_pol,3))
        absolute_loc = .false.
        if (present(absolute)) absolute_loc = absolute
        
        norm: do kd = 1,size(v_com,3)
            ! set up equidistant grid
            ierr = calc_eqd_grid(theta_eqd,min_theta_plot*pi,max_theta_plot*pi,&
                &excl_last=.true.)
            chckerr('')
            do id = 1,size(v_com,4)
                do ld = 1,size(v_com,5)
                    ! interpolate the geometric poloidal angle for outer points
                    do jd = 1,3,2
                        ! order
                        ierr = order_per_fun(theta(1:n_pol,jd,kd),&
                            &v_com(1:n_pol,jd,kd,id,ld),theta_ord,v_com_ord,&
                            &norm_disc_prec)
                        chckerr('')
                        
                        istat = spline(theta_ord,v_com_ord,theta_eqd,&
                            &v_com_interp(:,jd),ord=min(norm_disc_prec,3),&
                            &deriv=0)
                        deallocate(theta_ord,v_com_ord)
                        
                        if (istat.ne.0) then
                            call display_interp_warning(r)
                            v_com(:,2,kd,:,:) = 0._dp
                            cycle norm
                        end if
                    end do
                    
                    ! calculate difference and save in middle point
                    v_com_interp(:,2) = &
                        &(v_com_interp(:,3)-v_com_interp(:,1))
                    if (.not.absolute_loc) &                                    ! make it relative
                        &v_com_interp(:,2) = v_com_interp(:,2)/&
                        &sign(max(&
                        &tol_zero,abs(v_com_interp(:,3)+v_com_interp(:,1))),&
                        &v_com_interp(:,3)+v_com_interp(:,1))
                    
                    ! order
                    ierr = order_per_fun(theta_eqd,v_com_interp(:,2),&
                        &theta_ord,v_com_ord,norm_disc_prec)
                    chckerr('')
                    
                    ! interpolate back
                    istat = spline(theta_ord,v_com_ord,theta(:,2,kd),&
                        &v_com(:,2,kd,id,ld),min(norm_disc_prec,3),&
                        &deriv=0)
                    deallocate(theta_ord,v_com_ord)
                    
                    if (istat.ne.0) then
                        call display_interp_warning(r)
                        v_com(:,2,kd,:,:) = 0._dp
                        cycle norm
                    end if
                end do
            end do
        end do norm
    contains
        ! Displays warning if no interpolation possible.
        ! Uses variable kd from parent procedure.
        !> \private
        subroutine display_interp_warning(r)
            ! input / output
            real(dp), intent(in), optional :: r(:)                              ! normal positions for theta
            
            if (present(r)) then
                call writo('Cannot interpolate geometrical &
                    &theta at normal position '//&
                    &trim(r2str(r(kd))),warning=.true.)
            else
                call writo('Cannot interpolate geometrical &
                    &theta',warning=.true.)
            end if
            call lvl_ud(1)
            call writo('Are you sure the origin is chosen &
                &correctly?')
            call writo('Skipping this normal position')
            call lvl_ud(-1)
        end subroutine display_interp_warning
    end function calc_tor_diff_2d
    !> \private 0-D version
    integer function calc_tor_diff_0d(v_mag,theta,norm_disc_prec,absolute,r) &
        &result(ierr)
        character(*), parameter :: rout_name = 'calc_tor_diff_0D'
        
        ! input / output
        real(dp), intent(inout) :: v_mag(:,:,:)                                 !< magnitude <tt>(dim1,dim2,dim3)</tt>
        real(dp), intent(in) :: theta(:,:,:)                                    !< geometric poloidal angle
        integer, intent(in) :: norm_disc_prec                                   !< precision for normal derivatives
        logical, intent(in), optional :: absolute                               !< calculate absolute, not relative, difference
        real(dp), intent(in), optional :: r(:)                                  !< normal positions for theta
        
        ! local variable
        real(dp), allocatable :: v_com_loc(:,:,:,:,:)                           ! local copy of v_mag
        
        ! initialize ierr
        ierr = 0
        
        ! set up local copy of v_mag
        allocate(v_com_loc(size(v_mag,1),size(v_mag,2),size(v_mag,3),1,1))
        v_com_loc(:,:,:,1,1) = v_mag
        
        ! call 2-D version
        ierr = calc_tor_diff_2d(v_com_loc,theta,norm_disc_prec,&
            &absolute=absolute,r=r)
        chckerr('')
        
        ! copy
        v_mag = v_com_loc(:,:,:,1,1)
    end function calc_tor_diff_0d
    
    !> Calculates \f$X\f$, \f$Y\f$ and \f$Z\f$ on a grid \c grid_XYZ, determined
    !! through its E(quilibrium) coordinates.
    !!
    !! Furthermore, a  grid \c grid_eq  must be provided,  which is the  grid in
    !! which the  variables concerning \f$R\f$  and \f$Z\f$ are  tabulated, i.e.
    !! the full equilibrium grid in E(quilibrium) coordinates.
    !!
    !! Of this grid, however, only \c r_E  is used, and the rest ignored. It can
    !! therefore be provided without the angular part, i.e. by reconstructing it
    !! with a subset.
    !!
    !! If VMEC is the equilibrium model, this routine also optionally calculates
    !! \f$\lambda\f$ on the grid, as this  is also needed some times. for HELENA
    !! this variable is not used.
    !!
    !! Optionally, \f$R\f$ and \f$\lambda\f$ (only for VMEC) can be returned.
    !!
    !! \note
    !!  -# For VMEC, the trigonometric factors of \c grid_XYZ must be calculated
    !!  beforehand.
    !!  -# The normalization factor \c R_0  for length is taken into account and
    !!  the output is transformed back to unnormalized values.
    !!
    !! \return ierr
    integer function calc_xyz_grid(grid_eq,grid_XYZ,X,Y,Z,L,R) result(ierr)
        use num_vars, only: eq_style, use_normalization
        use num_utilities, only: round_with_tol
        use eq_vars, only: r_0
        
        character(*), parameter :: rout_name = 'calc_XYZ_grid'
        
        ! input / output
        type(grid_type), intent(in) :: grid_eq                                  !< equilibrium grid
        type(grid_type), intent(in) :: grid_xyz                                 !< grid for which to calculate \f$X\f$, \f$Y\f$, \f$Z\f$ and optionally \f$L\f$
        real(dp), intent(inout) :: x(:,:,:)                                     !< \f$X\f$ of grid
        real(dp), intent(inout) :: y(:,:,:)                                     !< \f$Y\f$ of grid
        real(dp), intent(inout) :: z(:,:,:)                                     !< \f$Z\f$ of grid
        real(dp), intent(inout), optional :: l(:,:,:)                           !< \f$\lambda\f$ of grid
        real(dp), intent(inout), optional :: r(:,:,:)                           !< \f$R\f$ of grid
        
        ! local variables
        character(len=max_str_ln) :: err_msg                                    ! error message
        
        ! initialize ierr
        ierr = 0
        
        ! test
        if (size(x,1).ne.grid_xyz%n(1) .or. size(x,2).ne.grid_xyz%n(2) .or. &
            &size(x,3).ne.grid_xyz%loc_n_r) then
            ierr = 1
            err_msg =  'X needs to have the correct dimensions'
            chckerr(err_msg)
        end if
        if (size(y,1).ne.grid_xyz%n(1) .or. size(y,2).ne.grid_xyz%n(2) .or. &
            &size(y,3).ne.grid_xyz%loc_n_r) then
            ierr = 1
            err_msg =  'Y needs to have the correct dimensions'
            chckerr(err_msg)
        end if
        if (size(z,1).ne.grid_xyz%n(1) .or. size(z,2).ne.grid_xyz%n(2) .or. &
            &size(z,3).ne.grid_xyz%loc_n_r) then
            ierr = 1
            err_msg =  'Z needs to have the correct dimensions'
            chckerr(err_msg)
        end if
        if (present(l)) then
            if (size(l,1).ne.grid_xyz%n(1) .or. size(l,2).ne.grid_xyz%n(2) &
                &.or. size(l,3).ne.grid_xyz%loc_n_r) then
                ierr = 1
                err_msg =  'L needs to have the correct dimensions'
                chckerr(err_msg)
            end if
        end if
        
        ! choose which equilibrium style is being used:
        !   1:  VMEC
        !   2:  HELENA
        select case (eq_style)
            case (1)                                                            ! VMEC
                ierr = calc_xyz_grid_vmec(grid_eq,grid_xyz,x,y,z,l=l,r=r)
                chckerr('')
            case (2)                                                            ! HELENA
                ierr = calc_xyz_grid_hel(grid_eq,grid_xyz,x,y,z,r=r)
                chckerr('')
        end select
        
        ! if normalized, return to physical value
        if (use_normalization) then
            x = x*r_0
            y = y*r_0
            z = z*r_0
            if (present(r)) r = r*r_0
        end if
    contains
        ! VMEC version
        !> \private
        integer function calc_xyz_grid_vmec(grid_eq,grid_XYZ,X,Y,Z,L,R) &
            &result(ierr)
            use num_vars, only: norm_disc_prec_eq
            use vmec_utilities, only: fourier2real
            use vmec_vars, only: r_v_c, r_v_s, z_v_c, z_v_s, l_v_c, l_v_s, &
                &mnmax_v, is_asym_v
            use num_utilities, only: spline
            
            character(*), parameter :: rout_name = 'calc_XYZ_grid_VMEC'
            
            ! input / output
            type(grid_type), intent(in) :: grid_eq                              ! equilibrium grid
            type(grid_type), intent(in) :: grid_xyz                             ! grid for which to calculate X, Y, Z and optionally L
            real(dp), intent(inout) :: x(:,:,:), y(:,:,:), z(:,:,:)             ! X, Y and Z of grid
            real(dp), intent(inout), optional :: l(:,:,:)                       ! lambda of grid
            real(dp), intent(inout), optional :: r(:,:,:)                       ! R of grid
            
            ! local variables
            integer :: id                                                       ! counter
            real(dp), allocatable :: r_v_c_int(:,:), r_v_s_int(:,:)             ! interpolated version of R_V_c and R_V_s
            real(dp), allocatable :: z_v_c_int(:,:), z_v_s_int(:,:)             ! interpolated version of Z_V_c and Z_V_s
            real(dp), allocatable :: l_v_c_int(:,:), l_v_s_int(:,:)             ! interpolated version of L_V_c and L_V_s
            real(dp), allocatable :: r_loc(:,:,:)                               ! R in Cylindrical coordinates
            character(len=max_str_ln) :: err_msg                                ! error message
            
            ! initialize ierr
            ierr = 0
            
            ! test
            if (.not.allocated(grid_xyz%trigon_factors)) then
                ierr = 1
                err_msg = 'trigon factors of grid_XYZ need to be allocated'
                chckerr(err_msg)
            end if
            
            ! set up interpolated R_V_c_int, ..
            allocate(r_v_c_int(mnmax_v,grid_xyz%loc_n_r))
            allocate(r_v_s_int(mnmax_v,grid_xyz%loc_n_r))
            allocate(z_v_c_int(mnmax_v,grid_xyz%loc_n_r))
            allocate(z_v_s_int(mnmax_v,grid_xyz%loc_n_r))
            if (present(l)) then
                allocate(l_v_c_int(mnmax_v,grid_xyz%loc_n_r))
                allocate(l_v_s_int(mnmax_v,grid_xyz%loc_n_r))
            end if
            
            ! interpolate VMEC tables
            do id = 1,mnmax_v
                ierr = spline(grid_eq%r_E,r_v_c(id,:,0),grid_xyz%loc_r_E,&
                    &r_v_c_int(id,:),ord=norm_disc_prec_eq)
                chckerr('')
                ierr = spline(grid_eq%r_E,r_v_s(id,:,0),grid_xyz%loc_r_E,&
                    &r_v_s_int(id,:),ord=norm_disc_prec_eq)
                chckerr('')
                ierr = spline(grid_eq%r_E,z_v_c(id,:,0),grid_xyz%loc_r_E,&
                    &z_v_c_int(id,:),ord=norm_disc_prec_eq)
                chckerr('')
                ierr = spline(grid_eq%r_E,z_v_s(id,:,0),grid_xyz%loc_r_E,&
                    &z_v_s_int(id,:),ord=norm_disc_prec_eq)
                chckerr('')
                if (present(l)) then
                    ierr = spline(grid_eq%r_E,l_v_c(id,:,0),grid_xyz%loc_r_E,&
                        &l_v_c_int(id,:),ord=norm_disc_prec_eq)
                    chckerr('')
                    ierr = spline(grid_eq%r_E,l_v_s(id,:,0),grid_xyz%loc_r_E,&
                        &l_v_s_int(id,:),ord=norm_disc_prec_eq)
                    chckerr('')
                end if
            end do
            
            ! allocate local R
            allocate(r_loc(grid_xyz%n(1),grid_xyz%n(2),grid_xyz%loc_n_r))
            
            ! inverse fourier transform with trigonometric factors
            ierr = fourier2real(r_v_c_int,r_v_s_int,grid_xyz%trigon_factors,&
                &r_loc,sym=[.true.,is_asym_v])
            chckerr('')
            ierr = fourier2real(z_v_c_int,z_v_s_int,grid_xyz%trigon_factors,z,&
                &sym=[is_asym_v,.true.])
            chckerr('')
            if (present(l)) then
                ierr = fourier2real(l_v_c_int,l_v_s_int,&
                    &grid_xyz%trigon_factors,l,sym=[is_asym_v,.true.])
                chckerr('')
            end if
            if (present(r)) r = r_loc
            
            ! transform cylindrical to cartesian
            ! (the geometrical zeta is equal to the VMEC zeta)
            x = r_loc*cos(grid_xyz%zeta_E)
            y = r_loc*sin(grid_xyz%zeta_E)
            
            ! deallocate
            deallocate(r_v_c_int,r_v_s_int,z_v_c_int,z_v_s_int)
            if (present(l)) deallocate(l_v_c_int,l_v_s_int)
            deallocate(r_loc)
        end function calc_xyz_grid_vmec
        
        ! HELENA version
        !> \private
        integer function calc_xyz_grid_hel(grid_eq,grid_XYZ,X,Y,Z,R) &
            &result(ierr)
            use helena_vars, only: r_h, z_h, chi_h, ias, nchi
            use num_vars, only: norm_disc_prec_eq
            use num_utilities, only: spline
            
            character(*), parameter :: rout_name = 'calc_XYZ_grid_HEL'
            
            ! input / output
            type(grid_type), intent(in) :: grid_eq                              ! equilibrium grid
            type(grid_type), intent(in) :: grid_xyz                             ! grid for which to calculate X, Y, Z and optionally L
            real(dp), intent(inout) :: x(:,:,:), y(:,:,:), z(:,:,:)             ! X, Y and Z of grid
            real(dp), intent(inout), optional :: r(:,:,:)                       ! R of grid
            
            ! local variables
            integer :: id, jd, kd                                               ! counters
            integer :: pmone                                                    ! plus or minus one
            integer :: bcs(2,3)                                                 ! boundary conditions (theta(even), theta(odd), r)
            real(dp) :: bcs_val(2,3)                                            ! values for boundary conditions
            real(dp), allocatable :: r_h_int(:,:), z_h_int(:,:)                 ! R and Z at interpolated normal value
            real(dp), allocatable :: r_loc(:,:,:)                               ! R in Cylindrical coordinates
            real(dp) :: theta_loc                                               ! local copy of theta_E
            
            ! initialize ierr
            ierr = 0
            
            ! set up boundary conditions
            if (ias.eq.0) then                                                  ! top-bottom symmetric
                bcs(:,1) = [1,1]                                                ! theta(even): zero first derivative
                bcs(:,2) = [2,2]                                                ! theta(odd): zero first derivative
            else
                bcs(:,1) = [-1,-1]                                              ! theta(even): periodic
                bcs(:,2) = [-1,-1]                                              ! theta(odd): periodic
            end if
            bcs(:,3) = [0,0]                                                    ! r: not-a-knot
            bcs_val = 0._dp
            
            ! allocate local R
            allocate(r_loc(grid_xyz%n(1),grid_xyz%n(2),grid_xyz%loc_n_r))
            
            ! set up interpolated R and Z
            allocate(r_h_int(nchi,grid_xyz%loc_n_r))
            allocate(z_h_int(nchi,grid_xyz%loc_n_r))
            
            ! interpolate HELENA output  R_H and Z_H for  every requested normal
            ! point
            do id = 1,nchi
                ierr = spline(grid_eq%r_E,r_h(id,:),grid_xyz%loc_r_E,&
                    &r_h_int(id,:),ord=norm_disc_prec_eq,bcs=bcs(:,3),&
                    &bcs_val=bcs_val(:,3))
                chckerr('')
                ierr = spline(grid_eq%r_E,z_h(id,:),grid_xyz%loc_r_E,&
                    &z_h_int(id,:),ord=norm_disc_prec_eq,bcs=bcs(:,3),&
                    &bcs_val=bcs_val(:,3))
                chckerr('')
            end do
            
            ! Note:  R_H and  Z_H  are not  adapted to  the  parallel grid,  but
            ! tabulated in the original HELENA poloidal grid.
            ! loop over normal points
            do kd = 1,grid_xyz%loc_n_r                                          ! loop over all normal points
                ! loop over toroidal points
                do jd = 1,grid_xyz%n(2)
                    ! interpolate at the requested poloidal points
                    do id = 1,grid_xyz%n(1)
                        ! initialize local theta
                        theta_loc = grid_xyz%theta_E(id,jd,kd)
                        
                        ! add or subtract 2pi to  the parallel angle until it is
                        ! at least 0 to get principal range 0..2pi
                        if (theta_loc.lt.0._dp) then
                            do while (theta_loc.lt.0._dp)
                                theta_loc = theta_loc + 2*pi
                            end do
                        else if (theta_loc.gt.2*pi) then
                            do while (theta_loc.gt.2._dp*pi)
                                theta_loc = theta_loc - 2*pi
                            end do
                        end if
                        
                        ! take into account possible symmetry
                        if (ias.eq.0 .and. theta_loc.gt.pi) then
                            theta_loc = 2*pi-theta_loc
                            pmone = -1
                        else
                            pmone = 1
                        end if
                        
                        ! interpolate  the HELENA  variables poloidally
                        ierr = spline(chi_h,r_h_int(:,kd),[theta_loc],&
                            &r_loc(id:id,jd,kd),ord=norm_disc_prec_eq,&
                            &bcs=bcs(:,1),bcs_val=bcs_val(:,1))                 ! even
                        chckerr('')
                        ierr = spline(chi_h,pmone*z_h_int(:,kd),[theta_loc],&
                            &z(id:id,jd,kd),ord=norm_disc_prec_eq,&
                            &bcs=bcs(:,2),bcs_val=bcs_val(:,2))                 ! odd
                        chckerr('')
                    end do
                end do
            end do
            if (present(r)) r = r_loc
            
            ! calculate X and Y, transforming cylindrical to cartesian
            ! (the geometrical zeta is the inverse of HELENA zeta)
            x = r_loc*cos(-grid_xyz%zeta_E)
            y = r_loc*sin(-grid_xyz%zeta_E)
            
            ! deallocate
            deallocate(r_loc)
        end function calc_xyz_grid_hel
    end function calc_xyz_grid
    
    !> Extend a grid angularly.
    !!
    !! This  is done  using  equidistant  variables of  \c  n_theta_plot and  \c
    !! n_zeta_plot angular and own \c loc_n_r points in F coordinates.
    !!
    !! Optionally,  the grid  can  also be  converted to  E  coordinates if  the
    !! equilibrium  grid is  provided.  This is  required  for many  operations,
    !! such  as  the  calculation  of   \f$X\f$,  \f$Y\f$  and  \f$Z\f$  through
    !! calc_XYZ_grid().
    !!
    !! \note For VMEC, it can be slow, as grid_utilities.coord_f2e() is used.
    !!
    !! \return ierr
    integer function extend_grid_f(grid_in,grid_ext,grid_eq,n_theta_plot,&
        &n_zeta_plot,lim_theta_plot,lim_zeta_plot) result(ierr)
        use num_vars, only: n_theta => n_theta_plot, n_zeta => n_zeta_plot, &
            &min_theta_plot, max_theta_plot, min_zeta_plot, max_zeta_plot
        use mpi_utilities, only: get_ser_var
        
        character(*), parameter :: rout_name = 'extend_grid_F'
        
        ! input / output
        type(grid_type), intent(in) :: grid_in                                  !< grid to be extended
        type(grid_type), intent(inout) :: grid_ext                              !< extended grid
        type(grid_type), intent(in), optional :: grid_eq                        !< equilibrium grid
        integer, intent(in), optional :: n_theta_plot                           !< number of poins in theta direction
        integer, intent(in), optional :: n_zeta_plot                            !< number of poins in zeta direction
        real(dp), intent(in), optional :: lim_theta_plot(2)                     !< limits in theta
        real(dp), intent(in), optional :: lim_zeta_plot(2)                      !< limits in zeta
        
        ! local variables
        type(grid_type) :: grid_ext_trim                                        ! trimmed grid_ext
        integer :: n_theta_plot_loc                                             ! local n_theta_plot
        integer :: n_zeta_plot_loc                                              ! local n_zeta_plot
        real(dp) :: lim_theta_plot_loc(2)                                       ! local limits of theta_plot
        real(dp) :: lim_zeta_plot_loc(2)                                        ! local limits of zeta_plot
        
        ! initialize ierr
        ierr = 0
        
        ! set up local variables
        n_theta_plot_loc = n_theta
        if (present(n_theta_plot)) n_theta_plot_loc = n_theta_plot
        n_zeta_plot_loc = n_zeta
        if (present(n_zeta_plot)) n_zeta_plot_loc = n_zeta_plot
        lim_theta_plot_loc = [min_theta_plot,max_theta_plot]
        if (present(lim_theta_plot)) lim_theta_plot_loc = lim_theta_plot
        lim_zeta_plot_loc = [min_zeta_plot,max_zeta_plot]
        if (present(lim_zeta_plot)) lim_zeta_plot_loc = lim_zeta_plot
        
        ! user ouput
        call writo('Theta limits: '//trim(r2strt(lim_theta_plot_loc(1)))//&
            &'pi .. '//trim(r2strt(lim_theta_plot_loc(2)))//'pi')
        call writo('Zeta limits: '//trim(r2strt(lim_zeta_plot_loc(1)))//&
            &'pi .. '//trim(r2strt(lim_zeta_plot_loc(2)))//'pi')
        
        ! creating  equilibrium  grid  for  the output  that  covers  the  whole
        ! geometry angularly in F coordinates
        ierr = grid_ext%init([n_theta_plot_loc,n_zeta_plot_loc,grid_in%n(3)],&
            &[grid_in%i_min,grid_in%i_max])
        chckerr('')
        grid_ext%r_F = grid_in%r_F
        grid_ext%loc_r_F = grid_in%loc_r_F
        ierr = calc_eqd_grid(grid_ext%theta_F,lim_theta_plot_loc(1)*pi,&
            &lim_theta_plot_loc(2)*pi,1)
        chckerr('')
        ierr = calc_eqd_grid(grid_ext%zeta_F,lim_zeta_plot_loc(1)*pi,&
            &lim_zeta_plot_loc(2)*pi,2)
        chckerr('')
        
        ! convert all F coordinates to E coordinates if requested
        if (present(grid_eq)) then
            call writo('convert F to E coordinates')
            call lvl_ud(1)
            
            ! set total r_E
            grid_ext%r_E = grid_in%r_E
            
            ! get local r_E and angular theta_E and zeta_E
            ierr = coord_f2e(grid_eq,&
                &grid_ext%loc_r_F,grid_ext%theta_F,grid_ext%zeta_F,&
                &grid_ext%loc_r_E,grid_ext%theta_E,grid_ext%zeta_E)
            chckerr('')
            
            ! trim external grid
            ierr = trim_grid(grid_ext,grid_ext_trim)
            chckerr('')
            
            ! clean up
            call grid_ext_trim%dealloc()
            
            call lvl_ud(-1)
        end if
    end function extend_grid_f
    
    !> Copy a grid A to a new grid B, that was not yet initialized.
    !!
    !! This  new grid  can contain  a  subsection of  the previous  grid in  all
    !! dimensions. It can also be a divided grid, by providing the limits
    !!
    !! \note
    !!  -# The normal  limits for the divided  grid should be given  in terms of
    !!  the normal dimension of \c grid B.
    !!  -# if the grids are not  purely normal, the procedure can currently only
    !!  handle the copying of a \c grid_A that is not divided.
    !!
    !! \return ierr
    integer function copy_grid(grid_A,grid_B,lims_B,i_lim) result(ierr)
        character(*), parameter :: rout_name = 'copy_grid'
        
        ! input / output
        class(grid_type), intent(in) :: grid_a                                  !< grid to be initialized
        class(grid_type), intent(inout) :: grid_b                               !< grid to be initialized
        integer, intent(in), optional :: lims_b(3,2)                            !< ranges for subset of grid
        integer, intent(in), optional :: i_lim(2)                               !< min. and max. local normal index
        
        ! local variables
        integer :: n_b(3)                                                       ! n of grid B
        integer :: lims_b_loc(3,2)                                              ! local lims_B
        integer :: i_lim_loc(2)                                                 ! local i_lim
        
        ! initialize ierr
        ierr = 0
        
        ! set up dimensions and of B
        lims_b_loc(1,:) = [1,grid_a%n(1)]
        lims_b_loc(2,:) = [1,grid_a%n(2)]
        lims_b_loc(3,:) = [1,grid_a%n(3)]
        if (present(lims_b)) lims_b_loc = lims_b
        n_b = lims_b_loc(:,2)-lims_b_loc(:,1)+1
        where (lims_b_loc(:,1).eq.lims_b_loc(:,2) &
            &.and. lims_b_loc(:,1).eq.[0,0,0]) n_b = 0
        i_lim_loc = lims_b_loc(3,:)
        if (present(i_lim)) i_lim_loc = lims_b_loc(3,1) - 1 + i_lim
        ! tests
        if (n_b(1).gt.grid_a%n(1) .or. n_b(2).gt.grid_a%n(2) .or. &
            &n_b(3).gt.grid_a%n(3)) then
            write(*,*) 'n_B' ,n_b
            write(*,*) 'grid_A', grid_a%n
            ierr = 1
            chckerr('lims_B is too large')
        end if
        if (i_lim_loc(1).gt.grid_a%n(3) .or. i_lim_loc(2).gt.grid_a%n(3)) then
            ierr = 1
            chckerr('normal limits too wide')
        end if
        
        ! allocate the new grid
        ierr = grid_b%init(n_b,i_lim)
        chckerr('')
        
        ! copy arrays, possibly subset
        grid_b%r_F = grid_a%r_F(lims_b_loc(3,1):lims_b_loc(3,2))
        grid_b%r_E = grid_a%r_E(lims_b_loc(3,1):lims_b_loc(3,2))
        grid_b%loc_r_F = grid_a%r_F(i_lim_loc(1):i_lim_loc(2))
        grid_b%loc_r_E = grid_a%r_E(i_lim_loc(1):i_lim_loc(2))
        if (n_b(1).ne.0 .and. n_b(2).ne.0) then
            if (grid_a%divided) then
                ierr = 1
                chckerr('grid_A cannot be divided')
            end if
            grid_b%theta_F = grid_a%theta_F(lims_b_loc(1,1):lims_b_loc(1,2),&
                &lims_b_loc(2,1):lims_b_loc(2,2),i_lim_loc(1):i_lim_loc(2))
            grid_b%theta_E = grid_a%theta_E(lims_b_loc(1,1):lims_b_loc(1,2),&
                &lims_b_loc(2,1):lims_b_loc(2,2),i_lim_loc(1):i_lim_loc(2))
            grid_b%zeta_F = grid_a%zeta_F(lims_b_loc(1,1):lims_b_loc(1,2),&
                &lims_b_loc(2,1):lims_b_loc(2,2),i_lim_loc(1):i_lim_loc(2))
            grid_b%zeta_E = grid_a%zeta_E(lims_b_loc(1,1):lims_b_loc(1,2),&
                &lims_b_loc(2,1):lims_b_loc(2,2),i_lim_loc(1):i_lim_loc(2))
        end if
    end function copy_grid
    
    !> Calculates volume integral on a 3D grid.
    !!
    !! Two angular and  one normal variable have  to be provided on a the grid.
    !!
    !! If the <tt>i</tt>'th dimension of the  grid is equal to one, the function
    !! to be integrated is assumed not to vary in this dimension.
    !!
    !! Furthermore, if \c i is 1 or 2, the corresponding <tt>i</tt>'th (angular)
    !! variable is the only variable that  is assumed to vary in that dimension.
    !! The other angular variable as well  as the normal variable are assumed to
    !! be constant like the function itself,  resulting in a factor \f$2 \pi\f$.
    !! However, if \c i is 3, an error is displayed as this does not represent a
    !! physical situation.
    !!
    !! A common case  through which to understand this is  the axisymmetric case
    !! where the first angular variable  \f$\theta\f$ varies in the dimensions 1
    !! and  3,  the second  angular  variable  \f$\zeta\f$  varies only  in  the
    !! dimension 2, and the normal variable only varies in the dimension 3.
    !!
    !! Alternatively, there  is the  case of a  grid-aligned set  of coordinates
    !! \f$\theta\f$ and  \f$\alpha\f$, where the first  dimension corresponds to
    !! the direction along the magnetic field line, the second to the geodesical
    !! direction and the third to the  normal direction. If the calculations for
    !! different  field  lines  are  decoupled,  the  variation  in  the  second
    !! dimension is not taken into account and no integration happens along it.
    !!
    !! Internally, the angular variables and  the normal variable are related to
    !! the  coordinates \f$\left(x,y,z\right)\f$  that correspond  to the  three
    !! dimensions. They thus  form a computational orthogonal grid  to which the
    !! original coordinates are related through the transformation of Jacobians:
    !!  \f[\mathcal{J}_{xyz} = \mathcal{J}_\text{F}
    !!      \frac{\partial r_\text{F}}{\partial z}
    !!      \left(\frac{\partial \gamma_1}{\partial x}
    !!      \frac{\gamma_2}{\partial y}-
    !!      \frac{\partial \gamma_1}{\partial y}
    !!      \frac{\gamma_2}{\partial x}\right) , \f]
    !! where \f$\gamma_1\f$  and \f$\gamma_2\f$  are \c ang_1  and \c  ang_2, so
    !! that the integral becomes
    !!  \f[ \sum_{x,y,z} \left[ f\left(x,y,z\right) \mathcal{J}_\text{F}
    !!      \frac{\partial r_\text{F}}{\partial z}
    !!      \left(\frac{\partial \gamma_1}{\partial x}
    !!      \frac{\gamma_2}{\partial y}-
    !!      \frac{\partial \gamma_1}{\partial y}
    !!      \frac{\gamma_2}{\partial x}\right)
    !!      \Delta_x \Delta_y \Delta_z \right] \f]
    !! where \f$\Delta_x\f$, \f$\Delta_y\f$ and \f$\Delta_z\f$ are all trivially
    !! equal to 1.
    !!
    !! The integrand  has to be  evaluated at the intermediate  positions inside
    !! the cells. This  is done by taking the average  of the \f$2^3=8\f$ points
    !! for  \c fj  = \f$f  J_\text{F}\f$ as  well as  the transformation  of the
    !! Jacobian to \f$\left(x,y,z\right)\f$ coordinates.
    !!
    !! \see See grid_vars.grid_type for a discussion on \c ang_1 and \c ang_2.
    !!
    !! \note
    !!  -# if the coordinates are independent,  this method is equivalent to the
    !!  repeated  numerical integration  using  the trapezoidal  method, \b  NOT
    !!  Simpson's 3/8 rule!
    !!  -# by setting \c debug_calc_int_vol, this  method can be compared to the
    !!  trapezoidal and simple method for  independent coordinates, again \b NOT
    !!  for Simpson's 3/8 rule!
    !!  -# The  Simpson's 3/8  rule could be  developed but it  is not  of great
    !!  importance.
    !!
    !! \return ierr
    integer function calc_int_vol(ang_1,ang_2,norm,J,f,f_int) result(ierr)
#if ldebug
        use num_vars, only: rank, n_procs
#endif
        
        character(*), parameter :: rout_name = 'calc_int_vol'
        
        ! input / output
        real(dp), intent(in) :: ang_1(:,:,:)                                    !< coordinate variable \f$\gamma_1\f$
        real(dp), intent(in) :: ang_2(:,:,:)                                    !< coordinate variable \f$\gamma_2\f$
        real(dp), intent(in) :: norm(:)                                         !< coordinate variable \f$r_\text{F}\f$
        real(dp), intent(in) :: j(:,:,:)                                        !< Jacobian
        complex(dp), intent(in) :: f(:,:,:,:)                                   !< input f(n_par,n_geo,n_r,size_X^2)
        complex(dp), intent(inout) :: f_int(:)                                  !< output integrated f
        
        ! local variables
        integer :: id, jd, kd, ld                                               ! counters
        integer :: nn_mod                                                       ! number of indices for V and V_int
        integer :: k_min                                                        ! minimum k
        integer :: dims(3)                                                      ! real dimensions
        real(dp), allocatable :: transf_j(:,:,:,:)                              ! comps. of transf. between J and Jxyz: theta_x, zeta_y, theta_y, zeta_x and r_F_z
        real(dp), allocatable :: transf_j_tot(:,:,:)                            ! transf. between J and Jxyz: (theta_x*zeta_y-theta_y*zeta_x)*r_F_z
        complex(dp), allocatable :: jf(:,:,:,:)                                 ! J*f
        character(len=max_str_ln) :: err_msg                                    ! error message
        integer :: i_a(3), i_b(3)                                               ! limits of building blocks of intermediary variables
        logical :: dim_1(2)                                                     ! whether the dimension is equal to one
#if ldebug
        character(len=max_str_ln), allocatable :: var_names(:,:)                ! names of variables
        complex(dp), allocatable :: f_int_alt(:)                                ! alternative calculation for output
        complex(dp), allocatable :: f_int_alt_1d(:,:)                           ! intermediary step in f_int_ALT
        complex(dp), allocatable :: f_int_alt_2d(:,:,:)                         ! intermediary step in f_int_ALT
        complex(dp), allocatable :: f_int_alt_alt(:)                            ! another alternative calculation for output
        integer :: loc_norm(2)                                                  ! local normal index
#endif
        
        ! initialize ierr
        ierr = 0
        
        ! set nn_mod
        nn_mod = size(f,4)
        
        ! tests
        if (size(f_int).ne.nn_mod) then
            ierr = 1
            err_msg = 'f and f_int need to have the same storage convention'
            chckerr(err_msg)
        end if
        if (size(ang_1).ne.size(ang_2) .or. size(ang_1,3).ne.size(norm)) then
            ierr = 1
            err_msg = 'The angular variables have to be compatible'
            chckerr(err_msg)
        end if
        if (size(norm).eq.1) then
            ierr = 1
            err_msg = 'The normal grid size has to be greater than one'
            chckerr(err_msg)
        end if
        
        ! set up dims
        dims = shape(ang_1)
        
        ! set up dim_1
        dim_1 = dims(1:2).eq.1
        
        ! set up Jf and transf_J
        allocate(jf(max(dims(1)-1,1),max(dims(2)-1,1),dims(3)-1,nn_mod))
        allocate(transf_j(max(dims(1)-1,1),max(dims(2)-1,1),dims(3)-1,5))
        allocate(transf_j_tot(max(dims(1)-1,1),max(dims(2)-1,1),dims(3)-1))
        
        ! set up k_min
        k_min = 1
        
        ! get intermediate Jf and components of transf_J:
        !   1: ang_1_x
        !   2: ang_2_y
        !   3: ang_1_y
        !   4: ang_2_x
        !   5: r_F_z
        ! Note: the are always 8 terms in the summation, hence the factor 0.125
        jf = 0._dp
        transf_j = 0._dp
        do kd = 1,2
            if (kd.eq.1) then
                i_a(3) = 1
                i_b(3) = dims(3)-1
            else
                i_a(3) = 2
                i_b(3) = dims(3)
            end if
            do jd = 1,2
                if (dim_1(2)) then
                    i_a(2) = 1
                    i_b(2) = dims(2)
                else
                    if (jd.eq.1) then
                        i_a(2) = 1
                        i_b(2) = dims(2)-1
                    else
                        i_a(2) = 2
                        i_b(2) = dims(2)
                    end if
                end if
                do id = 1,2
                    if (dim_1(1)) then
                        i_a(1) = 1
                        i_b(1) = dims(1)
                    else
                        if (id.eq.1) then
                            i_a(1) = 1
                            i_b(1) = dims(1)-1
                        else
                            i_a(1) = 2
                            i_b(1) = dims(1)
                        end if
                    end if
                    do ld = 1,nn_mod
                        ! Jf
                        jf(:,:,:,ld) = jf(:,:,:,ld) + 0.125_dp * &
                            &j(i_a(1):i_b(1),i_a(2):i_b(2),i_a(3):i_b(3)) * &
                            &f(i_a(1):i_b(1),i_a(2):i_b(2),i_a(3):i_b(3),ld)
                    end do
                    if (dim_1(1)) then
                        ! transf_J(1): ang_1_x
                        transf_j(:,:,:,1) = 2*pi
                        ! transf_J(4): ang_2_x
                        transf_j(:,:,:,4) = 0._dp
                    else
                        do ld = 1,dims(1)-1
                            ! transf_J(1): ang_1_x
                            transf_j(ld,:,:,1) = transf_j(ld,:,:,1) + &
                                &0.125_dp * ( &
                                &ang_1(ld+1,i_a(2):i_b(2),i_a(3):i_b(3)) - &
                                &ang_1(ld,i_a(2):i_b(2),i_a(3):i_b(3)) )
                            ! transf_J(4): ang_2_x
                            transf_j(ld,:,:,4) = transf_j(ld,:,:,4) + &
                                &0.125_dp * ( &
                                &ang_2(ld+1,i_a(2):i_b(2),i_a(3):i_b(3)) - &
                                &ang_2(ld,i_a(2):i_b(2),i_a(3):i_b(3)) )
                        end do
                    end if
                    if (dim_1(2)) then
                        ! transf_J(2): ang_2_y
                        transf_j(:,:,:,2) = 2*pi
                        ! transf_J(3): ang_1_y
                        transf_j(:,:,:,3) = 0._dp
                    else
                        do ld = 1,dims(2)-1
                            ! transf_J(2): ang_2_y
                            transf_j(:,ld,:,2) = transf_j(:,ld,:,2) + &
                                &0.125_dp * ( &
                                &ang_2(i_a(1):i_b(1),ld+1,i_a(3):i_b(3)) - &
                                &ang_2(i_a(1):i_b(1),ld,i_a(3):i_b(3)) )
                            ! transf_J(3): ang_1_y
                            transf_j(:,ld,:,3) = transf_j(:,ld,:,3) + &
                                &0.125_dp * ( &
                                &ang_1(i_a(1):i_b(1),ld+1,i_a(3):i_b(3)) - &
                                &ang_1(i_a(1):i_b(1),ld,i_a(3):i_b(3)) )
                        end do
                    end if
                    do ld = 1,dims(3)-1
                        ! transf_J(5): r_F_z
                        transf_j(:,:,ld,5) = transf_j(:,:,ld,5) + 0.125_dp * ( &
                            &norm(ld+1) - norm(ld) )
                    end do
                end do
            end do
        end do
        
        ! set up total transf_J
        transf_j_tot = (transf_j(:,:,:,1)*transf_j(:,:,:,2)-&
            &transf_j(:,:,:,3)*transf_j(:,:,:,4))*transf_j(:,:,:,5)
        
        ! integrate term  over ang_par_F for all equilibrium  grid points, using
        ! dx = dy = dz = 1
        do ld = 1,nn_mod
            f_int(ld) = sum(jf(:,:,:,ld)*transf_j_tot)
        end do
        
#if ldebug
        if (debug_calc_int_vol) then
            call writo('Testing whether volume integral is calculated &
                &correctly')
            call lvl_ud(1)
            
            allocate(var_names(nn_mod,2))
            var_names(:,1) = 'real part of integrand J f'
            var_names(:,2) = 'imaginary part of integrand J f'
            do ld = 1,nn_mod
                var_names(ld,1) = trim(var_names(ld,1))//' '//trim(i2str(ld))
                var_names(ld,2) = trim(var_names(ld,2))//' '//trim(i2str(ld))
            end do
            
            call plot_hdf5(var_names(:,1),'TEST_RE_Jf_'//trim(i2str(rank)),&
                &rp(jf))
            call plot_hdf5(var_names(:,2),'TEST_IM_Jf_'//trim(i2str(rank)),&
                &ip(jf))
            call plot_hdf5('transformation of Jacobians',&
                &'TEST_transf_J_'//trim(i2str(rank)),transf_j_tot)
            
            ! do alternative calculations
            ! assuming that coordinate i varies only in dimensions i!!!
            allocate(f_int_alt(nn_mod))
            allocate(f_int_alt_1d(size(f,3),nn_mod))
            allocate(f_int_alt_2d(size(f,2),size(f,3),nn_mod))
            allocate(f_int_alt_alt(nn_mod))
            f_int_alt = 0._dp
            f_int_alt_1d = 0._dp
            f_int_alt_2d = 0._dp
            f_int_alt_alt = 0._dp
            
            ! integrate in first coordinate
            do kd = 1,size(f,3)
                do jd = 1,size(f,2)
                    if (size(f,1).gt.1) then
                        do id = 2,size(f,1)
                            f_int_alt_2d(jd,kd,:) = f_int_alt_2d(jd,kd,:) + &
                                &0.5*(j(id,jd,kd)*f(id,jd,kd,:)+&
                                &j(id-1,jd,kd)*f(id-1,jd,kd,:))*&
                                &(ang_1(id,jd,kd)-ang_1(id-1,jd,kd))
                        end do
                    else
                        f_int_alt_2d(jd,kd,:) = j(1,jd,kd)*f(1,jd,kd,:)
                    end if
                end do
            end do
            
            ! integrate in second coordinate
            do kd = 1,size(f,3)
                if (size(f,2).gt.1) then
                    do jd = 2,size(f,2)
                        f_int_alt_1d(kd,:) = f_int_alt_1d(kd,:) + &
                            &0.5*(f_int_alt_2d(jd,kd,:)+&
                            &f_int_alt_2d(jd-1,kd,:))*&
                            &(ang_2(1,jd,kd)-ang_2(1,jd-1,kd))                  ! assuming that ang_2 is not dependent on dimension 1
                    end do
                else
                    f_int_alt_1d(kd,:) = f_int_alt_2d(1,kd,:)
                end if
            end do
            
            ! integrate in third coordinate
            if (size(f,3).gt.1) then
                do kd = 2,size(f,3)
                    f_int_alt(:) = f_int_alt(:) + &
                        &0.5*(f_int_alt_1d(kd,:)+f_int_alt_1d(kd-1,:))*&
                        &(norm(kd)-norm(kd-1))
                end do
            else
                f_int_alt = f_int_alt_1d(1,:)
            end if
            
            ! second alternative, first order integration
            loc_norm = [1,size(f,3)]
            if (rank.lt.n_procs-1) loc_norm(2) = loc_norm(2)-1                  ! no ghost region needed for this method
            do ld = 1,size(f,4)
                f_int_alt_alt(ld) = sum(j(:,:,loc_norm(1):loc_norm(2))*&
                    &f(:,:,loc_norm(1):loc_norm(2),ld))
            end do
            if (size(f,1).gt.1) f_int_alt_alt = f_int_alt_alt*&
                &(ang_1(2,1,1)-ang_1(1,1,1))
            if (size(f,2).gt.1) f_int_alt_alt = f_int_alt_alt*&
                &(ang_2(1,2,1)-ang_2(1,1,1))
            if (size(f,3).gt.1) f_int_alt_alt = f_int_alt_alt*&
                &(norm(2)-norm(1))
            
            call writo('Resulting integral: ',persistent=.true.)
            call lvl_ud(1)
            do ld = 1,nn_mod
                call writo('rank '//trim(i2str(rank))//' integral      '//&
                    &trim(i2str(ld))//': '//trim(c2str(f_int(ld))),&
                    &persistent=.true.)
                call writo('rank '//trim(i2str(rank))//'   alternative '//&
                    &trim(i2str(ld))//': '//trim(c2str(f_int_alt(ld))),&
                    &persistent=.true.)
                call writo('rank '//trim(i2str(rank))//'   other alt.  '//&
                    &trim(i2str(ld))//': '//trim(c2str(f_int_alt_alt(ld))),&
                    &persistent=.true.)
            end do
            call lvl_ud(-1)
            
            call writo('(Note that alternative calculations are only valid if &
                &independent coordinates!)')
            
            call lvl_ud(-1)
        end if
#endif
        
        ! deallocate local variables
        deallocate(jf,transf_j,transf_j_tot)
    end function calc_int_vol
    
    !> Trim a grid, removing any overlap between the different regions.
    !! 
    !! by default, the routine assumes a symmetric ghost region and cuts as many
    !! grid points from the end of the previous process as from the beginning of
    !! the next  process, but if the  number of overlapping grid  points is odd,
    !! the next process looses one more point.
    !!
    !! optionally, the trimmed  indices in the normal direction  can be provided
    !! in  \c  norm_id,  i.e.  the  indices in  the  old,  untrimmed  grid  that
    !! correspond to the start and end indices of the trimmed grid. E.g. if
    !!  - proc 0:  3 ... 25
    !!  - proc 1: 20 ... 50
    !!
    !! then the trimmed grid will be:
    !!  - proc 0:  3 ... 22
    !!  - proc 1: 23 ... 50
    !!
    !! which is shifted down by 2 to 
    !!  - proc 0:  1 ... 20
    !!  - proc 1: 21 ... 48
    !!
    !! in the trimmed  grid. The indices of the  previous step (3 & 22  and 23 &
    !! 50) are saved in \c norm_id.
    !!
    !! \return ierr
    integer function trim_grid(grid_in,grid_out,norm_id) result(ierr)
        use num_vars, only: n_procs, rank
        use mpi_utilities, only: get_ser_var
        
        character(*), parameter :: rout_name = 'trim_grid'
        
        ! input / output
        type(grid_type), intent(in) :: grid_in                                  !< input grid
        type(grid_type), intent(inout) :: grid_out                              !< trimmed grid
        integer, intent(inout), optional :: norm_id(2)                          !< normal indices corresponding to trimmed part
        
        ! local variables
        integer, allocatable :: tot_i_min(:)                                    ! i_min of grid of all processes
        integer, allocatable :: tot_i_max(:)                                    ! i_max of grid of all processes
        integer :: i_lim_out(2)                                                 ! i_lim of output grid
        integer :: n_out(3)                                                     ! n of output grid
        
        ! initialize ierr
        ierr = 0
        
        ! detect whether grid divided
        if (grid_in%divided) then
            ! get min_i's of the grid_in
            ierr = get_ser_var([grid_in%i_min],tot_i_min,scatter=.true.)
            chckerr('')
            
            ! get max_i's of the grid_in
            ierr = get_ser_var([grid_in%i_max],tot_i_max,scatter=.true.)
            chckerr('')
            
            ! set  i_lim of trimmed output  grid (not yet shifted  by first proc
            ! min)
            if (rank.gt.0) then
                i_lim_out(1) = max(tot_i_min(1),tot_i_min(rank+1)+&
                    &ceiling((tot_i_max(rank)-tot_i_min(rank+1)+1._dp)/2))
            else
                i_lim_out(1) = tot_i_min(1)
            end if
            if (rank.lt.n_procs-1) then
                i_lim_out(2) = min(tot_i_max(n_procs),tot_i_max(rank+1)-&
                    &floor((tot_i_max(rank+1)-tot_i_min(rank+2)+1._dp)/2))
            else
                i_lim_out(2) = tot_i_max(n_procs)
            end if
            
            ! limit to input limits (necessary if a process has only one point)
            i_lim_out(1) = max(min(i_lim_out(1),grid_in%i_max),grid_in%i_min)
            i_lim_out(2) = max(min(i_lim_out(2),grid_in%i_max),grid_in%i_min)
            
            ! get  min_i's and max_i's  of the grid_out,  not shifted by  min of
            ! first process
            ierr = get_ser_var([i_lim_out(1)],tot_i_min,scatter=.true.)
            chckerr('')
            ierr = get_ser_var([i_lim_out(2)],tot_i_max,scatter=.true.)
            chckerr('')
            
            ! set n of output grid
            n_out(1) = grid_in%n(1)
            n_out(2) = grid_in%n(2)
            n_out(3) = tot_i_max(n_procs)-tot_i_min(1)+1
            
            ! create new grid
            ierr = grid_out%init(n_out,i_lim_out-tot_i_min(1)+1,divided=.true.) ! limits shifted by min of first process
            chckerr('')
            
            ! recycle  i_lim_out  for  shifted  array  indices, set  norm_id  if
            ! requested
            i_lim_out = i_lim_out - grid_in%i_min + 1
            if (present(norm_id)) norm_id = i_lim_out
        else
            ! set n of output grid
            n_out = grid_in%n
            
            ! set unshifted i_lim_out and norm_id if requested
            i_lim_out = [grid_in%i_min,grid_in%i_max]
            if (present(norm_id)) norm_id = i_lim_out
            
            ! shift i_lim_out by min.
            i_lim_out = i_lim_out-i_lim_out(1)+1
            
            ! create new grid
            ierr = grid_out%init(n_out,i_lim_out)                               ! grid not divided
            chckerr('')
        end if
        
        ! copy local arrays
        if (grid_in%n(1).ne.0 .and. grid_in%n(2).ne.0) then                     ! only if 3D grid
            grid_out%theta_E = grid_in%theta_E(:,:,i_lim_out(1):i_lim_out(2))
            grid_out%zeta_E = grid_in%zeta_E(:,:,i_lim_out(1):i_lim_out(2))
            grid_out%theta_F = grid_in%theta_F(:,:,i_lim_out(1):i_lim_out(2))
            grid_out%zeta_F = grid_in%zeta_F(:,:,i_lim_out(1):i_lim_out(2))
        end if
        if (grid_in%divided) then                                               ! but if input grid divided, loc_r gets priority
            grid_out%loc_r_E = grid_in%loc_r_E(i_lim_out(1):i_lim_out(2))
            grid_out%loc_r_F = grid_in%loc_r_F(i_lim_out(1):i_lim_out(2))
        end if
        
        ! if divided, set total arrays
        if (grid_in%divided) then
            grid_out%r_E = grid_in%r_E(tot_i_min(1):tot_i_max(n_procs))
            grid_out%r_F = grid_in%r_F(tot_i_min(1):tot_i_max(n_procs))
        else
            grid_out%r_E = grid_in%r_E
            grid_out%r_F = grid_in%r_F
        end if
        
        ! if allocated, copy trigonometric factors
        if (allocated(grid_in%trigon_factors)) then
            allocate(grid_out%trigon_factors(&
                &size(grid_in%trigon_factors,1),&
                &size(grid_in%trigon_factors,2),&
                &size(grid_in%trigon_factors,3),&
                &grid_out%loc_n_r,2))
            grid_out%trigon_factors = &
                &grid_in%trigon_factors(:,:,:,i_lim_out(1):i_lim_out(2),:)
        end if
    end function trim_grid
    
    !> Untrims a trimmed grid by  introducing an asymmetric ghost regions at the
    !! right.
    !!
    !! The width of the ghost region has to be provided.
    !!
    !! \note
    !!  -# the ghosted grid should be deallocated (with dealloc_grid()).
    !!  -# the input grid has to be trimmed.
    !!
    !! \return ierr
    integer function untrim_grid(grid_in,grid_out,size_ghost) result(ierr)
        use num_vars, only: n_procs, rank
        use mpi_utilities, only: get_ghost_arr, get_ser_var
        
        character(*), parameter :: rout_name = 'untrim_grid'
        
        ! input / output
        type(grid_type), intent(in) :: grid_in                                  !< input grid
        type(grid_type), intent(inout) :: grid_out                              !< ghosted grid
        integer, intent(in) :: size_ghost                                       !< width of ghost region
        
        ! local variables
        integer :: i_lim_in(2)                                                  ! limits of input grid
        integer :: i_lim_out(2)                                                 ! limits of ghosted grid
        integer, allocatable :: tot_loc_n_r(:)                                  ! loc_n_r of all processes
        integer :: size_ghost_loc                                               ! local size_ghost
        
        ! initialize ierr
        ierr = 0
        
        ! get array sizes of all processes
        ierr = get_ser_var([grid_in%loc_n_r],tot_loc_n_r,scatter=.true.)
        chckerr('')
        
        ! set local size_ghost
        size_ghost_loc = min(size_ghost,minval(tot_loc_n_r)-1)
        
        ! set i_lim_in and i_lim_out
        i_lim_in = [grid_in%i_min,grid_in%i_max]
        i_lim_out = i_lim_in
        if (rank.lt.n_procs-1) i_lim_out(2) = i_lim_out(2)+size_ghost_loc
        
        ! create grid
        ierr = grid_out%init(grid_in%n,i_lim_out)
        chckerr('')
        
        ! set ghosted variables
        if (grid_in%n(1).ne.0 .and. grid_in%n(2).ne.0) then                     ! only if 3d grid
            grid_out%theta_e(:,:,1:grid_in%loc_n_r) = grid_in%theta_e
            grid_out%zeta_e(:,:,1:grid_in%loc_n_r) = grid_in%zeta_e
            grid_out%theta_f(:,:,1:grid_in%loc_n_r) = grid_in%theta_f
            grid_out%zeta_f(:,:,1:grid_in%loc_n_r) = grid_in%zeta_f
            ierr = get_ghost_arr(grid_out%theta_e,size_ghost_loc)
            chckerr('')
            ierr = get_ghost_arr(grid_out%zeta_e,size_ghost_loc)
            chckerr('')
            ierr = get_ghost_arr(grid_out%theta_f,size_ghost_loc)
            chckerr('')
            ierr = get_ghost_arr(grid_out%zeta_f,size_ghost_loc)
            chckerr('')
        end if
        if (grid_in%divided) then                                               ! but if input grid divided, loc_r gets priority
            grid_out%loc_r_e = grid_in%r_e(i_lim_out(1):i_lim_out(2))
            grid_out%loc_r_f = grid_in%r_f(i_lim_out(1):i_lim_out(2))
        end if
        grid_out%r_e = grid_in%r_e
        grid_out%r_f = grid_in%r_f
    end function untrim_grid
    
    !> Calculates  contra- and  covariant  components of  a  vector in  multiple
    !! coordinate systems.
    !!
    !! Chain of coordinate systems considered:
    !!  -# Flux F: \f$(\alpha,\psi,\theta)\f$,
    !!  -# Magnetic M: \f$(\psi,\theta,\zeta)\f$,
    !!      - H: \f$(\psi,\theta,\phi)\f$ if HELENA is used,
    !!      - V: \f$(r,\theta,\zeta)\f$ if VMEC is used,
    !!  -# Cylindrical C: \f$(R,\varphi,Z)\f$,
    !!  -# Cartesian X: \f$(X,Y,Z)\f$,
    !! as well  as the magnitude, starting  from the input in  Flux coordinates.
    !! Note that  the Cartesian components  in X,  Y and Z  can be used  to plot
    !! the  real  vector  and  that  covariant components  should  be  equal  to
    !! contravariant components in this coordinate system.
    !!
    !! Also, the fluxes can be calculated and plot by providing a \c base_name.
    !!
    !! By  default, the  cartesian  components  are returned,  but  this can  be
    !! indicated differently by providing \c max_transf.
    !!
    !! \note
    !!  -# Plots  for different Richardson  levels can  be combined to  show the
    !!  total grid by just plotting them all individually and sequentially.
    !!  -# The metric factors and  transformation matrices have to be allocated.
    !!  They can be calculated using the routines from eq_ops, for
    !!  <tt>deriv = [0,0,0]</tt>.
    !!  -# For VMEC, the trigonometric factors of \c grid_XYZ must be calculated
    !!  beforehand.  Optionally, by  providing \c  X, \c  Y and  \c Z,  the ones
    !!  calculated  in this  routine are  overwritten. This  is useful  for, for
    !!  example, slab geometries.
    !!  -# The  normalization factors are taken  into account and the  output is
    !!  transformed back to unnormalized values.
    !!
    !! \return ierr
    integer function calc_vec_comp(grid,grid_eq,eq_1,eq_2,v_com,norm_disc_prec,&
        &v_mag,base_name,max_transf,v_flux_tor,v_flux_pol,XYZ,compare_tor_pos) &
        &result(ierr)
        
        use eq_vars, only: eq_1_type, eq_2_type, max_flux_f
        use eq_utilities, only: calc_inv_met
        use num_vars, only: eq_jobs_lims, eq_job_nr, use_pol_flux_f, eq_style, &
            &use_normalization, rank, tol_zero, rz_0, &
            &compare_tor_pos_glob => compare_tor_pos
        use num_utilities, only: c, calc_int, order_per_fun, spline
        use eq_vars, only: r_0, b_0, psi_0
        use vmec_utilities, only: calc_trigon_factors
        use mpi_utilities, only: get_ser_var
        
        character(*), parameter :: rout_name = 'calc_vec_comp'
        
        ! input / output
        type(grid_type), intent(in) :: grid                                     !< grid on which vector is calculated
        type(grid_type), intent(in) :: grid_eq                                  !< grid on which equilibrium variables are calculated
        type(eq_1_type), intent(in) :: eq_1                                     !< flux equilibrium quantities
        type(eq_2_type), intent(in) :: eq_2                                     !< metric equilibrium variables
        real(dp), intent(inout) :: v_com(:,:,:,:,:)                             !< covariant and contravariant components <tt>(dim1,dim2,dim3,3,2)</tt>
        integer, intent(in) :: norm_disc_prec                                   !< precision for normal derivatives
        real(dp), intent(inout), optional :: v_mag(:,:,:)                       !< magnitude <tt>(dim1,dim2,dim3)</tt>
        character(len=*), intent(in), optional :: base_name                     !< base name for output plotting
        integer, intent(in), optional :: max_transf                             !< maximum transformation level (2: Magnetic, 3: Equilibrium, 4: Cylindrical, 5: Cartesian [def])
        real(dp), intent(inout), allocatable, optional :: v_flux_pol(:,:)       !< poloidal flux as function of normal coordinate for all toroidal positions
        real(dp), intent(inout), allocatable, optional :: v_flux_tor(:,:)       !< toroidal flux as function of normal coordinate for all poloidal positions
        real(dp), intent(in), optional :: xyz(:,:,:,:)                          !< \f$X\f$, \f$Y\f$ and \f$Z\f$ of grid
        logical, intent(in), optional :: compare_tor_pos                        !< compare toroidal positions
        
        ! local variables
        type(grid_type) :: grid_trim                                            ! trimmed plot grid
        integer :: max_transf_loc                                               ! local max_transf
        integer :: id, jd, kd, ld                                               ! counters
        integer :: plot_dim(4)                                                  ! dimensions of plot
        integer :: plot_offset(4)                                               ! local offset of plot
        integer :: c_loc                                                        ! local c
        integer :: tor_id(2)                                                    ! toroidal indices
        integer :: norm_id(2)                                                   ! untrimmed normal indices for trimmed grids
        logical :: cont_plot                                                    ! continued plot
        logical :: do_plot                                                      ! perform plotting
        logical :: compare_tor_pos_loc                                          ! local compare_tor_pos
        character(len=max_str_ln) :: description(3)                             ! description of plots
        character(len=max_str_ln) :: file_names(3)                              ! plot file names
        character(len=max_str_ln) :: var_names(3,2)                             ! variable names
        character(len=5) :: coord_names(3)                                      ! name of coordinates
        character(len=max_str_ln) :: err_msg                                    ! error message
        real(dp) :: norm_len                                                    ! normalization factor for lengths, to cancel the one introduced in "calc_XYZ_grid"
        real(dp), allocatable :: q_saf(:,:)                                     ! interpolated q_saf_FD and derivative
        real(dp), allocatable :: jac(:,:,:)                                     ! interpolated jac_FD
        real(dp), allocatable :: xyzr(:,:,:,:)                                  ! X, Y, Z and R of surface in cylindrical coordinates, untrimmed grid
        real(dp), allocatable :: theta_geo(:,:,:)                               ! geometrical theta
        real(dp), allocatable :: x(:,:,:,:), y(:,:,:,:), z(:,:,:,:)             ! copy of X, Y and Z, trimmed grid
        real(dp), allocatable :: v_temp(:,:,:,:,:)                              ! temporary variable for v
        real(dp), allocatable :: v_ser_temp(:)                                  ! temporary serial variable
        real(dp), allocatable :: v_ser_temp_int(:)                              ! temporary integrated serial variable
        real(dp), allocatable :: t_ba(:,:,:,:,:,:,:), t_ab(:,:,:,:,:,:,:)       ! transformation matrices from A to B
        real(dp), allocatable :: d1r(:,:,:), d2r(:,:,:)                         ! dR/dpsi, dR/dtheta in HELENA coords.
        real(dp), allocatable :: d1z(:,:,:), d2z(:,:,:)                         ! dZ/dpsi, dZ/dtheta in HELENA coords.
        
        ! initialize ierr
        ierr = 0
        
        ! user output
        call writo('Prepare calculation of vector components')
        call lvl_ud(1)
        
        ! set up maximum level up to which to transform and whether to plot
        max_transf_loc = 5
        if (present(max_transf)) then
            if (max_transf.ge.1 .and. max_transf.le.5) &
                &max_transf_loc = max_transf
        end if
        do_plot = .false.
        if (present(base_name)) then
            if (trim(base_name).ne.'') do_plot = .true.
        end if
        
        ! set up continued plot
        cont_plot = eq_job_nr.gt.1
        
        ! set up normalization factor
        norm_len = 1._dp
        if (use_normalization) norm_len = r_0
        
        ! trim plot grid
        ierr = trim_grid(grid,grid_trim,norm_id)
        chckerr('')
        
        ! set up X, Y, Z and R in grid
        allocate(xyzr(grid%n(1),grid%n(2),grid%loc_n_r,4))
        ierr = calc_xyz_grid(grid_eq,grid,xyzr(:,:,:,1),xyzr(:,:,:,2),&
            &xyzr(:,:,:,3),r=xyzr(:,:,:,4))
        chckerr('')
        
        ! set up local compare_tor_pos
        compare_tor_pos_loc = compare_tor_pos_glob
        if (present(compare_tor_pos)) compare_tor_pos_loc = compare_tor_pos
        
        ! get geometrical poloidal angle if comparing toroidal position
        if (compare_tor_pos_loc) then
            allocate(theta_geo(grid%n(1),grid%n(2),grid%loc_n_r))
            theta_geo = atan2(xyzr(:,:,:,3)-rz_0(2),xyzr(:,:,:,4)-rz_0(1))
            where (theta_geo.lt.0) theta_geo = theta_geo + 2*pi
        end if
        
        ! if plotting is required
        if (do_plot) then
            ! set up plot dimensions and local dimensions
            plot_dim = [grid_trim%n,3]
            plot_offset = [0,0,grid_trim%i_min-1,0]
            tor_id = [1,size(v_com,2)]
            if (compare_tor_pos_loc) then
                if (plot_dim(2).ne.3) then
                    ierr = 1
                    err_msg = 'When comparing toroidal positions, need to &
                        &have 3 toroidal points (one in the middle)'
                    chckerr(err_msg)
                end if
                plot_dim(2) = 1
                tor_id = [2,2]
            end if
            
            ! possibly modify if multiple equilibrium parallel jobs
            if (size(eq_jobs_lims,2).gt.1) then
                plot_dim(1) = eq_jobs_lims(2,size(eq_jobs_lims,2)) - &
                    &eq_jobs_lims(1,1) + 1
                plot_offset(1) = eq_jobs_lims(1,eq_job_nr) - 1
                if (compare_tor_pos_loc) then
                    ierr = 1
                    err_msg = 'When comparing toroidal positions, cannot have &
                        &multiple equilibrium jobs'
                    chckerr(err_msg)
                end if
            end if
            
            ! tests
            if (eq_style.eq.1 .and. .not.allocated(grid%trigon_factors)) then
                ierr = 1
                err_msg = 'trigonometric factors not allocated'
                chckerr(err_msg)
            end if
            
            ! user output
            if (compare_tor_pos_loc) call writo('Comparing toroidal positions')
            
            ! copy X, Y and Z to trimmed grid copies
            allocate(x(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,1))
            allocate(y(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,1))
            allocate(z(grid_trim%n(1),grid_trim%n(2),grid_trim%loc_n_r,1))
            if (present(xyz)) then
                x(:,:,:,1) = xyz(:,:,norm_id(1):norm_id(2),1)
                y(:,:,:,1) = xyz(:,:,norm_id(1):norm_id(2),2)
                z(:,:,:,1) = xyz(:,:,norm_id(1):norm_id(2),3)
            else
                x(:,:,:,1) = xyzr(:,:,norm_id(1):norm_id(2),1)
                y(:,:,:,1) = xyzr(:,:,norm_id(1):norm_id(2),2)
                z(:,:,:,1) = xyzr(:,:,norm_id(1):norm_id(2),3)
            end if
        end if
        
        ! set up temporal interpolated copy of  v, T_BA and T_AB
        allocate(v_temp(grid%n(1),grid%n(2),grid%loc_n_r,3,2))
        allocate(t_ba(grid%n(1),grid%n(2),grid%loc_n_r,9,0:0,0:0,0:0))
        allocate(t_ab(grid%n(1),grid%n(2),grid%loc_n_r,9,0:0,0:0,0:0))
        allocate(q_saf(grid%loc_n_r,0:1))
        if ((present(v_flux_tor) .or. present(v_flux_pol)) .and. &
            &.not.compare_tor_pos_loc) then
            allocate(jac(grid%n(1),grid%n(2),grid%loc_n_r))
            if (rank.eq.0 .and. .not.cont_plot) then
                if (present(v_flux_tor)) then
                    allocate(v_flux_tor(grid_trim%n(3),plot_dim(2)))
                    v_flux_tor = 0._dp
                end if
                if (present(v_flux_pol)) then
                    allocate(v_flux_pol(grid_trim%n(3),plot_dim(1)))
                    v_flux_pol = 0._dp
                end if
            end if
        end if
        
        call lvl_ud(-1)
        
        ! 1. Flux coordinates (alpha,psi,theta)
        call writo('Flux coordinates')
        call lvl_ud(1)
        if (present(v_mag)) then
            v_mag = 0._dp
            do id = 1,3
                v_mag = v_mag + v_com(:,:,:,id,1)*v_com(:,:,:,id,2)
            end do
            v_mag = sqrt(v_mag)
        end if
        
        ! save temporary copy, normalized
        v_temp = v_com
        
        ! set up plot variables
        if (do_plot .and. .not.compare_tor_pos_loc) then                        ! comparison only works for periodic quantities
            coord_names(1) = 'alpha'
            coord_names(2) = 'psi'
            if (use_pol_flux_f) then
                coord_names(3) = 'theta'
            else
                coord_names(3) = 'zeta'
            end if
            var_names = trim(base_name)
            do id = 1,3
                var_names(id,1) = trim(var_names(id,1))//'_sub_'//&
                    &trim(coord_names(id))
                var_names(id,2) = trim(var_names(id,2))//'_sup_'//&
                    &trim(coord_names(id))
            end do
            description(1) = 'covariant components of the magnetic field in &
                &Flux coordinates'
            description(2) = 'contravariant components of the magnetic field &
                &in Flux coordinates'
            description(3) = 'magnitude of the magnetic field in Flux &
                &coordinates'
            file_names(1) = trim(base_name)//'_F_sub'
            file_names(2) = trim(base_name)//'_F_sup'
            file_names(3) = trim(base_name)//'_F_mag'
            if (use_normalization) then
                v_com(:,:,:,1,1) = v_com(:,:,:,1,1) * r_0                       ! norm factor for e_alpha
                v_com(:,:,:,2,1) = v_com(:,:,:,2,1) / (r_0*b_0)                 ! norm factor for e_psi
                v_com(:,:,:,3,1) = v_com(:,:,:,3,1) * r_0                       ! norm factor for e_theta
                v_com(:,:,:,1,2) = v_com(:,:,:,1,2) / r_0                       ! norm factor for e^alpha
                v_com(:,:,:,2,2) = v_com(:,:,:,2,2) * (r_0*b_0)                 ! norm factor for e^psi
                v_com(:,:,:,3,2) = v_com(:,:,:,3,2) / r_0                       ! norm factor for e^theta
            end if
            do id = 1,2
                call plot_hdf5(var_names(:,id),trim(file_names(id)),&
                    &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,id),&
                    &tot_dim=plot_dim,loc_offset=plot_offset,&
                    &x=x(:,tor_id(1):tor_id(2),:,:),&
                    &y=y(:,tor_id(1):tor_id(2),:,:),&
                    &z=z(:,tor_id(1):tor_id(2),:,:),&
                    &cont_plot=cont_plot,descr=description(id))
            end do
            if (present(v_mag)) then
                call plot_hdf5(trim(base_name),trim(file_names(3)),&
                    &v_mag(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2)),&
                    &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3),&
                    &x=x(:,tor_id(1):tor_id(2),:,1),&
                    &y=y(:,tor_id(1):tor_id(2),:,1),&
                    &z=z(:,tor_id(1):tor_id(2),:,1),&
                    &cont_plot=cont_plot,descr=description(3))
            end if
        end if
        
        call lvl_ud(-1)
        
        ! 2.   Magnetic  coordinates   (phi,theta,zeta)   by  converting   using
        ! transformation matrices:
        !   (v_i)_M = T_M^F (v_i)_F
        !   (v^i)_M = (v^i)_F T_F^M
        ! where
        !             ( -q' theta 1 0 )
        !   T_MF    = (    -q     0 1 )
        !             (     1     0 0 )
        ! if the poloidal flux is used, or
        !             ( -q'/q zeta q 0 )
        !   T_MF    = (    -1      0 0 )
        !             (    1/q     0 1 )
        ! if it is the toroidal flux. Its inverse can be calculated as well.
        call writo('magnetic coordinates')
        call lvl_ud(1)
        t_ba = 0._dp
        t_ab = 0._dp
        do id = 0,1
            ierr = spline(grid_eq%loc_r_F,eq_1%q_saf_FD(:,id),grid%loc_r_F,&
                &q_saf(:,id),ord=norm_disc_prec)
            chckerr('')
        end do
        c_loc = c([1,1],.false.)
        if (use_pol_flux_f) then
            t_ba(:,:,:,c_loc,0,0,0) = -grid%theta_F
            do kd = 1,grid%loc_n_r
                t_ba(:,:,kd,c_loc,0,0,0) = t_ba(:,:,kd,c_loc,0,0,0)*q_saf(kd,1)
                t_ba(:,:,kd,c([2,1],.false.),0,0,0) = -q_saf(kd,0)
            end do
            t_ba(:,:,:,c([3,1],.false.),0,0,0) = 1._dp
            t_ba(:,:,:,c([1,2],.false.),0,0,0) = 1._dp
            t_ba(:,:,:,c([2,3],.false.),0,0,0) = 1._dp
        else
            t_ba(:,:,:,c_loc,0,0,0) = -grid%zeta_F
            do kd = 1,grid%loc_n_r
                t_ba(:,:,kd,c_loc,0,0,0) = &
                    &t_ba(:,:,kd,c_loc,0,0,0)*q_saf(kd,1)/q_saf(kd,0)
                t_ba(:,:,kd,c([2,1],.false.),0,0,0) = 1._dp/q_saf(kd,0)
                t_ba(:,:,kd,c([1,2],.false.),0,0,0) = q_saf(kd,0)
            end do
            t_ba(:,:,:,c([2,1],.false.),0,0,0) = -1._dp
            t_ba(:,:,:,c([3,3],.false.),0,0,0) = 1._dp
        end if
        ierr = calc_inv_met(t_ab,t_ba,[0,0,0])
        chckerr('')
        v_com = 0._dp
        do jd = 1,3
            do id = 1,3
                v_com(:,:,:,jd,1) = v_com(:,:,:,jd,1) + t_ba(:,:,:,&
                    &c([jd,id],.false.),0,0,0) * v_temp(:,:,:,id,1)
                v_com(:,:,:,jd,2) = v_com(:,:,:,jd,2) + t_ab(:,:,:,&
                    &c([id,jd],.false.),0,0,0) * v_temp(:,:,:,id,2)
            end do
        end do
        if (present(v_mag)) then
            v_mag = 0._dp
            do id = 1,3
                v_mag = v_mag + v_com(:,:,:,id,1)*v_com(:,:,:,id,2)
            end do
            v_mag = sqrt(v_mag)
        end if
        
        ! set up plot variables
        if (do_plot) then
            coord_names(1) = 'psi'
            coord_names(2) = 'theta'
            coord_names(3) = 'zeta'
            var_names = trim(base_name)
            do id = 1,3
                var_names(id,1) = trim(var_names(id,1))//'_sub_'//&
                    &trim(coord_names(id))
                var_names(id,2) = trim(var_names(id,2))//'_sup_'//&
                    &trim(coord_names(id))
            end do
            description(1) = 'covariant components of the magnetic field in &
                &Magnetic coordinates'
            description(2) = 'contravariant components of the magnetic field &
                &in Magnetic coordinates'
            description(3) = 'magnitude of the magnetic field in Magnetic &
                &coordinates'
            file_names(1) = trim(base_name)//'_M_sub'
            file_names(2) = trim(base_name)//'_M_sup'
            file_names(3) = trim(base_name)//'_M_mag'
            if (compare_tor_pos_loc) then
                do id = 1,3
                    file_names(id) = trim(file_names(id))//'_COMP'
                end do
            end if
            if (use_normalization) then
                v_com(:,:,:,1,1) = v_com(:,:,:,1,1) / (r_0*b_0)                 ! norm factor for e_psi
                v_com(:,:,:,2,1) = v_com(:,:,:,2,1) * r_0                       ! norm factor for e_theta
                v_com(:,:,:,3,1) = v_com(:,:,:,3,1) * r_0                       ! norm factor for e_zeta
                v_com(:,:,:,1,2) = v_com(:,:,:,1,2) * (r_0*b_0)                 ! norm factor for e^psi
                v_com(:,:,:,2,2) = v_com(:,:,:,2,2) / r_0                       ! norm factor for e^theta
                v_com(:,:,:,3,2) = v_com(:,:,:,3,2) / r_0                       ! norm factor for e^zeta
            end if
            if (compare_tor_pos_loc) then
                ierr = calc_tor_diff(v_com,theta_geo,norm_disc_prec)
                chckerr('')
            end if
            do id = 1,2
                call plot_hdf5(var_names(:,id),trim(file_names(id)),&
                    &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,id),&
                    &tot_dim=plot_dim,loc_offset=plot_offset,&
                    &x=x(:,tor_id(1):tor_id(2),:,:),&
                    &y=y(:,tor_id(1):tor_id(2),:,:),&
                    &z=z(:,tor_id(1):tor_id(2),:,:),&
                    &cont_plot=cont_plot,descr=description(id))
            end do
            if (present(v_mag)) then
                if (compare_tor_pos_loc) then
                    ierr = calc_tor_diff(v_mag,theta_geo,norm_disc_prec)
                    chckerr('')
                end if
                call plot_hdf5(trim(base_name),trim(file_names(3)),&
                    &v_mag(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2)),&
                    &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3),&
                    &x=x(:,tor_id(1):tor_id(2),:,1),&
                    &y=y(:,tor_id(1):tor_id(2),:,1),&
                    &z=z(:,tor_id(1):tor_id(2),:,1),&
                    &cont_plot=cont_plot,descr=description(3))
            end if
        end if
        
        if ((present(v_flux_tor) .or. present(v_flux_pol)) .and. &
            &.not.compare_tor_pos_loc) then
            ! tests
            if (maxval(grid%theta_F(grid%n(1),:,:)).lt.&
                &minval(grid%theta_F(1,:,:))) then
                call writo('theta of the grid monotomically decreases in Flux &
                    &quantities.',alert=.true.)
                call lvl_ud(1)
                call writo('This inverts the sign of the toroidal flux.')
                call writo('Remember that the grid limits are prescribed &
                    &in Flux quantities.')
                call lvl_ud(-1)
            end if
            if (maxval(grid%zeta_F(:,grid%n(2),:)).lt.&
                &minval(grid%zeta_F(:,1,:))) then
                call writo('zeta of the grid monotomically decreases in Flux &
                    &quantities.',alert=.true.)
                call lvl_ud(1)
                call writo('This inverts the sign of the poloidal flux.')
                call writo('Remember that the grid limits are prescribed &
                    &in Flux quantities.')
                if (eq_style.eq.1) call writo('For VMEC, these are inverse.')
                call lvl_ud(-1)
            end if
            if (grid_trim%r_F(grid_trim%n(3)).lt.grid_trim%r_F(1)) then
                call writo('r of the grid monotomically decreases in Flux &
                    &quantities.',alert=.true.)
                call lvl_ud(1)
                call writo('This inverts the sign of the fluxes.')
                call writo('Remember that the grid limits are prescribed &
                    &in Flux quantities.')
                call lvl_ud(-1)
            end if
            if (abs(minval(grid_trim%r_F)).gt.tol_zero) then
                call writo('r of the grid does not start at zero.',alert=.true.)
                call lvl_ud(1)
                call writo('This leaves out part of the fluxes.')
                call lvl_ud(-1)
            end if
            
            ! initialize temporary variable
            allocate(v_ser_temp_int(grid_trim%n(3)))
            
            ! set up plot variables
            var_names(1,2) = 'integrated poloidal flux of '//trim(base_name)
            var_names(2,2) = 'integrated toroidal flux of '//trim(base_name)
            file_names(1) = trim(base_name)//'_flux_pol_int'
            file_names(2) = trim(base_name)//'_flux_tor_int'
            
#if ldebug
            if (debug_calc_vec_comp) then
                ! user output
                call writo('Instead of calculating fluxes, returning &
                    &integrals in the angular variables.',alert=.true.)
                call lvl_ud(1)
                    call writo('Useful to check whether Maxwell holds.')
                    call writo('i.e. whether loop integral of J returns &
                        &B flux')
                    call writo('Note that the J-variables now refer to B and &
                        &vice-versa')
                    call writo('Don''t forget the contribution of the toroidal &
                        &field B_zeta on axis, times 2piR.')
                    call writo('Don''t forget the vacuum permeability!')
                call lvl_ud(-1)
            else
#endif
                ! multiply angular contravariant coordinates by Jacobian
                do jd = 1,grid_eq%n(2)
                    do id = 1,grid_eq%n(1)
                        ierr = spline(grid_eq%loc_r_F,&
                            &eq_2%jac_FD(id,jd,:,0,0,0),grid%loc_r_F,&
                            &jac(id,jd,:),ord=norm_disc_prec)
                        chckerr('')
                    end do
                end do
                if (use_normalization) jac = jac*r_0/b_0
                v_com(:,:,:,2,2) = v_com(:,:,:,2,2)*jac
                v_com(:,:,:,3,2) = v_com(:,:,:,3,2)*jac
                deallocate(jac)
#if ldebug
            end if
#endif
            
            ! integrate toroidal flux
            if (present(v_flux_tor) .and. grid_trim%n(1).gt.1 .and. &
                &.not.compare_tor_pos_loc) then                                 ! integrate poloidally and normally
                ! loop over all toroidal points
                do jd = 1,grid_trim%n(2)
#if ldebug
                    if (debug_calc_vec_comp) then
                        ! for  all  normal   points,  integrate  poloidally  the
                        ! covariant poloidal quantity
                        do kd = norm_id(1),norm_id(2)
                            ierr = calc_int(v_com(:,jd,kd,2,1),&
                                &grid%theta_F(:,jd,kd),v_com(:,jd,kd,2,2))      ! save in contravariant variable
                            chckerr('')
                        end do
                        
                        ! gather on master
                        ierr = get_ser_var(v_com(grid_trim%n(1),jd,&
                            &norm_id(1):norm_id(2),2,2),v_ser_temp)
                        chckerr('')
                        
                        ! update integrals if master
                        if (rank.eq.0) then
                            if (use_pol_flux_f) then
                                v_flux_tor(:,jd) = v_flux_tor(:,jd) + v_ser_temp
                            else
                                v_flux_tor(:,jd+plot_offset(1)) = v_ser_temp
                            end if
                        end if
                    else
#endif
                        ! for all normal points, integrate poloidally
                        do kd = norm_id(1),norm_id(2)
                            ierr = calc_int(v_com(:,jd,kd,3,2),&
                                &grid%theta_F(:,jd,kd),v_com(:,jd,kd,3,1))      ! save in covariant variable
                            chckerr('')
                        end do
                        
                        ! gather on master
                        ierr = get_ser_var(v_com(grid_trim%n(1),jd,&
                            &norm_id(1):norm_id(2),3,1),v_ser_temp)
                        chckerr('')
                        
                        ! integrate normally and update integrals if master
                        if (rank.eq.0) then
                            ierr = calc_int(v_ser_temp,&
                                &grid_trim%r_F(norm_id(1):),v_ser_temp_int)
                            chckerr('')
                            if (use_normalization) v_ser_temp_int = &
                                &v_ser_temp_int*psi_0
                            if (use_pol_flux_f) then
                                v_flux_tor(:,jd) = v_flux_tor(:,jd) + &
                                    &v_ser_temp_int
                            else
                                v_flux_tor(:,jd+plot_offset(1)) = v_ser_temp_int
                            end if
                        end if
#if ldebug
                    end if
#endif
                end do
                
                ! plot output
                if (rank.eq.0 .and. do_plot .and. &
                    &eq_job_nr.eq.size(eq_jobs_lims,2)) then
                    call print_ex_2d([var_names(2,2)],&
                        &trim(file_names(2)),v_flux_tor,&
                        &x=reshape(grid_trim%r_F(norm_id(1):)*2*pi/max_flux_f,&
                        &[grid_trim%n(3),1]),draw=.false.)
                    call draw_ex([var_names(2,2)],trim(file_names(2)),&
                        &plot_dim(2),1,.false.)
                end if
            end if
            
            ! integrate poloidal flux
            if (present(v_flux_pol) .and. grid_trim%n(2).gt.1 .and. &
                &.not.compare_tor_pos_loc) then                                 ! integrate toroidally and normally
                ! loop over all poloidal points
                do id = 1,grid_trim%n(1)
#if ldebug
                    if (debug_calc_vec_comp) then
                        ! for  all  normal   points,  integrate  toroidally  the
                        ! covariant toroidal quantity
                        do kd = norm_id(1),norm_id(2)
                            ierr = calc_int(v_com(id,:,kd,3,1),&
                                &grid%zeta_F(id,:,kd),v_com(id,:,kd,3,2))       ! save in contravariant variable
                            chckerr('')
                        end do
                        
                        ! gather on master
                        ierr = get_ser_var(v_com(id,grid_trim%n(2),&
                            &norm_id(1):norm_id(2),3,2),v_ser_temp)
                        chckerr('')
                        
                        ! update integrals if master
                        if (rank.eq.0) then
                            if (use_pol_flux_f) then
                                v_flux_pol(:,id+plot_offset(1)) = v_ser_temp
                            else
                                v_flux_pol(:,id) = v_flux_pol(:,id) + v_ser_temp
                            end if
                        end if
                    else
#endif
                        ! for all normal points, integrate toroidally
                        do kd = norm_id(1),norm_id(2)
                            ierr = calc_int(v_com(id,:,kd,2,2),&
                                &grid%zeta_F(id,:,kd),v_com(id,:,kd,2,1))       ! save in covariant variable
                            chckerr('')
                        end do
                        
                        ! gather on master
                        ierr = get_ser_var(v_com(id,grid_trim%n(2),&
                            &norm_id(1):norm_id(2),2,1),v_ser_temp)
                        chckerr('')
                        
                        ! integrate normally and update integrals if master
                        if (rank.eq.0) then
                            ierr = calc_int(v_ser_temp,&
                                &grid_trim%r_F(norm_id(1):),v_ser_temp_int)
                            chckerr('')
                            if (use_normalization) v_ser_temp_int = &
                                &v_ser_temp_int*psi_0
                            if (use_pol_flux_f) then
                                v_flux_pol(:,id+plot_offset(1)) = v_ser_temp_int
                            else
                                v_flux_pol(:,id) = v_flux_pol(:,id) + &
                                    &v_ser_temp_int
                            end if
                        end if
#if ldebug
                    end if
#endif
                end do
                
                ! plot output
                if (rank.eq.0 .and. do_plot .and. &
                    &eq_job_nr.eq.size(eq_jobs_lims,2)) then
                    call print_ex_2d([var_names(1,2)],&
                        &trim(file_names(1)),v_flux_pol,&
                        &x=reshape(grid_trim%r_F(norm_id(1):)*2*pi/max_flux_f,&
                        &[grid_trim%n(3),1]),draw=.false.)
                    call draw_ex([var_names(1,2)],trim(file_names(1)),&
                        &plot_dim(1),1,.false.)
                end if
            end if
            
            ! clean up
            deallocate(v_ser_temp)
            if (rank.eq.0) deallocate(v_ser_temp_int)
        end if
        
        call lvl_ud(-1)
        
        if (max_transf_loc.eq.2) return
        
        ! 3.   HELENA   coordinates   (psi,theta,zeta)   or   VMEC   coordinates
        ! (r,theta,phi), by converting using transformation matrices
        !   (v_i)_H = T_H^F (v_i)_F
        !   (v^i)_H = (v^i)_F T_F^H
        ! Note: This is  done directly from step 1; The results  from step 2 are
        ! skipped, i.e. v_temp is not overwritten after step 2.
        call writo('Equilibrium coordinates')
        call lvl_ud(1)
        do ld = 1,9
            do jd = 1,grid_eq%n(2)
                do id = 1,grid_eq%n(1)
                    ierr = spline(grid_eq%loc_r_F,&
                        &eq_2%T_FE(id,jd,:,ld,0,0,0),grid%loc_r_F,&
                        &t_ab(id,jd,:,ld,0,0,0),ord=norm_disc_prec)
                    chckerr('')
                    ierr = spline(grid_eq%loc_r_F,&
                        &eq_2%T_EF(id,jd,:,ld,0,0,0),grid%loc_r_F,&
                        &t_ba(id,jd,:,ld,0,0,0),ord=norm_disc_prec)
                    chckerr('')
                end do
            end do
        end do
        v_com = 0._dp
        do jd = 1,3
            do id = 1,3
                v_com(:,:,:,jd,1) = v_com(:,:,:,jd,1) + &
                    &t_ba(:,:,:,c([jd,id],.false.),0,0,0) * &
                    &v_temp(:,:,:,id,1)
                v_com(:,:,:,jd,2) = v_com(:,:,:,jd,2) + &
                    &t_ab(:,:,:,c([id,jd],.false.),0,0,0) * &
                    &v_temp(:,:,:,id,2)
            end do
        end do
        if (present(v_mag)) then
            v_mag = 0._dp
            do id = 1,3
                v_mag = v_mag + v_com(:,:,:,id,1)*v_com(:,:,:,id,2)
            end do
            v_mag = sqrt(v_mag)
        end if
        
        ! save temporary copy, normalized
        v_temp = v_com
        
        ! set up plot variables
        if (do_plot) then
            select case (eq_style)
                case (1)                                                        ! VMEC
                    coord_names(1) = 'r'
                    coord_names(2) = 'theta'
                    coord_names(3) = 'phi'
                    description(1) = 'covariant components of the magnetic &
                        &field in VMEC coordinates'
                    description(2) = 'contravariant components of the magnetic &
                        &field in VMEC coordinates'
                    description(3) = 'magnitude of the magnetic field in VMEC &
                        &coordinates'
                    file_names(1) = trim(base_name)//'_V_sub'
                    file_names(2) = trim(base_name)//'_V_sup'
                    file_names(3) = trim(base_name)//'_V_mag'
                case (2)                                                        ! HELENA
                    coord_names(1) = 'psi'
                    coord_names(2) = 'theta'
                    coord_names(3) = 'phi'
                    description(1) = 'covariant components of the magnetic &
                        &field in HELENA coordinates'
                    description(2) = 'contravariant components of the magnetic &
                        &field in HELENA coordinates'
                    description(3) = 'magnitude of the magnetic field in &
                        &HELENA coordinates'
                    file_names(1) = trim(base_name)//'_H_sub'
                    file_names(2) = trim(base_name)//'_H_sup'
                    file_names(3) = trim(base_name)//'_H_mag'
            end select
            if (compare_tor_pos_loc) then
                do id = 1,3
                    file_names(id) = trim(file_names(id))//'_COMP'
                end do
            end if
            var_names = trim(base_name)
            do id = 1,3
                var_names(id,1) = trim(var_names(id,1))//'_sub_'//&
                    &trim(coord_names(id))
                var_names(id,2) = trim(var_names(id,2))//'_sup_'//&
                    &trim(coord_names(id))
            end do
            if (use_normalization) then
                v_com(:,:,:,1,1) = v_com(:,:,:,1,1) * r_0                       ! norm factor for e_r or e_psi
                v_com(:,:,:,2,1) = v_com(:,:,:,2,1) * r_0                       ! norm factor for e_theta
                v_com(:,:,:,3,1) = v_com(:,:,:,3,1) * r_0                       ! norm factor for e_zeta or e_phi
                v_com(:,:,:,1,2) = v_com(:,:,:,1,2) / r_0                       ! norm factor for e^r or e^psi
                v_com(:,:,:,2,2) = v_com(:,:,:,2,2) / r_0                       ! norm factor for e^theta
                v_com(:,:,:,3,2) = v_com(:,:,:,3,2) / r_0                       ! norm factor for e^zeta or e^phi
            end if
            if (compare_tor_pos_loc) then
                ierr = calc_tor_diff(v_com,theta_geo,norm_disc_prec)
                chckerr('')
            end if
            do id = 1,2
                call plot_hdf5(var_names(:,id),trim(file_names(id)),&
                    &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,id),&
                    &tot_dim=plot_dim,loc_offset=plot_offset,&
                    &x=x(:,tor_id(1):tor_id(2),:,:),&
                    &y=y(:,tor_id(1):tor_id(2),:,:),&
                    &z=z(:,tor_id(1):tor_id(2),:,:),&
                    &cont_plot=cont_plot,descr=description(id))
            end do
            if (present(v_mag)) then
                if (compare_tor_pos_loc) then
                    ierr = calc_tor_diff(v_mag,theta_geo,norm_disc_prec)
                    chckerr('')
                end if
                call plot_hdf5(trim(base_name),trim(file_names(3)),&
                    &v_mag(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2)),&
                    &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3),&
                    &x=x(:,tor_id(1):tor_id(2),:,1),&
                    &y=y(:,tor_id(1):tor_id(2),:,1),&
                    &z=z(:,tor_id(1):tor_id(2),:,1),&
                    &cont_plot=cont_plot,descr=description(3))
            end if
        end if
        
        call lvl_ud(-1)
        
        if (max_transf_loc.eq.3) return
        
        ! 4.   cylindrical    coordinates   (R,phi,Z)   by    converting   using
        ! transformation matrices:
        !   (v_i)_C = T_C^E (v_i)_E
        !   (v^i)_C = (v^i)_E T_E^C
        ! where for VMEC (E=V), the transformation matrix T_VC is tabulated, and
        ! its inverse  is calculate,  and for  HELENA (E=H),  the transformation
        ! matrix T_HC is given by
        !             (    R'    0    Z'   )
        !   T_HC    = ( R_theta  0 Z_theta )
        !             (    0    -1    0    )
        ! and its inverse can be calculated as well.
        call writo('Cylindrical coordinates')
        call lvl_ud(1)
        t_ba = 0._dp
        t_ab = 0._dp
        select case (eq_style)
            case (1)                                                            ! VMEC
                do ld = 1,9
                    do jd = 1,grid_eq%n(2)
                        do id = 1,grid_eq%n(1)
                            ierr = spline(grid_eq%loc_r_F,&
                                &eq_2%T_VC(id,jd,:,ld,0,0,0),grid%loc_r_F,&
                                &t_ab(id,jd,:,ld,0,0,0),ord=norm_disc_prec)
                            chckerr('')
                        end do
                    end do
                end do
            case (2)                                                            ! HELENA
                ! setup D1R and D2R, and same thing for Z
                ! (use untrimmed grid, with ghost region)
                allocate(d1r(grid%n(1),grid%n(2),grid%loc_n_r))
                allocate(d2r(grid%n(1),grid%n(2),grid%loc_n_r))
                allocate(d1z(grid%n(1),grid%n(2),grid%loc_n_r))
                allocate(d2z(grid%n(1),grid%n(2),grid%loc_n_r))
                do jd = 1,grid%n(2)
                    do id = 1,grid%n(1)
                        ierr = spline(grid%loc_r_E,xyzr(id,jd,:,4)/norm_len,&
                            &grid%loc_r_E,d1r(id,jd,:),ord=norm_disc_prec,&
                            &deriv=1)
                        chckerr('')
                        ierr = spline(grid%loc_r_E,xyzr(id,jd,:,3)/norm_len,&
                            &grid%loc_r_E,d1z(id,jd,:),ord=norm_disc_prec,&
                            &deriv=1)
                        chckerr('')
                    end do
                end do
                do kd = 1,grid%loc_n_r
                    do jd = 1,grid%n(2)
                        ierr = spline(grid%theta_E(:,jd,kd),&
                            &xyzr(:,jd,kd,4)/norm_len,grid%theta_E(:,jd,kd),&
                            &d2r(:,jd,kd),ord=norm_disc_prec,deriv=1)
                        chckerr('')
                        ierr = spline(grid%theta_E(:,jd,kd),&
                            &xyzr(:,jd,kd,3)/norm_len,grid%theta_E(:,jd,kd),&
                            &d2z(:,jd,kd),ord=norm_disc_prec,deriv=1)
                        chckerr('')
                    end do
                end do
                t_ab(:,:,:,c([1,1],.false.),0,0,0) = d1r
                t_ab(:,:,:,c([1,3],.false.),0,0,0) = d1z
                t_ab(:,:,:,c([2,1],.false.),0,0,0) = d2r
                t_ab(:,:,:,c([2,3],.false.),0,0,0) = d2z
                t_ab(:,:,:,c([3,2],.false.),0,0,0) = -1._dp
                deallocate(d1r,d2r,d1z,d2z)
        end select
        ierr = calc_inv_met(t_ba,t_ab,[0,0,0])
        chckerr('')
        v_com = 0._dp
        do jd = 1,3
            do id = 1,3
                v_com(:,:,:,jd,1) = v_com(:,:,:,jd,1) + t_ba(:,:,:,&
                    &c([jd,id],.false.),0,0,0) * v_temp(:,:,:,id,1)
                v_com(:,:,:,jd,2) = v_com(:,:,:,jd,2) + t_ab(:,:,:,&
                    &c([id,jd],.false.),0,0,0) * v_temp(:,:,:,id,2)
            end do
        end do
        if (present(v_mag)) then
            v_mag = 0._dp
            do id = 1,3
                v_mag = v_mag + v_com(:,:,:,id,1)*v_com(:,:,:,id,2)
            end do
            v_mag = sqrt(v_mag)
        end if
        
        ! save temporary copy, normalized
        v_temp = v_com
        
        ! set up plot variables
        if (do_plot) then
            description(1) = 'covariant components of the magnetic field in &
                &cylindrical coordinates'
            description(2) = 'contravariant components of the magnetic field &
                &in cylindrical coordinates'
            description(3) = 'magnitude of the magnetic field in cylindrical &
                &coordinates'
            file_names(1) = trim(base_name)//'_C_sub'
            file_names(2) = trim(base_name)//'_C_sup'
            file_names(3) = trim(base_name)//'_C_mag'
            if (compare_tor_pos_loc) then
                do id = 1,3
                    file_names(id) = trim(file_names(id))//'_COMP'
                end do
            end if
            coord_names(1) = 'R'
            coord_names(2) = 'phi'
            coord_names(3) = 'Z'
            var_names = trim(base_name)
            do id = 1,3
                var_names(id,1) = trim(var_names(id,1))//'_sub_'//&
                    &trim(coord_names(id))
                var_names(id,2) = trim(var_names(id,2))//'_sup_'//&
                    &trim(coord_names(id))
            end do
            if (use_normalization) then
                !v_com(:,:,:,1,1) = v_com(:,:,:,1,1)                             ! norm factor for e_R
                v_com(:,:,:,2,1) = v_com(:,:,:,2,1) * r_0                       ! norm factor for e_phi
                !v_com(:,:,:,3,1) = v_com(:,:,:,3,1)                             ! norm factor for e_Z
                !v_com(:,:,:,1,2) = v_com(:,:,:,1,2)                             ! norm factor for e^R
                v_com(:,:,:,2,2) = v_com(:,:,:,2,2) / r_0                       ! norm factor for e^phi
                !v_com(:,:,:,3,2) = v_com(:,:,:,3,2)                             ! norm factor for e^Z
            end if
            if (compare_tor_pos_loc) then
                ierr = calc_tor_diff(v_com,theta_geo,norm_disc_prec)
                chckerr('')
            end if
            do id = 1,2
                call plot_hdf5(var_names(:,id),trim(file_names(id)),&
                    &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,id),&
                    &tot_dim=plot_dim,loc_offset=plot_offset,&
                    &x=x(:,tor_id(1):tor_id(2),:,:),&
                    &y=y(:,tor_id(1):tor_id(2),:,:),&
                    &z=z(:,tor_id(1):tor_id(2),:,:),&
                    &cont_plot=cont_plot,descr=description(id))
            end do
            if (present(v_mag)) then
                if (compare_tor_pos_loc) then
                    ierr = calc_tor_diff(v_mag,theta_geo,norm_disc_prec)
                    chckerr('')
                end if
                call plot_hdf5(trim(base_name),trim(file_names(3)),&
                    &v_mag(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2)),&
                    &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3),&
                    &x=x(:,tor_id(1):tor_id(2),:,1),&
                    &y=y(:,tor_id(1):tor_id(2),:,1),&
                    &z=z(:,tor_id(1):tor_id(2),:,1),&
                    &cont_plot=cont_plot,descr=description(3))
            end if
        end if
        
        call lvl_ud(-1)
        
        if (max_transf_loc.eq.4) return
        
        ! 5. Cartesian  coordinates (X,Y,Z)  by converting  using transformation
        ! matrices
        !    (e_R)    (   cos(phi)   sin(phi)  0 ) (e_X)
        !   (e_phi) = ( -R sin(phi) R cos(phi) 0 ) (e_Y)
        !    (e_Z)    (     0          0       1 ) (e_Z)
        ! with  phi  =  -  zeta_F.   For  the  transformation  of  contravariant
        ! components, the factor R becomes 1/R and the transpose is taken.
        call writo('Cartesian coordinates')
        call lvl_ud(1)
        t_ba = 0._dp
        t_ab = 0._dp
        t_ab(:,:,:,c([1,1],.false.),0,0,0) = cos(-grid%zeta_F)
        t_ab(:,:,:,c([1,2],.false.),0,0,0) = sin(-grid%zeta_F)
        t_ab(:,:,:,c([2,1],.false.),0,0,0) = -xyzr(:,:,:,4)/norm_len*&
            &sin(-grid%zeta_F)
        t_ab(:,:,:,c([2,2],.false.),0,0,0) = xyzr(:,:,:,4)/norm_len*&
            &cos(-grid%zeta_F)
        t_ab(:,:,:,c([3,3],.false.),0,0,0) = 1._dp
        ierr = calc_inv_met(t_ba,t_ab,[0,0,0])
        chckerr('')
        v_com = 0._dp
        do jd = 1,3
            do id = 1,3
                v_com(:,:,:,jd,1) = v_com(:,:,:,jd,1) + t_ba(:,:,:,&
                    &c([jd,id],.false.),0,0,0) * v_temp(:,:,:,id,1)
                v_com(:,:,:,jd,2) = v_com(:,:,:,jd,2) + t_ab(:,:,:,&
                    &c([id,jd],.false.),0,0,0) * v_temp(:,:,:,id,2)
            end do
        end do
        if (present(v_mag)) then
            v_mag = 0._dp
            do id = 1,3
                v_mag = v_mag + v_com(:,:,:,id,1)*v_com(:,:,:,id,2)
            end do
            v_mag = sqrt(v_mag)
        end if
        
        ! set up plot variables
        if (do_plot) then
            description(1) = 'covariant components of the magnetic field in &
                &cartesian coordinates'
            description(2) = 'contravariant components of the magnetic field &
                &in cartesian coordinates'
            description(3) = 'magnitude of the magnetic field in cartesian &
                &coordinates'
            file_names(1) = trim(base_name)//'_X_sub'
            file_names(2) = trim(base_name)//'_X_sup'
            file_names(3) = trim(base_name)//'_X_mag'
            if (compare_tor_pos_loc) then
                do id = 1,3
                    file_names(id) = trim(file_names(id))//'_COMP'
                end do
            end if
            coord_names(1) = 'X'
            coord_names(2) = 'Y'
            coord_names(3) = 'Z'
            var_names = trim(base_name)
            do id = 1,3
                var_names(id,1) = trim(var_names(id,1))//'_sub_'//&
                    &trim(coord_names(id))
                var_names(id,2) = trim(var_names(id,2))//'_sup_'//&
                    &trim(coord_names(id))
            end do
            if (compare_tor_pos_loc) then
                ierr = calc_tor_diff(v_com,theta_geo,norm_disc_prec)
                chckerr('')
            end if
            do id = 1,2
                call plot_hdf5(var_names(:,id),trim(file_names(id)),&
                    &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,id),&
                    &tot_dim=plot_dim,loc_offset=plot_offset,&
                    &x=x(:,tor_id(1):tor_id(2),:,:),&
                    &y=y(:,tor_id(1):tor_id(2),:,:),&
                    &z=z(:,tor_id(1):tor_id(2),:,:),&
                    &cont_plot=cont_plot,descr=description(id))
            end do
            if (present(v_mag)) then
                if (compare_tor_pos_loc) then
                    ierr = calc_tor_diff(v_mag,theta_geo,norm_disc_prec)
                    chckerr('')
                end if
                call plot_hdf5(trim(base_name),trim(file_names(3)),&
                    &v_mag(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2)),&
                    &tot_dim=plot_dim(1:3),loc_offset=plot_offset(1:3),&
                    &x=x(:,tor_id(1):tor_id(2),:,1),&
                    &y=y(:,tor_id(1):tor_id(2),:,1),&
                    &z=z(:,tor_id(1):tor_id(2),:,1),&
                    &cont_plot=cont_plot,descr=description(3))
            end if
            
            ! plot vector as well
            call plot_hdf5([trim(base_name)],trim(base_name)//'_vec',&
                &v_com(:,tor_id(1):tor_id(2),norm_id(1):norm_id(2),:,1),&
                &tot_dim=plot_dim,loc_offset=plot_offset,&
                &x=x(:,tor_id(1):tor_id(2),:,:),y=y(:,tor_id(1):tor_id(2),:,:),&
                &z=z(:,tor_id(1):tor_id(2),:,:),col=4,cont_plot=cont_plot,&
                &descr='magnetic field vector')
        end if
        
        call lvl_ud(-1)
        
        ! clean up
        call grid_trim%dealloc()
    end function  calc_vec_comp
    
    !> Calculates the local  number of parallel grid points  for this Richardson
    !! level, taking into account that it ould be half the actual number.
    !!
    !! \see calc_ang_grid_eq_b().
    !!
    !! \return ierr
    integer function calc_n_par_x_rich(n_par_X_rich,only_half_grid) result(ierr)
        use rich_vars, only: n_par_x
        
        character(*), parameter :: rout_name = 'calc_n_par_X_rich'
        
        ! input / output
        integer, intent(inout) :: n_par_x_rich                                  !< n_par_X for this Richardson level
        logical, intent(in), optional :: only_half_grid                         !< calculate only half grid with even points
        
        ! local variables
        character(len=max_str_ln) :: err_msg                                    ! error message
        
        ! initialize ierr
        ierr = 0
        
        ! possibly divide n_par_X by 2
        n_par_x_rich = n_par_x
        if (present(only_half_grid)) then
            if (only_half_grid) then
                if (mod(n_par_x,2).eq.1) then
                    n_par_x_rich = (n_par_x-1)/2
                else
                    ierr = 1
                    err_msg = 'Need odd number of points'
                    chckerr(err_msg)
                end if
            end if
        end if
    end function calc_n_par_x_rich
    
    !> calculates the  cosine and sine mode  numbers of a function  defined on a
    !! non-regular grid.
    !!
    !! If a plot name is provided, the modes are plotted.
    !!
    !! \note The  fundamental interval is  assumed to be \f$0\ldots  2\pi\f$ but
    !! there should be no overlap between the first and last point.
    !!
    !! \return ierr
    integer function nufft(x,f,f_F,plot_name) result(ierr)
        use num_utilities, only: order_per_fun, spline
        
        character(*), parameter :: rout_name = 'nufft'
        
        ! input / output
        real(dp), intent(in) :: x(:)                                            !< coordinate values
        real(dp), intent(in) :: f(:)                                            !< function values
        real(dp), intent(inout), allocatable :: f_f(:,:)                        !< Fourier modes for \f$\cos\f$ and \f$\sin\f$
        character(len=*), intent(in), optional :: plot_name                     !< name of possible plot
        
        ! local variables
        integer :: m_f                                                          ! nr. of modes
        integer :: n_x                                                          ! nr. of points
        integer :: id                                                           ! counter
        logical :: print_log = .false.                                          ! print log plot as well
        real(dp), allocatable :: x_loc(:)                                       ! local x, extended
        real(dp), allocatable :: work(:)                                        ! work array
        real(dp), allocatable :: f_loc(:)                                       ! local f, extended
        real(dp), allocatable :: f_int(:)                                       ! interpolated f, later Fourier modes
        character(len=max_str_ln) :: plot_title(2)                              ! name of plot
        
        ! tests
        if (size(x).ne.size(f)) then
            ierr = 1
            chckerr('x and f must have same size')
        end if
        if (maxval(x)-minval(x).gt.2*pi) then
            ierr = 1
            chckerr('x is not a fundamental interval')
        end if
        
        ! create local x and f that go from <0 to <2pi, with a royal overlap for
        ! interpolation
        ierr = order_per_fun(x,f,x_loc,f_loc,10)
        chckerr('')
        
        ! set local f and interpolate
        n_x = size(x)
        allocate(f_int(n_x))
        ierr = spline(x_loc,f_loc,[((id-1._dp)/n_x*2*pi,id=1,n_x)],f_int,&
            &ord=3)
        chckerr('')
        
        ! clean up
        deallocate(x_loc)
        deallocate(f_loc)
        
        !!do id = 1,n_x
            !!f_int(id) = -4._dp+&
                !!&3*cos((id-1._dp)/n_x*2*pi*1)+&
                !!&0.5*cos((id-1._dp)/n_x*2*pi*100)+&
                !!&1.5*cos((id-1._dp)/n_x*2*pi*101)+&
                !!&4*sin((id-1._dp)/n_x*2*pi*1)+&
                !!&2*sin((id-1._dp)/n_x*2*pi*50)+&
                !!&4.5*sin((id-1._dp)/n_x*2*pi*55)
        !!end do
        
        ! set up variables for fft
        m_f = (n_x-1)/2                                                         ! nr. of modes
        allocate(work(3*n_x+15))                                                ! from fftpack manual
        
        ! calculate fft
        call dffti(n_x,work)
        call dfftf(n_x,f_int,work)
        deallocate(work)
        
        ! rescale
        f_int(:) = f_int(:)*2/n_x
        f_int(1) = f_int(1)/2
        
        ! separate cos and sine
        if (allocated(f_f)) deallocate(f_f)
        allocate(f_f(m_f+1,2))
        f_f(1,1) = f_int(1)
        f_f(1,2) = 0._dp
        f_f(2:m_f+1,1) = f_int(2:2*m_f:2)
        f_f(2:m_f+1,2) = -f_int(3:2*m_f+1:2)                                    ! routine returns - sine factors
        
        ! output in plot if requested
        if (present(plot_name)) then
            plot_title = ['cos','sin']
            call print_ex_2d(plot_title,plot_name,f_f,draw=.false.)
            call draw_ex(plot_title,plot_name,2,1,.false.)
            if (print_log) then
                plot_title = ['cos [log]','sin [log]']
                call print_ex_2d(plot_title,trim(plot_name)//'_log',&
                    &log10(abs(f_f)),draw=.false.)
                call draw_ex(plot_title,trim(plot_name)//'_log',2,1,.false.)
            end if
        end if
    end function nufft
    
    !> finds smallest range that comprises a minimum and maximum value.
    subroutine find_compr_range(r_F,lim_r,lim_id)
        ! input / output
        real(dp), intent(in) :: r_f(:)                                          !< all values of coordinate
        real(dp), intent(in) :: lim_r(2)                                        !< limiting range
        integer, intent(inout) :: lim_id(2)                                     !< limiting indices
        
        ! local variables
        integer :: id                                                           ! counter
        integer :: n_r                                                          ! number of points in coordinate
        real(dp), parameter :: tol = 1.e-6                                      ! tolerance for grids
        
        ! set n_r
        n_r = size(r_f)
        
        lim_id = [1,n_r]                                                        ! initialize with full range
        if (r_f(1).lt.r_f(n_r)) then                                            ! ascending r_F
            do id = 1,n_r
                if (r_f(id).le.minval(lim_r)-tol) lim_id(1) = id                ! move lower limit up
                if (r_f(n_r-id+1).ge.maxval(lim_r)+tol) &
                    &lim_id(2) = n_r-id+1                                       ! move upper limit down
            end do
        else                                                                    ! descending r_F
            do id = 1,n_r
                if (r_f(id).ge.maxval(lim_r)+tol) lim_id(1) = id                ! move lower limit up
                if (r_f(n_r-id+1).le.minval(lim_r)-tol) &
                    &lim_id(2) = n_r-id+1                                       ! move upper limit down
            end do
        end if
    end subroutine find_compr_range
    
    !> Calculate arclength angle.
    !!
    !! This   angle   is   based   on   the   calculation   of   the   arclength
    !! from   the  start   of   the   grid,  and   then   normalizing  it   from
    !! \f$\min(\gamma)\ldots\max(\gamma)\f$, where \f$\gamma\f$  is the parallel
    !! angle.
    !!
    !! By default, the Flux variables are used,  but this can be changed with \c
    !! use_E.
    !!
    !! \note For  field-aligned grids  the projection is  taken on  the poloidal
    !! cross-section. For non-axisymmetric equilibria, this makes little sense.
    !!
    !! \return ierr
    integer function calc_arc_angle(grid,eq_1,eq_2,arc,use_E) result(ierr)
        use eq_vars, only: eq_1_type, eq_2_type, &
            &r_0
        use num_vars, only: use_normalization, use_pol_flux_f
        use num_utilities, only: c, calc_int
        
        character(*), parameter :: rout_name = 'calc_arc_angle'
        
        ! input / output
        type(grid_type), intent(in) :: grid                                     !< equilibrium grid
        type(eq_1_type), intent(in) :: eq_1                                     !< flux equilibrium variables
        type(eq_2_type), intent(in) :: eq_2                                     !< metric equilibrium variables
        real(dp), intent(inout), allocatable :: arc(:,:,:)                      !< arclength angle
        logical, intent(in), optional :: use_e                                  !< use E variables instead of F
        
        ! local variables
        integer :: jd, kd                                                       ! counters
        logical :: use_e_loc                                                    ! local use_E
        
        ! initialize ierr
        ierr = 0
        
        ! set local use_E
        use_e_loc = .false.
        if (present(use_e)) use_e_loc = use_e
        
        ! calculate arclength: int dtheta |e_theta|
        allocate(arc(grid%n(1),grid%n(2),grid%loc_n_r))
        do kd = 1,grid%loc_n_r
            do jd = 1,grid%n(2)
                if (use_e_loc) then
                    ierr = calc_int(&
                        &eq_2%g_E(:,jd,kd,c([2,2],.true.),0,0,0)**(0.5),&
                        &grid%theta_E(:,jd,kd),arc(:,jd,kd))
                    chckerr('')
                    if (use_normalization) arc(:,jd,kd) = arc(:,jd,kd) * r_0    ! not really necessary as we take fractional arc length
                    arc(:,jd,kd) = grid%theta_E(1,jd,kd) + &
                        &arc(:,jd,kd) / arc(grid%n(1),jd,kd) * &
                        &(grid%theta_E(grid%n(1),jd,kd)-grid%theta_E(1,jd,kd))
                else
                    if (use_pol_flux_f) then
                        ! e_arc = e_theta - q e_alpha
                        ierr = calc_int((&
                            &eq_2%g_FD(:,jd,kd,c([3,3],.true.),0,0,0) - &
                            &eq_2%g_FD(:,jd,kd,c([1,3],.true.),0,0,0) * &
                            &2*eq_1%q_saf_FD(kd,0) + &
                            &eq_2%g_FD(:,jd,kd,c([1,1],.true.),0,0,0) * &
                            &eq_1%q_saf_FD(kd,0)**2)**(0.5),&
                            &grid%theta_F(:,jd,kd),arc(:,jd,kd))
                        chckerr('')
                    else
                        ! e_arc = -e_alpha
                        ierr = calc_int(&
                            &(-eq_2%g_FD(:,jd,kd,c([1,1],.true.),0,0,0))&
                            &**(0.5),grid%theta_F(:,jd,kd),arc(:,jd,kd))
                        chckerr('')
                    end if
                    if (use_normalization) arc(:,jd,kd) = arc(:,jd,kd) * r_0    ! not really necessary as we take fractional arc length
                    arc(:,jd,kd) = grid%theta_F(1,jd,kd) + &
                        &arc(:,jd,kd) / arc(grid%n(1),jd,kd) * &
                        &(grid%theta_F(grid%n(1),jd,kd)-grid%theta_F(1,jd,kd))
                end if
            end do
        end do
    end function calc_arc_angle
end module grid_utilities