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

adding most based turbulence closure

parent 5aefb264
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 1
− 0
View file @ 7b7685ea
......@@ -48,6 +48,7 @@ set(SOURCES
obl_k_epsilon.f90
obl_pph.f90
obl_pph_dyn.f90
obl_most.f90
obl_fluxes.f90
obl_scm.f90
obl_main.f90
......
obl_config.f90
+ 6
− 0
View file @ 7b7685ea
......@@ -35,10 +35,12 @@ module obl_config
integer, parameter :: obl_model_pph = 0 !< pacanowski-philander
integer, parameter :: obl_model_pph_dyn = 1 !< pacanowski-philander (dynamic)
integer, parameter :: obl_model_k_epsilon = 2 !< k-epsilon
integer, parameter :: obl_model_most = 3 !< most scheme
character(len = 16), parameter :: obl_model_pph_tag = 'pph'
character(len = 16), parameter :: obl_model_pph_dyn_tag = 'pph-dyn'
character(len = 16), parameter :: obl_model_k_epsilon_tag = 'k-epsilon'
character(len = 16), parameter :: obl_model_most_tag = 'most'
contains
......@@ -97,6 +99,8 @@ contains
id = obl_model_pph_dyn
else if (trim(tag) == trim(obl_model_k_epsilon_tag)) then
id = obl_model_k_epsilon
else if (trim(tag) == trim(obl_model_most_tag)) then
id = obl_model_most
end if
end function
......@@ -113,6 +117,8 @@ contains
tag = obl_model_pph_dyn_tag
else if (id == obl_model_k_epsilon) then
tag = obl_model_k_epsilon_tag
else if (id == obl_model_most) then
tag = obl_model_most_tag
end if
end function
......
obl_main.f90
+ 14
− 1
View file @ 7b7685ea
......@@ -41,6 +41,11 @@ program obl_main
init_turbulence_closure_k_epsilon => init_turbulence_closure, &
advance_turbulence_closure_k_epsilon => advance_turbulence_closure, &
TKE_init, EPS_init, TKE_bc, EPS_bc
use obl_most, only: mostParamType, &
set_config_param_most => set_config_param, &
define_stability_functions_most => define_stability_functions, &
init_turbulence_closure_most => init_turbulence_closure, &
advance_turbulence_closure_most => advance_turbulence_closure
use obl_io_plt
use io
use io_metadata
......@@ -69,6 +74,7 @@ program obl_main
type(pphParamType) :: param_pph
type(pphDynParamType) :: param_pph_dyn
type(kepsilonParamType) :: param_k_epsilon
type(mostParamType) :: param_most
!< boundary conditions data
type(oblBcType) :: bc
......@@ -131,6 +137,7 @@ program obl_main
!< = obl_model_pph, pacanowski-philander
!< = obl_model_pph_dyn, pacanowski-philander (dynamic)
!< = obl_model_k_epsilon, k-epsilon
!< = obl_model_most, MOST scheme
obl_model_id = obl_model_k_epsilon
!< default output = 1, (netcdf)
......@@ -153,7 +160,7 @@ program obl_main
write(*, *) ' builtin setup [key] = kato (default) || papa-fluxes || papa-meteo || cbl || cyclone'
write(*, *) ' use configuration file [filename]'
write(*, *) ' --model [key]'
write(*, *) ' [key] = pph || pph-dyn || k-epsilon (default)'
write(*, *) ' [key] = pph || pph-dyn || k-epsilon (default) || most'
return
end if
......@@ -297,6 +304,9 @@ program obl_main
else if (obl_model_id.eq.obl_model_k_epsilon) then
call define_stability_functions_k_epsilon(param_k_epsilon, bc, grid)
call init_turbulence_closure_k_epsilon(param_k_epsilon, bc, grid)
else if (obl_model_id.eq.obl_model_most) then
call define_stability_functions_most(param_most, bc, grid)
call init_turbulence_closure_most(param_most, bc, grid)
endif
! ----------------------------------------------------------------------------
......@@ -326,6 +336,9 @@ program obl_main
else if (obl_model_id.eq.obl_model_k_epsilon) then
call define_stability_functions_k_epsilon(param_k_epsilon, bc, grid)
call advance_turbulence_closure_k_epsilon(param_k_epsilon, bc, grid, dt)
else if (obl_model_id.eq.obl_model_most) then
call define_stability_functions_most(param_most, bc, grid)
call advance_turbulence_closure_most(param_most, bc, grid, dt)
endif
! ----------------------------------------------------------------------------
......
obl_most.f90 0 → 100644
+ 374
− 0
View file @ 7b7685ea
#include "obl_def.fi"
module obl_most
!< @brief MOST based scheme
! --------------------------------------------------------------------------------
! modules used
! --------------------------------------------------------------------------------
#ifdef USE_CONFIG_PARSER
use iso_c_binding, only: C_NULL_CHAR
use config_parser
#endif
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 :: set_config_param
! --------------------------------------------------------------------------------
!> @brief MOST scheme parameters
type, public :: mostParamType
real :: Ri_g_cr = 0.20 !< critical gradient Richardson number [n/d]
!< stability functions in unstable case
real :: Ri_g_m = -0.20
real :: Ri_g_h = -1.0
real :: a_m = 16.0
real :: b_m = 1.26
real :: c_m = 8.38
real :: a_h = 16.0
real :: b_h = -28.86
real :: c_h = 98.96
!<
real :: Km_b = 0.0001 !< background eddy viscosity, [m**2 / s]
real :: Kh_b = 0.00005 !< background eddy diffusivity, [m**2 / s]
real :: kappa = 0.4 !< von Karman constant, [n/d]
real :: PrT0 = 1.0 !< neutral Prandtl number, [n/d]
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(mostParamType), 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(mostParamType), intent(in) :: param
type(oblBcType), intent(in) :: bc
type (gridDataType), intent(in) :: grid
! --------------------------------------------------------------------------------
call get_eddy_viscosity(Km, Ri_grad, S2, bc%U_dynH, param, grid)
call get_eddy_diffusivity(Kh, Ri_grad, S2, bc%U_dynH, param, grid)
end subroutine
! --------------------------------------------------------------------------------
subroutine advance_turbulence_closure(param, bc, grid, dt)
!< @brief advance turbulence closure
! ----------------------------------------------------------------------------
use obl_state
type(mostParamType), 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 get_eddy_viscosity(Km, Ri_grad, S2, bc%U_dynH, param, grid)
call get_eddy_diffusivity(Kh, Ri_grad, S2, bc%U_dynH, param, grid)
end subroutine
! --------------------------------------------------------------------------------
subroutine get_eddy_viscosity(Km, Ri_grad, S2, U_dynH, param, grid)
!< @brief calculate eddy viscosity on grid
! ----------------------------------------------------------------------------
type(mostParamType), intent(in) :: param
type (gridDataType), intent(in) :: grid
real, dimension(grid%cz), intent(in) :: Ri_grad !< Richardson gradient number
real, dimension(grid%cz), intent(in) :: S2 !< square of shear frequency [1 / s**2]
real, intent(in) :: U_dynH !< surface velocity scale [m/s]
real, dimension(grid%cz), intent(out) :: Km !< eddy viscosity, [m**2 / s]
real, dimension(grid%cz) :: depth
integer :: k
! --------------------------------------------------------------------------------
do k = 1, grid%cz
depth(k) = grid%zpos + grid%height - grid%z(k)
end do
call get_eddy_viscosity_vec(Km, Ri_grad, S2, U_dynH, param, depth, grid%cz)
end subroutine
subroutine get_eddy_viscosity_vec(Km, Ri_grad, S2, U_dynH, param, depth, n)
!< @brief calculate eddy viscosity
! ----------------------------------------------------------------------------
type(mostParamType), intent(in) :: param
integer, intent(in) :: n !< vector size
real, dimension(n), intent(in) :: Ri_grad !< Richardson gradient number
real, dimension(n), intent(in) :: S2 !< square of shear frequency [1 / s**2]
real, intent(in) :: U_dynH !< surface velocity scale [m/s]
real, dimension(n), intent(in) :: depth !< depth coordinates, [m]
real, dimension(n), intent(out) :: Km !< eddy viscosity, [m**2 / s]
real :: l_m
real :: g_m
integer :: k
! --------------------------------------------------------------------------------
do k = 1, n
if (Ri_grad(k) >= 0.0) then
!< stable case
Km(k) = param%Km_b
if (Ri_grad(k) <= param%Ri_g_cr) then
l_m = param%kappa * depth(k) * (1.0 - Ri_grad(k) / param%Ri_g_cr)
Km(k) = Km(k) + l_m * l_m * sqrt(S2(k))
endif
else
!< unstable case
if (Ri_grad(k) >= param%Ri_g_m) then
g_m = (1.0 - param%a_m * Ri_grad(k))**0.25
else
g_m = (param%b_m - param%c_m * Ri_grad(k))**(1.0 / 3.0)
endif
l_m = param%kappa * depth(k) * g_m
Km(k) = l_m * l_m * sqrt(max(S2(k), param%c_S2_min))
endif
end do
end subroutine
! --------------------------------------------------------------------------------
subroutine get_eddy_diffusivity(Kh, Ri_grad, S2, U_dynH, param, grid)
!< @brief calculate eddy diffusivity on grid
! ----------------------------------------------------------------------------
type(mostParamType), intent(in) :: param
type (gridDataType), intent(in) :: grid
real, dimension(grid%cz), intent(in) :: Ri_grad !< Richardson gradient number
real, dimension(grid%cz), intent(in) :: S2 !< square of shear frequency [1 / s**2]
real, intent(in) :: U_dynH !< surface velocity scale [m/s]
real, dimension(grid%cz), intent(out) :: Kh !< eddy diffusivity, [m**2 / s]
real, dimension(grid%cz) :: depth
integer :: k
! --------------------------------------------------------------------------------
do k = 1, grid%cz
depth(k) = grid%zpos + grid%height - grid%z(k)
end do
call get_eddy_diffusivity_vec(Kh, Ri_grad, S2, U_dynH, param, depth, grid%cz)
end subroutine
subroutine get_eddy_diffusivity_vec(Kh, Ri_grad, S2, U_dynH, param, depth, n)
!< @brief calculate eddy diffusivity
! ----------------------------------------------------------------------------
type(mostParamType), intent(in) :: param
integer, intent(in) :: n !< vector size
real, dimension(n), intent(in) :: Ri_grad !< Richardson gradient number
real, dimension(n), intent(in) :: S2 !< square of shear frequency [1 / s**2]
real, intent(in) :: U_dynH !< surface velocity scale [m/s]
real, dimension(n), intent(in) :: depth !< depth coordinates, [m]
real, dimension(n), intent(out) :: Kh !< eddy diffusivity, [m**2 / s]
real :: l_m, l_h
real :: g_m, g_h
integer :: k
! --------------------------------------------------------------------------------
do k = 1, n
if (Ri_grad(k) >= 0.0) then
!< stable case
Kh(k) = param%Kh_b
if (Ri_grad(k) <= param%Ri_g_cr) then
l_m = param%kappa * depth(k) * (1.0 - Ri_grad(k) / param%Ri_g_cr)
Kh(k) = Kh(k) + (1.0 / param%PrT0) * l_m * l_m * sqrt(S2(k))
endif
else
!< unstable case
if (Ri_grad(k) >= param%Ri_g_h) then
g_m = (1.0 - param%a_m * Ri_grad(k))**0.25
g_h = (1.0 - param%a_h * Ri_grad(k))**0.5
else
g_m = (param%b_m - param%c_m * Ri_grad(k))**(1.0 / 3.0)
g_h = (param%b_h - param%c_h * Ri_grad(k))**(1.0 / 3.0)
endif
l_m = param%kappa * depth(k) * g_m
l_h = param%kappa * depth(k) * g_h
Kh(k) = (1.0 / param%PrT0) * l_m * l_h * sqrt(max(S2(k), param%c_S2_min))
endif
end do
end subroutine
! --------------------------------------------------------------------------------
subroutine set_config_param(param, tag, ierr)
!< @brief set turbulence closure parameters
! ----------------------------------------------------------------------------
type(mostParamType), 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)//".Ri_g_cr"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".Ri_g_cr"//C_NULL_CHAR, param%Ri_g_cr, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".Ri_g_m"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".Ri_g_m"//C_NULL_CHAR, param%Ri_g_m, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".Ri_g_h"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".Ri_g_h"//C_NULL_CHAR, param%Ri_g_h, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".a_m"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".a_m"//C_NULL_CHAR, param%a_m, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".b_m"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".b_m"//C_NULL_CHAR, param%b_m, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".c_m"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".c_m"//C_NULL_CHAR, param%c_m, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".a_h"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".a_h"//C_NULL_CHAR, param%a_h, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".b_h"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".b_h"//C_NULL_CHAR, param%b_h, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".c_h"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".c_h"//C_NULL_CHAR, param%c_h, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".Km_b"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".Km_b"//C_NULL_CHAR, param%Km_b, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
endif
call c_config_is_varname(trim(tag)//".Kh_b"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float(trim(tag)//".Kh_b"//C_NULL_CHAR, param%Kh_b, 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment