Skip to content
Snippets Groups Projects
fo_stability_functions.f90 4.54 KiB
Newer Older
  • Learn to ignore specific revisions
  • Debolskiy Andrey's avatar
    Debolskiy Andrey committed
    ! 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, oc_flag)
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            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
    
            integer, intent(in):: oc_flag
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
    
            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(oc_flag > 0) then
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                    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
    
    
        !>esau byrkjedal 2007  +   viterbo 1999  unstable stratification
        subroutine esau_stab_functions( fnriuv,fnritq,rinum, &
                &   vdlm, vdlh, zkh, dz, oc_flag)
            implicit none
            real fnriuv,fnritq,rinum, vdlm, vdlh, zkh, dz
            integer oc_flag
            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
                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(oc_flag >= 1) 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 esau_stab_functions
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
    
    end module fo_stability_functions