#include "../includeF/sfx_def.fi"
module sfx_sheba
!< @brief SHEBA surface flux module
! modules used
! --------------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
use sfx_common
#endif
use sfx_data
use sfx_surface
use sfx_sheba_param
#if defined(INCLUDE_CUDA) || defined(INCLUDE_CXX)
use iso_c_binding, only: c_loc
use C_FUNC
#endif
! --------------------------------------------------------------------------------
! directives list
! --------------------------------------------------------------------------------
implicit none
private
! --------------------------------------------------------------------------------
! public interface
! --------------------------------------------------------------------------------
public :: get_surface_fluxes
public :: get_surface_fluxes_vec
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
type, public :: numericsType
integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness
end type
! --------------------------------------------------------------------------------
contains
#if defined(INCLUDE_CXX)
subroutine set_c_struct_sfx_sheba_param_values(sfx_model_param)
type (sfx_sheba_param), intent(inout) :: sfx_model_param
sfx_model_param%kappa = kappa
sfx_model_param%Pr_t_0_inv = Pr_t_0_inv
sfx_model_param%alpha_m = alpha_m
sfx_model_param%alpha_h = alpha_h
sfx_model_param%a_m = a_m
sfx_model_param%b_m = b_m
sfx_model_param%a_h = a_h
sfx_model_param%b_h = b_h
sfx_model_param%c_h = c_h
end subroutine set_c_struct_sfx_sheba_param_values
#endif
! --------------------------------------------------------------------------------
subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n)
!< @brief surface flux calculation for array data
!< @details contains C/C++ & CUDA interface
! ----------------------------------------------------------------------------
type (sfxDataVecType), intent(inout) :: sfx
type (meteoDataVecType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
integer, intent(in) :: n
! ----------------------------------------------------------------------------
! --- local variables
type (meteoDataType) meteo_cell
type (sfxDataType) sfx_cell
integer i
! ----------------------------------------------------------------------------
#if defined(INCLUDE_CXX)
type (meteoDataVecTypeC), pointer :: meteo_c !< meteorological data (input)
type (sfxDataVecTypeC), pointer :: sfx_c !< surface flux data (output)
type (sfx_sheba_param) :: model_param
type (sfx_surface_param) :: surface_param
type (sfx_sheba_numericsTypeC) :: numerics_c
type (sfx_phys_constants) :: phys_constants
numerics_c%maxiters_charnock = numerics%maxiters_charnock
phys_constants%Pr_m = Pr_m;
phys_constants%nu_air = nu_air;
phys_constants%g = g;
call set_c_struct_sfx_sheba_param_values(model_param)
call set_c_struct_sfx_surface_param_values(surface_param)
call set_meteo_vec_c(meteo, meteo_c)
call set_sfx_vec_c(sfx, sfx_c)
! call get_surface_fluxes_sheba(c_loc(sfx_c), c_loc(meteo_c), model_param, surface_param, numerics_c, phys_constants, n)
deallocate(meteo_c)
deallocate(sfx_c)
#else
do i = 1, n
meteo_cell = meteoDataType(&
h = meteo%h(i), &
U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
z0_m = meteo%z0_m(i))
call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
call push_sfx_data(sfx, sfx_cell, i)
end do
#endif
end subroutine get_surface_fluxes_vec
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
subroutine get_surface_fluxes(sfx, meteo, numerics)
!< @brief surface flux calculation for single cell
!< @details contains C/C++ interface
! ----------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
use ieee_arithmetic
#endif
type (sfxDataType), intent(out) :: sfx
type (meteoDataType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
! ----------------------------------------------------------------------------
! --- meteo derived datatype name shadowing
! ----------------------------------------------------------------------------
real :: h !< constant flux layer height [m]
real :: U !< abs(wind speed) at 'h' [m/s]
real :: dT !< difference between potential temperature at 'h' and at surface [K]
real :: Tsemi !< semi-sum of potential temperature at 'h' and at surface [K]
real :: dQ !< difference between humidity at 'h' and at surface [g/g]
real :: z0_m !< surface aerodynamic roughness (should be < 0 for water bodies surface)
! ----------------------------------------------------------------------------
! --- local variables
! ----------------------------------------------------------------------------
real z0_t !< thermal roughness [m]
real B !< = ln(z0_m / z0_t) [n/d]
real h0_m, h0_t !< = h / z0_m, h / z0_h [n/d]
real u_dyn0 !< dynamic velocity in neutral conditions [m/s]
real Re !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]
real zeta !< = z/L [n/d]
real Rib !< bulk Richardson number
real Udyn, Tdyn, Qdyn !< dynamic scales
real phi_m, phi_h !< stability functions (momentum) & (heat) [n/d]
real Km !< eddy viscosity coeff. at h [m^2/s]
real Pr_t_inv !< invese Prandt number [n/d]
real Cm, Ct !< transfer coeff. for (momentum) & (heat) [n/d]
integer surface_type !< surface type = (ocean || land)
#ifdef SFX_CHECK_NAN
real NaN
#endif
! ----------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
! --- checking if arguments are finite
if (.not.(is_finite(meteo%U).and.is_finite(meteo%Tsemi).and.is_finite(meteo%dT).and.is_finite(meteo%dQ) &
.and.is_finite(meteo%z0_m).and.is_finite(meteo%h))) then
NaN = ieee_value(0.0, ieee_quiet_nan) ! setting NaN
sfx = sfxDataType(zeta = NaN, Rib = NaN, &
Re = NaN, B = NaN, z0_m = NaN, z0_t = NaN, &
Rib_conv_lim = NaN, &
Cm = NaN, Ct = NaN, Km = NaN, Pr_t_inv = NaN)
return
end if
#endif
! --- shadowing names for clarity
U = meteo%U
Tsemi = meteo%Tsemi
dT = meteo%dT
dQ = meteo%dQ
h = meteo%h
z0_m = meteo%z0_m
! --- define surface type
if (z0_m < 0.0) then
surface_type = surface_ocean
else
surface_type = surface_land
end if
if (surface_type == surface_ocean) then
! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
! --- define relative height
h0_m = h / z0_m
endif
if (surface_type == surface_land) then
! --- define relative height
h0_m = h / z0_m
! --- define dynamic velocity in neutral conditions
u_dyn0 = U * kappa / log(h0_m)
end if
! --- define thermal roughness & B = log(z0_m / z0_h)
Re = u_dyn0 * z0_m / nu_air
call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
! --- define relative height [thermal]
h0_t = h / z0_t
! --- define Ri-bulk
Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2
! --- get the fluxes
! ----------------------------------------------------------------------------
call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10)
! ----------------------------------------------------------------------------
call get_phi(phi_m, phi_h, zeta)
! ----------------------------------------------------------------------------
! --- define transfer coeff. (momentum) & (heat)
Cm = 0.0
if (U > 0.0) then
Cm = Udyn / U
end if
Ct = 0.0
if (abs(dT) > 0.0) then
Ct = Tdyn / dT
end if
! --- define eddy viscosity & inverse Prandtl number
Km = kappa * Cm * U * h / phi_m
Pr_t_inv = phi_m / phi_h
! --- setting output
sfx = sfxDataType(zeta = zeta, Rib = Rib, &
Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
Rib_conv_lim = 0.0, &
Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
!< @brief get dynamic scales
! --------------------------------------------------------------------------------
subroutine get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
U, Tsemi, dT, dQ, z, z0_m, z0_t, beta, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: Udyn, Tdyn, Qdyn !< dynamic scales
real, intent(out) :: zeta !< = z/L
real, intent(in) :: U !< abs(wind speed) at z
real, intent(in) :: Tsemi !< semi-sum of temperature at z and at surface
real, intent(in) :: dT, dQ !< temperature & humidity difference between z and at surface
real, intent(in) :: z !< constant flux layer height
real, intent(in) :: z0_m, z0_t !< roughness parameters
real, intent(in) :: beta !< buoyancy parameter
integer, intent(in) :: maxiters !< maximum number of iterations
! ----------------------------------------------------------------------------
! --- local variables
real, parameter :: gamma = 0.61
real :: psi_m, psi_h
real :: psi0_m, psi0_h
real :: Linv
integer :: i
! ----------------------------------------------------------------------------
Udyn = kappa * U / log(z / z0_m)
Tdyn = kappa * dT * Pr_t_0_inv / log(z / z0_t)
Qdyn = kappa * dQ * Pr_t_0_inv / log(z / z0_t)
zeta = 0.0
! --- no wind
if (Udyn < 1e-5) return
Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn)
zeta = z * Linv
! --- near neutral case
if (Linv < 1e-5) return
do i = 1, maxiters
call get_psi(psi_m, psi_h, zeta)
call get_psi_mh(psi0_m, psi0_h, z0_m * Linv, z0_t * Linv)
Udyn = kappa * U / (log(z / z0_m) - (psi_m - psi0_m))
Tdyn = kappa * dT * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
Qdyn = kappa * dQ * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
if (Udyn < 1e-5) exit
Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn)
zeta = z * Linv
end do
end subroutine get_dynamic_scales
! --------------------------------------------------------------------------------
! stability functions
! --------------------------------------------------------------------------------
subroutine get_phi(phi_m, phi_h, zeta)
!< @brief stability functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: phi_m, phi_h !< stability functions
real, intent(in) :: zeta !< = z/L
! ----------------------------------------------------------------------------
if (zeta >= 0.0) then
phi_m = 1.0 + (a_m * zeta * (1.0 + zeta)**(1.0 / 3.0)) / (1.0 + b_m * zeta)
phi_h = 1.0 + (a_h * zeta + b_h * zeta * zeta) / (1.0 + c_h * zeta + zeta * zeta)
else
phi_m = (1.0 - alpha_m * zeta)**(-0.25)
phi_h = (1.0 - alpha_h * zeta)**(-0.5)
end if
end subroutine
! --------------------------------------------------------------------------------
! universal functions
! --------------------------------------------------------------------------------
subroutine get_psi(psi_m, psi_h, zeta)
!< @brief universal functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions
real, intent(in) :: zeta !< = z/L
! ----------------------------------------------------------------------------
! --- local variables
real :: x_m, x_h
real :: q_m, q_h
! ----------------------------------------------------------------------------
if (zeta >= 0.0) then
q_m = ((1.0 - b_m) / b_m)**(1.0 / 3.0)
q_h = sqrt(c_h * c_h - 4.0)
x_m = (1.0 + zeta)**(1.0 / 3.0)
x_h = zeta
psi_m = -3.0 * (a_m / b_m) * (x_m - 1.0) + 0.5 * (a_m / b_m) * q_m * (&
2.0 * log((x_m + q_m) / (1.0 + q_m)) - &
log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + &
2.0 * sqrt(3.0) * (&
atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - &
atan((2.0 - q_m) / (sqrt(3.0) * q_m))))
psi_h = -0.5 * b_h * log(1.0 + c_h * x_h + x_h * x_h) + &
((-a_h / q_h) + ((b_h * c_h) / (2.0 * q_h))) * (&
log((2.0 * x_h + c_h - q_h) / (2.0 * x_h + c_h + q_h)) - &
log((c_h - q_h) / (c_h + q_h)))
else
x_m = (1.0 - alpha_m * zeta)**(0.25)
x_h = (1.0 - alpha_h * zeta)**(0.25)
psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m)
psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h))
end if
end subroutine
subroutine get_psi_mh(psi_m, psi_h, zeta_m, zeta_h)
!< @brief universal functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions
real, intent(in) :: zeta_m, zeta_h !< = z/L
! ----------------------------------------------------------------------------
! --- local variables
real :: x_m, x_h
real :: q_m, q_h
! ----------------------------------------------------------------------------
if (zeta_m >= 0.0) then
q_m = ((1.0 - b_m) / b_m)**(1.0 / 3.0)
x_m = (1.0 + zeta_m)**(1.0 / 3.0)
psi_m = -3.0 * (a_m / b_m) * (x_m - 1.0) + 0.5 * (a_m / b_m) * q_m * (&
2.0 * log((x_m + q_m) / (1.0 + q_m)) - &
log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + &
2.0 * sqrt(3.0) * (&
atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - &
atan((2.0 - q_m) / (sqrt(3.0) * q_m))))
else
x_m = (1.0 - alpha_m * zeta_m)**(0.25)
psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m)
end if
if (zeta_h >= 0.0) then
q_h = sqrt(c_h * c_h - 4.0)
x_h = zeta_h
psi_h = -0.5 * b_h * log(1.0 + c_h * x_h + x_h * x_h) + &
((-a_h / q_h) + ((b_h * c_h) / (2.0 * q_h))) * (&
log((2.0 * x_h + c_h - q_h) / (2.0 * x_h + c_h + q_h)) - &
log((c_h - q_h) / (c_h + q_h)))
else
x_h = (1.0 - alpha_h * zeta_h)**(0.25)
psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h))
end if
end subroutine
! --------------------------------------------------------------------------------
end module sfx_sheba