do_write_avg_moments Subroutine

  subroutine do_write_avg_moments
    use diagnostics_moments, only: getmoms
    use mp, only: proc0
    use kt_grids, only: ntheta0, naky
    use species, only: nspec
    use fields_arrays, only: phinew, bparnew
    use gs2_io, only: nc_loop_moments
    use theta_grid, only: ntgrid
    implicit none
    integer :: it, is
    complex, dimension (-ntgrid:ntgrid,ntheta0,naky,nspec) :: ntot, density, &
         upar, tpar, tperp, qparflux, pperpj1, qpperpj1
    complex, dimension (ntheta0, nspec) :: ntot00, density00, upar00, tpar00, tperp00
    complex, dimension (ntheta0) :: phi00
    real, dimension (nspec) :: ntot2, ntot20, tpar2, tperp2
    real, dimension (ntheta0, naky, nspec) :: ntot2_by_mode, ntot20_by_mode
    real, dimension (ntheta0, naky, nspec) :: tpar2_by_mode, tperp2_by_mode

    ! We seem to call this routine a lot in one loop
    call getmoms (phinew, bparnew, ntot, density, upar, tpar, tperp, qparflux, pperpj1, qpperpj1)

    if (proc0) then
       do is = 1, nspec
          do it = 1, ntheta0
             call get_fldline_avg(ntot(-ntgrid:ntgrid,it,1,is),ntot00(it,is))
             call get_fldline_avg(density(-ntgrid:ntgrid,it,1,is),density00(it,is))
             call get_fldline_avg(upar(-ntgrid:ntgrid,it,1,is),upar00(it,is))
             call get_fldline_avg(tpar(-ntgrid:ntgrid,it,1,is),tpar00(it,is))
             call get_fldline_avg(tperp(-ntgrid:ntgrid,it,1,is),tperp00(it,is))
          end do
       end do

       do it = 1, ntheta0
          call get_fldline_avg(phinew(-ntgrid:ntgrid,it,1),phi00(it))
       end do

       do is = 1, nspec
          call get_vol_average (ntot(:,:,:,is), ntot(:,:,:,is), &
               ntot2(is), ntot2_by_mode(:,:,is))
          call get_vol_average (tpar(:,:,:,is), tpar(:,:,:,is), &
               tpar2(is), tpar2_by_mode(:,:,is))
          call get_vol_average (tperp(:,:,:,is), tperp(:,:,:,is), &
               tperp2(is), tperp2_by_mode(:,:,is))
       end do

       do is = 1, nspec
          call get_vol_average (ntot(0,:,:,is), ntot(0,:,:,is), &
               ntot20(is), ntot20_by_mode(:,:,is))
       end do

       call nc_loop_moments (get_netcdf_file_id(), nout, ntot2, ntot2_by_mode, ntot20, &
            ntot20_by_mode, phi00, ntot00, density00, upar00, &
            tpar00, tperp00, tpar2_by_mode, tperp2_by_mode)
    end if
  end subroutine do_write_avg_moments