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
! --------------------------------------------------------------------------------
! screen output parameters
integer, parameter :: nscreen = 1000
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
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
192
193
194
! --- 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(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)
#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_Rig(Ri_grad, N2, S2, nz)
call get_hml (hml, maxcellnumber, N2, dz, nz, z)
342
343
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
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
!> advancing time
i = i + 1
time_current = time_current + dt
if (mod(i, nscreen) == 0) then
write(*, '(a,g0)') ' mld = ', hml
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)') '-------------------------------------------------'
endif
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)
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
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
! 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)
! > removing grid data
call deallocate_grid(grid)
end program