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

k-epsilon model implementation update

parent e2b4de82
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 1
− 1
View file @ dc783f58
......@@ -34,7 +34,7 @@ set(SOURCES
obl_k_epsilon.f90
obl_pacanowski_philander.f90
obl_pacanowski_philander_plus.f90
obl_equation.f90
obl_scm.f90
obl_main.f90
io_metadata.f90
io.f90
......
obl_k_epsilon.f90
+ 20
− 4
View file @ dc783f58
......@@ -129,12 +129,14 @@ module obl_k_epsilon
end do
end subroutine get_Km_eps
subroutine solve_TKE_eq (Field, Kvisc, nz, dz, dt, P_TKE, B_TKE, eps)
subroutine solve_TKE_eq (Field, Km, nz, dz, dt, P_TKE, B_TKE, eps)
!< @brief solver for TKE equation (explicit)
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
......@@ -148,6 +150,8 @@ 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
gamma = 0.0
do k = 1, nz
......@@ -196,12 +200,14 @@ module obl_k_epsilon
end subroutine solve_TKE_eq
subroutine solve_eps_eq (Field, Kvisc, nz, dz, dt, P_TKE, B_TKE, TKE)
subroutine solve_eps_eq (Field, Km, nz, dz, dt, P_TKE, B_TKE, TKE)
!< @brief solver for epsilon equation (explicit)
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
......@@ -216,6 +222,8 @@ module obl_k_epsilon
real :: eps(nz) !< TKE dissipation rate, [W/kg]
real :: TKE(nz) !< TKE, [J/kg]
Kvisc(:) = Km(:) / Sch_eps
gamma = 0.0
do k = 1, nz
......@@ -263,12 +271,14 @@ module obl_k_epsilon
end subroutine solve_eps_eq
subroutine solve_TKE_eq_semiimplicit (Field, Kvisc, nz, dz, dt, P_TKE, B_TKE, eps)
subroutine solve_TKE_eq_semiimplicit (Field, Km, nz, dz, dt, P_TKE, B_TKE, eps)
!< @brief solver for TKE equation (semimplicit)
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
......@@ -282,6 +292,8 @@ 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
gamma = 0.0
do k = 1, nz
......@@ -331,12 +343,14 @@ module obl_k_epsilon
end subroutine solve_TKE_eq_semiimplicit
subroutine solve_eps_eq_semiimplicit (Field, Kvisc, nz, dz, dt, P_TKE, B_TKE, TKE)
subroutine solve_eps_eq_semiimplicit (Field, Km, nz, dz, dt, P_TKE, B_TKE, TKE)
!< @brief solver for epsilon equation (semiimplict)
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
......@@ -352,6 +366,8 @@ module obl_k_epsilon
real :: TKE(nz) !< TKE, [J/kg]
Kvisc(:) = Km(:) / Sch_eps
gamma = 0.0
do k = 1, nz
......
obl_main.f90
+ 4
− 14
View file @ dc783f58
......@@ -19,7 +19,7 @@ program obl_main
use obl_initial_conditions
use obl_forcing_and_boundary
use obl_rhs
use obl_equation
use obl_scm
use obl_turb_closure
use obl_diag
use obl_pacanowski_philander
......@@ -85,9 +85,6 @@ program obl_main
real :: mld !< mixed layer depth, [m]
real :: lab_mld !< theoretical mixed layer depth, [m]
real, allocatable :: Km_TKE(:) !eddy viscosity for momentum adjusted for Schmidt TKE number, [m**2 / s]
real, allocatable :: Km_eps(:) !eddy viscosity for momentum adjusted for Schmidt eps number, [m**2 / s]
!< just a temporary to hold value
real :: fvalue
......@@ -233,8 +230,6 @@ program obl_main
allocate(f_u_right(1:grid%cz))
allocate(f_v_right(1:grid%cz))
allocate(Km_TKE(1:grid%cz))
allocate(Km_eps(1:grid%cz))
......@@ -396,12 +391,10 @@ program obl_main
bc%u_momentum_fluxH, bc%v_momentum_fluxH, grid%cz)
call eps_bc(EPS, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, grid%cz, grid%dz)
call get_Km_TKE(Km_TKE, Km, grid%cz)
call get_Km_eps(Km_eps, Km, grid%cz)
!call solve_TKE_eq_semiimplicit (TKE, Km_TKE, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, eps)
call solve_TKE_eq (TKE, Km_TKE, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, EPS)
!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_eps, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, TKE)
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)
endif
! ----------------------------------------------------------------------------
......@@ -493,9 +486,6 @@ program obl_main
deallocate(f_u_right)
deallocate(f_v_right)
deallocate(Km_TKE)
deallocate(Km_eps)
!> removing time-dependent forcing data
call deallocate_tforcing(sensible_hflux_surf)
call deallocate_tforcing(latent_hflux_surf)
......
obl_equation.f90 → obl_scm.f90
+ 4
− 1
View file @ dc783f58
module obl_equation
module obl_scm
!< @brief solvers for general equations for scalars and vector
! modules used
use obl_math
use obl_state
use obl_boundary
implicit none
contains
subroutine solve_scalar_eq (Field, Kvisc, nz, dz, dt, flux_surf, flux_bot, f_right)
!< @brief solver for equation for scalar fields (temperature and salinity)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment