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

turbulence closure implementation update

parent c2d23107
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
obl_k_epsilon.f90
+ 7
− 6
View file @ 6663a5b1
......@@ -13,7 +13,7 @@ module obl_k_epsilon
!> @brief k-epsilon parameters
type, public :: kepsilonParamType
real :: C_mu = 0.09 !< momentum stability function constant
real :: PrT = 1.25 !< turbulent Prandtl number
real :: PrT0 = 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
......@@ -37,10 +37,11 @@ module obl_k_epsilon
type (gridDataType), intent(in) :: grid
! --------------------------------------------------------------------------------
call get_N2(N2, Rho, bc%Rho_dynH, bc%Rho_dyn0, 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_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%u_momentum_flux0, grid)
call get_Rig_vec(Ri_grad, N2, S2, grid%cz)
call get_Ri_gradient(Ri_grad, N2, S2, grid)
end subroutine
! --------------------------------------------------------------------------------
......@@ -55,8 +56,8 @@ module obl_k_epsilon
real, intent(in) :: dt
! --------------------------------------------------------------------------------
call get_TKE_production_vec(TKE_production, Km, S2, grid%cz)
call get_TKE_buoyancy_vec(TKE_buoyancy, Kh, N2, grid%cz)
call get_TKE_production(TKE_production, Km, S2, grid)
call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid)
call TKE_bc(TKE, &
bc%U_dyn0, bc%U_dynH, param, grid%cz)
......@@ -160,7 +161,7 @@ module obl_k_epsilon
integer :: k !< counter
do k = 1, nz
Kh(k) = Km(k) / param%PrT
Kh(k) = Km(k) / param%PrT0
end do
end subroutine get_Kh_k_eps
......
obl_pph.f90
+ 8
− 4
View file @ 6663a5b1
......@@ -26,6 +26,9 @@ module obl_pph
real :: n = 2.0 !< constant
real :: Km_b = 0.0001 !< background eddy viscosity, [m**2 / s], *: INMCM: 5 * 10**(-6)
real :: Kh_b = 0.00001 !< background eddy diffusivity, [m**2 / s]
real :: kappa = 0.4 !< von Karman constant, [n/d]
real :: PrT0 = 1.0 !< neutral Prandtl number, [n/d]
end type
contains
......@@ -43,10 +46,11 @@ module obl_pph
type (gridDataType), intent(in) :: grid
! --------------------------------------------------------------------------------
call get_N2(N2, Rho, bc%rho_dynH, bc%rho_dyn0, 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_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%u_momentum_flux0, grid)
call get_Rig_vec(Ri_grad, N2, S2, grid%cz)
call get_Ri_gradient(Ri_grad, N2, S2, grid)
end subroutine
! --------------------------------------------------------------------------------
......@@ -61,8 +65,8 @@ module obl_pph
real, intent(in) :: dt
! --------------------------------------------------------------------------------
call get_TKE_production_vec(TKE_production, Km, S2, grid%cz)
call get_TKE_buoyancy_vec(TKE_buoyancy, Kh, N2, grid%cz)
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, param, grid)
call get_eddy_diffusivity(Kh, Ri_grad, param, grid)
......
obl_pph_dyn.f90
+ 8
− 6
View file @ 6663a5b1
......@@ -16,8 +16,9 @@ module obl_pph_dyn
real :: n = 2.0 !< constant
real :: Km_b = 0.0001 !< background eddy viscosity, [m**2 / s], *: INMCM: 5 * 10**(-6)
real :: Kh_b = 0.00001 !< background eddy diffusivity, [m**2 / s]
real :: kappa = 0.4 !< von Karman constant
real :: PrT = 0.8 !< turbulent Prandtl number, *: probably redundant
real :: PrT0 = 1.0 !< turbulent Prandtl number, *: probably redundant
end type
contains
......@@ -33,10 +34,11 @@ module obl_pph_dyn
type (gridDataType), intent(in) :: grid
! --------------------------------------------------------------------------------
call get_N2(N2, Rho, bc%Rho_dynH, bc%Rho_dyn0, 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_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%u_momentum_flux0, grid)
call get_Rig_vec(Ri_grad, N2, S2, grid%cz)
call get_Ri_gradient(Ri_grad, N2, S2, grid)
end subroutine
! --------------------------------------------------------------------------------
......@@ -53,8 +55,8 @@ module obl_pph_dyn
real :: mld
! --------------------------------------------------------------------------------
call get_TKE_production_vec(TKE_production, Km, S2, grid%cz)
call get_TKE_buoyancy_vec(TKE_buoyancy, Kh, N2, grid%cz)
call get_TKE_production(TKE_production, Km, S2, grid)
call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid)
call get_mld(mld, N2, grid%dz, grid%cz)
......@@ -103,7 +105,7 @@ module obl_pph_dyn
integer :: k !< counter
do k = 1, nz
Kh(k) = Km(k) / param%PrT
Kh(k) = Km(k) / param%PrT0
end do
end subroutine
......
obl_turb_closure.f90
+ 159
− 102
View file @ 6663a5b1
......@@ -15,93 +15,131 @@ module obl_turb_closure
! public interface
! --------------------------------------------------------------------------------
public :: get_N2_vec, get_S2_vec, get_Rig_vec
public :: get_N2, get_S2
public :: get_N2_vec, get_S2_vec, get_Ri_gradient_vec
public :: get_N2, get_S2, get_Ri_gradient
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 :: S2_min = 1e-5 !1e-5
real, parameter :: c_S2_min = 1e-5 !1e-5
real, parameter :: kappa_ri = 0.4 !? где задавать каппу?
real, parameter :: PrT_ri = 0.8 !? прандль будет здесь публичный!!!
contains
subroutine get_N2_vec(N2, Rho, Rho_dyn_surf, Rho_dyn_bot, dz, nz)
! ----------------------------------------------------------------------------
subroutine get_N2(N2, Rho, rho_dyn0, rho_dynH, kappa, PrT0, grid)
!< @brief calculate square of buoyancy frequency
! ----------------------------------------------------------------------------
integer, intent(in) :: nz !< number of z-steps
real, intent(out), dimension(nz) :: N2 !< square of buoyancy frequency, [s**(-2)]
type (gridDataType), intent(in) :: grid
real, intent(in), dimension(nz) :: Rho !< density, [kg/m**3]
real, intent(in), dimension(nz) :: dz !< depth(z) step [m]
real, intent(in) :: Rho_dyn_surf, Rho_dyn_bot !< density, [kg/m**3]
real, intent(out) :: N2(grid%cz) !< square of buoyancy frequency, [s**(-2)]
real, intent(in) :: Rho(grid%cz) !< 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]
real, intent(in) :: PrT0 !< neutral turbulent Prandtl number, [n/d]
real :: dz_plus(nz), dz_minus(nz) !< half-sums for depth-step dz
real :: dRho_tmp_surf, dRho_tmp_bot
integer :: k
! ----------------------------------------------------------------------------
do k = 1, nz-1
dz_plus(k) = (dz(k)+dz(k+1))/2.
do k = 2, grid%cz - 1
N2(k) = - (g / Rho_ref) * (0.5 * &
((Rho(k+1) - Rho(k)) / grid%dzph(k) + (Rho(k) - Rho(k-1)) / grid%dzmh(k)))
end do
do k = 2, nz
dz_minus(k) = (dz(k)+dz(k-1))/2.
!< bottom layer
N2(1) = - (g / Rho_ref) * (1.0 / PrT0) * (rho_dyn0/ (kappa * grid%dz(1)/2.0))
!< surface layer
N2(grid%cz) = - (g / Rho_ref) * (1.0 / PrT0) * (rho_dynH/ (kappa * grid%dz(grid%cz)/2.0))
end subroutine
! ----------------------------------------------------------------------------
subroutine get_N2_vec(N2, Rho, rho_dyn0, rho_dynH, kappa, PrT0, dz, cz)
!< @brief calculate square of buoyancy frequency
! ----------------------------------------------------------------------------
integer, intent(in) :: cz !< number of z-steps
real, dimension(cz), intent(out) :: N2 !< square of buoyancy frequency, [s**(-2)]
real, dimension(cz), intent(in) :: Rho !< density, [kg/m**3]
real, dimension(cz), intent(in) :: dz !< grid steps, [m]
real, intent(in) :: rho_dyn0, rho_dynH !< density dynamic scales, [kg/m**3]
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]
integer :: k
! ----------------------------------------------------------------------------
do k = 1, cz - 1
dzph(k) = (dz(k) + dz(k+1)) / 2.0
end do
do k = 2, cz
dzmh(k) = (dz(k) + dz(k-1)) / 2.0
end do
do k = 2, nz-1
do k = 2, cz - 1
N2(k) = - g/Rho_ref * (0.5 * &
((Rho(k+1) - Rho(k)) / dz_plus(k) + (Rho(k) - Rho(k-1)) / dz_minus(k)))
((Rho(k+1) - Rho(k)) / dzph(k) + (Rho(k) - Rho(k-1)) / dzmh(k)))
end do
dRho_tmp_bot = 1.0 / PrT_ri * (Rho_dyn_bot/ (kappa_ri * dz(1)/2.0))
dRho_tmp_surf = 1.0 / PrT_ri * (Rho_dyn_surf/ (kappa_ri * dz(nz)/2.0))
N2(1) = - g/Rho_ref * dRho_tmp_bot
N2(nz) = - g/Rho_ref * dRho_tmp_surf
!< bottom layer
N2(1) = - (g / Rho_ref) * (1.0 / PrT0) * (rho_dyn0 / (kappa * dz(1)/2.0))
!< surface layer
N2(cz) = - (g / Rho_ref) * (1.0 / PrT0) * (rho_dynH / (kappa * dz(cz)/2.0))
end subroutine
subroutine get_N2(N2, Rho, Rho_dyn_surf, Rho_dyn_bot, grid)
!< @brief calculate square of buoyancy frequency
! ----------------------------------------------------------------------------
subroutine get_S2(S2, U, V, flux_u_surf, flux_u_bot, flux_v_surf, flux_v_bot, grid)
!< @brief calculate square of shear frequency
! ----------------------------------------------------------------------------
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, intent(in) :: Rho_dyn_surf, Rho_dyn_bot !< density, [kg/m**3]
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 :: dz_plus(grid%cz), dz_minus(grid%cz) !< half-sums for depth-step dz
real :: dRho_tmp_surf, dRho_tmp_bot
integer :: k !< counter
! ----------------------------------------------------------------------------
do k = 1, grid%cz - 1
dz_plus(k) = (grid%dz(k) + grid%dz(k + 1)) / 2.0
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.0
dz_minus(k) = (grid%dz(k)+grid%dz(k-1))/2.
end do
do k = 2, grid%cz - 1
N2(k) = - g/Rho_ref * (0.5 * &
((Rho(k+1) - Rho(k)) / dz_plus(k) + (Rho(k) - Rho(k-1)) / dz_minus(k)))
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
end do
dRho_tmp_bot = 1.0 / PrT_ri * (Rho_dyn_bot/ (kappa_ri * grid%dz(1)/2.0))
dRho_tmp_surf = 1.0 / PrT_ri * (Rho_dyn_surf/ (kappa_ri * grid%dz(grid%cz)/2.0))
N2(1) = - g/Rho_ref * dRho_tmp_bot
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
N2(grid%cz) = - g/Rho_ref * dRho_tmp_surf
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
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)
!< @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
......@@ -109,8 +147,10 @@ module obl_turb_closure
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
! ----------------------------------------------------------------------------
do k = 1, nz-1
dz_plus(k) = (dz(k)+dz(k+1))/2.
......@@ -134,84 +174,101 @@ module obl_turb_closure
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, S2_min, nz)
call limit_min_array(S2, c_S2_min, nz)
end subroutine
subroutine get_S2(S2, U, V, flux_u_surf, flux_u_bot, flux_v_surf, flux_v_bot, grid)
!< @brief calculate square of shear frequency
! ----------------------------------------------------------------------------
subroutine get_Ri_gradient(Ri_grad, N2, S2, grid)
!< @brief calculate Richardson gradient number
! ----------------------------------------------------------------------------
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 :: dz_plus(grid%cz), dz_minus(grid%cz) !< half-sums for depth-step dz
integer :: k !< counter
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)]
! ----------------------------------------------------------------------------
do k = 1, grid%cz - 1
dz_plus(k) = (grid%dz(k)+grid%dz(k+1))/2.
end do
call get_Ri_gradient_vec(Ri_grad, N2, S2, grid%cz)
end subroutine
do k = 2, grid%cz
dz_minus(k) = (grid%dz(k)+grid%dz(k-1))/2.
end do
! ----------------------------------------------------------------------------
subroutine get_Ri_gradient_vec(Ri_grad, N2, S2, 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)]
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
end do
integer :: k
! ----------------------------------------------------------------------------
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
do k = 1, cz
if (S2(k) < c_S2_min) then
Ri_grad(k) = N2(k) / c_S2_min
else
Ri_grad(k) = N2(k) / S2(k)
endif
enddo
end subroutine
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
! ----------------------------------------------------------------------------
subroutine get_TKE_production(P_TKE, Km, S2, grid)
!< @brief calculation of TKE production by shear
! ----------------------------------------------------------------------------
type (gridDataType), intent(in) :: grid
call limit_min_array(S2, S2_min, grid%cz)
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]
! ----------------------------------------------------------------------------
call get_TKE_production_vec(P_TKE, Km, S2, grid%cz)
end subroutine
subroutine get_Rig_vec(Ri_grad, N2, S2, nz)
!< @brief calculate Richardson gradient number
! ----------------------------------------------------------------------------
subroutine get_TKE_production_vec(P_TKE, Km, S2, cz)
!< @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]
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: Ri_grad(nz) !< Richardson gradient number
real, intent(in) :: N2(nz), S2(nz) !< square of buoyancy and shear frequency, [s**(-2)]
integer :: k !< counter
integer :: k
! ----------------------------------------------------------------------------
do k=1, nz
Ri_grad(k) = N2(k) / S2(k)
do k = 1, cz
P_TKE(k) = Km(k) * S2(k)
enddo
end subroutine
subroutine get_TKE_production_vec(P_TKE, Km, S2, nz)
!< @brief calculation of TKE production by shear
real, intent(in) :: S2(nz) !< shear frequency ^ 2, [s^-2]
real, intent(in) :: Km(nz) !< eddy viscosity for momentum, [m**2 / s]
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: P_TKE(nz) !< shear production of TKE, [m**2 / s**3]
integer :: k !< counter
! ----------------------------------------------------------------------------
subroutine get_TKE_buoyancy(B_TKE, Kh, N2, grid)
!< @brief calculation of TKE buoyancy production
! ----------------------------------------------------------------------------
type (gridDataType), intent(in) :: grid
do k = 1, nz
P_TKE(k) = Km(k) * S2(k)
enddo
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]
! ----------------------------------------------------------------------------
call get_TKE_buoyancy_vec(B_TKE, Kh, N2, grid%cz)
end subroutine
subroutine get_TKE_buoyancy_vec(B_TKE, Kh, N2, nz)
! ----------------------------------------------------------------------------
subroutine get_TKE_buoyancy_vec(B_TKE, Kh, N2, cz)
!< @brief calculation of TKE buoyancy production
real, intent(in) :: N2(nz) !< buoyancy frequency ^ 2, [s^-2]
real, intent(in) :: Kh(nz) !< eddy duffusivity for tracers, [m**2 / s]
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
integer :: k !< counter
! ----------------------------------------------------------------------------
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]
integer :: k
! ----------------------------------------------------------------------------
do k = 1, nz
do k = 1, cz
B_TKE(k) = - Kh(k) * N2(k)
enddo
end subroutine
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment