Skip to content
Snippets Groups Projects
Commit 95c99b68 authored by Evgeny Mortikov's avatar Evgeny Mortikov
Browse files

k-epsilon module update

parent 6b751af8
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
obl_k_epsilon.f90
+ 53
− 37
View file @ 95c99b68
......@@ -9,104 +9,115 @@ module obl_k_epsilon
! directives list
implicit none
real, parameter :: C_mu = 0.09 !< momentum stability function constant
real, parameter :: PrT = 1.25 !< turbulent Prandtl number
real, parameter :: TKE_min = 0.000001 !< !< minimal valur for TKE, [J/kg]
real, parameter :: eps_min = TKE_min**(3.0/2.0) !0.0000001 !TKE_min**(3.0/2.0)!< minimal value for epsilon (TKE dissipation rate), [W/kg]
real, parameter :: kappa_k_eps = 0.4 !< von Karman constant
real, parameter :: C1_eps = 1.44 !< eps-eq.: production constant
real, parameter :: C2_eps = 1.92 !< eps-eq.: dissipation constant
real, parameter :: C3_eps = -0.40 !< eps-eq.: buoyancy constant (-0.4 stable stratification or 1.14 unstable stratification)
real, parameter :: Sch_TKE = 1.0 !< Schmidt number for TKE
real, parameter :: Sch_eps = 1.11 !1.11 !< Scmidt number for eps
!> @brief k-epsilon parameters
type, public :: kepsilonParamType
real :: C_mu = 0.09 !< momentum stability function constant
real :: PrT = 1.25 !< turbulent Prandtl number
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 :: 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 :: Sch_TKE = 1.0 !< Schmidt number for TKE
real :: Sch_eps = 1.11 !< Scmidt number for eps
end type
contains
subroutine TKE_init (TKE, nz)
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 = 2, nz-1
TKE(k) = 1.0 / (C_mu)**(1.0/2.0) * 0.001 * 0.001
TKE(k) = 1.0 / (param%C_mu)**(1.0/2.0) * 0.001 * 0.001
end do
end subroutine
subroutine eps_init (eps, nz, z)
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 = 2, nz-1
eps(k) = 0.001 * 0.001 * 0.001 / (kappa_k_eps * z) !not full z, but current z - must be fixed
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, cz)
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(C_mu)) * U_dyn0 * U_dyn0
TKE(cz) = (1.0 / sqrt(C_mu)) * U_dynH * U_dynH
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, dz, cz)
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) / (kappa_k_eps * 0.5 * dz(1))
EPS(cz) = (U_dynH * U_dynH * U_dynH) / (kappa_k_eps * 0.5 * dz(cz))
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_Km_k_eps (Km, TKE, eps, nz)
subroutine get_Km_k_eps (Km, TKE, eps, param, nz)
!< @brief calculate momentum transfer coefficient in k-epsilon closure
type(kepsilonParamType), intent(in) :: param
integer, intent(in) :: nz !< number of z-steps
real, intent(inout) :: TKE(nz), eps(nz)
real, intent(out) :: Km(nz) !eddy viscosity for momentum, [m**2 / s]
integer :: k !< counter
do k = 1, nz
if (eps(k) < eps_min) then
Km(k) = C_mu * TKE(k) * TKE(k) / eps_min
if (eps(k) < param%eps_min) then
Km(k) = param%C_mu * TKE(k) * TKE(k) / param%eps_min
else
Km(k) = C_mu * TKE(k) * TKE(k) / eps(k)
Km(k) = param%C_mu * TKE(k) * TKE(k) / eps(k)
endif
end do
end subroutine get_Km_k_eps
subroutine get_Kh_k_eps (Kh, Km, nz)
subroutine get_Kh_k_eps (Kh, Km, param, nz)
!< @brief calculate scalar transfer coefficient in k-epsilon closure
type(kepsilonParamType), intent(in) :: param
integer, intent(in) :: nz !< number of z-steps
real, intent(in) :: Km(nz) !eddy viscosity for momentum, [m**2 / s]
real, intent(out) :: Kh(nz) !eddy diffusivity for tracers, [m**2 / s]
integer :: k !< counter
do k = 1, nz
Kh(k) = Km(k) / PrT
Kh(k) = Km(k) / param%PrT
end do
end subroutine get_Kh_k_eps
subroutine solve_TKE_eq (Field, Km, nz, dz, dt, P_TKE, B_TKE, eps)
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]
......@@ -125,7 +136,7 @@ module obl_k_epsilon
real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
real :: eps(nz) !< TKE dissipation rate, [W/kg]
Kvisc(:) = Km(:) / Sch_TKE
Kvisc(:) = Km(:) / param%Sch_TKE
gamma = 0.0
......@@ -175,9 +186,11 @@ module obl_k_epsilon
end subroutine solve_TKE_eq
subroutine solve_eps_eq (Field, Km, nz, dz, dt, P_TKE, B_TKE, TKE)
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]
......@@ -197,12 +210,12 @@ module obl_k_epsilon
real :: eps(nz) !< TKE dissipation rate, [W/kg]
real :: TKE(nz) !< TKE, [J/kg]
Kvisc(:) = Km(:) / Sch_eps
Kvisc(:) = Km(:) / param%Sch_eps
gamma = 0.0
do k = 1, nz
f_right(k) = Field(k) / TKE(k) * (C1_eps * P_TKE(k) + C3_eps * B_TKE(k) - C2_eps * Field(k))
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
......@@ -246,9 +259,11 @@ module obl_k_epsilon
end subroutine solve_eps_eq
subroutine solve_TKE_eq_semiimplicit (Field, Km, nz, dz, dt, P_TKE, B_TKE, eps)
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]
......@@ -267,7 +282,7 @@ module obl_k_epsilon
real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
real :: eps(nz) !< TKE dissipation rate, [W/kg]
Kvisc(:) = Km(:) / Sch_TKE
Kvisc(:) = Km(:) / param%Sch_TKE
gamma = 0.0
......@@ -318,8 +333,9 @@ module obl_k_epsilon
end subroutine solve_TKE_eq_semiimplicit
subroutine solve_eps_eq_semiimplicit (Field, Km, nz, dz, dt, P_TKE, B_TKE, TKE)
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]
......@@ -341,12 +357,12 @@ module obl_k_epsilon
real :: TKE(nz) !< TKE, [J/kg]
Kvisc(:) = Km(:) / Sch_eps
Kvisc(:) = Km(:) / param%Sch_eps
gamma = 0.0
do k = 1, nz
f_right(k) = Field(k) / TKE(k) * (C1_eps * P_TKE(k) + C3_eps * B_TKE(k))
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
......@@ -369,7 +385,7 @@ module obl_k_epsilon
do k = 2, nz-1
a(k) = dt / dz(k) * Kvisc_minus(k)/dz_minus(k) * (1.0 - gamma)
b(k) = -1.0 - C2_eps * dt * Field(k)/TKE(k) - &
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)
......
obl_main.f90
+ 14
− 13
View file @ 95c99b68
......@@ -47,6 +47,7 @@ program obl_main
! turbulence closure parameters
type(pphParamType) :: param_pph
type(pphDynParamType) :: param_pph_dyn
type(kepsilonParamType) :: param_k_epsilon
!< surface forcing
type(timeForcingDataType) :: sensible_hflux_surf, latent_hflux_surf !< heat fluxes, [W/m^2]
......@@ -241,8 +242,8 @@ program obl_main
!initialization of TKE & eps in case of k-epsilon closure
if (closure_mode.eq.3 .or. closure_mode.eq.4) then
call TKE_init(TKE, grid%cz)
call eps_init(EPS, grid%cz, grid%height)
call TKE_init(TKE, param_k_epsilon, grid%cz)
call eps_init(EPS, param_k_epsilon, grid%cz, grid%height)
endif
!< setting atmospheric forcing
......@@ -382,20 +383,20 @@ program obl_main
else if (closure_mode.eq.4) then
call TKE_bc(TKE, &
bc%U_dyn0, bc%U_dynH, grid%cz)
bc%U_dyn0, bc%U_dynH, param_k_epsilon, grid%cz)
call EPS_bc(EPS, &
bc%U_dyn0, bc%U_dynH, grid%dz, grid%cz)
call get_Km_k_eps (Km, TKE, EPS, grid%cz)
call get_Kh_k_eps (Kh, Km, grid%cz)
bc%U_dyn0, bc%U_dynH, param_k_epsilon, grid%dz, grid%cz)
call get_Km_k_eps (Km, TKE, EPS, param_k_epsilon, grid%cz)
call get_Kh_k_eps (Kh, Km, param_k_epsilon, grid%cz)
call TKE_bc(TKE, &
bc%U_dyn0, bc%U_dynH, grid%cz)
bc%U_dyn0, bc%U_dynH, param_k_epsilon, grid%cz)
call eps_bc(EPS, &
bc%U_dyn0, bc%U_dynH, grid%dz, grid%cz)
!call solve_TKE_eq_semiimplicit (TKE, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, eps)
call solve_TKE_eq (TKE, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, EPS)
call limit_min_array(TKE, TKE_min, grid%cz)
call solve_eps_eq_semiimplicit (EPS, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, TKE)
call limit_min_array(EPS, eps_min, grid%cz)
bc%U_dyn0, bc%U_dynH, param_k_epsilon, grid%dz, grid%cz)
!call solve_TKE_eq_semiimplicit (TKE, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, eps, param_k_epsilon)
call solve_TKE_eq (TKE, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, EPS, param_k_epsilon)
call limit_min_array(TKE, param_k_epsilon%TKE_min, grid%cz)
call solve_eps_eq_semiimplicit (EPS, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, TKE, param_k_epsilon)
call limit_min_array(EPS, param_k_epsilon%eps_min, grid%cz)
endif
! ----------------------------------------------------------------------------
......
obl_scm.f90
+ 1
− 1
View file @ 95c99b68
......@@ -227,7 +227,7 @@ module obl_scm
! --------------------------------------------------------------------------------
do k = 1, cz
U_rhs(k) = - f * (V(k) - Vg)
U_rhs(k) = f * (V(k) - Vg)
V_rhs(k) = - f * (U(k) - Vg)
end do
end subroutine
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment