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_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
! 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
real :: value
! --------------------------------------------------------------------------------
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
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
187
188
189
190
191
! --- 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)
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)
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
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)
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
! 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