!!******************************************************************************
!!
!!  This file is part of the AMUN source code, a program to perform
!!  Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!!  adaptive mesh.
!!
!!  Copyright (C) 2008-2021 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License as published by
!!  the Free Software Foundation, either version 3 of the License, or
!!  (at your option) any later version.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: STATISTICS
!!
!!  This module provides subroutines to calculate and store mesh and variable
!!  statistics.
!!
!!******************************************************************************
!
module statistics

  implicit none

! MODULE PARAMETERS:
! =================
!
!   munit - a file handler for the mesh statistics file;
!   cunit - a file handler for the conserved variable integrals file;
!   sunit - a file handler for the variable statistics file;
!   eunit - a file handler for the integration errors file;
!
  integer(kind=4), save :: munit = 10
  integer(kind=4), save :: cunit = 11
  integer(kind=4), save :: sunit = 12
  integer(kind=4), save :: eunit = 13

!   iintd - the number of steps between subsequent intervals;
!
  integer(kind=4), save :: iintd = 1

! the mesh coverage and efficiency factors
!
  real(kind=8), save :: fcv = 1.0d+00, fef = 1.0d+00

! the format strings
!
  character(len=32), save :: mfmt1, mfmt2, mfmt3
  character(len=32), save :: efmt

  private

  public :: initialize_statistics, finalize_statistics
  public :: store_statistics

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
  contains
!
!===============================================================================
!!
!!***  PUBLIC SUBROUTINES  *****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine INITIALIZE_STATISTICS:
! -------------------------------
!
!   Subroutine initializes module INTEGRALS by preparing the statistics files.
!
!   Arguments:
!
!     verbose - the verbose flag;
!     nrun    - the number of resumed run;
!     status  - the subroutine call status;
!
!===============================================================================
!
  subroutine initialize_statistics(verbose, nrun, status)

    use coordinates, only : toplev, ncells, bcells, nghosts, domain_base_dims
    use equations  , only : pvars, nf
    use evolution  , only : error_control
    use mpitools   , only : master, nprocs
    use parameters , only : get_parameter

    implicit none

    logical, intent(in)  :: verbose
    integer, intent(in)  :: nrun
    integer, intent(out) :: status

    character(len=32) :: fname, stmp
    logical           :: append, flag
    integer           :: l

!-------------------------------------------------------------------------------
!
    status = 0

! only master process writes the files to the disk
!
    if (master) then

! process parameters
!
      call get_parameter("statistics_interval", iintd)
      iintd = max(1, iintd)

      stmp = "off"
      call get_parameter("statistics_append", stmp)
      select case(trim(stmp))
      case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
        append = .true.
      case default
        append = .false.
      end select

! prepare the mesh statistics file
!
      if (append) then
        write(fname,"('mesh-statistics.dat')")
        inquire(file=fname, exist=flag)
        flag = flag .and. nrun > 1
      else
        write(fname,"('mesh-statistics_',i2.2,'.dat')") nrun
        flag = .false.
      end if

      if (flag) then
#ifdef __INTEL_COMPILER
        open(newunit=munit, file=fname, form='formatted', status='old',        &
                                        position='append', buffered='yes')
#else /* __INTEL_COMPILER */
        open(newunit=munit, file=fname, form='formatted', status='old',        &
                                        position='append')
#endif /* __INTEL_COMPILER */
        write(munit,"('#')")
      else
#ifdef __INTEL_COMPILER
        open(newunit=munit, file=fname, form='formatted',                      &
                                        status='replace', buffered='yes')
#else /* __INTEL_COMPILER */
        open(newunit=munit, file=fname, form='formatted', status='replace')
#endif /* __INTEL_COMPILER */
      end if

! write the header of the mesh statistics file
!
      write(munit,"('#',5x,'step',4x,'time',11x,'coverage',7x,'efficiency',"// &
                  "5x,'blocks',5x,'leafs',1x,'|',1x," //                       &
                  "'block distribution per level')", advance="no")
#ifdef MPI
      write(stmp,"(a,i0,a)") "(", max(1, 10 * toplev - 28), "x,'|',1x,a)"
      write(munit,stmp) "block distribution per process"
#endif /* MPI */
      write(munit,"('#',75x,'|')", advance="no")
      do l = 1, toplev
        write(munit,"(1x,i9)", advance="no") l
      end do
#ifdef MPI
      write(stmp,"(a,i0,a)") "(", max(1, 10 * (3 - toplev)), "x,'|')"
      write(munit,stmp, advance="no")
      do l = 0, nprocs - 1
        write(munit,"(1x,i9)", advance="no") l
      end do
#endif /* MPI */
      write(munit,"('')")

! calculate the coverage and efficiency factors
!
      fcv = 1.0d+00 / (product(domain_base_dims) * 2**(NDIMS*(toplev - 1)))
      fef = 1.0d+00
      do l = 1, NDIMS
        fef = fef * (domain_base_dims(l) * ncells * 2**(toplev - 1)            &
                          + 2 * nghosts) / bcells
      end do

! prepare the format strings for mesh statistics
!
      mfmt1 = "(1x,i9,3(1x,1es14.6e3),2(1x,i9))"
      write (mfmt2,"(a,i0,a)") '(2x,', toplev, '(1x,i9))'
#ifdef MPI
      write(stmp,"(a,i0,a)") "(", max(2, 10 * (3 - toplev) + 1), "x,"
      write(mfmt3,"(a,i0,a)") trim(stmp), nprocs, '(1x,i9))'
#endif /* MPI */

! prepare the conserved variable integrals file
!
      if (append) then
        write(fname, "('variable-integrals.dat')")
        inquire(file=fname, exist=flag)
        flag = flag .and. nrun > 1
      else
        write(fname, "('variable-integrals_',i2.2,'.dat')") nrun
        flag = .false.
      end if

      if (flag) then
#ifdef __INTEL_COMPILER
        open(newunit=cunit, file=fname, form='formatted', status='old',        &
                                        position='append', buffered='yes')
#else /* __INTEL_COMPILER */
        open(newunit=cunit, file=fname, form='formatted', status='old',        &
                                        position='append')
#endif /* __INTEL_COMPILER */
        write(cunit,"('#')")
      else
#ifdef __INTEL_COMPILER
        open(newunit=cunit, file=fname, form='formatted',                      &
                                        status='replace', buffered='yes')
#else /* __INTEL_COMPILER */
        open(newunit=cunit, file=fname, form='formatted', status='replace')
#endif /* __INTEL_COMPILER */
      end if

      write(cunit,"('#',a8,13(1x,a18))") 'step', 'time', 'dt'                  &
                                       , 'mass', 'momx', 'momy', 'momz'        &
                                       , 'ener', 'ekin', 'emag', 'eint'        &
                                       , 'einj', 'erat', 'arms'
      write(cunit,"('#')")

! prepare the variable statistics file
!
      if (append) then
        write(fname, "('variable-statistics.dat')")
        inquire(file=fname, exist=flag)
        flag = flag .and. nrun > 1
      else
        write(fname, "('variable-statistics_',i2.2,'.dat')") nrun
        flag = .false.
      end if

      if (flag) then
#ifdef __INTEL_COMPILER
        open(newunit=sunit, file=fname, form='formatted', status='old',        &
                                        position='append', buffered='yes')
#else /* __INTEL_COMPILER */
        open(newunit=sunit, file=fname, form='formatted', status='old',        &
                                        position='append')
#endif /* __INTEL_COMPILER */
        write(sunit,"('#')")
      else
#ifdef __INTEL_COMPILER
        open(newunit=sunit, file=fname, form='formatted',                      &
                                        status='replace', buffered='yes')
#else /* __INTEL_COMPILER */
        open(newunit=sunit, file=fname, form='formatted', status='replace')
#endif /* __INTEL_COMPILER */
      end if

      write(sunit,"('#',a8,23(1x,a18))") 'step', 'time'                        &
                                 , 'mean(dens)', 'min(dens)', 'max(dens)'      &
                                 , 'mean(pres)', 'min(pres)', 'max(pres)'      &
                                 , 'mean(velo)', 'min(velo)', 'max(velo)'      &
                                 , 'mean(magn)', 'min(magn)', 'max(magn)'      &
                                 , 'mean(bpsi)', 'min(bpsi)', 'max(bpsi)'      &
                                 , 'mean(Mson)', 'min(Mson)', 'max(Mson)'      &
                                 , 'mean(Malf)', 'min(Malf)', 'max(Malf)'
      write(sunit,"('#')")

! prepare the integration errors file
!
      if (error_control) then

        if (append) then
          write(fname, "('integration-errors.dat')")
          inquire(file=fname, exist=flag)
          flag = flag .and. nrun > 1
        else
          write(fname, "('integration-errors_',i2.2,'.dat')") nrun
          flag = .false.
        end if

! check if the file exists; if not, create a new one, otherwise move to the end
!
        if (flag) then
#ifdef __INTEL_COMPILER
          open(newunit=eunit, file=fname, form='formatted', status='old',      &
                                          position='append', buffered='yes')
#else /* __INTEL_COMPILER */
          open(newunit=eunit, file=fname, form='formatted', status='old',      &
                                          position='append')
#endif /* __INTEL_COMPILER */
          write(cunit,"('#')")
        else
#ifdef __INTEL_COMPILER
          open(newunit=eunit, file=fname, form='formatted',                    &
                                          status='replace', buffered='yes')
#else /* __INTEL_COMPILER */
          open(newunit=eunit, file=fname, form='formatted', status='replace')
#endif /* __INTEL_COMPILER */
        end if

        write(stmp,"(i9)") nf + 4
        efmt = "('#',a8," // trim(adjustl(stmp)) // "(1x,a18))"
        write(eunit,efmt) 'step', 'time', 'dt_cfl', 'dt_err',                  &
                                          'max_err', pvars(1:nf)
        write(eunit,"('#')")

        efmt = "(i9," // trim(adjustl(stmp)) // "(1x,1es18.8e3))"

      end if ! error_control

    end if ! master

!-------------------------------------------------------------------------------
!
  end subroutine initialize_statistics
!
!===============================================================================
!
! subroutine FINALIZE_STATISTICS:
! -----------------------------
!
!   Subroutine finalizes module STATISTICS.
!
!   Arguments:
!
!     status  - an integer flag for error return value;
!
!===============================================================================
!
  subroutine finalize_statistics(status)

    use evolution, only : error_control
    use mpitools , only : master

    implicit none

    integer, intent(out) :: status

!-------------------------------------------------------------------------------
!
    status = 0

    if (master) then
      close(munit)
      close(cunit)
      close(sunit)
      if (error_control) close(eunit)
    end if

!-------------------------------------------------------------------------------
!
  end subroutine finalize_statistics
!
!===============================================================================
!
! subroutine STORE_STATISTICS:
! --------------------------
!
!   Subroutine calculates the processes the variable integrals and statistics,
!   mesh statistics and integration errors, collects them from all processes
!   and stores in the corresponding files.
!
!===============================================================================
!
  subroutine store_statistics()

    use blocks     , only : block_leaf, block_data
    use blocks     , only : list_leaf, list_data
    use blocks     , only : get_mblocks, get_nleafs
    use coordinates, only : ni => ncells, nb, ne
    use coordinates, only : advol, voli
    use coordinates, only : toplev
    use equations  , only : idn, ipr, ivx, ivz, ibx, iby, ibz, ibp
    use equations  , only : ien, imx, imy, imz
    use equations  , only : magnetized, adiabatic_index, csnd
    use equations  , only : errors
    use evolution  , only : error_control
    use evolution  , only : step, time, dt, dth, dte
    use forcing    , only : einj, rinj, arms
    use helpers    , only : print_message, flush_and_sync
    use mesh       , only : changed
    use mpitools   , only : master, nprocs
#ifdef MPI
    use mpitools       , only : reduce_minimum, reduce_maximum, reduce_sum
#endif /* MPI */
    use workspace      , only : resize_workspace, work, work_in_use

    implicit none

    logical, save :: first = .true.
    integer       :: n, l, u, nk, nl, nm, status
    real(kind=8)  :: dvol, dvolh

    type(block_leaf), pointer :: pleaf
    type(block_data), pointer :: pdata

    real(kind=8), dimension(:,:,:), pointer, save :: vel, mag, sqd, tmp

    integer         , parameter :: narr = 16

    real(kind=8), dimension(narr) :: inarr, avarr, mnarr, mxarr

    real(kind=8)    , parameter :: eps = epsilon(1.0d+00)
    real(kind=8)    , parameter :: big = huge(1.0d+00)

    integer(kind=4), dimension(toplev) :: ldist
#ifdef MPI
    integer(kind=4), dimension(nprocs) :: cdist
#endif /* MPI */

    character(len=*), parameter :: loc = 'INTEGRALS:store_statistics()'

!-------------------------------------------------------------------------------
!
! process and store the mesh statistics only on the master node
!
    if (master) then
      if (changed) then

! get the new numbers of meta blocks and leafs
!
        nm = get_mblocks()
        nl = get_nleafs()

! determine the block distribution per level and per process
!
        ldist(:) = 0
#ifdef MPI
        cdist(:) = 0
#endif /* MPI */
        pleaf => list_leaf
        do while(associated(pleaf))
          l = pleaf%meta%level
          ldist(l) = ldist(l) + 1
#ifdef MPI
          n = pleaf%meta%process+1
          cdist(n) = cdist(n) + 1
#endif /* MPI */
          pleaf => pleaf%next
        end do

! store the block statistics
!
        write (munit,mfmt1,advance='no') step, time, fcv * nl, fef / nl, nm, nl
        write (munit,mfmt2,advance='no') ldist(:)
#ifdef MPI
        write (munit,mfmt3,advance='no') cdist(:)
#endif /* MPI */
        write(munit,"('')")

        call flush_and_sync(munit)

        changed = .false.

      end if ! number of blocks or leafs changed
    end if ! master

! return if the storage interval was not reached
!
    if (mod(step, iintd) > 0) return

    if (first) then
      n = 4 * ni**NDIMS

      call resize_workspace(n, status)
      if (status /= 0) then
        call print_message(loc, "Could not resize the workspace!")
        go to 100
      end if

#if NDIMS == 3
      nk = ni
#else /* NDIMS == 3 */
      nk = 1
#endif /* NDIMS == 3 */

      n = ni**NDIMS
      l = 1
      u = n
      vel(1:ni,1:ni,1:nk) => work(l:u)
      l = l + n
      u = u + n
      mag(1:ni,1:ni,1:nk) => work(l:u)
      l = l + n
      u = u + n
      sqd(1:ni,1:ni,1:nk) => work(l:u)
      l = l + n
      u = u + n
      tmp(1:ni,1:ni,1:nk) => work(l:u)

      first = .false.
    end if

! reset the integrals array
!
    inarr(:) =   0.0d+00
    avarr(:) =   0.0d+00
    mnarr(:) =   big
    mxarr(:) = - big

! reset some statistics if they are not used
!
    if (ipr < 1) then
      mnarr(2) = 0.0d+00
      mxarr(2) = 0.0d+00
    end if
    if (.not. magnetized) then
      mnarr(4) = 0.0d+00
      mxarr(4) = 0.0d+00
      mnarr(5) = 0.0d+00
      mxarr(5) = 0.0d+00
      mnarr(7) = 0.0d+00
      mxarr(7) = 0.0d+00
    end if

    if (work_in_use) &
      call print_message(loc, "Workspace is being used right now! " //         &
                              "Corruptions can occur!")

    work_in_use = .true.

! associate the pointer with the first block on the data block list
!
    pdata => list_data

! iterate over all data blocks on the list
!
    do while(associated(pdata))

! obtain the volume elements for the current block
!
      dvol  = advol(pdata%meta%level)
      dvolh = 0.5d+00 * dvol

! sum up density and momenta components
!
#if NDIMS == 3
      inarr(1) = inarr(1) + sum(pdata%u(idn,nb:ne,nb:ne,nb:ne)) * dvol
      inarr(2) = inarr(2) + sum(pdata%u(imx,nb:ne,nb:ne,nb:ne)) * dvol
      inarr(3) = inarr(3) + sum(pdata%u(imy,nb:ne,nb:ne,nb:ne)) * dvol
      inarr(4) = inarr(4) + sum(pdata%u(imz,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
      inarr(1) = inarr(1) + sum(pdata%u(idn,nb:ne,nb:ne,  :  )) * dvol
      inarr(2) = inarr(2) + sum(pdata%u(imx,nb:ne,nb:ne,  :  )) * dvol
      inarr(3) = inarr(3) + sum(pdata%u(imy,nb:ne,nb:ne,  :  )) * dvol
      inarr(4) = inarr(4) + sum(pdata%u(imz,nb:ne,nb:ne,  :  )) * dvol
#endif /* NDIMS == 3 */

! sum up total energy
!
      if (ien > 0) then
#if NDIMS == 3
        inarr(5) = inarr(5) + sum(pdata%u(ien,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
        inarr(5) = inarr(5) + sum(pdata%u(ien,nb:ne,nb:ne,  :  )) * dvol
#endif /* NDIMS == 3 */
      end if

! sum up kinetic energy
!
#if NDIMS == 3
      inarr(6) = inarr(6) + sum((pdata%u(imx,nb:ne,nb:ne,nb:ne)**2             &
                               + pdata%u(imy,nb:ne,nb:ne,nb:ne)**2             &
                               + pdata%u(imz,nb:ne,nb:ne,nb:ne)**2)            &
                               / pdata%u(idn,nb:ne,nb:ne,nb:ne)) * dvolh
#else /* NDIMS == 3 */
      inarr(6) = inarr(6) + sum((pdata%u(imx,nb:ne,nb:ne,  :  )**2             &
                               + pdata%u(imy,nb:ne,nb:ne,  :  )**2             &
                               + pdata%u(imz,nb:ne,nb:ne,  :  )**2)            &
                               / pdata%u(idn,nb:ne,nb:ne,  :  )) * dvolh
#endif /* NDIMS == 3 */

! sum up magnetic energy
!
      if (magnetized) then
#if NDIMS == 3
        inarr(7) = inarr(7) + sum(pdata%u(ibx,nb:ne,nb:ne,nb:ne)**2            &
                                + pdata%u(iby,nb:ne,nb:ne,nb:ne)**2            &
                                + pdata%u(ibz,nb:ne,nb:ne,nb:ne)**2) * dvolh
#else /* NDIMS == 3 */
        inarr(7) = inarr(7) + sum(pdata%u(ibx,nb:ne,nb:ne,  :  )**2            &
                                + pdata%u(iby,nb:ne,nb:ne,  :  )**2            &
                                + pdata%u(ibz,nb:ne,nb:ne,  :  )**2) * dvolh
#endif /* NDIMS == 3 */
      end if

! sum up the injected energy and injection rate
!
      inarr( 9) = einj
      inarr(10) = rinj

! get average, minimum and maximum values of density
!
#if NDIMS == 3
      tmp(:,:,:) = pdata%q(idn,nb:ne,nb:ne,nb:ne)
#else /* NDIMS == 3 */
      tmp(:,:,:) = pdata%q(idn,nb:ne,nb:ne,  :  )
#endif /* NDIMS == 3 */
      avarr(1) = avarr(1) + sum(tmp(:,:,:)) * dvol
      mnarr(1) = min(mnarr(1), minval(tmp(:,:,:)))
      mxarr(1) = max(mxarr(1), maxval(tmp(:,:,:)))

! get average, minimum and maximum values of pressure
!
      if (ipr > 0) then
#if NDIMS == 3
        tmp(:,:,:) = pdata%q(ipr,nb:ne,nb:ne,nb:ne)
#else /* NDIMS == 3 */
        tmp(:,:,:) = pdata%q(ipr,nb:ne,nb:ne,  :  )
#endif /* NDIMS == 3 */
        avarr(2) = avarr(2) + sum(tmp(:,:,:)) * dvol
        mnarr(2) = min(mnarr(2), minval(tmp(:,:,:)))
        mxarr(2) = max(mxarr(2), maxval(tmp(:,:,:)))
      end if

! get average, minimum and maximum values of velocity amplitude
!
#if NDIMS == 3
      vel(:,:,:) = sqrt(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne)**2, 1))
#else /* NDIMS == 3 */
      vel(:,:,:) = sqrt(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,  :  )**2, 1))
#endif /* NDIMS == 3 */
      avarr(3) = avarr(3) + sum(vel(:,:,:)) * dvol
      mnarr(3) = min(mnarr(3), minval(vel(:,:,:)))
      mxarr(3) = max(mxarr(3), maxval(vel(:,:,:)))

! get average, minimum and maximum values of magnetic field amplitude, and
! divergence potential
!
      if (magnetized) then
#if NDIMS == 3
        mag(:,:,:) = sqrt(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne)**2, 1))
#else /* NDIMS == 3 */
        mag(:,:,:) = sqrt(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,  :  )**2, 1))
#endif /* NDIMS == 3 */
        avarr(4) = avarr(4) + sum(mag(:,:,:)) * dvol
        mnarr(4) = min(mnarr(4), minval(mag(:,:,:)))
        mxarr(4) = max(mxarr(4), maxval(mag(:,:,:)))

#if NDIMS == 3
        tmp(:,:,:) = pdata%q(ibp,nb:ne,nb:ne,nb:ne)
#else /* NDIMS == 3 */
        tmp(:,:,:) = pdata%q(ibp,nb:ne,nb:ne,  :  )
#endif /* NDIMS == 3 */
        avarr(5) = avarr(5) + sum(tmp(:,:,:)) * dvol
        mnarr(5) = min(mnarr(5), minval(tmp(:,:,:)))
        mxarr(5) = max(mxarr(5), maxval(tmp(:,:,:)))
      end if

! calculate square root of density
!
#if NDIMS == 3
      sqd(:,:,:) = sqrt(pdata%q(idn,nb:ne,nb:ne,nb:ne))
#else /* NDIMS == 3 */
      sqd(:,:,:) = sqrt(pdata%q(idn,nb:ne,nb:ne,  :  ))
#endif /* NDIMS == 3 */

! get average, minimum and maximum values of sonic Mach number
!
      if (ipr > 0) then
        tmp(:,:,:) = sqd(:,:,:) * vel(:,:,:) /                                 &
#if NDIMS == 3
                     sqrt(adiabatic_index * pdata%q(ipr,nb:ne,nb:ne,nb:ne))
#else /* NDIMS == 3 */
                     sqrt(adiabatic_index * pdata%q(ipr,nb:ne,nb:ne,  :  ))
#endif /* NDIMS == 3 */
      else
        tmp(:,:,:) = vel(:,:,:) / csnd
      end if
      avarr(6) = avarr(6) + sum(tmp(:,:,:)) * dvol
      mnarr(6) = min(mnarr(6), minval(tmp(:,:,:)))
      mxarr(6) = max(mxarr(6), maxval(tmp(:,:,:)))

! get average, minimum and maximum values of Alfvénic Mach number
!
      if (magnetized) then
        tmp(:,:,:) = sqd(:,:,:) * vel(:,:,:) / max(eps, mag(:,:,:))
        avarr(7) = avarr(7) + sum(tmp(:,:,:)) * dvol
        mnarr(7) = min(mnarr(7), minval(tmp(:,:,:)))
        mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:)))
      end if

! associate the pointer with the next block on the list
!
      pdata => pdata%next

    end do ! data blocks

    work_in_use = .false.

#ifdef MPI
! sum the integral array from all processes
!
    call reduce_sum(inarr(1:narr))

! reduce average, minimum and maximum values
!
    call reduce_sum(avarr(1:narr))
    call reduce_minimum(mnarr(1:narr))
    call reduce_maximum(mxarr(1:narr))
#endif /* MPI */

! calculate the internal energy
!
    if (ien > 0) inarr(8) = inarr(5) - inarr(6) - inarr(7)

! normalize the averages by the volume of domain
!
    avarr(:) = avarr(:) * voli

! write down the integrals and statistics to appropriate files
!
    if (master) then
      write(cunit,"(i9,13(1x,1es18.8e3))") step, time, dt, inarr(1:10), arms
      write(sunit,"(i9,23(1x,1es18.8e3))") step, time                          &
                                            , avarr(1), mnarr(1), mxarr(1)     &
                                            , avarr(2), mnarr(2), mxarr(2)     &
                                            , avarr(3), mnarr(3), mxarr(3)     &
                                            , avarr(4), mnarr(4), mxarr(4)     &
                                            , avarr(5), mnarr(5), mxarr(5)     &
                                            , avarr(6), mnarr(6), mxarr(6)     &
                                            , avarr(7), mnarr(7), mxarr(7)

      call flush_and_sync(cunit)
      call flush_and_sync(sunit)

      if (error_control) then
        write(eunit,efmt) step, time, dth, dte, maxval(errors(:)), errors(:)
        call flush_and_sync(eunit)
      end if
    end if

    100 continue

!-------------------------------------------------------------------------------
!
  end subroutine store_statistics

!===============================================================================
!
end module