module sfx_z0m_all_surface
!< @brief surface roughness parameterizations
use sfx_phys_const
implicit none
public :: get_dynamic_roughness_ch
public :: get_dynamic_roughness_map
public :: get_dynamic_roughness_ow
public :: get_dynamic_roughness_fetch
public :: get_dynamic_roughness_and
public :: get_dynamic_roughness_coast1
public :: get_dynamic_roughness_coast2
public :: get_dynamic_roughness_coast3
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
real, parameter, private :: kappa = 0.40 !< von Karman constant [n/d]
! --------------------------------------------------------------------------------
!< 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
real, parameter :: gamma_min = 0.01
real, parameter :: gamma_max = 0.11
real, parameter :: f_c = 100
real, parameter :: eps = 1
! --------------------------------------------------------------------------------
contains
! charnock roughness definition
! --------------------------------------------------------------------------------
subroutine get_dynamic_roughness_ch(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
!write(*,*) 'ch', z0_m, u_dyn0
end subroutine
! --------------------------------------------------------------------------------
subroutine get_dynamic_roughness_ow(z0_m, u_dyn0, U, h, maxiters)
!Owen 1964
! ----------------------------------------------------------------------------
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 :: b1, b2, Cuz, betta_u, nu_m, C_z0,c
real :: f
integer :: i, j
! ----------------------------------------------------------------------------
Uc=U
C_z0=0.007
betta_u=0.111
nu_m=0.0000133
b1=log(h*g/C_z0)
b2=betta_u*nu_m*g/C_z0
Cuz=25.0
do i = 1, maxiters
f = c1_charnock - 2.0 * log(Uc)
c = (f + 2.0 * log(Cuz)) / kappa
Cuz=(1.0/kappa)*(b1+log(U/Cuz)-log(b2+(U/Cuz)*(U/Cuz)))
if(Cuz==0.0) exit
z0_m=h*exp(-kappa*Cuz)
end do
u_dyn0 = Uc / c
!write(*,*) 'ow', z0_m, u_dyn0
end subroutine
! --------------------------------------------------------------------------------
subroutine get_dynamic_roughness_fetch(z0_m, u_dyn0, U, depth, 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) :: U !< abs(wind speed) [m/s]
real, intent(in) :: depth !< depth [m]
real, intent(in) :: h !< constant flux layer height [m]
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
real :: A_lake, B_lake, gamma_c, fetch, c1_charnock_lake, c2_charnock_lake
integer :: i, j
! ----------------------------------------------------------------------------
Uc = U
a = 0.0
b = 25.0
c_min = log(h_charnock) / kappa
fetch = 25.0 * depth !25.0 * depth
!< z0 = Re_visc_min * (nu / u_dyn) + gamma_c * (u_dyn^2 / g)
!< gamma_c = gamma_min + (gamma_max - gamma_min) * exp(-min(A_lake, B_lake))
!< А_lake = (fetch * g / U^2)^(1/3) / f_c
!< B_lake = eps (sqrt(depth * g)/U)
do i = 1, maxiters
A_lake = ((fetch * g / (U)**2)**(1/3)) / f_c
B_lake = eps * (sqrt(depth * g)/U)
gamma_c = gamma_min + (gamma_max - gamma_min) * exp(-min(A_lake, B_lake))
!write(*,*) A_lake
!write(*,*) B_lake
c1_charnock_lake = log(h_charnock * (g / gamma_c))
c2_charnock_lake = Re_visc_min * nu_air * c1_charnock_lake
f = c1_charnock_lake - 2.0 * log(Uc)
do j = 1, maxiters
c = (f + 2.0 * log(b)) / kappa
if (U <= 8.0e0) a = log(1.0 + c2_charnock_lake * ((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
subroutine get_dynamic_roughness_map(z0_m, u_dyn0, U, h, z0m_map)
! ----------------------------------------------------------------------------
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) :: z0m_map !< aerodynamic roughness from map[m]
real, intent(in) :: U !< abs(wind speed) [m/s]
! ----------------------------------------------------------------------------
real :: h0_m
z0_m=z0m_map
h0_m = h / z0_m
u_dyn0 = U * kappa / log(h0_m)
!write(*,*) 'map', z0_m, u_dyn0
end subroutine
! --------------------------------------------------------------------------------
subroutine get_dynamic_roughness_and(z0_m, u_dyn0, U, h, z0m_map)
! ----------------------------------------------------------------------------
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) :: z0m_map !< aerodynamic roughness from map[m]
real, intent(in) :: U !< abs(wind speed) [m/s]
! ----------------------------------------------------------------------------
real :: h0_m, u_star_prev, nu, g
real :: tolerance
integer :: max_iterations, iter
nu = 1.7e-5
g = 9.81
u_dyn0 = 0.2
tolerance = 1.0e-5
max_iterations = 10
do iter = 1, max_iterations
u_star_prev = u_dyn0
z0_m = (0.135 * nu / u_dyn0) + (0.035 * u_dyn0**2 / g) * &
(1.0 + exp(-((u_dyn0 - 0.18) / 0.1)**2))
h0_m = h / z0_m
u_dyn0 = U * kappa / log(h0_m)
if (abs(u_dyn0 - u_star_prev) < tolerance) exit
end do
u_dyn0 = U * kappa / log(h0_m)
end subroutine
! --------------------------------------------------------------------------------
subroutine get_dynamic_roughness_coast1(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
! ----------------------------------------------------------------------------
! Taylor&Yelland formulation
! --- local variables
real :: Uc, c, c_min, a, b ! wind speed at h_charnock [m/s] & Cd^(-1)
real :: pi=3.14159265
real :: hs, Tp, Lp
integer :: i, j
! ----------------------------------------------------------------------------
Uc = U
a = 0
b = 0
c_min = log(h_charnock) / kappa
do i = 1, maxiters
hs = 0.0248*(Uc**2.) !hs is the significant wave height
Tp = 0.729*max(Uc,0.1) !Tp dominant wave period
Lp = g*Tp**2/(2*pi) !Lp is the wavelength of the dominant wave
do j = 1, maxiters
z0_m = 1200.*hs*(hs/Lp)**4.5
c = log(h_charnock / z0_m) / kappa
if (Uc <= 8.0e0) then
a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
c = max(c - a, c_min)
b = c
z0_m = h_charnock * exp(-kappa*c)
end if
end do
z0_m = max( z0_m, 1.27e-7) !These max/mins were suggested by
z0_m = min( z0_m, 2.85e-3) !Davis et al. (2008)
Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
end do
u_dyn0 = Uc * kappa / log (h_charnock / z0_m)
! c = Uc/u_dyn0
! write (*,*) 'out1', u_dyn0
end subroutine
! Hwang
subroutine get_dynamic_roughness_coast2(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, c, c1, a, b, c_min ! wind speed at h_charnock [m/s] & Cd^(-1)
integer :: i, j
! ----------------------------------------------------------------------------
Uc = U
a = 0
b = 0
c_min = log(h_charnock) / kappa
do i = 1, maxiters
c1 = 0.016 * Uc**2
do j = 1, maxiters
c = 1.0 / (0.01*sqrt(8.058 + 0.967 * Uc - c1))
if (Uc <= 8.0e0) then
a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
c = max(c - a, c_min)
b = c
end if
end do
z0_m = h_charnock * exp(-kappa*c)
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
! write (*,*) 'out1', u_dyn0
! --------------------------------------------------------------------------------
end subroutine
! Zhao et al 2015 10^3*z0=15.6*u_dyn0**2/g + 10**(-2)
subroutine get_dynamic_roughness_coast3(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, u_dyn, z0_m0 ! wind speed at h_charnock [m/s]
real :: c, c_min, a, b
real :: f1
integer :: i, j
! ----------------------------------------------------------------------------
Uc = U
u_dyn = U / 28.0
a = 0
b = 0
z0_m0=0.082
c_min = log(h_charnock) / kappa
if (Uc >8.00 .AND. Uc <= 35.0e0) then
do i = 1, maxiters
f1 = 15.6*u_dyn**2/g
z0_m0 = 0.001 * (f1 + 0.01)
c = log(h_charnock / z0_m0) / kappa
end do
end if
if (Uc <= 8.0e0) then
do j = 1, maxiters
a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
c = max(c - a, c_min)
b = c
z0_m0 = h_charnock * exp(-c * kappa)
end do
else
z0_m0 = 2.82
end if
z0_m = max(z0_m0, 0.000015e0)
Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
! --- define dynamic velocity in neutral conditions
! u_dyn0 = Uc / c
u_dyn0 = Uc * kappa / log (h_charnock / z0_m)
end subroutine
end module sfx_z0m_all_surface