#include "../includeF/sfx_def.fi"
module sfx_esm
!< @brief main Earth System Model surface flux module
! modules used
! --------------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
use sfx_common
#endif
use sfx_data
use sfx_surface
use sfx_esm_param
#if defined(INCLUDE_CXX)
use iso_c_binding, only: C_LOC, C_PTR
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_convection = 10 !< maximum (actual) number of iterations in convection
integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness
end type
! --------------------------------------------------------------------------------
contains
! --------------------------------------------------------------------------------
#if defined(INCLUDE_CXX)
subroutine set_c_struct_sfx_esm_param_values(sfx_model_param)
type (sfx_esm_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%Pr_t_inf_inv = Pr_t_inf_inv
sfx_model_param%alpha_m = alpha_m
sfx_model_param%alpha_h = alpha_h
sfx_model_param%alpha_h_fix = alpha_h_fix
sfx_model_param%beta_m = beta_m
sfx_model_param%beta_h = beta_h
sfx_model_param%Rib_max = Rib_max
end subroutine set_c_struct_sfx_esm_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), target :: meteo_c !< meteorological data (input)
type (sfxDataVecTypeC), target :: sfx_c !< surface flux data (output)
type(C_PTR) :: meteo_c_ptr, sfx_c_ptr
type (sfx_esm_param) :: model_param
type (sfx_surface_param) :: surface_param
type (sfx_esm_numericsTypeC) :: numerics_c
type (sfx_phys_constants) :: phys_constants
numerics_c%maxiters_convection = numerics%maxiters_convection
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_esm_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)
meteo_c_ptr = C_LOC(meteo_c)
sfx_c_ptr = C_LOC(sfx_c)
call get_surface_fluxes_esm(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n)
#else
do i = 1, n
#ifdef SFX_FORCE_DEPRECATED_ESM_CODE
#else
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)
#endif
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 zeta_conv_lim !< z/L critical value for matching free convection limit [n/d]
real Rib_conv_lim !< Ri-bulk critical value for matching free convection limit [n/d]
real f_m_conv_lim !< stability function (momentum) value in free convection limit [n/d]
real f_h_conv_lim !< stability function (heat) value in free convection limit [n/d]
real psi_m, psi_h !< universal functions (momentum) & (heat) [n/d]
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)
real fval !< just a shortcut for partial calculations
#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
! --- define free convection transition zeta = z/L value
call get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, &
h0_m, h0_t, B)
! --- get the fluxes
! ----------------------------------------------------------------------------
if (Rib > 0.0) then
! --- stable stratification block
! --- restrict bulk Ri value
! *: note that this value is written to output
Rib = min(Rib, Rib_max)
call get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
fval = beta_m * zeta
phi_m = 1.0 + fval
phi_h = 1.0/Pr_t_0_inv + fval
else if (Rib < Rib_conv_lim) then
! --- strong instability block
call get_psi_convection(psi_m, psi_h, zeta, Rib, &
zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, h0_m, h0_t, B, numerics%maxiters_convection)
fval = (zeta_conv_lim / zeta)**(1.0/3.0)
phi_m = fval / f_m_conv_lim
phi_h = fval / (Pr_t_0_inv * f_h_conv_lim)
else if (Rib > -0.001) then
! --- nearly neutral [-0.001, 0] block
call get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
zeta = 0.0
phi_m = 1.0
phi_h = 1.0 / Pr_t_0_inv
else
! --- weak & semistrong instability block
call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection)
phi_m = (1.0 - alpha_m * zeta)**(-0.25)
phi_h = 1.0 / (Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta))
end if
! ----------------------------------------------------------------------------
! --- define transfer coeff. (momentum) & (heat)
Cm = kappa / psi_m
Ct = kappa / psi_h
! --- 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 = Rib_conv_lim, &
Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
! convection universal functions shortcuts
! --------------------------------------------------------------------------------
function f_m_conv(zeta)
! ----------------------------------------------------------------------------
real :: f_m_conv
real, intent(in) :: zeta
! ----------------------------------------------------------------------------
f_m_conv = (1.0 - alpha_m * zeta)**0.25
end function f_m_conv
function f_h_conv(zeta)
! ----------------------------------------------------------------------------
real :: f_h_conv
real, intent(in) :: zeta
! ----------------------------------------------------------------------------
f_h_conv = (1.0 - alpha_h * zeta)**0.5
end function f_h_conv
! --------------------------------------------------------------------------------
! universal functions
! --------------------------------------------------------------------------------
subroutine get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
!< @brief universal functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions
real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h
real, intent(in) :: B !< = log(z0_m / z0_h)
! ----------------------------------------------------------------------------
psi_m = log(h0_m)
psi_h = log(h0_t) / Pr_t_0_inv
!*: this looks redundant z0_t = z0_m in case |B| ~ 0
if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv
end subroutine
subroutine get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
!< @brief universal functions (momentum) & (heat): stable case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions [n/d]
real, intent(out) :: zeta !< = z/L [n/d]
real, intent(in) :: Rib !< bulk Richardson number [n/d]
real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d]
real, intent(in) :: B !< = log(z0_m / z0_h) [n/d]
! ----------------------------------------------------------------------------
! --- local variables
real :: Rib_coeff
real :: psi0_m, psi0_h
real :: phi
real :: c
! ----------------------------------------------------------------------------
psi0_m = alog(h0_m)
psi0_h = B / psi0_m
Rib_coeff = beta_m * Rib
c = (psi0_h + 1.0) / Pr_t_0_inv - 2.0 * Rib_coeff
zeta = psi0_m * (sqrt(c**2 + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / &
(2.0 * beta_m * (1.0 - Rib_coeff))
phi = beta_m * zeta
psi_m = psi0_m + phi
psi_h = (psi0_m + B) / Pr_t_0_inv + phi
end subroutine
subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters)
!< @brief universal functions (momentum) & (heat): semi-strong convection case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions [n/d]
real, intent(out) :: zeta !< = z/L [n/d]
real, intent(in) :: Rib !< bulk Richardson number [n/d]
real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d]
real, intent(in) :: B !< = log(z0_m / z0_h) [n/d]
integer, intent(in) :: maxiters !< maximum number of iterations
! --- local variables
real :: zeta0_m, zeta0_h
real :: f0_m, f0_h
real :: f_m, f_h
integer :: i
! ----------------------------------------------------------------------------
psi_m = log(h0_m)
psi_h = log(h0_t)
if (abs(B) < 1.0e-10) psi_h = psi_m
zeta = Rib * Pr_t_0_inv * psi_m**2 / psi_h
do i = 1, maxiters + 1
zeta0_m = zeta / h0_m
zeta0_h = zeta / h0_t
if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
f_m = (1.0 - alpha_m * zeta)**0.25e0
f_h = sqrt(1.0 - alpha_h_fix * zeta)
f0_m = (1.0 - alpha_m * zeta0_m)**0.25e0
f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h)
f0_m = max(f0_m, 1.000001e0)
f0_h = max(f0_h, 1.000001e0)
psi_m = log((f_m - 1.0e0)*(f0_m + 1.0e0)/((f_m + 1.0e0)*(f0_m - 1.0e0))) + 2.0e0*(atan(f_m) - atan(f0_m))
psi_h = log((f_h - 1.0e0)*(f0_h + 1.0e0)/((f_h + 1.0e0)*(f0_h - 1.0e0))) / Pr_t_0_inv
! *: don't update zeta = z/L at last iteration
if (i == maxiters + 1) exit
zeta = Rib * psi_m**2 / psi_h
end do
end subroutine
subroutine get_psi_convection(psi_m, psi_h, zeta, Rib, &
zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, &
h0_m, h0_t, B, maxiters)
!< @brief universal functions (momentum) & (heat): fully convective case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions [n/d]
real, intent(out) :: zeta !< = z/L [n/d]
real, intent(in) :: Rib !< bulk Richardson number [n/d]
real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d]
real, intent(in) :: B !< = log(z0_m / z0_h) [n/d]
integer, intent(in) :: maxiters !< maximum number of iterations
real, intent(in) :: zeta_conv_lim !< convective limit zeta
real, intent(in) :: f_m_conv_lim, f_h_conv_lim !< universal function shortcuts in limiting case
! ----------------------------------------------------------------------------
! --- local variables
real :: zeta0_m, zeta0_h
real :: f0_m, f0_h
real :: p_m, p_h
real :: a_m, a_h
real :: c_lim, f
integer :: i
! ----------------------------------------------------------------------------
p_m = 2.0 * atan(f_m_conv_lim) + log((f_m_conv_lim - 1.0) / (f_m_conv_lim + 1.0))
p_h = log((f_h_conv_lim - 1.0) / (f_h_conv_lim + 1.0))
zeta = zeta_conv_lim
do i = 1, maxiters + 1
zeta0_m = zeta / h0_m
zeta0_h = zeta / h0_t
if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
f0_m = (1.0 - alpha_m * zeta0_m)**0.25
f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h)
a_m = -2.0*atan(f0_m) + log((f0_m + 1.0)/(f0_m - 1.0))
a_h = log((f0_h + 1.0)/(f0_h - 1.0))
c_lim = (zeta_conv_lim / zeta)**(1.0 / 3.0)
f = 3.0 * (1.0 - c_lim)
psi_m = f / f_m_conv_lim + p_m + a_m
psi_h = (f / f_h_conv_lim + p_h + a_h) / Pr_t_0_inv
! *: don't update zeta = z/L at last iteration
if (i == maxiters + 1) exit
zeta = Rib * psi_m**2 / psi_h
end do
end subroutine
! --------------------------------------------------------------------------------
! convection limit definition
! --------------------------------------------------------------------------------
subroutine get_convection_lim(zeta_lim, Rib_lim, f_m_lim, f_h_lim, &
h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: zeta_lim !< limiting value of z/L
real, intent(out) :: Rib_lim !< limiting value of Ri-bulk
real, intent(out) :: f_m_lim, f_h_lim !< limiting values of universal functions shortcuts
real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d]
real, intent(in) :: B !< = log(z0_m / z0_h) [n/d]
! ----------------------------------------------------------------------------
! --- local variables
real :: psi_m, psi_h
real :: f_m, f_h
real :: c
! ----------------------------------------------------------------------------
! --- define limiting value of zeta = z / L
c = (Pr_t_inf_inv / Pr_t_0_inv)**4
zeta_lim = (2.0 * alpha_h - c * alpha_m - &
sqrt((c * alpha_m)**2 + 4.0 * c * alpha_h * (alpha_h - alpha_m))) / (2.0 * alpha_h**2)
f_m_lim = f_m_conv(zeta_lim)
f_h_lim = f_h_conv(zeta_lim)
! --- universal functions
f_m = zeta_lim / h0_m
f_h = zeta_lim / h0_t
if (abs(B) < 1.0e-10) f_h = f_m
f_m = (1.0 - alpha_m * f_m)**0.25
f_h = sqrt(1.0 - alpha_h_fix * f_h)
psi_m = 2.0 * (atan(f_m_lim) - atan(f_m)) + alog((f_m_lim - 1.0) * (f_m + 1.0)/((f_m_lim + 1.0) * (f_m - 1.0)))
psi_h = alog((f_h_lim - 1.0) * (f_h + 1.0)/((f_h_lim + 1.0) * (f_h - 1.0))) / Pr_t_0_inv
! --- bulk Richardson number
Rib_lim = zeta_lim * psi_h / (psi_m * psi_m)
end subroutine
! --------------------------------------------------------------------------------
end module sfx_esm