Newer
Older
Evgeny Mortikov
committed
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
module sfx_roughness
!> @brief surface roughness parameterizations
! modules used
! --------------------------------------------------------------------------------
use sfx_phys_const
! --------------------------------------------------------------------------------
! directives list
! --------------------------------------------------------------------------------
implicit none
private
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
real, parameter :: kappa = 0.40 !> von Karman constant [n/d]
! --------------------------------------------------------------------------------
! public interface
! --------------------------------------------------------------------------------
public :: get_charnock_roughness
! --------------------------------------------------------------------------------
!> Charnock parameters
!> z0 = Re_visc_min * (nu / u_dyn) + gamma_c * (u_dyn^2 / g)
! --------------------------------------------------------------------------------
real, parameter :: gamma_c = 0.0144
real, parameter :: Re_visc_min = 0.111
real, parameter :: h_charnock = 10.0
real, parameter :: c1_charnock = log(h_charnock * (g / gamma_c))
real, parameter :: c2_charnock = Re_visc_min * nu_air * c1_charnock
! --------------------------------------------------------------------------------
contains
! charnock roughness definition
! --------------------------------------------------------------------------------
subroutine get_charnock_roughness(z0_m, u_dyn0, U, h, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_m !> aerodynamic roughness [m]
real, intent(out) :: u_dyn0 !> dynamic velocity in neutral conditions [m/s]
real, intent(in) :: h !> constant flux layer height [m]
real, intent(in) :: U !> abs(wind speed) [m/s]
integer, intent(in) :: maxiters !> maximum number of iterations
! ----------------------------------------------------------------------------
! --- local variables
real :: Uc ! wind speed at h_charnock [m/s]
real :: a, b, c, c_min
real :: f
integer :: i, j
! ----------------------------------------------------------------------------
Uc = U
a = 0.0
b = 25.0
c_min = log(h_charnock) / kappa
do i = 1, maxiters
f = c1_charnock - 2.0 * log(Uc)
do j = 1, maxiters
c = (f + 2.0 * log(b)) / kappa
! looks like the check should use U10 instead of U
! but note that a1 should be set = 0 in cycle beforehand
if (U <= 8.0e0) a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
c = max(c - a, c_min)
b = c
end do
z0_m = h_charnock * exp(-c * kappa)
z0_m = max(z0_m, 0.000015e0)
Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
end do
! --- define dynamic velocity in neutral conditions
u_dyn0 = Uc / c
end subroutine
! --------------------------------------------------------------------------------
end module sfx_roughness