advance_nonlinear_term Subroutine

  subroutine advance_nonlinear_term(g_state, istep, phi, apar, bpar, &
       dt, fields_func, source_func)
    use dist_fn_arrays, only: g_adjust, to_g_gs2, from_g_gs2
    use mp, only: mp_abort, get_mp_times, proc0
    use job_manage, only: time_message
    use nonlinear_terms, only: time_add_explicit_terms, time_add_explicit_terms_mpi
    use run_parameters, only: has_phi, has_apar, has_bpar
    use array_utils, only: copy
    implicit none
    integer, intent(in) :: istep
    complex, dimension (:,:,:), intent(in out) :: g_state
    complex, dimension (:,:,:), intent(in) :: phi, apar, bpar
    real, intent(in) :: dt
    complex, dimension (:,:,:), allocatable :: phinew, aparnew, bparnew
    real :: mp_total_after, mp_total
    procedure(invert_field_func) :: fields_func
    procedure(source_term_func) :: source_func

    call time_message(.false., time_add_explicit_terms, 'Explicit terms')
    call get_mp_times(total_time = mp_total)
    call reset_step_statistics

    ! Copy current fields to local copies
    allocate(phinew, mold = phi)
    allocate(aparnew, mold = apar)
    allocate(bparnew, mold = bpar)
    if (has_phi) call copy(phi, phinew)
    if (has_apar) call copy(apar, aparnew)
    if (has_bpar) call copy(bpar, bparnew)

    ! If we want to advance h rather than g then adjust
    ! here to get h. Note we rely on the call site also
    ! providing the correct fields_func to calculate fields
    ! from h rather than g. A different option would be
    ! to require the call site to always pass a way to calculate
    ! from g and a way from h and we then decide here which to use.
    if (advance_nonadiabatic_dfn) then
       call g_adjust(g_state, phinew, bparnew, direction = from_g_gs2)
    end if

    ! Call requested method for advancing nonlinear term
    select case(split_method_switch)
    case(split_method_ab3)
       call advance_nonlinear_term_ab3(g_state, phinew, aparnew, bparnew, dt, &
            fields_func, source_func)
    case(split_method_backwards_euler)
       call advance_nonlinear_term_beuler(g_state, phinew, aparnew, bparnew, dt, &
            fields_func, source_func)
    case(split_method_rk)
       call advance_nonlinear_term_rk(g_state, phinew, aparnew, bparnew, dt, &
            fields_func, source_func)
    case(split_method_picard)
       call advance_nonlinear_term_picard(g_state, phinew, aparnew, bparnew, dt, &
            fields_func, source_func)
    case default
       call mp_abort("Invalid split_method_switch", .true.)
    end select

    if (proc0 .and. show_statistics) call report_step_statistics(istep)

    ! If we advanced h then we need to reconstruct g from this
    if (advance_nonadiabatic_dfn) then
       call g_adjust(g_state, phi, bpar, direction = to_g_gs2)
    end if

    ! Now that we've advanced the explicit terms by dt and set the
    ! current g_state to this new result. Note that we do not update the
    ! fields, instead resetting them to their original values. The
    ! update to the fields is calculated as a part of the implicit
    ! field solve and including the NL contribution here leads to
    ! "double counting".

    call time_message(.false., time_add_explicit_terms, 'Explicit terms')
    call get_mp_times(total_time = mp_total_after)
    time_add_explicit_terms_mpi = time_add_explicit_terms_mpi + (mp_total_after - mp_total)
  end subroutine advance_nonlinear_term