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