noniterative sheba scheme for stable surface layer

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

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

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]'

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

srcF/sfx_sheba_noit_param.f90

+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

srcF/sfx_sheba_noniterative.f90

+#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