estimate_error Subroutine

  subroutine estimate_error (app1, app2, kmax_local, errest, erridx, trap)
    use kt_grids, only: naky, ntheta0
    use theta_grid, only: ntgrid
    use le_grids, only: ng2, jend, il_is_trapped
    implicit none
    complex, dimension (-ntgrid:,:,:), intent (in) :: app1
    complex, dimension (-ntgrid:,:,:,:), intent (in) :: app2
    real, dimension (:, :), intent (in) :: kmax_local
    logical, intent (in), optional :: trap
    real, dimension (:), intent (out) :: errest
    integer, dimension (:), intent (out) :: erridx
    integer :: ik, it, ig, ipt, end_idx
    integer :: igmax, ikmax, itmax, gpcnt
    real :: gdsum, gdmax, gpavg, gnsum, gsmax, gpsum, gptmp

    igmax = 0; ikmax = 0; itmax = 0
    gdsum = 0.0; gdmax = 0.0; gpavg = 0.0; gnsum = 0.0; gsmax = 1.0
    do ik = 1, naky
       do it = 1, ntheta0
          do ig = -ntgrid, ntgrid
             if (.not. present(trap) .or. il_is_trapped(jend(ig))) then
                gpcnt = 0; gpsum = 0.0

                if (present(trap)) then
                   end_idx = jend(ig)-ng2
                else
                   end_idx = size(app2, 4)
                end if

                do ipt=1,end_idx
                   gptmp = kmax_local(it, ik) * abs(app1(ig, it, ik) - app2(ig, it, ik, ipt))
                   gpsum = gpsum + gptmp
                   gpcnt = gpcnt + 1
                end do

                gpavg = gpsum/gpcnt

                if (gpavg > gdmax) then
                   igmax = ig
                   ikmax = ik
                   itmax = it
                   gdmax = gpavg
                   gsmax = kmax_local(it, ik) * abs(app1(ig, it, ik))
                end if

                gnsum = gnsum + gpavg
                gdsum = gdsum + kmax_local(it, ik) * abs(app1(ig, it, ik))
             end if
          end do
       end do
    end do

    gdmax = gdmax / gsmax

    erridx(1) = igmax
    erridx(2) = ikmax
    erridx(3) = itmax
    errest(1) = gdmax

    ! Guard against divide by zero. This can commonly occur on the first time step
    ! in runs with fapar /= 0 due to our initial conditions typically ensuring that
    ! the initial apar is identically zero.
    if (abs(gdsum) < epsilon(0.0)) then
       errest(2) = 0.0
    else
       errest(2) = gnsum / gdsum
    end if
  end subroutine estimate_error