Newer
Older
program obl_main
!< @brief main program for calculations for ocean boundary layer
use obl_time_and_space
use obl_initial_conditions
use obl_forcing_and_boundary
use obl_rhs
use obl_equation
use obl_richardson
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
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 :: 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 :: 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
real, allocatable :: Rho(:) !< density, [kg / m**3]
real :: Rho_dyn_surf, Rho_dyn_bot !< density dynamic, [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]
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
real, allocatable :: Theta_write(:,:)
real, allocatable :: Salin_write(:,:)
real, allocatable :: U_write(:,:)
real, allocatable :: V_write(:,:)
real, allocatable :: Rho_write(:,:)
real, allocatable :: N2_write(:,:)
real, allocatable :: S2_write(:,:)
real, allocatable :: Ri_grad_write(:,:)
real, allocatable :: Kh_write(:,:)
real, allocatable :: Km_write(:,:)
real, allocatable :: hml_write(:)
real, allocatable :: maxcellnumber_write(:)
real, allocatable :: Theta_C_write(:,:)
real, allocatable :: lab_hml_write(:)
closure_mode = 4 !< 1 - pacanowski-philander, 2 - pacanowski-philander+, 3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
call set_dz_dt(dz, dt, nz, nt, t, z)
call set_Time_Space_Forcing(nf, df)
allocate(Theta(1:nz))
allocate(Salin(1:nz))
allocate(U(1:nz))
allocate(V(1:nz))
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(Km(1:nz))
allocate(Kh(1:nz))
allocate(Rho(1:nz))
allocate(N2(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))
allocate(V_write(nz, nt))
allocate(Rho_write(nz, nt))
allocate(N2_write(nz, nt))
allocate(S2_write(nz, nt))
allocate(Ri_grad_write(nz, nt))
allocate(Kh_write(nz, nt))
allocate(Km_write(nz, nt))
allocate(hml_write(nt))
allocate(maxcellnumber_write(nt))
allocate(Theta_C_write(nz, nt))
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 (199, file= 'output_Daria/hml.txt', status ='replace')
!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 k = 1, nz
! Km(k) = 2.0 * 0.01 !100.0
! Kh(k) = 2.0 * 0.01 !100.0
!enddo
do while (time_current < time_end )
if (closure_mode.ne.5) then
call get_value_interpolate(flux_heat_surf_res, flux_heat_surf, Times_flux_heat_surf, nf, time_current)
call get_value_interpolate(flux_heat_bot_res, flux_heat_bot, Times_flux_heat_bot, nf, time_current)
call get_value_interpolate(flux_sal_surf_res, flux_sal_surf, Times_flux_sal_surf, nf, time_current)
call get_value_interpolate(flux_sal_bot_res, flux_sal_bot, Times_flux_sal_bot, nf, time_current)
call get_value_interpolate(flux_u_surf_res, flux_u_surf, Times_flux_u_surf, nf, time_current)
call get_value_interpolate(flux_v_surf_res, flux_v_surf, Times_flux_v_surf, nf, time_current)
call get_value_interpolate(flux_u_bot_res, flux_u_bot, Times_flux_u_surf, nf, time_current)
call get_value_interpolate(flux_v_bot_res, flux_v_bot, Times_flux_v_bot, nf, time_current)
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)
call get_S2 (S2, U, V, flux_u_surf_res, flux_v_surf_res, flux_u_bot_res, flux_v_bot_res, dz, nz)
call get_Rig(Ri_grad, N2, S2, nz)
call get_hml (hml, maxcellnumber, N2, dz, nz, z)
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
call get_lab_hml(lab_hml, flux_u_surf_res, flux_v_surf_res, time_current, z)
if (closure_mode.eq.1) then
call get_Km(Km, Ri_grad, nz)
call get_Kh(Kh, Ri_grad, nz)
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, hml, 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)
!do k = 1, nz
! Theta (k) = Theta(k) + Theta_0
!end do
Theta_write(:, i) = Theta_dev
Salin_write(:, i) = Salin_dev
U_write(:, i) = U
V_write(:, i) = V
N2_write(:, i) = N2
S2_write(:, i) = S2
Ri_grad_write(:, i) = Ri_grad
Kh_write(:, i) = Kh
Km_write(:, i) = Km
hml_write(i) = hml
maxcellnumber_write(i) = maxcellnumber
!write(199,*) time_current, hml, lab_hml
!write(20,*) time_current, Theta_dev(nz) + Theta_ref
i = i+1
time_current = time_current+dt
enddo
Theta_C_write = Theta_write - 273.15
Theta_write = Theta_write + Theta_ref
Salin_write = Salin_write + Salin_ref
!Rho_write = Rho_write + Rho_ref
! Writing 2D arrays (variables with depth and time)
call write_2d_real_nc(Theta_write, 'output.nc', meta_theta)
call write_2d_real_nc(Salin_write, 'output.nc', meta_salin)
call write_2d_real_nc(U_write, 'output.nc', meta_u)
call write_2d_real_nc(V_write, 'output.nc', meta_v)
call write_2d_real_nc(Rho_write, 'output.nc', meta_rho)
call write_2d_real_nc(N2_write, 'output.nc', meta_n2)
call write_2d_real_nc(S2_write, 'output.nc', meta_s2)
call write_2d_real_nc(Ri_grad_write, 'output.nc', meta_ri_grad)
call write_2d_real_nc(Kh_write, 'output.nc', meta_kh)
call write_2d_real_nc(Km_write, 'output.nc', meta_km)
call write_2d_real_nc(Theta_C_write, 'output.nc', meta_theta_C)
! Writing 1D arrays (scalar variables over time)
call write_1d_real_nc(hml_write, 'output.nc', meta_hml)
call write_1d_real_nc(maxcellnumber_write, 'output.nc', meta_maxcellnumber)
!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)
call write_1d_real_nc(lab_hml_write, 'output.nc', meta_lab_hml)
! Writing 2D arrays (variables with depth and time) to ASCII files
call write_2d_real_ascii(Theta_write, 'output_Theta.txt', meta_theta)
call write_2d_real_ascii(Salin_write, 'output_Salin.txt', meta_salin)
call write_2d_real_ascii(U_write, 'output_U.txt', meta_u)
call write_2d_real_ascii(V_write, 'output_V.txt', meta_v)
call write_2d_real_ascii(Rho_write, 'output_Rho.txt', meta_rho)
call write_2d_real_ascii(N2_write, 'output_N2.txt', meta_n2)
call write_2d_real_ascii(S2_write, 'output_S2.txt', meta_s2)
call write_2d_real_ascii(Ri_grad_write, 'output_Ri_grad.txt', meta_ri_grad)
call write_2d_real_ascii(Kh_write, 'output_Kh.txt', meta_kh)
call write_2d_real_ascii(Km_write, 'output_Km.txt', meta_km)
! Writing 1D arrays (scalar variables over time) to ASCII files
call write_1d_real_ascii(hml_write, 'output_hml.txt', meta_hml)
call write_1d_real_ascii(lab_hml_write, 'output_lab_hml.txt', meta_lab_hml)
call write_1d_real_ascii(maxcellnumber_write, 'output_maxcellnumber.txt', meta_maxcellnumber)
!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)
!close(199)
!close(20)
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
! deallocation
deallocate(dz)
deallocate(Theta)
deallocate(Theta_dev)
deallocate(Salin)
deallocate(Salin_dev)
deallocate(U)
deallocate(V)
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(Km)
deallocate(Kh)
deallocate(Rho)
deallocate(Rho_dev)
deallocate(N2)
deallocate(S2)
deallocate(Ri_grad)
deallocate(TKE)
deallocate(eps)
deallocate(P_TKE)
deallocate(B_TKE)
deallocate(Km_TKE)
deallocate(Km_eps)
deallocate(Theta_write)
deallocate(Salin_write)
deallocate(U_write)
deallocate(V_write)
deallocate(Rho_write)
deallocate(N2_write)
deallocate(S2_write)
deallocate(Ri_grad_write)
deallocate(Kh_write)
deallocate(Km_write)
deallocate(hml_write)
deallocate(maxcellnumber_write)
deallocate(Theta_C_write)
deallocate(lab_hml_write)
end program