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_pph, only: pphParamType, &
define_stability_functions_pph => define_stability_functions, &
advance_turbulence_closure_pph => advance_turbulence_closure
use obl_pph_dyn, only: pphDynParamType, &
define_stability_functions_pph_dyn => define_stability_functions, &
advance_turbulence_closure_pph_dyn => advance_turbulence_closure
use obl_k_epsilon, only: kepsilonParamType, &
define_stability_functions_k_epsilon => define_stability_functions, &
advance_turbulence_closure_k_epsilon => advance_turbulence_closure, &
TKE_init, EPS_init, TKE_bc, EPS_bc
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
#ifdef USE_SFX
type(meteoDataType) :: meteo
type(sfxDataType) :: sfx
type(numericsType_most) :: sfx_numerics
#endif
type(pphParamType) :: param_pph
type(pphDynParamType) :: param_pph_dyn
!< surface forcing
type(timeForcingDataType) :: sensible_hflux_surf, latent_hflux_surf !< heat fluxes, [W/m^2]
type(timeForcingDataType) :: salin_flux_surf !< salinity flux, [PSU*m/s]
type(timeForcingDataType) :: tau_x_surf, tau_y_surf !< momentum flux, [N/m**2]
type(timeForcingDataType) :: Ua, Va !< air wind, [m/s]
type(timeForcingDataType) :: Ta !< air temperature, [K]
type(timeForcingDataType) :: Pa !< air pressure, [Pa]
type(timeForcingDataType) :: Qa !< air humidity, [kg/kg]
type(timeForcingDataType) :: RHa !< air relative humidity, [%]
type(timeForcingDataType) :: sw_flux_surf !< shortwave radiation flux IN, [W/m^2]
type(timeForcingDataType) :: lw_in_surf, lw_net_surf !< longwave radiation flux IN/NET, [W/m^2]
integer :: is_meteo_setup
!< bottom forcing
type(timeForcingDataType) :: hflux_bot !< heat flux, [W/m^2]
type(timeForcingDataType) :: salin_flux_bot !< salinity flux, [PSU*m/s]
type(timeForcingDataType) :: tau_x_bot, tau_y_bot !< momentum flux, [N/m**2]
!< boundary conditions data
type(oblBcType) :: bc
real :: domain_height
integer :: grid_cz
!< output
type(oblOutputStx) :: scm_output
real :: dt !< time step [s]
integer :: i, k !< counters
integer :: status, num !< for file input/output
real :: time_begin, time_end, time_current !< time for output
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 = 3 ! 1 -- netcdf, 2 -- ascii, 3 -- tecplot
integer, parameter :: obl_setup = 1 ! 1 - Kato-Phillips, 2 - Papa station (fluxes), 3 - Papa station (meteo)
real, parameter :: Cd0 = 0.001 ! bottom drag coefficient
real, parameter :: lw_albedo = 0.03
real, parameter :: surface_emissivity = 0.97
real :: mld !< mixed layer depth, [m]
real :: lab_mld !< theoretical mixed layer depth, [m]
!< just a temporary to hold value
real :: fvalue
! 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
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
! --- 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, domain_height, 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, grid_cz, 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, time_end, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
end if
call c_config_is_varname("time.begin"//C_NULL_CHAR, status)
call c_config_get_float("time.begin"//C_NULL_CHAR, time_begin, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
end if
call c_config_is_varname("time.dt"//C_NULL_CHAR, status)
if (status /= 0) then
call c_config_get_float("time.dt"//C_NULL_CHAR, dt, status)
if (status == 0) then
ierr = 1 ! signal ERROR
return
end if
end if
#endif
endif
enddo
!< setting model time
! ----------------------------------------------------------------------------
if (obl_setup == 1) then
time_begin = 0.0
time_end = 300.0 * 3600.0
dt = 1.0
domain_height = 100.0
grid_cz = 32
endif
if ((obl_setup == 2).or.(obl_setup == 3)) then
domain_height = 128.0
grid_cz = 32
#endif
time_current = time_begin
! ----------------------------------------------------------------------------
! setting grid -- just an example
! > [zpos, height, cz, gcz]
call set_uniform_grid(grid, 0.0, domain_height, grid_cz)
! debug grid print
call print_grid(grid)
! initialize scm
call init_scm_vec(grid%cz)
!< setting phys
! ----------------------------------------------------------------------------
if (obl_setup == 1) then
f = 0.0
endif
if ((obl_setup == 2).or.(obl_setup == 3)) then
f = 1.116 * 1e-4
a_band_ratio = 0.67
a_extinction_coeff = 1.0
b_extinction_coeff = 1.0 / 17.0
endif
! ----------------------------------------------------------------------------
!< initialization of main fields
! ----------------------------------------------------------------------------
if (obl_setup == 1) then
call set_linear_profile(Theta, 330.0, 0.3, grid)
call set_const_profile(Salin, 35.0, grid)
endif
if ((obl_setup == 2).or.(obl_setup == 3)) then
call set_external_profile(Theta, 'PAPA_06_2017/Theta.dat', grid)
call set_external_profile(Salin, 'PAPA_06_2017/Salin.dat', grid)
endif
Theta_dev = Theta - Theta_ref
Salin_dev = Salin - Salin_ref
call get_density(Rho, Theta_dev, Salin_dev, grid%cz)
call set_const_profile(U, 0.0, grid)
call set_const_profile(V, 0.0, grid)
!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, param_k_epsilon, grid%cz)
call eps_init(EPS, param_k_epsilon, grid%cz, grid%height)
! ----------------------------------------------------------------------------
!< setting atmospheric forcing
! ----------------------------------------------------------------------------
if (obl_setup == 1) then
call set_const_tforcing(sensible_hflux_surf, 0.0)
call set_const_tforcing(latent_hflux_surf, 0.0)
call set_const_tforcing(salin_flux_surf, 0.0)
call set_const_tforcing(tau_x_surf, 0.1)
call set_const_tforcing(tau_y_surf, 0.0)
call set_const_tforcing(sw_flux_surf, 0.0)
call set_const_tforcing(lw_net_surf, 0.0)
endif
if (obl_setup == 2) then
call set_external_tforcing(sensible_hflux_surf, 'PAPA_06_2017/sensible_hflux.dat')
call set_external_tforcing(latent_hflux_surf, 'PAPA_06_2017/latent_hflux.dat')
call set_const_tforcing(salin_flux_surf, 0.0)
call set_external_tforcing(tau_x_surf, 'PAPA_06_2017/tau-x.dat')
call set_external_tforcing(tau_y_surf, 'PAPA_06_2017/tau-y.dat')
call set_external_tforcing(sw_flux_surf, 'PAPA_06_2017/SW-down.dat')
call set_external_tforcing(lw_in_surf, 'PAPA_06_2017/LW-down.dat')
!< normalize time in external forcing: hrs -> sec
call normalize_time_tforcing(sensible_hflux_surf, 3600.0)
call normalize_time_tforcing(latent_hflux_surf, 3600.0)
call normalize_time_tforcing(tau_x_surf, 3600.0)
call normalize_time_tforcing(tau_y_surf, 3600.0)
call normalize_time_tforcing(sw_flux_surf, 3600.0)
call normalize_time_tforcing(lw_in_surf, 3600.0)
341
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
if (obl_setup == 3) then
is_meteo_setup = 1
call set_external_tforcing(Ua, 'PAPA_06_2017/u-wind.dat')
call set_external_tforcing(Va, 'PAPA_06_2017/v-wind.dat')
call set_const_tforcing(salin_flux_surf, 0.0)
call set_external_tforcing(Ta, 'PAPA_06_2017/Tair.dat')
call set_external_tforcing(Pa, 'PAPA_06_2017/Pair.dat')
call set_external_tforcing(RHa, 'PAPA_06_2017/RHair.dat')
call set_external_tforcing(sw_flux_surf, 'PAPA_06_2017/SW-down.dat')
call set_external_tforcing(lw_in_surf, 'PAPA_06_2017/LW-down.dat')
!< normalize time in external forcing: hrs -> sec
call normalize_time_tforcing(Ua, 3600.0)
call normalize_time_tforcing(Va, 3600.0)
call normalize_time_tforcing(Ta, 3600.0)
call normalize_time_tforcing(Pa, 3600.0)
call normalize_time_tforcing(RHa, 3600.0)
call normalize_time_tforcing(sw_flux_surf, 3600.0)
call normalize_time_tforcing(lw_in_surf, 3600.0)
endif
! ----------------------------------------------------------------------------
!< setting bottom forcing
! ----------------------------------------------------------------------------
call set_const_tforcing(hflux_bot, 0.0)
call set_const_tforcing(salin_flux_bot, 0.0)
call set_const_tforcing(tau_x_bot, 0.0)
call set_const_tforcing(tau_y_bot, 0.0)
! ----------------------------------------------------------------------------
!< setting reference data
! ----------------------------------------------------------------------------
N2_0 = 0.00044
! ----------------------------------------------------------------------------
status = 0
num = 0
i=1
do while (time_current < time_end )
! ----------------------------------------------------------------------------
!< define fluxes & dynamic scales [surface]
! ----------------------------------------------------------------------------
if (is_meteo_setup == 0) then
!< heat flux
bc%heat_fluxH = 0.0
call get_value_tforcing(fvalue, time_current, sensible_hflux_surf)
bc%heat_fluxH = bc%heat_fluxH + fvalue
call get_value_tforcing(fvalue, time_current, latent_hflux_surf)
bc%heat_fluxH = bc%heat_fluxH + fvalue
if (lw_net_surf%num > 0) then
call get_value_tforcing(fvalue, time_current, lw_net_surf)
bc%heat_fluxH = bc%heat_fluxH + fvalue
else if (lw_in_surf%num > 0) then
call get_value_tforcing(fvalue, time_current, lw_in_surf)
bc%heat_fluxH = bc%heat_fluxH + (1.0 - lw_albedo) * fvalue
!< adding LW outgoing flux to balance
block
real :: Theta_surf
real :: lw_out_surf
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
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
Theta_surf = Theta_dev(grid%cz) + Theta_ref
lw_out_surf = surface_emissivity * stefan_boltzmann_const * &
Theta_surf * Theta_surf * Theta_surf * Theta_surf
bc%heat_fluxH = bc%heat_fluxH - lw_out_surf
end block
endif
!< kinematic heat flux = F / (rho_ref * cp)
bc%heat_fluxH = bc%heat_fluxH / (Rho_ref * cp_water)
!< momentum flux
call get_value_tforcing(bc%u_momentum_fluxH, time_current, tau_x_surf)
call get_value_tforcing(bc%v_momentum_fluxH, time_current, tau_y_surf)
!< kinematic momentum flux = tau / rho_ref
bc%u_momentum_fluxH = bc%u_momentum_fluxH / Rho_ref
bc%v_momentum_fluxH = bc%v_momentum_fluxH / Rho_ref
else if (is_meteo_setup == 1) then
#ifdef USE_SFX
block
real :: Ua_t, Va_t, Ta_t, Pa_t, Qa_t, RHa_t, e_sat_a
real :: Theta_surf, Q_surf, e_sat
real :: rho_a
real :: Lv
rho_a = 1.22
Lv = 2.5 * 1e6
call get_value_tforcing(Ua_t, time_current, Ua)
call get_value_tforcing(Va_t, time_current, Va)
call get_value_tforcing(Ta_t, time_current, Ta)
call get_value_tforcing(Pa_t, time_current, Pa)
if (Qa%num > 0) then
call get_value_tforcing(Qa_t, time_current, Qa)
else if (RHa%num > 0) then
call get_value_tforcing(RHa_t, time_current, RHa)
call get_water_saturation_vapor_pressure_in_kPa(e_sat_a, Ta_t)
call get_specific_humidity(Qa_t, e_sat_a, Pa_t / 1000.0, R_d, R_v)
Qa_t = Qa_t * (RHa_t / 100.0)
endif
Theta_surf = Theta_dev(grid%cz) + Theta_ref
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
496
497
498
499
500
501
call get_water_saturation_vapor_pressure_in_kPa(e_sat, Theta_surf)
call get_specific_humidity(Q_surf, e_sat, Pa_t / 1000.0, R_d, R_v)
meteo%U = sqrt(Ua_t * Ua_t + Va_t * Va_t)
meteo%dT = Ta_t - Theta_surf
meteo%Tsemi = 0.5 * (Ta_t + Theta_surf)
meteo%dQ = Qa_t - Q_surf
meteo%h = 2.0
meteo%z0_m = -1.0
call get_surface_fluxes_most(sfx, meteo, sfx_numerics)
!< heat flux
bc%heat_fluxH = &
rho_a * cp_air * sfx%Cm * sfx%Ct * meteo%U * meteo%dT + &
rho_a * Lv * sfx%Cm * sfx%Ct * meteo%U * meteo%dQ
if (lw_net_surf%num > 0) then
call get_value_tforcing(fvalue, time_current, lw_net_surf)
bc%heat_fluxH = bc%heat_fluxH + fvalue
else if (lw_in_surf%num > 0) then
call get_value_tforcing(fvalue, time_current, lw_in_surf)
bc%heat_fluxH = bc%heat_fluxH + (1.0 - lw_albedo) * fvalue
!< adding LW outgoing flux to balance
block
real :: Theta_surf
real :: lw_out_surf
Theta_surf = Theta_dev(grid%cz) + Theta_ref
lw_out_surf = surface_emissivity * stefan_boltzmann_const * &
Theta_surf * Theta_surf * Theta_surf * Theta_surf
bc%heat_fluxH = bc%heat_fluxH - lw_out_surf
end block
endif
!< kinematic heat flux = F / (rho_ref * cp)
bc%heat_fluxH = bc%heat_fluxH / (Rho_ref * cp_water)
!< momentum flux
bc%u_momentum_fluxH = rho_a * sfx%Cm * sfx%Cm * meteo%U * Ua_t
bc%v_momentum_fluxH = rho_a * sfx%Cm * sfx%Cm * meteo%U * Va_t
!< kinematic momentum flux = tau / rho_ref
bc%u_momentum_fluxH = bc%u_momentum_fluxH / Rho_ref
bc%v_momentum_fluxH = bc%v_momentum_fluxH / Rho_ref
end block
#endif
endif
call get_value_tforcing(bc%salin_fluxH, time_current, salin_flux_surf)
!< shortwave radiation flux [W/m^2]
call get_value_tforcing(bc%sw_fluxH, time_current, sw_flux_surf)
call get_u_dynamic(bc%U_dynH, bc%u_momentum_fluxH, bc%v_momentum_fluxH)
call get_rho_dynamic(bc%rho_dynH, &
bc%U_dynH, bc%heat_fluxH, bc%salin_fluxH)
! ----------------------------------------------------------------------------
!< define fluxes & dynamic scales [bottom]
! ----------------------------------------------------------------------------
!< heat flux
call get_value_tforcing(bc%heat_flux0, time_current, hflux_bot)
!< kinematic heat flux = F / (rho_ref * cp)
bc%heat_flux0 = bc%heat_flux0 / (Rho_ref * cp_water)
call get_value_tforcing(bc%salin_flux0, time_current, salin_flux_bot)
call get_value_tforcing(bc%u_momentum_flux0, time_current, tau_x_bot)
call get_value_tforcing(bc%v_momentum_flux0, time_current, tau_y_bot)
!< kinematic momentum flux = tau / rho_ref
bc%u_momentum_flux0 = bc%u_momentum_flux0 / Rho_ref
bc%v_momentum_flux0 = bc%v_momentum_flux0 / Rho_ref
!< U* def.:
call get_u_dynamic(bc%U_dyn0, bc%u_momentum_flux0, bc%v_momentum_flux0)
call get_rho_dynamic(bc%rho_dyn0, &
bc%U_dyn0, bc%heat_flux0, bc%salin_flux0)
! ----------------------------------------------------------------------------
!< advance turbulence closure
! ----------------------------------------------------------------------------
call define_stability_functions_pph(param_pph, bc, grid)
call advance_turbulence_closure_pph(param_pph, bc, grid, dt)
call define_stability_functions_pph_dyn(param_pph_dyn, bc, grid)
call advance_turbulence_closure_pph_dyn(param_pph_dyn, bc, grid, dt)
call define_stability_functions_k_epsilon(param_k_epsilon, bc, grid)
call advance_turbulence_closure_k_epsilon(param_k_epsilon, bc, grid, dt)
! ----------------------------------------------------------------------------
! ----------------------------------------------------------------------------
! ----------------------------------------------------------------------------
call get_mld(mld, N2, grid%dz, grid%cz)
call get_mld_ref(lab_mld, bc%U_dynH, N2_0, time_current, grid%height)
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 push_state_vec(scm_output, grid%cz)
call get_mld(mld, N2, grid%dz, grid%cz)
call get_mld_ref(lab_mld, bc%U_dynH, N2_0, time_current, grid%height)
call push_value_to_tseries(scm_output%mld, mld)
call push_value_to_tseries(scm_output%mld_ref, lab_mld)
call push_value_to_tseries(scm_output%tau_x, bc%u_momentum_fluxH * Rho_ref)
call push_value_to_tseries(scm_output%tau_y, bc%v_momentum_fluxH * Rho_ref)
call push_value_to_tseries(scm_output%Theta_surf, Theta_dev(grid%cz) + Theta_ref)
call push_value_to_tseries(scm_output%time, time_current / 3600.0)
enddo
if (output_mode.eq.1) then
write(*, *) ' >> writing netcdf output ...'
call write_netcdf(scm_output, grid%z)
if (output_mode.eq.2) then
write(*, *) ' >> writing ascii output ...'
endif
if (output_mode.eq.3) then
write(*, *) ' >> writing tecplot output ...'
call write_tecplot(scm_output, grid%z)
!> deallocate scm
call deallocate_scm_vec
!> removing time-dependent forcing data
call deallocate_tforcing(sensible_hflux_surf)
call deallocate_tforcing(latent_hflux_surf)
call deallocate_tforcing(salin_flux_surf)
call deallocate_tforcing(tau_x_surf)
call deallocate_tforcing(tau_y_surf)
call deallocate_tforcing(hflux_bot)
call deallocate_tforcing(salin_flux_bot)
call deallocate_tforcing(tau_x_bot)
call deallocate_tforcing(tau_y_bot)
call deallocate_tforcing(Ua)
call deallocate_tforcing(Va)
call deallocate_tforcing(Ta)
call deallocate_tforcing(Pa)
call deallocate_tforcing(Qa)
call deallocate_tforcing(RHa)
call deallocate_tforcing(sw_flux_surf)
call deallocate_tforcing(lw_in_surf)
call deallocate_tforcing(lw_net_surf)
! > removing grid data
call deallocate_grid(grid)