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
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) :: 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]
!< 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
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
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
! --- 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) then
time_begin = 0.0
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)
call Theta_init (Theta, Theta_dev, grid%cz, grid%dz)
call Salin_init(Salin, Salin_dev, grid%cz, grid%dz)
call get_density(Rho, Theta_dev, Salin_dev, grid%cz)
call U_init(U, grid%cz)
call V_init(V, grid%cz)
!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)
! ----------------------------------------------------------------------------
!< 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)
! ----------------------------------------------------------------------------
!open (20, file= 'output_Daria/surf_temp.txt', status ='replace')
status = 0
num = 0
do while (time_current < time_end )
! ----------------------------------------------------------------------------
!< define fluxes & dynamic scales [surface]
! ----------------------------------------------------------------------------
!< heat flux
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
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)
call get_value_tforcing(bc%salin_fluxH, time_current, salin_flux_surf)
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
!< 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(sw_flux_surf)
call deallocate_tforcing(lw_in_surf)
call deallocate_tforcing(lw_net_surf)
! > removing grid data
call deallocate_grid(grid)