set_grids Subroutine

  subroutine set_grids
    use species, only: nspec, spec, f0_maxwellian
    use theta_grid, only: ntgrid, nbset, bset, eps_trapped
    use file_utils, only: open_output_file, close_output_file
    use mp, only: proc0
    implicit none
    integer :: is
    logical :: has_maxwellian_species
    integer :: unit
    allocate (speed(negrid,nspec),speed_maxw(negrid),w_maxw(negrid),energy_maxw(negrid))

    w_maxw = 0.0
    energy_maxw = 0.0
    speed_maxw = 0.0
    has_maxwellian_species = .false. 

    speed = sqrt(energy)

    do is = 1,nspec
      if  (spec(is)%f0type == f0_maxwellian) then
         has_maxwellian_species = .true.
         w_maxw(:) = w(:,is)
         energy_maxw(:) = energy(:,is)
         speed_maxw(:) = speed(:,is)
         exit
      end if
    end do
    if (.not. has_maxwellian_species) write(*,*) &
        'Warning; no maxwellian species; collisions will fail'

    if (trapped_particles .and. eps_trapped > epsilon(0.0)) then
       nlambda = ng2+nbset
       lmax = nlambda-1
    else
       nlambda = ng2
       lmax = nlambda
    end if
    nxi = max(2*nlambda-1, 2*ng2)
    allocate (al(nlambda))
    allocate (wl(-ntgrid:ntgrid,nlambda))
    allocate (jend(-ntgrid:ntgrid))
    allocate (forbid(-ntgrid:ntgrid,nlambda))
    allocate (is_ttp(-ntgrid:ntgrid,nlambda))
    allocate (is_bounce_point(-ntgrid:ntgrid,nlambda))
    allocate (can_be_ttp(nlambda))
    if (grid_has_trapped_particles()) then
       al(ng2+1:nlambda) = 1.0/bset
    end if
    call lgridset

    if (proc0) then
       call open_output_file(unit, '.vspace_integration_error')
       call report_velocity_integration_error_estimate(unit)
       call close_output_file(unit)
    end if

  end subroutine set_grids