Newer
Older
program obl_main
!< @brief main program for calculations for ocean boundary layer
#ifdef USE_CONFIG_PARSER
use iso_c_binding, only: C_NULL_CHAR
use config_parser
#endif
use obl_time_and_space
use obl_initial_conditions
use obl_forcing_and_boundary
use obl_rhs
use obl_equation
use obl_pacanowski_philander
use obl_pacanowski_philander_plus
use obl_k_epsilon
use io
use io_metadata
#ifdef USE_SFX
use sfx_data, only: meteoDataType, sfxDataType
use sfx_most, only: get_surface_fluxes_most => get_surface_fluxes, &
#endif
#ifdef USE_CONFIG_PARSER
!use vertical_mixing, default = off
!use vermix
implicit none
! turbulence closure parameters
type(pacanowskiParamType) :: param_pacanowski
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 :: Flux_heat_surf(:), Flux_heat_bot(:) !< heat fluxes, [K*m/s]
real, allocatable :: Times_flux_Heat_surf(:), Times_flux_Heat_bot(:)
real, allocatable :: Flux_sal_surf(:), Flux_sal_bot(:) !< salt fluxes, [PSU*m/s] ???
real, allocatable :: Times_flux_Sal_surf(:), Times_flux_Sal_bot(:)
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 :: Times_flux_u_surf(:), Times_flux_u_bot(:)
real, allocatable :: flux_v_surf(:), flux_v_bot(:) !< [(m/s)**2]
real, allocatable :: Times_flux_v_surf(:), Times_flux_v_bot(:)
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
integer, parameter :: output_mode = 1 ! 1 -- netcdf, 2 -- ascii, 3 -- tecplot
real :: Rho_dyn_surf, Rho_dyn_bot !< density dynamic, [kg / m**3]
real :: mld !< mixed layer depth, [m]
real :: lab_mld !< theoretical mixed layer depth, [m]
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]
real :: flux_heat_surf_res, flux_heat_bot_res, flux_sal_surf_res, flux_sal_bot_res
real :: flux_u_surf_res, flux_u_bot_res, flux_v_surf_res, flux_v_bot_res
! command line arguments
! --------------------------------------------------------------------------------
integer :: num_args
character(len = 128) :: arg
character(len = 128), parameter :: arg_key_help = '--help'
character(len = 128), parameter :: arg_key_config = "--config"
integer :: ierr
! --------------------------------------------------------------------------------
! screen output parameters
integer, parameter :: nscreen = 1000
! file output parameters
integer, parameter :: noutput = 60
closure_mode = 4 !< 1 - pacanowski-philander, 2 - pacanowski-philander+, 3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
! --- command line arguments processing
num_args = command_argument_count()
do i = 1, num_args
call get_command_argument(i, arg)
if (trim(arg) == trim(arg_key_help)) then
write(*, *) ' obl model, usage:'
write(*, *) ' --help'
write(*, *) ' print usage options'
write(*, *) ' --config [filename]'
write(*, *) ' use configuration file'
return
end if
if (trim(arg) == trim(arg_key_config)) then
if (i == num_args) then
write(*, *) ' FAILURE! > missing configuration file [key] argument'
ierr = 1 ! signal ERROR
return
end if
call get_command_argument(i + 1, arg)
#ifdef USE_CONFIG_PARSER
call c_config_run(trim(arg)//C_NULL_CHAR, status)
if (status == 0) then
write(*, *) ' FAILURE! > unable to parse configuration file: ', trim(arg)
ierr = 1 ! signal ERROR
return
end if
call c_config_is_varname("domain.depth"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float("domain.depth"//C_NULL_CHAR, z, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
end if
call c_config_is_varname("grid.cz"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_int("grid.cz"//C_NULL_CHAR, nz, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
end if
call c_config_is_varname("time.end"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float("time.end"//C_NULL_CHAR, t, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
end if
call c_config_is_varname("time.n"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_int("time.n"//C_NULL_CHAR, nt, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
end if
#endif
endif
enddo
call set_Time_Space(t, nt, z, nz)
#endif
call set_dz_dt(dz, dt, nz, nt, t, z)
call set_Time_Space_Forcing(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)
allocate(Times_flux_Heat_surf(1:nf))
allocate(flux_Heat_bot(1:nf))
allocate(Times_flux_Heat_bot(1:nf))
allocate(flux_Sal_surf(1:nf))
allocate(Times_flux_Sal_surf(1:nf))
allocate(flux_Sal_bot(1:nf))
allocate(Times_flux_Sal_bot(1:nf))
allocate(Tau_x_surf(1:nf))
allocate(Tau_y_surf(1:nf))
allocate(Tau_x_bot(1:nf))
allocate(Tau_y_bot(1:nf))
allocate(flux_u_surf(1:nf))
allocate(Times_flux_u_surf(1:nf))
allocate(flux_u_bot(1:nf))
allocate(Times_flux_u_bot(1:nf))
allocate(flux_v_surf(1:nf))
allocate(Times_flux_v_surf(1:nf))
allocate(flux_v_bot(1:nf))
allocate(Times_flux_v_bot(1:nf))
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(TKE(1:nz))
allocate(eps(1:nz))
allocate(P_TKE(1:nz))
allocate(B_TKE(1:nz))
allocate(Km_TKE(1:nz))
call Theta_init (Theta, Theta_dev, nz, 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 set_Flux_heat_surf (Flux_heat_surf, Times_flux_heat_surf, nf, df)
call set_Flux_heat_bot (Flux_heat_bot, Times_flux_heat_bot, nf, df)
call set_Flux_sal_surf (Flux_sal_surf, Times_flux_sal_surf, nf, df)
call set_Flux_sal_bot (Flux_sal_bot, Times_flux_sal_bot, nf, df)
call set_Tau_x_surf(Tau_x_surf, flux_u_surf, Times_flux_u_surf, nf, df)
call set_Tau_y_surf(Tau_y_surf, flux_v_surf, Times_flux_v_surf, nf, df)
call set_Tau_x_bot(Tau_x_bot, flux_u_bot, Times_flux_u_bot, nf, df)
call set_Tau_y_bot(Tau_y_bot, flux_v_bot, Times_flux_v_bot, nf, df)
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 (20, file= 'output_Daria/surf_temp.txt', status ='replace')
status = 0
num = 0
time_begin = 0.0
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 while (time_current < time_end )
! ----------------------------------------------------------------------------
call get_value_interpolate(flux_heat_surf_res, time_current, flux_heat_surf, Times_flux_heat_surf, nf)
call get_value_interpolate(flux_heat_bot_res, time_current, flux_heat_bot, Times_flux_heat_bot, nf)
call get_value_interpolate(flux_sal_surf_res, time_current, flux_sal_surf, Times_flux_sal_surf, nf)
call get_value_interpolate(flux_sal_bot_res, time_current, flux_sal_bot, Times_flux_sal_bot, nf)
call get_value_interpolate(flux_u_surf_res, time_current, flux_u_surf, Times_flux_u_surf, nf)
call get_value_interpolate(flux_v_surf_res, time_current, flux_v_surf, Times_flux_v_surf, nf)
call get_value_interpolate(flux_u_bot_res, time_current, flux_u_bot, Times_flux_u_surf, nf)
call get_value_interpolate(flux_v_bot_res, time_current, flux_v_bot, Times_flux_v_bot, nf)
call get_dyn_velocity(u_dynH, flux_u_surf_res, flux_v_surf_res)
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 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_eddy_viscosity(Km, Ri_grad, param_pacanowski, grid)
call get_eddy_diffusivity(Kh, Ri_grad, param_pacanowski, grid)
call get_TKE_generation(P_TKE, Km, S2, nz)
call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
else if (closure_mode.eq.2) then
call get_Km_plus(Km, Ri_grad, flux_u_surf_res, flux_v_surf_res, mld, nz, dz, z)
call get_Kh_plus(Kh, Km, nz)
call get_TKE_generation(P_TKE, Km, S2, nz)
call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
else if (closure_mode.eq.4) then
call TKE_bc(TKE, flux_u_surf_res, flux_v_surf_res, nz)
call eps_bc(eps, flux_u_surf_res, flux_v_surf_res, nz, dz)
call get_Km_k_eps (Km, TKE, eps, nz)
call get_Kh_k_eps (Kh, Km, nz)
call get_TKE_generation(P_TKE, Km, S2, nz)
call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
call TKE_bc(TKE, flux_u_surf_res, flux_v_surf_res, nz)
call eps_bc(eps, flux_u_surf_res, flux_v_surf_res, 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, P_TKE, B_TKE, eps)
call solve_TKE_eq (TKE, Km_TKE, nz, dz, dt, P_TKE, B_TKE, eps)
call limit_min_array(TKE, TKE_min, nz)
call solve_eps_eq_semiimplicit (eps, Km_eps, nz, dz, dt, P_TKE, B_TKE, TKE)
call limit_min_array(eps, eps_min, nz)
endif
call solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, flux_heat_surf_res, flux_heat_bot_res, f_heat_right)
call solve_scalar_eq (Salin_dev, Kh, nz, dz, dt, flux_sal_surf_res, flux_sal_bot_res, f_sal_right)
call solve_vector_eq (U, Km, nz, dz, dt, flux_u_surf_res, flux_u_bot_res, f_u_right)
call solve_vector_eq (V, Km, nz, dz, dt, flux_v_surf_res, flux_v_bot_res, f_v_right)
!else if (closure_mode.eq.5) then
! call solve_state_eq (Rho, Theta_dev, Salin_dev, nz)
! call vermix_init(rho, u, v, Tau_x_surf(i), Tau_y_surf(i), Ri_grad, Kh, Km)
! call solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right)
! call solve_scalar_eq (Salin_dev, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_right)
! call solve_vector_eq(U, Km, nz, dz, dt, flux_u_surf(i), flux_u_bot(i), f_u_right)
! call solve_vector_eq(V, Km, nz, dz, dt, flux_v_surf(i), flux_v_bot(i), f_v_right)
call get_mld(mld, N2, dz, nz, z)
call get_lab_mld(lab_mld, u_dynH, N2_0, time_current, z)
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)') '-------------------------------------------------'
call get_mld(mld, N2, dz, nz, z)
call get_lab_mld(lab_mld, u_dynH, N2_0, time_current, z)
call push_value_to_tseries(output_mld, mld)
call push_value_to_tseries(output_mld_ref, lab_mld)
call push_value_to_tseries(output_time, time_current / 3600.0)
enddo
output_Theta%data(:,1:output_Theta%num) = output_Theta%data(:,1:output_Theta%num) + Theta_ref
output_Salin%data(:,1:output_Salin%num) = output_Salin%data(:,1:output_Salin%num) + Salin_ref
output_TinC%data(:,1:output_TinC%num) = output_TinC%data(:,1:output_TinC%num) + Theta_ref - 273.15
if (output_mode.eq.1) then
write(*, *) ' >> writing netcdf output ...'
! Writing 2D arrays (variables with depth and time)
call write_netcdf
! Writing 1D arrays (scalar variables over time)
!call write_1d_real_nc(Tau_x_surf, 'output.nc', meta_tau_u)
!call write_1d_real_nc(Tau_y_surf, 'output.nc', meta_tau_v)
endif
if (output_mode.eq.2) then
write(*, *) ' >> writing ascii output ...'
! Writing 2D arrays (variables with depth and time) to ASCII files
call write_ascii
! Writing 1D arrays (scalar variables over time) to ASCII files
!call write_1d_real_ascii(Tau_x_surf, 'output_tau_u.txt', meta_tau_u)
!call write_1d_real_ascii(Tau_y_surf, 'output_tau_v.txt', meta_tau_v)
endif
if (output_mode.eq.3) then
write(*, *) ' >> writing tecplot output ...'
! time slice output
call write_tecplot(grid%z, output_time%data)
! Writing 1D arrays (scalar variables over time) to Tecplot files
!call write_1d_real_ascii(Tau_x_surf, 'output_tau_u.txt', meta_tau_u)
!call write_1d_real_ascii(Tau_y_surf, 'output_tau_v.txt', meta_tau_v)
endif
deallocate(dz)
deallocate(flux_Heat_surf)
deallocate(Times_flux_Heat_surf)
deallocate(flux_Heat_bot)
deallocate(Times_flux_Heat_bot)
deallocate(flux_Sal_surf)
deallocate(Times_flux_Sal_surf)
deallocate(flux_Sal_bot)
deallocate(Times_flux_Sal_bot)
deallocate(Tau_x_surf)
deallocate(Tau_y_surf)
deallocate(Tau_x_bot)
deallocate(Tau_y_bot)
deallocate(flux_u_surf)
deallocate(Times_flux_u_surf)
deallocate(flux_u_bot)
deallocate(Times_flux_u_bot)
deallocate(flux_v_surf)
deallocate(Times_flux_v_surf)
deallocate(flux_v_bot)
deallocate(Times_flux_v_bot)
deallocate(f_Heat_right)
deallocate(f_Sal_right)
deallocate(f_u_right)
deallocate(f_v_right)
deallocate(TKE)
deallocate(eps)
deallocate(P_TKE)
deallocate(B_TKE)
deallocate(Km_TKE)
!> removing time slice data
call output_cleanup
! > removing grid data
call deallocate_grid(grid)
end program