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

esm, most, most_snow bug fixed

parent bc6b2d07
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 35376 additions and 33401 deletions
CMakeLists.txt
+ 0
− 2
View file @ 54ab84cb
......@@ -85,8 +85,6 @@ set(SOURCES_F
srcF/sfx_phys_const.f90
srcF/sfx_z0m_all_surface.f90
srcF/sfx_z0t_all_surface.f90
srcF/sfx_z0m_all.f90
srcF/sfx_z0t_all.f90
srcF/sfx_surface.f90
srcF/sfx_most.f90
srcF/sfx_most_param.f90
......
data/MOSAiC.txt
+ 35199
− 33327
View file @ 54ab84cb
This diff is collapsed.
makefile
+ 1
− 1
View file @ 54ab84cb
......@@ -10,7 +10,7 @@ ifeq ($(COMPILER),gnu)
FC = gfortran
endif
OBJ_F90 = sfx_phys_const.o sfx_common.o sfx_io.o sfx_data.o sfx_z0t_all_surface.o sfx_z0m_all_surface.o sfx_surface.o sfx_config.o sfx_log_param.o sfx_log.o sfx_most_param.o sfx_most.o sfx_most_snow_param.o sfx_most_snow.o sfx_sheba_param.o sfx_sheba.o sfx_esm_param.o sfx_esm.o sfx_sheba_noit_param.o sfx_sheba_noniterative.o sfx_run.o sfx_main.o
OBJ_F90 = sfx_phys_const.o sfx_common.o sfx_io.o sfx_data.o sfx_z0t_all_surface.o sfx_z0m_all_surface.o sfx_z0m_all.o sfx_z0t_all.o sfx_surface.o sfx_config.o sfx_log_param.o sfx_log.o sfx_most_param.o sfx_most.o sfx_most_snow_param.o sfx_most_snow.o sfx_sheba_param.o sfx_sheba.o sfx_esm_param.o sfx_esm.o sfx_sheba_noit_param.o sfx_sheba_noniterative.o sfx_run.o sfx_main.o
OBJ_F =
OBJ = $(OBJ_F90) $(OBJ_F)
......
srcF/sfx_config.f90
+ 4
− 4
View file @ 54ab84cb
......@@ -174,8 +174,8 @@ contains
dataset%depth = 10.0
if (id == dataset_mosaic) then
dataset%h = 3.8
dataset%z0_m = 0.01
dataset%h = 6.0
dataset%z0_m = 0.0078
else if (id == dataset_irgason) then
dataset%h = 9.2
dataset%z0_m = 0.01
......@@ -189,7 +189,7 @@ contains
dataset%z0_m = -1.0
else if (id == dataset_papa) then
dataset%h = 10.0
dataset%surface = surface_ocean
dataset%surface = surface_lake
dataset%z0_m = -1.0
else if (id == dataset_toga) then
! *: check & fix values
......@@ -197,7 +197,7 @@ contains
dataset%surface = surface_ocean
dataset%z0_m = -1.0
end if
write (*,*) dataset%surface, 'dataset%surface'
end subroutine
function get_dataset_filename(id) result(filename)
......
srcF/sfx_esm.f90
+ 10
− 5
View file @ 54ab84cb
......@@ -209,7 +209,7 @@ contains
integer surface_type !< surface type = (ocean || land)
real fval !< just a shortcut for partial calculations
real z0_m1
#ifdef SFX_CHECK_NAN
real NaN
#endif
......@@ -235,22 +235,27 @@ contains
dT = meteo%dT
dQ = meteo%dQ
h = meteo%h
z0_m = meteo%z0_m
z0_m1 = meteo%z0_m
depth = meteo%depth
lai = meteo%lai
surface_type=meteo%surface_type
!write (*,*) surface_type, 'esm'
call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
forest_z0m_id, usersf_z0m_id, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m1, z0m_id)
call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
forest_z0t_id, usersf_z0t_id, z0t_id)
Re = u_dyn0 * z0_m / nu_air
!write (*,*) surface_type, z0_m, z0m_id, z0t_id, z0_t, 'hi3'
!write (*,*) u_dyn0, 'u_dyn0'
call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
! --- define relative height
......
srcF/sfx_most.f90
+ 3
− 3
View file @ 54ab84cb
......@@ -117,7 +117,7 @@ contains
real Cm, Ct !< transfer coeff. for (momentum) & (heat) [n/d]
integer surface_type !< surface type = (ocean || land)
real z0_m1
#ifdef SFX_CHECK_NAN
real NaN
......@@ -144,7 +144,7 @@ contains
dT = meteo%dT
dQ = meteo%dQ
h = meteo%h
z0_m = meteo%z0_m
z0_m1 = meteo%z0_m
depth = meteo%depth
lai = meteo%lai
surface_type=meteo%surface_type
......@@ -152,7 +152,7 @@ contains
call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
forest_z0m_id, usersf_z0m_id, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m1, z0m_id)
call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
forest_z0t_id, usersf_z0t_id, z0t_id)
......
srcF/sfx_most_snow.f90
+ 55
− 11
View file @ 54ab84cb
......@@ -118,14 +118,15 @@ contains
real Pr_t_inv !< invese Prandt number [n/d]
real Cm, Ct !< transfer coeff. for (momentum) & (heat) [n/d]
real z0_m1
! --- local additional variables
! ----------------------------------------------------------------------------
!real phi_m, phi_m
real hfx, mfx
real S_mean, Lsnow
real z0_s, h_salt
!real S_salt
integer surface_type !< surface type = (ocean || land)
......@@ -158,7 +159,7 @@ contains
dT = meteo%dT
dQ = meteo%dQ
h = meteo%h
z0_m = meteo%z0_m
z0_m1 = meteo%z0_m
depth = meteo%depth
lai = meteo%lai
surface_type=meteo%surface_type
......@@ -167,7 +168,7 @@ contains
call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
forest_z0m_id, usersf_z0m_id, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m1, z0m_id)
call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
forest_z0t_id, usersf_z0t_id, z0t_id)
......@@ -175,7 +176,7 @@ contains
Re = u_dyn0 * z0_m / nu_air
call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
! --- define relative height
h0_m = h / z0_m
! --- define relative height [thermal]
......@@ -251,8 +252,8 @@ contains
real :: psi_m, psi_h
real :: psi0_m, psi0_h
real :: Linv
real :: w_snow, sigma_m
integer :: i
real :: w_snow, sigma_m, S_salt, d_s, a, b
integer :: i, n
! ----------------------------------------------------------------------------
......@@ -280,19 +281,20 @@ contains
Qdyn = kappa * dQ * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
if (Udyn < 1e-5) exit
Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn)
zeta = z * Linv
!S_salt =0.0004
call S_mean_fun(S_salt, Udyn)
if (Udyn>u_thsnow) then
call simpson_integration(d_s, 0.00000886, 0.0001, 1000);
call get_S_mean(S_mean, S_salt, h_s, z)
call get_sigma_m(sigma_m, rho_air, rho_s)
call get_w_snow(w_snow, sigma_m, g, d_s, nu_air)
Linv=Linv*((1-S_mean)/(1+sigma_m*S_mean))+(g*w_snow*sigma_m*S_mean/(Udyn**3.0))/(1+sigma_m*S_mean)
zeta = z * Linv
Lsnow=1/Linv
!write(*,*) S_mean, sigma_m, w_snow
write(*,*) S_mean, sigma_m, w_snow
!pause
!stop
......@@ -410,5 +412,47 @@ contains
end subroutine
! --------------------------------------------------------------------------------
subroutine S_mean_fun(S_salt, Udyn)
real, intent(in) :: Udyn
real, intent(out) :: S_salt
real S_salt1
if (Udyn <= u_thsnow) then
S_salt1 = 0.0
else
S_salt1 = (Udyn*Udyn - u_thsnow*u_thsnow) / (3.25 * Udyn * 9.8 * 0.08436 * Udyn**1.27)
end if
S_salt = S_salt1 / (S_salt1 + rho_s / rho_air)
end subroutine
subroutine simpson_integration(d_s, a, b, n)
real, intent(in) :: a, b
integer, intent(in) :: n
real, intent(out) :: d_s
real h, sum
integer i
h = (b - a) / n
sum=(exp(-0.5 * a * a)/sqrt(2 * 3.14))+(exp(-0.5 * b * b)/sqrt(2 * 3.14))
do i = 1, n - 1, 2
sum = sum + 4 * (exp(-0.5 * (a + i * h) * (a + i * h))/sqrt(2 * 3.14))
end do
do i = 2, n - 2, 2
sum = sum + 2 * (exp(-0.5 * (a + i * h) * (a + i * h))/sqrt(2 * 3.14))
end do
d_s = sum * (h / 3)
end subroutine
end module sfx_most_snow
\ No newline at end of file
srcF/sfx_most_snow_param.f90
+ 3
− 3
View file @ 54ab84cb
......@@ -29,10 +29,10 @@ module sfx_most_snow_param
!< snow parameters
real, parameter :: rho_s =900
real, parameter :: d_s=0.0000886
real, parameter :: h_s=0.07
!real, parameter :: d_s=0.00000886
real, parameter :: h_s=0.05
real, parameter :: u_thsnow=0.25
real, parameter :: S_salt=0.0004
!real, parameter :: S_salt=0.000004
real, parameter :: rho_air=1.2
end module sfx_most_snow_param
srcF/sfx_run.f90
+ 4
− 3
View file @ 54ab84cb
......@@ -89,9 +89,9 @@ contains
write(*, '(a,g0)') ' h = ', dataset%h
write(*, '(a,g0)') ' z0(m) = ', dataset%z0_m
write(*, '(a,g0)') ' z0(h) = ', dataset%z0_h
write(*, '(a,g0)') ' z0(h) = ', dataset%lai
write(*, '(a,g0)') ' z0(h) = ', dataset%depth
write(*, '(a,g0)') ' z0(h) = ', dataset%surface_type
write(*, '(a,g0)') ' lai = ', dataset%lai
write(*, '(a,g0)') ' depth = ', dataset%depth
write(*, '(a,g0)') ' surface_type = ', dataset%surface_type
! --- define number of elements
......@@ -430,6 +430,7 @@ contains
!< save nmax if previously set
nmax = dataset%nmax
call set_dataset(dataset, id)
dataset%nmax = nmax
call c_config_is_varname("dataset.filename"//C_NULL_CHAR, status)
......
srcF/sfx_surface.f90
+ 22
− 17
View file @ 54ab84cb
......@@ -48,8 +48,9 @@ module sfx_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_map_id = 3
integer, public, parameter :: z0m_user = 4
integer, public, parameter :: z0m_and = 5
integer, public, parameter :: z0t_kl_water = 0
......@@ -71,6 +72,8 @@ module sfx_surface
character(len = 16), parameter :: z0m_tag_ow = 'owen'
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 :: z0t_tag_kl_water = 'kl_water'
character(len = 16), parameter :: z0t_tag_kl_land = 'kl_land'
......@@ -87,16 +90,16 @@ module sfx_surface
integer, public, parameter :: ocean_z0m_id = z0m_ch !< ocean surface
integer, public, parameter :: land_z0m_id = z0m_map !< land 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 :: 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
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_map_id !< user surface
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 :: snow_z0t_id = z0t_kl_land !< snow covered surface
integer, public, parameter :: forest_z0t_id = z0t_du !< forest csurface
integer, public, parameter :: usersf_z0t_id = z0t_mix !< user surface
......@@ -134,7 +137,7 @@ module sfx_surface
!real, parameter :: B_max_lake = 8.0
contains
! surface type definition
! --------------------------------------------------------------------------------
function get_surface_id(tag) result(id)
......@@ -195,7 +198,7 @@ contains
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
z0m_id = z0m_map_id
else if (trim(tag_z0m) == trim(z0m_tag_user)) then
z0m_id = z0m_user
end if
......@@ -214,7 +217,7 @@ contains
tag_z0m = z0m_tag_fe
else if (z0m_id == z0m_ow) then
tag_z0m = z0m_tag_ow
else if (z0m_id == z0m_map) then
else if (z0m_id == z0m_map_id) then
tag_z0m = z0m_tag_map
else if (z0m_id == z0m_user) then
tag_z0m = z0m_tag_user
......@@ -252,7 +255,7 @@ contains
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)
......@@ -347,7 +350,7 @@ contains
else if (surface_type == surface_user) then
z0m_id = usersf_z0m_id
end if
!write (*,*) z0m_id, surface_type
end subroutine
! ----------------------------------------------------------------------------
......@@ -366,19 +369,21 @@ contains
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
else if (z0m_id == z0m_map_id) then
call get_dynamic_roughness_map(z0_m, u_dyn0, U, h, z0m_map)
else if (z0m_id == z0m_and) then
call get_dynamic_roughness_and(z0_m, u_dyn0, U, h, z0m_map)
else if (z0m_id == z0m_user) then
write(*, *) 'z0m_user'
end if
end subroutine
! --------------------------------------------------------------------------------
......@@ -413,7 +418,7 @@ subroutine get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t
else if (surface_type == surface_user) then
z0t_id = usersf_z0t_id
end if
end subroutine
!---------------------------------------------------------------
......@@ -433,7 +438,7 @@ subroutine get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t
! ---------------------------------------------------------------------------
real Czm !< proportionality coefficient z0_t =Czm*z0_m
Czm=1
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
......@@ -459,7 +464,7 @@ subroutine get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t
else if (z0t_id == z0t_user) then
write(*, *) 'z0t_user'
end if
end subroutine
......
srcF/sfx_z0m_all.f90
+ 15
− 14
View file @ 54ab84cb
......@@ -167,34 +167,35 @@ module sfx_z0m_all
! --------------------------------------------------------------------------------
! ----------------------------------------------------------------------------
subroutine get_dynamic_roughness_definition(surface_type, id_ocean, id_land, id_snow, id_lake, &
id_forest, id_user, z0m_id)
subroutine get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, snow_z0m_id, lake_z0m_id, &
forest_z0m_id, usersf_z0m_id, z0m_id)
! ----------------------------------------------------------------------------
real, intent(out) :: 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
real, intent(in) :: ocean_z0m_id
real, intent(in) :: land_z0m_id
real, intent(in) :: snow_z0m_id
real, intent(in) :: lake_z0m_id
real, intent(in) :: forest_z0m_id
real, intent(in) :: usersf_z0m_id
! ---------------------------------------------------------------------------
Write (*,*) 'get_dynamic_roughness_definition'
if (surface_type == surface_ocean) then
z0m_id = id_ocean
z0m_id = ocean_z0m_id
else if (surface_type == surface_land) then
z0m_id = id_land
z0m_id = land_z0m_id
else if (surface_type == surface_snow) then
z0m_id = id_snow
z0m_id = snow_z0m_id
else if (surface_type == surface_lake) then
z0m_id = id_lake
z0m_id = lake_z0m_id
else if (surface_type == surface_forest) then
z0m_id = id_forest
z0m_id = forest_z0m_id
else if (surface_type == surface_user) then
z0m_id = id_user
z0m_id = usersf_z0m_id
end if
end subroutine
......
srcF/sfx_z0m_all_surface.f90
+ 53
− 3
View file @ 54ab84cb
......@@ -3,7 +3,11 @@ module sfx_z0m_all_surface
use sfx_phys_const
implicit none
public
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
! --------------------------------------------------------------------------------
......@@ -117,7 +121,7 @@ module sfx_z0m_all_surface
u_dyn0 = Uc / c
! write(*,*) z0_m, 'oewn'
end subroutine
! --------------------------------------------------------------------------------
......@@ -196,11 +200,57 @@ subroutine get_dynamic_roughness_map(z0_m, u_dyn0, U, h, z0m_map)
z0_m=z0m_map
h0_m = h / z0_m
u_dyn0 = U * kappa / log(h0_m)
write(*,*) z0_m, 'map'
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
! --------------------------------------------------------------------------------
......
srcF/sfx_z0t_all.f90
+ 1
− 3
View file @ 54ab84cb
......@@ -9,7 +9,6 @@ module sfx_z0t_all
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
......@@ -218,7 +217,7 @@ module sfx_z0t_all
write(*, *) 'z0t_user'
end if
write (*,*) z0_t, B, z0_m, Re, u_dyn, Czm, LAI, z0t_id
end subroutine
! --------------------------------------------------------------------------------
......@@ -258,7 +257,6 @@ module sfx_z0t_all
else if (surface_type == surface_user) then
usersf_z0t_id = id_user
end if
end subroutine
end module sfx_z0t_all
srcF/sfx_z0t_all_surface.f90
+ 6
− 5
View file @ 54ab84cb
......@@ -46,11 +46,10 @@ module sfx_z0t_all_surface
! *: 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
! --------------------------------------------------------------------------------
......@@ -272,12 +271,14 @@ module sfx_z0t_all_surface
B=alog(-0.56*(4.0*(Re)**(0.5)-3.4)) !Repina, 2023
z0_t = z0_m*(-0.56*(4.0*(Re)**(0.5)-3.4))
B=alog(z0_m/z0_t)
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
!z0_t = z0_m / exp(B)
!write(*,*) z0_t, B, Re, 're'
end subroutine
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment