! Copyright (c) 2013 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2013-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2013-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2013-2014 Simon Zimny <s.zimny@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) 2017 Sindhuja Budaraju <nagasai.budaraju@student.uni-siegen.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.
! ****************************************************************************** !
!> author: Kannan Masilamani
!! This module contains data type and modules related to musubi
!! lattice to physical unit convertion and vice versa. \n
!! physics data type is global for all scheme, it is defined
!! in the following format: \n
!!```lua
!! physics = { dt = dt_phy, -- physical time step size
!!             rho0 = rho0_phy, -- reference density
!!             temp0 = t_phy -- reference temperature}
!!```
!! Of these quantities, dt is mandetory for conversion.
!! Others can be omitted, thus use default values.
!!
!! To add a new conversion factor, one has to do the following:
!! 1. add this new factor into mus_convertFac_type
!! 2. add the defination of factor inside routine mus_set_convFac
!! 3. add this new factor into routine mus_physics_out
!!
module mus_physics_module

  ! include treelm modules
  use env_module,          only: rk
  use tem_aux_module,      only: tem_abort
  use treelmesh_module,    only: treelmesh_type
  use tem_geometry_module, only: tem_ElemSizeLevel
  use tem_logging_module,  only: logUnit, tem_toStr
  use tem_tools_module,    only: tem_horizontalSpacer

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

  implicit none

  private

  !> Reference coulomb is set to fundamental electrical charge
  real(kind=rk), parameter, public :: coulomb_ref = 1.60217657e-19_rk

  !> Reference mole is set to inverse of Avogadro's constant
  real(kind=rk), parameter, public :: mole_ref = 1e-23_rk / 6.02214129

  !> the boltzmann constant J K^-1
  real(kind=rk), parameter, public :: k_b = 1.38064852e-23_rk

  ! Faraday constant C/mol
  real(kind=rk), parameter, public :: faraday = 96485.3365_rk
  
  ! Ideal gas constant N m / (mol K)
  real(kind=rk), parameter, public :: gasConst_R = 8.3144621_rk

  public :: mus_physics_type
  public :: mus_load_physics
  public :: mus_create_funcStr
  public :: mus_physics_out
  public :: mus_physics_out_conv
  public :: mus_convertFac_type
  public :: mus_set_convFac
  public :: set_values_by_levels
  public :: mus_physics_dump2outUnit
  public :: mus_set_scaleFac

  !> This type contains the converstion factor for derived variables
  !! from lattice to physical.\n
  !! Their inverses can be used to convert physical to lattice units
  !! use reference density to parmeterize kg and reference mole density
  !! to parmeterize mol
  type mus_convertFac_type
    !> length (m) = dx
    real(kind=rk) :: length
    !> time (s) = dt
    real(kind=rk) :: time
    !> velocity(m/s) = dx/dt
    real(kind=rk) :: vel
    !> kinematic viscosity(m^2/s) = dx^2/dt
    real(kind=rk) :: visc
    !> Dynamic viscosity (Pa s) = kg/m/s
    real(kind=rk) :: viscDyna
    !> acceleration(m/s^2) = dx/dt^2
    real(kind=rk) :: accel
    !> Force(N)(kg m/s^2) = rho0*dx^4/dt^2
    real(kind=rk) :: force
    !> Force per unit volume (N/m^3)(kg/s^2/m^2) = rho0*dx/dt^2
    real(kind=rk) :: body_force
    !> Pressure(N/m^2)(kg/m/s^2) = rho0*dx^2/dt^2
    real(kind=rk) :: press
    !> Strain Rate (1/s) = 1/dt
    real(kind=rk) :: strainRate
    !> Energy (N-m) (kg*m^2/s^2) = rho0*dx^5/dt^2
    real(kind=rk) :: energy
    !> mole density(mol/m^3) = mole0/dx^3
    !real(kind=rk) :: moleDens
    !> Charge density (C/m^3) = Coulomb0/dx^3
    real(kind=rk) :: chargeDens
    !> Current density (C/s/m^2) = Coulomb0/dt/dx^2
    real(kind=rk) :: currentDens
    !> mole flux(mol/m^2/s) = moleDen0*dx/dt
    real(kind=rk) :: moleFlux
    !> mass flux(kg/m^2/s) = rho0*dx/dt
    real(kind=rk) :: flux
    !> diffusivity(m^2/s) = dx^2/dt
    real(kind=rk) :: diffusivity
    !> faraday (C/mol) = coulomb0/moleDens0/dx^3
    real(kind=rk) :: faraday
    !> gas constant (J/mol/K) (N m/mol/K) (kg m^2/s^2/mol/K) 
    !! = rho0*dx^5/dt^2/mole0/temp0
    real(kind=rk) :: gasConst
    !> Potential (V) (kg m^2/(C*S^2))
    real(kind=rk) :: potential
  end type mus_convertFac_type

  !> This type contains the basic SI units used to convert
  !! lattice to physical unit and vice versa.
  !! keep reference mass density, mole density and molecular weight
  !! same for all levels.
  type mus_physics_type
    !> needed to check if physics table is defined
    logical :: active = .false.
    !> reference length - discretization size of the coarsest level
    !! SI unit - meter
    real(kind=rk) :: dx = -1.0_rk
    real(kind=rk), allocatable :: dxLvl(:)
    !> reference time - time discretization for discretization size 
    !! of the coarsest level
    !! SI unit - seconds
    real(kind=rk) :: dt = -1.0_rk
    real(kind=rk), allocatable :: dtLvl(:)
    !> reference physical mass density
    !! SI unit - kg/m^3
    real(kind=rk) :: rho0 = -1.0_rk
    !> reference physical mole density
    !! SI unit - mol/m^3
    real(kind=rk) :: moleDens0 = -1.0_rk
    !> reference molecular weight
    !! SI unit - kg/mol
    real(kind=rk) :: molWeight0 = -1.0_rk
    !> reference temperature
    !! SI unit - Kelvin
    real(kind=rk) :: temp0 = -1.0_rk
    !> reference fundamental electrical charge
    !! SI unit - Coulomb
    real(kind=rk) :: coulomb0 = -1.0_rk
    !> mole is defined by inverse of Avogadro Constant 
    !! Avogadro Constant = 6.02214129e23 [1/mol]
    real(kind=rk) :: mole0 = -1.0_rk
    !> reference mass in kg derived from density or moleweight
    !! SI unit :: kg
    real(kind=rk) :: mass0 = -1.0_rk

    !> Level-wise conversion factor for derived variables
    !! size: minLevel:maxLevel
    !! allocated in mus_load_physics
    !! \todo KM: conversion factor should not be level-dependent.
    !! it should be same for all levels, the lattice dx and dt for each level
    !! must be considered to scale variables in multilevel.
    !! Implemented force, visc, etc using dtL according to formula
    !! Introduced lattice speed variable: dx/dt for each level
    !! it should be same for all level for acoustic scaling and different
    !! for diffusive scaling
    type( mus_convertFac_type ), allocatable :: fac(:)

    !> Pressure (strain rate) over level scale factor.
    !! This factor is meant to convert pressure in LB unit on source level to
    !! the required pressure on target level.
    !! It is mainly used in interpolation routine.
    !! It is allocated as: allocate(pFac( minLevel:maxLevel, minLevel:maxLevel))
    !! It is allocated and initialized in routine: mus_set_scaleFac
    !! How to use it in the code:
    !! pTargetLevel = pSourceLevel * pFac( sourceLevel, targetLevel )
    real(kind=rk), allocatable :: pFac(:,:)

    !> Velocity over level scale factor
    !! Its usage is the same as pressure scale factor
    real(kind=rk), allocatable :: vFac(:,:)

    !> Strain rate over level scale factor
    !! Its usage is the same as pressure scale factor
    real(kind=rk), allocatable :: sFac(:,:)

  end type mus_physics_type

  contains

  ! **************************************************************************** !
  !> This routine loads the physics table from musubi config file
  !!
  subroutine mus_load_physics( me, conf, tree, scaleFactor, dtRef, dxRef )
    ! ---------------------------------------------------------------------------
    !> physics type
    type( mus_physics_type ), intent(out) :: me
    !> flu state
    type( flu_State ) :: conf
    !> global treelm mesh
    type( treelmesh_type), intent(in) :: tree
    !> scaling factor: diffusive -> 4; acoustic -> 2
    integer, intent(in) :: scaleFactor
    !> reference time step if none 
    real(kind=rk), optional, intent(in) :: dtRef
    !> reference spacestep if none 
    real(kind=rk), optional, intent(in) :: dxRef
    ! ---------------------------------------------------------------------------
    integer :: thandle
    integer :: iError
    real(kind=rk) :: dt_Ref, dx_ref ! reference step if none
    ! ---------------------------------------------------------------------------
    if( present( dtRef )) then
      dt_ref = dtRef
    else
      dt_ref = 1._rk
    end if

    if( present( dxRef )) then
      dx_ref = dxRef
    else
      dx_ref = 1._rk
    end if

    call tem_horizontalSpacer(fUnit = logUnit(1))
    write(logUnit(1),*) ' Loading physics table ...'
    call aot_table_open( L=conf, thandle=thandle, key='physics' )

    if ( thandle > 0 ) then
      me%active = .true.
      ! Found a physics table.
      ! Activate the conversion
      write(logUnit(1), *) ' Physics table is defined.'
      write(logUnit(5), *) ' All values are expressed in physical units.'

      ! reference dx is always the coarsest level in the tree
      me%dx = tem_ElemSizeLevel( tree, tree%global%minLevel )

      ! load dt
      call aot_get_val( L = conf, thandle = thandle, key = 'dt',               &
        &               val = me%dt, ErrCode = iError )
      if (btest(iError, aoterr_Fatal)) then
        write(logUnit(1),*)'FATAL Error occured, while retrieving dt.'
        if (btest(iError, aoterr_NonExistent)) then
          write(logUnit(1),*)'Time step definition not found. Setting to dt=1'
          me%dt = dt_ref
        end if
        if (btest(iError, aoterr_WrongType)) then
          write(logUnit(1),*)'Variable has wrong type!'
          write(logUnit(1),*)'STOPPING'
          call tem_abort()
        endif
      end if

      ! define mole before loading density because if molecular weight is defined
      ! to compute reference mass than we need reference mole.
      ! try to load reference mole density, if not defined then set reference
      ! mole to inverse of avogadro's constant
      call aot_get_val( L = conf, thandle = thandle, key = 'moleDens0',        &
        &               val = me%moleDens0, ErrCode = iError )
      if (btest(iError, aoterr_Fatal)) then
        write(logUnit(7),*) 'No value given for moleDens0.'
        if (btest(iError, aoterr_WrongType)) then
          write(logUnit(1),*)'Error! moleDens0 has wrong type!'
          call tem_abort()
        end if
      end if

      ! reference mole density not defined so try to load reference mole
      if (btest(iError, aoterr_NonExistent)) then
        write(logUnit(7),*)'WARNING: Reference moleDens0 is not found. '
        write(logUnit(7),*)'Loading reference mole0(mol):'
        call aot_get_val( L = conf, thandle = thandle, key = 'mole0',          &
          &               val = me%mole0, ErrCode = iError )
        if (btest(iError, aoterr_Fatal)) then
          write(logUnit(7),*)'FATAL Error occured, while retrieving mole0.'
          if (btest(iError, aoterr_WrongType)) then
            write(logUnit(1),*)'Variable has wrong type!'
            call tem_abort()
          end if
        end if
        ! reference mole is also not defined set inverse of Avogadro's constant
        if (btest(iError, aoterr_NonExistent)) then
          write(logUnit(7),*)'WARNING: Reference mole0 is not found. '
          write(logUnit(7),*)"Setting reference mole to inverse of Avogadro's" &
            &                //"constant."
          me%mole0 = mole_ref
        endif
        ! derive moleDens0 from mole0
        me%moleDens0 = me%mole0 / me%dx**3
      else ! if defined moleDens0 derive mole from moledensity
        me%mole0 = me%moleDens0 * me%dx**3
      endif

      ! try to load reference density, If not defined then try to load reference
      ! molecular weight, if not defined then try to load reference mass in 'kg'
      ! if that also is not defined prompt an error message
      ! load rho0
      call aot_get_val( L = conf, thandle = thandle, key = 'rho0',             &
        &               val = me%rho0, ErrCode = iError )
      if (btest(iError, aoterr_Fatal)) then
        write(logUnit(3),*) 'No value given for rho0.'
        if (btest(iError, aoterr_WrongType)) then
          write(logUnit(1),*) 'rho0 has wrong type!'
          call tem_abort()
        end if
      end if

      ! reference mass density not defined so try to load reference
      ! molecular weight
      if (btest(iError, aoterr_NonExistent)) then
        write(logUnit(1),*)'WARNING: Reference mass density rho0 is not found.'
        write(logUnit(1),*)'Loading reference molecular weight '               &
          &              //'molWeight0(kg/mol):'
        call aot_get_val( L = conf, thandle = thandle, key = 'molWeight0',     &
          &               val = me%molWeight0, ErrCode = iError )
        if (btest(iError, aoterr_Fatal)) then
          write(logUnit(3),*)'No value given for molWeight0.'
          if (btest(iError, aoterr_WrongType)) then
            write(logUnit(1),*) 'molWeight has wrong type!'
            call tem_abort()
          end if
        end if
        ! if reference molecular weight is not defined than load reference mass
        if (btest(iError, aoterr_NonExistent)) then
          write(logUnit(5),*)'WARNING: Reference molecular weight molWeight0 ' &
          &                //'is not found. '
          write(logUnit(5),*)'Loading reference mass0(kg):'
          call aot_get_val( L = conf, thandle = thandle, key = 'mass0',        &
            &               val = me%mass0, ErrCode = iError )
          if (btest(iError, aoterr_Fatal)) then
            write(logUnit(3),*)'No value given for mass0'
            if (btest(iError, aoterr_WrongType)) then
              write(logUnit(1),*)'Variable has wrong type!'
              call tem_abort()
            end if
          end if
          ! if reference mass is also not defined prompt an error message
          if (btest(iError, aoterr_NonExistent)) then
            write(logUnit(1),*) 'ERROR: Unable to obtain reference mass0(kg)'
            write(logUnit(1),*) "Solution: Provide reference density 'rho0' "  &
              &               //"in kg/m^3"
            write(logUnit(1),*) "or reference molecular weight 'molWeight0' "   &
              &               //"in kg/mol"
            write(logUnit(1),*) "or reference mass 'mass0' in kg"
            call tem_abort()
          else 
          ! reference mass is defined derive density and molWeight0
            me%rho0 = me%mass0/me%dx**3
            ! same as rho0/moleDens0
            ! molecular Weight in independent of dx
            me%molWeight0 = me%mass0/me%mole0
          endif
        else
        ! reference molecular weight is defined derive mass and density
          me%mass0 = me%molWeight0 * me%mole0
          me%rho0 = me%mass0/me%dx**3
        endif
      else
        ! reference density is defined derive mass and molWeigh
        me%mass0 = me%rho0 * me%dx**3
        ! same as rho0/moleDens0
        ! molecular Weight in independent of dx
        me%molWeight0 = me%mass0 / me%mole0
      endif

      ! try to load coloumb, if not defined set to fundamental electrical charge
      call aot_get_val( L = conf, thandle = thandle, key = 'coulomb0',         &
        &               val = me%coulomb0, ErrCode = iError )
      if (btest(iError, aoterr_Fatal)) then
        write(logUnit(3),*)'No value given for coulomb0'
        if (btest(iError, aoterr_WrongType)) then
          write(logUnit(1),*) 'coulomb0 has wrong type!'
          call tem_abort()
        end if
      end if
      if (btest(iError, aoterr_NonExistent)) then
        write(logUnit(4),*) 'WARNING: Reference coulomb0 is not found. '
        write(logUnit(4),*) "Setting reference coulomb0 to fundamental" &
          &                 //" electrical charge."
        me%coulomb0 = coulomb_ref
        !me%coulomb0 = me%rho0*me%dx**5 / me%dt**2
      endif

      ! load temp0
      call aot_get_val(L = conf, thandle = thandle, key = 'temp0',             &
        &              val = me%temp0, default = 1.0_rk, ErrCode = iError)
      if (btest(iError, aoterr_WrongType)) then
        write(logUnit(1),*)'temp0 has wrong type! Stopping!'
        call tem_abort()
      end if
    else
      ! No physics table defined.
      me%dt         = dt_ref
      me%dx         = dx_ref
      me%rho0       = 1._rk
      me%moleDens0  = 1._rk
      me%molWeight0 = 1._rk
      me%mass0      = 1._rk
      me%mole0      = 1._rk
      me%temp0      = 1._rk
      me%coulomb0   = 1._rk
      write(logUnit(1),*) ' Physics table is not defined.'
      write(logUnit(5),*) ' All values are expressed in Lattice Units.'
    end if

    call aot_table_close( L=conf, thandle=thandle )

    ! Assign dx and dt for each level
    allocate(me%dxLvl( tree%global%minLevel:tree%global%maxLevel ))
    allocate(me%dtLvl( tree%global%minLevel:tree%global%maxLevel ))
    allocate(  me%fac( tree%global%minLevel:tree%global%maxLevel ))

    ! dx and dt is set according scaling type
    me%dxLvl( tree%global%minLevel:tree%global%maxLevel ) =  &
      & set_values_by_levels( me%dx, tree%global%minLevel,  &
      &                              tree%global%maxLevel, 2 )

    me%dtLvl( tree%global%minLevel:tree%global%maxLevel ) =  &
      & set_values_by_levels( me%dt, tree%global%minLevel,             &
      &                              tree%global%maxLevel, scaleFactor )

    ! compute and store conversion factors in converstionFac type
    call mus_set_convFac( me = me, minLevel = tree%global%minLevel, &
      &                            maxLevel = tree%global%maxLevel )
    ! set scale factor
    call mus_set_scaleFac( me, tree%global%minLevel, tree%global%maxLevel )

    call mus_physics_dump2outUnit( me, logUnit(4), tree%global%minLevel,  &
      &                           tree%global%maxLevel )

    call tem_horizontalSpacer(fUnit = logUnit(1))

  end subroutine mus_load_physics
! ****************************************************************************** !


! ****************************************************************************** !
  !> This routine computed conversion factors for lattice to physical units.
  !! inverse of this factors can be used to convert from physical to lattice
  !! units.\n
  !! use reference density to parmeterize kg and reference mole density
  !! to parmeterize mol.\n
  !! Multiply these factors with the LB quantity to get the physical quantity
  !! Divide the physical quantity by these factors to get the LB units.
  subroutine mus_set_convFac( me, minLevel, maxLevel )
    ! ---------------------------------------------------------------------------
    type( mus_physics_type ), intent(inout) :: me !< physics type
    integer, intent(in) :: minLevel
    integer, intent(in) :: maxLevel
    ! ---------------------------------------------------------------------------
    integer :: iLevel
    ! ---------------------------------------------------------------------------

    do iLevel = minLevel, maxLevel
      ! length
      me%fac( iLevel )%length = me%dxLvl( iLevel )
      ! time
      me%fac( iLevel )%time = me%dtLvl( iLevel )
      ! velocity
      me%fac( iLevel )%vel = me%dxLvl( iLevel )/me%dtLvl( iLevel )
      ! kinematic viscosity
      me%fac( iLevel )%visc = me%dxLvl( iLevel )**2/me%dtLvl( iLevel )
      ! dynamic viscosity
      me%fac( iLevel )%viscDyna = me%rho0 * me%fac( iLevel )%visc
      ! acceleration
      me%fac( iLevel )%accel = me%dxLvl( iLevel )/me%dtLvl( iLevel )**2
      ! force
      me%fac( iLevel )%force = me%rho0*me%dxLvl( iLevel )**4                   &
        &                    / me%dtLvl( iLevel )**2
      ! body_force
      me%fac( iLevel )%body_force = me%rho0*me%dxLvl( iLevel )                 &
        &                         / me%dtLvl( iLevel )**2
      ! pressure or shear stress
      me%fac( iLevel )%press = me%rho0 * me%dxLvl( iLevel )**2                 &
        &                    / me%dtLvl( iLevel )**2
      ! strain rate
      me%fac( iLevel )%strainRate = 1._rk / me%dtLvl( iLevel )
      ! Energy
      me%fac( iLevel )%energy = me%rho0*me%dxLvl( iLevel )**5                  &
        &                     / me%dtLvl( iLevel )**2
      ! Mass
      !me%fac( iLevel )%mass = me%rho0*me%dxLvl( iLevel )**3
      ! Molar mass or Molecular weight
      !me%fac( iLevel )%molWeigh = me%fac( iLevel )%mass/me%mole0
      ! mole density
      !me%fac( iLevel )%moleDens = me%mole0/me%dxLvl( iLevel )**3
      ! mole flux (mol/m^2/s) (moleDens*velocity)
      me%fac( iLevel )%moleFlux = me%moleDens0*me%dxLvl( iLevel ) &
        &                       / me%dtLvl( iLevel )
      ! Charge density (C/m^3) = Coulomb0/dx^3
      me%fac( iLevel )%chargeDens = me%coulomb0 / me%dxLvl( iLevel )**3
      ! Current density (C/s/m^2) = Coulomb0/dt/dx^2
      me%fac( iLevel )%currentDens = me%coulomb0 / me%dxLvl( iLevel )**2 &
        &                          / me%dtLvl( iLevel )
      ! mass flux (kg/m^2/s)  (density*velocity)
      me%fac( iLevel )%flux = me%rho0*me%dxLvl( iLevel )/me%dtLvl( iLevel )
      !diffusivity
      me%fac( iLevel )%diffusivity = me%dxLvl( iLevel )**2/me%dtLvl( iLevel )
      !faraday
      me%fac( iLevel )%faraday = me%coulomb0/me%mole0
      !gas constant
      me%fac( iLevel )%gasConst = me%rho0*me%dxLvl( iLevel )**5     &
        &                       / me%dtLvl( iLevel )**2 / me%mole0  &
        &                       / me%temp0
      ! Potential
      me%fac( iLevel )%potential = me%rho0*me%dxLvl( iLevel )**5  &
        &                     / me%dtLvl( iLevel )**2 / me%coulomb0

    end do

  end subroutine mus_set_convFac
! ****************************************************************************** !


! ****************************************************************************** !
  !> This routine creates musubi specific lua function to compute dx and dt.
  !!
  subroutine mus_create_funcStr( fun_str )
    ! ---------------------------------------------------------------------------
    !> This string contains lua functions to compute dt from visocosity or
    !! velocity
    character(len=*) :: fun_str
    ! ---------------------------------------------------------------------------
    character(len=256) :: dtVel, dtOmega, dtVisc, dxLevel, dxnL, omega_dt,     &
      &                   vel_dt
    ! ---------------------------------------------------------------------------
    write(dxLevel,*) 'function getdxFromLevel(t) '//                           &
      &              ' return t.len_bnd/2^t.level end'

    write(dxnL,*) 'function getdxFromnL(t)'//                                  &
      &           ' return t.len_p/t.nL end'

    write(dtVel,*) 'function getdtFromVel(t) return t.dx*t.u_l/t.u_p end'

    write(dtOmega,*) 'function getdtFromOmega(t) '//                           &
      &              ' return ((1.0/t.omega - 0.5)/3.0)*t.dx*t.dx/t.nu_p end'

    write(dtVisc,*) ' function getdtFromVisc(t)'//                             &
      &             ' return t.dx*t.dx*t.nu_l/t.nu_p end'

    write(omega_dt,*) 'function getOmegaFromdt(t) '//                          &
      &               ' return 1.0/(3.0*t.nu_p*t.dt/(t.dx*t.dx) + 0.5) end'

    write(vel_dt, *) 'function getVelFromdt(t) '//                             &
      &              ' return t.u_p*t.dt/t.dx end'

    write(fun_str,*) trim(dtVel)//' '//trim(dtOmega)//' '//trim(dtVisc)        &
      & //' '//trim(dxLevel)//' '//trim(dxnL)//' '//trim(omega_dt)             &
      & //' '//trim(vel_dt)

  end subroutine mus_create_funcStr
! ****************************************************************************** !


! ****************************************************************************** !
  !> This routine write reference physics parameters into solver specific 
  !! string in lua format.
  !!
  !! This dumped table is loaded back using mus_load_physics
  subroutine mus_physics_out( me, conf )
    ! ---------------------------------------------------------------------------
    type( mus_physics_type ), intent(in) :: me !< physics type
    type( aot_out_type ) :: conf
    ! ---------------------------------------------------------------------------

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

    call aot_out_val( put_conf = conf, vname = 'dt', val = me%dt )
    call aot_out_val( put_conf = conf, vname = 'rho0', val = me%rho0 )
    call aot_out_val( put_conf = conf, vname = 'moleDens0', val = me%moleDens0 )
    call aot_out_val( put_conf = conf, vname = 'molWeight0', &
      &               val = me%molWeight0 )
    call aot_out_val( put_conf = conf, vname = 'temp0', val = me%temp0 )
    call aot_out_val( put_conf = conf, vname = 'coulomb0', val = me%coulomb0 )

    call aot_out_close_table( put_conf = conf )

  end subroutine mus_physics_out
! ****************************************************************************** !


! ****************************************************************************** !
  !> This routine write physics convert factor into solver specific string in
  !! lua format.
  !! use reference density to parmeterize kg and reference mole density
  !! to parmeterize mol.
  subroutine mus_physics_out_conv( me, conf, minLevel, maxLevel )
    ! ---------------------------------------------------------------------------
    type( mus_physics_type ), intent(in) :: me !< physics type
    type( aot_out_type ) :: conf
    integer, intent(in) :: minLevel
    integer, intent(in) :: maxLevel
    ! ---------------------------------------------------------------------------
    integer :: iLevel
    ! ---------------------------------------------------------------------------

    call aot_out_open_table( put_conf = conf, tname = 'physics' )
    do iLevel = minLevel, maxLevel
      call aot_out_open_table( put_conf = conf )
      call aot_out_val( put_conf = conf, vname = 'length',                     &
        &               val = me%fac(iLevel)%length )
      call aot_out_val( put_conf = conf, vname = 'time',                       &
        &               val = me%fac(iLevel)%time )
      call aot_out_val( put_conf = conf, vname = 'density',                    &
        &               val = me%rho0 )
      call aot_out_val( put_conf = conf, vname = 'vel',                        &
        &               val = me%fac(iLevel)%vel )
      call aot_out_val( put_conf = conf, vname = 'visc',                       &
        &               val = me%fac(iLevel)%visc )
      call aot_out_val( put_conf = conf, vname = 'viscDyna',                   &
        &               val = me%fac(iLevel)%viscDyna )
      call aot_out_val( put_conf = conf, vname = 'accel',                      &
        &               val = me%fac(iLevel)%accel )
      call aot_out_val( put_conf = conf, vname = 'force',                      &
        &               val = me%fac(iLevel)%force )
      call aot_out_val( put_conf = conf, vname = 'body_force',                 &
        &               val = me%fac(iLevel)%body_force )
      call aot_out_val( put_conf = conf, vname = 'press',                      &
        &               val = me%fac(iLevel)%press )
      call aot_out_val( put_conf = conf, vname = 'strainRate',                 &
        &               val = me%fac(iLevel)%strainRate )
      call aot_out_val( put_conf = conf, vname = 'energy',                     &
        &               val = me%fac(iLevel)%energy )
      call aot_out_val( put_conf = conf, vname = 'temp',                       &
        &               val = me%temp0 )
!      call aot_out_val( put_conf = conf, vname = 'mass',                       &
!        &               val = me%fac(iLevel)%mass )
!      call aot_out_val( put_conf = conf, vname = 'mole',                       &
!        &               val = me%mole0 )
      call aot_out_val( put_conf = conf, vname = 'moleDens',                   &
        &               val = me%moleDens0 )
      call aot_out_val( put_conf = conf, vname = 'molWeight',                  &
        &               val = me%molWeight0 )
      call aot_out_val( put_conf = conf, vname = 'coulomb',                    &
        &               val = me%coulomb0 )
      call aot_out_val( put_conf = conf, vname = 'moleFlux',                   &
        &               val = me%fac(iLevel)%moleFlux )
      call aot_out_val( put_conf = conf, vname = 'chargeDens',                 &
        &               val = me%fac(iLevel)%chargeDens )
      call aot_out_val( put_conf = conf, vname = 'currentDens',                &
        &               val = me%fac(iLevel)%currentDens )
      call aot_out_val( put_conf = conf, vname = 'flux',                       &
        &               val = me%fac(iLevel)%flux )
      call aot_out_val( put_conf = conf, vname = 'diffusivity',                &
        &               val = me%fac(iLevel)%diffusivity )
      call aot_out_val( put_conf = conf, vname = 'faraday',                    &
        &               val = me%fac(iLevel)%faraday )
      call aot_out_val( put_conf = conf, vname = 'gasConst',                   &
        &               val = me%fac(iLevel)%gasConst )
      call aot_out_close_table( put_conf = conf )
    end do

    call aot_out_close_table( put_conf = conf )

  end subroutine mus_physics_out_conv
! ****************************************************************************** !

  ! **************************************************************************** !
  pure function set_values_by_levels( valMinLevel, minLevel, maxLevel, scaleFac )&
    &                                                       result( values )
    ! ---------------------------------------------------------------------------
    real(kind=rk), intent(in) :: valMinLevel !< value at min level
    integer,       intent(in) :: minLevel
    integer,       intent(in) :: maxLevel
    integer,       intent(in) :: scaleFac !< scale factor between levels
    real(kind=rk) :: values(minLevel:maxLevel) !< return value
    ! ---------------------------------------------------------------------------
    integer :: iLevel
    ! ---------------------------------------------------------------------------

    values( minLevel:maxLevel ) = valMinLevel / &
      & real([(scaleFac**(iLevel-minLevel), iLevel=minLevel, maxLevel)],rk)

  end function set_values_by_levels
  ! **************************************************************************** !

  ! **************************************************************************** !
  subroutine mus_physics_dump2outUnit( me, outUnit, minLevel, maxLevel )
    ! ---------------------------------------------------------------------------
    type( mus_physics_type ), intent(in) :: me
    integer, intent(in) :: outUnit
    integer, intent(in) :: minLevel
    integer, intent(in) :: maxLevel
    ! ---------------------------------------------------------------------------
    integer :: iLevel

    write(outUnit,'(A)')       'Reference physical quantities:'
    write(outUnit,"(A)") '  rho0       = '//trim(tem_toStr(me%rho0))
    write(outUnit,"(A)") '  mass0      = '//trim(tem_toStr(me%mass0))
    write(outUnit,"(A)") '  mole0      = '//trim(tem_toStr(me%mole0))
    write(outUnit,"(A)") '  moleDens0  = '//trim(tem_toStr(me%moleDens0))
    write(outUnit,"(A)") '  molWeight0 = '//trim(tem_toStr(me%molWeight0))
    write(outUnit,"(A)") '  Temp0      = '//trim(tem_toStr(me%temp0))
    write(outUnit,"(A)") '  coulomb0   = '//trim(tem_toStr(me%coulomb0))

    do iLevel = minLevel, maxLevel
      write(outUnit,"(A,I0)")    'Conversion factors on level: ', iLevel
      write(outUnit,"(A)") '  dx            = '//trim(tem_toStr(me%dxLvl( iLevel )))
      write(outUnit,"(A)") '  dt            = '//trim(tem_toStr(me%dtLvl( iLevel )))
      write(outUnit,"(A)") '  pressure      = '//trim(tem_toStr(me%fac( iLevel )%press))
      write(outUnit,"(A)") '  velocity      = '//trim(tem_toStr(me%fac( iLevel )%vel))
      write(outUnit,"(A)") '  kine. visc.   = '//trim(tem_toStr(me%fac( iLevel )%visc))
      write(outUnit,"(A)") '  dyn. visc.    = '//trim(tem_toStr(me%fac( iLevel )%viscDyna))
      write(outUnit,"(A)") '  acceleration  = '//trim(tem_toStr(me%fac( iLevel )%accel))
      write(outUnit,"(A)") '  force         = '//trim(tem_toStr(me%fac( iLevel )%force))
      write(outUnit,"(A)") '  body force    = '//trim(tem_toStr(me%fac( iLevel )%body_force))
      write(outUnit,"(A)") '  strain rate   = '//trim(tem_toStr(me%fac( iLevel )%strainRate))
      write(outUnit,"(A)") '  energy        = '//trim(tem_toStr(me%fac( iLevel )%energy))
      write(outUnit,"(A)") '  moleFlux      = '//trim(tem_toStr(me%fac( iLevel )%moleFlux))
      write(outUnit,"(A)") '  chargeDens    = '//trim(tem_toStr(me%fac( iLevel )%chargeDens))
      write(outUnit,"(A)") '  currentDens   = '//trim(tem_toStr(me%fac( iLevel )%currentDens))
      write(outUnit,"(A)") '  massflux      = '//trim(tem_toStr(me%fac( iLevel )%flux))
      write(outUnit,"(A)") '  diffusivity   = '//trim(tem_toStr(me%fac( iLevel )%diffusivity))
      write(outUnit,"(A)") '  faraday       = '//trim(tem_toStr(me%fac( iLevel )%faraday))
      write(outUnit,"(A)") '  gasConst      = '//trim(tem_toStr(me%fac( iLevel )%gasConst))
      write(outUnit,"(A)") '  potential     = '//trim(tem_toStr(me%fac( iLevel )%potential))
      write(outUnit,"(A)") ''
    end do

    write(outUnit,"(A)")   'Over level scale factor: '
    do iLevel = minLevel, maxLevel
      write(outUnit,"(A)") '    pressure: '//trim(tem_toStr(me%pFac(iLevel,minLevel:maxLevel),','))
      write(outUnit,"(A)") '    velocity: '//trim(tem_toStr(me%vFac(iLevel,minLevel:maxLevel),','))
      write(outUnit,"(A)") ' strain rate: '//trim(tem_toStr(me%sFac(iLevel,minLevel:maxLevel),','))
    end do

  end subroutine mus_physics_dump2outUnit
  ! **************************************************************************** !

  ! **************************************************************************** !
  subroutine mus_set_scaleFac( me, minLevel, maxLevel )
    ! ---------------------------------------------------------------------------
    type( mus_physics_type ), intent(inout) :: me
    integer, intent(in) :: minLevel
    integer, intent(in) :: maxLevel
    ! ---------------------------------------------------------------------------
    integer :: iLevel
    ! ---------------------------------------------------------------------------

    allocate( me%pFac( minLevel:maxLevel, minLevel:maxLevel ) )
    allocate( me%vFac( minLevel:maxLevel, minLevel:maxLevel ) )
    allocate( me%sFac( minLevel:maxLevel, minLevel:maxLevel ) )

    do iLevel = minLevel, maxLevel
      me%pFac( :, iLevel ) = me%fac(:)%press      / me%fac(iLevel)%press
      me%vFac( :, iLevel ) = me%fac(:)%vel        / me%fac(iLevel)%vel
      me%sFac( :, iLevel ) = me%fac(:)%strainRate / me%fac(iLevel)%strainRate
    end do

  end subroutine mus_set_scaleFac
  ! **************************************************************************** !

end module mus_physics_module
! ****************************************************************************** !
