#include "obl_def.fi"
program obl_main
!< @brief ocean boundary layer (OBL) main
! --------------------------------------------------------------------------------
#ifdef USE_CONFIG_PARSER
use iso_c_binding, only: C_NULL_CHAR
use config_parser
#endif
! modules used
! --------------------------------------------------------------------------------
use obl_grid
use obl_state
use obl_state_eq
use obl_common_phys
use obl_tforcing
use obl_tslice
use obl_tseries
use obl_output
use obl_init
use obl_fluxes
use obl_scm
use obl_config
use obl_diag
use obl_bc
use obl_pph, only: pphParamType, &
set_config_param_pph => set_config_param, &
define_stability_functions_pph => define_stability_functions, &
init_turbulence_closure_pph => init_turbulence_closure, &
advance_turbulence_closure_pph => advance_turbulence_closure
use obl_pph_dyn, only: pphDynParamType, &
set_config_param_pph_dyn => set_config_param, &
define_stability_functions_pph_dyn => define_stability_functions, &
init_turbulence_closure_pph_dyn => init_turbulence_closure, &
advance_turbulence_closure_pph_dyn => advance_turbulence_closure
use obl_k_epsilon, only: kepsilonParamType, &
set_config_param_k_epsilon => set_config_param, &
define_stability_functions_k_epsilon => define_stability_functions, &
init_turbulence_closure_k_epsilon => init_turbulence_closure, &
advance_turbulence_closure_k_epsilon => advance_turbulence_closure, &
TKE_init, EPS_init, TKE_bc, EPS_bc
use obl_most, only: mostParamType, &
set_config_param_most => set_config_param, &
define_stability_functions_most => define_stability_functions, &
init_turbulence_closure_most => init_turbulence_closure, &
advance_turbulence_closure_most => advance_turbulence_closure
use obl_io_plt
use io
use io_metadata
! directives list
! --------------------------------------------------------------------------------
implicit none
! model data
! --------------------------------------------------------------------------------
integer :: obl_config_id ! --- OBL builtin configs
! = obl_config_kato: Kato-Phillips
! = obl_config_papa_fluxes: Papa-station (fluxes)
! = obl_config_papa_meteo: Papa-station (meteo)
! = obl_config_cbl: convective boundary layer (Willis exp.)
! = obl_config_cyclone: cyclone setup
! = obl_config_papa_long_fluxes: Papa-station 2 years (fluxes)
! = obl_config_papa_long_meteo: Papa-station 2 years (meteo)
integer :: obl_model_id ! --- OBL model def.
! = obl_model_pph: pacanowski-philander
! = obl_model_pph_dyn: pacanowski-philander (dynamic)
! = obl_model_k_epsilon: k-epsilon
type(gridDataType) :: grid
!< turbulence closure parameters
type(pphParamType) :: param_pph
type(pphDynParamType) :: param_pph_dyn
type(kepsilonParamType) :: param_k_epsilon
type(mostParamType) :: param_most
!< boundary conditions data
type(oblBcType) :: bc
!< output
type(oblOutputStx) :: scm_output
!< time
real*8 :: time_begin_real8, time_end_real8, time_current_real8
real :: time_begin, time_end, time_current
real :: dt
integer :: time_index
!< screen output parameters
integer, parameter :: nscreen = 1000
!< file output parameters
integer :: output_mode ! --- OBL output mode
! = 1 -- netcdf
! = 2 -- ascii
! = 3 -- tecplot
integer, parameter :: noutput = 60
!< additional variables & parameters
real, parameter :: Cd0 = 0.001 ! bottom drag coefficient
real :: mld, mld_ref ! mixed layer depth (model & reference), [m]
real :: eld, eld_ref ! entrainment layer depth (model & reference), [m]
real, parameter :: N2_ref = 0.000044 ! reference N**2 (used in mld ref. estimate), [1/s**2]
! = 0.000044 (default Kato)
! = 2.0 * 1e-4 (default CBL)
real, parameter :: Bsurf_ref = 0.52 * 1e-7 ! reference buoyancy flux (used in eld ref. estimate)
! = 0.52 * 1e-7 (default CBL)
! command line arguments & configuration file variables
! --------------------------------------------------------------------------------
integer :: num_args
character(len = 128) :: arg
character(len = 128), parameter :: arg_key_help = '--help'
character(len = 128), parameter :: arg_key_config = "--config"
character(len = 128), parameter :: arg_key_model = "--model"
integer :: i, status
integer :: ierr
! --------------------------------------------------------------------------------
!< default config = obl_config_kato (Kato-Phiilps)
!< possible values:
!< = obl_config_kato
!< = obl_config_papa_fluxes
!< = obl_config_papa_meteo
!< = obl_config_cbl
!< = obl_config_cyclone
!< = obl_config_papa_long_fluxes
!< = obl_config_papa_long_meteo
obl_config_id = obl_config_kato
!< default closure = 4, k-epsilon
!< poissible values:
!< = obl_model_pph, pacanowski-philander
!< = obl_model_pph_dyn, pacanowski-philander (dynamic)
!< = obl_model_k_epsilon, k-epsilon
!< = obl_model_most, MOST scheme
obl_model_id = obl_model_k_epsilon
!< default output = 1, (netcdf)
!< possible values:
!< = 1, netcdf
!< = 2, ascii
!< = 3, tecplot
output_mode = 3
! --- 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 [key] || [filename]'
write(*, *) ' builtin setup [key] = kato (default) || papa-fluxes || papa-meteo || cbl || cyclone'
write(*, *) ' use configuration file [filename]'
write(*, *) ' --model [key]'
write(*, *) ' [key] = pph || pph-dyn || k-epsilon (default) || most'
return
end if
if (trim(arg) == trim(arg_key_config)) then
if (i == num_args) then
write(*, *) ' FAILURE! > missing config [key] argument'
ierr = 1 ! signal ERROR
return
end if
call get_command_argument(i + 1, arg)
obl_config_id = get_obl_config_id(arg)
if (obl_config_id == -1) then
#ifdef USE_CONFIG_PARSER
!< try reading configuration file
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
!< forcing configuration file setup
obl_config_id = 999
#else
write(*, *) ' FAILURE! > unknown config [key]: ', trim(arg)
ierr = 1 ! signal ERROR
return
#endif
end if
else if (trim(arg) == trim(arg_key_model)) then
if (i == num_args) then
write(*, *) ' FAILURE! > missing model [key] argument'
ierr = 1 ! signal ERROR
return
end if
call get_command_argument(i + 1, arg)
obl_model_id = get_obl_model_id(arg)
if (obl_model_id == -1) then
write(*, *) ' FAILURE! > unknown model [key]: ', trim(arg)
ierr = 1 ! signal ERROR
return
end if
endif
enddo
!< setting grid
! ----------------------------------------------------------------------------
call set_grid(grid, obl_config_id, ierr)
if (ierr /= 0) then
write(*, *) ' FAILURE! > unable to set grid '
return
endif
!< debug grid print
call print_grid(grid)
! ----------------------------------------------------------------------------
!< setting model time
! ----------------------------------------------------------------------------
call set_time(time_begin, time_end, dt, obl_config_id, ierr)
if (ierr /= 0) then
write(*, *) ' FAILURE! > unable to set time '
return
endif
! ----------------------------------------------------------------------------
!< allocating state vector
! ----------------------------------------------------------------------------
call allocate_state_vec(grid%cz)
! ----------------------------------------------------------------------------
!< initialize scm
! ----------------------------------------------------------------------------
call init_scm_vec(grid%cz)
! ----------------------------------------------------------------------------
!< setting phys
! ----------------------------------------------------------------------------
call set_phys(obl_config_id, ierr)
if (ierr /= 0) then
write(*, *) ' FAILURE! > unable to set phys parameters '
return
endif
! ----------------------------------------------------------------------------
!< setting forcing
! ----------------------------------------------------------------------------
call set_forcing(obl_config_id, ierr)
if (ierr /= 0) then
write(*, *) ' FAILURE! > unable to set forcing '
return
endif
! ----------------------------------------------------------------------------
!< setting turbulence closure parameters
! ----------------------------------------------------------------------------
call set_config_param_pph(param_pph, "turbulence_model.pph", ierr)
if (ierr /= 0) then
write(*, *) ' FAILURE! > unable to set ''pph'' parameters '
return
endif
call set_config_param_pph_dyn(param_pph_dyn, "turbulence_model.pph_dyn", ierr)
if (ierr /= 0) then
write(*, *) ' FAILURE! > unable to set ''pph-dyn'' parameters '
return
endif
call set_config_param_k_epsilon(param_k_epsilon, "turbulence_model.k_epsilon", ierr)
if (ierr /= 0) then
write(*, *) ' FAILURE! > unable to set ''k-epsilon'' parameters '
return
endif
! ----------------------------------------------------------------------------
!< initialization of main fields
! ----------------------------------------------------------------------------
call set_initial_conditions(grid, obl_config_id, ierr)
if (ierr /= 0) then
write(*, *) ' FAILURE! > unable to set initial conditions '
return
endif
Theta_dev = Theta - Theta_ref
Salin_dev = Salin - Salin_ref
call get_density(Rho, Theta_dev, Salin_dev, grid%cz)
! ----------------------------------------------------------------------------
!< initialization of turbulence closure
! ----------------------------------------------------------------------------
!< define fluxes & dynamic scales at startup
call advance_surface_fluxes(bc, time_begin, grid)
call advance_bottom_fluxes(bc, time_begin, grid)
if (obl_model_id.eq.obl_model_pph) then
call define_stability_functions_pph(param_pph, bc, grid)
call init_turbulence_closure_pph(param_pph, bc, grid)
else if (obl_model_id.eq.obl_model_pph_dyn) then
call define_stability_functions_pph_dyn(param_pph_dyn, bc, grid)
call init_turbulence_closure_pph_dyn(param_pph_dyn, bc, grid)
else if (obl_model_id.eq.obl_model_k_epsilon) then
call define_stability_functions_k_epsilon(param_k_epsilon, bc, grid)
call init_turbulence_closure_k_epsilon(param_k_epsilon, bc, grid)
else if (obl_model_id.eq.obl_model_most) then
call define_stability_functions_most(param_most, bc, grid)
call init_turbulence_closure_most(param_most, bc, grid)
endif
! ----------------------------------------------------------------------------
time_begin_real8 = time_begin
time_current_real8 = time_begin
time_end_real8 = time_end
time_current = time_begin
time_index = 1
do while (time_current < time_end )
! ----------------------------------------------------------------------------
!< define fluxes & dynamic scales [surface]
! ----------------------------------------------------------------------------
call advance_surface_fluxes(bc, time_current, grid)
! ----------------------------------------------------------------------------
!< define fluxes & dynamic scales [bottom]
! ----------------------------------------------------------------------------
call advance_bottom_fluxes(bc, time_current, grid)
! ----------------------------------------------------------------------------
!< advance turbulence closure
! ----------------------------------------------------------------------------
if (obl_model_id.eq.obl_model_pph) then
call define_stability_functions_pph(param_pph, bc, grid)
call advance_turbulence_closure_pph(param_pph, bc, grid, dt)
else if (obl_model_id.eq.obl_model_pph_dyn) then
call define_stability_functions_pph_dyn(param_pph_dyn, bc, grid)
call advance_turbulence_closure_pph_dyn(param_pph_dyn, bc, grid, dt)
else if (obl_model_id.eq.obl_model_k_epsilon) then
call define_stability_functions_k_epsilon(param_k_epsilon, bc, grid)
call advance_turbulence_closure_k_epsilon(param_k_epsilon, bc, grid, dt)
else if (obl_model_id.eq.obl_model_most) then
call define_stability_functions_most(param_most, bc, grid)
call advance_turbulence_closure_most(param_most, bc, grid, dt)
endif
! ----------------------------------------------------------------------------
!< advance single-column model eq.
! ----------------------------------------------------------------------------
call advance_scm_eq(bc, grid, dt)
! ----------------------------------------------------------------------------
!> advance time
! ----------------------------------------------------------------------------
time_index = time_index + 1
time_current_real8 = time_current_real8 + dt
time_current = time_current + dt
! ----------------------------------------------------------------------------
!> advance screen output
! ----------------------------------------------------------------------------
if (mod(time_index, nscreen) == 0) then
call get_mld(mld, N2, grid%dz, grid%cz)
call get_eld(eld, TKE_buoyancy, grid%dz, grid%cz)
! *: add finite check
write(*, '(a,g0,a,g0)') ' mld = ', mld, ' eld = ', eld
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
! ----------------------------------------------------------------------------
!> advance file output
! ----------------------------------------------------------------------------
if (mod(time_index, noutput) == 0) then
call push_state_vec(scm_output, grid%cz)
call get_mld(mld, N2, grid%dz, grid%cz)
call get_mld_ref(mld_ref, bc%U_dynH, N2_ref, time_current, grid%height)
call get_eld(eld, TKE_buoyancy, grid%dz, grid%cz)
call get_eld_ref(eld_ref, N2_ref, Bsurf_ref, time_current, grid%height)
call push_value_to_tseries(scm_output%mld, mld)
call push_value_to_tseries(scm_output%mld_ref, mld_ref)
call push_value_to_tseries(scm_output%eld, eld)
call push_value_to_tseries(scm_output%eld_ref, eld_ref)
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)
endif
! ----------------------------------------------------------------------------
enddo
!> writing file output
! ----------------------------------------------------------------------------
if (output_mode.eq.1) then
write(*, *) ' >> writing netcdf output ...'
call write_netcdf(scm_output, grid%z)
endif
if (output_mode.eq.2) then
write(*, *) ' >> writing ascii output ...'
call write_ascii(scm_output, grid%z)
endif
if (output_mode.eq.3) then
write(*, *) ' >> writing tecplot output ...'
call write_tecplot(scm_output, grid%z)
endif
! ----------------------------------------------------------------------------
!> model cleanup
! ----------------------------------------------------------------------------
!> deallocate state
call deallocate_state_vec
!> deallocate scm
call deallocate_scm_vec
!> deallocate time-dependent forcing
call deallocate_fluxes_data
!> removing time slice data
call output_cleanup(scm_output)
! > removing grid data
call deallocate_grid(grid)
! ----------------------------------------------------------------------------
end program