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

adding grid module

parent c15bc10e
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 1
− 0
View file @ 7b120d70
......@@ -12,6 +12,7 @@ enable_language(Fortran)
set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/modules/)
set(SOURCES
obl_grid.f90
obl_time_and_space.f90
obl_common_phys_parameters.f90
obl_mathematics.f90
......
obl_def.fi
+ 1
− 0
View file @ 7b120d70
......@@ -2,3 +2,4 @@
!#define OBL_EXCLUDE_NETCDF
!#define USE_CONFIG_PARSER
!#define OBL_USE_GRID_DATATYPE
obl_grid.f90 0 → 100644
+ 112
− 0
View file @ 7b120d70
module obl_grid
!< @brief obl grid def.
! --------------------------------------------------------------------------------
! modules used
! --------------------------------------------------------------------------------
! directives list
implicit none
private
! public interface
! --------------------------------------------------------------------------------
public :: set_uniform_grid
public :: deallocate_grid
public :: print_grid
! --------------------------------------------------------------------------------
!> @brief grid datatype
type, public :: gridDataType
real :: zpos, height !< domain in [zpos, zpos + height]
integer :: cz !< number of cells
real, allocatable :: z(:) !< grid cell center coordinates [m]
real, allocatable :: dz(:) !< grid cell sizes [m]
end type
contains
! --------------------------------------------------------------------------------
subroutine set_uniform_grid(grid, zpos, height, cz)
!> @brief setting uniform grid
! ----------------------------------------------------------------------------
type (gridDataType), intent(inout) :: grid
real, intent(in) :: zpos, height
integer, intent(in) :: cz
integer :: k
real :: dz
real, dimension(cz + 1) :: edge
! ----------------------------------------------------------------------------
dz = height / cz
do k = 1, cz + 1
edge(k) = zpos + (k - 1) * dz
end do
call init_grid(grid, edge, zpos, height, cz)
end subroutine set_uniform_grid
! --------------------------------------------------------------------------------
subroutine init_grid(grid, edge, zpos, height, cz)
!> @brief init grid data
! ----------------------------------------------------------------------------
type (gridDataType), intent(inout) :: grid
real, intent(in) :: zpos, height
integer, intent(in) :: cz
real, dimension(cz + 1) :: edge
integer :: k
! ----------------------------------------------------------------------------
grid%zpos = zpos
grid%height = height
grid%cz = cz
allocate(grid%z(grid%cz))
allocate(grid%dz(grid%cz))
do k = 1, grid%cz
grid%z(k) = 0.5 * (edge(k) + edge(k + 1))
grid%dz(k) = edge(k + 1) - edge(k)
end do
end subroutine init_grid
subroutine deallocate_grid(grid)
!> @brief free grid data
! ----------------------------------------------------------------------------
type (gridDataType), intent(inout) :: grid
! ----------------------------------------------------------------------------
if (grid%cz > 0) then
deallocate(grid%z)
deallocate(grid%dz)
end if
end subroutine deallocate_grid
subroutine print_grid(grid)
!> @brief debug print
! ----------------------------------------------------------------------------
type (gridDataType), intent(in) :: grid
integer :: k
! ----------------------------------------------------------------------------
write(*, '(a,g0)') ' zpos = ', grid%zpos
write(*, '(a,g0)') ' height = ', grid%height
write(*, '(a,g0)') ' cz = ', grid%cz
do k = 1, grid%cz
write(*, '(a,g0)') ' z = ', grid%z(k)
end do
end subroutine print_grid
end module
\ No newline at end of file
obl_main.f90
+ 35
− 5
View file @ 7b120d70
......@@ -9,6 +9,7 @@ program obl_main
#endif
! modules used
use obl_grid
use obl_time_and_space
use obl_initial_conditions
use obl_forcing_and_boundary
......@@ -34,6 +35,8 @@ program obl_main
! directives list
implicit none
type (gridDataType) :: grid
real :: t !< time, [s]
integer :: nt !< number of time steps
integer :: nf !< number of time steps for forcing
......@@ -97,10 +100,10 @@ program obl_main
character(len = 128), parameter :: arg_key_config = "--config"
integer :: ierr
real :: value
! --------------------------------------------------------------------------------
! screen output parameters
integer, parameter :: nscreen = 1000
! arrays for output
real, allocatable :: Theta_write(:,:)
......@@ -199,6 +202,13 @@ program obl_main
call set_Time_Space_Forcing(nf, df)
write (*,*) t, dt, nt, z, nz, nf, df
! setting grid -- just an example
! > [zpos, height, cz, gcz]
call set_uniform_grid(grid, 0.0, z, nz)
! debug grid print
call print_grid(grid)
! allocation
allocate(Theta(1:nz))
allocate(Theta_dev(1:nz))
......@@ -317,8 +327,16 @@ program obl_main
call get_rho_dyn(Rho_dyn_surf, flux_u_surf_res, flux_v_surf_res, flux_heat_surf_res, flux_sal_surf_res)
call get_rho_dyn(Rho_dyn_bot, flux_u_bot_res, flux_v_bot_res, flux_heat_bot_res, flux_sal_bot_res)
call solve_state_eq (Rho, Theta_dev, Salin_dev, nz)
#ifndef OBL_USE_GRID_DATATYPE
call get_N2(N2, Rho, Rho_dyn_surf, Rho_dyn_bot, dz, nz)
#else
call get_N2_on_grid(N2, Rho, Rho_dyn_surf, Rho_dyn_bot, grid)
#endif
#ifndef OBL_USE_GRID_DATATYPE
call get_S2(S2, U, V, flux_u_surf_res, flux_v_surf_res, flux_u_bot_res, flux_v_bot_res, dz, nz)
#else
call get_S2_on_grid(S2, U, V, flux_u_surf_res, flux_v_surf_res, flux_u_bot_res, flux_v_bot_res, grid)
#endif
call get_Rig(Ri_grad, N2, S2, nz)
call get_hml (hml, maxcellnumber, N2, dz, nz, z)
call get_lab_hml(lab_hml, flux_u_surf_res, flux_v_surf_res, time_current, z)
......@@ -383,8 +401,17 @@ program obl_main
!write(199,*) time_current, hml, lab_hml
!write(20,*) time_current, Theta_dev(nz) + Theta_ref
!> advancing time
i = i + 1
time_current = time_current + dt
if (mod(i, nscreen) == 0) then
write(*, '(a,g0)') ' mld = ', hml
write(*, '(a,g0)') ' Theta(surface) = ', Theta_dev(grid%cz) + Theta_ref
write(*, '(a,g0,a,g0,a)') ' current time = ', time_current / 3600.0, ' HRS [ ', &
(time_current / time_end) * 100.0, '% ]'
write(*, '(a)') '-------------------------------------------------'
endif
enddo
Theta_C_write = Theta_write - 273.15
......@@ -494,5 +521,8 @@ program obl_main
deallocate(Theta_C_write)
deallocate(lab_hml_write)
! > removing grid data
call deallocate_grid(grid)
end program
\ No newline at end of file
obl_pacanowski_philander.f90
+ 1
− 0
View file @ 7b120d70
......@@ -15,6 +15,7 @@ module obl_pacanowski_philander
real, parameter :: Kh_b = 0.00001 !< background eddy diffusivity for tracers, [m**2 / s]
contains
subroutine get_Km(Km, Ri_grad, nz)
!< @brief calculate momentum transfer coefficient in P-Ph closure
integer, intent(in) :: nz !< number of z-steps
......
obl_richardson.f90
+ 75
− 0
View file @ 7b120d70
......@@ -2,6 +2,7 @@ module obl_richardson
!< @brief calculate Richardson gradient number
! modules used
use obl_grid
use obl_mathematics
use obl_common_phys_parameters
use obl_state_eq
......@@ -52,6 +53,41 @@ module obl_richardson
end subroutine
subroutine get_N2_on_grid(N2, Rho, Rho_dyn_surf, Rho_dyn_bot, grid)
!< @brief calculate square of buoyancy 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 :: 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
end do
do k = 2, grid%cz
dz_minus(k) = (grid%dz(k) + grid%dz(k - 1)) / 2.0
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)))
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
N2(grid%cz) = - g/Rho_ref * dRho_tmp_surf
end subroutine
subroutine get_S2(S2, U, V, flux_u_surf, flux_u_bot, flux_v_surf, flux_v_bot, dz, nz)
!< @brief calculate square of shear frequency
......@@ -91,6 +127,45 @@ module obl_richardson
end subroutine
subroutine get_S2_on_grid(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) :: 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
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
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
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
call limit_min_array(S2, S2_min, grid%cz)
end subroutine
subroutine get_Rig(Ri_grad, N2, S2, nz)
!< @brief calculate Richardson gradient number
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment