eq_utilities.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------!
!> Numerical utilities related to equilibrium variables.
!------------------------------------------------------------------------------!
module eq_utilities
#include <PB3D_macros.h>
    use str_utilities
    use output_ops
    use messages
    use num_vars, only: pi, dp, max_str_ln, max_deriv
    use grid_vars, only: grid_type
    use eq_vars, only: eq_1_type, eq_2_type
#if ldebug
    use num_utilities, only: check_deriv
#endif
    
    implicit none
    private
    public calc_inv_met, calc_g, transf_deriv, calc_f_derivs, do_eq, eq_info, &
        &print_info_eq, calc_memory_eq
#if ldebug
    public debug_calc_inv_met_ind
#endif
    
    ! global variables
#if ldebug
    !> \ldebug
    logical :: debug_calc_inv_met_ind = .false.                                 !< plot debug information for calc_inv_met_ind()
#endif
    
    ! interfaces
    
    !> \public  Transforms  derivatives  of  the  equilibrium  quantities  in  E
    !! coordinates to derivatives in the F coordinates.
    !!
    !! \return ierr
    interface calc_f_derivs
        !> \public
        module procedure calc_f_derivs_1
        !> \public
        module procedure calc_f_derivs_2
    end interface
    
    !> \public Calculate
    !!      \f$D_1^{m_1} D_2^{m_2} D_3^{m_3} X\f$
    !! from
    !!      \f$D_1^{i_1} D_2^{i_2} D_3^{i_3} X\f$
    !! and
    !!      \f$D_1^{j_1} D_2^{j_2} D_3^{j_3} Y\f$
    !! where \f$XY=1\f$ and \f$i,j = 0\ldots m\f$, according to \cite Weyens3D.
    !!
    !! 
    !! \f$D_1^{m_1} D_2^{m_2} D_3^{m_3} \f$ is defined as
    !!      \f$\left(\frac{\partial}{\partial u^1}\right)^{m_1}
    !!      \left(\frac{\partial}{\partial u^2}\right)^{m_2}
    !!      \left(\frac{\partial}{\partial u^3}\right)^{m_3} \f$
    !!
    !! \note It is assumed that the lower order derivatives have been calculated
    !! already. If not, the results will be incorrect.
    !!
    !! \return ierr
    interface calc_inv_met
        !> \public
        module procedure calc_inv_met_ind
        !> \public
        module procedure calc_inv_met_arr
        !> \public
        module procedure calc_inv_met_ind_0d
        !> \public
        module procedure calc_inv_met_arr_0d
    end interface
    
    !> \public Calculates derivatives in a  coordinate system B from derivatives
    !! in  a coordinates  system  A,  making use  of  the transformation  matrix
    !! \f$\overline{\text{T}}_\text{B}^\text{A}\f$.
    !!
    !! The routine works by exchanging the  derivatives in the coordinates B for
    !! derivatives in coordinates A using the formula
    !!  \f[\mathbf{D}_\text{B}^m X = \overline{\text{T}}_\text{B}^\text{A}
    !!      \mathbf{D}^1_\text{A} \left(\mathbf{D}^{m-1}_\text{B} X\right) , \f]
    !! where \f$\mathbf{D}^m\f$  is a tensor  of rank \f$m\f$ that  contains all
    !! the derivatives of total rank \f$m\f$. For example,
    !!  \f[\mathbf{D}^1 = \vec{D} =
    !!      \left(\begin{array}{c}\frac{\partial}{\partial u^1} \\
    !!      \frac{\partial}{\partial u^2} \\ \frac{\partial}{\partial u^3}
    !!      \end{array}\right)  \f]
    !! 
    !! This  is done  for the  derivatives in  each of  the coordinates  B until
    !! degree 0 is reached.
    !!
    !! Furthermore, each of these degrees of derivatives in coordinates B can be
    !! be derived optionally  in the original coordinate system  A, which yields
    !! the formula:
    !!  \f[ \mathbf{D}^p_\text{A} \left(\mathbf{D}^m_\text{B} X \right) =
    !!      \sum_q \left(\begin{array}{c}p\\q\end{array}\right)
    !!      \mathbf{D}^q_\text{A}
    !!      \left(\overline{\text{T}}_\text{B}^\text{A}\right)
    !!      \left(\mathbf{D}^{p-q}_\text{A} \mathbf{D}^1_\text{A}\right)
    !!      \left( \mathbf{D}^{m-1}_\text{B} X \right)\f]
    !!
    !! This way, ultimately the desired derivatives  in the coordinates B can be
    !! obtained  recursively from  the lower  orders  in the  coordinates B  and
    !! higher orders in the coordinates A.
    !!
    !! For example:
    !!  - For order \f$M\f$, the formula can be used with \f$p=0\f$.
    !!  -   For   this,   it    is   necessary   that   \f$\mathbf{D}^1_\text{A}
    !!  \mathbf{D}^{m-1}_\text{B} X\f$ be precalculated, which can be done using
    !!  the formula with \f$p=1\f$ and \f$m=M-1\f$.
    !!  -   For   this,   it    is   necessary   that   \f$\mathbf{D}^1_\text{A}
    !!  \mathbf{D}^{M-1}_\text{B}     X\f$      and     \f$\mathbf{D}^2_\text{A}
    !!  \mathbf{D}^{M-1}_\text{B} X\f$ are precalculated, which can be done with
    !!  \f$p=1\f$ and \f$m=M-1\f$ as well as \f$p=2\f$ and \f$m=M-1\f$.
    !!  - etc.
    !!
    !! \see \cite Weyens3D for more detailed information.
    !!
    !! \return ierr
    interface transf_deriv
        !> \public
        module procedure transf_deriv_3_ind
        !> \public
        module procedure transf_deriv_3_arr
        !> \public
        module procedure transf_deriv_3_arr_2d
        !> \public
        module procedure transf_deriv_1_ind
    end interface
    
contains
    !> \private individual matrix version
    integer function calc_inv_met_ind(X,Y,deriv) result(ierr)
        use num_utilities, only: calc_inv, calc_mult, c, conv_mat
#if ldebug
        use num_vars, only: n_procs
#endif
        
        character(*), parameter :: rout_name = 'calc_inv_met_ind'
        
        ! input / output
        real(dp), intent(inout) :: x(1:,1:,1:,1:,0:,0:,0:)                      !< X
        real(dp), intent(in) :: y(1:,1:,1:,1:,0:,0:,0:)                         !< Y
        integer, intent(in) :: deriv(:)                                         !< derivatives
        
        ! local variables
        integer :: r, t, z                                                      ! counters for derivatives
        integer :: m1, m2, m3                                                   ! alias for deriv(i)
        integer :: dims(3)                                                      ! dimensions of coordinates
        real(dp) :: bin_fac(3)                                                  ! binomial factor for norm., pol. and tor. sum
        real(dp), allocatable :: dum1(:,:,:,:)                                  ! dummy variable representing metric matrix for some derivatives
        real(dp), allocatable :: dum2(:,:,:,:)                                  ! dummy variable representing metric matrix for some derivatives
        integer :: nn                                                           ! size of matrices
        character(len=max_str_ln) :: err_msg                                    ! error message
#if ldebug
        real(dp), allocatable :: xy(:,:,:,:)                                    ! to check if XY is indeed unity
#endif
        
        ! initialize ierr
        ierr = 0
        
        ! set nn
        nn = size(x,4)
        
        ! tests
        if (size(y,4).ne.nn) then
            ierr = 1
            err_msg = 'Both matrices X and Y have to be either in full or &
                &lower diagonal storage'
            chckerr(err_msg)
        end if
        
        ! set mi
        m1 = deriv(1)
        m2 = deriv(2)
        m3 = deriv(3)
        
        ! set dims
        dims = [size(x,1),size(x,2),size(x,3)]
        
        ! initialize the requested derivative to zero
        x(:,:,:,:,m1,m2,m3) = 0._dp
        
        ! set up dummy variables
        allocate(dum1(dims(1),dims(2),dims(3),9))
        allocate(dum2(dims(1),dims(2),dims(3),9))
        
        ! initialize dum2
        dum2 = 0._dp
        
        ! calculate terms
        if (m1.eq.0 .and. m2.eq.0 .and. m3.eq.0) then                           ! direct inverse
            ierr = calc_inv(x(:,:,:,:,0,0,0),y(:,:,:,:,0,0,0),3)
            chckerr('')
            
#if ldebug
            ! check if X*Y is unity indeed if debugging
            if (debug_calc_inv_met_ind) then
                if (n_procs.gt.1) call writo('Debugging of calc_inv_met_ind &
                    &should be done with one processor only',warning=.true.)
                call writo('checking if X*Y = 1')
                call lvl_ud(1)
                allocate(xy(dims(1),dims(2),dims(3),9))
                ierr = calc_mult(x(:,:,:,:,0,0,0),y(:,:,:,:,0,0,0),xy,3)
                chckerr('')
                do r = 1,3
                    do t = 1,3
                        call writo(&
                            &trim(r2strt(minval(xy(:,1,:,c([r,t],.false.)))))//&
                            &' < element ('//trim(i2str(t))//','//&
                            &trim(i2str(r))//') < '//&
                            &trim(r2strt(maxval(xy(:,1,:,c([r,t],.false.))))))
                    end do
                end do
                deallocate(xy)
                call lvl_ud(-1)
            end if
#endif
        else                                                                    ! calculate using the other inverses
            do z = 0,m3                                                         ! derivs in third coord
                if (z.eq.0) then                                                ! first term in sum
                    bin_fac(3) = 1.0_dp
                else
                    bin_fac(3) = bin_fac(3)*(m3-(z-1))/z
                end if
                do t = 0,m2                                                     ! derivs in second coord
                    if (t.eq.0) then                                            ! first term in sum
                        bin_fac(2) = 1.0_dp
                    else
                        bin_fac(2) = bin_fac(2)*(m2-(t-1))/t
                    end if
                    do r = 0,m1                                                 ! derivs in first coord
                        if (r.eq.0) then                                        ! first term in sum
                            bin_fac(1) = 1.0_dp
                        else
                            bin_fac(1) = bin_fac(1)*(m1-(r-1))/r
                        end if
                        if (z+t+r .lt. m1+m2+m3) then                           ! only add if not all r,t,z are equal to m1,m2,m3
                            ! multiply X(r,t,z) and Y(m1-r,m2-t,m3-z)
                            ierr = calc_mult(x(:,:,:,:,r,t,z),&
                                &y(:,:,:,:,m1-r,m2-t,m3-z),dum1,3)
                            chckerr('')
                            ! add result, multiplied by bin_fac, to dum2
                            dum2 = dum2 - product(bin_fac)*dum1
                        end if
                    end do
                end do
            end do
            
            ! right-multiply by X(0,0,0) and save in dum1
            ierr = calc_mult(dum2,x(:,:,:,:,0,0,0),dum1,3)
            chckerr('')
            ! save result in X(m1,m2,m3), converting if necessary
            ierr = conv_mat(dum1,x(:,:,:,:,m1,m2,m3),3)
            chckerr('')
            
            ! deallocate local variables
            deallocate(dum1,dum2)
        end if
    end function calc_inv_met_ind
    !> \private individual scalar version
    integer function calc_inv_met_ind_0d(X,Y,deriv) result(ierr)
        character(*), parameter :: rout_name = 'calc_inv_met_ind_0D'
        
        ! input / output
        real(dp), intent(inout) :: x(1:,1:,1:,0:,0:,0:)                         !< X
        real(dp), intent(in) :: y(1:,1:,1:,0:,0:,0:)                            !< Y
        integer, intent(in) :: deriv(:)                                         !< derivatives
        
        ! local variables
        integer :: r, t, z                                                      ! counters for derivatives
        integer :: m1, m2, m3                                                   ! alias for deriv(i)
        real(dp) :: bin_fac(3)                                                  ! binomial factor for norm., pol. and tor. sum
        
        ! initialize ierr
        ierr = 0
        
        m1 = deriv(1)
        m2 = deriv(2)
        m3 = deriv(3)
        
        ! initialize the requested derivative
        x(:,:,:,m1,m2,m3) = 0.0_dp
        
        ! calculate terms
        if (m1.eq.0 .and. m2.eq.0 .and. m3.eq.0) then                           ! direct inverse
            x(:,:,:,0,0,0) = 1._dp/y(:,:,:,0,0,0)
        else                                                                    ! calculate using the other inverses
            do z = 0,m3                                                         ! derivs in third coord
                if (z.eq.0) then                                                ! first term in sum
                    bin_fac(3) = 1.0_dp
                else
                    bin_fac(3) = bin_fac(3)*(m3-(z-1))/z
                    ! ADD SOME ERROR CHECKING HERE !!!!
                    chckerr('')
                end if
                do t = 0,m2                                                     ! derivs in second coord
                    if (t.eq.0) then                                            ! first term in sum
                        bin_fac(2) = 1.0_dp
                    else
                        bin_fac(2) = bin_fac(2)*(m2-(t-1))/t
                    end if
                    do r = 0,m1                                                 ! derivs in first coord
                        if (r.eq.0) then                                        ! first term in sum
                            bin_fac(1) = 1.0_dp
                        else
                            bin_fac(1) = bin_fac(1)*(m1-(r-1))/r
                        end if
                        if (z+t+r .lt. m1+m2+m3) then                           ! only add if not all r,t,z are equal to m1,m2,m3
                            x(:,:,:,m1,m2,m3) = x(:,:,:,m1,m2,m3)-&
                                &bin_fac(1)*bin_fac(2)*bin_fac(3)* &
                                &x(:,:,:,r,t,z)*y(:,:,:,m1-r,m2-t,m3-z)
                        end if
                    end do
                end do
            end do
            
            ! right-multiply by X(0,0,0)
            x(:,:,:,m1,m2,m3) = x(:,:,:,m1,m2,m3)*x(:,:,:,0,0,0)
        end if
    end function calc_inv_met_ind_0d
    !< \private array matrix version
    integer function calc_inv_met_arr(X,Y,deriv) result(ierr)
        character(*), parameter :: rout_name = 'calc_inv_met_arr'
        
        ! input / output
        real(dp), intent(inout) :: x(1:,1:,1:,1:,0:,0:,0:)                      !< X
        real(dp), intent(in) :: y(1:,1:,1:,1:,0:,0:,0:)                         !< Y
        integer, intent(in) :: deriv(:,:)                                       !< derivatives
        
        ! local variables
        integer :: id
        
        ! initialize ierr
        ierr = 0
        
        do id = 1, size(deriv,2)
            ierr = calc_inv_met_ind(x,y,deriv(:,id))
            chckerr('')
        end do
    end function calc_inv_met_arr
    !< \private array scalar version
    integer function calc_inv_met_arr_0d(X,Y,deriv) result(ierr)
        character(*), parameter :: rout_name = 'calc_inv_met_arr_0D'
        
        ! input / output
        real(dp), intent(inout) :: x(1:,1:,1:,0:,0:,0:)                         !< X
        real(dp), intent(in) :: y(1:,1:,1:,0:,0:,0:)                            !< Y
        integer, intent(in) :: deriv(:,:)                                       !< derivatives
        
        ! local variables
        integer :: id
        
        ! initialize ierr
        ierr = 0
        
        do id = 1, size(deriv,2)
            ierr = calc_inv_met_ind_0d(x,y,deriv(:,id))
            chckerr('')
        end do
    end function calc_inv_met_arr_0d
 
    !> Calculate  the metric coefficients in  a coordinate system B  ! using the
    !metric  coefficients in  a coordinate  system  A and  the !  transformation
    !matrix \f$\overline{\text{T}}_\text{B}^\text{A}\f$.
    !!
    !! This is  based on the  general treatment  using the formula  described in
    !! \cite Weyens3D :
    !!  \f[ \vec{D}_1^{m_1} \vec{D}_2^{m_2} \vec{D}_3^{m_3} g_\text{B}
    !!      = \sum_1 \sum_2 \sum_3 C_1 C_2 C_3 \vec{D}^{i_1,j_1,k_1}
    !!      \overline{\text{T}}_\text{B}^\text{A}
    !!      \vec{D}^{i_2,j_2,k_2} g_\text{A} \vec{D}^{i_3,j_3,k_3}
    !!      \left(\overline{\text{T}}_\text{B}^\text{A}\right)^T \ ,
    !!  \f]
    !! with
    !!  \f[ k_l = m_l-i_l-j_l \ , \f]
    !! and with  the \f$\sum_i\f$  double summations, with  the first  one going
    !! from \f$0\f$ to \f$m_i\f$ and the second one from \f$m_j\f$ to \f$0\f$.
    !! 
    !! The coefficients \f$C_l\f$ are calculated as
    !!  \f[
    !!      C_l = \left(\begin{array}{c}m_l\\i_l,j_l,k_l\end{array}\right)
    !!      \equiv \frac{m_l!}{i_l!j_l!(m_l-i_l-j_l)!} \ .
    !!  \f]
    !!
    !! \note It is assumed that the lower order derivatives have been calculated
    !! already. If not, the results will be incorrect. This is not checked.
    integer function calc_g(g_A,T_BA,g_B,deriv,max_deriv) result(ierr)
        use num_utilities, only: calc_mult, conv_mat
        
        character(*), parameter :: rout_name = 'calc_g'
        
        ! input / output
        real(dp), intent(in) :: g_a(:,:,:,:,0:,0:,0:)                           !< \f$g_\text{A}\f$
        real(dp), intent(in) :: t_ba(:,:,:,:,0:,0:,0:)                          !< \f$\overline{\text{T}}_\text{B}^\text{A}\f$
        real(dp), intent(inout) :: g_b(:,:,:,:,0:,0:,0:)                        !< \f$g_\text{B}\f$
        integer, intent(in) :: deriv(3)                                         !< derivatives
        integer, intent(in) :: max_deriv                                        !< maximum derivative
        
        ! local variables
        integer :: i1, j1, i2, j2, i3, j3, k1, k2, k3                           ! counters
        integer :: m1, m2, m3                                                   ! alias for deriv
        real(dp), allocatable :: c1(:), c2(:), c3(:)                            ! coeff. for (il)
        real(dp), allocatable :: dum1(:,:,:,:)                                  ! dummy variable representing metric matrix for some derivatives
        real(dp), allocatable :: dum2(:,:,:,:)                                  ! dummy variable representing metric matrix for some derivatives
        
        ! initialize ierr
        ierr = 0
        
#if ldebug
        ! check the derivatives requested
        ierr = check_deriv(deriv,max_deriv,'calc_g')
        chckerr('')
#endif
        
        ! set ml
        m1 = deriv(1)
        m2 = deriv(2)
        m3 = deriv(3)
        
        ! set up dummies
        allocate(dum1(size(g_a,1),size(g_a,2),size(g_a,3),9))
        allocate(dum2(size(g_a,1),size(g_a,2),size(g_a,3),9))
        
        ! initialize the requested derivative to zero
        g_b(:,:,:,:,m1,m2,m3) = 0.0_dp
        
        ! calculate the terms in the summation
        d1: do j1 = 0,m1                                                        ! derivatives in first coordinate
            call calc_c(m1,j1,c1)                                               ! calculate coeff. C1
            do i1 = m1-j1,0,-1
                d2: do j2 = 0,m2                                                ! derivatives in second coordinate
                    call calc_c(m2,j2,c2)                                       ! calculate coeff. C2
                    do i2 = m2-j2,0,-1
                        d3: do j3 = 0,m3                                        ! derivatives in third coordinate
                            call calc_c(m3,j3,c3)                               ! calculate coeff. C3
                            do i3 = m3-j3,0,-1
                                ! set up ki
                                k1 = m1 - j1 - i1
                                k2 = m2 - j2 - i2
                                k3 = m3 - j3 - i3
                                ! calculate term to add to the summation
                                ierr = calc_mult(g_a(:,:,:,:,j1,j2,j3),&
                                    &t_ba(:,:,:,:,k1,k2,k3),dum1,3,&
                                    &transp=[.false.,.true.])
                                chckerr('')
                                ierr = calc_mult(t_ba(:,:,:,:,i1,i2,i3),&
                                    &dum1,dum2,3)
                                chckerr('')
                                ! convert result to lower-diagonal storage
                                ierr = conv_mat(dum2,dum1(:,:,:,1:6),3)
                                chckerr('')
                                g_b(:,:,:,:,m1,m2,m3) = g_b(:,:,:,:,m1,m2,m3) &
                                    &+c1(i1)*c2(i2)*c3(i3)*dum1(:,:,:,1:6)
                            end do
                        end do d3
                    end do
                end do d2
            end do
        end do d1
        
        ! deallocate local variables
        deallocate(dum1,dum2)
    contains
        ! Calculate the coeff. C  at the all i values for  the current value for
        ! j, optionally making use of the coeff. C at previous value for j:
        ! - If it is the very first value (0,m), then it is just equal to 1
        ! - If  it is  the first value  in the current  i-summation, then  it is
        ! calculated from the first coeff. of the previous i-summation:
        !   C(j,m) = C(j-1,m) * (m-j+1)/j
        ! - If it is not the first  value in the current i-summation, then it is
        ! calculated from the previous value of the current i-summation
        !   C(j,i) = C(j,i+1) * (i+1)/(m-i-j)
        ! - If  it is  the last  value in  the current  i-summation, then  it is
        ! divided by 2 if (m-j) is even (see [ADD REF])
        !> \private
        subroutine calc_c(m,j,C)
            ! input / output
            integer, intent(in) :: m, j
            real(dp), intent(inout), allocatable :: c(:)
            
            ! local variables
            integer :: i                                                        ! counter
            real(dp) :: cm_prev
            
            ! if j>0, save the value Cm_prev 
            if (j.gt.0) cm_prev = c(m-j+1)
            
            ! allocate C
            if (allocated(c)) deallocate (c)
            allocate(c(0:m-j))
            
            ! first value i = m-j
            if (j.eq.0) then                                                    ! first coeff., for j = 0
                c(m) = 1.0_dp
            else                                                                ! first coeff., for j > 0 -> need value Cm_prev from previous array
                c(m-j) = (m-j+1.0_dp)/j * cm_prev
            end if
            
            ! all other values i = m-1 .. 0
            do i = m-j-1,0,-1
                c(i) = (i+1.0_dp)/(m-i-j) * c(i+1)
            end do
        end subroutine
    end function calc_g
    
    !> \private 3-D scalar version with one derivative
    integer recursive function transf_deriv_3_ind(x_a,t_ba,x_b,max_deriv,&
        &deriv_B,deriv_A_input) result(ierr)
        use num_utilities, only: add_arr_mult, c
        
        character(*), parameter :: rout_name = 'transf_deriv_3_ind'
        
        ! input / output
        real(dp), intent(in) :: x_a(1:,1:,1:,0:,0:,0:)                          !< variable and derivs. in coord. system A
        real(dp), intent(in) :: t_ba(1:,1:,1:,1:,0:,0:,0:)                      !< transf. mat. and derivs. between coord. systems A and B
        real(dp), intent(inout) :: x_b(1:,1:,1:)                                !< requested derivs. of variable in coord. system B
        integer, intent(in) :: max_deriv                                        !< maximum degrees of derivs.
        integer, intent(in) :: deriv_b(:)                                       !< derivs. in coord. system B
        integer, intent(in), optional :: deriv_a_input(:)                       !< derivs. in coord. system A (optional)
        
        ! local variables
        integer :: id, jd, kd, ld                                               ! counters
        integer :: deriv_id_b                                                   ! holds the deriv. currently calculated
        integer, allocatable :: deriv_a(:)                                      ! holds either deriv_A or 0
        real(dp), allocatable :: x_b_x(:,:,:,:,:,:)                             ! X_B and derivs. D^p-q_A with extra 1 exchanged deriv. D_A
        integer, allocatable :: deriv_a_x(:)                                    ! holds A derivs. for exchanged X_B_x
        integer, allocatable :: deriv_b_x(:)                                    ! holds B derivs. for exchanged X_B_x
        
        ! initialize ierr
        ierr = 0
        
        ! setup deriv_A
        allocate(deriv_a(size(deriv_b)))
        if (present(deriv_a_input)) then
            deriv_a = deriv_a_input
        else
            deriv_a = 0
        end if
        
#if ldebug
        ! check the derivatives requested
        ! (every B deriv. needs all the A derivs. -> sum(deriv_B) needed)
        ierr = check_deriv(deriv_a + sum(deriv_b),max_deriv,'transf_deriv')
        chckerr('')
#endif
        
        ! detect first deriv. in the B coord. system that needs to be exchanged,
        ! with derivs in the A coord. system going from B coord. 1 to B coord. 2
        ! to B coord. 3
        ! if equal to  zero on termination, this means that  no more exchange is
        ! necessary
        deriv_id_b = 0
        do id = 1,3
            if (deriv_b(id).gt.0) then
                deriv_id_b = id
                exit
            end if
        end do
        
        ! calculate the  derivative in coord.  deriv_id_B of coord. system  B of
        ! one order lower than requested here
        ! check if we have reached the zeroth derivative in coord. B
        if (deriv_id_b.eq.0) then                                               ! return the function with its required A derivs.
            x_b = x_a(:,:,:,deriv_a(1),deriv_a(2),deriv_a(3))
        else                                                                    ! apply the formula to obtain X_B using exchanged derivs.
            ! allocate  the  exchanged  X_B_x  and  the  derivs.  in  A  and  B,
            ! corresponding to D^m-1_B X
            allocate(x_b_x(size(x_b,1),size(x_b,2),size(x_b,3),&
                &0:deriv_a(1),0:deriv_a(2),0:deriv_a(3)))
            allocate(deriv_a_x(size(deriv_a)))
            allocate(deriv_b_x(size(deriv_b)))
            
            ! initialize X_B
            x_b = 0.0_dp
            
            ! loop over the 3 terms in the sum due to the transf. mat.
            do id = 1,3
                ! calculate X_B_x for orders 0 up to deriv_A
                do jd = 0,deriv_a(1)
                    do kd = 0,deriv_a(2)
                        do ld = 0,deriv_a(3)
                            ! calculate the  derivs. in  A for X_B_x:  for every
                            ! component in the sum due  to the transf. mat., the
                            ! deriv. is one order  higher than [jd,kd,ld] at the
                            ! corresponding coord. id
                            deriv_a_x = [jd,kd,ld]
                            deriv_a_x(id) = deriv_a_x(id) + 1
                            
                            ! calculate the derivs. in B  for X_B_x: this is one
                            ! order lower in the coordinate deriv_id_B
                            deriv_b_x = deriv_b
                            deriv_b_x(deriv_id_b) = deriv_b_x(deriv_id_b) - 1
                            
                            ! recursively call the subroutine again to calculate
                            ! X_B_x for the current derivatives
                            ierr = transf_deriv_3_ind(x_a,t_ba,&
                                &x_b_x(:,:,:,jd,kd,ld),max_deriv,deriv_b_x,&
                                &deriv_a_x)
                            chckerr('')
                        end do
                    end do
                end do
                
                ! use X_B_x at this term in summation due to the transf. mat. to
                ! update X_B using the formula
                ierr = add_arr_mult(&
                    &t_ba(:,:,:,c([deriv_id_b,id],.false.),0:,0:,0:),&
                    &x_b_x(:,:,:,0:,0:,0:),x_b,deriv_a)
                chckerr('')
            end do
        end if
    end function transf_deriv_3_ind
    !> \private 3-D scalar version with multiple derivatives
    integer function transf_deriv_3_arr(X_A,T_BA,X_B,max_deriv,derivs) &
        &result(ierr)
        character(*), parameter :: rout_name = 'transf_deriv_3_arr'
        
        ! input / output
        real(dp), intent(in) :: x_a(1:,1:,1:,0:,0:,0:)                          !< variable and derivs. in coord. system A
        real(dp), intent(in) :: t_ba(1:,1:,1:,1:,0:,0:,0:)                      !< transf. mat. and derivs. between coord. systems A and B
        real(dp), intent(inout) :: x_b(1:,1:,1:,0:,0:,0:)                       !< requested derivs. of variable in coord. system B
        integer, intent(in) :: max_deriv                                        !< maximum degrees of derivs.
        integer, intent(in) :: derivs(:,:)                                      !< series of derivs. (in coordinate system B)
        
        ! local variables
        integer :: id                                                           ! counter
        
        ! initialize ierr
        ierr = 0
        
        do id = 1, size(derivs,2)
            ierr = transf_deriv_3_ind(x_a,t_ba,&
                &x_b(:,:,:,derivs(1,id),derivs(2,id),derivs(3,id)),&
                &max_deriv,derivs(:,id))
            chckerr('')
        end do
    end function transf_deriv_3_arr
    !> \private 3-D matrix version  with multiple derivatives
    integer function transf_deriv_3_arr_2d(X_A,T_BA,X_B,max_deriv,derivs) &
        &result(ierr)
        use num_utilities, only: is_sym, c
        
        character(*), parameter :: rout_name = 'transf_deriv_3_arr_2D'
        
        ! input / output
        real(dp), intent(in) :: x_a(1:,1:,1:,1:,0:,0:,0:)                       !< variable and derivs. in coord. system A
        real(dp), intent(in) :: t_ba(1:,1:,1:,1:,0:,0:,0:)                      !< transf. mat. and derivs. between coord. systems A and B
        real(dp), intent(inout) :: x_b(1:,1:,1:,1:,0:,0:,0:)                    !< requested derivs. of variable in coord. system B
        integer, intent(in) :: max_deriv                                        !< maximum degrees of derivs.
        integer, intent(in) :: derivs(:,:)                                      !< series of derivs. (in coordinate system B)
        
        ! local variables
        integer :: id, kd                                                       ! counters
        character(len=max_str_ln) :: err_msg                                    ! error message
        logical :: sym                                                          ! sym
        
        ! initialize ierr
        ierr = 0
        
        ! test whether the dimensions of X_A and X_B agree
        if (size(x_a,4).ne.size(x_b,4)) then
            err_msg = 'X_A and X_B need to have the same sizes'
            ierr = 1
            chckerr(err_msg)
        end if
        
        ierr = is_sym(3,size(x_a,4),sym)
        chckerr('')
        
        do kd = 1,size(x_a,4)
            do id = 1, size(derivs,2)
                ierr = transf_deriv_3_ind(x_a(:,:,:,kd,:,:,:),t_ba,&
                    &x_b(:,:,:,kd,derivs(1,id),derivs(2,id),derivs(3,id)),&
                    &max_deriv,derivs(:,id))
                chckerr('')
            end do
        end do
    end function transf_deriv_3_arr_2d
    !> \private 1-D scalar version with one derivative
    integer recursive function transf_deriv_1_ind(x_a,t_ba,x_b,max_deriv,&
        &deriv_B,deriv_A_input) result(ierr)
        use num_utilities, only: add_arr_mult
        
        character(*), parameter :: rout_name = 'transf_deriv_1_ind'
        
        ! input / output
        real(dp), intent(in) :: x_a(1:,0:)                                      !< variable and derivs. in coord. system A
        real(dp), intent(inout) :: x_b(1:)                                      !< requested derivs. of variable in coord. system B
        real(dp), intent(in) :: t_ba(1:,0:)                                     !< transf. mat. and derivs. between coord. systems A and B
        integer, intent(in), optional :: deriv_a_input                          !< derivs. in coord. system A (optional)
        integer, intent(in) :: deriv_b                                          !< derivs. in coord. system B
        integer, intent(in) :: max_deriv                                        !< maximum degrees of derivs.
        
        ! local variables
        integer :: jd                                                           ! counters
        integer :: deriv_a                                                      ! holds either deriv_A or 0
        real(dp), allocatable :: x_b_x(:,:)                                     ! X_B and derivs. D^p-q_A with extra 1 exchanged deriv. D_A
        integer :: deriv_a_x                                                    ! holds A derivs. for exchanged X_B_x
        integer :: deriv_b_x                                                    ! holds B derivs. for exchanged X_B_x
        
        ! initialize ierr
        ierr = 0
        
        ! setup deriv_A
        if (present(deriv_a_input)) then
            deriv_a = deriv_a_input
        else
            deriv_a = 0
        end if
        
#if ldebug
        ! check the derivatives requested
        ! (every B deriv. needs all the A derivs. -> sum(deriv_B) needed)
        ierr = check_deriv([deriv_a+deriv_b,0,0],max_deriv,'transf_deriv')
        chckerr('')
#endif
        
        ! calculate the  derivative in coord.  deriv_id_B of coord. system  B of
        ! one order lower than requested here
        ! check if we have reached the zeroth derivative in coord. B
        if (deriv_b.eq.0) then                                                  ! return the function with its required A derivs.
            x_b = x_a(:,deriv_a)
        else                                                                    ! apply the formula to obtain X_B using exchanged derivs.
            ! allocate  the  exchanged  X_B_x  and  the  derivs.  in  A  and  B,
            ! corresponding to D^m-1_B X
            allocate(x_b_x(size(x_b,1),0:deriv_a))
            
            ! initialize X_B
            x_b = 0.0_dp
            
            ! calculate X_B_x for orders 0 up to deriv_A
            do jd = 0,deriv_a
                ! calculate the derivs.  in A for X_B_x: for  every component in
                ! the  sum due  to the  transf. mat.,  the deriv.  is one  order
                ! higher than jd
                deriv_a_x = jd+1
                
                ! calculate the derivs. in B for X_B_x: this is one order lower
                deriv_b_x = deriv_b-1
                
                ! recursively call  the subroutine again to  calculate X_B_x for
                ! the current derivatives
                ierr = transf_deriv_1_ind(x_a,t_ba,x_b_x(:,jd),max_deriv,&
                    &deriv_b_x,deriv_a_input=deriv_a_x)
                chckerr('')
            end do
            
            ! use X_B_x to calculate X_B using the formula
            ierr = add_arr_mult(t_ba(:,0:),x_b_x(:,0:),x_b,[deriv_a,0,0])
            chckerr('')
        end if
    end function transf_deriv_1_ind
    
    !> \private flux version
    integer function calc_f_derivs_1(grid_eq,eq) result(ierr)
        use num_vars, only: eq_style, max_deriv, use_pol_flux_f
        use num_utilities, only: derivs, c, fac
        use eq_vars, only: max_flux_e
        
        character(*), parameter :: rout_name = 'calc_F_derivs_1'
        
        ! input / output
        type(grid_type), intent(in) :: grid_eq                                  !< equilibrium grid
        type(eq_1_type), intent(inout) :: eq                                    !< flux equilibrium variables
        
        ! local variables
        integer :: id
        real(dp), allocatable :: t_fe_loc(:,:)                                  ! T_FE(2,1)
        
        ! initialize ierr
        ierr = 0
        
        ! Transform depending on equilibrium style being used:
        !   1:  VMEC
        !   2:  HELENA
        select case (eq_style)
            case (1)                                                            ! VMEC
                ! user output
                call writo('Transform flux equilibrium quantities to flux &
                    &coordinates')
                
                call lvl_ud(1)
                
                ! set up local T_FE
                allocate(t_fe_loc(grid_eq%loc_n_r,0:max_deriv+1))
                if (use_pol_flux_f) then
                    do id = 0,max_deriv+1
                        t_fe_loc(:,id) = 2*pi/max_flux_e*eq%q_saf_E(:,id)
                    end do
                else
                    t_fe_loc(:,0) = -2*pi/max_flux_e
                    t_fe_loc(:,1:max_deriv+1) = 0._dp
                end if
                
                ! Transform the  derivatives in E coordinates  to derivatives in
                ! the F coordinates.
                do id = 0,max_deriv+1
                    ! pres_FD
                    ierr = transf_deriv(eq%pres_E,t_fe_loc,eq%pres_FD(:,id),&
                        &max_deriv+1,id)
                    chckerr('')
                    
                    ! flux_p_FD
                    ierr = transf_deriv(eq%flux_p_E,t_fe_loc,&
                        &eq%flux_p_FD(:,id),max_deriv+1,id)
                    chckerr('')
                    
                    ! flux_t_FD
                    ierr = transf_deriv(eq%flux_t_E,t_fe_loc,&
                        &eq%flux_t_FD(:,id),max_deriv+1,id)
                    chckerr('')
                    eq%flux_t_FD(:,id) = -eq%flux_t_FD(:,id)                    ! multiply by plus minus one
                    
                    ! q_saf_FD
                    ierr = transf_deriv(eq%q_saf_E,t_fe_loc,eq%q_saf_FD(:,id),&
                        &max_deriv+1,id)
                    chckerr('')
                    eq%q_saf_FD(:,id) = -eq%q_saf_FD(:,id)                      ! multiply by plus minus one
                    
                    ! rot_t_FD
                    ierr = transf_deriv(eq%rot_t_E,t_fe_loc,eq%rot_t_FD(:,id),&
                        &max_deriv+1,id)
                    chckerr('')
                    eq%rot_t_FD(:,id) = -eq%rot_t_FD(:,id)                      ! multiply by plus minus one
                end do
                
                call lvl_ud(-1)
            case (2)                                                            ! HELENA
                ! E derivatives are equal to F derivatives
                eq%flux_p_FD = eq%flux_p_E
                eq%flux_t_FD = eq%flux_t_E
                eq%pres_FD = eq%pres_E
                eq%q_saf_FD = eq%q_saf_E
                eq%rot_t_FD = eq%rot_t_E
        end select
    end function calc_f_derivs_1
    !> \private metric version
    integer function calc_f_derivs_2(eq) result(ierr)
        use num_vars, only: eq_style
        use num_utilities, only: derivs, c
        
        character(*), parameter :: rout_name = 'calc_F_derivs_2'
        
        ! input / output
        type(eq_2_type), intent(inout) :: eq                                    !< metric equilibrium variables
        
        ! local variables
        integer :: id
        integer :: pmone                                                        ! plus or minus one
        
        ! initialize ierr
        ierr = 0
        
        ! Set up pmone, depending on equilibrium style being used
        !   1:  VMEC
        !   2:  HELENA
        select case (eq_style)
            case (1)                                                            ! VMEC
                pmone = -1                                                      ! conversion VMEC LH -> RH coord. system
            case (2)                                                            ! HELENA
                pmone = 1
        end select
        
        ! user output
        call writo('Transform metric equilibrium quantities to flux &
            &coordinates')
        
        call lvl_ud(1)
        
        ! Transform the  derivatives in E coordinates  to derivatives in        
        ! the F coordinates.
        do id = 0,max_deriv
            ! g_FD
            ierr = transf_deriv(eq%g_F,eq%T_FE,eq%g_FD,max_deriv,derivs(id))
            chckerr('')
            
            ! h_FD
            ierr = transf_deriv(eq%h_F,eq%T_FE,eq%h_FD,max_deriv,derivs(id))
            chckerr('')
            
            ! jac_FD
            ierr = transf_deriv(eq%jac_F,eq%T_FE,eq%jac_FD,max_deriv,&
                &derivs(id))
            chckerr('')
        end do
        
        call lvl_ud(-1)
    end function calc_f_derivs_2
    
    !> Calculate memory in MB necessary for variables in equilibrium job.
    !!
    !! The size of  these variables is equal to the  product of the non-parallel
    !! dimensions (e.g. \c n_geo x \c loc_n_r), times the number of variables.
    !! 
    !! The latter should be:
    !!  -  PB3D:  only take  into  account  (\f$2\cdot6+1 =  13\f$)  equilibrium
    !!  variables \c g_FD, \c h_FD and  \c jac_FD, as the perturbation variables
    !!  are  divided in  jobs occupying  the remaning  space. These  equilibrium
    !!  variables are tabulated on the  equilibrium grid. Note that they contain
    !!  derivatives in extra dimensions, so that their size should be multiplied
    !!  by (\c max_deriv+1)^3.
    !!  - POST: take  into account these 13 equilibrium variables,  as well as 4
    !!  variables \c U and \c DU, with double size due to being complex, and the
    !!  additional dimension equal to \c n_mod_X, but without derivatives.
    !!
    !! \return ierr
    integer function calc_memory_eq(arr_size,n_par,mem_size) result(ierr)
        use iso_c_binding
        
        character(*), parameter :: rout_name = 'calc_memory_eq'
        
        ! input / output
        integer, intent(in) :: arr_size                                         !< size of part of X array
        integer, intent(in) :: n_par                                            !< number of parallel points
        real(dp), intent(inout) :: mem_size                                     !< total size
        
        ! local variables
        integer(C_SIZE_T) :: dp_size                                            ! size of dp
        
        ! initialize ierr
        ierr = 0
        
        call lvl_ud(1)
        
        ! get size of real variable
        dp_size = sizeof(1._dp)
        
        ! set memory size for metric equilibrium variables
        mem_size = 13._dp*arr_size
        mem_size = mem_size*n_par*(1+max_deriv)**3
        mem_size = mem_size*dp_size
        
        ! convert B to MB
        mem_size = mem_size*1.e-6_dp
        
        ! apply 50% safety factor (empirical)
        mem_size = mem_size*1.5_dp
        
        ! test overflow
        if (mem_size.lt.0) then
            ierr = 1
            chckerr('Overflow occured')
        end if
        
        call lvl_ud(-1)
    end function calc_memory_eq
    
    !> If this equilibrium job should be done, also increment \c eq_job_nr.
    logical function do_eq()
        use num_vars, only: eq_jobs_lims, eq_job_nr, rich_restart_lvl, &
            &jump_to_sol
        use rich_vars, only: rich_lvl
        
        ! increment equilibrium job nr.
        eq_job_nr = eq_job_nr + 1
        
        if (rich_lvl.eq.rich_restart_lvl .and. jump_to_sol) then                ! jumping to solution for restart level
            if (eq_job_nr.eq.1) then
                do_eq = .true.
            else
                do_eq = .false.
            end if
        else
            if (eq_job_nr.le.size(eq_jobs_lims,2)) then                         ! normal behavior, not yet at last equilibrium job
                do_eq = .true.
            else                                                                ! normal behavior, at last equilibrium job
                do_eq = .false.
            end if
        end if
    end function do_eq
    
    !> Returns string  with possible extension  with equilibrium job as  well as
    !! parallel job, or nothing if only one level and one parallel job.
    elemental character(len=max_str_ln) function eq_info()
        use num_vars, only: eq_jobs_lims, eq_job_nr
        
        if (size(eq_jobs_lims,2).gt.1) then
            eq_info = ' for Equilibrium job '//trim(i2str(eq_job_nr))//&
                &' of '//trim(i2str(size(eq_jobs_lims,2)))
        else
            eq_info = ''
        end if
    end function eq_info
    
    !> Prints information for equilibrium parallel job.
    subroutine print_info_eq(n_par_X_rich)
        use num_vars, only: eq_job_nr, eq_jobs_lims
        
        ! input / output
        integer, intent(in) :: n_par_x_rich                                     !< number of parallel points in this Richardson level
        
        ! user output
        if (size(eq_jobs_lims,2).gt.1) then
            call writo('but parallel limits for this job only ['//&
                &trim(i2str(eq_jobs_lims(1,eq_job_nr)))//'..'//&
                &trim(i2str(eq_jobs_lims(2,eq_job_nr)))//&
                &'] of [1..'//trim(i2str(n_par_x_rich))//']')
        end if
    end subroutine print_info_eq
end module eq_utilities