! Copyright (c) 2022 Kannan Masilamani <kannan.masilaman@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 contains turbulent wall function type and routines to calculate
!! friction velocity and stream-wise velocity component.
module mus_turb_wallFunc_module

  ! include treelm modules
  use env_module,               only: rk, long_k, labelLen
  use tem_aux_module,           only: tem_abort
  use tem_logging_module,       only: logUnit
  use tem_float_module,         only: operator(.fle.)

  ! include aotus modules
  use aotus_module,     only: flu_State, aoterr_Fatal, aoterr_NonExistent, &
    &                         aoterr_WrongType, aot_get_val

  implicit none

  private

  public :: mus_turb_wallFunc_type
  public :: mus_load_turb_wallFunc

  !> Contains friction velocity and turbulent viscosity on boundary elements on
  !! each level
  type mus_turb_wallFunc_data_type
    !> Turbulent viscosity on boundary element computed using mixing length
    !! formulation in lattice unit
    !! nu_t = (vonKarman*distToBnd)**2 * |du/dy|
    real(kind=rk), allocatable :: tVisc(:)

    !> Friction velocity on first neighbor in normal direction in lattice unit
    !! computed from wall model
    real(kind=rk), allocatable :: velTau(:)

    !> Distance to boundary from first fluid in normal direction
    !! in lattice unit.
    real(kind=rk), allocatable :: distToBnd(:)

    !> Distance to boundary from first fluid neighbor in normal direction
    !! in lattice unit.
    real(kind=rk), allocatable :: neighDistToBnd(:)
  end type mus_turb_wallFunc_data_type

  !> Contains function pointers to compute fuction velocity and stream-wise
  !! velocity component
  type mus_turb_wallFunc_type
    !> is true if wall function is active
    logical :: isActive = .false.

    !> Wall model function
    character(len=labelLen) :: wall_func

    !> Nonlinear solver type
    character(len=labelLen) :: nonlinear_solver

    !> Von-Karman constant. Default = 0.4_rk
    real(kind=rk) :: vonKarman = 0.4_rk

    !> Use vanDriest damping function to damp turbulent viscosity
    logical :: useVanDriest = .true.

    !> Contains data computed in turbulent wall bc routine on each level
    type(mus_turb_wallFunc_data_type), allocatable :: dataOnLvl(:)

    !> Function pointer to compute friction velocity
    procedure(mus_proc_calcFricVel), pointer, nopass :: calcFricVel => null()

    !> Function pointer to compute strean-wise velocity component
    procedure(mus_proc_calcStreamWiseVel), pointer, nopass :: &
      & calcStreamWiseVel => null()
  end type mus_turb_wallFunc_type

  !> Interface definition for the turbulent wall bc routines
  abstract interface
    !> This abstract interface defines the interface to calculate turbulent wall
    !! friction velocity from given velocity and distance to boundary.
    !! All inputs and output are in lattice units.
    subroutine mus_proc_calcFricVel(velTau, velSW, distToBnd, &
      &                             viscKine, nElems)
      import :: rk
      !> Friction velocity computed from wall model.
      !! it is inout to provide velTau from previous timestep as initial velTau
      !! for fixed-point or Newton iteration solver
      real(kind=rk), intent(inout) :: velTau(:)
      !> Stream-wise velocity component from which friction velocity is computed
      real(kind=rk), intent(in) :: velSW(:)
      !> Distance to the boundary in the discrete normal direction
      real(kind=rk), intent(in) :: distToBnd(:)
      !> Kinematic viscosity
      real(kind=rk), intent(in) :: viscKine(:)
      !> Number of elements in input and output arrays
      integer, intent(in) :: nElems
    end subroutine mus_proc_calcFricVel

    !> This abstract interface defines the interface to calculate stream-wise
    !! velocity component from friction velocity and distance to boundary.
    !! All inputs and output are in lattice units.
    subroutine mus_proc_calcStreamWiseVel(velSW, velTau, distToBnd, &
      &                                   viscKine, nElems)
      import :: rk
      !> Stream-wise velocity component from wall model
      real(kind=rk), intent(out) :: velSW(:)
      !> Friction velocity computd from wall model
      real(kind=rk), intent(in) :: velTau(:)
      !> Distance to the boundary in the discrete normai direction
      real(kind=rk), intent(in) :: distToBnd(:)
      !> Kinematic viscosity
      real(kind=rk), intent(in) :: viscKine(:)
      !> Number of elements in input and output arrays
      integer, intent(in) :: nElems
    end subroutine mus_proc_calcStreamWiseVel
  end interface

  !! Constant parameters for Implicit equation solver
  integer, parameter :: imEq_nIter = 1000
  real(kind=rk), parameter :: imEq_tol = 1e-10

  !! Constant parameters for Reichardt's law
  real(kind=rk), parameter :: reivonKA = 0.4_rk
  real(kind=rk), parameter :: oneOvervonKA = 1.0_rk/reivonKA
  real(kind=rk), parameter :: reiC  = 7.8_rk
  real(kind=rk), parameter :: reiB1 = 11.0_rk
  real(kind=rk), parameter :: reiB2 = 3.0_rk

  !! constants for Werner & Wengle model
  !! http://trace-portal.de/userguide/trace/page_wallFunctionWernerWengle.html
  real(kind=rk), parameter :: ww_C_m = 8.3_rk
  real(kind=rk), parameter :: ww_m = 1.0_rk/7.0_rk
  real(kind=rk), parameter :: ww_uLmt_yPlus = ww_C_m**(1.0_rk/(1.0_rk-ww_m))
  real(kind=rk), parameter :: ww_uLmt_vel = WW_uLmt_yPlus**2.0_rk
  real(kind=rk), parameter :: ww_onePlusM = 1.0_rk+ww_m
  real(kind=rk), parameter :: ww_onePlusMInv = 1.0_rk/ww_onePlusM
  real(kind=rk), parameter :: ww_mDivOnePlusM = ww_m/ww_onePlusM

  !! Werner & Wengle two layer model does not cover buffer layer so
  !! three different layer model from Schmidt is used to compute stream-wise
  !! velocity at first element.
  !! Viscous sublayer < 5
  !! buffer layer > 5 and < 30
  !! log layer > 30
  real(kind=rk), parameter :: sc_lLmt = 5.0_rk ! lower limit of buffer layer
  real(kind=rk), parameter :: sc_uLmt = 30.0_rk ! upper limit for buffer layer
  real(kind=rk), parameter :: sc_vonKarman = 0.4_rk
  real(kind=rk), parameter :: sc_a = ((1.0_rk/sc_vonKarman)*log(30.0_rk) &
    &                              + 0.2_rk)/log(6.0_rk)
  real(kind=rk), parameter :: sc_b = 5.0_rk-sc_a*log(5.0_rk)

contains

  ! -------------------------------------------------------------------------- !
  !> This routine loads wall model and nonlinear solver type for nonlinear
  !! equation
  subroutine mus_load_turb_wallFunc(me, conf, parent)
    ! -------------------------------------------------------------------- !
    !> Turbulent wall model type to fill assign wall model
    type(mus_turb_wallFunc_type), intent(inout) :: me
    !> lua flu state
    type( flu_state ) :: conf
    !> bc parent handle
    integer, intent(in) :: parent
    ! -------------------------------------------------------------------- !
    integer :: iError
    ! -------------------------------------------------------------------- !
    write(logUnit(10), "(A)") '  Loading turbulent wall_function ...'
    me%isActive = .true.
    ! wall model function name
    call aot_get_val(L       = conf,            &
      &              thandle = parent,          &
      &              key     = 'wall_function', &
      &              val     = me%wall_func,    &
      &              default = 'musker',        &
      &              ErrCode = iError           )

    if (btest(iError, aoterr_Fatal)) then
      write(logUnit(1),*)'FATAL Error occured, while retrieving wall_function:'
      if (btest(iError, aoterr_WrongType)) then
         write(logUnit(1),*)'Variable has wrong type!'
         write(logUnit(1),*)'STOPPING'
         call tem_abort()
      end if
    end if

    write(logUnit(1),*) '  wall function: ', trim(me%wall_func)

    ! Get nonlinear solver for wall model other then power_law
    if (trim(me%wall_func) /= 'power_law') then
      call aot_get_val(L       = conf,                &
        &              thandle = parent,              &
        &              key     = 'nonlinear_solver',  &
        &              val     = me%nonlinear_solver, &
        &              default = 'fixed_point',       &
        &              ErrCode = iError               )

      if (btest(iError, aoterr_Fatal)) then
        write(logUnit(1),*)'FATAL Error occured, while retrieving ' &
          &              //'nonlinear_solver:'
        if (btest(iError, aoterr_WrongType)) then
           write(logUnit(1),*)'Variable has wrong type!'
           write(logUnit(1),*)'STOPPING'
           call tem_abort()
        end if
      end if
      write(logUnit(1),*) '  nonlinear solver: ', trim(me%nonlinear_solver)
    end if

    ! Von-Karman constant
    call aot_get_val(L       = conf,         &
      &              thandle = parent,       &
      &              key     = 'von_karman', &
      &              val     = me%vonKarman, &
      &              default = 0.4_rk,       &
      &              ErrCode = iError        )

    if (btest(iError, aoterr_Fatal)) then
      write(logUnit(1),*)'FATAL Error occured, while retrieving ' &
        &              //'von_karman:'
      if (btest(iError, aoterr_WrongType)) then
         write(logUnit(1),*)'Variable has wrong type!'
         write(logUnit(1),*)'STOPPING'
         call tem_abort()
      end if
    end if
    write(logUnit(1),*) '  von Karman: ', me%vonKarman

    ! Use van Driest damping function
    call aot_get_val(L       = conf,             &
      &              thandle = parent,           &
      &              key     = 'use_van_driest', &
      &              val     = me%useVanDriest,  &
      &              default = .true.,           &
      &              ErrCode = iError            )

    if (btest(iError, aoterr_Fatal)) then
      write(logUnit(1),*)'FATAL Error occured, while retrieving ' &
        &              //'use_van_driest:'
      if (btest(iError, aoterr_WrongType)) then
         write(logUnit(1),*)'Variable has wrong type!'
         write(logUnit(1),*)'STOPPING'
         call tem_abort()
      end if
    end if
    write(logUnit(1),*) '  use van Driest: ', me%useVanDriest

    ! Assign turbulent wall model function pointer for friction and streamwise
    ! velocity calculation
    call assign_turb_wallFunc(me)

  end subroutine mus_load_turb_wallFunc
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> Assign wall models pointers to calculate friction velocity and stream-wise
  !! velocity component
  subroutine assign_turb_wallFunc(me)
    ! -------------------------------------------------------------------- !
    !> Turbulent wall model type to fill assign wall model
    type(mus_turb_wallFunc_type), intent(inout) :: me
    ! -------------------------------------------------------------------- !
    ! -------------------------------------------------------------------- !
    select case (trim(me%wall_func))
    case ('musker')
      me%calcStreamWiseVel => streamWiseVel_musker
      select case (trim(me%nonlinear_solver))
      case ('fixed_point')
        me%calcFricVel => fricVel_musker_fixedPoint
      case ('newton')
        me%calcFricVel => fricVel_musker_newton
      case default
        call tem_abort('Error: Unsupported nonlinear solver for wall function')
      end select
    case ('reichardt')
      me%calcStreamWiseVel => streamWiseVel_reichardt
      select case (trim(me%nonlinear_solver))
      case ('fixed_point')
        me%calcFricVel => fricVel_reichardt_fixedPoint
      case ('newton')
        me%calcFricVel => fricVel_reichardt_newton
      case default
        call tem_abort('Error: Unsupported nonlinear solver for wall function')
      end select
    case ('power_law')
      me%calcStreamWiseVel => streamWiseVel_powerLaw
      me%calcFricVel => fricVel_powerLaw
    case default
      call tem_abort("Error: Unknown turbulent wall function")
    end select

  end subroutine assign_turb_wallFunc
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This routine computes friction velocity from Musker wall model profile
  !! using fixed-point iteration method
  subroutine fricVel_musker_fixedPoint(velTau, velSW, distToBnd, viscKine, &
    &                                  nElems)
    ! -------------------------------------------------------------------- !
    !> Friction velocity computed from wall model.
    !! it is inout to provide velTau from previous timestep as initial velTau
    !! for fixed-point or Newton iteration solver
    real(kind=rk), intent(inout) :: velTau(:)
    !> Stream-wise velocity component from which friction velocity is computed
    real(kind=rk), intent(in) :: velSW(:)
    !> Distance to the boundary in the discrete normai direction
    real(kind=rk), intent(in) :: distToBnd(:)
    !> Kinematic viscosity
    real(kind=rk), intent(in) :: viscKine(:)
    !> Number of elements in input and output arrays
    integer, intent(in) :: nElems
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: error
    integer :: iElem, iter
    real(kind=rk) :: velTau_old, velTau_new, yPlus, yPlus_fac
    ! -------------------------------------------------------------------- !
    do iElem = 1, nElems
      ! Fixed point iteration
      iter = 0
      velTau_old = velTau(iElem)
      yPlus_fac = distToBnd(iElem) / viscKine(iElem)
      fixed_point: do
        yPlus = velTau_old * yPlus_fac
        velTau_new = velSW(iElem) / ( 5.424_rk                    &
          &        * atan((2.0_rk*yPlus - 8.15_rk)/16.7_rk)       &
          &        + log10( (yPlus + 10.6_rk)**9.6_rk             &
          &            / ( yPlus**2.0_rk - 8.15_rk*yPlus          &
          &                + 86.0_rk)**2.0_rk ) - 3.5072790194_rk )
        error = abs(velTau_new - velTau_old)
        velTau_old = velTau_new
        iter = iter+1
        if ((error .fle. imEq_tol)    &
          & .or. (iter >= imEq_nIter)) exit fixed_point
      end do fixed_point
      velTau(iElem) = velTau_new
    end do

  end subroutine fricVel_musker_fixedPoint
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This routine computes friction velocity from Musker wall model profile
  !! using Newton iteration method
  subroutine fricVel_musker_newton(velTau, velSW, distToBnd, viscKine, &
    &                              nElems)
    ! -------------------------------------------------------------------- !
    !> Friction velocity computed from wall model.
    !! it is inout to provide velTau from previous timestep as initial velTau
    !! for fixed-point or Newton iteration solver
    real(kind=rk), intent(inout) :: velTau(:)
    !> Stream-wise velocity component from which friction velocity is computed
    real(kind=rk), intent(in) :: velSW(:)
    !> Distance to the boundary in the discrete normai direction
    real(kind=rk), intent(in) :: distToBnd(:)
    !> Kinematic viscosity
    real(kind=rk), intent(in) :: viscKine(:)
    !> Number of elements in input and output arrays
    integer, intent(in) :: nElems
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: error
    integer :: iElem, iter
    real(kind=rk) :: velTau_old, velTau_new, yPlus, yPlus_fac
    real(kind=rk) :: log_num, log_den, tan_num, fx, dfdx
    ! -------------------------------------------------------------------- !
    do iElem = 1, nElems
      ! Fixed point iteration
      iter = 0
      velTau_old = velTau(iElem)
      yPlus_fac = distToBnd(iElem) / viscKine(iElem)
      fixed_point: do
        yPlus = velTau_old * yPlus_fac
        velTau_new = velSW(iElem) / ( 5.424_rk                    &
          &        * atan((2.0_rk*yPlus - 8.15_rk)/16.7_rk)       &
          &        + log10( (yPlus + 10.6_rk)**9.6_rk             &
          &            / ( yPlus**2.0_rk - 8.15_rk*yPlus          &
          &                + 86.0_rk)**2.0_rk ) - 3.5072790194_rk )

        log_num = (yPlus + 10.6_rk)**9.6_rk
        log_den = (yPlus**2.0_rk - 8.15_rk * yPlus + 86.0_rk)**2.0_rk
        tan_num = 2.0_rk*yPlus - 8.15_rk
        ! Function of u_tau
        fx = velTau_old*(5.424_rk*atan(tan_num/16.7_rk) &
          & + log10(log_num/log_den) - 3.5072790194_rk)-velSW(iElem)

        ! derivative of function w.r.t u_tau
        dfdx = velTau_old * ( ( 9.6_rk*yPlus_fac*(yPlus+10.6_rk)**8.6_rk      &
          &  / log_den - (2_rk*log_num)*(2_rk*yPlus-8.15_rk*yPlus_Fac)        &
          &  / log_den**1.5_rk ) * ( log_den / (log(10.0_rk)*log_num) )       &
          & + 0.649581*yPlus_fac/(0.00358564_rk*(2.0_rk*yPlus-8.15)**2 + 1) ) &
          & + log(log_num/log_den)/log(10._rk)                                &
          & + 5.424_rk*atan(0.0598802_rk*(tan_num))-3.5072790194_rk

        velTau_new = velTau_old - (fx/dfdx)

        error = abs(velTau_new - velTau_old)
        velTau_old = velTau_new
        iter = iter+1
        if ((error .fle. imEq_tol)    &
          & .or. (iter >= imEq_nIter)) exit fixed_point
      end do fixed_point
      velTau(iElem) = velTau_new
    end do
  end subroutine fricVel_musker_newton
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This routines computes streamWise velocity component from friction velocity
  !! and distance to boundary using Musker profile.
  subroutine streamWiseVel_musker(velSW, velTau, distToBnd, &
    &                             viscKine, nElems)
    ! -------------------------------------------------------------------- !
    !> Stream-wise velocity component from wall model
    real(kind=rk), intent(out) :: velSW(:)
    !> Friction velocity computd from wall model
    real(kind=rk), intent(in) :: velTau(:)
    !> Distance to the boundary in the discrete normai direction
    real(kind=rk), intent(in) :: distToBnd(:)
    !> Kinematic viscosity
    real(kind=rk), intent(in) :: viscKine(:)
    !> Number of elements in input and output arrays
    integer, intent(in) :: nElems
    ! -------------------------------------------------------------------- !
    integer :: iElem
    real(kind=rk) :: yPlus
    ! -------------------------------------------------------------------- !
    do iElem = 1, nElems
      yPlus = distToBnd(iElem) * velTau(iElem) / viscKine(iElem)

      velSW(iElem) = velTau(iElem) * ( 5.424_rk                   &
        &          * atan((2.0_rk*yPlus - 8.15_rk)/16.7_rk)       &
        &          + log10( (yPlus + 10.6_rk)**9.6_rk             &
        &              / ( yPlus**2.0_rk - 8.15_rk*yPlus          &
        &                  + 86.0_rk)**2.0_rk ) - 3.5072790194_rk )
    end do
  end subroutine streamWiseVel_musker
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This routine computes friction velocity from reichardt wall model profile
  !! using fixed-point iteration method
  subroutine fricVel_reichardt_fixedPoint(velTau, velSW, distToBnd, viscKine, &
    &                                     nElems)
    ! -------------------------------------------------------------------- !
    !> Friction velocity computed from wall model.
    !! it is inout to provide velTau from previous timestep as initial velTau
    !! for fixed-point or Newton iteration solver
    real(kind=rk), intent(inout) :: velTau(:)
    !> Stream-wise velocity component from which friction velocity is computed
    real(kind=rk), intent(in) :: velSW(:)
    !> Distance to the boundary in the discrete normai direction
    real(kind=rk), intent(in) :: distToBnd(:)
    !> Kinematic viscosity
    real(kind=rk), intent(in) :: viscKine(:)
    !> Number of elements in input and output arrays
    integer, intent(in) :: nElems
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: error
    integer :: iElem, iter
    real(kind=rk) :: velTau_old, velTau_new, yPlus, yPlus_fac
    ! -------------------------------------------------------------------- !
    do iElem = 1, nElems
      ! Fixed point iteration
      iter = 0
      velTau_old = velTau(iElem)
      yPlus_fac = distToBnd(iElem) / viscKine(iElem)
      fixed_point: do
        yPlus = velTau_old * yPlus_fac
        velTau_new =  velSW(iElem) / ( oneOvervonKA              &
           &               * log( 1.0_rk + reivonKA*yPlus )      &
           &               + reiC * ( 1.0_rk-exp(-yPlus/reiB1)   &
           &               - (yPlus/reiB1) * exp(-yPlus/reiB2) ) )

        error = abs(velTau_new - velTau_old)
        velTau_old = velTau_new
        iter = iter+1
        if ((error .fle. imEq_tol)    &
          & .or. (iter >= imEq_nIter)) exit fixed_point
      end do fixed_point
      velTau(iElem) = velTau_new
    end do

  end subroutine fricVel_reichardt_fixedPoint
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This routine computes friction velocity from reichardt wall model profile
  !! using Newton iteration method
  subroutine fricVel_reichardt_newton(velTau, velSW, distToBnd, viscKine, &
    &                                 nElems)
    ! -------------------------------------------------------------------- !
    !> Friction velocity computed from wall model.
    !! it is inout to provide velTau from previous timestep as initial velTau
    !! for fixed-point or Newton iteration solver
    real(kind=rk), intent(inout) :: velTau(:)
    !> Stream-wise velocity component from which friction velocity is computed
    real(kind=rk), intent(in) :: velSW(:)
    !> Distance to the boundary in the discrete normai direction
    real(kind=rk), intent(in) :: distToBnd(:)
    !> Kinematic viscosity
    real(kind=rk), intent(in) :: viscKine(:)
    !> Number of elements in input and output arrays
    integer, intent(in) :: nElems
    ! -------------------------------------------------------------------- !
    ! -------------------------------------------------------------------- !
    call tem_abort("Error: fricVel_reichardt_newton not implemented yet")
  end subroutine fricVel_reichardt_newton
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This routines computes streamWise velocity component from friction velocity
  !! and distance to boundary using reichardt profile.
  subroutine streamWiseVel_reichardt(velSW, velTau, distToBnd, &
    &                                viscKine, nElems)
    ! -------------------------------------------------------------------- !
    !> Stream-wise velocity component from wall model
    real(kind=rk), intent(out) :: velSW(:)
    !> Friction velocity computd from wall model
    real(kind=rk), intent(in) :: velTau(:)
    !> Distance to the boundary in the discrete normai direction
    real(kind=rk), intent(in) :: distToBnd(:)
    !> Kinematic viscosity
    real(kind=rk), intent(in) :: viscKine(:)
    !> Number of elements in input and output arrays
    integer, intent(in) :: nElems
    ! -------------------------------------------------------------------- !
    integer :: iElem
    real(kind=rk) :: yPlus
    ! -------------------------------------------------------------------- !
    do iElem = 1, nElems
      yPlus = distToBnd(iElem) * velTau(iElem) / viscKine(iElem)
      velSW(iElem) = velTau(iElem) * ( oneOvervonKA      &
           &       * log( 1.0_rk + reivonKA*yPlus )      &
           &       + reiC * ( 1.0_rk-exp(-yPlus/reiB1)   &
           &       - (yPlus/reiB1) * exp(-yPlus/reiB2) ) )
    end do
  end subroutine streamWiseVel_reichardt
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This routine computes friction velocity from explicit two equation
  !! Werner and Wengle wall model.
  !! Haussmann, Marc; BARRETO, Alejandro CLARO; KOUYI, Gislain LIPEME;
  !! Rivière, Nicolas; Nirschl, Hermann; Krause, Mathias J. (2019):
  !! Large-eddy simulation coupled with wall models for turbulent channel
  !! flows at high Reynolds numbers with a lattice Boltzmann method
  !! — Application to Coriolis mass flowmeter. In Computers & Mathematics
  !! with Applications 78 (10), pp. 3285–3302. DOI: 10.1016/j.camwa.2019.04.033.
  !!
  !! The explicit power-law in terms of friction velocity given in Eq. 33 in
  !! S. Wilhelm, J. Jacob, and P. Sagaut, "An explicit power-law-based wall 
  !! model for lattice Boltzmann method–Reynolds-averaged numerical simulations 
  !! of the flow around airfoils", Physics of Fluids 30, 065111 (2018) 
  !! https://doi.org/10.1063/1.5031764
  !!
  !! The model constants are chosen according to Werner and Wengle:
  !! Wengle, H. and Werner, H. (1993) ‘Large-eddy Simulation of Turbulent Flow
  !! Over Sharp-edged Obstacles in a Plate Channel’, (1985), pp. 192–199.
  subroutine fricVel_powerLaw(velTau, velSW, distToBnd, viscKine, &
    &                             nElems)
    ! -------------------------------------------------------------------- !
    !> Friction velocity computed from wall model.
    !! it is inout to provide velTau from previous timestep as initial velTau
    !! for fixed-point or Newton iteration solver
    real(kind=rk), intent(inout) :: velTau(:)
    !> Stream-wise velocity component from which friction velocity is computed
    real(kind=rk), intent(in) :: velSW(:)
    !> Distance to the boundary in the discrete normai direction
    real(kind=rk), intent(in) :: distToBnd(:)
    !> Kinematic viscosity
    real(kind=rk), intent(in) :: viscKine(:)
    !> Number of elements in input and output arrays
    integer, intent(in) :: nElems
    ! -------------------------------------------------------------------- !
    integer :: iElem
    real(kind=rk) :: visc_div_dist, fac
    ! -------------------------------------------------------------------- !
    do iElem = 1, nElems
      visc_div_dist = viscKine(iElem) / distToBnd(iElem)
      if (velSW(iElem) <= ww_uLmt_vel * visc_div_dist) then
        ! y_plus below 11.81
        velTau(iElem) = sqrt(velSW(iElem) * visc_div_dist)
      else
        ! power-law
        fac = (visc_div_dist)**ww_mDivOnePlusM
        velTau(iElem) = fac * (velSW(iElem) / ww_C_m)**ww_onePlusMInv
      end if
    end do

  end subroutine fricVel_powerLaw
  ! -------------------------------------------------------------------------- !

  ! -------------------------------------------------------------------------- !
  !> This routines computes streamWise velocity component from friction velocity
  !! and distance to boundary using two equation Werner and Wengle profile.
  !! Ref to following paper for Schmitt three layer equations.
  !! Haussmann, M. et al. (2019) ‘Large-eddy simulation coupled with wall models
  !! for turbulent channel flows at high Reynolds numbers with a lattice
  !! Boltzmann method — Application to Coriolis mass flowmeter’, Computers &
  !! Mathematics with Applications. Elsevier Ltd, 78(10), pp. 3285–3302.
  subroutine streamWiseVel_powerLaw(velSW, velTau, distToBnd, &
    &                                   viscKine, nElems)
    ! -------------------------------------------------------------------- !
    !> Stream-wise velocity component from wall model
    real(kind=rk), intent(out) :: velSW(:)
    !> Friction velocity computd from wall model
    real(kind=rk), intent(in) :: velTau(:)
    !> Distance to the boundary in the discrete normai direction
    real(kind=rk), intent(in) :: distToBnd(:)
    !> Kinematic viscosity
    real(kind=rk), intent(in) :: viscKine(:)
    !> Number of elements in input and output arrays
    integer, intent(in) :: nElems
    ! -------------------------------------------------------------------- !
    integer :: iElem
    real(kind=rk) :: yPlus
    ! -------------------------------------------------------------------- !
    do iElem = 1, nElems
      yPlus = distToBnd(iElem) * velTau(iElem) / viscKine(iElem)
      if (yPlus >= sc_uLmt) then
        ! Interial layer, use powerlaw profile from werner and wengle
        velSW(iElem) = (ww_C_m * ( yPlus )**ww_m) * velTau(iElem)
      else if ( yPlus >= sc_lLmt .and. yPlus < sc_uLmt) then
        ! Buffer layer, use logarithmic profile
        velSW(iElem) = (sc_a*log(yPlus) + sc_b) * velTau(iElem)
      else
        ! viscous sublayer, use linear profile.
        velSW(iElem) = yPlus * velTau(iElem)
      end if
    end do
  end subroutine streamWiseVel_powerLaw
  ! -------------------------------------------------------------------------- !

end module mus_turb_wallFunc_module
