avg_hk Subroutine

  subroutine avg_hk (hk, hk_hist, istep, navg)
    use mp, only: proc0
    use species, only: nspec
    implicit none
    !> Heating diagnostics as a function of \((\theta, k_y)\) at the current timestep
    type (heating_diagnostics), dimension(:,:), intent(in out) :: hk
    !> Heating diagnostics as a function of \((\theta, k_y)\) over the last [[gs2_diagnostics:navg]] timesteps
    type (heating_diagnostics), dimension(:,:,0:), intent(in out) :: hk_hist
    integer, intent (in) :: istep, navg
    integer :: is, i, m, n, mmax, nmax

    mmax = size(hk, 1)
    nmax = size(hk, 2)

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

  end subroutine avg_hk