! Created by Andrey Debolskiy on 24.01.2025.
module fo_stability_functions
    !< Constants to for stability functions
    !< louis stability functions
    real, parameter:: yb = 5.e0
    real, parameter:: yc = 5.e0
    real, parameter:: yd = 5.e0
    real, parameter:: ye = yb*yc*sqrt(1.e0/3.e0)
    real, parameter::ricr = 2.e0/(3.e0*yd)
    real, parameter:: vlinfm = 120.0e0
    real, parameter:: vlinfh = vlinfm * sqrt(3.e0*yd/2.e0)


    !    security parameters
    !    du2min is a minimal value to calculate wind shear
    !   cormin is a minimal value of coriolis parameter to calculate the
    !   pbl extention
    real, parameter::  du2min = 1.e-02
    real, parameter::  cormin = 5.e-05
    contains

    !>
    subroutine louis_stab_functions(fnriuv,fnritq,rinum,vdlm, &
            & vdlh,zkh, rolim)
        implicit none
        real fnriuv,fnritq
        real, intent(in):: rinum
        real, intent(in):: vdlm,vdlh, zkh
        real scf, ucfm, ucfh
        !real yb, yc, yd, ye, ricr
        real, intent(in):: rolim

        if(rinum.ge.0.e0) then
            scf = sqrt(1.e0 + yd*rinum)
            fnriuv = 1.e0/(1.e0 + 2.e0*yb * rinum / scf)
            fnritq = 1.e0/(1.e0 + 3.e0*yb * rinum * scf)
        else
            ucfm = 1.e0/(1.e0 + ye * (vdlm/zkh)**2 * sqrt(abs(rinum)))
            ucfh = 1.e0/(1.e0 + ye * (vdlh/zkh)**2 * sqrt(abs(rinum)))
            fnriuv = 1.e0 - 2.e0*yb * rinum * ucfm
            if(rolim.gt.0.99.and.rolim.lt.1.5) then
                fnritq = 1.e0 - 3.0*3.e0*yb * rinum * ucfh    !enhanced mixing for unstable stratification over land
            else
                fnritq = 1.e0 - 3.e0*yb * rinum * ucfh
            end if
        end if
    end subroutine louis_stab_functions

    !>
    !>esau byrkjedal 2007  +   viterbo 1999  unstable stratification
    subroutine esau_vit_stab_functions( fnriuv,fnritq,rinum, &
            &   vdlm,zkh, dz)
        implicit none
        real fnriuv,fnritq,rinum, vdlm, zkh, dz
        real scf, ucfm, ucfh
        real yb, yc, yd, ye, ricr
        real gh
        real, parameter:: p_m = 21.0
        real, parameter:: q_m = 0.005
        real, parameter:: p_h = 10.0
        real, parameter:: q_h = 0.0012
        real, parameter:: bb = 5.0
        real, parameter:: cc = 5.0

        yb = 5.e0
        yc = 5.e0
        yd = 5.e0
        ye = yb*yc*sqrt(1.e0/3.e0)
        ricr = 2.e0/(3.e0*yd)


        if(rinum.ge.0.e0) then
            scf = sqrt(1.e0 + yd*rinum)
            fnriuv = 1.0/((1.0 + p_m * rinum)*(1.0 + p_m * rinum)) &
                    &   + q_m * sqrt(rinum)
            fnritq = 1.0/((1.0+p_h*rinum)*(1.0+p_h*rinum)*(1.0+p_h*rinum)) &
                    &  + q_h
        else
            gh = ((vdlm * vdlm)/(dz * sqrt(abs(dz * zkh)))) &
                    &     * ((1 + dz / zkh)**(1/3) - 1)**(3/2)

            fnriuv =  1.0 - (2.0 * bb * rinum) &
                    &  / (1.0 + 3.0 * bb * cc * gh * sqrt(abs(rinum)))

            fnritq = 1.0 - (2.0 * bb * rinum) &
                    &  / (1.0 + 3.0 * bb * cc * gh * sqrt(abs(rinum)))
        end if
    end subroutine esau_vit_stab_functions


end module fo_stability_functions