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

scm dynamics implementation update

parent dc783f58
No related branches found
No related tags found
No related merge requests found
......@@ -28,7 +28,6 @@ set(SOURCES
obl_initial_conditions.f90
obl_forcing_and_boundary.f90
obl_boundary.f90
obl_rhs.f90
obl_turb_closure.f90
obl_diag.f90
obl_k_epsilon.f90
......
......@@ -18,7 +18,6 @@ program obl_main
use obl_output
use obl_initial_conditions
use obl_forcing_and_boundary
use obl_rhs
use obl_scm
use obl_turb_closure
use obl_diag
......@@ -71,8 +70,6 @@ program obl_main
real :: time_begin, time_end, time_current !< time for output
real, allocatable :: f_Heat_right(:), f_Sal_right(:), f_u_right(:), f_v_right(:) !< RHS
integer :: closure_mode !< closure type:
!1 - pacanowski-philander, 2 - pacanowski-philander+,
!3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
......@@ -224,11 +221,8 @@ program obl_main
! allocation
call allocate_state_vec(grid%cz)
allocate(f_Heat_right(1:grid%cz))
allocate(f_Sal_right(1:grid%cz))
allocate(f_u_right(1:grid%cz))
allocate(f_v_right(1:grid%cz))
! initialize scm
call init_scm_vec(grid%cz)
......@@ -281,12 +275,6 @@ program obl_main
call set_const_tforcing(tau_y_bot, 0.0)
! ----------------------------------------------------------------------------
!set right-hand side
call set_f_Heat_right (f_Heat_right, grid%cz)
call set_f_Sal_right (f_Sal_right, grid%cz)
call set_f_u_right (f_v_right, V, grid%cz)
call set_f_v_right (f_v_right, U, grid%cz)
!open (199, file= 'output_Daria/mld.txt', status ='replace')
!open (20, file= 'output_Daria/surf_temp.txt', status ='replace')
......@@ -401,15 +389,7 @@ program obl_main
!< advance dynamics
! ----------------------------------------------------------------------------
call solve_scalar_eq (Theta_dev, Kh, grid%cz, grid%dz, dt, &
bc%heat_fluxH, bc%heat_flux0, f_heat_right)
call solve_scalar_eq (Salin_dev, Kh, grid%cz, grid%dz, dt, &
bc%salin_fluxH, bc%salin_flux0, f_sal_right)
call solve_state_eq(Rho, Theta_dev, Salin_dev, grid%cz)
call solve_vector_eq (U, Km, grid%cz, grid%dz, dt, &
bc%u_momentum_fluxH, bc%u_momentum_flux0, f_u_right)
call solve_vector_eq (V, Km, grid%cz, grid%dz, dt, &
bc%v_momentum_fluxH, bc%v_momentum_flux0, f_v_right)
call advance_scm(bc, grid, dt)
! ----------------------------------------------------------------------------
!> advance time
......@@ -478,13 +458,11 @@ program obl_main
!> deallocation
!> deallocate state
call deallocate_state_vec
deallocate(f_Heat_right)
deallocate(f_Sal_right)
deallocate(f_u_right)
deallocate(f_v_right)
!> deallocate scm
call deallocate_scm_vec
!> removing time-dependent forcing data
call deallocate_tforcing(sensible_hflux_surf)
......
module obl_rhs
!< @brief calculate RHS for general equations
! modules used
implicit none
real, parameter :: f_cor = 0 !< Coriolis parameter
real, parameter :: Ug = 0 !< u-geostrofic сurrent, [m/s]
real, parameter :: Vg = 0 !< v-geostrofic сurrent, [m/s]
contains
subroutine set_f_Heat_right (f_Heat_right, nz)
!< @brief set RHS for Theta equation
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: f_Heat_right(nz) !< RHS
integer :: k !< counter
do k = 1, nz
f_Heat_right(k) = 0
end do
end subroutine
subroutine set_f_Sal_right (f_Sal_right, nz)
!< @brief set RHS for Salinity equation
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: f_Sal_right(nz) !< RHS
integer :: k !< counter
do k = 1, nz
f_Sal_right(k) = 0
end do
end subroutine
subroutine set_f_u_right (f_u_right, V, nz)
!< @brief set RHS for U equation
integer, intent(in) :: nz !< number of z-steps
real, intent(in) :: V(nz) !< V, [m/s]
real, intent(out) :: f_u_right(nz) !< RHS
integer :: k !< counter
do k = 1, nz
f_u_right(k) = - f_cor * (V(k) - Vg)
end do
end subroutine
subroutine set_f_v_right (f_v_right, U, nz)
!< @brief set RHS for V equation
integer, intent(in) :: nz !< number of z-steps
real, intent(in) :: U(nz) !< U, [m/s]
real, intent(out) :: f_v_right(nz) !< RHS
integer :: k !< counter
do k = 1, nz
f_v_right(k) = - f_cor * (U(k) - Vg)
end do
end subroutine
end module
\ No newline at end of file
......@@ -4,12 +4,173 @@ module obl_scm
! modules used
use obl_math
use obl_state
use obl_grid
use obl_boundary
implicit none
real, allocatable :: Theta_rhs(:)
real, allocatable :: Salin_rhs(:)
real, allocatable :: U_rhs(:), V_rhs(:)
real :: f !< Coriolis parameter [1/s]
real :: Ug !< u-geostrophic сurrent, [m/s]
real :: Vg !< v-geostrophic сurrent, [m/s]
contains
! --------------------------------------------------------------------------------
subroutine init_scm_vec(cz)
!> @brief initialize single-column model
! ----------------------------------------------------------------------------
integer, intent(in) :: cz
! ----------------------------------------------------------------------------
f = 0.0
Ug = 0.0
Vg = 0.0
allocate(Theta_rhs(cz))
allocate(Salin_rhs(cz))
allocate(U_rhs(cz), V_rhs(cz))
end subroutine init_scm_vec
! --------------------------------------------------------------------------------
subroutine deallocate_scm_vec
!> @brief free single-column model data
! ----------------------------------------------------------------------------
! ----------------------------------------------------------------------------
deallocate(Theta_rhs)
deallocate(Salin_rhs)
deallocate(U_rhs, V_rhs)
end subroutine deallocate_scm_vec
! --------------------------------------------------------------------------------
subroutine advance_scm(bc, grid, dt)
!> @brief single-column model time advancement
! ----------------------------------------------------------------------------
type(oblBcType), intent(in) :: bc
type(gridDataType), intent(in) :: grid
real, intent(in) :: dt
! ----------------------------------------------------------------------------
call advance_temperature_eq_scm(bc, grid, dt)
call advance_salinity_eq_scm(bc, grid, dt)
call advance_state_eq_scm(bc, grid, dt)
call advance_dynamics_eq_scm(bc, grid, dt)
end subroutine advance_scm
! --------------------------------------------------------------------------------
subroutine advance_temperature_eq_scm(bc, grid, dt)
!> @brief single-column model time advancement
! ----------------------------------------------------------------------------
type(oblBcType), intent(in) :: bc
type(gridDataType), intent(in) :: grid
real, intent(in) :: dt
! ----------------------------------------------------------------------------
call set_temperature_eq_rhs(Theta_rhs, grid%cz)
call solve_scalar_eq(Theta_dev, Kh, grid%cz, grid%dz, dt, &
bc%heat_fluxH, bc%heat_flux0, Theta_rhs)
end subroutine advance_temperature_eq_scm
! --------------------------------------------------------------------------------
subroutine advance_salinity_eq_scm(bc, grid, dt)
!> @brief single-column model time advancement
! ----------------------------------------------------------------------------
type(oblBcType), intent(in) :: bc
type(gridDataType), intent(in) :: grid
real, intent(in) :: dt
! ----------------------------------------------------------------------------
call set_salinity_eq_rhs(Salin_rhs, grid%cz)
call solve_scalar_eq (Salin_dev, Kh, grid%cz, grid%dz, dt, &
bc%salin_fluxH, bc%salin_flux0, Salin_rhs)
end subroutine advance_salinity_eq_scm
! --------------------------------------------------------------------------------
subroutine advance_state_eq_scm(bc, grid, dt)
!> @brief single-column model time advancement
! ----------------------------------------------------------------------------
type(oblBcType), intent(in) :: bc
type(gridDataType), intent(in) :: grid
real, intent(in) :: dt
! ----------------------------------------------------------------------------
call solve_state_eq(Rho, Theta_dev, Salin_dev, grid%cz)
end subroutine advance_state_eq_scm
! --------------------------------------------------------------------------------
subroutine advance_dynamics_eq_scm(bc, grid, dt)
!> @brief single-column model time advancement
! ----------------------------------------------------------------------------
type(oblBcType), intent(in) :: bc
type(gridDataType), intent(in) :: grid
real, intent(in) :: dt
! ----------------------------------------------------------------------------
call set_dynamics_eq_rhs(U_rhs, V_rhs, U, V, f, grid%cz)
call solve_vector_eq (U, Km, grid%cz, grid%dz, dt, &
bc%u_momentum_fluxH, bc%u_momentum_flux0, U_rhs)
call solve_vector_eq (V, Km, grid%cz, grid%dz, dt, &
bc%v_momentum_fluxH, bc%v_momentum_flux0, V_rhs)
end subroutine advance_dynamics_eq_scm
subroutine set_temperature_eq_rhs(Theta_rhs, nz)
!< @brief set RHS for Theta equation
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: Theta_rhs(nz) !< RHS
integer :: k !< counter
do k = 1, nz
Theta_rhs(k) = 0
end do
end subroutine
subroutine set_salinity_eq_rhs(Salin_rhs, nz)
!< @brief set RHS for Salinity equation
integer, intent(in) :: nz !< number of z-steps
real, intent(out) :: Salin_rhs(nz) !< RHS
integer :: k !< counter
do k = 1, nz
Salin_rhs(k) = 0
end do
end subroutine
subroutine set_dynamics_eq_rhs(U_rhs, V_rhs, U, V, f, nz)
!< @brief set RHS for U equation
integer, intent(in) :: nz !< number of z-steps
real, intent(in) :: U(nz) !< V, [m/s]
real, intent(in) :: V(nz) !< V, [m/s]
real, intent(in) :: f
real, intent(out) :: U_rhs(nz) !< RHS
real, intent(out) :: V_rhs(nz) !< RHS
integer :: k !< counter
do k = 1, nz
U_rhs(k) = - f * (V(k) - Vg)
V_rhs(k) = - f * (U(k) - Vg)
end do
end subroutine
! --------------------------------------------------------------------------------
subroutine solve_scalar_eq (Field, Kvisc, nz, dz, dt, flux_surf, flux_bot, f_right)
!< @brief solver for equation for scalar fields (temperature and salinity)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment