do_write_heating Subroutine

  subroutine do_write_heating(t, file_unit, file_unit2, h)
    use species, only: nspec
    implicit none
    real, intent(in) :: t
    integer, intent(in) :: file_unit, file_unit2
    !> Heating diagnostics
    type(heating_diagnostics), intent (in) :: h

    integer :: is

    !
    ! For case with two species:
    !
    ! Column     Item               
    !   1        time              
    !   2        Energy              
    !   3        dEnergy/dt            
    !   4        J_ant.E             
    !   5        [h_(i+1)*h_*]/2 * C[h_(i+1)] * T_0 for species 1
    !   6        [h_(i+1)*h_*]/2 * C[h_(i+1)] * T_0 for species 2
    !   7       -[h H(h) * T_0]_1
    !   8       -[h H(h) * T_0]_2
    !   9       -[h C(h) * T_0]_1 
    !  10       -[h C(h) * T_0]_2
    !  11        [h w_* h]_1
    !  12        [h w_* h]_2
    !  13        [h * (q dchi/dt - dh/dt * T0)]_1
    !  14        [h * (q dchi/dt - dh/dt * T0)]_2
    !  15      sum (h C(h) * T_0)  in total, as in 5, 6      
    !  16     -sum (h H(h) * T_0)      
    !  17     -sum (h C(h) * T_0)   
    !  18      sum (h w_* h)  
    !  19      sum [h (q dchi/dt - dh/dt * T0)]
    !  20      3 + 4 + 18 + 19
    !  21      (k_perp A)**2
    !  22      B_par**2
    !  23      df_1 ** 2
    !  24      df_2 ** 2
    !  25      h_1 ** 2
    !  26      h_2 ** 2
    !  27      Phi_bar_1 ** 2
    !  28      Phi_bar_2 ** 2
    !
    !
    ! For case with one species:
    !
    ! Column     Item               
    !   1        time              
    !   2        Energy              
    !   3        dEnergy/dt            
    !   4        J_ant.E             
    !   5        [h_(i+1)*h_*]/2 * C[h_(i+1)] * T_0 
    !   6       -[h H(h) * T_0]
    !   7       -[h C(h) * T_0]
    !   8        [h w_* h]
    !   9        [h * (q dchi/dt - dh/dt * T0)]_1
    !  10      sum (h C(h) * T_0)  in total, as in 5, 6      
    !  11     -sum (h H(h) * T_0)      
    !  12     -sum (h C(h) * T_0)   
    !  13      sum (h w_* h)  
    !  14      sum [h (q dchi/dt - dh/dt * T0)]
    !  15      3 + 4 + 9 + 10
    !  16      (k_perp A)**2
    !  17      B_par**2
    !  18      df ** 2
    !  19      h ** 2
    !  20      Phi_bar ** 2
    
    write (unit=file_unit, fmt="(28es12.4)") t,h % energy,  &
         h % energy_dot, h % antenna, h % imp_colls, h % hypercoll, h % collisions, &
         h % gradients, h % heating, sum(h % imp_colls), sum(h % hypercoll), sum(h % collisions), &
         sum(h % gradients), sum(h % heating),sum(h%heating)+h%antenna+sum(h%gradients)+h%energy_dot, &
         h % eapar, h % ebpar, h % delfs2(:),  h % hs2(:), h % phis2(:)

    do is=1,nspec
      write (unit=file_unit2, fmt="(15es12.4)") t,h % energy,  &
           h % energy_dot, h % antenna, h % imp_colls(is), h % hypercoll(is), h % collisions(is), &
           h % gradients(is), h % heating(is), &
           h % eapar, h % ebpar, h % delfs2(is),  h % hs2(is), h % phis2(is), real(is)
      write (unit=file_unit2, fmt=*)
    end do
    write (unit=file_unit2, fmt=*)

    flush (file_unit)
    flush (file_unit2)
  end subroutine do_write_heating