diff --git a/srcF/sfx_z0m_all.f90 b/srcF/sfx_z0m_all.f90 new file mode 100644 index 0000000000000000000000000000000000000000..89fcaf3df9d8a4e24fc115342e1821309fd9c2b7 --- /dev/null +++ b/srcF/sfx_z0m_all.f90 @@ -0,0 +1,208 @@ +module sfx_z0m_all + + !< @brief surface thermal roughness parameterizations for all type surface + + use z0m_all_surface + + + implicit none + public :: get_dynamic_roughness_all + public :: get_dynamic_roughness_definition + + + integer, public, parameter :: surface_ocean = 0 !< ocean surface + integer, public, parameter :: surface_land = 1 !< land surface + integer, public, parameter :: surface_lake = 2 !< lake surface + integer, public, parameter :: surface_snow = 3 !< snow covered surface + integer, public, parameter :: surface_forest = 4 !< forest csurface + integer, public, parameter :: surface_user = 5 !< user surface + + + integer, public, parameter :: z0m_ch = 0 + integer, public, parameter :: z0m_fe = 1 + integer, public, parameter :: z0m_ow = 2 + integer, public, parameter :: z0m_map = 3 + integer, public, parameter :: z0m_user = 4 + + + + character(len = 16), parameter :: surface_ocean_tag = 'ocean' + character(len = 16), parameter :: surface_land_tag = 'land' + character(len = 16), parameter :: surface_lake_tag = 'lake' + character(len = 16), parameter :: surface_snow_tag = 'snow' + character(len = 16), parameter :: surface_forest_tag = 'forest' + character(len = 16), parameter :: surface_user_tag = 'user' + + character(len = 16), parameter :: z0m_tag_ch = 'charnock' + character(len = 16), parameter :: z0m_tag_fe = 'fetch' + character(len = 16), parameter :: z0m_tag_ow = 'owen' + character(len = 16), parameter :: z0m_tag_map = 'map' + character(len = 16), parameter :: z0m_tag_user = 'user' + + + integer, public, parameter :: ocean_z0m_id = z0m_ch !< ocean surface + integer, public, parameter :: land_z0m_id = z0m_map !< land surface + integer, public, parameter :: lake_z0m_id = z0m_fe !< lake surface + integer, public, parameter :: snow_z0m_id = z0m_ow !< snow covered surface + integer, public, parameter :: forest_z0m_id = z0m_map !< forest csurface + integer, public, parameter :: usersf_z0m_id = z0m_map !< user surface + + + contains + + + ! surface type definition + ! -------------------------------------------------------------------------------- + function get_surface_id(tag) result(id) + implicit none + character(len=*), intent(in) :: tag + integer :: id + + id = - 1 + if (trim(tag) == trim(surface_ocean_tag)) then + id = surface_ocean + else if (trim(tag) == trim(surface_land_tag)) then + id = surface_land + else if (trim(tag) == trim(surface_lake_tag)) then + id = surface_lake + else if (trim(tag) == trim(surface_snow_tag)) then + id = surface_snow + end if + + end function + + function get_surface_tag(id) result(tag) + implicit none + integer :: id + character(len=:), allocatable :: tag + + tag = 'undefined' + if (id == surface_ocean) then + tag = surface_ocean_tag + else if (id == surface_land) then + tag = surface_land_tag + else if (id == surface_lake) then + tag = surface_lake_tag + else if (id == surface_snow) then + tag = surface_snow_tag + end if + + end function + + + ! -------------------------------------------------------------------------------- + ! surface type definition + ! -------------------------------------------------------------------------------- + function get_surface_z0m_id(tag_z0m) result(z0m_id) + implicit none + character(len=*), intent(in) :: tag_z0m + integer :: z0m_id + + z0m_id = - 1 + if (trim(tag_z0m) == trim(z0m_tag_ch)) then + z0m_id = z0m_ch + else if (trim(tag_z0m) == trim(z0m_tag_fe)) then + z0m_id = z0m_fe + else if (trim(tag_z0m) == trim(z0m_tag_ow)) then + z0m_id = z0m_ow + else if (trim(tag_z0m) == trim(z0m_tag_map)) then + z0m_id = z0m_map + else if (trim(tag_z0m) == trim(z0m_tag_user)) then + z0m_id = z0m_user + end if + + end function + + function get_surface_z0m_tag(z0m_id) result(tag_z0m) + implicit none + integer :: z0m_id + character(len=:), allocatable :: tag_z0m + + tag_z0m = 'undefined' + if (z0m_id == z0m_ch) then + tag_z0m = z0m_tag_ch + else if (z0m_id == z0m_fe) then + tag_z0m = z0m_tag_fe + else if (z0m_id == z0m_ow) then + tag_z0m = z0m_tag_ow + else if (z0m_id == z0m_map) then + tag_z0m = z0m_tag_map + else if (z0m_id == z0m_user) then + tag_z0m = z0m_tag_user + end if + end function + + + + ! ---------------------------------------------------------------------------- + subroutine get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h,& + maxiters, z0m_map, z0m_id) + ! ---------------------------------------------------------------------------- + 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] + real, intent(in) :: z0m_map !< aerodynamic roughness from map [m] + integer, intent(in) :: maxiters !< maximum number of iterations + + integer, intent(in) :: z0m_id + ! --------------------------------------------------------------------------- + + + if (z0m_id == z0m_ch) then + call get_dynamic_roughness_ch(z0_m, u_dyn0, U, h, maxiters) + else if (z0m_id == z0m_fe) then + call get_dynamic_roughness_fetch(z0_m, u_dyn0, U, depth, h, maxiters) + else if (z0m_id == z0m_ow) then + call get_dynamic_roughness_ow(z0_m, u_dyn0, U, h, maxiters) + else if (z0m_id == z0m_map) then + call get_dynamic_roughness_map(z0_m, u_dyn0, U, h, z0m_map) + else if (z0m_id == z0m_user) then + write(*, *) 'z0m_user' + end if + + end subroutine + ! -------------------------------------------------------------------------------- + + ! ---------------------------------------------------------------------------- + subroutine get_dynamic_roughness_definition(surface_type, id_ocean, id_land, id_snow, id_lake, & + id_forest, id_user, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, & + forest_z0m_id, usersf_z0m_id) + ! ---------------------------------------------------------------------------- + real, intent(out) :: ocean_z0m_id + real, intent(out) :: land_z0m_id + real, intent(out) :: lake_z0m_id + real, intent(out) :: snow_z0m_id + real, intent(out) :: forest_z0m_id + real, intent(out) :: usersf_z0m_id + + real, intent(in) :: surface_type + + real, intent(in) :: id_ocean + real, intent(in) :: id_land + real, intent(in) :: id_snow + real, intent(in) :: id_lake + real, intent(in) :: id_forest + real, intent(in) :: id_user + + ! --------------------------------------------------------------------------- + + if (surface_type == surface_ocean) then + ocean_z0m_id = id_ocean + else if (surface_type == surface_land) then + land_z0m_id = id_land + else if (surface_type == surface_snow) then + snow_z0m_id = id_snow + else if (surface_type == surface_lake) then + lake_z0m_id = id_lake + else if (surface_type == surface_forest) then + forest_z0m_id = id_forest + else if (surface_type == surface_user) then + usersf_z0m_id = id_user + end if + + end subroutine + +end module sfx_z0m_all diff --git a/srcF/sfx_z0t_all.f90 b/srcF/sfx_z0t_all.f90 new file mode 100644 index 0000000000000000000000000000000000000000..8f615b145bb998a1455afbc62914041dea867202 --- /dev/null +++ b/srcF/sfx_z0t_all.f90 @@ -0,0 +1,264 @@ +module sfx_z0t_all + + !< @brief surface thermal roughness parameterizations for all type surface + + use z0t_all_surface + + + implicit none + public :: get_thermal_roughness_all + public :: get_thermal_roughness_definition + + + integer, public, parameter :: surface_ocean = 0 !< ocean surface + integer, public, parameter :: surface_land = 1 !< land surface + integer, public, parameter :: surface_lake = 2 !< lake surface + integer, public, parameter :: surface_snow = 3 !< snow covered surface + integer, public, parameter :: surface_forest = 4 !< forest surface + integer, public, parameter :: surface_user = 5 !< user surface + + + integer, public, parameter :: z0t_kl_water = 0 + integer, public, parameter :: z0t_kl_land = 1 + integer, public, parameter :: z0t_re = 2 + integer, public, parameter :: z0t_zi = 3 + integer, public, parameter :: z0t_ca = 4 + integer, public, parameter :: z0t_cz = 5 + integer, public, parameter :: z0t_br = 6 + integer, public, parameter :: z0t_ot = 7 + integer, public, parameter :: z0t_du = 8 + integer, public, parameter :: z0t_zm = 9 + integer, public, parameter :: z0t_mix = 10 + integer, public, parameter :: z0t_user = 11 + + + + character(len = 16), parameter :: surface_ocean_tag = 'ocean' + character(len = 16), parameter :: surface_land_tag = 'land' + character(len = 16), parameter :: surface_lake_tag = 'lake' + character(len = 16), parameter :: surface_snow_tag = 'snow' + character(len = 16), parameter :: surface_forest_tag = 'forest' + character(len = 16), parameter :: surface_user_tag = 'user' + + character(len = 16), parameter :: z0t_tag_kl_water = 'kl_water' + character(len = 16), parameter :: z0t_tag_kl_land = 'kl_land' + character(len = 16), parameter :: z0t_tag_re = 're' + character(len = 16), parameter :: z0t_tag_zi = 'zi' + character(len = 16), parameter :: z0t_tag_ca = 'ca' + character(len = 16), parameter :: z0t_tag_cz = 'cz' + character(len = 16), parameter :: z0t_tag_br = 'br' + character(len = 16), parameter :: z0t_tag_ot = 'ot' + character(len = 16), parameter :: z0t_tag_du = 'du' + character(len = 16), parameter :: z0t_tag_zm = 'zm' + character(len = 16), parameter :: z0t_tag_mix = 'mix' + character(len = 16), parameter :: z0t_tag_user = 'zt_user' + + integer, public, parameter :: ocean_z0t_id = z0t_kl_water !< ocean surface + integer, public, parameter :: land_z0t_id = z0t_kl_land !< land surface + integer, public, parameter :: lake_z0t_id = z0t_re !< lake surface + integer, public, parameter :: snow_z0t_id = z0t_ca !< snow covered surface + integer, public, parameter :: forest_z0t_id = z0t_du !< forest csurface + integer, public, parameter :: usersf_z0t_id = z0t_mix !< user surface + + + contains + + + ! surface type definition + ! -------------------------------------------------------------------------------- + function get_surface_id(tag) result(id) + implicit none + character(len=*), intent(in) :: tag + integer :: id + + id = - 1 + if (trim(tag) == trim(surface_ocean_tag)) then + id = surface_ocean + else if (trim(tag) == trim(surface_land_tag)) then + id = surface_land + else if (trim(tag) == trim(surface_lake_tag)) then + id = surface_lake + else if (trim(tag) == trim(surface_snow_tag)) then + id = surface_snow + end if + + end function + + function get_surface_tag(id) result(tag) + implicit none + integer :: id + character(len=:), allocatable :: tag + + tag = 'undefined' + if (id == surface_ocean) then + tag = surface_ocean_tag + else if (id == surface_land) then + tag = surface_land_tag + else if (id == surface_lake) then + tag = surface_lake_tag + else if (id == surface_snow) then + tag = surface_snow_tag + end if + + end function + + + ! -------------------------------------------------------------------------------- + ! surface type definition + ! -------------------------------------------------------------------------------- + function get_surface_z0t_id(tag_z0t) result(z0t_id) + implicit none + character(len=*), intent(in) :: tag_z0t + integer :: z0t_id + + z0t_id = - 1 + if (trim(tag_z0t) == trim(z0t_tag_kl_water)) then + z0t_id = z0t_kl_water + else if (trim(tag_z0t) == trim(z0t_tag_kl_land)) then + z0t_id = z0t_kl_land + else if (trim(tag_z0t) == trim(z0t_tag_re)) then + z0t_id = z0t_re + else if (trim(tag_z0t) == trim(z0t_tag_zi)) then + z0t_id = z0t_zi + else if (trim(tag_z0t) == trim(z0t_tag_ca)) then + z0t_id = z0t_ca + else if (trim(tag_z0t) == trim(z0t_tag_cz)) then + z0t_id = z0t_cz + else if (trim(tag_z0t) == trim(z0t_tag_br)) then + z0t_id = z0t_br + else if (trim(tag_z0t) == trim(z0t_tag_ot)) then + z0t_id = z0t_ot + else if (trim(tag_z0t) == trim(z0t_tag_du)) then + z0t_id = z0t_du + else if (trim(tag_z0t) == trim(z0t_tag_zm)) then + z0t_id = z0t_zm + else if (trim(tag_z0t) == trim(z0t_tag_mix)) then + z0t_id = z0t_mix + else if (trim(tag_z0t) == trim(z0t_tag_user)) then + z0t_id = z0t_user + end if + + end function + + function get_surface_z0t_tag(z0t_id) result(tag_z0t) + implicit none + integer :: z0t_id + character(len=:), allocatable :: tag_z0t + + tag_z0t = 'undefined' + if (z0t_id == z0t_kl_water) then + tag_z0t = z0t_tag_kl_water + else if (z0t_id == z0t_kl_land) then + tag_z0t = z0t_tag_kl_land + else if (z0t_id == z0t_re) then + tag_z0t = z0t_tag_re + else if (z0t_id == z0t_zi) then + tag_z0t = z0t_tag_zi + else if (z0t_id == z0t_ca) then + tag_z0t = z0t_tag_ca + else if (z0t_id == z0t_cz) then + tag_z0t = z0t_tag_cz + else if (z0t_id == z0t_br) then + tag_z0t = z0t_tag_br + else if (z0t_id == z0t_ot) then + tag_z0t = z0t_tag_ot + else if (z0t_id == z0t_du) then + tag_z0t = z0t_tag_du + else if (z0t_id == z0t_zm) then + tag_z0t = z0t_tag_zm + else if (z0t_id == z0t_mix) then + tag_z0t = z0t_tag_mix + else if (z0t_id == z0t_user) then + tag_z0t = z0t_tag_user + end if + end function + + + + ! ---------------------------------------------------------------------------- + subroutine get_thermal_roughness_all(z0_t, B, & + z0_m, Re, u_dyn, Czm, LAI, z0t_id) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: Re !< roughness Reynolds number [n/d] + real, intent(in) :: LAI !< leaf-area index + real, intent(in) :: u_dyn !< dynamic velocity [m/s] + real, intent(in) :: Czm !< proportionality coefficient z0_t =Czm*z0_m + + integer, intent(in) :: z0t_id + ! --------------------------------------------------------------------------- + + + if (z0t_id == z0t_kl_water) then + call get_thermal_roughness_kl_water(z0_t, B, z0_m, Re) + else if (z0t_id == z0t_kl_land) then + call get_thermal_roughness_kl_land(z0_t, B, z0_m, Re) + else if (z0t_id == z0t_re) then + call get_thermal_roughness_re(z0_t, B, z0_m, Re) + else if (z0t_id == z0t_zi) then + call get_thermal_roughness_zi(z0_t, B, z0_m, Re) + else if (z0t_id == z0t_ca) then + call get_thermal_roughness_ca(z0_t, B, z0_m, Re) + else if (z0t_id == z0t_cz) then + call get_thermal_roughness_cz(z0_t, B, z0_m, Re) + else if (z0t_id == z0t_br) then + call get_thermal_roughness_br(z0_t, B, z0_m, Re) + else if (z0t_id == z0t_ot) then + call get_thermal_roughness_ot(z0_t, B, z0_m, Re) + else if (z0t_id == z0t_du) then + call get_thermal_roughness_du(z0_t, B, z0_m, u_dyn, LAI) + else if (z0t_id == z0t_zm) then + call get_thermal_roughness_zm(z0_t, B, z0_m, Czm) + else if (z0t_id == z0t_mix) then + call get_thermal_roughness_mix(z0_t, B, z0_m, u_dyn, Re) + else if (z0t_id == z0t_user) then + write(*, *) 'z0t_user' + end if + + + end subroutine + ! -------------------------------------------------------------------------------- + + ! ---------------------------------------------------------------------------- + subroutine get_thermal_roughness_definition(surface_type, id_ocean, id_land, id_snow, id_lake, & + id_forest, id_user, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, & + forest_z0t_id, usersf_z0t_id) + ! ---------------------------------------------------------------------------- + real, intent(out) :: ocean_z0t_id + real, intent(out) :: land_z0t_id + real, intent(out) :: lake_z0t_id + real, intent(out) :: snow_z0t_id + real, intent(out) :: forest_z0t_id + real, intent(out) :: usersf_z0t_id + + real, intent(in) :: surface_type + + real, intent(in) :: id_ocean + real, intent(in) :: id_land + real, intent(in) :: id_snow + real, intent(in) :: id_lake + real, intent(in) :: id_forest + real, intent(in) :: id_user + + ! --------------------------------------------------------------------------- + + if (surface_type == surface_ocean) then + ocean_z0t_id = id_ocean + else if (surface_type == surface_land) then + land_z0t_id = id_land + else if (surface_type == surface_snow) then + snow_z0t_id = id_snow + else if (surface_type == surface_lake) then + lake_z0t_id = id_lake + else if (surface_type == surface_forest) then + forest_z0t_id = id_forest + else if (surface_type == surface_user) then + usersf_z0t_id = id_user + end if + + end subroutine + +end module sfx_z0t_all diff --git a/srcF/z0m_all_surface.f90 b/srcF/z0m_all_surface.f90 new file mode 100644 index 0000000000000000000000000000000000000000..46262d6e5798afda7f596b5607a1ab852cbc2f70 --- /dev/null +++ b/srcF/z0m_all_surface.f90 @@ -0,0 +1,207 @@ +module z0m_all_surface + !< @brief surface roughness parameterizations + + use sfx_phys_const + implicit none + public + ! -------------------------------------------------------------------------------- + + + + + ! -------------------------------------------------------------------------------- + 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 + + 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 + + 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) + + end subroutine + ! -------------------------------------------------------------------------------- + + + + + + end module z0m_all_surface \ No newline at end of file diff --git a/srcF/z0t_all_surface.f90 b/srcF/z0t_all_surface.f90 new file mode 100644 index 0000000000000000000000000000000000000000..ca3a8b4f8b1d9baa4a84c1c7698b2313e82a9e78 --- /dev/null +++ b/srcF/z0t_all_surface.f90 @@ -0,0 +1,285 @@ +module z0t_all_surface + !< @brief surface thermal roughness parameterizations + + implicit none + public + + + + ! -------------------------------------------------------------------------------- + real, parameter, private :: kappa = 0.40 !< von Karman constant [n/d] + real, parameter, private :: Pr_m = 0.71 !< molecular Prandtl number (air) [n/d] + !< Re fully roughness minimum value [n/d] + real, parameter :: Re_rough_min = 16.3 + !< roughness model coeff. [n/d] + !< --- transitional mode + !< B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2 + real, parameter :: B1_rough = 5.0 / 6.0 + real, parameter :: B2_rough = 0.45 + real, parameter :: B3_rough = kappa * Pr_m + !< --- fully rough mode (Re > Re_rough_min) + !< B = B4 * Re^(B2) + real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8) + real, parameter :: B_max_ocean = 8.0 + real, parameter :: B_max_land = 2.0 + + + contains + + + ! thermal roughness definition by Kazakov, Lykosov + ! -------------------------------------------------------------------------------- + subroutine get_thermal_roughness_kl_land(z0_t, B, & + z0_m, Re) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: Re !< roughness Reynolds number [n/d] + ! ---------------------------------------------------------------------------- + + !--- define B = log(z0_m / z0_t) + if (Re <= Re_rough_min) then + B = B1_rough * alog(B3_rough * Re) + B2_rough + else + ! *: B4 takes into account Re value at z' ~ O(10) z0 + B = B4_rough * (Re**B2_rough) + end if + + B = min(B, B_max_land) + + z0_t = z0_m / exp(B) + + end subroutine + ! -------------------------------------------------------------------------------- + +! -------------------------------------------------------------------------------- + subroutine get_thermal_roughness_kl_water(z0_t, B, & + z0_m, Re) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: Re !< roughness Reynolds number [n/d] + ! ---------------------------------------------------------------------------- + + !--- define B = log(z0_m / z0_t) + if (Re <= Re_rough_min) then + B = B1_rough * alog(B3_rough * Re) + B2_rough + else + ! *: B4 takes into account Re value at z' ~ O(10) z0 + B = B4_rough * (Re**B2_rough) + end if + + B = min(B, B_max_ocean) + + z0_t = z0_m / exp(B) + + end subroutine + + + ! thermal roughness definition by Chen, F., Zhang, Y., 2009. + ! -------------------------------------------------------------------------------- + subroutine get_thermal_roughness_cz(z0_t, B, & + z0_m, Re) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: Re !< roughness Reynolds number [n/d] + + + + B=(kappa*10.0**(-0.4*z0_m/0.07))*(Re**0.45) !Chen and Zhang + + ! --- define roughness [thermal] + z0_t = z0_m / exp(B) + + end subroutine + ! -------------------------------------------------------------------------------- + ! thermal roughness definition by Zilitinkevich, S., 1995. + ! -------------------------------------------------------------------------------- + subroutine get_thermal_roughness_zi(z0_t, B, & + z0_m, Re) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: Re !< roughness Reynolds number [n/d] + + + + + B=0.1*kappa*(Re**0.5) !6-Zilitinkevich + + + + ! --- define roughness [thermal] + z0_t = z0_m / exp(B) + + end subroutine + ! -------------------------------------------------------------------------------- + ! thermal roughness definition by Cahill, A.T., Parlange, M.B., Albertson, J.D., 1997. + ! It is better to use for dynamic surfaces such as sand + ! -------------------------------------------------------------------------------- + subroutine get_thermal_roughness_ca(z0_t, B, & + z0_m, Re) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: Re !< roughness Reynolds number [n/d] + + + + + B=2.46*(Re**0.25)-3.8 !4-Cahill et al. + + ! --- define roughness [thermal] + z0_t = z0_m / exp(B) + + end subroutine + + + ! -------------------------------------------------------------------------------- + ! thermal roughness definition by Brutsaert W., 2003. + ! -------------------------------------------------------------------------------- + subroutine get_thermal_roughness_br(z0_t, B, & + z0_m, Re) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: Re !< roughness Reynolds number [n/d] + + + + + B=2.46*(Re**0.25)-2.0 !Brutsaert + + ! --- define roughness [thermal] + z0_t = z0_m / exp(B) + + end subroutine + + + + ! -------------------------------------------------------------------------------- + ! thermal roughness definition by Owen P. R., Thomson W. R., 1963. + ! -------------------------------------------------------------------------------- + subroutine get_thermal_roughness_ot(z0_t, B, & + z0_m, Re) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: Re !< roughness Reynolds number [n/d] + + + + B=kappa*(Re**0.45) !Owen P. R., Thomson W. R. + + ! --- define roughness [thermal] + z0_t = z0_m / exp(B) + + end subroutine + ! -------------------------------------------------------------------------------- + ! thermal roughness definition by Duynkerke P. G., 1992. + !It is better to use for surfaces wiht forest + ! -------------------------------------------------------------------------------- + subroutine get_thermal_roughness_du(z0_t, B, & + z0_m, u_dyn, LAI) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: u_dyn !< dynamic velocity [m/s] + real, intent(in) :: LAI !< leaf-area index + + + + + B=(13*u_dyn**0.4)/LAI+0.85 !Duynkerke P. G., 1992. + + ! --- define roughness [thermal] + z0_t = z0_m / exp(B) + + end subroutine + ! -------------------------------------------------------------------------------- + ! thermal roughness definition z0_t = C*z0_m + ! -------------------------------------------------------------------------------- + subroutine get_thermal_roughness_zm(z0_t, B, & + z0_m, Czm) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: Czm !< proportionality coefficient + + + + z0_t =Czm*z0_m + B=log(z0_m / z0_t) + end subroutine + ! -------------------------------------------------------------------------------- + ! thermal roughness definition by Chen and Zhang and Zilitinkevich + ! -------------------------------------------------------------------------------- + subroutine get_thermal_roughness_mix(z0_t, B, & + z0_m, u_dyn, Re) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: u_dyn !< dynamic velocity [m/s] + real, intent(in) :: Re !< roughness Reynolds number [n/d] + + + real, parameter :: u_dyn_th=0.17 !< dynamic velocity treshhold [m/s] + + if (u_dyn <= u_dyn_th) then + B=0.1*kappa*(Re**0.5) !Zilitinkevich + else + B=(kappa*10.0**(-0.4*z0_m/0.07))*(Re**0.45) !Chen and Zhang + end if + + + ! --- define roughness [thermal] + z0_t = z0_m / exp(B) + + end subroutine + +! -------------------------------------------------------------------------------- + subroutine get_thermal_roughness_re(z0_t, B, & + z0_m, Re) + ! ---------------------------------------------------------------------------- + real, intent(out) :: z0_t !< thermal roughness [m] + real, intent(out) :: B !< = log(z0_m / z0_t) [n/d] + + real, intent(in) :: z0_m !< aerodynamic roughness [m] + real, intent(in) :: Re !< roughness Reynolds number [n/d] + + + + + + B=alog(-0.56*(4.0*(Re)**(0.5)-3.4)) !Repina, 2023 + + + ! --- define roughness [thermal] + z0_t = z0_m / exp(B) + + end subroutine + + + +end module z0t_all_surface