advance_nonlinear_term_ab3 Subroutine

  subroutine advance_nonlinear_term_ab3(g_state, phinew, aparnew, bparnew, dt, &
       fields_func, source_func)
    use theta_grid, only: ntgrid
    use gs2_time, only: get_adams_bashforth_coefficients
    use gs2_time, only: code_dt_min, dt_not_set, code_dt_prev1, code_dt_prev2
    use mp, only: mp_abort, max_allreduce
    use gs2_layouts, only: g_lo
    use array_utils, only: copy, zero_array
    implicit none
    complex, dimension (-ntgrid:, :, g_lo%llim_proc:), intent(in out) :: g_state
    complex, dimension (-ntgrid:,:,:), intent (in out) :: phinew, aparnew, bparnew
    !> Arrays to hold the nonlinear source history similar to gexp_1/2/3
    complex, dimension(:, :, :), allocatable :: gnl_1, gnl_2, gnl_3
    real, intent(in) :: dt
    real :: internal_time, time_step, max_vel, half_time
    real, dimension(:), allocatable :: ab_coeffs
    integer :: counter, iglo
    procedure(invert_field_func) :: fields_func
    procedure(source_term_func) :: source_func

    ! Initialise internal variables
    internal_time = 0.0
    counter = 1

    allocate(gnl_1(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc)) ; call zero_array(gnl_1)
    allocate(gnl_2(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc)) ; call zero_array(gnl_2)
    allocate(gnl_3(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc)) ; call zero_array(gnl_3)

    ! Reset time step sizes to reflect fact we're restarting the explicit scheme
    code_dt_prev1 = dt_not_set
    code_dt_prev2 = dt_not_set

    ! Keep taking explicit steps until we have advanced by dt
    do while (internal_time < dt)
       ! Get the current nonlinear source term and cfl limit
       call source_func(g_state, gnl_1, phinew, aparnew, bparnew, &
            max_vel, need_to_adjust = .not. advance_nonadiabatic_dfn, &
            calculate_cfl_limit = .true.)
       nrhs_calls = nrhs_calls + 1

       ! Get the maximum velocity across all processors
       ! and convert to time_step estimate.
       call max_allreduce(max_vel)
       time_step = 1.0 / max_vel
       time_step = time_step * time_step_safety_factor

       ! If the cfl time step is too small then abort
       if (time_step < code_dt_min) call mp_abort("Time step has become too small in advance_nonlinear_term.", .true.)

       ! Now make sure we're don't overshoot the target time by limiting time step.
       time_step = min(time_step, dt - internal_time)

       ! Get the Adams-Bashforth coefficients for the current collection of time step sizes
       ab_coeffs = get_adams_bashforth_coefficients(target_dt = time_step)

       half_time = 0.5 * time_step

       ! Advance solution by time_step
       select case(counter)
       case(1)
          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(iglo) &
          !$OMP SHARED(g_lo, g_state, half_time, ab_coeffs, gnl_1) &
          !$OMP SCHEDULE(static)
          do iglo = g_lo%llim_proc, g_lo%ulim_proc
             g_state(:, :, iglo) = g_state(:, :, iglo) + half_time *( &
                  ab_coeffs(1) * gnl_1(:, :, iglo)       &
                  )
          end do
          !$OMP END PARALLEL DO
       case(2)
          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(iglo) &
          !$OMP SHARED(g_lo, g_state, half_time, ab_coeffs, gnl_1, gnl_2) &
          !$OMP SCHEDULE(static)
          do iglo = g_lo%llim_proc, g_lo%ulim_proc
             g_state(:, :, iglo) = g_state(:, :, iglo) + half_time *( &
                  ab_coeffs(1) * gnl_1(:, :, iglo) +     &
                  ab_coeffs(2) * gnl_2(:, :, iglo)       &
                  )
          end do
          !$OMP END PARALLEL DO
       case default
          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(iglo) &
          !$OMP SHARED(g_lo, g_state, half_time, ab_coeffs, gnl_1, gnl_2, gnl_3) &
          !$OMP SCHEDULE(static)
          do iglo = g_lo%llim_proc, g_lo%ulim_proc
             g_state(:, :, iglo) = g_state(:, :, iglo) + half_time *( &
                  ab_coeffs(1) * gnl_1(:, :, iglo) +     &
                  ab_coeffs(2) * gnl_2(:, :, iglo) +     &
                  ab_coeffs(3) * gnl_3(:, :, iglo)       &
                  )
          end do
          !$OMP END PARALLEL DO
       end select

       ! Calculate the fields consistent with the new g
       call calculate_fields(fields_func, g_state, phinew, aparnew, bparnew)

       ! Increment internal state
       internal_time = internal_time + time_step
       counter = counter + 1
       nsuccessful_steps = nsuccessful_steps + 1

       ! Shift time history along one point
       code_dt_prev2 = code_dt_prev1
       code_dt_prev1 = time_step
       call copy(gnl_2, gnl_3)
       call copy(gnl_1, gnl_2)
    end do

  end subroutine advance_nonlinear_term_ab3