get_omegaavg Subroutine

  subroutine get_omegaavg (istep, exit, omegaavg, debopt)
    use kt_grids, only: ntheta0, naky
    use gs2_time, only: code_dt
    use optionals, only: get_option_with_default
    implicit none
    !> The current timestep
    integer, intent (in) :: istep
    !> Should the simulation exit?
    !> FIXME: This could be `intent(out)` only
    logical, intent (in out) :: exit
    !> The time-averaged complex frequency
    complex, dimension (:,:), intent (out) :: omegaavg
    !> If true, write some debug messages
    logical, intent (in), optional :: debopt
    logical :: debug
    debug = get_option_with_default(debopt, .false.)
    if (debug) write(6,*) "get_omeaavg: allocate domega"
    if (.not. allocated(domega)) allocate(domega(navg,ntheta0,naky))
    if (debug) write(6,*) "get_omeaavg: start"
    if (debug) write(6,*) "get_omeaavg: set omegahist"

    ! Get and store the current instantaneous omega
    omegahist(mod(istep, navg), :, :) = calculate_instantaneous_omega(igomega, omegatol)

    if (debug) write(6,*) "get_omeaavg: set omegahist at istep = 0"
    !During initialisation fieldnew==field but floating error can lead to finite omegahist
    !Force omegahist=0 here to avoid erroneous values.
    !Could think about forcing omegahist=0 where abs(omegahist)<tol
    if(istep==0) omegahist(:,:,:)=0.0

    if (debug) write(6,*) "get_omeaavg: set omegaavg"
    omegaavg = sum(omegahist/real(navg),dim=1)
    if (debug) write(6,*) "get_omegaavg: omegaavg=",omegaavg
    if (istep > navg) then
       domega = spread(omegaavg,1,navg) - omegahist
       if (all(sqrt(sum(abs(domega)**2/real(navg),dim=1)) &
            .le. min(abs(omegaavg),1.0)*omegatol)) &
            then
          if (write_ascii) write (out_unit, "('*** omega converged')")
          exit = .true. .and. exit_when_converged
       end if

       if (any(abs(omegaavg)*code_dt > omegatinst)) then
          if (write_ascii) write (out_unit, "('*** numerical instability detected')") 
          exit = .true.
       end if
    end if
    if (debug) write(6,*) "get_omegaavg: done"
  end subroutine get_omegaavg