! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2011-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2011 Jan Hueckelheim <j.hueckelheim@grs-sim.de>
! Copyright (c) 2012-2014, 2016-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2012-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2012, 2021 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2012 Sathish Krishnan P S <s.krishnan@grs-sim.de>
! Copyright (c) 2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016 Philipp Otte <otte@mathcces.rwth-aachen.de>
! Copyright (c) 2018 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2019, 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2021 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@dlr.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! **************************************************************************** !
!> This module keeps all information about the fluid
!!
!! In this module, all infos about the fluid is collected.
!! The file type [[mus_fluid_type]] contains all relevant information.
!! This includes physical parameters such as the viscosity and reference
!! density. Also, LBM-specific parameters such as the relaxation rates are
!! defined.
!!
module mus_fluid_module
  use mpi
  ! include treelm modules
  use env_module,           only: rk
  use tem_tools_module,     only: tem_horizontalSpacer
  use tem_aux_module,       only: tem_abort
  use tem_param_module,     only: cs2, cs2inv, rho0, div2_3
  use tem_logging_module,   only: logUnit
  use tem_time_module,      only: tem_time_type
  use tem_stencil_module,   only: tem_stencilHeader_type
  use tem_general_module,   only: tem_general_type
  use tem_spacetime_fun_module, only: tem_load_spacetime
  use tem_construction_module,  only: tem_levelDesc_type
  use tem_grow_array_module,    only: init

  ! include aotus modules
  use aotus_module,     only: flu_State,          &
    &                         aot_get_val,        &
    &                         aoterr_NonExistent, &
    &                         aoterr_Fatal
  use aot_table_module, only: aot_table_open, aot_table_close
  use aot_out_module,   only: aot_out_type, aot_out_val, aot_out_open_table,   &
    &                         aot_out_close_table

  ! include musubi modules
  use mus_physics_module,       only: mus_physics_type
  use mus_nonNewtonian_module,  only: mus_nNwtn_type, mus_nNwtn_load, &
    &                                 mus_assign_nNwtnVisc_ptr,       &
    &                                 mus_nNwtn_dump2outUnit, mus_nNwtn_save2lua
  use mus_pdf_module,            only: pdf_data_type
  use mus_turbulence_module,     only: mus_turbulence_type,  &
    &                                  mus_load_turbulence,  &
    &                                  mus_init_turbulenceData
  use mus_turb_viscosity_module, only: mus_assign_turbVisc_ptr
  use mus_scheme_header_module,  only: mus_scheme_header_type
  use mus_relaxationParam_module, only: mus_viscosity_type,              &
    &                                   mus_init_relaxParam,             &
    &                                   mus_calcOmegaFromVisc,           &
    &                                   mus_update_relaxParamFromViscSTfun
  use mus_mrtRelaxation_module,   only: mus_proc_mrt,                      &
    &                                   mus_assign_mrt_ptr
  use mus_cumulantInit_module,    only: cumulant_omega_check

  implicit none

  private

  public :: mus_fluid_type
  public :: mus_load_fluid
  public :: mus_fluid_save2lua
  public :: mus_init_fluid
  public :: mus_fluid_cleanup

  !> collection of properties of the fluid
  type mus_fluid_type

    logical :: active = .false.       !< is this object a fluid?

    !> Magic value for TRT collision model
    !! Lambda = ( 1/omega_+ - 0.5 ) * ( 1/omega_- - 0.5 )
    real(kind=rk) :: lambda = 0.25_rk

    !> level-wise bulk omegas, used as relaxation in mrt model
    !! allocated in mus_init_fluid
    real(kind=rk), allocatable :: omegaBulkLvl(:)
    !> Level wise bulk viscosity in lattice
    real(kind=rk), allocatable :: viscBulkLvl(:)

    !> Contains information for turbulence model
    type(mus_turbulence_type) :: turbulence

    !> nonNewtonian fluid parameter
    type(mus_nNwtn_type) :: nNwtn

    !> function pointer to get MRT diagonal relaxation matrix
    procedure(mus_proc_mrt), nopass, pointer :: mrtPtr => null()

    !> kinematic viscosity
    !! \todo KM: implement interpolation routine for constant viscosity
    type(mus_viscosity_type) :: viscKine

    real(kind=rk) :: viscBulk_phy  !< physical bulk viscosity

    real(kind=rk) :: force(3)        !< forcing term

    ! also works around ICE in Intel 19.1
    ! HRR
    real(kind=rk) :: HRR_sigma = 0.98_rk

    ! DRT
    real(kind=rk) :: DRT_tauN = 0.70_rk

    ! Cumulant extended
    real(kind=rk) :: omega_Cum(10)
    real(kind=rk) :: omega_Lim(3)

  end type mus_fluid_type


contains


! **************************************************************************** !
  !> Read in the fluid property from the LUA parameter file
  subroutine mus_load_fluid( me, conf, parent, minLevel, physics, schemeHeader )
    !--------------------------------------------------------------------------
    !> fluid type
    type( mus_fluid_type ),intent(out) :: me
    !> lua state
    type( flu_state )              :: conf
    !> parent handle
    integer, intent(in), optional  :: parent
    !> global pdf info
    integer, intent(in) :: minLevel
    !> physics type to convert physics to lattice unit or vice versa
    type( mus_physics_type ), intent(in) :: physics
    !> identifier of the scheme
    type( mus_scheme_header_type ), intent(in) :: schemeHeader
    ! --------------------------------------------------------------------------
    integer :: fluid_table, ii
    integer :: iError
    integer :: vecError(3)
    real(kind=rk) :: omegaBulk_def ! default omegaBulk for incompressible model
    character(len=6) :: omega_name = "omega_"
    character(len=8) :: omega_string
    character(len=10) :: omega_lim_name = "omega_lim_"
    character(len=11) :: omega_lim_string
    real(kind=rk) :: omegaCum_def
    ! --------------------------------------------------------------------------


    ! if fluid informations in scheme table parentHandle /= 0
    if( present(parent) ) then
      call aot_table_open( L       = conf,                                     &
        &                  parent  = parent,                                   &
        &                  thandle = fluid_table,                              &
        &                  key     = 'fluid' )
    else
      call aot_table_open( L=conf, thandle = fluid_table, key = 'fluid' )
    end if

    if ( fluid_table == 0 ) then
      write(logUnit(1),*)'No fluid table defined'
      call tem_abort()
    endif

    write(logUnit(1),*) 'Loading fluid informations'
    ! Activating the current fluid object
    me%active = .true.

    ! load kinematic viscosity as spacetime function
    call tem_load_spacetime( me      = me%viscKine%STfun,      &
      &                      conf    = conf,                   &
      &                      parent  = fluid_table,            &
      &                      key     = 'kinematic_viscosity',  &
      &                      errCode = iError                  )

    if (iError /= 0) then
      call tem_abort('ERROR: Unable to load "kinematic_viscosity"' &
        &          //'from fluid table')
    end if

    ! Load bulk viscosity as scalar and convert it to LB
    call aot_get_val( L       = conf,             &
      &               thandle = fluid_table,      &
      &               key     = 'bulk_viscosity', &
      &               val     = me%viscBulk_phy,  &
      &               ErrCode = iError            )

    ! If not provided, use default value only for incompressible model
    if (btest(iError, aoterr_NonExistent)) then
      if (trim(schemeHeader%kind) == 'fluid') then
        if (trim(schemeHeader%layout) == 'd2q9') then
          omegaBulk_def = 1.63_rk
          write(logUnit(1),*) 'Warning: Using default omegaBulk=1.63'
        else if (trim(schemeHeader%layout) == 'd3q27') then
          omegaBulk_def = 1.54_rk
          write(logUnit(1),*) 'Warning: Using default omegaBulk=1.54'
        else
          call tem_abort('ERROR: "bulk_viscosity" is not provided in '// &
            &            'fluid table')
        end if
      else
        ! using default value from
        ! D’Humières, D., Ginzburg, I., Krafczyk, M., Lallemand, P.,
        ! & Luo, L.-S. (2002).
        ! Multiple-relaxation-time lattice Boltzmann models in three
        ! dimensions. Philosophical Transactions. Series A, Mathematical,
        ! Physical, and Engineering Sciences, 360(1792), 437–51.
        select case(trim(schemeHeader%layout))
        case ('d2q9')
          ! default value for d2q9 is from
          ! Chen, S., Peng, C., Teng, Y., Wang, L. P., & Zhang, K. (2016).
          ! Improving lattice Boltzmann simulation of moving particles in a
          ! viscous flow using local grid refinement. Computers and Fluids,
          omegaBulk_def = 1.63_rk
        case ('d3q15')
          omegaBulk_def = 1.60_rk
        case ('d3q27')
          omegaBulk_def = 1.54_rk
        case default
          omegaBulk_def = 1.19_rk
        end select

        write(logUnit(1),*) 'Warning: Using default omegaBulk:', omegaBulk_def
      end if

      if (trim(schemeHeader%layout) == 'd2q9') then
        !! "Theory of the lattice Boltzmann method: Dispersion, dissipation,
        !! isotropy, Galilean invariance, and stability", Pierre Lallemand and
        !! Li-Shi Luo, Phys. Rev. E 61, 2000.
        me%viscBulk_phy = ((8._rk - 12._rk * cs2) / 24._rk) &
          &             * (1.0_rk/omegaBulk_def - 0.5_rk )  &
          &             * physics%fac(minLevel)%visc
      else
        me%viscBulk_phy = ((5._rk - 9._rk * cs2) / 9._rk)  &
          &             * (1.0_rk/omegaBulk_def - 0.5_rk ) &
          &             * physics%fac(minLevel)%visc
      endif

    endif

    ! load lambda (magic value for trt model)
    if (schemeHeader%relaxation(1:3) == 'trt') then
      call aot_get_val(L       = conf,        &
        &              thandle = fluid_table, &
        &              key     = 'lambda',    &
        &              val     = me%lambda,   &
        &              default = 0.25_rk,     &
        &              ErrCode = iError       )
    end if

    ! load forces
    call aot_get_val( L       = conf,                &
      &               thandle = fluid_table,         &
      &               key     = 'force',             &
      &               val     = me%force,            &
      &               default = [0._rk,0._rk,0._rk], &
      &               ErrCode = vecError             )

    ! convert force to lattice
    me%force = me%force / physics%fac(minLevel)%force

    ! load sigma for HRR_bgk
    if (trim(schemeHeader%relaxation) == 'hrr_bgk') then
      call aot_get_val(L       = conf,        &
        &              thandle = fluid_table, &
        &              key     = 'hrr_sigma',     &
        &              val     = me%HRR_sigma,&
        &              default = 0.98_rk,     &
        &              ErrCode = iError       )

      if (btest(iError, aoterr_fatal)) then
        call tem_abort('ERROR: wrong data type for sigma_HRR ' &
            &          // 'fluid table')
      end if
    end if
    ! load tauN for DRT_bgk
    if (trim(schemeHeader%relaxation) == 'drt_bgk') then
      call aot_get_val(L       = conf,        &
        &              thandle = fluid_table, &
        &              key     = 'drt_taun',     &
        &              val     = me%DRT_tauN,&
        &              default = 0.70_rk,     &
        &              ErrCode = iError       )

      if (btest(iError, aoterr_fatal)) then
        call tem_abort('ERROR: wrong data type for DRT tauN ' &
            &          // 'fluid table')
      end if
    end if

    ! load omega vec for cumulant and limiter
    if (trim(schemeHeader%relaxation) == 'cumulant_extended' .or. &
      & trim(schemeHeader%relaxation) == 'cumulant_extended_generic') then
      me%omega_Cum = -1._rk

      call aot_get_val(L       = conf,            &
        &              thandle = fluid_table,     &
        &              key     = 'omega_2',       &
        &              val     = me%omega_Cum(2), &
        &              default = -1._rk,          &
        &              ErrCode = iError           )

      if (me%omega_Cum(2) < 0._rk) then
        write(logUnit(1),*) 'Cumulant_extended is active'
        write(logUnit(1),*) '  Warning: Using omega_2 = omegaBulk'
      end if

      if (btest(iError, aoterr_fatal)) then
        call tem_abort('ERROR: wrong data type for omega_2 '// &
            &            'fluid table')
      end if

      do ii = 3, 10
        if (ii < 10) then
          write (omega_string, "(A6,I1)") omega_name, ii
        else
          write (omega_string, "(A6,I2)") omega_name, ii
        endif

        if (ii < 6) then
          omegaCum_def = -1._rk
        else
          omegaCum_def = 1._rk
        endif

        call aot_get_val(L       = conf,               &
          &              thandle = fluid_table,        &
          &              key     = trim(omega_string), &
          &              val     = me%omega_Cum(ii),   &
          &              default = omegaCum_def,       &
          &              ErrCode = iError              )

        if (btest(iError, aoterr_fatal)) then
          call tem_abort('ERROR: wrong data type for ' // trim(omega_string) &
              & // ' cumulant fluid table')
        end if

      enddo

      do ii = 1, 3
        write (omega_lim_string, "(A10,I1)") omega_lim_name, ii

        call aot_get_val(L       = conf,                       &
          &              thandle = fluid_table,                &
          &              key     = trim(omega_lim_string),     &
          &              val     = me%omega_Lim(ii),           &
          &              default = 1.e-2_rk,                   &
          &              ErrCode = iError                      )

        if (btest(iError, aoterr_fatal)) then
          call tem_abort('ERROR: wrong data type for ' &
            &            // trim(omega_lim_string) &
            &            // ' cumulant fluid table')
        end if

      enddo

    end if

    ! load nonNewtonian fluid feature
    call mus_nNwtn_load( me     = me%nNwtn,   &
      &                  conf   = conf,       &
      &                  parent = fluid_table )

    ! load turbulence
    call mus_load_turbulence( me     = me%turbulence, &
      &                       conf   = conf,          &
      &                       parent = fluid_table    )

    call aot_table_close( L=conf, thandle=fluid_table )
    call tem_horizontalSpacer(fUnit = logUnit(1))

  end subroutine mus_load_fluid
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> This routine initilizes fluid visocity and relaxation paramters for each
  !! level
  subroutine mus_init_fluid(me, physics, schemeHeader, minLevel, maxLevel,    &
    &                       levelDesc, pdf, stencil, general, tNow)
    ! --------------------------------------------------------------------------
    !> fluid type
    type(mus_fluid_type), intent(inout) :: me
    !> physics type to convert physics to lattice unit or vice versa
    type(mus_physics_type), intent(in) :: physics
    !> scheme header
    type(mus_scheme_header_type), intent(in) :: schemeHeader
    !> min and max level
    integer, intent(in) :: minLevel, maxLevel
    !> level descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc(minLevel:maxLevel)
    !> pdf info with neigh array for all levels
    type(pdf_data_type), intent(in) :: pdf(minLevel:maxLevel)
    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> general type contains communication pattern and proc info
    type(tem_general_type), intent(in)      :: general
    !> current simulation time
    type(tem_time_type),intent(in) :: tNow
    ! --------------------------------------------------------------------------
    integer :: iLevel
    ! --------------------------------------------------------------------------
    write(logUnit(1),'(A)') 'Initialize fluid'
    ! allocate array to store kinematic viscosity and bulk viscosity
    allocate(me%viscKine%dataOnLvl(minLevel:maxLevel))
    do iLevel = minLevel, maxLevel
      call init(me%viscKine%dataOnLvl(iLevel), pdf(iLevel)%nElems_local)
      me%viscKine%dataOnLvl(iLevel)%val = -100.0_rk
    end do
    allocate(me%viscBulkLvl(minLevel:maxLevel))
    allocate(me%omegaBulkLvl(minLevel:maxLevel))

    ! New initialization
    ! allocate relaxation parameter array to size nSolve
    call mus_init_relaxParam(omLvl    = me%viscKine%omLvl,  &
      &                      minLevel = minLevel,           &
      &                      maxLevel = maxLevel,           &
      &                      nElems   = pdf(:)%nElems_local )

    ! assign function pointer
    call mus_assign_mrt_ptr(me%mrtPtr, schemeHeader)

    ! intialize kinematic omega
    write(logUnit(3),'(A)') '  Set relaxation parameter from kinematic ' &
      &                     //'viscosity'
    do iLevel = minLevel, maxLevel
      call mus_update_relaxParamFromViscSTfun(             &
        & omega       = me%viscKine%omLvl(iLevel)%val,     &
        & visc        = me%viscKine%dataOnLvl(iLevel)%val, &
        & viscSTfun   = me%viscKine%STfun,                 &
        & viscRef     = physics%fac(iLevel)%visc,          &
        & nElems      = pdf(iLevel)%nElems_local,          &
        & baryOfTotal = levelDesc(iLevel)%baryOfTotal,     &
        & tNow        = tNow                               )

      ! Initialize bulk omega
      me%viscBulkLvl(iLevel) = me%viscBulk_phy / physics%fac(iLevel)%visc
      ! for 2D the formula is different!!!!!!
      !write(logUnit(3),'(A,I0)') '   nDims: ', stencil%nDims
      if (stencil%nDims == 2) then
        ! for 2d the formula of omega from bulk is different!
        ! see ref: "Theory of the Lattice Boltzmann Method: Dispersion,
        ! Dissipation, Isotropy, Galilean Invariance, and Stability"
        ! Luo 2000, PHYSICAL REVIEW E
        me%omegaBulkLvl(iLevel) = 1.0_rk / (24.0_rk * me%viscBulkLvl(iLevel) &
          & / (8._rk - 12._rk * cs2) + 0.5_rk )
      else
        ! see ref: "A D3Q27 multiple-relaxation-time lattice Boltzmann method
        ! for turbulent flows", Suga 2015, Computers&Fluids
        !me%omegaBulkLvl(iLevel) = 2.0_rk                                    &
        !  &                     / (9.0_rk * me%viscBulkLvl(iLevel) + 1.0_rk )
        me%omegaBulkLvl(iLevel) = 1.0_rk / (9.0_rk * me%viscBulkLvl(iLevel) &
          & / (5._rk - 9._rk * cs2) + 0.5_rk )
      endif

    end do

    ! Assign function pointer for nonNewtonian model
    if (me%nNwtn%active) then
      write(logUnit(3),'(A)') '  Assign function to compute viscosity ' &
        &                     //'for non-Newtonian model'
      call mus_assign_nNwtnVisc_ptr(me%nNwtn, schemeHeader)
    end if

    if (me%turbulence%active) then
      write(logUnit(3),'(A)') '  Assign function to compute viscosity ' &
        &                     //'for turbulence model'
      ! assign function pointer to compute turbulence viscosity
      call mus_assign_turbVisc_ptr(me%turbulence, schemeHeader)

      ! Initialize communication buffer for viscosity to compute
      ! velocity gradient for turbulence model
      ! This step must be done after construct_connectivity and
      ! init_levelBuffers
      allocate( me%turbulence%dataOnLvl( minLevel:maxLevel ) )
      do iLevel = minLevel, maxLevel
        call mus_init_turbulenceData(                     &
          & me         = me%turbulence%dataOnLvl(iLevel), &
          !& turbConfig = me%turbulence%config,            &
          & levelDesc  = levelDesc(iLevel),               &
          & pattern    = general%commPattern,             &
          & nSize      = pdf(iLevel)%nSize                )
      end do !iLevel

    end if

    ! Dump fluid information
    call mus_fluid_dump( me, minLevel, maxLevel, physics, pdf(:)%nElems_solve, &
      &                  general, logUnit(1), schemeHeader )

  end subroutine mus_init_fluid
  ! ************************************************************************** !

  ! ************************************************************************** !
  subroutine mus_fluid_dump( me, minLevel, maxLevel, physics, nSolve, general, &
    &                        outUnit, schemeHeader )
    ! --------------------------------------------------------------------------
    type( mus_fluid_type ), intent(inout) :: me !< fluid type
    !> minlevel and maxlevel
    integer, intent(in) :: minLevel, maxLevel
    !> physics type to convert physics to lattice unit or vice versa
    type( mus_physics_type ), intent(in) :: physics
    !> number of elements to solve per level (fluid+ghost)
    integer, intent(in) :: nSolve(minLevel:maxLevel)
    !> general type
    type(tem_general_type), intent(in) :: general
    integer, intent(in) :: outUnit
    !> scheme header
    type(mus_scheme_header_type), intent(in) :: schemeHeader
    ! --------------------------------------------------------------------------
    integer :: iLevel, iError
    real(kind=rk) :: visc_min, visc_max, viscRef, om_min, om_max
    real(kind=rk) :: glob_visc_min, glob_visc_max
    ! --------------------------------------------------------------------------

    write(outUnit,"(A)")       'Fluid properties:'
    write(outUnit,"(A,F10.7)") '  Reference density rho0: ', rho0
    write(outUnit,"(A)") '  Kinematic viscosity STfun kind: ' // &
      &                        trim(me%viscKine%STfun%fun_kind)

    do iLevel = minLevel, maxLevel
      ! get global min and max viscosity
      visc_min = minval(me%viscKine%dataOnLvl(iLevel)%val(1:nSolve(iLevel)))
      visc_max = maxval(me%viscKine%dataOnLvl(iLevel)%val(1:nSolve(iLevel)))

      call mpi_reduce( visc_min, glob_visc_min, 1, mpi_double_precision, &
                       mpi_min, 0, general%proc%comm, ierror    )
      call mpi_reduce( visc_max, glob_visc_max, 1, mpi_double_precision, &
                       mpi_max, 0, general%proc%comm, ierror    )
      viscRef = physics%fac(iLevel)%visc
      om_min = mus_calcOmegaFromVisc(glob_visc_min)
      om_max = mus_calcOmegaFromVisc(glob_visc_max)
      write(outUnit,'(A,I0)') 'Kinematic viscosity, (min, max) on level: ', &
        &                        iLevel
      write(outUnit,'(A,F10.7,A,F10.7,A)') '  Physical: (', &
        & glob_visc_min*viscRef, ',', glob_visc_max*viscRef, ')'
      write(outUnit,'(A,F10.7,A,F10.7,A)') '  Lattice: (', glob_visc_min, &
        &                                     ',', glob_visc_max, ')'
      write(outUnit,'(A,F10.7,A,F10.7,A)') '  Kinematic omega (min,max):(', &
        &                                     om_min,',', om_max, ')'

      write(outUnit,"(A)") ' viscBulk lattice and omega: '
      write(outUnit,"(A,F10.7)") '  viscBulk:  ', me%viscBulkLvl(iLevel)
      write(outUnit,"(A,F10.7)") '  omegaBulk: ', me%omegaBulkLvl(iLevel)

    end do

    write(outUnit, "(A,F10.7)") 'Magic lambda for TRT: ', me%lambda

    write(outUnit, "(A,F10.7)") 'HRR sigma: ', me%HRR_sigma

    write(outUnit, "(A,F10.7)") 'DRT tauN: ', me%DRT_tauN

    if (trim(schemeHeader%relaxation) == 'cumulant_extended' .or. &
      & trim(schemeHeader%relaxation) == 'cumulant_extended_generic') then
      write(outUnit, "(A)") 'omega Cumulant: -1.000 means adjust during runtime'
      do iLevel = 1, 10
        write(outUnit, "(A,I2,A,F10.7)") '  ', iLevel,' = ', me%omega_Cum(iLevel)
      enddo
      write(outUnit, "(A)") 'omega limiter Cumulant: 1.e10 means unlimited'
      do iLevel = 1, 3
        write(outUnit, "(A,I2,A,F10.7)") '  ', iLevel,' = ', me%omega_Lim(iLevel)
      enddo
      
      ! check Cumulant_extended range of admissibility for omegas
      write(logUnit(1),*) 'Checking admissibility of omegas for parametrized Cumulant scheme'

      do iLevel = minLevel, maxLevel
        call cumulant_omega_check( omegaVisc = me%viscKine%omLvl(iLevel)%val(:), &
          &                        omegaBulk = me%omegaBulkLvl(iLevel),          &
          &                        omegaIn   = me%omega_Cum(:),                  &
          &                        nSolve    = nSolve(iLevel),                   &
          &                        level     = iLevel                            )
      enddo
    end if

    if( maxval( abs( me%force )) > 0._rk )  &
      & write(outUnit,*) ' Forcing:             ', real(me%force)

    ! Dump nonNewtonian parameters to outUnit
    call mus_nNwtn_dump2outUnit( me%nNwtn, outUnit )

  end subroutine mus_fluid_dump
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> write fluid prop into a lua file
  !!
  subroutine mus_fluid_save2lua( me, conf )
    ! ---------------------------------------------------------------------------
    !> single fluid type
    type( mus_fluid_type ), intent(in) :: me
    type( aot_out_type ) :: conf
    ! ---------------------------------------------------------------------------

    call aot_out_open_table( put_conf = conf, tname = 'fluid' )

    call aot_out_val( put_conf = conf,             &
      &               vname    = 'bulk_viscosity', &
      &               val      = me%viscBulk_phy   )
    if (trim(me%viscKine%STfun%fun_kind)=='const') then
      call aot_out_val( put_conf = conf,                      &
        &               vname    = 'kine_viscosity',          &
        &               val      = me%viscKine%STfun%const(1) )
    end if

    if ( me%nNwtn%active ) call mus_nNwtn_save2lua( me%nNwtn, conf )

    call aot_out_close_table( put_conf = conf )

  end subroutine mus_fluid_save2lua
  ! **************************************************************************** !

  ! ************************************************************************** !
  !> This routines act as a destructor for fluid type
  subroutine mus_fluid_cleanup( me )
    ! ---------------------------------------------------------------------------
    !> single fluid type
    type( mus_fluid_type ), intent(inout) :: me
    ! ---------------------------------------------------------------------------

    if (allocated(me%viscKine%dataOnLvl)) then
      deallocate(me%viscKine%dataOnLvl)
    end if
    if (allocated(me%viscKine%omLvl)) then
      deallocate(me%viscKine%omLvl)
    end if
    if (allocated(me%viscBulkLvl)) then
      deallocate(me%viscBulkLvl)
    end if
    if (allocated(me%omegaBulkLvl)) then
      deallocate(me%omegaBulkLvl)
    end if
    if (me%turbulence%active) then
      if (allocated(me%turbulence%dataOnLvl)) then
        deallocate(me%turbulence%dataOnLvl)
      end if
    end if

  end subroutine mus_fluid_cleanup
  ! ************************************************************************** !

end module mus_fluid_module
! ****************************************************************************** !
