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

code documentation

parent 50c5d069
No related branches found
No related tags found
No related merge requests found
module obl_equation
!< @brief solvers for general equations for scalars and vector
! modules used
use obl_mathematics
implicit none
......@@ -20,7 +22,7 @@ module obl_equation
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!!!!!
real :: flux_surf, flux_bot !< fluxes
gamma = 0.0
......@@ -86,7 +88,7 @@ module obl_equation
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!!!!!
real :: flux_surf, flux_bot !< fluxes
gamma = 0.0
......
program obl_main
!< @brief main program for calculations for ocean boundary layer
! modules used
use obl_time_and_space
use obl_initial_conditions
use obl_forcing_and_boundary
......@@ -9,51 +11,53 @@ program obl_main
use obl_pacanowski_philander
use obl_pacanowski_philander_plus
use obl_k_epsilon
use io
use io_metadata
!use vertical_mixing
use vermix
! directives list
implicit none
real :: t !Time
integer :: nt !Number of steps in t (time)
real :: z !Depth
integer :: nz !Number of steps in z (depth)
real :: dt !Timestep
real, allocatable :: dz(:) !z-step
integer :: i, k
integer :: status, num
real :: time_begin, time_end, time_current !time for output
real, allocatable :: Theta(:)
real, allocatable :: Salin(:)
real, allocatable :: U(:)
real, allocatable :: V(:)
real, allocatable :: Km(:)
real, allocatable :: Kh(:)
real, allocatable :: Flux_heat_surf(:), Flux_heat_bot(:)
real, allocatable :: Flux_sal_surf(:), Flux_sal_bot(:)
real, allocatable :: Tau_x_surf(:), Tau_y_surf(:), Tau_x_bot(:), Tau_y_bot(:)
real, allocatable :: flux_u_surf(:), flux_u_bot(:)
real, allocatable :: flux_v_surf(:), flux_v_bot(:)
real, allocatable :: f_Heat_right(:)
real, allocatable :: f_Sal_right(:)
real, allocatable :: f_u_right(:)
real, allocatable :: f_v_right(:)
integer :: closure_mode
real, allocatable :: Rho(:), N2(:), S2(:), Ri_grad(:)
real :: hml, lab_hml
integer :: maxcellnumber
real, allocatable :: TKE(:), eps(:), P_TKE(:), B_TKE(:)
real, allocatable :: Km_TKE(:), Km_eps(:)
real :: t !< time, [s]
integer :: nt !< number of time steps
real :: z !< depth, [m]
integer :: nz !< number of steps in z (depth), [m]
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
real, allocatable :: Theta(:) !< temperature, [K]
real, allocatable :: Salin(:) !< salinity, [PSU]
real, allocatable :: U(:) !< U-velocity [m/s]
real, allocatable :: V(:) !< V-velocity [m/s]
real, allocatable :: Km(:) !< eddy viscosity for momentum, [m**2 / s]
real, allocatable :: Kh(:) !eddy diffusivity for tracers, [m**2 / s]
real, allocatable :: Flux_heat_surf(:), Flux_heat_bot(:) !< heat fluxes, [K*m/s]
real, allocatable :: Flux_sal_surf(:), Flux_sal_bot(:) !< salt fluxes, [PSU*m/s] ???
real, allocatable :: Tau_x_surf(:), Tau_y_surf(:), Tau_x_bot(:), Tau_y_bot(:) !< [N/m**2]
real, allocatable :: flux_u_surf(:), flux_u_bot(:) !< [(m/s)**2]
real, allocatable :: flux_v_surf(:), flux_v_bot(:) !< [(m/s)**2]
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
real, allocatable :: Rho(:) !< density, [kg / m**3]
real, allocatable :: N2(:) !< square of buoyancy frequency, [s**(-2)]
real, allocatable :: S2(:) !< square of shear frequency, [s**(-2)]
real, allocatable :: Ri_grad(:) !< Richardson gradient number
real :: hml !< mixed layer depth, [m]
real :: lab_hml !< theoretical mixed layer depth, [m]
integer :: maxcellnumber !< index of cell with maximum N2
real, allocatable :: TKE(:) !< TKE, [J/kg]
real, allocatable :: eps(:) !< TKE dissipation rate, [W/kg]
real, allocatable :: P_TKE(:) !< shear production of TKE, [m**2 / s**3]
real, allocatable :: B_TKE(:) !< buoyancy production of TKE, [m**2 / s**3]
real, allocatable :: Km_TKE(:) !eddy viscosity for momentum adjusted for Schmidt TKE number, [m**2 / s]
real, allocatable :: Km_eps(:) !eddy viscosity for momentum adjusted for Schmidt eps number, [m**2 / s]
! arrays for output
real, allocatable :: Theta_write(:,:)
real, allocatable :: Salin_write(:,:)
real, allocatable :: U_write(:,:)
......@@ -71,14 +75,13 @@ program obl_main
closure_mode = 1 ! 1 - pacanowski-philander, 2 - pacanowski-philander+, 3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
! set time, depth, steps and number of space
call set_Time_Space(t, nt, z, nz)
allocate(dz(1:nz))
call set_dz_dt(dz, dt, nz, nt, t, z)
write (*,*) t, dt, nt, z, nz
! allocation
allocate(Theta(1:nz))
allocate(Salin(1:nz))
allocate(U(1:nz))
......@@ -95,29 +98,22 @@ program obl_main
allocate(flux_u_bot(1:nt))
allocate(flux_v_surf(1:nt))
allocate(flux_v_bot(1:nt))
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(Km(1:nz))
allocate(Kh(1:nz))
allocate(Rho(1:nz))
allocate(N2(1:nz))
allocate(S2(1:nz))
allocate(Ri_grad(1:nz)) !
allocate(TKE(1:nz)) !
allocate(eps(1:nz)) !
allocate(P_TKE(1:nz)) !
allocate(B_TKE(1:nz)) !
allocate(Km_TKE(1:nz)) !
allocate(Km_eps(1:nz)) !
allocate(Ri_grad(1:nz))
allocate(TKE(1:nz))
allocate(eps(1:nz))
allocate(P_TKE(1:nz))
allocate(B_TKE(1:nz))
allocate(Km_TKE(1:nz))
allocate(Km_eps(1:nz))
allocate(Theta_write(nz, nt))
allocate(Salin_write(nz, nt))
allocate(U_write(nz, nt))
......@@ -133,42 +129,28 @@ program obl_main
allocate(Theta_C_write(nz, nt))
allocate(lab_hml_write(nt))
! initialization of main fields
call Theta_init (Theta, nz, dz)
call Salin_init(Salin, nz, dz)
call U_init(U, nz)
call V_init(V, nz)
!set boundary conditions
call set_Flux_heat_surf (Flux_heat_surf, nt)
call set_Flux_heat_bot (Flux_heat_bot, nt)
call set_Flux_sal_surf (Flux_sal_surf, nt)
call set_Flux_sal_bot (Flux_sal_bot, nt)
call set_Tau_x_surf(Tau_x_surf, flux_u_surf, nt)
call set_Tau_y_surf(Tau_y_surf, flux_v_surf, nt)
call set_Tau_x_bot(Tau_x_bot, flux_u_bot, nt)
call set_Tau_y_bot(Tau_y_bot, flux_v_bot, nt)
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)
open (199, file= 'output_Daria/hml.txt', status ='replace')
open (20, file= 'output_Daria/surf_temp.txt', status ='replace')
status = 0
......@@ -178,20 +160,14 @@ program obl_main
time_end = t
time_current = time_begin
!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)
endif
do k = 1, nz
Km(k) = 10.0**(-1.0)
Kh(k) = 10.0**(-1.0)
end do
i=1
do while (time_current < time_end ) !.and. i<=nt
if (closure_mode.eq.1) then
call solve_state_eq (Rho, Theta, Salin, nz)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment