diff --git a/CMakeLists.txt b/CMakeLists.txt
index f2d4ef2d87474b1f0c9a93b69ec46842bd77127f..1ea93b8d37358c65720167f2e88ec66082241da7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -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
diff --git a/obl_config.f90 b/obl_config.f90
index 24efc529eda2c1520fef9fe58e6f00ff27d655fc..cf81293d674c98ceae89b6da1115975cb0595216 100644
--- a/obl_config.f90
+++ b/obl_config.f90
@@ -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
diff --git a/obl_main.f90 b/obl_main.f90
index a9786c7e72c4c2b14b8a9251697f33890a4f24f6..9a63c122b6fc243c572e45d7ba498117fc460d27 100644
--- a/obl_main.f90
+++ b/obl_main.f90
@@ -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
! ----------------------------------------------------------------------------
diff --git a/obl_most.f90 b/obl_most.f90
new file mode 100644
index 0000000000000000000000000000000000000000..d493864f63f6c4a5efdadaf63171760cc085e363
--- /dev/null
+++ b/obl_most.f90
@@ -0,0 +1,374 @@
+#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