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

turbulence common module update

parent 6663a5b1
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 @ 16c558b3
......@@ -42,7 +42,7 @@ set(SOURCES
obl_state_eq.f90
obl_init.f90
obl_bc.f90
obl_turb_closure.f90
obl_turb_common.f90
obl_diag.f90
obl_k_epsilon.f90
obl_pph.f90
......
obl_k_epsilon.f90
+ 6
− 4
View file @ 16c558b3
......@@ -5,7 +5,7 @@ module obl_k_epsilon
use obl_math
use obl_grid
use obl_bc
use obl_turb_closure
use obl_turb_common
! directives list
implicit none
......@@ -22,6 +22,8 @@ module obl_k_epsilon
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
real :: c_S2_min = 1e-5 !< min shear frequency, [(1/s)**2]
end type
contains
......@@ -39,9 +41,9 @@ module obl_k_epsilon
call get_N2(N2, Rho, bc%Rho_dyn0, bc%Rho_dynH, &
param%kappa, param%PrT0, grid)
call get_S2(S2, U, V, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%u_momentum_flux0, grid)
call get_Ri_gradient(Ri_grad, N2, S2, 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
! --------------------------------------------------------------------------------
......
obl_pph.f90
+ 6
− 4
View file @ 16c558b3
......@@ -4,7 +4,7 @@ module obl_pph
! modules used
use obl_grid
use obl_bc
use obl_turb_closure
use obl_turb_common
! directives list
implicit none
......@@ -29,6 +29,8 @@ module obl_pph
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
......@@ -48,9 +50,9 @@ module obl_pph
call get_N2(N2, Rho, bc%rho_dyn0, bc%rho_dynH, &
param%kappa, param%PrT0, grid)
call get_S2(S2, U, V, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%u_momentum_flux0, grid)
call get_Ri_gradient(Ri_grad, N2, S2, 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
! --------------------------------------------------------------------------------
......
obl_pph_dyn.f90
+ 6
− 4
View file @ 16c558b3
......@@ -5,7 +5,7 @@ module obl_pph_dyn
use obl_grid
use obl_bc
use obl_diag
use obl_turb_closure
use obl_turb_common
! directives list
implicit none
......@@ -19,6 +19,8 @@ module obl_pph_dyn
real :: kappa = 0.4 !< von Karman constant
real :: PrT0 = 1.0 !< turbulent Prandtl number, *: probably redundant
real :: c_S2_min = 1e-5 !< min shear frequency, [(1/s)**2]
end type
contains
......@@ -36,9 +38,9 @@ module obl_pph_dyn
call get_N2(N2, Rho, bc%Rho_dyn0, bc%Rho_dynH, &
param%kappa, param%PrT0, grid)
call get_S2(S2, U, V, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%u_momentum_flux0, grid)
call get_Ri_gradient(Ri_grad, N2, S2, 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
! --------------------------------------------------------------------------------
......
obl_turb_closure.f90 → obl_turb_common.f90
+ 81
− 96
View file @ 16c558b3
module obl_turb_closure
!< @brief calculate turbulent parameters
module obl_turb_common
!< @brief common turbulence model subroutines
! --------------------------------------------------------------------------------
! modules used
! --------------------------------------------------------------------------------
use obl_grid
use obl_math
use obl_common_phys
use obl_state_eq
! directives list
implicit none
......@@ -15,19 +12,13 @@ module obl_turb_closure
! public interface
! --------------------------------------------------------------------------------
public :: get_N2_vec, get_S2_vec, get_Ri_gradient_vec
public :: get_N2, get_S2, get_Ri_gradient
public :: get_N2_vec, get_S2_vec, get_Ri_gradient_vec
public :: get_TKE_production, get_TKE_buoyancy
public :: get_TKE_production_vec, get_TKE_buoyancy_vec
! --------------------------------------------------------------------------------
real, parameter :: Ri_grad_min = 0.0000001 !< Richardson gradient number mainimum value
real, parameter :: Ri_grad_max = 10 !< Richardson gradient number maximum value
real, parameter :: c_S2_min = 1e-5 !1e-5
real, parameter :: kappa_ri = 0.4 !? где задавать каппу?
real, parameter :: PrT_ri = 0.8 !? прандль будет здесь публичный!!!
contains
......@@ -35,10 +26,13 @@ module obl_turb_closure
subroutine get_N2(N2, Rho, rho_dyn0, rho_dynH, kappa, PrT0, grid)
!< @brief calculate square of buoyancy frequency
! ----------------------------------------------------------------------------
use obl_state_eq
use obl_common_phys
type (gridDataType), intent(in) :: grid
real, intent(out) :: N2(grid%cz) !< square of buoyancy frequency, [s**(-2)]
real, intent(in) :: Rho(grid%cz) !< density, [kg/m**3]
real, dimension(grid%cz), intent(out) :: N2 !< square of buoyancy frequency, [s**(-2)]
real, dimension(grid%cz), intent(in) :: Rho !< density, [kg/m**3]
real, intent(in) :: rho_dyn0, rho_dynH !< density dynamic scales, [kg/m**3]
real, intent(in) :: kappa !< von Karman constant, [n/d]
......@@ -63,6 +57,9 @@ module obl_turb_closure
subroutine get_N2_vec(N2, Rho, rho_dyn0, rho_dynH, kappa, PrT0, dz, cz)
!< @brief calculate square of buoyancy frequency
! ----------------------------------------------------------------------------
use obl_state_eq
use obl_common_phys
integer, intent(in) :: cz !< number of z-steps
real, dimension(cz), intent(out) :: N2 !< square of buoyancy frequency, [s**(-2)]
......@@ -73,7 +70,7 @@ module obl_turb_closure
real, intent(in) :: kappa !< von Karman constant, [n/d]
real, intent(in) :: PrT0 !< neutral turbulent Prandtl number, [n/d]
real :: dzph(cz), dzmh(cz) !< half-sums of grid steps, [m]
real, dimension(cz) :: dzph, dzmh !< half-sums of grid steps, [m]
integer :: k
! ----------------------------------------------------------------------------
......@@ -97,107 +94,95 @@ module obl_turb_closure
end subroutine
! ----------------------------------------------------------------------------
subroutine get_S2(S2, U, V, flux_u_surf, flux_u_bot, flux_v_surf, flux_v_bot, grid)
subroutine get_S2(S2, U, V, U_dyn0, U_dynH, kappa, grid)
!< @brief calculate square of shear frequency
! ----------------------------------------------------------------------------
type (gridDataType), intent(in) :: grid
real, intent(out) :: S2(grid%cz) !< square of shear frequency, [s**(-2)]
real, intent(in) :: U(grid%cz), V(grid%cz) !< currents, [m/s]
real :: dU_tmp(grid%cz), dV_tmp(grid%cz) !< temporary for currents
real, intent(in) :: flux_u_surf, flux_u_bot !< [(m/s)**2]
real, intent(in) :: flux_v_surf, flux_v_bot !< [(m/s)**2]
real, dimension(grid%cz), intent(out) :: S2 !< square of shear frequency, [s**(-2)]
real, dimension(grid%cz), intent(in) :: U, V !< currents velocity, [m/s]
real :: dz_plus(grid%cz), dz_minus(grid%cz) !< half-sums for depth-step dz
integer :: k !< counter
! ----------------------------------------------------------------------------
real, intent(in) :: U_dyn0, U_dynH !< dynamic velocity, [m/s]
real, intent(in) :: kappa !< von Karman constant, [n/d]
do k = 1, grid%cz - 1
dz_plus(k) = (grid%dz(k)+grid%dz(k+1))/2.
end do
do k = 2, grid%cz
dz_minus(k) = (grid%dz(k)+grid%dz(k-1))/2.
end do
real :: Ugrad, Vgrad
integer :: k
! ----------------------------------------------------------------------------
do k = 2, grid%cz - 1
dU_tmp(k) = 0.5 * ((U(k+1) - U(k)) / (dz_plus(k)) + (U(k) - U(k-1)) / dz_minus(k))
dV_tmp(k) = 0.5 * ((V(k+1) - V(k)) / (dz_plus(k)) + (V(k) - V(k-1))/ dz_minus(k))
S2(k) = dU_tmp(k)**2 + dV_tmp(k)**2
Ugrad = 0.5 * ((U(k+1) - U(k)) / (grid%dzph(k)) + (U(k) - U(k-1)) / grid%dzmh(k))
Vgrad = 0.5 * ((V(k+1) - V(k)) / (grid%dzph(k)) + (V(k) - V(k-1)) / grid%dzmh(k))
S2(k) = Ugrad**2 + Vgrad**2
end do
dU_tmp(1) = sqrt(flux_u_bot) / (kappa_ri * grid%dz(1)/2.0)
dV_tmp(1) = sqrt(flux_v_bot) / (kappa_ri * grid%dz(1)/2.0)
S2(1) = dU_tmp(1)**2.0 + dV_tmp(1)**2.0
dU_tmp(grid%cz) = sqrt(flux_v_surf) / (kappa_ri * grid%dz(grid%cz)/2.0)
dV_tmp(grid%cz) = sqrt(flux_v_surf) / (kappa_ri * grid%dz(grid%cz)/2.0)
S2(grid%cz) = dU_tmp(grid%cz)**2.0 + dV_tmp(grid%cz)**2.0
!< bottom layer
S2(1) = (U_dyn0 / (kappa * grid%dz(1)/2.0))**2
!< surface layer
S2(grid%cz) = (U_dynH / (kappa * grid%dz(grid%cz)/2.0))**2
call limit_min_array(S2, c_S2_min, grid%cz)
end subroutine
! ----------------------------------------------------------------------------
subroutine get_S2_vec(S2, U, V, flux_u_surf, flux_u_bot, flux_v_surf, flux_v_bot, dz, nz)
subroutine get_S2_vec(S2, U, V, U_dyn0, U_dynH, kappa, dz, cz)
!< @brief calculate square of shear frequency
! ----------------------------------------------------------------------------
real, intent(out) :: S2(nz) !< square of shear frequency, [s**(-2)]
real, intent(in) :: U(nz), V(nz) !< currents, [m/s]
real :: dU_tmp(nz), dV_tmp(nz) !< temporary for currents
integer, intent(in) :: nz !< number of z-steps
real, intent(in) :: dz(nz) !< depth(z) step [m]
real, intent(in) :: flux_u_surf, flux_u_bot !< [(m/s)**2]
real, intent(in) :: flux_v_surf, flux_v_bot !< [(m/s)**2]
real :: dz_plus(nz), dz_minus(nz) !< half-sums for depth-step dz
integer :: k !< counter
integer, intent(in) :: cz !< number of z-steps
real, dimension(cz), intent(out) :: S2 !< square of shear frequency, [s**(-2)]
real, dimension(cz), intent(in) :: U, V !< currents velocity, [m/s]
real, dimension(cz), intent(in) :: dz !< grid steps [m]
real, intent(in) :: U_dyn0, U_dynH !< dynamic velocity, [m/s]
real, intent(in) :: kappa !< von Karman constant, [n/d]
real, dimension(cz) :: dzph, dzmh !< half-sums of grid steps, [m]
real :: Ugrad, Vgrad
integer :: k
! ----------------------------------------------------------------------------
do k = 1, nz-1
dz_plus(k) = (dz(k)+dz(k+1))/2.
do k = 1, cz - 1
dzph(k) = (dz(k) + dz(k+1)) / 2.0
end do
do k = 2, nz
dz_minus(k) = (dz(k)+dz(k-1))/2.
do k = 2, cz
dzmh(k) = (dz(k) + dz(k-1)) / 2.0
end do
do k = 2, nz-1
dU_tmp(k) = 0.5 * ((U(k+1) - U(k)) / (dz_plus(k)) + (U(k) - U(k-1)) / dz_minus(k))
dV_tmp(k) = 0.5 * ((V(k+1) - V(k)) / (dz_plus(k)) + (V(k) - V(k-1))/ dz_minus(k))
S2(k) = dU_tmp(k)**2 + dV_tmp(k)**2
do k = 2, cz - 1
Ugrad = 0.5 * ((U(k+1) - U(k)) / (dzph(k)) + (U(k) - U(k-1)) / dzmh(k))
Vgrad = 0.5 * ((V(k+1) - V(k)) / (dzph(k)) + (V(k) - V(k-1)) / dzmh(k))
S2(k) = Ugrad**2 + Vgrad**2
end do
dU_tmp(1) = sqrt(flux_u_bot) / (kappa_ri * dz(1)/2.0)
dV_tmp(1) = sqrt(flux_v_bot) / (kappa_ri * dz(1)/2.0)
S2(1) = dU_tmp(1)**2.0 + dV_tmp(1)**2.0
dU_tmp(nz) = sqrt(flux_v_surf) / (kappa_ri * dz(nz)/2.0)
dV_tmp(nz) = sqrt(flux_v_surf) / (kappa_ri * dz(nz)/2.0)
S2(nz) = dU_tmp(nz)**2.0 + dV_tmp(nz)**2.0
call limit_min_array(S2, c_S2_min, nz)
!< bottom layer
S2(1) = (U_dyn0 / (kappa * dz(1)/2.0))**2
!< surface layer
S2(cz) = (U_dynH / (kappa * dz(cz)/2.0))**2
end subroutine
! ----------------------------------------------------------------------------
subroutine get_Ri_gradient(Ri_grad, N2, S2, grid)
subroutine get_Ri_gradient(Ri_grad, N2, S2, c_S2_min, grid)
!< @brief calculate Richardson gradient number
! ----------------------------------------------------------------------------
type (gridDataType), intent(in) :: grid
real, intent(out) :: Ri_grad(grid%cz) !< Richardson gradient number
real, intent(in) :: N2(grid%cz), S2(grid%cz) !< square of buoyancy and shear frequency, [s**(-2)]
real, dimension(grid%cz), intent(out) :: Ri_grad !< Richardson gradient number
real, dimension(grid%cz), intent(in) :: N2, S2 !< square of buoyancy and shear frequency, [s**(-2)]
real, intent(in) :: c_S2_min !< min S2 value
! ----------------------------------------------------------------------------
call get_Ri_gradient_vec(Ri_grad, N2, S2, grid%cz)
call get_Ri_gradient_vec(Ri_grad, N2, S2, c_S2_min, grid%cz)
end subroutine
! ----------------------------------------------------------------------------
subroutine get_Ri_gradient_vec(Ri_grad, N2, S2, cz)
subroutine get_Ri_gradient_vec(Ri_grad, N2, S2, c_S2_min, cz)
!< @brief calculate Richardson gradient number
! ----------------------------------------------------------------------------
integer, intent(in) :: cz !< number of z-steps
real, intent(out) :: Ri_grad(cz) !< Richardson gradient number
real, intent(in) :: N2(cz), S2(cz) !< square of buoyancy and shear frequency, [s**(-2)]
real, dimension(cz), intent(out) :: Ri_grad !< Richardson gradient number
real, dimension(cz), intent(in) :: N2, S2 !< square of buoyancy and shear frequency, [s**(-2)]
real, intent(in) :: c_S2_min !< min S2 value
integer :: k
! ----------------------------------------------------------------------------
......@@ -217,9 +202,9 @@ module obl_turb_closure
! ----------------------------------------------------------------------------
type (gridDataType), intent(in) :: grid
real, intent(in) :: S2(grid%cz) !< shear frequency ^ 2, [s^-2]
real, intent(in) :: Km(grid%cz) !< eddy viscosity for momentum, [m**2 / s]
real, intent(out) :: P_TKE(grid%cz) !< shear production of TKE, [m**2 / s**3]
real, dimension(grid%cz), intent(in) :: S2 !< shear frequency ^ 2, [s^-2]
real, dimension(grid%cz), intent(in) :: Km !< eddy viscosity for momentum, [m**2 / s]
real, dimension(grid%cz), intent(out) :: P_TKE !< shear production of TKE, [m**2 / s**3]
! ----------------------------------------------------------------------------
call get_TKE_production_vec(P_TKE, Km, S2, grid%cz)
......@@ -230,9 +215,9 @@ module obl_turb_closure
!< @brief calculation of TKE production by shear
! ----------------------------------------------------------------------------
integer, intent(in) :: cz !< number of z-steps
real, intent(in) :: S2(cz) !< shear frequency ^ 2, [s^-2]
real, intent(in) :: Km(cz) !< eddy viscosity for momentum, [m**2 / s]
real, intent(out) :: P_TKE(cz) !< shear production of TKE, [m**2 / s**3]
real, dimension(cz), intent(in) :: S2 !< shear frequency ^ 2, [s^-2]
real, dimension(cz), intent(in) :: Km !< eddy viscosity for momentum, [m**2 / s]
real, dimension(cz), intent(out) :: P_TKE !< shear production of TKE, [m**2 / s**3]
integer :: k
! ----------------------------------------------------------------------------
......@@ -248,9 +233,9 @@ module obl_turb_closure
! ----------------------------------------------------------------------------
type (gridDataType), intent(in) :: grid
real, intent(in) :: N2(grid%cz) !< buoyancy frequency ^ 2, [s^-2]
real, intent(in) :: Kh(grid%cz) !< eddy duffusivity for tracers, [m**2 / s]
real, intent(out) :: B_TKE(grid%cz) !< buoyancy production of TKE, [m**2 / s**3]
real, dimension(grid%cz), intent(in) :: N2 !< buoyancy frequency ^ 2, [s^-2]
real, dimension(grid%cz), intent(in) :: Kh !< eddy duffusivity for tracers, [m**2 / s]
real, dimension(grid%cz), intent(out) :: B_TKE !< buoyancy production of TKE, [m**2 / s**3]
! ----------------------------------------------------------------------------
call get_TKE_buoyancy_vec(B_TKE, Kh, N2, grid%cz)
......@@ -261,9 +246,9 @@ module obl_turb_closure
!< @brief calculation of TKE buoyancy production
! ----------------------------------------------------------------------------
integer, intent(in) :: cz !< number of z-steps
real, intent(in) :: N2(cz) !< buoyancy frequency ^ 2, [s^-2]
real, intent(in) :: Kh(cz) !< eddy duffusivity for tracers, [m**2 / s]
real, intent(out) :: B_TKE(cz) !< buoyancy production of TKE, [m**2 / s**3]
real, dimension(cz), intent(in) :: N2 !< buoyancy frequency ^ 2, [s^-2]
real, dimension(cz), intent(in) :: Kh !< eddy duffusivity for tracers, [m**2 / s]
real, dimension(cz), intent(out) :: B_TKE !< buoyancy production of TKE, [m**2 / s**3]
integer :: k
! ----------------------------------------------------------------------------
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment