#include "../includeF/sfx_def.fi"

module sfx_sheba_noniterative

    ! 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
    public :: get_psi_stable
    integer z0m_id
    integer z0t_id
    ! --------------------------------------------------------------------------------

    ! --------------------------------------------------------------------------------
    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
        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_noit_numericsType_C 
        integer(C_INT) :: maxiters_charnock 
        integer(C_INT) :: maxiters_convection
    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_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
        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_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

            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), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
            !write(*,*) 'get_flux', meteo%surface_type(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)
        real :: depth   
        real :: lai
        integer :: surface_type
        ! ----------------------------------------------------------------------------

        ! --- 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 Rib, Ri_sn         !< bulk Richardson number, snow Richardson number


        real zeta               !< = z/L [n/d]
        real Re                 !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]

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

        

        real fval               !< just a shortcut for partial calculations

        real w_snow, deltaS, h_salt
        real sigma_w, sigma_r
        real S_mean, S_salt, Linv
        real z0_m1
        integer i
#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_m1 = meteo%z0_m
        depth = meteo%depth
        lai = meteo%lai
        surface_type=meteo%surface_type
       
        
        !write (*,*) surface_type, 'esm'
        call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
        forest_z0m_id, usersf_z0m_id, ice_z0m_id, z0m_id)
        
        
        call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m1, z0m_id)
        
        call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
        forest_z0t_id, usersf_z0t_id, ice_z0t_id, z0t_id)
        
        
        Re = u_dyn0 * z0_m / nu_air

        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
        ! --- define relative height
            h0_m = h / z0_m
        ! --- 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)
        
        !  --- define snow consentration
        !call get_sigma(sigma_r, sigma_w, rho_air, rho_s)
        !call get_w_snow(w_snow, sigma_w, g, d_s, nu_air)
        !call get_h_salt(h_salt, u_dyn0)
        !call get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt)
        !call get_S_mean(S_mean,  S_salt, h_salt, h, w_snow, u_dyn0)
                      
        !deltaS=S_salt-S_mean
    
        !Ri_sn = (g * sigma_r * deltaS * h) / U**2
                !Ri_sn=0.0
        
        ! --- 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)
        !do i = 1, numerics%maxiters_convection  
        !if (surface_type == surface_snow) then    
            !write(*,*) 'RIsnow', Rib, Ri_sn   
         !   Rib=Rib+Ri_sn 
        !endif

            call get_zeta(zeta, Rib, h, z0_m, z0_t)
            !write(*,*) 'get_zeta', zeta, h
            
            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)
            S_mean=0.0
            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))
            !write(*,*) 'sfx_before_snow', Udyn, zeta, S_mean

            if (surface_type==3.or.surface_type==6.) then
                if (Udyn>u_thsnow) then
            
            call get_sigma(sigma_r, sigma_w, rho_air, rho_s)
            call get_w_snow(w_snow, sigma_w, g, d_s, nu_air)
            call get_h_salt(h_salt, Udyn)
            call get_S_salt(S_salt, Udyn, u_thsnow, g, h_salt)
            call get_S_mean(S_mean,  S_salt, h_salt, h, w_snow, Udyn)
                      
            Linv=Linv*((1-S_mean)/(1+sigma_w*S_mean))+(g*w_snow*sigma_w*S_mean/(Udyn**3.0))/(1+sigma_w*S_mean)
            zeta = h * Linv 
            !write(*,*) 'sfx_snow1', Udyn, Linv, zeta, S_mean
                  
                 
             
            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))

            write(*,*) 'sfx_snow2', Udyn, zeta, S_mean, Linv
            endif 
            endif
        
           ! Ri_sn = (g * sigma_r * deltaS * h) / U**2
            
            !end do
            
        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)
        !write(*,*) 'Smean_0ut', S_mean
        sfx = sfxDataType(zeta = zeta, Rib = Rib, &
            Re = Linv, B = B, z0_m = z0_m, z0_t = z0_t, &
            Rib_conv_lim = S_mean, &
            Cm = Cm, Ct = Ct, Km = S_salt, Pr_t_inv = Udyn)

    end subroutine get_surface_fluxes
    ! --------------------------------------------------------------------------------
    !--------------------------------------snow functions-----------------------------
    subroutine get_sigma(sigma_r,  sigma_w, rho_air, rho_s)
        !< @brief function for 
        ! ----------------------------------------------------------------------------
        real, intent(out) :: sigma_r, sigma_w  !< s
        real, intent(in) ::  rho_air, rho_s            
        ! ----------------------------------------------------------------------------
        sigma_r  = ((rho_s/rho_air)-1)
        sigma_w  = (rho_s-rho_air)/rho_air
    end subroutine

    subroutine get_h_salt(h_salt, u_dyn0)
        ! ----------------------------------------------------------------------------
        real, intent(out) :: h_salt   
        real, intent(in) :: u_dyn0
        ! ----------------------------------------------------------------------------
        h_salt  = 0.08436*u_dyn0**1.27
    end subroutine

    
    subroutine get_S_mean(S_mean, S_salt, h_salt, z, w_snow, u_dyn0)
        !< @brief function for snow consentration
        ! ----------------------------------------------------------------------------
        real, intent(out) :: S_mean   !< snow consentration
        real, intent(in) ::  S_salt, h_salt, z, w_snow, u_dyn0         
        ! ----------------------------------------------------------------------------
        S_mean  = S_salt *  (z/h_salt)**(-w_snow/(kappa*u_dyn0))
        !write(*,*) 'Smean_func', S_mean, u_dyn0
    end subroutine

    subroutine get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt)
        !< @brief function for snow consentration
        ! ----------------------------------------------------------------------------
        real, intent(out) :: S_salt   !< snow consentration
        real, intent(in) ::  u_dyn0,u_thsnow, g, h_salt         
        ! ----------------------------------------------------------------------------
        real qs
        if (u_dyn0>u_thsnow) then
            qs  = (u_dyn0*u_dyn0-u_thsnow*u_thsnow)/(Csn*u_dyn0*g*h_salt)
            S_salt=qs/(qs+rho_s/rho_air)
            !write(*,*) 'S_salt', S_salt
        else
            S_salt=0.0
        endif
    end subroutine

    subroutine get_w_snow(w_snow, sigma_w, g, d_s, nu_air)
        !< @brief function for smow velosity
        ! ----------------------------------------------------------------------------
        real, intent(out) :: w_snow   !< 
        real, intent(in) ::  sigma_w, g, d_s, nu_air            
        ! ----------------------------------------------------------------------------
        w_snow  = (sigma_w * g * d_s * d_s) / (18.0 * nu_air);

    end subroutine
         !-----------------------------------------------------------------------------




    subroutine get_zeta(zeta, Rib, h, z0_m, z0_t)
    real,intent(out) :: zeta
    real,intent(in) :: Rib, h, z0_m, z0_t

    real :: Ribl, C1, A1, A2, lne, lnet
    real :: psi_m, psi_h, psi0_m, psi0_h

        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)

    end subroutine get_zeta

    ! 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 get_psi_stable


    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