module obl_k_epsilon
    !< @brief k-epsilon scheme module
    ! --------------------------------------------------------------------------------
    
    ! modules used
    ! --------------------------------------------------------------------------------
#ifdef USE_CONFIG_PARSER
    use iso_c_binding, only: C_NULL_CHAR
    use config_parser
#endif

    use obl_math
    use obl_grid
    use obl_bc
    use obl_turb_common

    ! directives list
    ! --------------------------------------------------------------------------------
    implicit none
    private

    ! public interface
    ! --------------------------------------------------------------------------------
    public :: init_turbulence_closure
    public :: advance_turbulence_closure
    public :: define_stability_functions

    public :: get_eddy_viscosity, get_eddy_diffusivity
    public :: get_eddy_viscosity_vec, get_eddy_diffusivity_vec

    public :: TKE_init, EPS_init
    public :: TKE_bc, EPS_bc

    public :: set_config_param
    ! --------------------------------------------------------------------------------

    !> @brief k-epsilon parameters
    type, public :: kepsilonParamType
        real :: C_mu = 0.09                     !< momentum stability function constant
        real :: C1_eps = 1.44                   !< eps-eq.: production constant
        real :: C2_eps = 1.92                   !< eps-eq.: dissipation constant
        real :: C3_eps = -0.40                  !< eps-eq.: buoyancy constant (-0.4 stable stratification or 1.14 unstable stratification)
        real :: ScT_TKE = 1.0                   !< Schmidt number for TKE
        real :: ScT_EPS = 1.11                  !< Scmidt number for eps

        real :: TKE_min = 0.000001              !< minimal valur for TKE, [J/kg]
        real :: EPS_min = (0.000001)**(3.0/2.0) !< minimal value for epsilon (TKE dissipation rate), [W/kg]

        real :: kappa = 0.4                     !< von Karman constant
        real :: PrT0 = 1.25                     !< turbulent Prandtl number

        real :: c_S2_min = 1e-5     !< min shear frequency, [(1/s)**2]
    end type

    contains


    ! --------------------------------------------------------------------------------
    subroutine define_stability_functions(param, bc, grid)
        !< @brief advance stability functions: N**2, S**2, Ri-gradient
        ! ----------------------------------------------------------------------------
        use obl_state

        type(kepsilonParamType), intent(in) :: param
        type(oblBcType), intent(in) :: bc
        type (gridDataType), intent(in) :: grid
        ! --------------------------------------------------------------------------------

        call get_N2(N2, Rho, bc%Rho_dyn0, bc%Rho_dynH, &
            param%kappa, param%PrT0, grid)
        call get_S2(S2, U, V, bc%U_dyn0, bc%U_dynH, &
            param%kappa, grid)
        call get_Ri_gradient(Ri_grad, N2, S2, param%c_S2_min, grid)
    end subroutine

    ! --------------------------------------------------------------------------------
    subroutine init_turbulence_closure(param, bc, grid)
        !< @brief advance turbulence closure
        ! ----------------------------------------------------------------------------
        use obl_state

        type(kepsilonParamType), intent(in) :: param
        type(oblBcType), intent(in) :: bc
        type (gridDataType), intent(in) :: grid
        ! --------------------------------------------------------------------------------

        call TKE_init(TKE, param, grid%cz)
        call EPS_init(EPS, param, grid%cz, grid%height)

        call TKE_bc(TKE, &
            bc%U_dyn0, bc%U_dynH, param, grid%cz)
        call EPS_bc(EPS, &
            bc%U_dyn0, bc%U_dynH, param, grid%dz, grid%cz)

        call get_eddy_viscosity(Km, TKE, EPS, param, grid)
        call get_eddy_diffusivity(Kh, Km, param, grid)

    end subroutine

    ! --------------------------------------------------------------------------------
    subroutine advance_turbulence_closure(param, bc, grid, dt)
        !< @brief advance turbulence closure
        ! ----------------------------------------------------------------------------
        use obl_state

        type(kepsilonParamType), intent(in) :: param
        type(oblBcType), intent(in) :: bc
        type (gridDataType), intent(in) :: grid
        real, intent(in) :: dt
        ! --------------------------------------------------------------------------------

        call get_TKE_production(TKE_production, Km, S2, grid)
        call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid)

        call TKE_bc(TKE, &
            bc%U_dyn0, bc%U_dynH, param, grid%cz)
        call EPS_bc(EPS, &
            bc%U_dyn0, bc%U_dynH, param, grid%dz, grid%cz)
        call get_eddy_viscosity(Km, TKE, EPS, param, grid)
        call get_eddy_diffusivity(Kh, Km, param, grid)
        call TKE_bc(TKE, &
            bc%U_dyn0, bc%U_dynH, param, grid%cz)
        call eps_bc(EPS, &
            bc%U_dyn0, bc%U_dynH, param, grid%dz, grid%cz)
        !call solve_TKE_eq_semiimplicit (TKE, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, eps, param)
        call solve_TKE_eq (TKE, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, EPS, param)
        call limit_min_array(TKE, param%TKE_min, grid%cz)
        call solve_eps_eq_semiimplicit (EPS, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, TKE, param)
        call limit_min_array(EPS, param%eps_min, grid%cz) 

    end subroutine

    subroutine TKE_init(TKE, param, nz)
        !< @brief TKE initialization
        type(kepsilonParamType), intent(in) :: param
        integer, intent(in) :: nz !< number of z-steps
        real, intent(out) :: TKE(nz) !< TKE, [J/kg]
        integer :: k !< counter

        do k = 1, nz
            TKE(k) = 1.0 / (param%C_mu)**(1.0/2.0) * 0.001 * 0.001
        end do     
    end subroutine

    subroutine EPS_init(eps, param, nz, z)
        !< @brief TKE dissipation rate initialization
        type(kepsilonParamType), intent(in) :: param
        integer, intent(in) :: nz  !< number of z-steps
        real, intent(in) :: z !< depth, [m]
        real, intent(out) :: eps(nz) !< TKE dissipation rate, [W/kg]
        integer :: k !< counter

        do k = 1, nz
            eps(k) = 0.001 * 0.001 * 0.001 / (param%kappa * z) !not full z, but current z - must be fixed
        end do     
    end subroutine

   
    ! ----------------------------------------------------------------------------
    subroutine TKE_bc(TKE, U_dyn0, U_dynH, param, cz) 
        !< @brief TKE boundary conditions
        ! ----------------------------------------------------------------------------
        type(kepsilonParamType), intent(in) :: param
        real, intent(in) :: U_dyn0, U_dynH      !< dynamic velocity, [(m/s)]
        integer, intent(in) :: cz               !< grid size
        
        real, intent(out) :: TKE(cz)            !< TKE, [(m/s)**2]
        ! ----------------------------------------------------------------------------

        TKE(1) = (1.0 / sqrt(param%C_mu)) * U_dyn0 * U_dyn0
        TKE(cz) = (1.0 / sqrt(param%C_mu)) * U_dynH * U_dynH

    end subroutine TKE_bc

    ! ----------------------------------------------------------------------------
    subroutine EPS_bc(EPS, U_dyn0, U_dynH, param, dz, cz)
        !< @brief epsilon boundary conditions
        ! ----------------------------------------------------------------------------
        type(kepsilonParamType), intent(in) :: param
        real, intent(in) :: U_dyn0, U_dynH      !< dynamic velocity, [(m/s)]
        integer, intent(in) :: cz               !< grid size
        real, intent(in) :: dz(cz)              !< grid steps

        real, intent(out) :: EPS(cz)            !< dissipation rate, [m^2/s^3]

        EPS(1) = (U_dyn0 * U_dyn0 * U_dyn0) / (param%kappa * 0.5 * dz(1))
        EPS(cz) = (U_dynH * U_dynH * U_dynH) / (param%kappa * 0.5 * dz(cz)) 

    end subroutine EPS_bc

    ! --------------------------------------------------------------------------------
    subroutine get_eddy_viscosity(Km, TKE, EPS, param, grid)
        !< @brief calculate eddy viscosity
        ! ----------------------------------------------------------------------------
        type(kepsilonParamType), intent(in) :: param
        type (gridDataType), intent(in) :: grid
        
        real, dimension(grid%cz), intent(out) :: Km          !< eddy viscosity, [m**2 / s] 
        real, dimension(grid%cz), intent(in) :: TKE, EPS     !<  
        
        ! ----------------------------------------------------------------------------

        call get_eddy_viscosity_vec(Km, TKE, EPS, param, grid%cz)
    end subroutine get_eddy_viscosity 

    subroutine get_eddy_viscosity_vec(Km, TKE, EPS, param, cz)
        !< @brief calculate eddy viscosity
        ! ----------------------------------------------------------------------------
        type(kepsilonParamType), intent(in) :: param
        
        integer, intent(in) :: cz                       !< grid size
        real, dimension(cz), intent(out) :: Km          !< eddy viscosity, [m**2 / s] 
        real, dimension(cz), intent(in) :: TKE, EPS     !<  
        
        integer :: k
        ! ----------------------------------------------------------------------------

        do k = 1, cz
            if (eps(k) < param%eps_min) then
                Km(k) = param%C_mu * TKE(k) * TKE(k) / param%eps_min
            else
                Km(k) = param%C_mu * TKE(k) * TKE(k) / eps(k)
            endif
        end do
    end subroutine get_eddy_viscosity_vec 

    ! --------------------------------------------------------------------------------
    subroutine get_eddy_diffusivity(Kh, Km, param, grid)
        !< @brief calculate eddy diffusivity
        ! ----------------------------------------------------------------------------
        type(kepsilonParamType), intent(in) :: param
        type (gridDataType), intent(in) :: grid
         
        real, dimension(grid%cz), intent(out) :: Kh      !< eddy diffusivity, [m**2 / s]
        real, dimension(grid%cz), intent(in) :: Km       !< eddy viscosity, [m**2 / s]
        ! ----------------------------------------------------------------------------

        call get_eddy_diffusivity_vec(Kh, Km, param, grid%cz)
    end subroutine get_eddy_diffusivity

    subroutine get_eddy_diffusivity_vec(Kh, Km, param, cz)
        !< @brief calculate eddy diffusivity
        ! ----------------------------------------------------------------------------
        type(kepsilonParamType), intent(in) :: param
        
        integer, intent(in) :: cz                   !< grid size 
        real, dimension(cz), intent(out) :: Kh      !< eddy diffusivity, [m**2 / s]
        real, dimension(cz), intent(in) :: Km       !< eddy viscosity, [m**2 / s]
        
        integer :: k
        ! ----------------------------------------------------------------------------

        do k = 1, cz
            Kh(k) = Km(k) / param%PrT0
        end do
    end subroutine get_eddy_diffusivity_vec

    ! --------------------------------------------------------------------------------
    subroutine solve_TKE_eq (Field, Km, nz, dz, dt, P_TKE, B_TKE, eps, param)
        !< @brief solver for TKE equation (explicit)

        type(kepsilonParamType), intent(in) :: param

        integer, intent (in) :: nz !< number of steps in z (depth), [m]
	    real, intent (in) :: dt !< time step [s]
        real, intent (in) :: dz(nz) !< depth(z) step [m] 
        real, intent(in) :: Km(nz)

        real :: Kvisc(nz) !< turbulent transfer coefficient [m^2/s]
        real :: G(nz) !< explicit part  
        real :: a(nz), b(nz), c(nz), f(nz) !< coefficients for tridiagonal matrix
        real :: Field(nz) !< value of unknown field in the point
        integer :: k !< counter   
        real :: Kvisc_plus(nz), Kvisc_minus(nz) !<half-sums for turbulent transfer coefficient Kvisc
        real :: dz_plus(nz), dz_minus(nz) !<half-sums for depth-step dz
        real :: f_right(nz) !< right-hand side
        real :: gamma  !< variable determining the explicitness/implicitness of a scheme (from 0.0 to 1.0)
        real :: P_TKE(nz) !< shear production of TKE, [m**2 / s**3]
        real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
        real :: eps(nz) !< TKE dissipation rate, [W/kg]    

        Kvisc(:) = Km(:) / param%ScT_TKE

            gamma = 0.0

            do k = 1, nz
                f_right(k) = P_TKE(k) + B_TKE(k) - eps(k) 
            end do  
            
            do k = 1, nz-1
                Kvisc_plus(k) = (Kvisc(k)+Kvisc(k+1))/2.0   
                dz_plus(k) =  (dz(k)+dz(k+1))/2.0   
            end do    
    
            do k = 2, nz
                Kvisc_minus(k) = (Kvisc(k)+Kvisc(k-1))/2.0  
                dz_minus(k) = (dz(k)+dz(k-1))/2.0 
            end do  
    
            do k = 2, nz-1
                G(k) = Kvisc_plus(k)/dz_plus(k) * (Field(k+1) - Field(k)) &
                - Kvisc_minus(k)/dz_minus(k) * (Field(k) - Field(k-1))
            end do  
    
                G(1) = Kvisc_plus(1)/dz_plus(1) * (Field(2) - Field(1))
                G(nz) = - Kvisc_minus(nz)/dz_minus(nz) * (Field(nz) - Field(nz-1))
                       
            do k = 2, nz-1
                a(k) = dt / dz(k) * Kvisc_minus(k)/dz_minus(k) * (1.0 - gamma)
                b(k) = -1.0 - dt / dz(k) * Kvisc_plus(k)/dz_plus(k)  * (1.0 - gamma) &
                    -  dt / dz(k) * Kvisc_minus(k)/dz_minus(k)  * (1.0 - gamma)     
                c(k) = dt / dz(k) * Kvisc_plus(k)/dz_plus(k) * (1.0 - gamma)
                f(k) = - Field(k) - f_right(k) * dt  &
                    - dt / dz(k) * G(k) * gamma

                if (k == 2) then
                    f(k) = f(k) - a(k) * Field(k-1)
                endif
                if (k == nz-1) then
                    f(k) = f(k) - c(k) * Field(k+1)
                endif
                if (k == 2) a(k) = 0.0
                if (k == nz-1) c(k) = 0.0
            
            end do  
    
            call tma (Field(2:nz-1), a(2:nz-1), b(2:nz-1), c(2:nz-1), f(2:nz-1), nz - 2) 					
        
    end subroutine solve_TKE_eq


    subroutine solve_eps_eq (Field, Km, nz, dz, dt, P_TKE, B_TKE, TKE, param)
        !< @brief solver for epsilon equation (explicit)

        type(kepsilonParamType), intent(in) :: param

        integer, intent (in) :: nz !< number of steps in z (depth), [m]
	    real, intent (in) :: dt !< time step [s]
        real, intent (in) :: dz(nz) !< depth(z) step [m] 
        real, intent(in) :: Km(nz)

        real :: Kvisc(nz) !< turbulent transfer coefficient [m^2/s]
        real :: G(nz) !< explicit part  
        real :: a(nz), b(nz), c(nz), f(nz) !< coefficients for tridiagonal matrix
        real :: Field(nz) !< value of unknown field in the point
        integer :: k !< counter   
        real :: Kvisc_plus(nz), Kvisc_minus(nz) !<half-sums for turbulent transfer coefficient Kvisc
        real :: dz_plus(nz), dz_minus(nz) !<half-sums for depth-step dz
        real :: f_right(nz) !< right-hand side
        real :: gamma  !< variable determining the explicitness/implicitness of a scheme (from 0.0 to 1.0)
        real :: P_TKE(nz) !< shear production of TKE, [m**2 / s**3]
        real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
        real :: eps(nz) !< TKE dissipation rate, [W/kg]    
        real :: TKE(nz) !< TKE, [J/kg]

        Kvisc(:) = Km(:) / param%ScT_EPS
                 
            gamma = 0.0

            do k = 1, nz
                f_right(k) = Field(k) / TKE(k) * (param%C1_eps * P_TKE(k) + param%C3_eps * B_TKE(k) - param%C2_eps * Field(k))  
            end do  
            
            do k = 1, nz-1
                Kvisc_plus(k) = (Kvisc(k)+Kvisc(k+1))/2.0   
                dz_plus(k) =  (dz(k)+dz(k+1))/2.0   
            end do    
    
            do k = 2, nz
                Kvisc_minus(k) = (Kvisc(k)+Kvisc(k-1))/2.0 
                dz_minus(k) = (dz(k)+dz(k-1))/2.0 
            end do  
    
            do k = 2, nz-1
                G(k) = Kvisc_plus(k)/dz_plus(k) * (Field(k+1) - Field(k)) &
                    - Kvisc_minus(K)/dz_minus(k) * (Field(k) - Field(k-1))
            end do  
    
                G(1) = Kvisc_plus(1)/dz_plus(1) * (Field(2) - Field(1))
                G(nz) = - Kvisc_minus(nz)/dz_minus(nz) * (Field(nz) - Field(nz-1))
                       
            do k = 2, nz-1
                a(k) = dt / dz(k) * Kvisc_minus(k)/dz_minus(k) * (1.0 - gamma)
                b(k) = -1.0 - dt / dz(k) * Kvisc_plus(k)/dz_plus(k)  * (1.0 - gamma) &
                    -  dt / dz(k) * Kvisc_minus(k)/dz_minus(k)  * (1.0 - gamma)     
                c(k) = dt / dz(k) * Kvisc_plus(k)/dz_plus(k) * (1.0 - gamma)
                f(k) = - Field(k) - f_right(k) * dt &
                        -  dt / dz(k) * G(k) * gamma

                if (k == 2) then
                    f(k) = f(k) - a(k) * Field(k-1)
                endif
                if (k == nz-1) then
                    f(k) = f(k) - c(k) * Field(k+1)
                endif
                if (k == 2) a(k) = 0.0
                if (k == nz-1) c(k) = 0.0
                
            end do  
    
            call tma (Field(2:nz-1), a(2:nz-1), b(2:nz-1), c(2:nz-1), f(2:nz-1), nz - 2) 	    				
            
        end subroutine solve_eps_eq

        subroutine solve_TKE_eq_semiimplicit (Field, Km, nz, dz, dt, P_TKE, B_TKE, eps, param)
            !< @brief solver for TKE equation (semimplicit)

            type(kepsilonParamType), intent(in) :: param

            integer, intent (in) :: nz !< number of steps in z (depth), [m]
	        real, intent (in) :: dt !< time step [s]
            real, intent (in) :: dz(nz) !< depth(z) step [m] 
            real, intent(in) :: Km(nz)

            real :: Kvisc(nz) !< turbulent transfer coefficient [m^2/s]
            real :: G(nz) !< explicit part  
            real :: a(nz), b(nz), c(nz), f(nz) !< coefficients for tridiagonal matrix
            real :: Field(nz) !< value of unknown field in the point
            integer :: k !< counter   
            real :: Kvisc_plus(nz), Kvisc_minus(nz) !<half-sums for turbulent transfer coefficient Kvisc
            real :: dz_plus(nz), dz_minus(nz) !<half-sums for depth-step dz
            real :: f_right(nz) !< right-hand side
            real :: gamma  !< variable determining the explicitness/implicitness of a scheme (from 0.0 to 1.0)
            real :: P_TKE(nz) !< shear production of TKE, [m**2 / s**3]
            real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
            real :: eps(nz) !< TKE dissipation rate, [W/kg]    
            
            Kvisc(:) = Km(:) / param%ScT_TKE

                gamma = 0.0

                do k = 1, nz
                    f_right(k) = P_TKE(k) + B_TKE(k)
                end do  
            
                do k = 1, nz-1
                    Kvisc_plus(k) = (Kvisc(k)+Kvisc(k+1))/2.0   
                    dz_plus(k) =  (dz(k)+dz(k+1))/2.0   
                end do    
    
                do k = 2, nz
                    Kvisc_minus(k) = (Kvisc(k)+Kvisc(k-1))/2.0  
                    dz_minus(k) = (dz(k)+dz(k-1))/2.0 
                end do  
    
                do k = 2, nz-1
                    G(k) = Kvisc_plus(k)/dz_plus(k) * (Field(k+1) - Field(k)) &
                        - Kvisc_minus(k)/dz_minus(k) * (Field(k) - Field(k-1))
                end do  
    
                    G(1) = Kvisc_plus(1)/dz_plus(1) * (Field(2) - Field(1))
                    G(nz) = - Kvisc_minus(nz)/dz_minus(nz) * (Field(nz) - Field(nz-1))
                       
                do k = 2, nz-1
                    a(k) = dt / dz(k) * Kvisc_minus(k)/dz_minus(k) * (1.0 - gamma)
                    b(k) = -1.0 - dt * eps(k)/Field(k) - dt / dz(k) * &
                            Kvisc_plus(k)/dz_plus(k) * (1.0 - gamma) &
                            -  dt / dz(k) * Kvisc_minus(k)/dz_minus(k) &
                             * (1.0 - gamma)         
                    c(k) = dt / dz(k) * Kvisc_plus(k)/dz_plus(k) * (1.0 - gamma)
                    f(k) = - Field(k) - f_right(k) * dt  &
                            - dt / dz(k) * G(k) * gamma

                    if (k == 2) then
                        f(k) = f(k) - a(k) * Field(k-1)
                    endif
                    if (k == nz-1) then
                        f(k) = f(k) - c(k) * Field(k+1)
                    endif
                    if (k == 2) a(k) = 0.0
                    if (k == nz-1) c(k) = 0.0
            
                end do  
    
                call tma (Field(2:nz-1), a(2:nz-1), b(2:nz-1), c(2:nz-1), f(2:nz-1), nz - 2) 					
            
        end subroutine solve_TKE_eq_semiimplicit

        subroutine solve_eps_eq_semiimplicit (Field, Km, nz, dz, dt, P_TKE, B_TKE, TKE, param)
	        !< @brief solver for epsilon equation (semiimplict)
            type(kepsilonParamType), intent(in) :: param

            integer, intent (in) :: nz !< number of steps in z (depth), [m]
            real, intent (in) :: dt !< time step [s]
            real, intent (in) :: dz(nz) !< depth(z) step [m] 
            real, intent(in) :: Km(nz)

            real :: Kvisc(nz) !< turbulent transfer coefficient [m^2/s]
            real :: G(nz) !< explicit part  
            real :: a(nz), b(nz), c(nz), f(nz) !< coefficients for tridiagonal matrix
            real :: Field(nz) !< value of unknown field in the point
            integer :: k !< counter   
            real :: Kvisc_plus(nz), Kvisc_minus(nz) !<half-sums for turbulent transfer coefficient Kvisc
            real :: dz_plus(nz), dz_minus(nz) !<half-sums for depth-step dz
            real :: f_right(nz) !< right-hand side
            real :: gamma  !< variable determining the explicitness/implicitness of a scheme (from 0.0 to 1.0)
            real :: P_TKE(nz) !< shear production of TKE, [m**2 / s**3]
            real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
            real :: eps(nz) !< TKE dissipation rate, [W/kg]    
            real :: TKE(nz) !< TKE, [J/kg]

        
            Kvisc(:) = Km(:) / param%ScT_EPS

                gamma = 0.0

                do k = 1, nz
                    f_right(k) = Field(k) / TKE(k) * (param%C1_eps * P_TKE(k) + param%C3_eps * B_TKE(k))  
                end do  
            
                do k = 1, nz-1
                    Kvisc_plus(k) = (Kvisc(k)+Kvisc(k+1))/2.0   
                    dz_plus(k) =  (dz(k)+dz(k+1))/2.0   
                end do    
            
                do k = 2, nz
                    Kvisc_minus(k) = (Kvisc(k)+Kvisc(k-1))/2.0 
                    dz_minus(k) = (dz(k)+dz(k-1))/2.0
                end do  
    
                do k = 2, nz-1
                    G(k) = Kvisc_plus(k)/dz_plus(k) * (Field(k+1) - Field(k)) &
                        - Kvisc_minus(k)/dz_minus(k) * (Field(k) - Field(k-1))
                end do  
    
                    G(1) = Kvisc_plus(1)/dz_plus(1) * (Field(2) - Field(1))
                    G(nz) = - Kvisc_minus(nz)/dz_minus(nz) * (Field(nz) - Field(nz-1))
                           
                do k =  2, nz-1
                    a(k) = dt / dz(k) * Kvisc_minus(k)/dz_minus(k) * (1.0 - gamma)
                    b(k) = -1.0 - param%C2_eps * dt * Field(k)/TKE(k) - &
                        dt / dz(k) * Kvisc_plus(k)/dz_plus(k)  * (1.0 - gamma) &
                        -  dt / dz(k) * Kvisc_minus(k)/dz_minus(k)  * (1.0 - gamma)    
                    c(k) = dt / dz(k) * Kvisc_plus(k)/dz_plus(k) * (1.0 - gamma)
                    f(k) = - Field(k) - f_right(k) * dt &
                        -  dt / dz(k) * G(k) * gamma

                    if (k == 2) then
                        f(k) = f(k) - a(k) * Field(k-1)
                    endif
                    if (k == nz-1) then
                        f(k) = f(k) - c(k) * Field(k+1)
                    endif
                    if (k == 2) a(k) = 0.0
                    if (k == nz-1) c(k) = 0.0
            end do  

            call tma (Field(2:nz-1), a(2:nz-1), b(2:nz-1), c(2:nz-1), f(2:nz-1), nz - 2) 	    				
            
        end subroutine solve_eps_eq_semiimplicit

    ! --------------------------------------------------------------------------------
    subroutine set_config_param(param, tag, ierr)
        !< @brief set turbulence closure parameters
        ! ----------------------------------------------------------------------------
        type(kepsilonParamType), intent(out) :: param
        integer, intent(out) :: ierr

        character(len = *), intent(in) :: tag

        integer :: status
        ! --------------------------------------------------------------------------------

        ierr = 0        ! = OK

#ifdef USE_CONFIG_PARSER
        call c_config_is_varname(trim(tag)//".C_mu"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".C_mu"//C_NULL_CHAR, param%C_mu, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".C1_eps"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".C1_eps"//C_NULL_CHAR, param%C1_eps, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".C2_eps"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".C2_eps"//C_NULL_CHAR, param%C2_eps, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".C3_eps"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".C3_eps"//C_NULL_CHAR, param%C3_eps, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".ScT_TKE"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".ScT_TKE"//C_NULL_CHAR, param%ScT_TKE, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".ScT_EPS"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".ScT_EPS"//C_NULL_CHAR, param%ScT_EPS, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".TKE_min"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".TKE_min"//C_NULL_CHAR, param%TKE_min, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".EPS_min"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".EPS_min"//C_NULL_CHAR, param%EPS_min, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".kappa"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".kappa"//C_NULL_CHAR, param%kappa, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".PrT0"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".PrT0"//C_NULL_CHAR, param%PrT0, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".c_S2_min"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".c_S2_min"//C_NULL_CHAR, param%c_S2_min, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
#else
        !> *: just skipping setup without configuration file
#endif

    end subroutine

end module