avg_h Subroutine

  subroutine avg_h (h, h_hist, istep, navg)
    use mp, only: proc0
    use species, only: nspec
    implicit none
    !> The heating diagnostics at the current timestep
    type (heating_diagnostics), intent(in out) :: h
    !> The last [[gs2_diagnostics:navg]] timesteps of heating diagnostics
    type (heating_diagnostics), dimension (0:), intent(in out) :: h_hist
    integer, intent (in) :: istep, navg
    integer :: is, i

    if (proc0) then
       if (navg > 1) then
          if (istep > 1) then
             h_hist(mod(istep,navg)) % energy     = h % energy
             h_hist(mod(istep,navg)) % energy_dot = h % energy_dot
             h_hist(mod(istep,navg)) % antenna    = h % antenna
             h_hist(mod(istep,navg)) % eapar      = h % eapar
             h_hist(mod(istep,navg)) % ebpar      = h % ebpar
             h_hist(mod(istep,navg)) % delfs2     = h % delfs2
             h_hist(mod(istep,navg)) % hs2        = h % hs2
             h_hist(mod(istep,navg)) % phis2      = h % phis2
             h_hist(mod(istep,navg)) % hypervisc  = h % hypervisc
             h_hist(mod(istep,navg)) % hyperres   = h % hyperres
             h_hist(mod(istep,navg)) % hypercoll  = h % hypercoll
             h_hist(mod(istep,navg)) % collisions = h % collisions
             h_hist(mod(istep,navg)) % imp_colls  = h % imp_colls
             h_hist(mod(istep,navg)) % gradients  = h % gradients
!             h_hist(mod(istep,navg)) % curvature  = h % curvature
             h_hist(mod(istep,navg)) % heating    = h % heating
          end if
          
          if (istep >= navg) then
             call zero_htype(h)
             do i=0,navg-1
                h % energy     = h%energy     + h_hist(i) % energy / real(navg)
                h % energy_dot = h%energy_dot + h_hist(i) % energy_dot / real(navg)
                h % antenna    = h%antenna    + h_hist(i) % antenna / real(navg)
                h % eapar      = h%eapar      + h_hist(i) % eapar    / real(navg)
                h % ebpar      = h%ebpar      + h_hist(i) % ebpar    / real(navg)
                do is = 1,nspec
                   h % delfs2(is)     = h%delfs2(is)     + h_hist(i) % delfs2(is) / real(navg)
                   h % hs2(is)        = h%hs2(is)        + h_hist(i) % hs2(is) / real(navg)
                   h % phis2(is)      = h%phis2(is)      + h_hist(i) % phis2(is) / real(navg)
                   h % hypervisc(is)  = h%hypervisc(is)  + h_hist(i) % hypervisc(is) / real(navg)
                   h % hyperres(is)   = h%hyperres(is)   + h_hist(i) % hyperres(is) / real(navg)
                   h % hypercoll(is)  = h%hypercoll(is)  + h_hist(i) % hypercoll(is) / real(navg)
                   h % collisions(is) = h%collisions(is) + h_hist(i) % collisions(is) / real(navg)
                   h % imp_colls(is) = h%imp_colls(is) + h_hist(i) % imp_colls(is) / real(navg)
                   h % gradients(is)  = h%gradients(is)  + h_hist(i) % gradients(is) / real(navg)
!                  h % curvature(is)  = h%curvature(is)  + h_hist(i) % curvature(is) / real(navg)
                   h % heating(is)    = h%heating(is)    + h_hist(i) % heating(is) / real(navg)
                end do
             end do
          end if
       end if
    end if

  end subroutine avg_h