Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
! 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