diff --git a/srcF/sfx_config.f90 b/srcF/sfx_config.f90
index b2a26202b9e09e683f3d73ce471d8b5bb80cbcc7..68df2c56de33194ccfce95cc22bedf691b7eff7b 100644
--- a/srcF/sfx_config.f90
+++ b/srcF/sfx_config.f90
@@ -167,8 +167,8 @@ contains
dataset%filename = get_dataset_filename(id)
dataset%nmax = 0
- dataset%surface = surface_ice
- dataset%surface_type = surface_ice
+ dataset%surface = surface_lake
+ dataset%surface_type = surface_lake
dataset%z0_h = -1.0
dataset%lai = 1.0
dataset%depth = 10.0
@@ -194,7 +194,7 @@ contains
else if (id == dataset_toga) then
! *: check & fix values
dataset%h = 15.0
- dataset%surface = surface_ocean
+ dataset%surface = surface_lake
dataset%z0_m = -1.0
end if
write (*,*) dataset%surface, 'dataset%surface'
diff --git a/srcF/sfx_surface.f90 b/srcF/sfx_surface.f90
index 8f4ebe64f9a9b336499f3d419dfce0640efb98c5..ed09eb1e7d712d4ef051bc2cb99e94845f455570 100644
--- a/srcF/sfx_surface.f90
+++ b/srcF/sfx_surface.f90
@@ -53,7 +53,9 @@ module sfx_surface
integer, public, parameter :: z0m_map_id = 3
integer, public, parameter :: z0m_user = 4
integer, public, parameter :: z0m_and = 5
-
+ integer, public, parameter :: z0m_cst1 = 6
+ integer, public, parameter :: z0m_cst2 = 7
+ integer, public, parameter :: z0m_cst3 = 8
integer, public, parameter :: z0t_kl_water = 0
integer, public, parameter :: z0t_kl_land = 1
@@ -75,6 +77,9 @@ module sfx_surface
character(len = 16), parameter :: z0m_tag_map = 'map'
character(len = 16), parameter :: z0m_tag_user = 'user'
character(len = 16), parameter :: z0m_tag_and = 'andreas'
+ character(len = 16), parameter :: z0m_tag_cst1 = 'coast1'
+ character(len = 16), parameter :: z0m_tag_cst2 = 'coast2'
+ character(len = 16), parameter :: z0m_tag_cst3 = 'coast3'
character(len = 16), parameter :: z0t_tag_kl_water = 'kl_water'
@@ -93,7 +98,7 @@ module sfx_surface
integer, public, parameter :: ocean_z0m_id = z0m_ch !< ocean surface
integer, public, parameter :: land_z0m_id = z0m_map_id !< land surface
- integer, public, parameter :: lake_z0m_id = z0m_fe !< lake surface
+ integer, public, parameter :: lake_z0m_id = z0m_map_id !< lake surface
integer, public, parameter :: snow_z0m_id = z0m_map_id !< snow covered surface
integer, public, parameter :: forest_z0m_id = z0m_map_id !< forest csurface
integer, public, parameter :: usersf_z0m_id = z0m_ch !< user surface
@@ -101,12 +106,12 @@ module sfx_surface
integer, public, parameter :: ocean_z0t_id = z0t_kl_water !< ocean surface
- integer, public, parameter :: land_z0t_id = z0t_mix !< land surface
- integer, public, parameter :: lake_z0t_id = z0t_cz !< lake surface
- integer, public, parameter :: snow_z0t_id = z0t_zi !< snow covered surface
- integer, public, parameter :: forest_z0t_id = z0t_ot !< forest csurface
+ integer, public, parameter :: land_z0t_id = z0t_kl_land !< land surface
+ integer, public, parameter :: lake_z0t_id = z0t_kl_land !< lake surface
+ integer, public, parameter :: snow_z0t_id = z0t_kl_land !< snow covered surface
+ integer, public, parameter :: forest_z0t_id = z0t_kl_land !< forest csurface
integer, public, parameter :: usersf_z0t_id = z0t_kl_water !< user surface
- integer, public, parameter :: ice_z0t_id = z0t_zi !< user surface
+ integer, public, parameter :: ice_z0t_id = z0t_kl_land !< user surface
! --------------------------------------------------------------------------------
real, parameter, private :: kappa = 0.40 !< von Karman constant [n/d]
@@ -210,6 +215,12 @@ contains
z0m_id = z0m_map_id
else if (trim(tag_z0m) == trim(z0m_tag_user)) then
z0m_id = z0m_user
+ else if (trim(tag_z0m) == trim(z0m_tag_cst1)) then
+ z0m_id = z0m_cst1
+ else if (trim(tag_z0m) == trim(z0m_tag_cst2)) then
+ z0m_id = z0m_cst2
+ else if (trim(tag_z0m) == trim(z0m_tag_cst3)) then
+ z0m_id = z0m_cst3
end if
end function
@@ -230,6 +241,12 @@ contains
tag_z0m = z0m_tag_map
else if (z0m_id == z0m_user) then
tag_z0m = z0m_tag_user
+ else if (z0m_id == z0m_cst1) then
+ tag_z0m = z0m_tag_cst1
+ else if (z0m_id == z0m_cst2) then
+ tag_z0m = z0m_tag_cst2
+ else if (z0m_id == z0m_cst3) then
+ tag_z0m = z0m_tag_cst3
end if
end function
@@ -393,6 +410,12 @@ contains
call get_dynamic_roughness_and(z0_m, u_dyn0, U, h, z0m_map)
else if (z0m_id == z0m_user) then
write(*, *) 'z0m_user'
+ else if (z0m_id == z0m_cst1) then
+ call get_dynamic_roughness_coast1(z0_m, u_dyn0, U, h, maxiters)
+ else if (z0m_id == z0m_cst2) then
+ call get_dynamic_roughness_coast2(z0_m, u_dyn0, U, h, maxiters)
+ else if (z0m_id == z0m_cst3) then
+ call get_dynamic_roughness_coast3(z0_m, u_dyn0, U, h, maxiters)
end if
end subroutine
diff --git a/srcF/sfx_z0m_all_surface.f90 b/srcF/sfx_z0m_all_surface.f90
index 4d6f2edea2d7cd0a67b9a5db64ec2c9ceacd5646..1f8b69c635eb91278ba7b1ad4e76733722d2cd67 100644
--- a/srcF/sfx_z0m_all_surface.f90
+++ b/srcF/sfx_z0m_all_surface.f90
@@ -8,6 +8,9 @@ module sfx_z0m_all_surface
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
! --------------------------------------------------------------------------------
@@ -252,6 +255,139 @@ subroutine get_dynamic_roughness_map(z0_m, u_dyn0, U, h, z0m_map)
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
\ No newline at end of file