fm_init_next_field_points Subroutine

  subroutine fm_init_next_field_points(self,field,pts_remain,kwork_filter,ifl)
    use kt_grids, only: naky, ntheta0
    use theta_grid, only: ntgrid
    use dist_fn, only: timeadv
    use fields_arrays, only: phi, apar, bpar, phinew, aparnew, bparnew
    use dist_fn_arrays, only: fieldeq, fieldeqa, fieldeqp
    use mp, only: barrier
    use array_utils, only: zero_array
    implicit none
    class(fieldmat_type), intent(inout) :: self
    complex, dimension(-ntgrid:ntgrid,ntheta0,naky), intent(inout) :: field
    integer, intent(inout) :: pts_remain
    integer, intent(in) :: ifl
    logical, dimension(ntheta0, naky), intent(inout) :: kwork_filter
    integer :: ik, is, iex, ic, iex_init, ig, it
    integer :: it_cur, ig_cur
    integer :: pts_set_here
    logical :: found_ig, found_it

    !AJ This is a nasty barrier, but it is in here as on ARCHER we have a problem
    !AJ that we get deadlock on large core counts without it.  This has been investigated 
    !AJ and looks like a problem with the Cray MPI library, but it's hard to get to the bottom of.
    !AJ The problem does not occur on small core counts (we see it between 2000 and 3000 cores)
    call barrier()
    
    pts_set_here=0
    !Loop over all ik
    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ik, is, ic, iex_init, iex, ig, it) &
    !$OMP SHARED(self, kwork_filter, field, ntgrid) &
    !$OMP REDUCTION(+:pts_set_here) &
    !$OMP SCHEDULE(static)
    do ik = 1, self%naky
       !Technically should be able to skip ik values which are not local (though could be empty)
       !/would have to remember to decrease pts_remain

       !If we've finished this ky set the 
       !filter and move on to next ik
       if (self%kyb(ik)%initdone) then
          kwork_filter(:, ik) = .true.
          cycle
       end if

       !Now look at all supercells with this ik
       do is = 1, self%kyb(ik)%nsupercell
          !If we've finished with this supercell
          !then set the filter and move on to next
          if (self%kyb(ik)%supercells(is)%initdone) then
             !Set filter for all it on supercell
             do ic = 1, self%kyb(ik)%supercells(is)%ncell
                kwork_filter(self%kyb(ik)%supercells(is)%cells(ic)%it_ind, ik) = .true.
             end do
             cycle
          end if

          !//If we get to this point then we know that we have at least
          !one grid point on supercell for which we've not made a delta-fn

          !Now look for the first point that we've not initialised
          iex_init = 0
          do iex = 1, self%kyb(ik)%supercells(is)%nextend
             !If we've already done this point then cycle
             if (self%kyb(ik)%supercells(is)%initialised(iex)) cycle
             !Else this is the first point so store it and exit
             iex_init = iex
             exit
          end do

          !If this is the last point then mark supercell as done
          !but don't set filter yet, this will be automatically
          !applied on next loop
          if (iex_init == self%kyb(ik)%supercells(is)%nextend) then
             self%kyb(ik)%supercells(is)%initdone = .true.
          endif

          !Mark this point as initialised
          self%kyb(ik)%supercells(is)%initialised(iex_init) = .true.

          !Increment the number of points handled in this call
          pts_set_here = pts_set_here + 1

          !Get the ig and it indices corresponding to this point
          call self%kyb(ik)%supercells(is)%iex_to_ig(iex_init,ig)
          call self%kyb(ik)%supercells(is)%iex_to_it(iex_init,it)

          !Initialise field
          field(ig, it, ik) = 1.0

          !Set the duplicate boundary point if we're at the left
          !boundary of this call and we are not the leftmost cell in
          !the chain.
          if ((ig == -ntgrid) .and. &
               (it /= self%kyb(ik)%supercells(is)%it_leftmost)) then
             !Get the left connection
             it = self%kyb(ik)%supercells(is)%get_left_it(it)
             field(ntgrid, it, ik) = 1.0
          end if
       end do !Ends loop over supercells

       !Decide if we've done this kyb
       self%kyb(ik)%initdone = all(self%kyb(ik)%supercells%initdone)

    end do !Ends loop over ik
    !$OMP END PARALLEL DO

    ! Update the number of points remaining
    pts_remain = pts_remain - pts_set_here

    !Now do the time advance
    call timeadv(phi,apar,bpar,phinew,aparnew,bparnew,0)
    !Now get the field equations
    call self%getfieldeq_nogath(phinew,aparnew,bparnew,fieldeq,fieldeqa,fieldeqp)

    !Now we need to store the field equations in the appropriate places
    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ik, ig_cur, it_cur, found_ig, found_it, it, is) &
    !$OMP SHARED(self, field, ntgrid, ntheta0, ifl, fieldeq, fieldeqa, fieldeqp) &
    !$OMP SCHEDULE(static)
    do ik = 1, self%naky
       !If no data then skip ik -- because we don't skip setting
       !the field in the first loop of this routine but we do skip
       !zeroing it out here on processors which don't have any interest
       !in this ky, we need to zero the field after this loop to
       !ensure it has been fully reset to 0 on all procs. We might
       !consider skipping the setting of the field in the first loop
       !instead, although this could mean the field is different on
       !different processors (which should be fine) and might complicate
       !how we track the number of points we need to visit etc.
       if (self%kyb(ik)%is_empty) cycle

       !Loop until we've done all supercells
       ! Note we dont use gs2_max_abs here as we are already inside
       ! an OpenMP region.
       do while (maxval(abs(field(: , : ,ik))) > 0)
          !Reset variables
          ig_cur = -ntgrid-1
          it_cur = 0
          found_ig = .false.
          found_it = .false.
          
          !Look for first entry which was set to a delta-fn
          theta0: do it = 1, ntheta0
             ! Find the first theta grid point equal to 1.0
             ig_cur = findloc(field(:, it, ik), cmplx(1.0, 0.0), dim = 1)
             found_ig = ig_cur > 0
             ! If we didn't find a theta grid value move on to
             ! next it
             if (.not. found_ig) cycle
             ! Adjust ig_cur to correct index range (findloc
             ! assumes arrays are 1 indexed)
             ig_cur = ig_cur - (ntgrid + 1)
             if (found_ig) then
                it_cur = it
                found_it = .true.
                exit theta0
             end if
          end do theta0

          !Now check we've found a point, else cycle
          !NOTE: If this triggers we're probably in trouble
          !as nothing will be different on the next iteration
          !and hence we will be stuck in an infinite loop.
          !As such this should really trigger an abort.
          if (.not. found_it) cycle

          !Zero out the field at this point so we ignore
          !it on future iterations
          field(ig_cur, it_cur, ik)=0.0

          !Find the appropriate supercell
          is = self%kyb(ik)%is_from_it(it_cur)

          !Also zero out the boundary point if needed
          if((ig_cur == -ntgrid) .and. &
               (it_cur /= self%kyb(ik)%supercells(is)%it_leftmost)) then
             it = self%kyb(ik)%supercells(is)%get_left_it(it_cur)
             field(ntgrid, it, ik) = 0.0
          end if

          !Now store data
          call self%kyb(ik)%store_fq(fieldeq, fieldeqa, fieldeqp, ifl, it_cur, ig_cur)
       end do !Ends loop over points
    end do !Ends loop over ky
    !$OMP END PARALLEL DO

    ! See comment in above loop for why this is needed
    call zero_array(field)

  end subroutine fm_init_next_field_points