calculate_omega Subroutine

  subroutine calculate_omega (gnostics)
    use fields_arrays, only: phi, apar, bpar, phinew, aparnew, bparnew
    use gs2_time, only: code_dt
    use constants, only: zi
    use diagnostics_config, only: diagnostics_type
    use run_parameters, only: ieqzip
    use nonlinear_terms, only: nonlin
    use mp, only: proc0
    use warning_helpers, only: is_zero
    implicit none
    type(diagnostics_type), intent(inout) :: gnostics
    real :: fac
    integer :: j

    ! Avoid calculation of omega in nonlinear calculations.
    ! Note we rely on omega_average and related arrays being
    ! initialised elsewhere.
    if (nonlin .and. .not. any(ieqzip)) then
       return
    end if

    if (debug) write(6,*) "get_omeaavg: start"

    j = gnostics%igomega
    if (is_zero(gnostics%omegatol)) then
       fac = 1.0
    else
       fac = 1000 / abs(gnostics%omegatol)
    end if
    
    !<DD>The logic below was originally
    !    where (abs(phinew(j,:,:)+aparnew(j,:,:)+bparnew(j,:,:)) < epsilon(0.0) &
    !           .or. abs(phi(j,:,:)+apar(j,:,:)+bpar(j,:,:)) < epsilon(0.0))
    !This results in not calculating the frequency whenever the fields drop below epsilon(0.0) [~1d-16 for DP]
    !What we really want to do is prevent dividing by too small a number. 
    ! Calculate omega
    where (abs(phinew(j,:,:)+aparnew(j,:,:)+bparnew(j,:,:)) < tiny(0.0)*fac &
         .or. abs(phi(j,:,:)+apar(j,:,:)+bpar(j,:,:)) < tiny(0.0)*fac)
       omegahist(mod(gnostics%istep,gnostics%navg),:,:) = 0.0
    elsewhere
       omegahist(mod(gnostics%istep,gnostics%navg),:,:) &
            = log((phinew(j,:,:) + aparnew(j,:,:) + bparnew(j,:,:)) &
            /(phi(j,:,:)   + apar(j,:,:)    + bpar(j,:,:)))*zi/code_dt
    end where
    
    ! Calculate omega_average
    omega_average = sum(omegahist/real(gnostics%navg),dim=1)
    ! Copy the results to gnostics
    gnostics%current_results%omega_average = omega_average
    
    if (debug) write(6,*) "calculate_omega: omega_average=",omega_average
    
    ! Check if converged by comparing the difference between
    ! omega_average and all values in omegahist, to the value of omega_average
    if (gnostics%istep > gnostics%navg) then
       domega = spread(omega_average,1,gnostics%navg) - omegahist
       if (all(sqrt(sum(abs(domega)**2/real(gnostics%navg),dim=1)) &
            .le. min(abs(omega_average),1.0)*gnostics%omegatol)) &
       then
          if (proc0 .and. gnostics%ascii_files%write_to_out) write (gnostics%ascii_files%out, "('*** omega converged')")
          if (gnostics%exit_when_converged) gnostics%exit = .true.
       end if
       
       if (any(abs(omega_average)*code_dt > gnostics%omegatinst)) then
          !In old diagnostics this would only be written by proc0 and it would
          !go to gnostics%ascii_files%out rather than screen. Should we make this
          !consistent here?
          write (*, "('*** numerical instability detected')") 
          gnostics%exit = .true.
       end if
    end if
    
    if (debug) write(6,*) "calculate_omega: done"
  end subroutine calculate_omega