kmesh_shell_reconstruct Subroutine

  subroutine kmesh_shell_reconstruct(kmesh_input, num_kpts, multi, dnn, nnshell, bweight)
    !================================================
    !
    !!  Include more shells to calculate higher-order finite difference: 2b, 3b, ... Nb shells
    !!  Note: some shells are overwritten
    !================================================
    use w90_types, only: kmesh_input_type, print_output_type

    implicit none

    ! arguments
    type(kmesh_input_type), intent(inout) :: kmesh_input
    integer, intent(in) :: num_kpts
    integer, intent(inout) :: multi(max(kmesh_input%search_shells, 6*kmesh_input%higher_order_n))   ! the number of bvectors in the shell
    real(kind=dp), intent(inout) :: dnn(max(kmesh_input%search_shells, 6*kmesh_input%higher_order_n))
    integer, intent(inout) :: nnshell(num_kpts, max(kmesh_input%search_shells, 6*kmesh_input%higher_order_n))
    real(kind=dp), intent(inout) :: bweight(kmesh_input%max_shells_h)

    ! local variables
    real(kind=dp) :: bweight_temp, fact
    integer :: shell, order, loop_j, temp_multi(kmesh_input%num_shells), temp_nnshell(num_kpts, kmesh_input%num_shells)
    real(kind=dp) :: temp_dnn(kmesh_input%num_shells)

    ! update new shells (after simplify the first-order shell list, e.g. 1,4,6, ... -> 1,2,3,...)
    do shell = 1, kmesh_input%num_shells
      temp_multi(shell) = multi(kmesh_input%shell_list(shell))
      temp_nnshell(:, shell) = nnshell(:, kmesh_input%shell_list(shell))
      temp_dnn(shell) = dnn(kmesh_input%shell_list(shell))
    end do

    do shell = 1, kmesh_input%num_shells
      kmesh_input%shell_list(shell) = shell
      multi(shell) = temp_multi(shell)
      nnshell(:, shell) = temp_nnshell(:, shell)
      dnn(shell) = temp_dnn(shell)
    end do

    do order = 2, kmesh_input%higher_order_n
      do shell = 1, kmesh_input%num_shells
        kmesh_input%shell_list((order - 1)*kmesh_input%num_shells + shell) = &
          (order - 1)*kmesh_input%num_shells + shell
        multi((order - 1)*kmesh_input%num_shells + shell) = multi(shell)
        nnshell(:, (order - 1)*kmesh_input%num_shells + shell) = nnshell(:, shell)
        dnn((order - 1)*kmesh_input%num_shells + shell) = dnn(shell)*order
      end do
    end do

    ! calculate new bweights w_b, w_2b, ..., w_Nb
    do shell = 1, kmesh_input%num_shells
      bweight_temp = bweight(shell)
      do order = 1, kmesh_input%higher_order_n
        fact = 1.0_dp/REAL(order**2, DP)
        do loop_j = 1, kmesh_input%higher_order_n
          if (loop_j == order) cycle
          fact = (fact*REAL(loop_j**2, DP))/REAL(loop_j**2 - order**2, DP)
        end do
        bweight(kmesh_input%num_shells*(order - 1) + shell) = bweight_temp*fact
      end do
    end do

    kmesh_input%num_shells = kmesh_input%num_shells*kmesh_input%higher_order_n

    return

  end subroutine kmesh_shell_reconstruct