estimate_error Function

  real function estimate_error(g1, g2, g3) result(error)
    use gs2_time, only: get_adams_bashforth_coefficients, time_history_available
    use mp, only: max_allreduce
    use array_utils, only: gs2_max_abs
    implicit none
    complex, dimension(:, :, :), intent(in) :: g1, g2, g3
    complex, dimension(:, :, :), allocatable :: source_order_3, source_order_2
    real, dimension(2) :: intermediate_error
    real, dimension(:), allocatable :: coeff_3rd, coeff_2nd

    allocate(coeff_3rd(min(time_history_available(), 3)))
    allocate(coeff_2nd(min(time_history_available(), 2)))

    ! Check that we actually have 3rd and 2nd order coefficients.
    ! This could fail on the first two time steps of a simulation, for
    ! example, where we don't have enough time history to allow a
    ! 2nd/3rd order calculation. For now just return an error of
    ! error_target times cfl squared. This is chosen such the target
    ! timestep remains the current timestep.
    if ( (size(coeff_3rd) /= 3) .or. (size(coeff_2nd) /= 2) ) then
       error = error_target *cfl*cfl
       return
    end if

    ! Calculate the coefficients required for variable time step schemes
    coeff_3rd = get_adams_bashforth_coefficients(maximum_order = 3)
    coeff_2nd = get_adams_bashforth_coefficients(maximum_order = 2)

    ! Calculate the source terms using 3rd and 2nd order methods
    allocate(source_order_3, mold = g1)
    allocate(source_order_2, mold = g1)
    source_order_3 = coeff_3rd(1) * g1 + coeff_3rd(2) * g2 + coeff_3rd(3) * g3
    source_order_2 = coeff_2nd(1) * g1 + coeff_2nd(2) * g2

    ! Here we find the maximum absolute error
    intermediate_error(1) = gs2_max_abs(source_order_3 - source_order_2)

    ! Now we find the maximum of source_order_2 to scale the absolute
    ! error with.  Note we use this rather than source_order_3 as the
    ! higher order method has a smaller region of stability so is
    ! perhaps more likely to blow up.
    intermediate_error(2) = gs2_max_abs(source_order_2)

    ! Reduce over all processors
    call max_allreduce(intermediate_error)

    ! Now we form a crude relative error by dividing by the maximum
    ! source value. We do this rather than forming the true relative
    ! error (i.e. dividing the entire absolute difference by the source)
    ! to avoid issues with dividing by small numbers. This will tend to
    ! underestimate the real relative error.
    error = intermediate_error(1) / intermediate_error(2)
  end function estimate_error