init_connected_bc Subroutine

  subroutine init_connected_bc
    use kt_grids, only: naky, ntheta0
    use le_grids, only: mixed_wfb
    use gs2_layouts, only: g_lo
    use array_utils, only: zero_array
    use mp, only: proc0
    use file_utils, only: error_unit
    implicit none
    integer :: jshift0, n_links_max

    if (connectinit) return
    connectinit = .true.

    ! Skip setup of linked bc data in cases without linked boundaries.
    ! We might want to remove this check so that we can calculate
    ! the linked bc data for testing etc.
    if (boundary_option_switch /= boundary_option_linked) return

    ! Get the kx index spacing between connections
    jshift0 = compute_jshift0()

    allocate (itleft(ntheta0, naky), itright(ntheta0, naky))
    call compute_itleft_and_itright(jshift0, itleft, itright)

    allocate (connections(g_lo%llim_proc:g_lo%ulim_alloc))
    call compute_connections(itleft, itright, connections)

    allocate (save_h(2, g_lo%llim_proc:g_lo%ulim_alloc))
    call set_save_h(connections, save_h)

    call count_links(itleft, itright, l_links, r_links, n_links)
    n_links_max = maxval(n_links)

    ! wfb -- mixed treatment requires some extra links to be communicated
    if (n_links_max > 0 .and. mixed_wfb) n_links_max = n_links_max + 3

    ! Could we move this call to init_bc? Would have to check that it returns
    ! the expected values when l_links = r_links = 0
    call setup_fields_implicit_class_arrays(l_links, r_links, M_class, N_class, i_class)

    ! now set up communication pattern:
    ! excluding wfb
    call setup_connected_bc_redistribute(l_links, r_links, n_links_max, &
         links_p, links_h, no_connections)

    ! If no connections then don't need to do any more work here. One
    ! might expect us to be able to set no_connections = n_links_max == 0
    ! rather than having to call setup_connected_bc_redistribute to
    ! determine this. This would allow a little bit of work to be skipped.
    ! Note that no_connections is a global quantity (i.e. agreed value
    ! on all processors). We may wish to introduce a local version as well.
    ! At scale each processor may only be responsible for one pitch angle,
    ! for example. If this is a trapped pitch angle then it will not have
    ! any work to do associated with the boundaries yet it will still go
    ! through the allocations etc. associated with setting up linked boundaries
    if (no_connections) then
       if (proc0) write(error_unit(), &
            '("No connections found when setting up linked boundaries.")')
       return
    end if

    ! take care of wfb
    if (mixed_wfb) then
       call setup_connected_bc_redistribute_mixed_wfb(l_links, r_links, n_links_max, &
            wfb_p, wfb_h)
    end if

    ! Setup redistributes for passing the incoming boundary to the previous linked
    ! cell. Used to generalise the linked boundary conditions to allow gnew /= 0 at
    ! the linked location.
    if (start_from_previous_solution) then
       call setup_pass_incoming_boundary_to_connections(l_links, r_links, incoming_links)
    end if

    ! n_links_max is typically 2 * number of cells in largest supercell
    ! Note this is potentially quite wasteful in terms of storage. In flux
    ! tube setup there is likely only one supercell with the maximal number
    ! of links and many more with fewer (or zero) links yet we allocate the
    ! full length for _all_ iglo points. For the longest supercell we expect
    ! n_links ~ 2 * ntheta0/jtwist - 1 and hence n_links_max is ~ 2*( 1 + ntheta/jtwist)
    ! For simulations with high nx and/or low jtwist n_links_max could be >> ntheta
    ! and therefore g_adj may dominate memory consumption in some cases.
    allocate (g_adj (n_links_max, 2, g_lo%llim_proc:g_lo%ulim_alloc))
    call zero_array(g_adj)

  end subroutine init_connected_bc