advance_nonlinear_term_picard Subroutine

  subroutine advance_nonlinear_term_picard(g_state, phinew, &
       aparnew, bparnew, dt, fields_func, source_func)
    use theta_grid, only: ntgrid
    use gs2_time, only: code_dt_min
    use mp, only: mp_abort, proc0, max_allreduce
    use gs2_layouts, only: g_lo
    use array_utils, only: copy, gs2_max_abs, 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
    complex, dimension(:, :, :), allocatable :: gnl_1
    real, intent(in) :: dt
    real :: internal_time, time_step, max_vel, cfl_limit
    integer :: counter, ipic, iglo
    logical, parameter :: debug = .false.
    real, dimension(4) :: errors
    real, save :: time_step_last = -1
    procedure(invert_field_func) :: fields_func
    procedure(source_term_func) :: source_func
    complex, dimension(:, :, :), allocatable :: g_start
    integer, parameter :: npic_limit = 15, npic_low_limit = 3
    real, parameter :: time_step_factor = 2.0
    logical, parameter :: check_cfl = .true.
    logical :: converged
    real :: half_time

    ! Initialise internal variables
    if (time_step_last < 0) time_step_last = dt
    internal_time = 0.0
    counter = 1
    time_step = time_step_last

    allocate(g_start(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc))
    allocate(gnl_1(-ntgrid:ntgrid, 2, g_lo%llim_proc:g_lo%ulim_alloc))
    call zero_array(g_start) ; call zero_array(gnl_1)
    ! Keep taking explicit steps until we have advanced by dt
    do while (internal_time < dt)
       call copy(g_state, g_start)
       half_time = 0.5 * time_step

       if(proc0.and.debug) print*,'Taking step with size',time_step,'at time',internal_time
       ipic = 1

       call calculate_fields(fields_func, g_state, phinew, aparnew, bparnew)

       ! Get the current nonlinear source term
       call source_func(g_state, gnl_1, &
            phinew, aparnew, bparnew, &
            max_vel, need_to_adjust = .not. advance_nonadiabatic_dfn, &
            calculate_cfl_limit = check_cfl)
       nrhs_calls = nrhs_calls + 1

       errors(1) = 0.0
       errors(2) = gs2_max_abs(gnl_1) * half_time
       errors(3) = gs2_max_abs(g_start)
       errors(4) = max_vel
       call max_allreduce(errors)

       ! Check for CFL limit - we can simply rescale the change in g and carry on
       ! so here we simply rescale the effective time step,
       cfl_limit = 1.0 / errors(4)
       if (time_step > cfl_limit .and. check_cfl) then
          ! Drop the effective time step a little from the limit
          cfl_limit = 0.9 * cfl_limit
          if(proc0.and.debug) print*,'Dropping time step due to cfl from',time_step,'to',cfl_limit
          errors(2) = errors(2) * cfl_limit / time_step
          time_step = cfl_limit
          half_time = time_step * 0.5
       end if

       converged = picard_converged(errors)

       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(iglo) &
       !$OMP SHARED(g_state, g_start, gnl_1, half_time, g_lo)&
       !$OMP SCHEDULE(static)
       do iglo = g_lo%llim_proc, g_lo%ulim_proc
          g_state(:, :, iglo) = g_start(:, :, iglo) + gnl_1(:, :, iglo) * half_time
       end do
       !$OMP END PARALLEL DO

       do while(.not. converged)
          ! Note the first step (i.e. call before this loop) is included in
          ! ipic here (i.e. it counts towards our approach towards npic_limit).
          ipic = ipic + 1
          if (ipic > npic_limit) exit

          call calculate_fields(fields_func, g_state, phinew, aparnew, bparnew)

          ! Get the current nonlinear source term
          call source_func(g_state, gnl_1, &
               phinew, aparnew, bparnew, &
               max_vel, need_to_adjust = .not. advance_nonadiabatic_dfn, &
               calculate_cfl_limit = check_cfl)
          nrhs_calls = nrhs_calls + 1

          errors(1) = errors(2)
          errors(2) = gs2_max_abs(gnl_1) * half_time
          errors(4) = max_vel
          call max_allreduce(errors)
          cfl_limit = 1.0 / errors(4)

          converged = picard_converged(errors)

          if(proc0.and.debug) print*,'ipic',ipic,converged,'|',errors

          if (time_step > cfl_limit .and. .not. converged) exit

          !$OMP PARALLEL DO DEFAULT(none) &
          !$OMP PRIVATE(iglo) &
          !$OMP SHARED(g_state, g_start, gnl_1, half_time, g_lo)&
          !$OMP SCHEDULE(static)
          do iglo = g_lo%llim_proc, g_lo%ulim_proc
             g_state(:, :, iglo) = g_start(:, :, iglo) + gnl_1(:, :, iglo) * half_time
          end do
          !$OMP END PARALLEL DO
       end do

       if (.not. converged) then
          ! If we didn't converge in the limit then reset state
          ! to the start of the step and reduce the time step.
          call copy(g_start, g_state)
          time_step = time_step / time_step_factor
          nfailed_steps = nfailed_steps + 1

          ! If the 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.)
       else
          internal_time = internal_time + time_step
          nsuccessful_steps = nsuccessful_steps + 1
          ! If we converged quickly then increase the time step
          if (ipic <= npic_low_limit) time_step = time_step * time_step_factor
       end if

       !Now update the time step, capped to dt
       time_step = min(time_step, dt)

       ! Make sure we don't overshoot the target time by limiting time step.
       if(internal_time < dt) then
          time_step = min(time_step, dt - internal_time)
          counter = counter + 1
       end if

    end do

    if(proc0.and.debug) print*,'Completed Picard step in',counter,'steps to',internal_time,'last time step',time_step

    time_step_last = time_step
  contains
    pure logical function picard_converged(errors) result(converged)
      implicit none
      real, dimension(4), intent(in) :: errors
      converged = &
           !Change in solution less than abs tol
           (errors(2) < absolute_tolerance) .or.  &
           ! Change in the change, less than abs tol
           (abs(errors(2) - errors(1)) < absolute_tolerance) .or. &
           ! Change is less than the relative tolerance
           (errors(2) < relative_tolerance * errors(3)) .or. &
           ! Change in the change less than the relative tolerance
           (abs(errors(2) - errors(1)) < relative_tolerance * errors(3))
    end function picard_converged

  end subroutine advance_nonlinear_term_picard