diff --git a/CMakeLists.txt b/CMakeLists.txt
index 39e636a6c21b711c709aa4864de7c52b3f2ccadf..8dcb6779bf338a2fbc79d440745086aa2ba7b7a4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -88,6 +88,7 @@ set(SOURCES_F
srcF/sfx_most.f90
srcF/sfx_most_param.f90
srcF/sfx_sheba.f90
+ srcF/sfx_sheba_noniterative.f90
srcF/sfx_sheba_param.f90
srcF/sfx_fc_wrapper.F90
srcF/sfx_api_inmcm.f90
diff --git a/srcF/sfx_config.f90 b/srcF/sfx_config.f90
index 99f7b22e57e362426e34c9f20a0275bc9ddef4cb..3fde2a0995da10b521591349c02afe7ba8e98d57 100644
--- a/srcF/sfx_config.f90
+++ b/srcF/sfx_config.f90
@@ -19,11 +19,13 @@ module sfx_config
integer, parameter :: model_log = 1 !< LOG simplified model
integer, parameter :: model_most = 2 !< MOST model
integer, parameter :: model_sheba = 3 !< SHEBA model
+ integer, parameter :: model_sheba_noit = 4 !< SHEBA model noniterative
character(len = 16), parameter :: model_esm_tag = 'esm'
character(len = 16), parameter :: model_log_tag = 'log'
character(len = 16), parameter :: model_most_tag = 'most'
character(len = 16), parameter :: model_sheba_tag = 'sheba'
+ character(len = 16), parameter :: model_sheba_noit_tag = 'sheba_noit'
!> @brief dataset enum def.
integer, parameter :: dataset_mosaic = 1 !< MOSAiC campaign
@@ -68,6 +70,8 @@ contains
id = model_most
else if (trim(tag) == trim(model_sheba_tag)) then
id = model_sheba
+ else if (trim(tag) == trim(model_sheba_noit_tag)) then
+ id = model_sheba_noit
end if
end function
@@ -86,6 +90,8 @@ contains
tag = model_most_tag
else if (id == model_sheba) then
tag = model_sheba_tag
+ else if (id == model_sheba_noit) then
+ tag = model_sheba_noit_tag
end if
end function
diff --git a/srcF/sfx_run.f90 b/srcF/sfx_run.f90
index 03d2a9ee4b8ba3ef668abf64af31cedb66001823..eec3f6ea7eb0196907298ea8892b6a21113a3482 100644
--- a/srcF/sfx_run.f90
+++ b/srcF/sfx_run.f90
@@ -40,6 +40,9 @@ contains
use sfx_sheba, only: &
get_surface_fluxes_vec_sheba => get_surface_fluxes_vec, &
numericsType_sheba => numericsType
+ use sfx_sheba_noniterative, only: &
+ get_surface_fluxes_vec_sheba_noit => get_surface_fluxes_vec, &
+ numericsType_sheba_noit => numericsType
! --------------------------------------------------------------------------------
character(len=*), intent(in) :: filename_out
@@ -59,6 +62,7 @@ contains
type(numericsType_log) :: numerics_log !< surface flux module (LOG) numerics parameters
type(numericsType_most) :: numerics_most !< surface flux module (MOST) numerics parameters
type(numericsType_sheba) :: numerics_sheba !< surface flux module (SHEBA) numerics parameters
+ type(numericsType_sheba_noit) :: numerics_sheba_noit !< surface flux module (SHEBA) numerics parameters
integer :: num !< number of 'cells' in input
! --------------------------------------------------------------------------------
@@ -147,6 +151,8 @@ contains
call get_surface_fluxes_vec_most(sfx, meteo, numerics_most, num)
else if (model == model_sheba) then
call get_surface_fluxes_vec_sheba(sfx, meteo, numerics_sheba, num)
+ else if (model == model_sheba_noit) then
+ call get_surface_fluxes_vec_sheba_noit(sfx, meteo, numerics_sheba_noit, num)
end if
@@ -233,7 +239,7 @@ contains
write(*, *) ' --help'
write(*, *) ' print usage options'
write(*, *) ' --model [key]'
- write(*, *) ' key = esm (default) || log || most || sheba'
+ write(*, *) ' key = esm (default) || log || most || sheba || sheba_noit'
write(*, *) ' --dataset [key]'
write(*, *) ' key = mosaic (default) || irgason || sheba'
write(*, *) ' = lake || papa || toga || user [filename] [param]'
diff --git a/srcF/sfx_sheba_noniterative.f90 b/srcF/sfx_sheba_noniterative.f90
new file mode 100644
index 0000000000000000000000000000000000000000..dad747ae09fef45ede28a4053159db27e0b1433f
--- /dev/null
+++ b/srcF/sfx_sheba_noniterative.f90
@@ -0,0 +1,514 @@
+#include "../includeF/sfx_def.fi"
+
+module sfx_sheba_noniterative
+ !< @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_CXX)
+ use iso_c_binding, only: C_LOC, C_PTR, C_INT, C_FLOAT
+ use C_FUNC
+#endif
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ private
+ ! --------------------------------------------------------------------------------
+
+ ! public interface
+ ! --------------------------------------------------------------------------------
+ public :: get_surface_fluxes
+ public :: get_surface_fluxes_vec
+ public :: get_psi
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ type, public :: numericsType
+ integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness
+ end type
+ ! --------------------------------------------------------------------------------
+
+#if defined(INCLUDE_CXX)
+ type, BIND(C), public :: sfx_sheba_param_C
+ real(C_FLOAT) :: kappa
+ real(C_FLOAT) :: Pr_t_0_inv
+
+ real(C_FLOAT) :: alpha_m
+ real(C_FLOAT) :: alpha_h
+ real(C_FLOAT) :: a_m
+ real(C_FLOAT) :: b_m
+ real(C_FLOAT) :: a_h
+ real(C_FLOAT) :: b_h
+ real(C_FLOAT) :: c_h
+ end type
+
+ type, BIND(C), public :: sfx_sheba_numericsType_C
+ integer(C_INT) :: maxiters_charnock
+ end type
+
+ INTERFACE
+ SUBROUTINE c_sheba_compute_flux(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C, &
+ name="c_sheba_compute_flux")
+ use sfx_data
+ use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR
+ Import :: sfx_sheba_param_C, sfx_sheba_numericsType_C
+ implicit none
+ INTEGER(C_INT) :: grid_size
+ type(C_PTR), value :: sfx
+ type(C_PTR), value :: meteo
+ type(sfx_sheba_param_C) :: model_param
+ type(sfx_surface_param) :: surface_param
+ type(sfx_sheba_numericsType_C) :: numerics
+ type(sfx_phys_constants) :: constants
+ END SUBROUTINE c_sheba_compute_flux
+ END INTERFACE
+#endif
+
+contains
+
+#if defined(INCLUDE_CXX)
+ subroutine set_c_struct_sfx_sheba_param_values(sfx_model_param)
+ type (sfx_sheba_param_C), 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), 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_sheba_param_C) :: model_param
+ type (sfx_surface_param) :: surface_param
+ type (sfx_sheba_numericsType_C) :: 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)
+
+ meteo_c_ptr = C_LOC(meteo_c)
+ sfx_c_ptr = C_LOC(sfx_c)
+
+ call c_sheba_compute_flux(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n)
+#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
+ ! ----------------------------------------------------------------------------
+ if(Rib > 0)then
+ call get_dynamic_scales_noniterative(Udyn, Tdyn, Qdyn, zeta, &
+ U, dT, dQ, h, z0_m, z0_t, Rib)
+ else
+ call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
+ U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10)
+ end if
+
+ ! ----------------------------------------------------------------------------
+
+ 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
+
+
+ subroutine get_dynamic_scales_noniterative(Udyn, Tdyn, Qdyn, zeta, &
+ U, dT, dQ, z, z0_m, z0_t, Rib)
+ ! ----------------------------------------------------------------------------
+ real, parameter :: gamma = 2.91, zeta_a = 3.6
+
+ 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) :: 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) :: Rib !< bulk Richardson number
+
+ ! ----------------------------------------------------------------------------
+
+ ! --- local variables
+ real :: psi_m, psi_h
+ real :: psi0_m, psi0_h
+ real :: C1,A1,A2,lne,lnet,Ribl
+ ! ----------------------------------------------------------------------------
+
+ Ribl = (Rib*Pr_t_0_inv) * (1 - z0_t / z) / ((1 - z0_m / z)**2)
+
+ call get_psi(psi_m, psi_h, zeta_a)
+ call get_psi_mh(psi0_m, psi0_h, zeta_a * z0_m / z, zeta_a * z0_t / z)
+
+ lne = log(z/z0_m)
+ lnet = log(z/z0_t)
+ C1 = (lne**2)/lnet
+ A1 = ((lne - psi_m + psi0_m)**(2*(gamma-1))) &
+& / ((zeta_a**(gamma-1))*((lnet-(psi_h-psi0_h)*Pr_t_0_inv)**(gamma-1)))
+ A2 = ((lne - psi_m + psi0_m)**2) / (lnet-(psi_h-psi0_h)*Pr_t_0_inv) - C1
+
+ zeta = C1 * Ribl + A1 * A2 * (Ribl**gamma)
+
+ call get_psi(psi_m, psi_h, zeta)
+ call get_psi_mh(psi0_m, psi0_h, zeta * z0_m / z, zeta * z0_t /z)
+
+ 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))
+
+
+ end subroutine get_dynamic_scales_noniterative
+ ! --------------------------------------------------------------------------------
+
+ ! 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_noniterative
\ No newline at end of file