diff --git a/CMakeLists.txt b/CMakeLists.txt
index 39e636a6c21b711c709aa4864de7c52b3f2ccadf..9ed5741ad1ed7d36301778922b5eab91f69451c9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -89,6 +89,8 @@ set(SOURCES_F
srcF/sfx_most_param.f90
srcF/sfx_sheba.f90
srcF/sfx_sheba_param.f90
+ srcF/sfx_sheba_noniterative.f90
+ srcF/sfx_sheba_noit_param.f90
srcF/sfx_fc_wrapper.F90
srcF/sfx_api_inmcm.f90
srcF/sfx_api_term.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.f90 b/srcF/sfx_sheba.f90
index c33e1a5144fb948ed18a070c835cbdd8319f3f60..4fd6142bb7ddd9a16e553d54e5b075a6e43d6f36 100644
--- a/srcF/sfx_sheba.f90
+++ b/srcF/sfx_sheba.f90
@@ -33,6 +33,7 @@ module sfx_sheba
! --------------------------------------------------------------------------------
type, public :: numericsType
+ integer :: maxiters_convection = 10 !< maximum (actual) number of iterations
integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness
end type
! --------------------------------------------------------------------------------
@@ -40,7 +41,8 @@ module sfx_sheba
#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) ::
+ Pr_t_0_inv
real(C_FLOAT) :: alpha_m
real(C_FLOAT) :: alpha_h
@@ -256,7 +258,7 @@ contains
! --- get the fluxes
! ----------------------------------------------------------------------------
call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
- U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10)
+ U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), numerics%maxiters_convection)
! ----------------------------------------------------------------------------
call get_phi(phi_m, phi_h, zeta)
@@ -343,6 +345,9 @@ contains
end do
end subroutine get_dynamic_scales
+
+
+
! --------------------------------------------------------------------------------
! stability functions
diff --git a/srcF/sfx_sheba_noit_param.f90 b/srcF/sfx_sheba_noit_param.f90
new file mode 100644
index 0000000000000000000000000000000000000000..6461dde50e1fa347b7f1cfe86934cdc22b3be970
--- /dev/null
+++ b/srcF/sfx_sheba_noit_param.f90
@@ -0,0 +1,40 @@
+module sfx_sheba_noit_param
+ !< @brief SHEBA surface flux model parameters
+ !< @details all in SI units
+
+ ! modules used
+ ! --------------------------------------------------------------------------------
+ use sfx_phys_const
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ ! --------------------------------------------------------------------------------
+
+
+ !< von Karman constant [n/d]
+ real, parameter :: kappa = 0.40
+ !< inverse Prandtl (turbulent) number in neutral conditions [n/d]
+ real, parameter :: Pr_t_0_inv = 1.15
+ !< inverse Prandtl (turbulent) number in free convection [n/d]
+ real, parameter :: Pr_t_inf_inv = 3.5
+
+ real, parameter :: Rib_max = 0.5
+
+ !< stability function coeff. (unstable)
+ real, parameter :: alpha_m = 16.0
+ real, parameter :: alpha_h = 16.0
+ real, parameter :: alpha_h_fix = 16.0
+
+ !< stability function coeff. (stable)
+ real, parameter :: a_m = 5.0
+ real, parameter :: b_m = a_m / 6.5
+ real, parameter :: a_h = 5.0
+ real, parameter :: b_h = 5.0
+ real, parameter :: c_h = 3.0
+
+ real, parameter :: gamma = 2.91, zeta_a = 3.6 ! for stable psi
+
+
+end module sfx_sheba_noit_param
diff --git a/srcF/sfx_sheba_noniterative.f90 b/srcF/sfx_sheba_noniterative.f90
new file mode 100644
index 0000000000000000000000000000000000000000..19d01aa62da6109ffdf7831015a05fbd78be5489
--- /dev/null
+++ b/srcF/sfx_sheba_noniterative.f90
@@ -0,0 +1,588 @@
+#include "../includeF/sfx_def.fi"
+
+module sfx_sheba_noniterative
+ !< @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_sheba_noit_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
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ 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
+ ! --------------------------------------------------------------------------------
+
+#if defined(INCLUDE_CXX)
+ type, BIND(C), public :: sfx_sheba_noit_param_C
+ real(C_FLOAT) :: kappa
+ real(C_FLOAT) :: Pr_t_0_inv
+ real(C_FLOAT) :: Pr_t_inf_inv
+
+ real(C_FLOAT) :: alpha_m
+ real(C_FLOAT) :: alpha_h
+ real(C_FLOAT) :: alpha_h_fix
+ real(C_FLOAT) :: Rib_max
+ real(C_FLOAT) :: gamma
+ real(C_FLOAT) :: zeta_a
+ end type
+
+ type, BIND(C), public :: sfx_sheba_noit_numericsType_C
+ integer(C_INT) :: maxiters_convection
+ integer(C_INT) :: maxiters_charnock
+ end type
+
+ INTERFACE
+ SUBROUTINE c_sheba_noit_compute_flux(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C, &
+ name="c_sheba_noit_compute_flux")
+ use sfx_data
+ use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR
+ Import :: sfx_sheba_noit_param_C, sfx_sheba_noit_numericsType_C
+ implicit none
+ integer(C_INT) :: grid_size
+ type(C_PTR), value :: sfx
+ type(C_PTR), value :: meteo
+ type(sfx_sheba_noit_param_C) :: model_param
+ type(sfx_surface_noit_param) :: surface_param
+ type(sfx_sheba_noit_numericsType_C) :: numerics
+ type(sfx_phys_constants) :: constants
+ END SUBROUTINE c_sheba_noit_compute_flux
+ END INTERFACE
+#endif
+
+contains
+
+ ! --------------------------------------------------------------------------------
+#if defined(INCLUDE_CXX)
+ subroutine set_c_struct_sfx_sheba_noit_param_values(sfx_model_param)
+ type (sfx_sheba_noit_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%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%Rib_max = Rib_max
+ sfx_model_param%gamma = gamma
+ sfx_model_param%zeta_a = zeta_a
+ end subroutine set_c_struct_sfx_sheba_noit_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_noit_param_C) :: model_param
+ type (sfx_surface_param) :: surface_param
+ type (sfx_sheba_noit_numericsType_C) :: 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_sheba_noit_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_noit_compute_flux(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n)
+#else
+ do i = 1, n
+#ifdef SFX_FORCE_DEPRECATED_sheba_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 psi0_m, psi0_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]
+
+ real Udyn, Tdyn
+
+ integer surface_type !< surface type = (ocean || land)
+
+ real fval !< just a shortcut for partial calculations
+
+ real :: C1,A1,A2,lne,lnet,Ribl
+
+#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)
+
+ Ribl = (Rib*Pr_t_0_inv) * (1 - z0_t / h) / ((1 - z0_m / h)**2)
+
+ call get_psi_stable(psi_m, psi_h, zeta_a, zeta_a)
+ call get_psi_stable(psi0_m, psi0_h, zeta_a * z0_m / h, zeta_a * z0_t / h)
+
+ lne = log(h/z0_m)
+ lnet = log(h/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_stable(psi_m, psi_h, zeta, zeta)
+ call get_psi_stable(psi0_m, psi0_h, zeta * z0_m / h, zeta * z0_t / h)
+
+ 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)
+
+ Udyn = kappa * U / (log(h / z0_m) - (psi_m - psi0_m))
+ Tdyn = kappa * dT * Pr_t_0_inv / (log(h / z0_t) - (psi_h - psi0_h))
+
+ 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)
+ if(Rib > 0)then
+ 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
+ else
+ Cm = kappa / psi_m
+ Ct = kappa / psi_h
+ 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 = 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_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
+ ! ----------------------------------------------------------------------------
+
+
+ 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))))
+ 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)))
+
+ 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_sheba_noniterative
\ No newline at end of file