! See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
!> author: Kannan Masilamani
!! author: Simon Zimny
!! This module contains information about all fields like fluid, 
!! species, temperature etc. This field type will be used for
!! multispecies and passive scalar transport.
!!
!! [mus_field_prop]: @ref mus_field_prop_module "mus_field_prop_module"
!!
module mus_field_module

  ! include treelm modules
  use env_module,               only: labelLen, rk
  use tem_param_module,         only: cs2inv, cs2, div1_2
  use tem_math_module,          only: invert_matrix
  use tem_time_module,          only: tem_time_type
  use tem_bc_prop_module,       only: tem_bc_prop_type
  use tem_debug_module,         only: dbgUnit
  use tem_tools_module,         only: tem_horizontalSpacer
  use tem_aux_module,           only: tem_abort, tem_CheckLabel
  use tem_ini_condition_module, only: tem_ini_condition_type, tem_load_ic
  use tem_varSys_module,        only: tem_varSys_type
  use tem_varMap_module,        only: tem_possible_variable_type
  use tem_restart_module,       only: tem_restart_type
  use tem_logging_module,       only: logUnit
  use tem_temporal_module,      only: tem_temporal_for
  use tem_construction_module,  only: tem_levelDesc_type
  use tem_stencil_module,       only: tem_stencilHeader_type
  use tem_stringKeyValuePair_module, only: init, truncate

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

  ! include musubi modules
  use mus_bc_header_module,     only: boundary_type, mus_load_bc,              &
    &                                 glob_boundary_type, &
    &                                 check_solid_in_bc, rearrange_bc_elems
  use mus_scheme_header_module, only: mus_scheme_header_type
  use mus_scheme_layout_module, only: mus_scheme_layout_type
  use mus_field_prop_module,    only: mus_field_prop_type, mus_field_prop_out, &
    &                                 mus_load_field_prop
  use mus_species_module,       only: compute_molWeightRatio,                  &
    &                                 compute_bulkViscOmega
  use mus_physics_module,       only: mus_physics_type
  use mus_mixture_module,       only: mus_mixture_type, mus_load_mixture
  use mus_source_var_module,    only: mus_source_type, mus_load_source_var
  use mus_nernstPlanck_module,  only: mus_nernstPlanck_type, mus_load_nernstPlanck

  implicit none

  private

  public :: mus_field_type
  public :: mus_load_fields, mus_fields_out
  public :: mus_load_fieldBaseInfos
  public :: setParameters_multispecies
  public :: remove_solid_in_bc
  public :: mus_check_allWall
  public :: mus_field_getSymmetricBCs


  !> This type contains all information on fields with ic and bc
  !! Example fields: fluid, species etc.
  !! Each field contains one initial condition and array of 
  !! boundary conditions
  !!
  type mus_field_type
    !> field label. Should be unique for each field
    character(len=labelLen) :: label
    !> lbm parameters (fluid and species) for field
    type( mus_field_prop_type ) :: fieldProp
    !> array of field boundary types for each field
    !! size: #BCs in the boundary_condition table
    !! allocated in mus_load_bc
    type( boundary_type ), allocatable :: bc(:)
    !> initialization case, one initial condition for each field
    type( tem_ini_condition_type ) :: ic
    !> field source applied only to current field
    type( mus_source_type ) :: source
    !> An instance of restart type
    type( tem_restart_type ) :: restart
  end type mus_field_type

  !> Interface for dumping a single field or a set of fields in a file in lua 
  !! format.
  !!
  interface mus_fields_out
    module procedure mus_fields_out_vec
    module procedure mus_field_out_scal
  end interface mus_fields_out

contains

! ****************************************************************************** !
  !> Subroutine to load the field table from the lua configuration file.
  !!
  !! If field table is not defined than load bc, ic, fluid, species from scheme
  !! table. If scheme table is not defined than load field variables from
  !! config parent.
  !!
  subroutine mus_load_fields( me, varSys, nFields, mixture, nernstPlanck,   &
    &                         bc_prop, conf, parent, minLevel, maxLevel,    &
    &                         schemeHeader, poss_srcVar, physics, rank,     &
    &                         scaling, layout, isMusHvs )
    ! --------------------------------------------------------------------------
    !> array of field type
    type(mus_field_type), intent(inout) :: me(:)
    !> Global variable system required to append annoymous source and 
    !! boundary variables
    type(tem_varSys_type), intent(inout) :: varSys
    integer, intent(in) :: nFields !< number of fields defined in lua file
    !> contains mixture information
    type(mus_mixture_type), intent(out) :: mixture
    !> contains solvent information
    type(mus_nernstPlanck_type), intent(out) :: nernstPlanck
    !> boundary data from mesh
    type(tem_bc_prop_type), intent(in):: bc_prop
    !> flu state
    type(flu_State), intent(inout) :: conf
    !> global pdf info
    integer, intent(in) :: minLevel, maxLevel
    !> identifier of the scheme
    type( mus_scheme_header_type ), intent(in) :: schemeHeader
    !> parent handle if scheme table is defined
    integer, intent(in), optional :: parent
    !> possible source variables
    type(tem_possible_variable_type), intent(in) :: poss_srcVar
    !> physics type to convert physics to lattice unit or vice versa
    type( mus_physics_type ), intent(in) :: physics
    !> the current rank
    integer, intent(in) :: rank
    !> scaling type
    character(len=labelLen), intent(in) :: scaling
    !> fluid stencil info
    type(mus_scheme_layout_type), intent(in) :: layout
    !> Logic to not to load tracking and variable table if this routine 
    !! is called from mus_hvs_config_load.
    !! Default is False
    logical, optional, intent(in) :: isMusHvs
    ! ---------------------------------------------------------------------------
    ! counter variables
    integer :: iField
    ! aotus handles
    integer :: field_handle, field_sub_handle
    ! ---------------------------------------------------------------------------
    call tem_horizontalSpacer(fUnit = logUnit(1))
    write(logUnit(1),*) 'Loading the fields ...' 
    write(logUnit(1),*) 'Number of fields to load: ', nFields 
    call aot_table_open( L       = conf,         &
      &                  parent  = parent,       &
      &                  thandle = field_handle, &
      &                  key     = 'field'       ) 

    if ( field_handle == 0 ) then
      if (present(parent)) then
        ! field table is not defined load field variable from scheme 
        write(logUnit(1),*) 'No field table defined.'          &
          &              // 'Loading field data from scheme' 

        call mus_load_field_single( me           = me(1),        &
          &                         varSys       = varSys,       &
          &                         nFields      = nFields,      &
          &                         bc_prop      = bc_prop,      &
          &                         conf         = conf,         &
          &                         minLevel     = minLevel,     &
          &                         maxLevel     = maxLevel,     &
          &                         parent       = parent,       &
          &                         poss_srcVar  = poss_srcVar,  &
          &                         physics      = physics,      &
          &                         schemeHeader = schemeHeader, &
          &                         scaling      = scaling,      &
          &                         rank         = rank,         &
          &                         layout       = layout,       &
          &                         isMusHvs     = isMusHvs      )

      ! scheme table is not defined load field variable from config
      else
        write(logUnit(1),*) 'No field table and scheme table are defined.' &
          &              // 'Loading field data from config' 

        call mus_load_field_single( me           = me(1),        &
          &                         varSys       = varSys,       &
          &                         nFields      = nFields,      &
          &                         bc_prop      = bc_prop,      &
          &                         conf         = conf,         &
          &                         minLevel     = minLevel,     &
          &                         maxLevel     = maxLevel,     &
          &                         poss_srcVar  = poss_srcVar,  &
          &                         physics      = physics,      &
          &                         schemeHeader = schemeHeader, &
          &                         scaling      = scaling,      &
          &                         rank         = rank,         &
          &                         layout       = layout,       &
          &                         isMusHvs     = isMusHvs      )
      endif ! present parent?
    else ! check whether field is a single table or multiple table

      call aot_table_open( L       = conf,                                     &
        &                  parent  = field_handle,                             & 
        &                  thandle = field_sub_handle,                         &
        &                  pos     = 1 )
      if ( field_sub_handle == 0 ) then
        ! field is a single table
        call aot_table_close( L=conf, thandle=field_sub_handle )
        write(logUnit(1),*) 'Field is a single table' 

        call mus_load_field_single( me           = me(1),        &
          &                         varSys       = varSys,       &
          &                         nFields      = nFields,      &
          &                         bc_prop      = bc_prop,      &
          &                         conf         = conf,         &
          &                         parent       = field_handle, &
          &                         minLevel     = minLevel,     &
          &                         maxLevel     = maxLevel,     &
          &                         poss_srcVar  = poss_srcVar,  &
          &                         physics      = physics,      &
          &                         schemeHeader = schemeHeader, &
          &                         scaling      = scaling,      &
          &                         rank         = rank,         &
          &                         layout       = layout,       &
          &                         isMusHvs     = isMusHvs      )

      else ! field is multiple table
        call aot_table_close( L=conf, thandle=field_sub_handle )
        write(logUnit(1),*) 'Field is a multiple table'
        do iField = 1, nFields
          call aot_table_open( L       = conf,                                 &
            &                  parent  = field_handle,                         &
            &                  thandle = field_sub_handle,                     &
            &                  pos     = iField )

          write(logUnit(1),*) 'Properties for Field ', iField

          call tem_horizontalSpacer(fUnit = logUnit(1))
          call mus_load_field_single( me           = me(iField),       &
            &                         varSys       = varSys,           &
            &                         nFields      = nFields,          &
            &                         bc_prop      = bc_prop,          &
            &                         conf         = conf,             &
            &                         minLevel     = minLevel,         &
            &                         maxLevel     = maxLevel,         &
            &                         parent       = field_sub_handle, &
            &                         poss_srcVar  = poss_srcVar,      &
            &                         physics      = physics,          &
            &                         schemeHeader = schemeHeader,     &
            &                         scaling      = scaling,          &
            &                         rank         = rank,             &
            &                         layout       = layout,           &
            &                         isMusHvs     = isMusHvs          )

          call aot_table_close( L=conf, thandle=field_sub_handle )
        end do ! number of fields
      end if ! if field table is present
    end if

    call aot_table_close( L=conf, thandle=field_handle )

    !load scheme specific field
    select case(trim(schemeHeader%kind))
    case('lbm', 'lbm_nNwtn', 'lbm_incomp', 'lbm_incomp_nNwtn', 'passive_scalar', &
      & 'isotherm_acEq', 'poisson','poisson_boltzmann_linear', &
      & 'poisson_boltzmann_nonlinear')
      !nFields must be = 1 scheme other than multi-species and nernst_planck
      if( nFields > 1 ) then
        write(logUnit(1),*) ' ERROR: Number of fields defined for '
        write(logUnit(1),*) '        '//trim(schemeHeader%kind)//' scheme: >1. '
        write(logUnit(1),*) '        Specify only one field. Aborting'
        call tem_abort()
      endif
    case('multi-species_gas','multi-species_liquid')
      !nFields must be > 1 for multispecies
      if( nFields == 1 ) then
        write(logUnit(1),*) ' ERROR: Number of fields(species) defined for '
        write(logUnit(1),*) '        multi-species scheme = 1. '
        write(logUnit(1),*) '        Specify nFields > 1. Aborting'
        call tem_abort()
      endif

      !compute molecular weight ratios incase of multispecies simulation
      call compute_molWeightRatio(                              &
        &  molWeights    = me(:)%fieldProp%species%molWeight,   &
        &  molWeigRatios = me(:)%fieldProp%species%molWeigRatio )

      !load mixture table and set bulk viscosity relaxation time
      ! for each species which is depends on molecular weight ratios
      call mus_load_mixture( me           = mixture,      &
        &                    conf         = conf,         &
        &                    parent       = parent,       &
        &                    schemeHeader = schemeHeader, & 
        &                    minLevel     = minLevel,     &
        &                    maxLevel     = maxLevel,     &
        &                    physics      = physics,      &
        &                    nFields      = nFields       )

      !compute bulk viscosity omega for each field
      write(logUnit(1),*)'  Bulk omega for each species:'
      do iField = 1,nFields
        write(logUnit(1),*) '   species ', iField
        call compute_bulkViscOmega( species = me(iField)%fieldProp%species,    &
          &                         bulkvisc = mixture%bulk_viscosityLB,       &
          &                         bulkviscLvl = mixture%relaxLvl(:)%bulkVisc,&
          &                         minLevel = minLevel, &
          &                         maxLevel = maxLevel )
      end do
    case('nernst_planck')
      write(logUnit(1),"(A)") ' Loading properties for nersnt_planck:'
      call mus_load_nernstPlanck(  me         = nernstPlanck, &
        &                          conf       = conf,         &
        &                          parent     = parent,       &
        &                          minLevel   = minLevel,     &
        &                          physics    = physics       )
    case default
      write(logUnit(1),*) 'The selected scheme kind is unknown '//           &
        &             trim(schemeHeader%kind)
      call tem_abort()
    end select

  end subroutine mus_load_fields
! ****************************************************************************** !


! ****************************************************************************** !
  !> load a single field table
  !! In includes:
  !!   load field property
  !!   load source variables
  !!   load boundary defination
  !!   load immersed boundary method
  !!   load initial condition defination and its property
  !!
  subroutine mus_load_field_single( me, varSys, nFields, bc_prop, conf,        &
    &                               parent, minLevel, maxLevel, poss_srcVar,   &
    &                               physics, schemeHeader, rank, scaling,      &
    &                               layout, isMusHvs )
    ! --------------------------------------------------------------------------
    !> field type
    type(mus_field_type), intent(inout) :: me
    !> Global variable system required to append annoymous source and 
    !! boundary variables
    type(tem_varSys_type), intent(inout) :: varSys
    !> number of fields defined in lua file
    integer, intent(in) :: nFields 
    !> boundary data from mesh
    type(tem_bc_prop_type), intent(in):: bc_prop
    !> flu state
    type(flu_State), intent(inout) :: conf
    !> parent handle if scheme table is defined
    integer, intent(in), optional :: parent
    !> global pdf info
    integer, intent(in) :: minLevel, maxLevel
    !> possible source variables
    type(tem_possible_variable_type), intent(in) :: poss_srcVar
    !> 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
    !> the current rank
    integer, intent(in) :: rank
    !> scaling type
    character(len=labelLen), intent(in) :: scaling
    !> fluid stencil info
    type(mus_scheme_layout_type), intent(in) :: layout
    !> Logic to not to load tracking and variable table if this routine 
    !! is called from mus_hvs_config_load.
    !! Default is False
    logical, optional, intent(in) :: isMusHvs
    ! --------------------------------------------------------------------------
    character(len=labelLen), allocatable :: ic_states(:)
    logical :: isMusHvs_loc
    integer :: IC_nVars
    !> Error code of loading ic variables
    integer, allocatable :: iError(:)
    ! --------------------------------------------------------------------------
    if (present(isMusHvs)) then
      isMusHvs_loc = isMusHvs
    else
      isMusHvs_loc = .false.
    end if

    write(logUnit(1),"(A)") ' Loading field: '//trim(me%label)
    ! load field properties
    call mus_load_field_prop( me           = me%fieldProp, &
      &                       nFields      = nFields,      &
      &                       conf         = conf,         &
      &                       minLevel     = minLevel,     &
      &                       maxLevel     = maxLevel,     &
      &                       parent       = parent,       &
      &                       physics      = physics,      &
      &                       scaling      = scaling,      &
      &                       cs_lattice   = layout%cs,    &
      &                       schemeHeader = schemeHeader  )

    ! Load Boundary, initial condition only for solver
    if (.not. isMusHvs_loc) then

      ! Read initial condition
      call mus_set_ic_states( schemeHeader%kind, ic_states, IC_nVars )
      allocate( iError( size(ic_states) ) )

      call tem_load_ic( me        = me%ic,               &
        &               conf      = conf,                &
        &               parent    = parent,              &
        &               key       = 'initial_condition', &
        &               ErrCode   = iError,              &
        &               StateName = ic_states            )

      if ( any( iError(1:IC_nVars) /= 0 ) ) then
        write(logUnit(1), "(A)") " Not all the KEY variables in initial condition" &
          &                      //"are well defined!"
        call tem_abort()
      end if

      deallocate( ic_states )
      deallocate( iError    )

      ! load boundary
      call mus_load_bc( me      = me%bc,          &
        &               bc_prop = bc_prop,        &
        &               conf    = conf,           &
        &               parent  = parent,         &
        &               varSys  = varSys,         &
        &               stencil = layout%fStencil )

      ! load the field specific source variables
      call mus_load_source_var( me       = me%source,   &
        &                       possVars = poss_srcVar, &
        &                       conf     = conf,        &
        &                       parent   = parent,      &
        &                       key      = 'source',    &
        &                       varSys   = varSys       )

    else
      ! allocate arrays for mus_harvesting
      allocate(me%bc(0))
    end if

    call tem_horizontalSpacer(fUnit = logUnit(1))

  end subroutine mus_load_field_single
! ****************************************************************************** !


  ! ************************************************************************ !
  !> This routine returns nFields and field labels from config file.
  !! It is required to initialize variable system.
  !! labels are loaded only if field table is present else default
  !! is set to empty string.
  subroutine mus_load_fieldBaseInfos( me, nFields, parent, conf )
    ! ---------------------------------------------------------------------------
    !> array of field type
    type( mus_field_type ), allocatable, intent(out) :: me(:)
    !> number of fields defined in lua file
    integer, intent(out) :: nFields
    !> parent handle if scheme table is defined
    integer, intent(in), optional :: parent
    !> flu state
    type(flu_State), intent(inout) :: conf
    ! ---------------------------------------------------------------------------
    integer :: iError, iField
    integer :: field_handle, field_sub_handle
    ! ---------------------------------------------------------------------------
    call tem_horizontalSpacer(fUnit = logUnit(1))
    write(logUnit(1),*) 'Loading the field base info ...' 

    call aot_table_open( L       = conf,         &
      &                  parent  = parent,       &
      &                  thandle = field_handle, &
      &                  key     = 'field'       ) 

    ! allocate the fields
    nFields = 1
    allocate(me(nFields))
    me(nFields)%label = ''

    ! load label only if field table is present 
    if ( field_handle /= 0 ) then
      ! check whether field is a single table or multiple table
      call aot_table_open( L       = conf,                                     &
        &                  parent  = field_handle,                             & 
        &                  thandle = field_sub_handle,                         &
        &                  pos     = 1 )

      if ( field_sub_handle == 0 ) then
        ! field is a single table
        call aot_get_val( L       = conf,         &
           &              thandle = field_handle, &
           &              key     = 'label',      &
           &              val     = me(1)%label,  &
           &              default = '',           &
           &              ErrCode = iError        )

        if( trim(me(1)%label) /= '')                                           &
          & me(1)%label = trim( me(1)%label )//'_'
        write(logUnit(1),*) 'Field label: ', trim(me(1)%label)

      else ! field is multiple table
        call aot_table_close( L=conf, thandle=field_sub_handle )
        nFields = aot_table_length( L=conf, thandle=field_handle )
        write(logUnit(1),*) 'Number of fields: ', nFields 

        ! reallocate the fields
        deallocate( me )
        allocate( me( nFields ))

        do iField = 1, nFields
          call aot_table_open( L       = conf,                                 &
            &                  parent  = field_handle,                         & 
            &                  thandle = field_sub_handle,                     &
            &                  pos     = iField )

          call aot_get_val( L       = conf,             &
            &               thandle = field_sub_handle, &
            &               key     = 'label',          &
            &               val     = me(iField)%label, &
            &               default = '',               &
            &               ErrCode = iError            )

          if (btest(iError, aoterr_NonExistent)) then
            write(logUnit(1),*) 'Error: field label is not specified' 
            write(logUnit(1),*) 'field label is neccessary when nFields>1' 
            call tem_abort()
          end if

          ! add underscore here to differentiate state variables and derived 
          ! variables of each field
          if( trim(me(iField)%label) /= '' )                                   &
            &          me(iField)%label = trim( me(iField)%label )//'_'

          ! Check, if the names of the fields are unique. If not, add an
          ! increasing number
          call tem_checkLabel( label = me(:)%label, nLabels = iField )

          call aot_table_close( L=conf, thandle=field_sub_handle )
          write(logUnit(1),*) iField, 'Field label: ', trim(me(iField)%label)
        end do !number of fields

      end if ! single field  
    else
      write(logUnit(1),*) 'No field table defined.'
      write(logUnit(1),*) 'Assuming single field and label = "".'
    end if
    call tem_horizontalSpacer(fUnit = logUnit(1))

  end subroutine mus_load_fieldBaseInfos
  ! ************************************************************************ !


! ****************************************************************************** !
  !> write array of fields into a lua file
  !!
  subroutine mus_fields_out_vec( me, conf, schemeHeader )
    ! ---------------------------------------------------------------------------
    !> array of field type
    type( mus_field_type ), intent(in) :: me(:)
    !> aotus out type
    type( aot_out_type ), intent(inout) :: conf
    !> identifier of the scheme
    type( mus_scheme_header_type ), intent(in) :: schemeHeader
    ! ---------------------------------------------------------------------------
    integer :: iField ! counter variable
    ! ---------------------------------------------------------------------------

    call aot_out_val( put_conf = conf,                                         &
      &               val      = size(me),                                     &
      &               vname    = 'nFields' )

    if (size(me) > 1) then
      call aot_out_open_table(put_conf = conf, tname = 'field')
      do iField = 1, size(me)
        call mus_field_out_scal( me           = me( iField ), &
          &                      conf         = conf,         &
          &                      schemeHeader = schemeHeader, & 
          &                      level        = 1             )
      end do
      call aot_out_close_table(put_conf = conf)
    else
      call mus_field_out_scal( me           = me( 1 ),      &
        &                      conf         = conf,         &
        &                      schemeHeader = schemeHeader, & 
        &                      level        = 0             )
    end if    
  end subroutine mus_fields_out_vec
! ****************************************************************************** !


! ****************************************************************************** !
  !> write single field into a lua file
  !!
  subroutine mus_field_out_scal( me, conf, schemeHeader, level )
    ! ---------------------------------------------------------------------------
    !> single field type
    type( mus_field_type ), intent(in) :: me
    !> aotus out type
    type( aot_out_type ), intent(inout) :: conf
    !> identifier of the scheme
    type( mus_scheme_header_type ), intent(in) :: schemeHeader
    !> To dump field with or without key
    integer, optional, intent(in) :: level
    ! ---------------------------------------------------------------------------
    integer :: level_loc, pos
    ! ---------------------------------------------------------------------------
    if (present(level)) then
      level_loc = level
    else
      level_loc = 0
    end if

    if( level_loc == 0) then
      call aot_out_open_table( put_conf = conf, tname = 'field' )
    else
      call aot_out_open_table( put_conf = conf )
    end if

    if (len(trim(me%label)) > 0) then
      ! Remove underscore from field name before dumping
      pos = INDEX(trim(me%label),'_', back=.true.)
      call aot_out_val( put_conf = conf,                   &
        &               vname    = 'label',                &
        &               val      = trim(me%label(1:pos-1)) )
    end if

    call mus_field_prop_out( me           = me%fieldProp,                      &
      &                      conf         = conf,                              &
      &                      schemeHeader = schemeHeader )
    call aot_out_close_table(put_conf = conf)

  end subroutine mus_field_out_scal
! ****************************************************************************** !


! ****************************************************************************** !
  !> Set ic states labels by scheme kind
  !!
  subroutine mus_set_ic_states( scheme_kind, ic_states, IC_nVars )
    ! --------------------------------------------------------------------------
    character(len=labelLen),  intent(in) :: scheme_kind
    character(len=labelLen), allocatable :: ic_states(:)
    !> Number of initial condition variables required to initialize state
    integer, intent(out) :: IC_nVars 
    ! --------------------------------------------------------------------------

    ! load scheme specific field
    select case( trim(scheme_kind) )
    case('lbm', 'lbm_incomp')
      IC_nVars = 4
      allocate(ic_states(10))
      ic_states = (/ 'pressure ', 'velocityX', 'velocityY', 'velocityZ',&
        &            'Sxx      ', 'Syy      ', 'Szz      ', 'Sxy      ',&
        &            'Syz      ', 'Sxz      ' /)
    case('lbm_nNwtn', 'lbm_incomp_nNwtn')
      IC_nVars = 4
      allocate(ic_states(11))
      ic_states = (/ 'pressure ', 'velocityX', 'velocityY', 'velocityZ',&
        &            'Sxx      ', 'Syy      ', 'Szz      ', 'Sxy      ',&
        &            'Syz      ', 'Sxz      ', 'omega    '/)
    case('passive_scalar')
      IC_nVars = 4
      allocate(ic_states(4))
      ic_states = (/ 'pressure ', 'velocityX', 'velocityY', 'velocityZ' /)
    case('multi-species_gas')
      IC_nVars = 4
      allocate( ic_states(4) )
      ic_states = (/ 'pressure ', 'velocityX', 'velocityY',&
        &            'velocityZ' /)
    case('multi-species_liquid', 'nernst_planck')
      IC_nVars = 4
      allocate( ic_states(4) )
      ic_states = (/ 'mole_fraction', 'velocityX    ', 'velocityY    ', &
        &            'velocityZ    ' /)
    case('isotherm_acEq')
      IC_nVars = 4
      allocate(ic_states(10))
      ic_states = (/ 'pressure ', 'velocityX', 'velocityY', 'velocityZ',&
        &            'Sxx      ', 'Syy      ', 'Szz      ', 'Sxy      ',&
        &            'Syz      ', 'Sxz      ' /)
    case('poisson', 'poisson_boltzmann_linear', 'poisson_boltzmann_nonlinear')
      IC_nVars = 1
      allocate(ic_states(1))
      ic_states = (/ 'potential' /)
    case default
      write(logUnit(1),*) 'The selected scheme kind model is not '&
         &              //'supported to initialize ic state names: '&
         &                 //trim( scheme_kind )
      call tem_abort()
    end select

  end subroutine mus_set_ic_states
! ****************************************************************************** !

! ****************************************************************************** !
  !> Set parameters for multispecies
  subroutine setParameters_multispecies( field, nFields, mixture, header, layout, &
    &                                    iLevel, tNow )
    ! --------------------------------------------------------------------------
    integer,                        intent(in) :: nFields
    type( mus_field_type ),      intent(inout) :: field(nFields)
    type( mus_scheme_layout_type ), intent(in) :: layout
    type( mus_scheme_header_type ), intent(in) :: header
    type( mus_mixture_type ),    intent(inout) :: mixture
    integer,                        intent(in) :: iLevel
    !> solver general info
    type( tem_time_type ),       intent(in) :: tNow
    ! --------------------------------------------------------------------------
    real(kind=rk) :: omega_diff ! diffusive relaxation parameter
    real(kind=rk) :: omega_kine ! kinematic shear relaxation parameter
    real(kind=rk) :: omega_bulk ! bulk relaxation parameter
    real(kind=rk) :: fac ! ramping multiplication factor
    integer :: iField, iDir
    real(kind=rk), dimension( layout%fStencil%QQ, layout%fStencil%QQ ) :: &
      &                                       identity, tmpMatrix, tmpMatrixInv
    ! --------------------------------------------------------------------------
!write(dbgUnit,*) 'setting omega paramter for multi-species'
    ! get the temporal multiplication factor
    fac = tem_temporal_for( temporal = mixture%omega_ramping, &
      &                     time     = tNow )
    !! Relaxation parameter for each level
    !! kine_viscosity = dt*c^2*(1/omega - 0.5)
    !!    =>   omega = dt / (viscosity/c^2 + dt/2)
    !! omega for each level is stored at fluid%omLvl
    omega_diff  =  fac * mixture%relaxLvl( iLevel )%paramB * cs2
!write(dbgUnit,*) 'omega_diff ', omega_diff
    mixture%relaxLvl( iLevel )%omega_diff = omega_diff


    ! if relaxation is mrt set kinematic and bulk viscosity omega
    if( header%relaxation(1:3) == 'mrt' ) then
      do iField = 1, nFields
        if( .not. allocated( field( iField )%fieldProp%species          &
          &          %mrt( iLevel )%s_mrt))                             &
          & allocate( field( iField )%fieldProp%species                 &
          &          %mrt( iLevel )%s_mrt( layout%fStencil%QQ,          &
          &                                layout%fStencil%QQ))

        if( .not. allocated( field( iField )%fieldProp%species            &
          &                                 %mrt( iLevel )%omegaMoments)) &
          & allocate( field( iField )%fieldProp%species                   &
          &          %mrt( iLevel )%omegaMoments( layout%fStencil%QQ,     &
          &                                       layout%fStencil%QQ)     )

        if( .not. allocated( field( iField )%fieldProp%species             &
          &                                 %mrt( iLevel )%omegaMomForce)) &
          & allocate( field( iField )%fieldProp%species                    &
          &          %mrt( iLevel )%omegaMomForce( layout%fStencil%QQ,     &
          &                                        layout%fStencil%QQ)     )

        field( iField )%fieldProp%species%mrt( iLevel )       &
          &                              %omegaMoments = 0.0_rk
        field( iField )%fieldProp%species%mrt( iLevel )       &
          &                              %omegaMomForce = 0.0_rk

        tmpMatrix    = 0.0_rk
        tmpMatrixInv = 0.0_rk
        Identity     = 0.0_rk
        do iDir = 1, layout%fStencil%QQ
          Identity(iDir, iDir) = 1.0_rk
        end do

        ! viscous omega
        omega_kine = mixture%relaxLvl( iLevel )%omega_kine * fac
        ! bulk omega
        omega_bulk = field(iField)%fieldProp%species%omBulkLvl(iLevel)  &
          &        * fac
!write(*,*) 'omega_diff ', omega_diff 
!write(*,*) 'omega_kine ', omega_kine        
!write(*,*) 'omega_bulk ', omega_bulk    
        ! Set the relaxation parameters (only the non-relaxing modes)
        ! \ref Multiple-relaxation-time lattice Boltzmann scheme for 
        !      homogeneous mixture flows with external force- Pietro Asinari
        ! \ref Jens ICCMES paper
        ! initialize all entries in relaxation matrix
        field( iField )%fieldProp%species%mrt( iLevel )%s_mrt(:,:) = 0.0_rk
        select case(trim(header%layout))
        case('d2q9' )
!write(*,*) 'omega_diaCenter ', (omega_bulk + omega_kine )/2.0_rk
!write(*,*) 'omega_nondia ', (omega_bulk - omega_kine )/2.0_rk
          ! diffusivity omegas
          do iDir=2,3
            field( iField )%fieldProp%species%mrt( iLevel )             &
              & %s_mrt( iDir, iDir ) = omega_diff
          end do    

          do iDir=4,5
            field( iField )%fieldProp%species%mrt( iLevel )             &
              & %s_mrt( iDir, iDir ) = ( omega_bulk + omega_kine ) * div1_2
          end do

          field( iField )%fieldProp%species%mrt( iLevel )               &
              & %s_mrt( 5, 4 ) = ( omega_bulk - omega_kine ) * div1_2

          field( iField )%fieldProp%species%mrt( iLevel )               &
              & %s_mrt( 4, 5 ) = ( omega_bulk - omega_kine ) * div1_2

          field( iField )%fieldProp%species%mrt( iLevel )               &
            & %s_mrt( 6, 6 ) = omega_kine

          do iDir=7,9
            field( iField )%fieldProp%species%mrt( iLevel )             &
              & %s_mrt( iDir, iDir ) = mixture%omega_hom
          end do

        case('d3q19')
!write(*,*) 'omega_diaCenter ', (omega_bulk + 2.0_rk * omega_kine )/3.0_rk
!write(*,*) 'omega_nondia ', (omega_bulk - omega_kine )/3.0_rk
          do iDir=2,4
            field( iField )%fieldProp%species%mrt( iLevel )             &
              & %s_mrt( iDir, iDir ) = omega_diff
          end do

          do iDir=5,7
            field( iField )%fieldProp%species%mrt( iLevel )             &
              & %s_mrt( iDir, iDir ) = ( omega_bulk + 2.0_rk * omega_kine ) / 3.0_rk
          end do

          field( iField )%fieldProp%species%mrt( iLevel )%s_mrt( 5, 6 ) &
            & = ( omega_bulk - omega_kine ) / 3.0_rk
          field( iField )%fieldProp%species%mrt( iLevel )%s_mrt( 5, 7 ) &
            & = ( omega_bulk - omega_kine ) / 3.0_rk

          field( iField )%fieldProp%species%mrt( iLevel )%s_mrt( 6, 5 ) &
            & = ( omega_bulk - omega_kine ) / 3.0_rk
          field( iField )%fieldProp%species%mrt( iLevel )%s_mrt( 6, 7 ) &
            & = ( omega_bulk - omega_kine ) / 3.0_rk

          field( iField )%fieldProp%species%mrt( iLevel )%s_mrt( 7, 5 ) &
            & = ( omega_bulk - omega_kine ) / 3.0_rk
          field( iField )%fieldProp%species%mrt( iLevel )%s_mrt( 7, 6 ) &
            & = ( omega_bulk - omega_kine ) / 3.0_rk

          do iDir=8,10 
            field( iField )%fieldProp%species%mrt( iLevel )             &
              & %s_mrt( iDir, iDir ) = omega_kine
          end do

          do iDir=11,19
            field( iField )%fieldProp%species%mrt( iLevel )             &
              & %s_mrt( iDir, iDir ) = mixture%omega_hom
          end do
        end select

        !omega moments = (M^-1 RelaxMat M)*(I+(M^-1 RelMat M)/2)^-1
        tmpMatrix = matmul( matmul( layout%moment%toPdf%A,              &
          & field( iField )%fieldProp%species%mrt( iLevel )%s_mrt ),    &
          & layout%moment%toMoments%A )

        tmpMatrixInv = invert_matrix( Identity + tmpMatrix * 0.5_rk )

        field( iField )%fieldProp%species%mrt( iLevel )%omegaMoments =  &
          & matmul(tmpMatrix, tmpMatrixInv)
          
        field( iField )%fieldProp%species%mrt( iLevel )      &
          &                      %omegaMomForce = tmpMatrixInv

      end do ! iField
    end if ! mrt?
!stop
  end subroutine setParameters_multispecies
! ****************************************************************************** !

! ****************************************************************************** !
!> First check count number of valid elements (non-solid) in each BC.
!! Then rearrange BC elements list so it contains only valid elements.
!! Update fields%bc%elemLvl%stencilPos fields%bc%elemLvl%posInNghElems accordingly.
!!
  subroutine remove_solid_in_bc( minLevel, maxLevel, nBCs, nFields, &
    &                            levelPointer, levelDesc, globBC, fields )
    ! ---------------------------------------------------------------------------
    integer, intent(in) :: minLevel, maxLevel, nBCs, nFields
    !> Level pointer
    integer, intent(in) :: levelPointer(:)
    !> Level Descriptor
    type( tem_levelDesc_type ), intent(in) :: levelDesc(minLevel:maxLevel)
    !>
    type( glob_boundary_type ) :: globBC( nBCs )
    !>
    type( mus_field_type ) :: fields( nFields )
    ! ---------------------------------------------------------------------------
    !> number of valid BC elements
    integer :: nValid(minLevel:maxLevel)
    !> positions of only valid elements (non-solid) in bc_elems_type
    integer, allocatable :: posInBCElem(:,:)
    integer :: iBC, iField, iLevel, iElem
    logical :: allWall
    ! ---------------------------------------------------------------------------

    write(dbgUnit(3),*) 'Get into remove_solid_in_bc routine'

    do iBC = 1, nBCs

      !! \todo: do this for both SBB and LIBB?
      allWall = .true.
      do iField = 1, nFields
        if( fields(iField)%bc(iBC)%BC_kind(1:4) /= 'wall' ) then
        ! if( trim(fields(iField)%bc(iBC)%BC_kind) /= 'wall' ) then
          allWall = .false.
          exit
        end if
      end do

      if ( .not. allWall ) then

        write(logUnit(3),"(A)") 'Try to remove solid from BC: '//trim(globBC(iBC)%label)

        ! --------------------------------------------------------------------------------
        allocate( posInBCElem( maxval(globBC(iBC)%nElems), minLevel:maxLevel ) )

        ! check number of valid elements in each BC
        call check_solid_in_bc( minLevel, maxLevel, levelPointer, levelDesc,   &
          &                     globBC(iBC)%nElems, globBC(iBC)%elemLvl,       &
          &                     nValid, posInBCElem )

        ! remove solid in globBC%elemLvl
        call rearrange_bc_elems( minLevel, maxLevel, nValid, posInBCElem, &
          &                      globBC(iBC)%nElems, globBC(iBC)%elemLvl )

        globBC(iBC)%nElems(:) = nValid(:)

        do iField = 1, nFields
          ! re-arrange stencilPos, posInNghElems only if this BC requires neighbors
          if ( fields(iField)%bc(iBC)%nNeighs > 0 ) then
            do iLevel = minLevel, maxLevel

              ! ----------------------------------------------------------------
              do iElem = 1, nValid( iLevel )
                fields(iField)%bc(iBC)%elemLvl(iLevel)%stencilPos( iElem ) = &
                  & fields(iField)%bc(iBC)%elemLvl(iLevel)%stencilPos( posInBCElem(iElem, iLevel) )

                fields(iField)%bc(iBC)%elemLvl(iLevel)%posInNghElems( iElem ) = &
                  & fields(iField)%bc(iBC)%elemLvl(iLevel)%posInNghElems( posInBCElem(iElem, iLevel) )
              end do ! iElem = 1, nValid
              ! ----------------------------------------------------------------

            end do ! iLevel
          end if ! nNeighs > 0
        end do ! iField

        deallocate( posInBCElem )

      end if ! not allWall

    end do ! iBC

    write(dbgUnit(3),*) 'Leave out of remove_solid_in_bc routine'
    write(dbgUnit(3),*) ''

  end subroutine remove_solid_in_bc
! ****************************************************************************** !

! ****************************************************************************** !
  !> Check if a BC is wall or symmetry for all fields
  pure function mus_check_allWall( nFields, fields, iBC ) result ( allWall )
    ! ---------------------------------------------------------------------------
    integer,                intent(in) :: iBC, nFields
    type( mus_field_type ), intent(in) :: fields( nFields )
    logical :: allWall
    ! ---------------------------------------------------------------------------
    integer :: iField
    ! ---------------------------------------------------------------------------

    allWall = .true.
    do iField = 1, nFields
      if( trim(fields(iField)%bc(iBC)%BC_kind) /= 'wall' .or. &
        & trim(fields(iField)%bc(iBC)%BC_kind) /= 'symmetry' ) then
        allWall = .false.
        exit
      end if ! not wall
    end do  ! iField

  end function mus_check_allWall
! ****************************************************************************** !

! ****************************************************************************** !
  !> This routine checks for the existence of symmetric boundaries and 
  !! returns the boundary IDs which are defined as symmetry
  subroutine mus_field_getSymmetricBCs(symmetricBCs, nSymBCs, nBCs, nFields, &
    & field)
    ! ---------------------------------------------------------------------------
    !> number of boundary conditions
    integer, intent(in) :: nBCs
    !> Symmetric boundary ids
    integer, intent(out) :: symmetricBCs(nBCs)
    !> Number of symmetric boundary conditions
    integer, intent(out) :: nSymBCs
    !> number of fields
    integer, intent(in) :: nFields
    !> all fields to access their boundary definitions
    type(mus_field_type), intent(in) :: field(nFields)
    ! ---------------------------------------------------------------------------
    integer :: iField, iBC
    logical :: isSymmetry(nFields)
    ! ---------------------------------------------------------------------------
    nSymBCs = 0
    do iBC = 1, nBCs
      isSymmetry = .false.
      do iField = 1, nFields
        if (trim(field(iField)%bc(iBC)%BC_kind) == 'symmetry') then
          isSymmetry(iField) = .true.
        end if
        if ( any(isSymmetry) ) then
          if (.not. all(isSymmetry) ) then
            call tem_abort('Error: Not all fields boundary label:'&
              &          //trim(field(iField)%bc(iBC)%label)//' are symmetry!')
          end if
          nSymBCs = nSymBCs + 1
          symmetricBCs(nSymBCs) = iBC
        end if
      end do  
    end do

  end subroutine mus_field_getSymmetricBCs  
! ****************************************************************************** !

end module mus_field_module
! ****************************************************************************** !
