Skip to content
Snippets Groups Projects
Commit d2bc97ad authored by Виктория Суязова's avatar Виктория Суязова Committed by Anna Shestakova
Browse files

added coast1,coast2

parent 53ad7eb5
No related branches found
No related tags found
No related merge requests found
......@@ -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'
......
......@@ -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
......
......@@ -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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment