sol_utilities.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------!
!> Numerical utilities related to the solution vectors.
!------------------------------------------------------------------------------!
module sol_utilities
#include <PB3D_macros.h>
#include <wrappers.h>
    use str_utilities
    use output_ops
    use messages
    use num_vars, only: dp, iu, max_str_ln, pi
    use grid_vars, only: grid_type
    use eq_vars, only: eq_1_type, eq_2_type
    use x_vars, only: x_1_type, modes_type
    use sol_vars, only: sol_type
 
    implicit none
    private
    public calc_xuq, calc_tot_sol_vec, calc_loc_sol_vec, interp_v_spline
#if ldebug
    public debug_calc_xuq_arr, debug_interp_v_spline
#endif
    
    ! global variables
#if ldebug
    logical :: debug_calc_xuq_arr = .false.                                     !< plot debug information for calc_XUQ_arr \ldebug
    logical :: debug_interp_v_spline = .false.                                  !< debug information for interp_v_spline \ldebug
#endif
    
    ! interfaces
    
    !> \public Calculates the normal (\f$\cdot_n\f$) or geodesic (\f$\cdot_g\f$)
    !! components   of   the   plasma  perturbation   \f$\vec{\xi}\f$   or   the
    !! magnetic  perturbation \f$\vec{Q}  = \nabla  \times \left(\vec{x}  \times
    !! \vec{B}\right)\f$.
    !!
    !! Which variables are plot depends on \c XUQ_style \cite Weyens3D :
    !!  - \c XUQ_style = 1: \f$X\f$ (supports parallel derivative)
    !!  - \c XUQ_style = 2: \f$U\f$ (supports parallel derivative)
    !!  - \c XUQ_style = 3: \f$Q_n\f$
    !!  - \c XUQ_style = 4: \f$Q_g\f$
    !!
    !! \f$X\f$ and \f$U\f$ support the calculation  of the parallel derivative.
    !!
    !! For \f$Q_n\f$ and \f$Q_g\f$, the metric  variables have to be provided as
    !! well.
    !!
    !! The perturbation grid is assumed to  have the same angular coordinates as
    !! the equilibrium grid, and the normal coordinates correspond to either the
    !! equilibrium grid (\c X_grid_style 1),  the solution grid (\c X_grid_style
    !! 2) or the enriched equilibrium grid (\c X_grid_style 3).
    !!
    !! The output grid,  furthermore, has the angular part  corresponding to the
    !! equilibrium grid, and the normal part to the solution grid.
    !!
    !! \return ierr
    interface calc_xuq
        !> \public
        module procedure calc_xuq_ind
        !> \public
        module procedure calc_xuq_arr
    end interface
 
contains
    !> \private (time) array version
    integer function calc_xuq_arr(mds_X,mds_sol,grid_eq,grid_X,grid_sol,eq_1,&
        &eq_2,X,sol,X_id,XUQ_style,time,XUQ,deriv) result(ierr)
        
        use num_vars, only: use_pol_flux_f, norm_disc_prec_sol, x_grid_style
        use num_utilities, only: con2dis, spline
        use x_utilities, only: trim_modes
#if ldebug
        use eq_vars, only: max_flux_f
        use num_utilities, only: calc_int
#endif
        
        character(*), parameter :: rout_name = 'calc_XUQ_arr'
        
        ! input / output
        type(modes_type), intent(in), target :: mds_x                           !< general modes variables in perturbation grid
        type(modes_type), intent(in), target :: mds_sol                         !< general modes variables in solution grid
        type(grid_type), intent(in) :: grid_eq                                  !< equilibrium grid
        type(grid_type), intent(in) :: grid_x                                   !< perturbation grid
        type(grid_type), intent(in) :: grid_sol                                 !< solution grid
        type(eq_1_type), intent(in) :: eq_1                                     !< flux equilibrium
        type(eq_2_type), intent(in) :: eq_2                                     !< metric equilibrium
        type(x_1_type), intent(in), target :: x                                 !< perturbation variables
        type(sol_type), intent(in) :: sol                                       !< solution variables
        integer, intent(in) :: x_id                                             !< nr. of Eigenvalue
        integer, intent(in) :: xuq_style                                        !< whether to calculate \f$X\f$, \f$U\f$, \f$Q_n\f$ or \f$Q_g\f$
        real(dp), intent(in) :: time(:)                                         !< time range
        complex(dp), intent(inout) :: xuq(:,:,:,:)                              !< \f$X\f$, \f$U\f$, \f$Q_n\f$ or \f$Q_g\f$
        logical, intent(in), optional :: deriv                                  !< return parallel derivative
        
        ! local variables
        integer :: m                                                            ! secondary mode number
        integer :: kd_tot                                                       ! kd in total grid
        integer :: ind_x                                                        ! index of mode in X tables
        integer :: ind_sol                                                      ! index of mode in sol tables
        integer :: n_mod_tot                                                    ! total number of modes
        integer :: n_t                                                          ! number of time points
        integer :: n_loc, m_loc                                                 ! local mode numbers
        integer :: id, jd, kd                                                   ! counters
        integer :: xuq_dims(4)                                                  ! dimensions of XUQ
        integer :: kdl_x(2), kdl_sol(2)                                         ! limits on normal index for a mode in X and sol grids
        integer :: id_lim_x(2), id_lim_sol(2)                                   ! limits on total modes
        integer, pointer :: sec_x_loc(:,:), sec_sol_loc(:,:)                    ! pointers to secondary mode variables
        real(dp), allocatable :: expon(:,:,:)                                   ! exponent in solution grid
        real(dp), allocatable :: par_fac(:)                                     ! multiplicative factor due to parallel derivative, without iu, in perturbation grid
        real(dp), allocatable :: jq(:)                                          ! iota or q, interpolated at perturbation grid
        real(dp), allocatable :: s(:,:,:), inv_j(:,:,:)                         ! S and 1/J, interpolated at perturbation grid
        real(dp), pointer :: r_x_loc(:), r_sol_loc(:)                           ! local r_X and r_sol for a mode
        real(dp), pointer :: theta_f(:,:,:), zeta_f(:,:,:)                      ! angles in solution normal grid
        complex(dp) :: sqrt_sol_val_norm                                        ! normalized sqrt(sol_val)
        complex(dp), allocatable :: fac_0(:,:,:), fac_1(:,:,:)                  ! factor to be multiplied with X and DX in perturbation grid
        complex(dp), allocatable :: xuq_loc(:,:)                                ! complex XUQ without time at a normal point
        complex(dp), allocatable :: sol_vec(:,:)                                ! total solution vector in solution grid
        complex(dp), allocatable :: sol_vec_tot(:,:,:)                          ! total solution vector in solution grid and derivatives
        complex(dp), allocatable :: dsol_vec(:,:)                               ! normal derivative of sol_vec in solution grid
        complex(dp), pointer :: u0(:,:,:), u1(:,:,:)                            ! U0 and U1 or parallel derivative
        logical :: extrap = .true.                                              ! whether extrapolation is used
        logical :: deriv_loc                                                    ! local copy of deriv
        character(len=max_str_ln) :: err_msg                                    ! error message
        
#if ldebug
        real(dp), allocatable :: sol_vec_alt(:)                                 ! sol_vec calculated from Dsol_vec in solution grid
#endif
        
        ! initialize ierr
        ierr = 0
        
        ! set up dimensions
        n_t = size(time)
        xuq_dims = shape(xuq)
        
        ! set up local copy of deriv
        deriv_loc = .false.
        if (present(deriv)) deriv_loc = deriv
        
        ! tests
        if (grid_eq%n(1).ne.grid_x%n(1) .or. grid_eq%n(2).ne.grid_x%n(2)) then
            ierr = 1
            err_msg = 'Grids need to be compatible'
            chckerr(err_msg)
        end if
        if (grid_eq%n(1).ne.xuq_dims(1) .or. grid_eq%n(2).ne.xuq_dims(2) .or. &
            &grid_sol%loc_n_r.ne.xuq_dims(3) .or. n_t.ne.xuq_dims(4)) then
            ierr = 1
            err_msg = 'XUQ needs to have the correct dimensions'
            chckerr(err_msg)
        end if
        
        ! limit input modes to output modes
        ! (X grid should comprise sol grid)
        ierr = trim_modes(mds_x,mds_sol,id_lim_x,id_lim_sol)
        chckerr('')
        sec_x_loc => mds_x%sec(id_lim_x(1):id_lim_x(2),:)
        sec_sol_loc => mds_sol%sec(id_lim_sol(1):id_lim_sol(2),:)
        
        n_mod_tot = size(sec_x_loc,1)
        
        if (size(sec_x_loc,1).ne.size(sec_sol_loc,1)) then
            ierr = 1
            err_msg = 'modes X and sol have inconsistent sizes'
            chckerr(err_msg)
        end if
        
        ! set up normalized sqrt(sol_val)
        sqrt_sol_val_norm = sqrt(sol%val(x_id))
        if (ip(sqrt_sol_val_norm).gt.0) sqrt_sol_val_norm = - sqrt_sol_val_norm ! exploding solution, not the decaying one
        sqrt_sol_val_norm = sqrt_sol_val_norm / abs(sqrt_sol_val_norm)          ! normalize
        
        ! allocate helper variables
        allocate(jq(xuq_dims(3)))
        allocate(s(xuq_dims(1),xuq_dims(2),xuq_dims(3)))
        if ((xuq_style.eq.3 .or. xuq_style.eq.4)) &
            &allocate(inv_j(xuq_dims(1),xuq_dims(2),xuq_dims(3)))
        allocate(xuq_loc(xuq_dims(1),xuq_dims(2)))
        allocate(dsol_vec(sol%n_mod,xuq_dims(3)))
        
        ! interpolate eq->sol
        if (use_pol_flux_f) then
            ierr = spline(grid_eq%loc_r_F,eq_1%q_saf_FD(:,0),&
                &grid_sol%loc_r_F,jq,ord=norm_disc_prec_sol)
            chckerr('')
        else
            ierr = spline(grid_eq%loc_r_F,eq_1%rot_t_FD(:,0),&
                &grid_sol%loc_r_F,jq,ord=norm_disc_prec_sol)
            chckerr('')
        end if
        do jd = 1,xuq_dims(2)
            do id = 1,xuq_dims(1)
                ierr = spline(grid_eq%loc_r_F,eq_2%S(id,jd,:),&
                    &grid_sol%loc_r_F,s(id,jd,:),ord=norm_disc_prec_sol)
                chckerr('')
                if ((xuq_style.eq.3 .or. xuq_style.eq.4)) then
                    ierr = spline(grid_eq%loc_r_F,&
                        &1._dp/eq_2%jac_FD(id,jd,:,0,0,0),grid_sol%loc_r_F,&
                        &inv_j(id,jd,:),ord=norm_disc_prec_sol)
                    chckerr('')
                end if
            end do
        end do
        
        ! initialize XUQ
        xuq = 0._dp
        
        ! set up normal derivative of X vec
        ! Note: In X_style 2 (fast), the mode  indices vary as a function as the
        ! normal coordinate.  All of  these contiguous  ranges are  tabulated in
        ! variable 'sec'.
        if (xuq_style.eq.2 .or. xuq_style.eq.4) then
            ! set up variables
            allocate(sol_vec_tot(n_mod_tot,grid_sol%loc_n_r,0:1))
            sol_vec_tot = 0._dp
            
            ! convert to total solution vector and derivate
            do id = 0,1
                ierr = calc_tot_sol_vec(mds_sol,grid_sol,sol%vec(:,:,x_id),&
                    &sol_vec,deriv=id)
                chckerr('')
                sol_vec_tot(:,:,id) = sol_vec
            end do
            
#if ldebug
            if (debug_calc_xuq_arr) then
                do m = 1,n_mod_tot
                    call writo('For mode '//trim(i2str(m))//', testing &
                        &whether Dsol_vec is correct by comparing its integral &
                        &with original sol_vec')
                    kdl_sol = sec_sol_loc(m,2:3)
                    allocate(sol_vec_alt(1:kdl_sol(2)-kdl_sol(1)+1))
                    sol_vec_alt = 0._dp
                    ierr = calc_int(rp(sol_vec_tot(m,kdl_sol(1):kdl_sol(2),1)),&
                        &grid_sol%loc_r_F(kdl_sol(1):kdl_sol(2)),sol_vec_alt)
                    chckerr('')
                    call print_ex_2d(['TEST_RE_Dsol_vec'],'',reshape(&
                        &[rp(sol_vec_tot(m,kdl_sol(1):kdl_sol(2),0)),&
                        &sol_vec_alt],[kdl_sol(2)-kdl_sol(1)+1,2]),x=reshape(&
                        &grid_sol%loc_r_F(kdl_sol(1):kdl_sol(2)),&
                        &[kdl_sol(2)-kdl_sol(1)+1,1])/max_flux_f*2*pi)
                    ierr = calc_int(ip(sol_vec_tot(m,kdl_sol(1):kdl_sol(2),1)),&
                        &grid_sol%loc_r_F(kdl_sol(1):kdl_sol(2)),sol_vec_alt)
                    chckerr('')
                    call print_ex_2d(['TEST_IM_Dsol_vec'],'',reshape(&
                        &[ip(sol_vec_tot(m,kdl_sol(1):kdl_sol(2),0)),&
                        &sol_vec_alt],[kdl_sol(2)-kdl_sol(1)+1,2]),x=reshape(&
                        &grid_sol%loc_r_F(kdl_sol(1):kdl_sol(2)),&
                        &[kdl_sol(2)-kdl_sol(1)+1,1])/max_flux_f*2*pi)
                    deallocate(sol_vec_alt)
                end do
            end if
#endif
            
            ! convert back to local solution vector
            ierr = calc_loc_sol_vec(mds_sol,grid_sol%i_min,sol_vec_tot(:,:,1),&
                &dsol_vec)
            chckerr('')
            
            ! clean up
            deallocate(sol_vec,sol_vec_tot)
        else
            dsol_vec = 0._dp
        end if
        
        ! interpolate or point angles
        select case (x_grid_style)
            case (1,3)                                                          ! equilibrium or enriched
                allocate(theta_f(xuq_dims(1),xuq_dims(2),xuq_dims(3)))
                allocate(zeta_f(xuq_dims(1),xuq_dims(2),xuq_dims(3)))
                do jd = 1,xuq_dims(2)
                    do id = 1,xuq_dims(1)
                        ierr = spline(grid_eq%loc_r_F,grid_eq%theta_F(id,jd,:),&
                            &grid_sol%loc_r_F,theta_f(id,jd,:))
                        chckerr('')
                        ierr = spline(grid_eq%loc_r_F,grid_eq%zeta_F(id,jd,:),&
                            &grid_sol%loc_r_F,zeta_f(id,jd,:))
                        chckerr('')
                    end do
                end do
            case (2)                                                            ! solution
                theta_f => grid_x%theta_F
                zeta_f => grid_x%zeta_F
        end select
        
        ! iterate over all secondary modes
        do m = 1,n_mod_tot
            ! test
            if (sec_x_loc(m,1).ne.sec_sol_loc(m,1)) then
                ierr = 1
                err_msg = 'For mode '//trim(i2str(m))//&
                    &', no consistency for modes in grid_X and grid_sol'
                chckerr(err_msg)
            end if
            
            ! set up and normal limits
            kdl_x = sec_x_loc(m,2:3)
            kdl_sol = sec_sol_loc(m,2:3)
            
            ! limit to grid ranges
            kdl_x(1) = max(kdl_x(1),grid_x%i_min)
            kdl_x(2) = min(kdl_x(2),grid_x%i_max)
            kdl_sol(1) = max(kdl_sol(1),grid_sol%i_min)
            kdl_sol(2) = min(kdl_sol(2),grid_sol%i_max)
            
            ! convert limits to local
            kdl_x = kdl_x - grid_x%i_min + 1
            kdl_sol = kdl_sol - grid_sol%i_min + 1
            
            ! cycle if mode does not exist anywhere
            if (kdl_x(1).gt.kdl_x(2) .or. kdl_sol(1).gt.kdl_sol(2)) cycle
            
            ! get perturbation and solution normal ranges for this mode number
            r_x_loc => grid_x%loc_r_F(kdl_x(1):kdl_x(2))
            r_sol_loc => grid_sol%loc_r_F(kdl_sol(1):kdl_sol(2))
            
            ! set up table indices
            ind_x = sec_x_loc(m,4)
            ind_sol = sec_sol_loc(m,4)
            
            ! Set up multiplicative factors fac_0 and fac_1 which are used in
            ! XUQ = fac_0*X + fac_1*DX with fac_0 and fac_1:
            !   1. fac_0 = 1 (i(nq-m) or i(n-iota m) for deriv)
            !      fac_1 = 0
            !   2. fac_0 = U_0 (DU_0 for deriv)
            !      fac_1 = U_1 (DU_1 for deriv)
            !   3. fac_0 = i(nq-m)/J or i(n-iota m)/J (no deriv)
            !      fac_1 = 0
            !   4. fac_0 = DU_0/J - S (no deriv)
            !      fac_1 = DU_1/J
            ! where par_fac = nq-m or n-iota m.
            ! These quantities are then multiplied by the exponent for each mode
            ! number.
            allocate(expon(xuq_dims(1),xuq_dims(2),kdl_sol(1):kdl_sol(2)))      ! normal subset of grid_sol%loc_r_F
            allocate(fac_0(xuq_dims(1),xuq_dims(2),kdl_sol(1):kdl_sol(2)))      ! normal subset of grid_sol%loc_r_F
            allocate(fac_1(xuq_dims(1),xuq_dims(2),kdl_sol(1):kdl_sol(2)))      ! normal subset of grid_sol%loc_r_F
            allocate(par_fac(kdl_sol(1):kdl_sol(2)))                            ! normal subset of grid_sol%loc_r_F
            
            ! calculate parallel multiplicative factor par_fac and exponent
            do kd = kdl_sol(1),kdl_sol(2)
                ! set up local mode numbers
                kd_tot = kd + grid_sol%i_min - 1
                n_loc = mds_sol%n(kd_tot,ind_sol)
                m_loc = mds_sol%m(kd_tot,ind_sol)
#if ldebug
                if (use_pol_flux_f .and. (m_loc.ne.sec_x_loc(m,1)) .or. &
                    &(.not.use_pol_flux_f .and. (n_loc.ne.sec_x_loc(m,1)))) then
                    ierr = 1
                    err_msg = 'mode numbers not consistent'
                    chckerr(err_msg)
                end if
#endif
                
                if (use_pol_flux_f) then
                    par_fac(kd) = n_loc*jq(kd) - m_loc
                else
                    par_fac(kd) = n_loc - m_loc*jq(kd)
                end if
                expon(:,:,kd) = n_loc*zeta_f(:,:,kd) - &
                    &m_loc*theta_f(:,:,kd)
            end do
            
            ! set up multiplicative factors in solution grid
            select case (xuq_style)
                case (1)                                                        ! calculating X
                    if (deriv_loc) then                                         ! parallel derivative
                        do kd = kdl_sol(1),kdl_sol(2)
                            fac_0(:,:,kd) = iu*par_fac(kd)
                        end do
                    else
                        fac_0 = 1._dp
                    end if
                    fac_1 = 0._dp
                case (2)                                                        ! calculating U
                    ! interpolate X-> sol or point
                    select case (x_grid_style)
                        case (1,3)                                              ! equilibrium or enriched
                            allocate(u0(xuq_dims(1),xuq_dims(2),&
                                &kdl_sol(1):kdl_sol(2)))
                            allocate(u1(xuq_dims(1),xuq_dims(2),&
                                &kdl_sol(1):kdl_sol(2)))
                            if (deriv_loc) then                                 ! parallel derivative
                                ierr = interp_v_spline(&
                                    &x%DU_0(:,:,kdl_x(1):kdl_x(2),ind_x),u0,&
                                    &r_x_loc,r_sol_loc,extrap)
                                chckerr('')
                                ierr = interp_v_spline(&
                                    &x%DU_1(:,:,kdl_x(1):kdl_x(2),ind_x),u1,&
                                    &r_x_loc,r_sol_loc,extrap)
                                chckerr('')
                            else
                                ierr = interp_v_spline(&
                                    &x%U_0(:,:,kdl_x(1):kdl_x(2),ind_x),u0,&
                                    &r_x_loc,r_sol_loc,extrap)
                                chckerr('')
                                ierr = interp_v_spline(&
                                    &x%U_1(:,:,kdl_x(1):kdl_x(2),ind_x),u1,&
                                    &r_x_loc,r_sol_loc,extrap)
                                chckerr('')
                            end if
                        case (2)                                                ! solution
                            if (deriv_loc) then                                 ! parallel derivative
                                u0 => x%DU_0(:,:,kdl_x(1):kdl_x(2),ind_x)
                                u1 => x%DU_1(:,:,kdl_x(1):kdl_x(2),ind_x)
                            else
                                u0 => x%U_0(:,:,kdl_x(1):kdl_x(2),ind_x)
                                u1 => x%U_1(:,:,kdl_x(1):kdl_x(2),ind_x)
                            end if
                    end select
                    
                    ! assign
                    fac_0 = u0
                    fac_1 = u1
                    
                    ! clean up
                    select case (x_grid_style)
                        case (1,3)                                              ! equilibrium or enriched
                            deallocate(u0,u1)
                        case (2)                                                ! solution
                            ! do nothing
                    end select
                case (3)                                                        ! calculating Qn
                    if (deriv_loc) then                                         ! parallel derivative
                        ierr = 1
                        err_msg = 'For XUQ_style = '//trim(i2str(xuq_style))//&
                            &' calculation of parallel derivative not supported'
                        chckerr('')
                    else
                        do kd = kdl_sol(1),kdl_sol(2)
                            fac_0(:,:,kd) = iu*par_fac(kd)*inv_j(:,:,kd)
                        end do
                    end if
                    fac_1 = 0._dp
                case (4)                                                        ! calculating Qg
                    ! interpolate X-> sol or point
                    err_msg = 'For XUQ_style = '//trim(i2str(xuq_style))//&
                        &' calculation of parallel derivative not supported'
                    chckerr('')
                    select case (x_grid_style)
                        case (1,3)                                              ! equilibrium or enriched
                            allocate(u0(xuq_dims(1),xuq_dims(2),&
                                &kdl_sol(1):kdl_sol(2)))
                            allocate(u1(xuq_dims(1),xuq_dims(2),&
                                &kdl_sol(1):kdl_sol(2)))
                            if (deriv_loc) then                                 ! parallel derivative
                                ierr = 1
                                chckerr(err_msg)
                            else
                                ierr = interp_v_spline(&
                                    &x%DU_0(:,:,kdl_x(1):kdl_x(2),ind_x),u0,&
                                    &r_x_loc,r_sol_loc,extrap)
                                chckerr('')
                                ierr = interp_v_spline(&
                                    &x%DU_1(:,:,kdl_x(1):kdl_x(2),ind_x),u1,&
                                    &r_x_loc,r_sol_loc,extrap)
                                chckerr('')
                            end if
                        case (2)                                                ! solution
                            if (deriv_loc) then                                 ! parallel derivative
                                ierr = 1
                                chckerr(err_msg)
                            else
                                u0 => x%DU_0(:,:,kdl_x(1):kdl_x(2),ind_x)
                                u1 => x%DU_1(:,:,kdl_x(1):kdl_x(2),ind_x)
                            end if
                    end select
                    
                    ! assign
                    fac_0 = u0*inv_j(:,:,kdl_sol(1):kdl_sol(2)) - &
                        &s(:,:,kdl_sol(1):kdl_sol(2))
                    fac_1 = u1*inv_j(:,:,kdl_sol(1):kdl_sol(2))
                    
                    ! clean up
                    select case (x_grid_style)
                        case (1,3)                                              ! equilibrium or enriched
                            deallocate(u0,u1)
                        case (2)                                                ! solution
                            ! do nothing
                    end select
                case default
                    err_msg = 'No XUQ style associated with '//&
                        &trim(i2str(xuq_style))
                    ierr = 1
                    chckerr(err_msg)
            end select
            
            ! iterate over all normal points in solution modes variables
            do kd = kdl_sol(1),kdl_sol(2)
                ! set up loc complex XUQ without time at this normal point
                xuq_loc = exp(iu*expon(:,:,kd)) * &
                    &(sol%vec(ind_sol,kd,x_id)*fac_0(:,:,kd) + &
                    &dsol_vec(ind_sol,kd)*fac_1(:,:,kd))
                
                ! iterate over time steps
                do id = 1,n_t
                    ! add current mode
                    xuq(:,:,kd,id) = xuq(:,:,kd,id) + xuq_loc * &
                        &exp(iu*sqrt_sol_val_norm*time(id)*2*pi)
                end do
            end do
            
            ! clean up
            deallocate(fac_0,fac_1,expon,par_fac)
        end do
        
        ! clean up
        select case (x_grid_style)
            case (1,3)                                                          ! equilibrium or enriched
                deallocate(theta_f,zeta_f)
            case (2)                                                            ! solution
                ! do nothing
        end select
        nullify(r_x_loc,r_sol_loc)
        nullify(theta_f,zeta_f)
        nullify(u0,u1)
        nullify(sec_x_loc,sec_sol_loc)
    end function calc_xuq_arr
    !> \private (time) individual version
    integer function calc_xuq_ind(mds_X,mds_sol,grid_eq,grid_X,grid_sol,eq_1,&
        &eq_2,X,sol,X_id,XUQ_style,time,XUQ,deriv) result(ierr)
        
        character(*), parameter :: rout_name = 'calc_XUQ_ind'
        
        ! input / output
        type(modes_type), intent(in) :: mds_x                                   !< general modes variables in perturbation grid
        type(modes_type), intent(in) :: mds_sol                                 !< general modes variables in solution grid
        type(grid_type), intent(in) :: grid_eq                                  !< equilibrium grid
        type(grid_type), intent(in) :: grid_x                                   !< perturbation grid
        type(grid_type), intent(in) :: grid_sol                                 !< solution grid
        type(eq_1_type), intent(in) :: eq_1                                     !< flux equilibrium
        type(eq_2_type), intent(in) :: eq_2                                     !< metric equilibrium
        type(x_1_type), intent(in) :: x                                         !< perturbation variables
        type(sol_type), intent(in) :: sol                                       !< solution variables
        integer, intent(in) :: x_id                                             !< nr. of Eigenvalue
        integer, intent(in) :: xuq_style                                        !< whether to calculate X, U, Qn or Qg
        real(dp), intent(in) :: time                                            !< time
        complex(dp), intent(inout) :: xuq(:,:,:)                                !< normal component of perturbation
        logical, intent(in), optional :: deriv                                  !< return parallel derivative
        
        ! local variables
        complex(dp), allocatable :: xuq_arr(:,:,:,:)                            ! dummy variable to use array version
        
        ! allocate array version of XUQ
        allocate(xuq_arr(size(xuq,1),size(xuq,2),size(xuq,3),1))
        
        ! call array version
        ierr = calc_xuq_arr(mds_x,mds_sol,grid_eq,grid_x,grid_sol,eq_1,eq_2,x,&
            &sol,x_id,xuq_style,[time],xuq_arr,deriv=deriv)
        chckerr('')
        
        ! copy array to individual XUQ
        xuq = xuq_arr(:,:,:,1)
        
        ! deallocate array version of XUQ
        deallocate(xuq_arr)
    end function calc_xuq_ind
    
    !> Calculate  the  total version  of  the  solution  vector from  the  local
    !! version.
    !!
    !! This is the  solution vector for all of the  possible mode numbers, which
    !! can be different from the local mode numbers for X style 2 (fast):
    !!  -  local: the  size of  the saved  perturbation and  solution variables,
    !!  prescribed  by the  user as  input:
    !!      - either directly through \c n_mod_X
    !!      -  or through  the limits  on the  modes using  \c min_sec_X  and \c
    !!      max_sec_X.
    !!  - total: all  the possible modes that can resonate  in the plasma, which
    !!  can be  different from the  local number for \c  X_style 2, as  the mode
    !!  numbers depend on the normal coordinate.
    !!
    !! If the output variable is not allocated, it is done here.
    !!
    !! Optionally  a  derivative  can  be requested,  depending  on  the  normal
    !! discretization order.
    !!
    !! \note Need to provide the grid in which mds is tabulated.
    !!
    !! \return ierr
    integer function calc_tot_sol_vec(mds,grid_sol,sol_vec_loc,sol_vec_tot,&
        &deriv) result(ierr)
        use num_vars, only: x_style, norm_disc_prec_sol
        use x_vars, only: n_mod_x
        use num_utilities, only: spline
        
        character(*), parameter :: rout_name = 'calc_tot_sol_vec'
        
        ! input / output
        type(modes_type), intent(in) :: mds                                     !< general modes variables
        type(grid_type), intent(in) :: grid_sol                                 !< solution grid
        complex(dp), intent(in) :: sol_vec_loc(:,:)                             !< solution vector for local resonating modes
        complex(dp), intent(inout), allocatable :: sol_vec_tot(:,:)             !< solution vector for all possible resonating modes
        integer, intent(in), optional :: deriv                                  !< order of derivative
        
        ! local variables
        character(len=max_str_ln) :: err_msg                                    ! error message
        integer :: n_mod_loc                                                    ! local number of modes that are used in the calculations
        integer :: n_mod_tot                                                    ! total number of modes
        integer :: ld                                                           ! counter
        integer :: kdl(2)                                                       ! kd limits
        integer :: ld_loc                                                       ! local ld
        integer :: n_r                                                          ! normal size of variables
        integer :: deriv_loc                                                    ! local derivative
        integer :: max_deriv                                                    ! maximum derivative
        integer :: min_range                                                    ! minimum range for derivative
        complex(dp), allocatable :: dsol_vec_loc(:)                             ! local Dsol_vec_tot in solution grid
        
        ! initialize ierr
        ierr = 0
        
        ! set n_r
        n_r = size(sol_vec_loc,2)
        
        ! set local derivative
        deriv_loc = 0
        if (present(deriv)) deriv_loc = deriv
        
        ! test
        select case (norm_disc_prec_sol)
            case (1)                                                            ! linear
                max_deriv = 1
                min_range = 2
            case (2)                                                            ! Akita hermite
                max_deriv = 1
                min_range = 4
            case (3)                                                            ! spline
                max_deriv = 3
                min_range = 4
            case default
                ierr = 1
                err_msg = 'normal discretization precision for solution &
                    &cannot be higher than 3 or lower than 0'
                chckerr(err_msg)
        end select
        if (deriv_loc.lt.0 .or. deriv_loc.gt.max_deriv) then
            ierr = 1
            err_msg = 'deriv must be between 0 and '//trim(i2str(max_deriv))//&
                &' for normal discretization '//trim(i2str(norm_disc_prec_sol))
            chckerr(err_msg)
        endif
        
        ! operations depending on X style
        select case (x_style)
            case (1)                                                            ! prescribed
                ! test allocation
                if (allocated(sol_vec_tot)) then
                    if (size(sol_vec_tot,1).ne.size(sol_vec_loc,1) .or. &
                        &size(sol_vec_tot,2).ne.size(sol_vec_loc,2)) then
                        ierr = 1
                        err_msg = 'Local and total solution vectors need to &
                            &be of the same length'
                        chckerr(err_msg)
                    end if
                else
                    allocate(sol_vec_tot(size(sol_vec_loc,1),n_r))
                end if
                
                ! copy or derivate
                if (deriv_loc.eq.0) then                                        ! copy
                    sol_vec_tot = sol_vec_loc
                else                                                            ! derivate
                    do ld = 1,size(sol_vec_tot,1)
                        ierr = spline(grid_sol%r_F,sol_vec_loc(ld,:),&
                            &grid_sol%r_F,sol_vec_tot(ld,:),&
                            &ord=norm_disc_prec_sol,deriv=deriv_loc)
                        chckerr('')
                    end do
                end if
            case (2)                                                            ! fast
                ! set local and total nr. of modes
                n_mod_loc = n_mod_x
                n_mod_tot = size(mds%sec,1)
                
                ! test allocation
                if (allocated(sol_vec_tot)) then
                    if (size(sol_vec_tot,1).ne.n_mod_tot .or. &
                        &size(sol_vec_tot,2).ne.size(sol_vec_loc,2)) then
                        ierr = 1
                        err_msg = 'Total solution vectors needs to be of &
                            &correct length '//trim(i2str(n_mod_tot))
                        chckerr(err_msg)
                    end if
                else
                    allocate(sol_vec_tot(n_mod_tot,n_r))
                end if
                
                ! track all possible modes and copy or derivate
                sol_vec_tot = 0._dp
                do ld = 1,size(mds%sec,1)
                    ld_loc = mds%sec(ld,1)-minval(mds%sec(:,1),1)+1
                    kdl = [mds%sec(ld,2),mds%sec(ld,3)]-grid_sol%i_min+1        ! normal range in solution vector
                    kdl = max(1,min(kdl,n_r))                                   ! limit to sol_vec_loc range
                    if (deriv_loc.eq.0) then
                        sol_vec_tot(ld_loc,kdl(1):kdl(2)) = &
                            &sol_vec_loc(mds%sec(ld,4),kdl(1):kdl(2))
                    else
                        allocate(dsol_vec_loc(kdl(2)-kdl(1)+1))
                        dsol_vec_loc = 0._dp
                        
                        if (kdl(2)-kdl(1)+1.ge.min_range) then                  ! only calculate if enough points
                            ierr = spline(grid_sol%r_F(kdl(1):kdl(2)),&
                                &sol_vec_loc(mds%sec(ld,4),kdl(1):kdl(2)),&
                                &grid_sol%r_F(kdl(1):kdl(2)),dsol_vec_loc,&
                                &ord=norm_disc_prec_sol,deriv=deriv_loc)
                            chckerr('')
                        end if
                        sol_vec_tot(ld_loc,kdl(1):kdl(2)) = dsol_vec_loc
                        
                        deallocate(dsol_vec_loc)
                    end if
                end do
        end select
    end function calc_tot_sol_vec
    
    !> Calculate  the  local version  of  the  solution  vector from  the  total
    !! version.
    !!
    !! \see See calc_tot_sol_vec().
    !!
    !! \return ierr
    integer function calc_loc_sol_vec(mds,i_min,sol_vec_tot,sol_vec_loc) &
        &result(ierr)
        use num_vars, only: x_style
        use x_vars, only: n_mod_x
        
        character(*), parameter :: rout_name = 'calc_loc_sol_vec'
        
        ! input / output
        type(modes_type), intent(in) :: mds                                     !< general modes variables
        integer, intent(in) :: i_min                                            !< \c i_min of grid in which variables are tabulated
        complex(dp), intent(in) :: sol_vec_tot(:,:)                             !< solution vector for all possible resonating modes
        complex(dp), intent(inout), allocatable :: sol_vec_loc(:,:)             !< solution vector for local resonating modes
        
        ! local variables
        character(len=max_str_ln) :: err_msg                                    ! error message
        integer :: n_mod_loc                                                    ! local number of modes that are used in the calculations
        integer :: n_mod_tot                                                    ! total number of modes
        integer :: ld                                                           ! counter
        integer :: kdl(2)                                                       ! kd limits
        integer :: ld_loc                                                       ! local ld
        integer :: n_r                                                          ! normal size of variables
        
        ! set n_r
        n_r = size(sol_vec_tot,2)
        
        ! initialize ierr
        ierr = 0
        
        ! operations depending on X style
        select case (x_style)
            case (1)                                                            ! prescribed
                ! test allocation
                if (allocated(sol_vec_loc)) then
                    if (size(sol_vec_loc,1).ne.size(sol_vec_tot,1) .or. &
                        &size(sol_vec_loc,2).ne.size(sol_vec_tot,2)) then
                        ierr = 1
                        err_msg = 'Local and total solution vectors need to &
                            &be of the same length'
                        chckerr(err_msg)
                    end if
                else
                    allocate(sol_vec_loc(size(sol_vec_tot,1),n_r))
                end if
                sol_vec_loc = sol_vec_tot
            case (2)                                                            ! fast
                ! set local and total nr. of modes
                n_mod_loc = n_mod_x
                n_mod_tot = size(mds%sec,1)
                
                ! test allocation
                if (allocated(sol_vec_loc)) then
                    if (size(sol_vec_loc,1).ne.n_mod_loc .or. &
                        &size(sol_vec_loc,2).ne.size(sol_vec_tot,2)) then
                        ierr = 1
                        err_msg = 'Total solution vectors needs to be of &
                            &correct length '//trim(i2str(n_mod_loc))
                        chckerr(err_msg)
                    end if
                else
                    allocate(sol_vec_loc(n_mod_loc,n_r))
                end if
                
                ! track all possible modes and copy
                sol_vec_loc = 0._dp
                do ld = 1,size(mds%sec,1)
                    ld_loc = mds%sec(ld,1)-minval(mds%sec(:,1),1)+1
                    kdl = [mds%sec(ld,2),mds%sec(ld,3)]-i_min+1                 ! range in sol_vec_loc
                    kdl = max(1,min(kdl,n_r))                                   ! limit to sol_vec_loc range
                    sol_vec_loc(mds%sec(ld,4),kdl(1):kdl(2)) = &
                        &sol_vec_tot(ld_loc,kdl(1):kdl(2))
                end do
        end select
    end function calc_loc_sol_vec
    
    !> Interpolation for a quantity V using splines.
    !!
    !! Optionally, a variable 'ivs_stat' is  returned that indicates whether the
    !! interpolation  was a  plain copy  (1), a  linear interpolation  (2) or  a
    !! spline interpolation (3).
    integer function interp_v_spline(V_i,V_o,r_i,r_o,extrap,ivs_stat) &
        &result(ierr)
        
        use num_utilities, only: spline
        
        character(*), parameter :: rout_name = 'interp_V_spline'
        
        ! input / output
        complex(dp), intent(in) :: v_i(:,:,:)                                   ! input variable
        complex(dp), intent(out) :: v_o(:,:,:)                                  ! output variable
        real(dp), intent(in) :: r_i(:)                                          ! input r
        real(dp), intent(in) :: r_o(:)                                          ! output r
        logical, intent(in) :: extrap                                           ! extrapolation possible
        integer, intent(out), optional :: ivs_stat                              ! which method was used (1: copy, 2: linear, 3: spline)
        
        ! local variables
        integer :: id, jd, kd                                                   ! counters
        character(len=max_str_ln) :: err_msg                                    ! error message
        
        ! initialize ierr
        ierr = 0
        
#if ldebug
        ! test sizes
        if (size(r_i).ne.size(v_i,3)) then
            ierr = 1
            err_msg = 'r_i and V_i are not compatible'
            chckerr(err_msg)
        end if
        if (size(r_o).ne.size(v_o,3)) then
            ierr = 1
            err_msg = 'r_o and V_o are not compatible'
            chckerr(err_msg)
        end if
        if (size(v_i,1).ne.size(v_o,1) .or. size(v_i,2).ne.size(v_o,2)) then
            ierr = 1
            err_msg = 'V_i and V_o are not compatible'
            chckerr(err_msg)
        end if
#endif
        
        ! interpolate depending on normal size
        select case (size(r_i))
            case (1)                                                            ! copy directly
                do kd = 1,size(v_o,3)
                    v_o(:,:,kd) = v_i(:,:,1)
                end do
                if (present(ivs_stat)) ivs_stat = 1
            case (2)                                                            ! manual linear (EZspline needs >= 3 pts)
                do kd = 1,size(v_o,3)
                    v_o(:,:,kd) = v_i(:,:,1) + &
                        &(v_i(:,:,2) - v_i(:,:,1)) / (r_i(2) - r_i(1)) * &
                        &(r_o(kd) - r_i(1))
                end do
                if (present(ivs_stat)) ivs_stat = 2
            case (3)                                                            ! linear via EZspline
                do jd = 1,size(v_i,2)
                    do id = 1,size(v_i,1)
                        ierr = spline(r_i,v_i(id,jd,:),r_o,v_o(id,jd,:),&
                            &ord=1,extrap=extrap)
                        chckerr('')
                    end do
                end do
                if (present(ivs_stat)) ivs_stat = 2
            case (4:)                                                           ! spline
                do jd = 1,size(v_i,2)
                    do id = 1,size(v_i,1)
                        ierr = spline(r_i,v_i(id,jd,:),r_o,v_o(id,jd,:),&
                            &ord=3,extrap=extrap)
                        chckerr('')
                    end do
                end do
                if (present(ivs_stat)) ivs_stat = 3
#if ldebug
            case default
                ierr = 1
                chckerr('Need more than 0 points')
#endif
        end select
        
#if ldebug
        ! visualize
        do jd = 1,size(v_i,2)
            do id = 1,size(v_i,1)
                if (debug_interp_v_spline) then
                    call plot_comp('real part',r_i,rp(v_i(id,jd,:)),&
                        &r_o,rp(v_o(id,jd,:)))
                    call plot_comp('imaginary part',r_i,ip(v_i(id,jd,:)),&
                        &r_o,ip(v_o(id,jd,:)))
                end if
            end do
        end do
    contains
        !> \private Plot to compare interpolation to original variables
        subroutine plot_comp(plot_name,x_i,y_i,x_o,y_o)
            ! input / output
            character(len=*), intent(in) :: plot_name                           ! name of plot
            real(dp), intent(in) :: x_i(:)                                      ! x input
            real(dp), intent(in) :: y_i(:)                                      ! y input
            real(dp), intent(in) :: x_o(:)                                      ! x output
            real(dp), intent(in) :: y_o(:)                                      ! y output
            
            ! local variables
            real(dp), allocatable :: x_plot(:,:)                                ! x for plot
            real(dp), allocatable :: y_plot(:,:)                                ! y for plot
            character(len=max_str_ln) :: var_names(2)                           ! variable names
            
            allocate(x_plot(max(size(x_i),size(x_o)),2))
            allocate(y_plot(max(size(x_i),size(x_o)),2))
            
            x_plot(1:size(x_i),1) = x_i
            x_plot(size(x_i):size(x_plot,1),1) = x_i(size(x_i))
            x_plot(1:size(x_o),2) = x_o
            x_plot(size(x_o):size(x_plot,1),2) = x_o(size(x_o))
            
            y_plot(1:size(x_i),1) = y_i(:)
            y_plot(size(x_i):size(y_plot,1),1) = y_i(size(x_i))
            y_plot(1:size(x_o),2) = y_o(:)
            y_plot(size(x_o):size(y_plot,1),2) = y_o(size(x_o))
            
            var_names(1) = 'input '//plot_name
            var_names(2) = 'output '//plot_name
            call print_ex_2d(var_names,'',y_plot,x=x_plot)
        end subroutine plot_comp
#endif
    end function interp_v_spline
end module sol_utilities