do_write_fyx Subroutine

  subroutine do_write_fyx (unit, phi, bpar)
    use mp, only: proc0, send, receive, barrier
    use gs2_layouts, only: il_idx, ig_idx, ik_idx, it_idx, is_idx, isign_idx, ie_idx
    use gs2_layouts, only: idx_local, proc_id, yxf_lo, accelx_lo, g_lo
    use le_grids, only: al, energy=>energy_maxw, forbid, negrid, nlambda
    use theta_grid, only: bmag, ntgrid
    use dist_fn_arrays, only: gnew, g_adjust, g_work, to_g_gs2, from_g_gs2
    use nonlinear_terms, only: accelerated
    use gs2_transforms, only: transform2, init_transforms
    use array_utils, only: zero_array, copy
#ifdef SHMEM
    use shm_mpi3, only: shm_alloc, shm_free
#endif
    integer, intent(in) :: unit
    !> Electrostatic potential and parallel magnetic field
    complex, dimension (-ntgrid:,:,:), intent (in) :: phi, bpar
    real, dimension (:,:), allocatable :: grs, gzf
#ifndef SHMEM
    real, dimension (:,:,:), allocatable :: agrs, agzf
#else
    real, dimension (:,:,:), pointer, contiguous :: agrs => null(), agzf => null()
    complex, dimension(:,:,:), pointer, contiguous :: g0_ptr => null()
#endif
    real, dimension (:), allocatable :: agp0, agp0zf
    real :: gp0, gp0zf
    integer :: ig, it, ik, il, ie, is, iyxlo, isign, iaclo, iglo, acc
    logical, save :: first = .true.

    if (accelerated) then
#ifndef SHMEM
       allocate (agrs(2*ntgrid+1, 2, accelx_lo%llim_proc:accelx_lo%ulim_alloc))
       allocate (agzf(2*ntgrid+1, 2, accelx_lo%llim_proc:accelx_lo%ulim_alloc))
#else
       call shm_alloc(agrs, [1, 2*ntgrid+1, 1, 2, accelx_lo%llim_proc, accelx_lo%ulim_alloc])
       call shm_alloc(agzf, [1, 2*ntgrid+1, 1, 2, accelx_lo%llim_proc, accelx_lo%ulim_alloc])
       call shm_alloc(g0_ptr, [1, 2*ntgrid+1, 1, 2, accelx_lo%llim_proc, accelx_lo%ulim_alloc])
#endif
       allocate (agp0(2), agp0zf(2))
       call zero_array(agrs); call zero_array(agzf); agp0 = 0.0; agp0zf = 0.0
    else
#ifdef SHMEM
       g0_ptr => g_work
#endif
       allocate (grs(yxf_lo%ny, yxf_lo%llim_proc:yxf_lo%ulim_alloc))
       allocate (gzf(yxf_lo%ny, yxf_lo%llim_proc:yxf_lo%ulim_alloc))
       call zero_array(grs); call zero_array(gzf); gp0 = 0.0; gp0zf = 0.0
    end if

    if (first) then
       if (proc0) then
          acc = 0
          if (accelerated) acc = 1
          write(unit,*) 2*negrid*nlambda, bmag(0), acc
       end if

       first = .false.
    end if

    call g_adjust (gnew, phi, bpar, direction = from_g_gs2)

#ifndef SHMEM
    call copy(gnew, g_work)
#else
    g0_ptr = gnew
#endif

#ifndef SHMEM
    if (accelerated) then
       call transform2 (g_work, agrs)
    else
       call transform2 (g_work, grs)
    end if

    call zero_array(g_work)
    do iglo=g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo,iglo)
       if (ik == 1) g_work(:,:,iglo) = gnew(:,:,iglo)
    end do

    call g_adjust (gnew, phi, bpar, direction = to_g_gs2)

    if (accelerated) then
       call transform2 (g_work, agzf)
    else
       call transform2 (g_work, gzf)
    end if

#else

    if (accelerated) then
       call transform2 (g0_ptr, agrs)
    else
       call transform2 (g0_ptr, grs)
    end if

    g0_ptr = 0.0
    do iglo=g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo,iglo)
       if (ik == 1) g0_ptr(:,:,iglo) = gnew(:,:,iglo)
    end do

    call g_adjust (gnew, phi, bpar, direction = to_g_gs2)

    if (accelerated) then
       call transform2 (g0_ptr, agzf)
    else
       call transform2 (g0_ptr, gzf)
    end if

#endif

    if (accelerated) then
       do iaclo=accelx_lo%llim_world, accelx_lo%ulim_world
          it = it_idx(accelx_lo, iaclo)
          ik = ik_idx(accelx_lo, iaclo)
          if (it == 1 .and. ik == 1) then
             il = il_idx(accelx_lo, iaclo)
             ig = 0
             if (.not. forbid(ig,il)) then
                ie = ie_idx(accelx_lo, iaclo)
                is = is_idx(accelx_lo, iaclo)

                if (proc0) then
                   if (.not. idx_local(accelx_lo, ik, it, il, ie, is)) then
                      call receive (agp0, proc_id(accelx_lo, iaclo))
                      call receive (agp0zf, proc_id(accelx_lo, iaclo))
                   else
                      agp0 = agrs(ig+ntgrid+1, :, iaclo)
                      agp0zf = agzf(ig+ntgrid+1, :, iaclo)
                   end if
                else if (idx_local(accelx_lo, ik, it, il, ie, is)) then
                   call send (agrs(ig+ntgrid+1, :, iaclo), 0)
                   call send (agzf(ig+ntgrid+1, :, iaclo), 0)
                end if

                if (proc0) then
                   write (unit, "(6(1x,e13.6))") energy(ie), al(il), &
                        agp0(1), agp0(2), agp0zf(1), agp0zf(2)
                end if
             end if
          end if
          call barrier
       end do
       deallocate(agp0, agp0zf)
#ifndef SHMEM
       deallocate(agrs, agzf)
#else
       call shm_free(agrs)
       call shm_free(agzf)
       call shm_free(g0_ptr)
#endif
    else
       do iyxlo=yxf_lo%llim_world, yxf_lo%ulim_world

          ig = ig_idx(yxf_lo, iyxlo)
          it = it_idx(yxf_lo, iyxlo)
          if (ig == 0 .and. it == 1) then
             il = il_idx(yxf_lo, iyxlo)
             if (.not. forbid(ig,il)) then
                ik = 1
                ie = ie_idx(yxf_lo, iyxlo)
                is = is_idx(yxf_lo, iyxlo)
                isign = isign_idx(yxf_lo, iyxlo)

                if (proc0) then
                   if (.not. idx_local(yxf_lo, ig, isign, it, il, ie, is)) then
                      call receive (gp0, proc_id(yxf_lo, iyxlo))
                      call receive (gp0zf, proc_id(yxf_lo, iyxlo))
                   else
                      gp0 = grs(ik, iyxlo)
                      gp0zf = gzf(ik, iyxlo)
                   end if
                else if (idx_local(yxf_lo, ig, isign, it, il, ie, is)) then
                   call send (grs(ik, iyxlo), 0)
                   call send (gzf(ik, iyxlo), 0)
                end if

                if (proc0) then
                   write (unit, "(4(1x,e13.6),i8)") energy(ie), al(il), &
                        gp0, gp0zf, isign
                end if
             end if
          end if
          call barrier
       end do
       deallocate (grs, gzf)
    end if

    if (proc0) flush (unit)

    if (proc0) then
       write(unit,*)
    end if

  end subroutine do_write_fyx