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

pacanowski-philander model implementation update

parent 01499f53
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 2
− 2
View file @ 8e5348ae
......@@ -31,8 +31,8 @@ set(SOURCES
obl_turb_closure.f90
obl_diag.f90
obl_k_epsilon.f90
obl_pacanowski_philander.f90
obl_pacanowski_philander_plus.f90
obl_pph.f90
obl_pph_dyn.f90
obl_scm.f90
obl_main.f90
io_metadata.f90
......
obl_main.f90
+ 8
− 7
View file @ 8e5348ae
......@@ -21,8 +21,8 @@ program obl_main
use obl_scm
use obl_turb_closure
use obl_diag
use obl_pacanowski_philander
use obl_pacanowski_philander_plus
use obl_pph
use obl_pph_dyn
use obl_k_epsilon
use obl_boundary
use obl_io_plt
......@@ -45,7 +45,8 @@ program obl_main
type(gridDataType) :: grid
! turbulence closure parameters
type(pacanowskiParamType) :: param_pacanowski
type(pphParamType) :: param_pph
type(pphDynParamType) :: param_pph_dyn
!< surface forcing
type(timeForcingDataType) :: sensible_hflux_surf, latent_hflux_surf !< heat fluxes, [W/m^2]
......@@ -371,13 +372,13 @@ program obl_main
call get_TKE_buoyancy_vec(TKE_buoyancy, Kh, N2, grid%cz)
if (closure_mode.eq.1) then
call get_eddy_viscosity(Km, Ri_grad, param_pacanowski, grid)
call get_eddy_diffusivity(Kh, Ri_grad, param_pacanowski, grid)
call get_eddy_viscosity(Km, Ri_grad, param_pph, grid)
call get_eddy_diffusivity(Kh, Ri_grad, param_pph, grid)
else if (closure_mode.eq.2) then
call get_Km_plus(Km, Ri_grad, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, mld, grid%cz, grid%dz, grid%height)
call get_Kh_plus(Kh, Km, grid%cz)
bc%u_momentum_fluxH, bc%v_momentum_fluxH, mld, param_pph_dyn, grid%cz, grid%dz, grid%height)
call get_Kh_plus(Kh, Km, param_pph_dyn, grid%cz)
else if (closure_mode.eq.4) then
call TKE_bc(TKE, &
......
obl_pacanowski_philander.f90 → obl_pph.f90
+ 29
− 23
View file @ 8e5348ae
module obl_pacanowski_philander
module obl_pph
!< @brief standard pacanowski-philander scheme
! modules used
......@@ -10,12 +10,12 @@ module obl_pacanowski_philander
! public interface
! --------------------------------------------------------------------------------
public :: get_eddy_viscosity_vec, get_eddy_diffusivity_vec
public :: get_eddy_viscosity, get_eddy_diffusivity
public :: get_eddy_viscosity_vec, get_eddy_diffusivity_vec
! --------------------------------------------------------------------------------
!> @brief Pacanowski-Philander parameters
type, public :: pacanowskiParamType
type, public :: pphParamType
real :: Km_0 = 2.0 * 0.01 !< neutral eddy viscosity, [m**2 / s] *: INMCM: 7.0 * 0.01
real :: Kh_0 = 2.0 * 0.01 !< neutral eddy diffusivity, [m**2 / s], *: INMCM: 5.0 * 0.01
real :: alpha = 5.0 !< constant
......@@ -26,60 +26,66 @@ module obl_pacanowski_philander
contains
! --------------------------------------------------------------------------------
subroutine get_eddy_viscosity(Km, Ri_grad, param, grid)
!< @brief calculate eddy viscosity on grid
! ----------------------------------------------------------------------------
type(pphParamType), intent(in) :: param
type (gridDataType), intent(in) :: grid
real, dimension(grid%cz), intent(in) :: Ri_grad !< Richardson gradient number
real, dimension(grid%cz), intent(out) :: Km !< eddy viscosity, [m**2 / s]
! --------------------------------------------------------------------------------
call get_eddy_viscosity_vec(Km, Ri_grad, param, grid%cz)
end subroutine
subroutine get_eddy_viscosity_vec(Km, Ri_grad, param, n)
!< @brief calculate eddy viscosity
! ----------------------------------------------------------------------------
type(pacanowskiParamType), intent(in) :: param
type(pphParamType), intent(in) :: param
integer, intent(in) :: n !< vector size
real, dimension(n), intent(in) :: Ri_grad !< Richardson gradient number
real, dimension(n), intent(out) :: Km !< eddy viscosity, [m**2 / s]
integer :: k
! --------------------------------------------------------------------------------
do k = 1, n
Km(k) = param%Km_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n) + param%Km_b
end do
end subroutine
subroutine get_eddy_viscosity(Km, Ri_grad, param, grid)
!< @brief calculate eddy viscosity on grid
! --------------------------------------------------------------------------------
subroutine get_eddy_diffusivity(Kh, Ri_grad, param, grid)
!< @brief calculate eddy diffusivity on grid
! ----------------------------------------------------------------------------
type(pacanowskiParamType), intent(in) :: param
type(pphParamType), intent(in) :: param
type (gridDataType), intent(in) :: grid
real, dimension(grid%cz), intent(in) :: Ri_grad !< Richardson gradient number
real, dimension(grid%cz), intent(out) :: Km !< eddy viscosity, [m**2 / s]
real, dimension(grid%cz), intent(out) :: Kh !< eddy diffusivity, [m**2 / s]
! --------------------------------------------------------------------------------
call get_eddy_viscosity_vec(Km, Ri_grad, param, grid%cz)
call get_eddy_diffusivity_vec(Kh, Ri_grad, param, grid%cz)
end subroutine
subroutine get_eddy_diffusivity_vec(Kh, Ri_grad, param, n)
!< @brief calculate eddy diffusivity
! ----------------------------------------------------------------------------
type(pacanowskiParamType), intent(in) :: param
type(pphParamType), intent(in) :: param
integer, intent(in) :: n !< vector size
real, dimension(n), intent(in) :: Ri_grad !< Richardson gradient number
real, dimension(n), intent(out) :: Kh !< eddy diffusivity, [m**2 / s]
integer :: k
! --------------------------------------------------------------------------------
do k = 1, n
Kh(k) = param%Kh_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n + 1.0) + param%Kh_b
end do
end subroutine
subroutine get_eddy_diffusivity(Kh, Ri_grad, param, grid)
!< @brief calculate eddy diffusivity on grid
! ----------------------------------------------------------------------------
type(pacanowskiParamType), intent(in) :: param
type (gridDataType), intent(in) :: grid
real, dimension(grid%cz), intent(in) :: Ri_grad !< Richardson gradient number
real, dimension(grid%cz), intent(out) :: Kh !< eddy diffusivity, [m**2 / s]
call get_eddy_diffusivity_vec(Kh, Ri_grad, param, grid%cz)
end subroutine
end module
obl_pacanowski_philander_plus.f90 → obl_pph_dyn.f90
+ 21
− 13
View file @ 8e5348ae
module obl_pacanowski_philander_plus
module obl_pph_dyn
!< @brief modified pacanowski-philander scheme
! modules used
......@@ -8,15 +8,22 @@ module obl_pacanowski_philander_plus
! directives list
implicit none
real, parameter :: alpha_turb = 5.0 !< constant
real, parameter :: n_turb = 2.0 !< constant
real, parameter :: Km_b = 0.0001 !< background eddy viscosity for momentum, [m**2 / s] - 5 * 10**(-6) in inmcm
real, parameter :: PrT = 0.8 !1.25 !< turbulent Prandtl number
real, parameter :: kappa_pph_plus = 0.4 !< von Karman constant
!> @brief Pacanowski-Philander dynamic parameters
type, public :: pphDynParamType
real :: alpha = 5.0 !< constant
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
end type
contains
subroutine get_Km_plus(Km, Ri_grad, flux_u_surf, flux_v_surf, mld, nz, dz, z)
subroutine get_Km_plus(Km, Ri_grad, flux_u_surf, flux_v_surf, mld, param, nz, dz, z)
!< @brief calculate momentum transfer coefficient in P-Ph+ closure
type(pphDynParamType), intent(in) :: param
real, intent(out) :: Km(nz) !< eddy viscosity for momentum, [m**2 / s]
integer, intent(in) :: nz !< number of z-steps
real, intent(in) :: Ri_grad(nz) !< Richardson gradient number
......@@ -39,22 +46,23 @@ module obl_pacanowski_philander_plus
if (mld < mld_min) mld_tmp = mld_min
if (mld > mld_max) mld_tmp = mld_max
Km_0 = flux_uv_surf * kappa_pph_plus * mld_tmp / 2.0
Km_0 = flux_uv_surf * param%kappa * mld_tmp / 2.0
do k = 1, nz
Km(k) = Km_0 / (1.0 + alpha_turb * Ri_grad(k))**(n_turb) + Km_b
Km(k) = Km_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n) + param%Km_b
end do
end subroutine
subroutine get_Kh_plus(Kh, Km, nz)
subroutine get_Kh_plus(Kh, Km, param, nz)
!< @brief calculate scalar transfer coefficient in in P-Ph+ closure
type(pphDynParamType), intent(in) :: param
integer, intent(in) :: nz !< number of z-steps
real, intent(in) :: Km(nz) !< eddy viscosity for momentum, [m**2 / s]
real, intent(out) :: Kh(nz) !eddy diffusivity for tracers, [m**2 / s]
integer :: k !< counter
do k = 1, nz
Kh(k) = Km(k) / PrT
Kh(k) = Km(k) / param%PrT
end do
end subroutine
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment