write_f Subroutine

  subroutine write_f (last)
    use mp, only: proc0, send, receive
    use file_utils, only: open_output_file, close_output_file
    use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx, il_idx, ie_idx
    use gs2_layouts, only: idx_local, proc_id
    use le_grids, only: al, energy, forbid, negrid, nlambda
    use egrid, only: zeroes, x0
    use theta_grid, only: bmag
    use dist_fn_arrays, only: gnew
    use species, only: nspec

    integer :: iglo, ik, it, is
    integer :: ie, il, ig
    integer, save :: unit
    real :: vpa, vpe
    complex, dimension(2) :: gtmp
    !> If true, also close the file
    logical, intent(in)  :: last 

    real, dimension (:,:), allocatable, save :: xpts
    real, dimension (:), allocatable, save :: ypts

    if (.not. allocated(xpts)) then
       allocate(xpts(negrid,nspec))
       allocate(ypts(nlambda))

    ! should really just get rid of xpts and ypts
       xpts(1:negrid-1,:) = zeroes(1:negrid-1,:)
       xpts(negrid,:) = x0(:)
       ypts = 0.0
       ! change argument of bmag depending on which theta you want to write out
       do il=1,nlambda
          if (1.0-al(il)*bmag(0) .gt. 0.0) ypts(il) = sqrt(1.0-al(il)*bmag(0))
       end do

       if (proc0) then
          call open_output_file (unit, ".dist")
          write(unit, *) negrid*nlambda
       end if
    endif

    do iglo = g_lo%llim_world, g_lo%ulim_world
       ik = ik_idx(g_lo, iglo) ; if (ik /= 1) cycle
       it = it_idx(g_lo, iglo) ; if (it /= 1) cycle
       is = is_idx(g_lo, iglo) ; if (is /= 1) cycle
       ie = ie_idx(g_lo, iglo) 
       ig = 0
       il = il_idx(g_lo, iglo)
       if (idx_local (g_lo, ik, it, il, ie, is)) then
          if (proc0) then 
             gtmp = gnew(ig,:,iglo)
          else
             call send (gnew(ig,:,iglo), 0)
          endif
       else if (proc0) then
          call receive (gtmp, proc_id(g_lo, iglo))
       endif
       if (proc0) then
          if (.not. forbid(ig,il)) then
             vpa = sqrt(energy(ie,is)*max(0.0, 1.0-al(il)*bmag(ig)))
             vpe = sqrt(energy(ie,is)*al(il)*bmag(ig))
             write (unit, "(8(1x,e12.6))") vpa, vpe, energy(ie,is), al(il), &
                  xpts(ie,is), ypts(il), real(gtmp(1)), real(gtmp(2))
          end if
       end if
    end do
    if (proc0) write (unit, *)
    if (last .and. proc0) call close_output_file (unit)
    
  end subroutine write_f