Skip to content
Snippets Groups Projects
Commit ab750afe authored by Daria Gladskikh's avatar Daria Gladskikh :ocean:
Browse files

preliminary structure of code with comments & rename cacl_obl -> obl_main

parent 222c05ca
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
obl_common_phys_parameters.f90
+ 4
− 4
View file @ ab750afe
module obl_common_phys_parameters
!< @brief common physical parameters and constants
implicit none
real, parameter :: g = 9.80655
real, parameter :: cp_water = 1003.5
real, parameter :: g = 9.80655 !< gravity acceleration [m / s^2]
real, parameter :: cp_water = 1003.5 !< specific heat water [J / (kg * K]
end module
\ No newline at end of file
obl_equation.f90
+ 30
− 28
View file @ ab750afe
module obl_equation
!< @brief solvers for general equations for scalars and vector
use obl_mathematics
implicit none
contains
subroutine solve_scalar_eq (Field, Kvisc, nz, dz, dt, flux_surf, flux_bot, f_right)
integer, intent (in) :: nz !Number of steps in z (depth)
real, intent (in) :: dt !Timestep
real, intent (in) :: dz(nz) !z-step
real :: Kvisc(nz) !Transfer coefficient
real :: G(nz) !Auxiliary variable
real :: a(nz), b(nz), c(nz), f(nz) !For tridiagonal matrix
real :: Field(nz) !Value in the point
integer :: k
real :: Kvisc_plus(nz), Kvisc_minus(nz) !Half-sum for K
real :: dz_plus(nz), dz_minus(nz) !Half-sum for dz
real :: f_right(nz) !hypothetic right-hand side
real :: gamma !Auxiliary variable
real :: flux_surf, flux_bot !fluxes
!< @brief solver for equation for scalar fields (temperature and salinity)
integer, intent (in) :: nz !< number of steps in z (depth), [m]
real, intent (in) :: dt !< time step [s]
real, intent (in) :: dz(nz) !< depth(z) step [m]
real :: Kvisc(nz) !< turbulent transfer coefficient [m^2/s]
real :: G(nz) !< explicit part
real :: a(nz), b(nz), c(nz), f(nz) !< coefficients for tridiagonal matrix
real :: Field(nz) !< value of unknown field in the point
integer :: k !< counter
real :: Kvisc_plus(nz), Kvisc_minus(nz) !<half-sums for turbulent transfer coefficient Kvisc
real :: dz_plus(nz), dz_minus(nz) !<half-sums for depth-step dz
real :: f_right(nz) !< right-hand side
real :: gamma !< variable determining the explicitness/implicitness of a scheme (from 0.0 to 1.0)
real :: flux_surf, flux_bot !fluxes!!!!!
gamma = 0.0
......@@ -71,20 +72,21 @@ module obl_equation
end subroutine solve_scalar_eq
subroutine solve_vector_eq (Field, Kvisc, nz, dz, dt, flux_surf, flux_bot, f_right)
integer, intent (in) :: nz !Number of steps in z (depth)
real, intent (in) :: dt !Timestep
real, intent (in) :: dz(nz) !z-step
real :: Kvisc(nz) !Transfer coefficient
real :: G(nz) !Auxiliary variable
real :: a(nz), b(nz), c(nz), f(nz) !For tridiagonal matrix
real :: Field(nz) !Value in the point
integer :: k
real :: Kvisc_plus(nz), Kvisc_minus(nz) !Half-sum for K
real :: dz_plus(nz), dz_minus(nz) !Half-sum for dz
real :: f_right(nz) !hypothetic right-hand side
real :: gamma !Auxiliary variable
real :: flux_surf, flux_bot !fluxes
!< @brief solver for equation for vectors (components of current velocity)
integer, intent (in) :: nz !< number of steps in z (depth), [m]
real, intent (in) :: dt !< time step [s]
real, intent (in) :: dz(nz) !< depth(z) step [m]
real :: Kvisc(nz) !< turbulent transfer coefficient [m^2/s]
real :: G(nz) !< explicit part
real :: a(nz), b(nz), c(nz), f(nz) !< coefficients for tridiagonal matrix
real :: Field(nz) !< value of unknown field in the point
integer :: k !< counter
real :: Kvisc_plus(nz), Kvisc_minus(nz) !<half-sums for turbulent transfer coefficient Kvisc
real :: dz_plus(nz), dz_minus(nz) !<half-sums for depth-step dz
real :: f_right(nz) !< right-hand side
real :: gamma !< variable determining the explicitness/implicitness of a scheme (from 0.0 to 1.0)
real :: flux_surf, flux_bot !fluxes!!!!!
gamma = 0.0
......
obl_forcing_and_boundary.f90
+ 66
− 39
View file @ ab750afe
module obl_forcing_and_boundary
!< @brief boundary conditions module
! modules used
use obl_common_phys_parameters
use obl_state_eq
! directives list
implicit none
contains
subroutine set_Flux_heat_bot (Flux_heat_bot, nt)
integer, intent(in) :: nt
real, intent(out) :: Flux_heat_bot(nt)
real :: Flux_heat_bot_tmp
integer :: i, status, num, set_Flux
!< @brief set bottom heat flux
integer, intent(in) :: nt !< number of timesteps
real, intent(out) :: Flux_heat_bot(nt) !< bottom heat flux, input: [W/m**2], output: [K*m/s] (?)
real :: Flux_heat_bot_tmp !< temporary bottom heat flux
integer :: i !< counter
integer :: status, num !< for file input
integer :: set_Flux !< how to set: 1 - const || 2 - function || 3 - file
set_Flux = 1
......@@ -20,13 +27,13 @@ module obl_forcing_and_boundary
num = 0
if (set_Flux.eq.1) then
do i = 1, nt
Flux_heat_bot(i) = 0 !flux_heat_bot_tmp;
Flux_heat_bot(i) = 0
end do
else if (set_Flux.eq.3) then
do while (status.eq.0)
read (444, *, iostat = status) Flux_heat_bot_tmp
num = num + 1
Flux_heat_bot(num) = Flux_heat_bot_tmp
Flux_heat_bot(num) = Flux_heat_bot_tmp/Rho_0/cp_water
if (num >= nt) then
exit
endif
......@@ -36,10 +43,14 @@ module obl_forcing_and_boundary
end subroutine
subroutine set_Flux_heat_surf (Flux_heat_surf, nt)
integer, intent(in) :: nt
real, intent(out) :: Flux_heat_surf(nt)
real :: Flux_heat_surf_tmp
integer :: i, status, num, set_Flux
!< @brief set surface heat flux
integer, intent(in) :: nt !< number of timesteps
real, intent(out) :: Flux_heat_surf(nt) !< surface heat flux, input: [W/m**2], output: [K*m/s] (?)
real :: Flux_heat_surf_tmp !< temporary surface heat flux
integer :: i !< counter
integer :: status, num !< for file input
integer :: set_Flux !< how to set: 1 - const || 2 - function || 3 - file
set_Flux = 1
......@@ -49,7 +60,7 @@ module obl_forcing_and_boundary
num = 0
if (set_Flux.eq.1) then
do i = 1, nt
Flux_heat_surf(i) = 0.0 !flux_heat_bot_tmp;
Flux_heat_surf(i) = 0.0
end do
else if (set_Flux.eq.3) then
do while (status.eq.0)
......@@ -65,11 +76,15 @@ module obl_forcing_and_boundary
end subroutine
subroutine set_Tau_x_surf(Tau_x_surf, flux_u_surf, nt)
!< @brief set surface u-momentum flux
integer, intent(in) :: nt
real, intent(out) :: Tau_x_surf(nt), flux_u_surf(nt)
real :: Tau_x_surf_tmp
integer :: i, status, num, set_Tau
integer, intent(in) :: nt !< number of timesteps
real, intent(out) :: Tau_x_surf(nt) !< [N/m**2]
real, intent(out) :: flux_u_surf(nt) !< [(m/s)**2] (?)
real :: Tau_x_surf_tmp !< temporary surface u-momentum flux
integer :: i !< counter
integer :: status, num !< for file input
integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file
set_Tau = 1
......@@ -99,17 +114,21 @@ module obl_forcing_and_boundary
close (188)
do i = 1, nt
flux_u_surf(i) = Tau_x_surf(i) / Rho_0
flux_u_surf(i) = Tau_x_surf(i) / Rho_0 ! sqrt (Tau/Rho) ??????
end do
end subroutine
subroutine set_Tau_y_surf(Tau_y_surf, flux_v_surf, nt)
!< @brief set surface v-momentum flux
integer, intent(in) :: nt
real, intent(out) :: Tau_y_surf(nt), flux_v_surf(nt)
real :: Tau_y_surf_tmp
integer :: i, status, num, set_Tau
integer, intent(in) :: nt !< number of timesteps
real, intent(out) :: Tau_y_surf(nt) !< [N/m**2]
real, intent(out) :: flux_v_surf(nt) !< [(m/s)**2] (?)
real :: Tau_y_surf_tmp !< temporary surface v-momentum flux
integer :: i !< counter
integer :: status, num !< for file input
integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file
set_Tau = 1
......@@ -145,11 +164,15 @@ module obl_forcing_and_boundary
subroutine set_Tau_x_bot(Tau_x_bot, flux_u_bot, nt)
!< @brief set bottom u-momentum flux
integer, intent(in) :: nt
real, intent(out) :: Tau_x_bot(nt), flux_u_bot(nt)
real :: Tau_x_bot_tmp
integer :: i, status, num, set_Tau
integer, intent(in) :: nt !< number of timesteps
real, intent(out) :: Tau_x_bot(nt) !< [N/m**2]
real, intent(out) :: flux_u_bot(nt) !< [(m/s)**2] (?)
real :: Tau_x_bot_tmp !< temporary bottom u-momentum flux
integer :: i !< counter
integer :: status, num !< for file input
integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file
set_Tau = 1
......@@ -158,17 +181,21 @@ module obl_forcing_and_boundary
end do
do i = 1, nt
flux_u_bot(i) = 0.0 !Tau_x_bot(i) / Rho_0
flux_u_bot(i) = Tau_x_bot(i) / Rho_0
end do
end subroutine
subroutine set_Tau_y_bot(Tau_y_bot, flux_v_bot, nt)
!< @brief set bottom v-momentum flux
integer, intent(in) :: nt
real, intent(out) :: Tau_y_bot(nt), flux_v_bot(nt)
real :: Tau_y_bot_tmp
integer :: i, status, num, set_Tau
integer, intent(in) :: nt !< number of timesteps
real, intent(out) :: Tau_y_bot(nt) !< [N/m**2]
real, intent(out) :: flux_v_bot(nt) !< [(m/s)**2] (?)
real :: Tau_y_bot_tmp !< temporary bottom v-momentum flux
integer :: i !< counter
integer :: status, num !< for file input
integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file
set_Tau = 1
......@@ -177,15 +204,16 @@ module obl_forcing_and_boundary
end do
do i = 1, nt
flux_v_bot(i) = 0.0 !Tau_y_bot(i) / Rho_0
flux_v_bot(i) = Tau_y_bot(i) / Rho_0
end do
end subroutine
subroutine set_Flux_sal_surf(Flux_sal_surf, nt)
!< @brief setsurface salinity flux
integer, intent(in) :: nt
real, intent(out) :: Flux_sal_surf(nt)
integer :: i
integer, intent(in) :: nt !< number of timesteps
real, intent(out) :: Flux_sal_surf(nt) ! [???]
integer :: i !< counter
do i = 1, nt
Flux_sal_surf(i) = 0.0
......@@ -193,16 +221,15 @@ module obl_forcing_and_boundary
end subroutine
subroutine set_Flux_sal_bot(Flux_sal_bot, nt)
!< @brief set bottom salinity flux
integer, intent(in) :: nt
real, intent(out) :: Flux_sal_bot(nt)
integer :: i
integer, intent(in) :: nt !< number of timesteps
real, intent(out) :: Flux_sal_bot(nt) ! [???]
integer :: i !< counter
do i = 1, nt
Flux_sal_bot(i) = 0.0
end do
end subroutine
end module
\ No newline at end of file
obl_initial_conditions.f90
+ 42
− 24
View file @ ab750afe
module obl_initial_conditions
!< @brief initial conditions module
! directives list
implicit none
real, parameter :: Theta_0 = 288.15
! reference temperature
real, parameter :: Theta_0 = 288.15 !< [K]
contains
subroutine Theta_init (Theta, nz, dz)
integer, intent(in) :: nz
real, intent(out) :: Theta(nz)
real, intent(in) :: dz(nz)
real :: Theta_start
integer :: k, status, num, init_Theta
!< @brief set initial conditions for temperature
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: Theta(nz) !< Theta [K]
real, intent(in) :: dz(nz) !< z-step [m]
real :: Theta_start !< temporary (for input from file)
integer :: k !< counter
integer :: status, num !< for file input
integer :: init_Theta !< how to set tau: 1 - const || 2 - function || 3 - file
init_Theta = 2
......@@ -22,7 +29,7 @@ module obl_initial_conditions
if (init_Theta.eq.1) then
do k = 1, nz
Theta(k) = 300.0 !Theta_start;
!Theta(k) = Theta(k) - Theta_0
Theta(k) = Theta(k) - Theta_0
end do
else if (init_Theta.eq.2) then
do k = 1, nz
......@@ -34,7 +41,7 @@ module obl_initial_conditions
read (222, *, iostat = status) Theta_start
num = num + 1
Theta(num) = Theta_start
!Theta(num) = Theta(num) - Theta_0
Theta(num) = Theta(num) - Theta_0
if (num >= nz) then
exit
endif
......@@ -44,14 +51,17 @@ module obl_initial_conditions
end subroutine
subroutine Salin_init (Salin, nz, dz)
!< @brief set initial conditions for salinity
integer, intent(in) :: nz
real, intent(out) :: Salin(nz)
real, intent(in) :: dz(nz)
real :: Salin_start
integer :: k, status, num, init_Salin
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: Salin(nz) !< Salinity [PSU]
real, intent(in) :: dz(nz) !< z-step [m]
real :: Salin_start !< temporary (for input from file)
integer :: k !< counter
integer :: status, num !< for file input
integer :: init_Salin !< how to set: 1 - const || 2 - function || 3 - file
init_Salin = 1 ! 1 - const, 2 - function(z), 3 - file
init_Salin = 1
open (777, file= 'Kato-Phillips/Salin_start.txt', status ='old') !Salin_PAPA.txt
!open (777, file= 'PAPA_06_2017/Salin_PAPA.txt', status ='old') !Salin_PAPA.txt
......@@ -80,12 +90,16 @@ module obl_initial_conditions
end subroutine
subroutine U_init (U, nz)
integer, intent(in) :: nz
real, intent(out) :: U(nz)
real :: U_start
integer :: k, status, num, init_U
!< @brief set initial conditions for U-velocity
init_U = 1 ! 1 - const, 2 - function(z), 3 - file
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: U(nz) !< U-velocity [m/s]
real :: U_start !< temporary (for input from file)
integer :: k !< counter
integer :: status, num !< for file input
integer :: init_U !< how to set: 1 - const || 2 - function || 3 - file
init_U = 1
open (888, file= 'Kato-Phillips/U_start.txt', status ='old')
status = 0
......@@ -112,12 +126,16 @@ module obl_initial_conditions
end subroutine
subroutine V_init (V, nz)
integer, intent(in) :: nz
real, intent(out) :: V(nz)
real :: V_start
integer :: k, status, num, init_V
!< @brief set initial conditions for V-velocity
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: V(nz) !< V-velocity [m/s]
real :: V_start !< temporary (for input from file)
integer :: k !< counter
integer :: status, num !< for file input
integer :: init_V !< how to set: 1 - const || 2 - function || 3 - file
init_V = 1 ! 1 - const, 2 - function(z), 3 - file
init_V = 1
open (999, file= 'Kato-Phillips/V_start.txt', status ='old')
status = 0
......
obl_k_epsilon.f90
+ 127
− 98
View file @ ab750afe
module obl_k_epsilon
!< @brief k-epsilon scheme module
! modules used
use obl_richardson
use obl_common_phys_parameters
use obl_mathematics
! directives list
implicit none
real, parameter :: C_mu = 0.09
real, parameter :: PrT = 0.8 !0.8 !1.25
real, parameter :: eps_min = 0.0000001
real, parameter :: TKE_min = 0.0001
real, parameter :: kappa_k_eps = 0.4
real, parameter :: C1_eps = 1.44
real, parameter :: C2_eps = 1.92
real, parameter :: C3_eps = -0.40
real, parameter :: Sch_TKE = 1.0
real, parameter :: Sch_eps = 1.11 !1.11
real, parameter :: C_mu = 0.09 !< momentum stability function constant
real, parameter :: PrT = 0.8 !< turbulent Prandtl number
real, parameter :: eps_min = 0.0000001 !< minimal value for epsilon (TKE dissipation rate), [W/kg]
real, parameter :: TKE_min = 0.0001 !< minimal valur for TKE, [J/kg]
real, parameter :: kappa_k_eps = 0.4 !< von Karman constant
real, parameter :: C1_eps = 1.44 !< eps-eq.: production constant
real, parameter :: C2_eps = 1.92 !< eps-eq.: dissipation constant
real, parameter :: C3_eps = -0.40 !< eps-eq.: buoyancy constant (-0.4 stable stratification or 1.14 unstable stratification)
real, parameter :: Sch_TKE = 1.0 !< Schmidt number for TKE
real, parameter :: Sch_eps = 1.11 !1.11 !< Scmidt number for eps
contains
subroutine TKE_init (TKE, nz)
integer, intent(in) :: nz
real, intent(out) :: TKE(nz)
integer :: k
!< @brief TKE initialization
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: TKE(nz) !< TKE, [J/kg]
integer :: k !< counter
do k = 2, nz-1
TKE(k) = 1.0 / (C_mu)**(1.0/2.0) * 0.001 * 0.001
......@@ -30,10 +34,11 @@ module obl_k_epsilon
end subroutine
subroutine eps_init (eps, nz, z)
integer, intent(in) :: nz
real, intent(in) :: z
real, intent(out) :: eps(nz)
integer :: k
!< @brief TKE dissipation rate initialization
integer, intent(in) :: nz !< number of z-steps
real, intent(in) :: z !< depth, [m]
real, intent(out) :: eps(nz) !< TKE dissipation rate, [W/kg]
integer :: k !< counter
do k = 2, nz-1
eps(k) = 0.001 * 0.001 * 0.001 / (kappa_k_eps * z) !not full z, but current z - must be fixed
......@@ -41,11 +46,12 @@ module obl_k_epsilon
end subroutine
subroutine get_TKE_generation (P_TKE, Km, S2, nz)
real, intent(in) :: S2(nz)
real, intent(in) :: Km(nz)
integer, intent(in) :: nz
real, intent(out) :: P_TKE(nz)
integer :: k
!< @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
do k = 1, nz
P_TKE(k) = Km(k) * S2(k)
......@@ -53,11 +59,12 @@ module obl_k_epsilon
end subroutine
subroutine get_TKE_buoyancy (B_TKE, Kh, N2, nz)
real, intent(in) :: N2(nz)
real, intent(in) :: Kh(nz)
integer, intent(in) :: nz
real, intent(out) :: B_TKE(nz)
integer :: k
!< @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
do k = 1, nz
B_TKE(k) = - Kh(k) * N2(k)
......@@ -65,10 +72,11 @@ module obl_k_epsilon
end subroutine
subroutine TKE_bc(TKE, flux_u_surf, flux_v_surf, nz)
real, intent(in) :: flux_u_surf, flux_v_surf
integer, intent(in) :: nz
real, intent(out) :: TKE(nz)
real :: flux_uv_surf
!< @brief TKE boundary conditions
real, intent(in) :: flux_u_surf, flux_v_surf !< fluxes (???)
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: TKE(nz) !< TKE, [J/kg]
real :: flux_uv_surf !< (????)
flux_uv_surf = (flux_u_surf**(2.0) + flux_v_surf**(2.0))**(1.0/4.0)
......@@ -79,11 +87,12 @@ module obl_k_epsilon
subroutine eps_bc(eps, flux_u_surf, flux_v_surf, nz, dz)
real, intent(in) :: flux_u_surf, flux_v_surf
integer, intent(in) :: nz
real, intent(in) :: dz(nz)
real, intent(out) :: eps(nz)
real :: flux_uv_surf
!< @brief epsilon boundary conditions
real, intent(in) :: flux_u_surf, flux_v_surf !< fluxes (???)
integer, intent(in) :: nz !< number of z-steps
real, intent(in) :: dz(nz) !< z-step
real, intent(out) :: eps(nz) !< epsilon ((TKE dissipation rate), [W/kg]
real :: flux_uv_surf !> ????
flux_uv_surf = (flux_u_surf**(2.0) + flux_v_surf**(2.0))**(1.0/4.0)
......@@ -93,10 +102,11 @@ module obl_k_epsilon
end subroutine eps_bc
subroutine get_Km_k_eps (Km, TKE, eps, nz)
integer, intent(in) :: nz
!< @brief calculate momentum transfer coefficient in k-epsilon closure
integer, intent(in) :: nz !< number of z-steps
real, intent(inout) :: TKE(nz), eps(nz)
real, intent(out) :: Km(nz)
integer :: k
real, intent(out) :: Km(nz) !eddy viscosity for momentum, [m**2 / s]
integer :: k !< counter
do k = 1, nz
call limit_min(eps(k), eps_min)
......@@ -106,10 +116,11 @@ module obl_k_epsilon
end subroutine get_Km_k_eps
subroutine get_Kh_k_eps (Kh, Km, nz)
integer, intent(in) :: nz
real, intent(in) :: Km(nz)
real, intent(out) :: Kh(nz)
integer :: k
!< @brief calculate scalar transfer coefficient in k-epsilon closure
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
......@@ -117,6 +128,7 @@ module obl_k_epsilon
end subroutine get_Kh_k_eps
subroutine get_Km_TKE(Km_TKE, Km, nz)
!< @brief calculate Km accounting Schmidt TKE number
integer, intent(in) :: nz
real, intent(out) :: Km_TKE(nz)
real, intent(in) :: Km(nz)
......@@ -128,6 +140,7 @@ module obl_k_epsilon
end subroutine get_Km_TKE
subroutine get_Km_eps (Km_eps, Km, nz)
!< @brief calculate Km accounting Schmidt eps number
integer, intent(in) :: nz
real, intent(out) :: Km_eps(nz)
real, intent(in) :: Km(nz)
......@@ -139,19 +152,23 @@ module obl_k_epsilon
end subroutine get_Km_eps
subroutine solve_TKE_eq (Field, Kvisc, nz, dz, dt, P_TKE, B_TKE, eps)
integer, intent (in) :: nz !Number of steps in z (depth)
real, intent (in) :: dt !Timestep
real, intent (in) :: dz(nz) !z-step
real :: Kvisc(nz) !Transfer coefficient
real :: G(nz) !Auxiliary variable
real :: a(nz), b(nz), c(nz), f(nz) !For tridiagonal matrix
real :: Field(nz) !Value in the point
integer :: k
real :: Kvisc_plus(nz), Kvisc_minus(nz) !Half-sum for K
real :: dz_plus(nz), dz_minus(nz) !Half-sum for dz
real :: f_right(nz), P_TKE(nz), B_TKE(nz), eps(nz)
real :: gamma !Auxiliary variable
!< @brief solver for TKE equation (explicit)
integer, intent (in) :: nz !< number of steps in z (depth), [m]
real, intent (in) :: dt !< time step [s]
real, intent (in) :: dz(nz) !< depth(z) step [m]
real :: Kvisc(nz) !< turbulent transfer coefficient [m^2/s]
real :: G(nz) !< explicit part
real :: a(nz), b(nz), c(nz), f(nz) !< coefficients for tridiagonal matrix
real :: Field(nz) !< value of unknown field in the point
integer :: k !< counter
real :: Kvisc_plus(nz), Kvisc_minus(nz) !<half-sums for turbulent transfer coefficient Kvisc
real :: dz_plus(nz), dz_minus(nz) !<half-sums for depth-step dz
real :: f_right(nz) !< right-hand side
real :: gamma !< variable determining the explicitness/implicitness of a scheme (from 0.0 to 1.0)
real :: P_TKE(nz) !< shear production of TKE, [m**2 / s**3]
real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
real :: eps(nz) !< TKE dissipation rate, [W/kg]
gamma = 0.0
......@@ -202,20 +219,24 @@ module obl_k_epsilon
subroutine solve_eps_eq (Field, Kvisc, nz, dz, dt, P_TKE, B_TKE, TKE)
integer, intent (in) :: nz !Number of steps in z (depth)
real, intent (in) :: dt !Timestep
real, intent (in) :: dz(nz) !z-step
real :: Kvisc(nz) !Transfer coefficient
real :: G(nz) !Auxiliary variable
real :: a(nz), b(nz), c(nz), f(nz) !For tridiagonal matrix
real :: Field(nz) !Value in the point
integer :: k
real :: Kvisc_plus(nz), Kvisc_minus(nz) !Half-sum for K
real :: dz_plus(nz), dz_minus(nz) !Half-sum for dz
real :: f_right(nz), P_TKE(nz), B_TKE(nz), TKE(nz) !hypothetic right-hand side
real :: C1_eps, C2_eps, C3_eps
real :: gamma !Auxiliary variable
!< @brief solver for epsilon equation (explicit)
integer, intent (in) :: nz !< number of steps in z (depth), [m]
real, intent (in) :: dt !< time step [s]
real, intent (in) :: dz(nz) !< depth(z) step [m]
real :: Kvisc(nz) !< turbulent transfer coefficient [m^2/s]
real :: G(nz) !< explicit part
real :: a(nz), b(nz), c(nz), f(nz) !< coefficients for tridiagonal matrix
real :: Field(nz) !< value of unknown field in the point
integer :: k !< counter
real :: Kvisc_plus(nz), Kvisc_minus(nz) !<half-sums for turbulent transfer coefficient Kvisc
real :: dz_plus(nz), dz_minus(nz) !<half-sums for depth-step dz
real :: f_right(nz) !< right-hand side
real :: gamma !< variable determining the explicitness/implicitness of a scheme (from 0.0 to 1.0)
real :: P_TKE(nz) !< shear production of TKE, [m**2 / s**3]
real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
real :: eps(nz) !< TKE dissipation rate, [W/kg]
real :: TKE(nz) !< TKE, [J/kg]
gamma = 0.0
......@@ -265,19 +286,23 @@ module obl_k_epsilon
end subroutine solve_eps_eq
subroutine solve_TKE_eq_semiimplicit (Field, Kvisc, nz, dz, dt, P_TKE, B_TKE, eps)
integer, intent (in) :: nz !Number of steps in z (depth)
real, intent (in) :: dt !Timestep
real, intent (in) :: dz(nz) !z-step
real :: Kvisc(nz) !Transfer coefficient
real :: G(nz) !Auxiliary variable
real :: a(nz), b(nz), c(nz), f(nz) !For tridiagonal matrix
real :: Field(nz) !Value in the point
integer :: k
real :: Kvisc_plus(nz), Kvisc_minus(nz) !Half-sum for K
real :: dz_plus(nz), dz_minus(nz) !Half-sum for dz
real :: f_right(nz), P_TKE(nz), B_TKE(nz), eps(nz)
real :: gamma !Auxiliary variable
!< @brief solver for TKE equation (semimplicit)
integer, intent (in) :: nz !< number of steps in z (depth), [m]
real, intent (in) :: dt !< time step [s]
real, intent (in) :: dz(nz) !< depth(z) step [m]
real :: Kvisc(nz) !< turbulent transfer coefficient [m^2/s]
real :: G(nz) !< explicit part
real :: a(nz), b(nz), c(nz), f(nz) !< coefficients for tridiagonal matrix
real :: Field(nz) !< value of unknown field in the point
integer :: k !< counter
real :: Kvisc_plus(nz), Kvisc_minus(nz) !<half-sums for turbulent transfer coefficient Kvisc
real :: dz_plus(nz), dz_minus(nz) !<half-sums for depth-step dz
real :: f_right(nz) !< right-hand side
real :: gamma !< variable determining the explicitness/implicitness of a scheme (from 0.0 to 1.0)
real :: P_TKE(nz) !< shear production of TKE, [m**2 / s**3]
real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
real :: eps(nz) !< TKE dissipation rate, [W/kg]
gamma = 0.0
......@@ -329,20 +354,24 @@ module obl_k_epsilon
end subroutine solve_TKE_eq_semiimplicit
subroutine solve_eps_eq_semiimplicit (Field, Kvisc, nz, dz, dt, P_TKE, B_TKE, TKE)
integer, intent (in) :: nz !Number of steps in z (depth)
real, intent (in) :: dt !Timestep
real, intent (in) :: dz(nz) !z-step
real :: Kvisc(nz) !Transfer coefficient
real :: G(nz) !Auxiliary variable
real :: a(nz), b(nz), c(nz), f(nz) !For tridiagonal matrix
real :: Field(nz) !Value in the point
integer :: k
real :: Kvisc_plus(nz), Kvisc_minus(nz) !Half-sum for K
real :: dz_plus(nz), dz_minus(nz) !Half-sum for dz
real :: f_right(nz), P_TKE(nz), B_TKE(nz), TKE(nz) !hypothetic right-hand side
real :: C1_eps, C2_eps, C3_eps
real :: gamma !Auxiliary variable
!< @brief solver for epsilon equation (semiimplict)
integer, intent (in) :: nz !< number of steps in z (depth), [m]
real, intent (in) :: dt !< time step [s]
real, intent (in) :: dz(nz) !< depth(z) step [m]
real :: Kvisc(nz) !< turbulent transfer coefficient [m^2/s]
real :: G(nz) !< explicit part
real :: a(nz), b(nz), c(nz), f(nz) !< coefficients for tridiagonal matrix
real :: Field(nz) !< value of unknown field in the point
integer :: k !< counter
real :: Kvisc_plus(nz), Kvisc_minus(nz) !<half-sums for turbulent transfer coefficient Kvisc
real :: dz_plus(nz), dz_minus(nz) !<half-sums for depth-step dz
real :: f_right(nz) !< right-hand side
real :: gamma !< variable determining the explicitness/implicitness of a scheme (from 0.0 to 1.0)
real :: P_TKE(nz) !< shear production of TKE, [m**2 / s**3]
real :: B_TKE(nz) !< buoyancy production of TKE, [m**2 / s**3]
real :: eps(nz) !< TKE dissipation rate, [W/kg]
real :: TKE(nz) !< TKE, [J/kg]
gamma = 0.0
......
calculate_obl.f90 → obl_main.f90
+ 1
− 1
View file @ ab750afe
program calculate_obl
program obl_main
use obl_time_and_space
use obl_initial_conditions
......
obl_mathematics.f90
+ 13
− 5
View file @ ab750afe
module obl_mathematics
!< @brief general mathematics module
! directives list
implicit none
contains
subroutine limit_min_array(Value, Value_min, nz)
integer, intent(in) :: nz
real, intent(inout) :: Value(nz)
real, intent(in) :: Value_min
integer :: k
!< @brief limitation function (minimum) for arrays
integer, intent(in) :: nz !< number of z-steps
real, intent(inout) :: Value(nz) !< value on current z-step
real, intent(in) :: Value_min !< minimal value
integer :: k !< counter
do k=1, nz
if (Value(k) < Value_min) then
......@@ -18,6 +21,7 @@ module obl_mathematics
end subroutine limit_min_array
subroutine limit_max_array(Value, Value_max, nz)
!< @brief limitation function (maximum) for arrays
integer, intent(in) :: nz
real, intent(inout) :: Value(nz)
real, intent(in) :: Value_max
......@@ -31,6 +35,7 @@ module obl_mathematics
end subroutine limit_max_array
subroutine limit_min(Value, Value_min)
!< @brief limitation function (minimum) for separate values
real, intent(inout) :: Value
real, intent(in) :: Value_min
......@@ -40,6 +45,7 @@ module obl_mathematics
end subroutine limit_min
subroutine limit_max(Value, Value_max)
!< @brief limitation function (maximum) for separate values
real, intent(inout) :: Value
real, intent(in) :: Value_max
......@@ -49,6 +55,7 @@ module obl_mathematics
end subroutine limit_max
subroutine tma(Field,a,b,c,f,nz)
!< @brief solver for tridiagonal matrix (Thomas algorithm)
integer, intent(in) :: nz
real, intent(in) :: a(nz), b(nz), c(nz), f(nz)
......@@ -73,6 +80,7 @@ module obl_mathematics
end subroutine tma
subroutine find_nan(Value, nz, k_nan)
!< @brief find NaN in array using NaN definition
integer, intent(in) :: nz
real, intent(in) :: Value(nz)
integer, intent(out) :: k_nan
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment