Skip to content
Snippets Groups Projects
Commit 3317d0d9 authored by Evgeny Mortikov's avatar Evgeny Mortikov
Browse files

adding fluxes module

parent 73f95bd8
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 1
− 0
View file @ 3317d0d9
......@@ -47,6 +47,7 @@ set(SOURCES
obl_k_epsilon.f90
obl_pph.f90
obl_pph_dyn.f90
obl_fluxes.f90
obl_scm.f90
obl_main.f90
io_metadata.f90
......
obl_common_phys.f90
+ 4
− 0
View file @ 3317d0d9
......@@ -22,6 +22,10 @@ module obl_common_phys
real, public, parameter :: R_d = 287.0 !< dry air gas constant [J / (K * kg)]
real, public, parameter :: R_v = 461.0 !< water vapor gas constant [J / (K * kg)]
real, public, parameter :: rho_air = 1.22 !< air reference density [kg / m**3]
real, public, parameter :: L_v = 2.5 * 1e6 !< latenth heat of vaporization [J/kg]
! --------------------------------------------------------------------------------
......
obl_fluxes.f90 0 → 100644
+ 278
− 0
View file @ 3317d0d9
module obl_fluxes
!< @brief obl fluxes calculation & setup
! --------------------------------------------------------------------------------
! TO DO:
! -- ***
! modules used
! --------------------------------------------------------------------------------
use obl_tforcing
#ifdef USE_SFX
use sfx_data, only: meteoDataType, sfxDataType
use sfx_most, only: get_surface_fluxes_most => get_surface_fluxes, &
numericsType_most => numericsType
#endif
! directives list
! --------------------------------------------------------------------------------
implicit none
private
! public interface
! --------------------------------------------------------------------------------
public :: advance_surface_fluxes
public :: advance_bottom_fluxes
public :: deallocate_fluxes_data
! --------------------------------------------------------------------------------
!> @brief setup mode [bool]
integer, public :: is_meteo_setup
!> @brief surface forcing
type(timeForcingDataType), public :: sensible_hflux_surf, latent_hflux_surf !< heat fluxes, [W/m^2]
type(timeForcingDataType), public :: salin_flux_surf !< salinity flux, [PSU*m/s]
type(timeForcingDataType), public :: tau_x_surf, tau_y_surf !< momentum flux, [N/m**2]
type(timeForcingDataType), public :: Ua, Va !< air wind, [m/s]
type(timeForcingDataType), public :: Ta !< air temperature, [K]
type(timeForcingDataType), public :: Pa !< air pressure, [Pa]
type(timeForcingDataType), public :: Qa !< air humidity, [kg/kg]
type(timeForcingDataType), public :: RHa !< air relative humidity, [%]
type(timeForcingDataType), public :: sw_flux_surf !< shortwave radiation flux IN, [W/m^2]
type(timeForcingDataType), public :: lw_in_surf, lw_net_surf !< longwave radiation flux IN/NET, [W/m^2]
!> @brief bottom forcing
type(timeForcingDataType), public :: hflux_bot !< heat flux, [W/m^2]
type(timeForcingDataType), public :: salin_flux_bot !< salinity flux, [PSU*m/s]
type(timeForcingDataType), public :: tau_x_bot, tau_y_bot !< momentum flux, [N/m**2]
real, parameter :: lw_albedo = 0.03
real, parameter :: surface_emissivity = 0.97
#ifdef USE_SFX
type(meteoDataType) :: meteo
type(sfxDataType) :: sfx
type(numericsType_most) :: sfx_numerics
#endif
contains
! --------------------------------------------------------------------------------
subroutine advance_surface_fluxes(bc, time_current, grid)
!< @brief set fluxes & dynamic scales b.c. (surface)
! ----------------------------------------------------------------------------
use obl_grid
use obl_state
use obl_bc
use obl_state_eq
use obl_common_phys
type(oblBcType), intent(out) :: bc
real, intent(in) :: time_current
type (gridDataType), intent(in) :: grid
real :: fvalue
real :: lw_net_t
! --------------------------------------------------------------------------------
!< define LW net radiation
lw_net_t = 0.0
if (lw_net_surf%num > 0) then
call get_value_tforcing(lw_net_t, time_current, lw_net_surf)
else if (lw_in_surf%num > 0) then
call get_value_tforcing(fvalue, time_current, lw_in_surf)
lw_net_t = (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
lw_net_t = lw_net_t - lw_out_surf
end block
endif
if (is_meteo_setup == 0) then
!< using 'flux' type forcing
! --------------------------------------------------------------------------------
!< 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
bc%heat_fluxH = bc%heat_fluxH + lw_net_t
!< 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)
else if (is_meteo_setup == 1) then
#ifdef USE_SFX
!< using 'meteo' type forcing
! --------------------------------------------------------------------------------
block
real :: Ua_t, Va_t, Ta_t, Pa_t, Qa_t, RHa_t, e_sat_a
real :: Theta_surf, Q_surf, e_sat
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)
!< define air specific humidity
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
!< define saturation specific humidity
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_air * cp_air * sfx%Cm * sfx%Ct * meteo%U * meteo%dT + &
rho_air * L_v * sfx%Cm * sfx%Ct * meteo%U * meteo%dQ
bc%heat_fluxH = bc%heat_fluxH + lw_net_t
!< momentum flux
bc%u_momentum_fluxH = rho_air * sfx%Cm * sfx%Cm * meteo%U * Ua_t
bc%v_momentum_fluxH = rho_air * sfx%Cm * sfx%Cm * meteo%U * Va_t
end block
#else
write(*, *) ' >> FAILURE! Meteo setup requires SFX'
stop
#endif
endif
!< kinematic heat flux = F / (rho_ref * cp)
bc%heat_fluxH = bc%heat_fluxH / (Rho_ref * cp_water)
!< 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
!< salinity flux
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)
!< U* def.:
call get_u_dynamic(bc%U_dynH, bc%u_momentum_fluxH, bc%v_momentum_fluxH)
!< rho* def.:
call get_rho_dynamic(bc%rho_dynH, &
bc%U_dynH, bc%heat_fluxH, bc%salin_fluxH)
end subroutine
! --------------------------------------------------------------------------------
subroutine advance_bottom_fluxes(bc, time_current, grid)
!< @brief set fluxes & dynamic scales b.c. (bottom)
! ----------------------------------------------------------------------------
use obl_grid
use obl_state
use obl_bc
use obl_state_eq
use obl_common_phys
type(oblBcType), intent(out) :: bc
real, intent(in) :: time_current
type (gridDataType), intent(in) :: grid
! --------------------------------------------------------------------------------
!< 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)
!< salinity flux
call get_value_tforcing(bc%salin_flux0, time_current, salin_flux_bot)
!< momentum flux
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)
!< rho* def.:
call get_rho_dynamic(bc%rho_dyn0, &
bc%U_dyn0, bc%heat_flux0, bc%salin_flux0)
end subroutine
! --------------------------------------------------------------------------------
subroutine deallocate_fluxes_data
!> @brief deallocate fluxes 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)
end subroutine deallocate_fluxes_data
end module
obl_main.f90
+ 5
− 220
View file @ 3317d0d9
......@@ -18,6 +18,7 @@ program obl_main
use obl_tseries
use obl_output
use obl_init
use obl_fluxes
use obl_scm
use obl_diag
use obl_bc
......@@ -34,11 +35,6 @@ program obl_main
use obl_io_plt
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, &
numericsType_most => numericsType
#endif
#ifdef USE_CONFIG_PARSER
use config_parser
#endif
......@@ -48,12 +44,6 @@ program obl_main
! directives list
implicit none
#ifdef USE_SFX
type(meteoDataType) :: meteo
type(sfxDataType) :: sfx
type(numericsType_most) :: sfx_numerics
#endif
type(gridDataType) :: grid
! turbulence closure parameters
......@@ -61,27 +51,6 @@ program obl_main
type(pphDynParamType) :: param_pph_dyn
type(kepsilonParamType) :: param_k_epsilon
!< 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
......@@ -108,16 +77,11 @@ program obl_main
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
real :: N2_0
! command line arguments
......@@ -392,165 +356,12 @@ program obl_main
!< 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
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
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
!< salinity flux
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)
!< U* def.:
call get_u_dynamic(bc%U_dynH, bc%u_momentum_fluxH, bc%v_momentum_fluxH)
!< rho* def.:
call get_rho_dynamic(bc%rho_dynH, &
bc%U_dynH, bc%heat_fluxH, bc%salin_fluxH)
call advance_surface_fluxes(bc, time_current, grid)
! ----------------------------------------------------------------------------
!< 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)
!< salinity flux
call get_value_tforcing(bc%salin_flux0, time_current, salin_flux_bot)
!< momentum flux
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)
!< rho* def.:
call get_rho_dynamic(bc%rho_dyn0, &
bc%U_dyn0, bc%heat_flux0, bc%salin_flux0)
call advance_bottom_fluxes(bc, time_current, grid)
! ----------------------------------------------------------------------------
!< advance turbulence closure
......@@ -634,34 +445,8 @@ program obl_main
!> 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)
!> deallocate time-dependent forcing
call deallocate_fluxes_data
!> removing time slice data
call output_cleanup(scm_output)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment