eq_utilities.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------!
!------------------------------------------------------------------------------!
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
 
    logical :: debug_calc_inv_met_ind = .false.                                 
#endif
    
    ! interfaces
    
    interface calc_f_derivs
        module procedure calc_f_derivs_1
        module procedure calc_f_derivs_2
    end interface
    
    interface calc_inv_met
        module procedure calc_inv_met_ind
        module procedure calc_inv_met_arr
        module procedure calc_inv_met_ind_0d
        module procedure calc_inv_met_arr_0d
    end interface
    
    interface transf_deriv
        module procedure transf_deriv_3_ind
        module procedure transf_deriv_3_arr
        module procedure transf_deriv_3_arr_2d
        module procedure transf_deriv_1_ind
    end interface
    
contains
    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:)
        real(dp), intent(in) :: y(1:,1:,1:,1:,0:,0:,0:)
        integer, intent(in) :: deriv(:)
        
        ! 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
    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:)
        real(dp), intent(in) :: y(1:,1:,1:,0:,0:,0:)
        integer, intent(in) :: deriv(:)
        
        ! 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
    
    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:)
        real(dp), intent(in) :: y(1:,1:,1:,1:,0:,0:,0:)
        integer, intent(in) :: deriv(:,:)
        
        ! 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
    
    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:)
        real(dp), intent(in) :: y(1:,1:,1:,0:,0:,0:)
        integer, intent(in) :: deriv(:,:)
        
        ! 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
 
    !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:)
        real(dp), intent(in) :: t_ba(:,:,:,:,0:,0:,0:)
        real(dp), intent(inout) :: g_b(:,:,:,:,0:,0:,0:)
        integer, intent(in) :: deriv(3)
        integer, intent(in) :: max_deriv
        
        ! 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])
        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
    
    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:)
        real(dp), intent(in) :: t_ba(1:,1:,1:,1:,0:,0:,0:)
        real(dp), intent(inout) :: x_b(1:,1:,1:)
        integer, intent(in) :: max_deriv
        integer, intent(in) :: deriv_b(:)
        integer, intent(in), optional :: deriv_a_input(:)
        
        ! 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
    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:)
        real(dp), intent(in) :: t_ba(1:,1:,1:,1:,0:,0:,0:)
        real(dp), intent(inout) :: x_b(1:,1:,1:,0:,0:,0:)
        integer, intent(in) :: max_deriv
        integer, intent(in) :: derivs(:,:)
        
        ! 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
    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:)
        real(dp), intent(in) :: t_ba(1:,1:,1:,1:,0:,0:,0:)
        real(dp), intent(inout) :: x_b(1:,1:,1:,1:,0:,0:,0:)
        integer, intent(in) :: max_deriv
        integer, intent(in) :: derivs(:,:)
        
        ! 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
    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:)
        real(dp), intent(inout) :: x_b(1:)
        real(dp), intent(in) :: t_ba(1:,0:)
        integer, intent(in), optional :: deriv_a_input
        integer, intent(in) :: deriv_b
        integer, intent(in) :: max_deriv
        
        ! 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
    
    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
        type(eq_1_type), intent(inout) :: eq
        
        ! 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
    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
        
        ! 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
    
    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
        integer, intent(in) :: n_par
        real(dp), intent(inout) :: mem_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
    
    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
    
    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
    
    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
        
        ! 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