get_verr Subroutine

  subroutine get_verr (errest, erridx, phi, bpar)
    use le_grids, only: integrate_species, eint_error, trap_error, lint_error, wdim
    use le_grids, only: ng2, nlambda, new_trap_int, grid_has_trapped_particles
    use theta_grid, only: ntgrid
    use kt_grids, only: ntheta0, naky, aky, akx
    use species, only: nspec, spec
    use dist_fn_arrays, only: gnew, aj0, vpa, g_adjust, g_work, to_g_gs2, from_g_gs2
    use run_parameters, only: has_phi, has_apar, beta
    use gs2_layouts, only: g_lo
    use collisions, only: adjust_vnmult
    use array_utils, only: zero_array
    implicit none
    !> Indices of maximum error.
    integer, dimension (5,3), intent (out) :: erridx
    !> Estimated error.
    real, dimension (5,2), intent (out) :: errest
    !> Electrostatic potential and parallel magnetic field
    complex, dimension (-ntgrid:,:,:), intent (in) :: phi, bpar
    integer :: ig, it, ik, iglo, isgn, ntrap
    complex, dimension (:,:,:), allocatable :: phi_app, apar_app
    complex, dimension (:,:,:,:), allocatable :: phi_e, phi_l, phi_t, apar_e, apar_l, apar_t
    real, dimension (:,:), allocatable :: kmax
    real, dimension (:), allocatable :: wgt
    real, dimension(2) :: errtmp
    integer, dimension(3) :: idxtmp
    logical, parameter :: trap_flag = .true.
    logical :: compute_trapped_error

    allocate(wgt(nspec))

    if (has_phi) then
       allocate(phi_app(-ntgrid:ntgrid,ntheta0,naky))
       allocate(phi_e(-ntgrid:ntgrid,ntheta0,naky,wdim))
       allocate(phi_l(-ntgrid:ntgrid,ntheta0,naky,ng2))
    end if

    if (has_apar) then
       allocate(apar_app(-ntgrid:ntgrid,ntheta0,naky))
       allocate(apar_e(-ntgrid:ntgrid,ntheta0,naky,wdim))
       allocate(apar_l(-ntgrid:ntgrid,ntheta0,naky,ng2))
    end if

    compute_trapped_error = grid_has_trapped_particles() .and. new_trap_int

    ! first call to g_adjust converts gyro-averaged dist. fn. (g)
    ! into nonadiabatic part of dist. fn. (h)
    call g_adjust (gnew, phi, bpar, direction = from_g_gs2)

    ! take gyro-average of h at fixed total position (not g.c. position)
    ! Note here phi_app and apar_app are effectively antot and antota that
    ! would be returned from a call to getan_from_dnf(gnew, antot, antota,...)
    if (has_phi) then
       do iglo = g_lo%llim_proc, g_lo%ulim_proc
          do isgn = 1, 2
             do ig = -ntgrid, ntgrid
                g_work(ig, isgn, iglo) = aj0(ig, iglo) * gnew(ig, isgn, iglo)
             end do
          end do
       end do

       wgt = spec%z*spec%dens
       call integrate_species (g_work, wgt, phi_app)

       ! integrates dist fn of each species over v-space
       ! after dropping an energy grid point and returns
       ! phi_e, which contains the integral approximations
       ! to phi for each point dropped
       call eint_error (g_work, wgt, phi_e)

       ! integrates dist fn of each species over v-space
       ! after dropping an untrapped lambda grid point and returns phi_l.
       ! phi_l contains ng2 approximations for the integral over lambda that
       ! come from dropping different pts from the gaussian quadrature grid
       call lint_error (g_work, wgt, phi_l)

       ! next loop gets error estimate for trapped particles, if there are any
       if (compute_trapped_error) then
          ntrap = nlambda - ng2
          allocate(phi_t(-ntgrid:ntgrid, ntheta0, naky, ntrap))
          call zero_array(phi_t)
          call trap_error (g_work, wgt, phi_t)
       end if

    end if

    if (has_apar) then
       do iglo = g_lo%llim_proc, g_lo%ulim_proc
          do isgn = 1, 2
             do ig = -ntgrid, ntgrid
                g_work(ig, isgn, iglo) = aj0(ig, iglo) * vpa(ig, isgn, iglo) * &
                     gnew(ig, isgn, iglo)
             end do
          end do
       end do

       wgt = 2.0 *beta * spec%z * spec%dens * sqrt(spec%temp / spec%mass)
       call integrate_species (g_work, wgt, apar_app)

       call eint_error (g_work, wgt, apar_e)
       call lint_error (g_work, wgt, apar_l)
       if (compute_trapped_error) then
          ntrap = nlambda - ng2
          allocate(apar_t(-ntgrid:ntgrid, ntheta0, naky, ntrap))
          call zero_array(apar_t)
          call trap_error (g_work, wgt, apar_t)
       end if
    end if
    deallocate (wgt)

    ! second call to g_adjust converts from h back to g
    call g_adjust (gnew, phi, bpar, direction = to_g_gs2)

    allocate (kmax(ntheta0, naky))
    do ik = 1, naky
       do it = 1, ntheta0
          kmax(it,ik) = max(akx(it),aky(ik))
       end do
    end do

    errest = 0.0
    erridx = 0

    if (has_phi) then

       call estimate_error (phi_app, phi_e, kmax, errtmp, idxtmp)
       errest(1,:) = errtmp
       erridx(1,:) = idxtmp

       call estimate_error (phi_app, phi_l, kmax, errtmp, idxtmp)
       errest(2,:) = errtmp
       erridx(2,:) = idxtmp

       if (compute_trapped_error) then
          call estimate_error (phi_app, phi_t, kmax, errtmp, idxtmp, trap_flag)
          errest(3,:) = errtmp
          erridx(3,:) = idxtmp
          deallocate (phi_t)
       end if
       deallocate(phi_app, phi_e, phi_l)
    end if

    ! No trapped errors for apar?
    if (has_apar) then

       call estimate_error (apar_app, apar_e, kmax, errtmp, idxtmp)
       errest(4,:) = errtmp
       erridx(4,:) = idxtmp

       call estimate_error (apar_app, apar_l, kmax, errtmp, idxtmp)
       errest(5,:) = errtmp
       erridx(5,:) = idxtmp

       deallocate(apar_app, apar_e, apar_l)
    end if

    call adjust_vnmult(errest, compute_trapped_error)
  end subroutine get_verr