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

grid and time setup update

parent 042f08b2
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 0
− 1
View file @ 0b9dfbef
......@@ -22,7 +22,6 @@ set(SOURCES
obl_tforcing.f90
obl_state.f90
obl_output.f90
obl_time_and_space.f90
obl_common_phys_parameters.f90
obl_math.f90
obl_state_eq.f90
......
obl_main.f90
+ 78
− 70
View file @ 0b9dfbef
......@@ -62,12 +62,11 @@ program obl_main
!< boundary conditions data
type(oblBcType) :: bc
real :: t !< time, [s]
integer :: nt !< number of time steps
real :: z !< depth, [m]
integer :: nz !< number of steps in z (depth), [m]
real :: domain_height
integer :: grid_cz
real :: dt !< time step [s]
real, allocatable :: dz(:) !< depth(z) step [m]
integer :: i, k !< counters
integer :: status, num !< for file input/output
real :: time_begin, time_end, time_current !< time for output
......@@ -186,46 +185,59 @@ program obl_main
endif
enddo
#ifndef USE_CONFIG_PARSER
! set time, depth, steps and number of space
call set_Time_Space(t, nt, z, nz)
#endif
allocate(dz(1:nz))
call set_dz_dt(dz, dt, nz, nt, t, z)
write (*,*) t, dt, nt, z, nz
!< setting model time
! ----------------------------------------------------------------------------
if (obl_setup == 1) then
time_begin = 0.0
time_end = 300.0 * 3600.0
time_current = time_begin
dt = 1.0
domain_height = 100.0
grid_cz = 32
endif
if (obl_setup == 2) then
time_begin = 0.0
time_end = 400.0 * 3600.0
time_current = time_begin
dt = 1.0
endif
! ----------------------------------------------------------------------------
! setting grid -- just an example
! > [zpos, height, cz, gcz]
call set_uniform_grid(grid, 0.0, z, nz)
call set_uniform_grid(grid, 0.0, domain_height, grid_cz)
! debug grid print
call print_grid(grid)
! allocation
call allocate_state_vec(nz)
call allocate_state_vec(grid%cz)
allocate(f_Heat_right(1:nz))
allocate(f_Sal_right(1:nz))
allocate(f_u_right(1:nz))
allocate(f_v_right(1:nz))
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))
allocate(Km_TKE(1:nz))
allocate(Km_eps(1:nz))
allocate(Km_TKE(1:grid%cz))
allocate(Km_eps(1:grid%cz))
! initialization of main fields
call Theta_init (Theta, Theta_dev, nz, dz)
call Theta_init (Theta, Theta_dev, grid%cz, grid%dz)
!call Theta_init (Theta, nz, dz)
call Salin_init(Salin, Salin_dev, nz, dz)
call U_init(U, nz)
call V_init(V, nz)
call Salin_init(Salin, Salin_dev, grid%cz, grid%dz)
call U_init(U, grid%cz)
call V_init(V, grid%cz)
!initialization of TKE & eps in case of k-epsilon closure
if (closure_mode.eq.3 .or. closure_mode.eq.4) then
call TKE_init(TKE, nz)
call eps_init(EPS, nz, z)
call TKE_init(TKE, grid%cz)
call eps_init(EPS, grid%cz, grid%height)
endif
!< setting atmospheric forcing
......@@ -262,10 +274,10 @@ program obl_main
! ----------------------------------------------------------------------------
!set right-hand side
call set_f_Heat_right (f_Heat_right, nz)
call set_f_Sal_right (f_Sal_right, nz)
call set_f_u_right (f_v_right, V, nz)
call set_f_v_right (f_v_right, U, nz)
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')
......@@ -274,9 +286,7 @@ program obl_main
status = 0
num = 0
time_begin = 0.0
time_end = t
time_current = time_begin
N2_0 = 0.0004
......@@ -340,62 +350,62 @@ program obl_main
bc%u_momentum_flux0, bc%v_momentum_flux0, bc%heat_flux0, bc%salin_flux0)
! ----------------------------------------------------------------------------
call solve_state_eq(Rho, Theta_dev, Salin_dev, nz)
call solve_state_eq(Rho, Theta_dev, Salin_dev, grid%cz)
#ifndef OBL_USE_GRID_DATATYPE
call get_N2(N2, Rho, bc%Rho_dynH, bc%rho_dyn0, dz, nz)
call get_N2(N2, Rho, bc%Rho_dynH, bc%rho_dyn0, grid%dz, grid%cz)
#else
call get_N2_on_grid(N2, Rho, bc%Rho_dynH, bc%Rho_dyn0, grid)
#endif
#ifndef OBL_USE_GRID_DATATYPE
call get_S2(S2, U, V, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%v_momentum_flux0, dz, nz)
bc%u_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%v_momentum_flux0, grid%dz, grid%cz)
#else
call get_S2_on_grid(S2, U, V, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%u_momentum_flux0, grid)
#endif
call get_Rig(Ri_grad, N2, S2, nz)
call get_Rig(Ri_grad, N2, S2, 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_TKE_generation(TKE_production, Km, S2, nz)
call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, nz)
call get_TKE_generation(TKE_production, Km, S2, grid%cz)
call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid%cz)
else if (closure_mode.eq.2) then
call get_Km_plus(Km, Ri_grad, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, mld, nz, dz, z)
call get_Kh_plus(Kh, Km, nz)
call get_TKE_generation(TKE_production, Km, S2, nz)
call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, nz)
bc%u_momentum_fluxH, bc%v_momentum_fluxH, mld, grid%cz, grid%dz, grid%height)
call get_Kh_plus(Kh, Km, grid%cz)
call get_TKE_generation(TKE_production, Km, S2, grid%cz)
call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid%cz)
else if (closure_mode.eq.4) then
call TKE_bc(TKE, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, nz)
bc%u_momentum_fluxH, bc%v_momentum_fluxH, grid%cz)
call eps_bc(EPS, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, nz, dz)
call get_Km_k_eps (Km, TKE, EPS, nz)
call get_Kh_k_eps (Kh, Km, nz)
call get_TKE_generation(TKE_production, Km, S2, nz)
call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, nz)
bc%u_momentum_fluxH, bc%v_momentum_fluxH, grid%cz, grid%dz)
call get_Km_k_eps (Km, TKE, EPS, grid%cz)
call get_Kh_k_eps (Kh, Km, grid%cz)
call get_TKE_generation(TKE_production, Km, S2, grid%cz)
call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid%cz)
call TKE_bc(TKE, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, nz)
bc%u_momentum_fluxH, bc%v_momentum_fluxH, grid%cz)
call eps_bc(EPS, &
bc%u_momentum_fluxH, bc%v_momentum_fluxH, nz, dz)
call get_Km_TKE(Km_TKE, Km, nz)
call get_Km_eps(Km_eps, Km, nz)
!call solve_TKE_eq_semiimplicit (TKE, Km_TKE, nz, dz, dt, TKE_production, TKE_buoyancy, eps)
call solve_TKE_eq (TKE, Km_TKE, nz, dz, dt, TKE_production, TKE_buoyancy, EPS)
call limit_min_array(TKE, TKE_min, nz)
call solve_eps_eq_semiimplicit (EPS, Km_eps, nz, dz, dt, TKE_production, TKE_buoyancy, TKE)
call limit_min_array(EPS, eps_min, nz)
bc%u_momentum_fluxH, bc%v_momentum_fluxH, grid%cz, grid%dz)
call get_Km_TKE(Km_TKE, Km, grid%cz)
call get_Km_eps(Km_eps, Km, grid%cz)
!call solve_TKE_eq_semiimplicit (TKE, Km_TKE, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, eps)
call solve_TKE_eq (TKE, Km_TKE, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, EPS)
call limit_min_array(TKE, TKE_min, grid%cz)
call solve_eps_eq_semiimplicit (EPS, Km_eps, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, TKE)
call limit_min_array(EPS, eps_min, grid%cz)
endif
call solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, &
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, nz, dz, dt, &
call solve_scalar_eq (Salin_dev, Kh, grid%cz, grid%dz, dt, &
bc%salin_fluxH, bc%salin_flux0, f_sal_right)
call solve_vector_eq (U, Km, nz, dz, dt, &
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, nz, dz, dt, &
call solve_vector_eq (V, Km, grid%cz, grid%dz, dt, &
bc%v_momentum_fluxH, bc%v_momentum_flux0, f_v_right)
!> advance time
......@@ -404,8 +414,8 @@ program obl_main
!> advance screen output
if (mod(i, nscreen) == 0) then
call get_mld(mld, N2, dz, nz, z)
call get_lab_mld(lab_mld, bc%U_dynH, N2_0, time_current, z)
call get_mld(mld, N2, grid%dz, grid%cz, grid%height)
call get_lab_mld(lab_mld, bc%U_dynH, N2_0, time_current, grid%height)
write(*, '(a,g0)') ' mld = ', mld
write(*, '(a,g0)') ' Theta(surface) = ', Theta_dev(grid%cz) + Theta_ref
......@@ -416,10 +426,10 @@ program obl_main
!> advance file output
if (mod(i, noutput) == 0) then
call push_state_vec(nz)
call push_state_vec(grid%cz)
call get_mld(mld, N2, dz, nz, z)
call get_lab_mld(lab_mld, bc%U_dynH, N2_0, time_current, z)
call get_mld(mld, N2, grid%dz, grid%cz, grid%height)
call get_lab_mld(lab_mld, bc%U_dynH, N2_0, time_current, grid%height)
call push_value_to_tseries(output_mld, mld)
call push_value_to_tseries(output_mld_ref, lab_mld)
......@@ -467,8 +477,6 @@ program obl_main
!> deallocation
call deallocate_state_vec
deallocate(dz)
deallocate(f_Heat_right)
deallocate(f_Sal_right)
deallocate(f_u_right)
......
obl_time_and_space.f90 deleted 100644 → 0
+ 0
− 66
View file @ 042f08b2
module obl_time_and_space
!< @brief set time and space
! directives list
implicit none
contains
subroutine set_Time_Space(t, nt, z, nz)
!< @brief set time, depth and number of steps
real, intent(out) :: t !< time, [s]
integer, intent(out) :: nt !< number of time steps
real, intent(out) :: z !< depth, [m]
integer, intent(out) :: nz !< number of steps in z (depth)
integer :: status, num !< for file input
integer :: i !< counter
!open (1, file= 'PAPA_06_2017/time_space_for_calc_PAPA.txt', status ='old')
open (1, file= 'Kato-Phillips/time_space_for_calc.txt', status ='old')
status = 0
do while (status.eq.0)
read (1, *, iostat = status) t, nt, z, nz
end do
close (1)
end subroutine
subroutine set_dz_dt(dz, dt, nz, nt, t, z)
!< @brief set time and depth steps
integer, intent(in) :: nz, nt !< number of depth and time steps
real, intent(in) :: z !< depth, [m]
real, intent(in) :: t !< time, [s]
real, intent(out) :: dz(nz) !< depth(z) step [m]
real, intent(out) :: dt !< time step [s]
integer :: status, num !< for file input
real :: dz_start !< temporary dz, [m]
integer :: init_dz !< how to set: 1 - const || 3 - file
integer :: k !< counter
dt = t/nt
init_dz = 1 ! 1 - const, 2 - function(z), 3 - file
!open (33, file= 'Kato-Phillips/dz.txt', status ='old') !'dz_PAPA.txt'
open (33, file= 'PAPA_06_2017/dz_PAPA.txt', status ='old') !'dz_PAPA.txt'
status = 0
num = 0
if (init_dz.eq.1) then
do k = 1, nz
dz(k) = z / nz
end do
else if (init_dz.eq.3) then
do while (status.eq.0)
read (33, *, iostat = status) dz_start
num = num + 1
dz(num) = dz_start
if (num >= nz) then
exit
endif
end do
endif
close (33)
end subroutine
end module
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment