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

sheba noit snow add

parent b3038599
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
srcF/sfx_sheba_noit_param.f90
+ 10
− 1
View file @ 7fd94a0b
......@@ -36,5 +36,14 @@ module sfx_sheba_noit_param
real, parameter :: gamma = 2.91, zeta_a = 3.6 ! for stable psi
!< snow parameters
real, parameter :: rho_s =900
real, parameter :: d_s=0.0000886
real, parameter :: h_salt_const=0.05
real, parameter :: u_thsnow=0.35
real, parameter :: S_salt_const=0.000004
real, parameter :: rho_air=1.2
real, parameter :: Csn=3.25
end module sfx_sheba_noit_param
srcF/sfx_sheba_noniterative.f90
+ 110
− 30
View file @ 7fd94a0b
......@@ -28,6 +28,8 @@ module sfx_sheba_noniterative
public :: get_surface_fluxes
public :: get_surface_fluxes_vec
public :: get_psi_stable
integer z0m_id
integer z0t_id
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
......@@ -184,6 +186,8 @@ contains
real :: Tsemi !< semi-sum of potential temperature at 'h' and at surface [K]
real :: dQ !< difference between humidity at 'h' and at surface [g/g]
real :: z0_m !< surface aerodynamic roughness (should be < 0 for water bodies surface)
real :: depth
real :: lai
! ----------------------------------------------------------------------------
! --- local variables
......@@ -193,10 +197,11 @@ contains
real h0_m, h0_t !< = h / z0_m, h / z0_h [n/d]
real u_dyn0 !< dynamic velocity in neutral conditions [m/s]
real Re !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]
real Rib, Ri_sn !< bulk Richardson number, snow Richardson number
real zeta !< = z/L [n/d]
real Rib !< bulk Richardson number
real Re !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]
real zeta_conv_lim !< z/L critical value for matching free convection limit [n/d]
real Rib_conv_lim !< Ri-bulk critical value for matching free convection limit [n/d]
......@@ -220,7 +225,10 @@ contains
real fval !< just a shortcut for partial calculations
real w_snow, deltaS, h_salt
real sigma_w, sigma_r
real S_mean, S_salt
real z0_m1
#ifdef SFX_CHECK_NAN
real NaN
#endif
......@@ -246,44 +254,54 @@ 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
! --- define surface type
if (z0_m < 0.0) then
surface_type = surface_ocean
else
surface_type = surface_land
end if
if (surface_type == surface_ocean) then
! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
! --- define relative height
h0_m = h / z0_m
endif
if (surface_type == surface_land) then
! --- define relative height
h0_m = h / z0_m
! --- define dynamic velocity in neutral conditions
u_dyn0 = U * kappa / log(h0_m)
end if
!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, ice_z0m_id, 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, ice_z0t_id, z0t_id)
! --- define thermal roughness & B = log(z0_m / z0_h)
Re = u_dyn0 * z0_m / nu_air
call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
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]
h0_t = h / z0_t
! --- define snow consentration
call get_sigma(sigma_r, sigma_w, rho_air, rho_s)
call get_w_snow(w_snow, sigma_w, g, d_s, nu_air)
call get_h_salt(h_salt, u_dyn0)
call get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt)
call get_S_mean(S_mean, S_salt, h_salt, h, w_snow, u_dyn0)
deltaS=S_salt-S_mean
! --- define Ri-bulk
Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2
Ri_sn = (g * sigma_r * deltaS * h) / U**2
!write(*,*) 'ri', S_salt, h_salt, h, w_snow, u_dyn0
! --- define free convection transition zeta = z/L value
call get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, &
h0_m, h0_t, B)
! --- get the fluxes
! ----------------------------------------------------------------------------
write(*,*) 'surface_type', surface_type
if (Rib > 0.0) then
! --- stable stratification block
......@@ -291,6 +309,11 @@ contains
! *: note that this value is written to output
! Rib = min(Rib, Rib_max)
if (surface_type == surface_snow) then
write(*,*) 'snow', Ri_sn, Rib
Rib=Rib+Ri_sn
endif
call get_zeta(zeta, Rib, h, z0_m, z0_t)
call get_psi_stable(psi_m, psi_h, zeta, zeta)
......@@ -357,6 +380,63 @@ contains
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
!--------------------------------------snow functions-----------------------------
subroutine get_sigma(sigma_r, sigma_w, rho_air, rho_s)
!< @brief function for
! ----------------------------------------------------------------------------
real, intent(out) :: sigma_r, sigma_w !< s
real, intent(in) :: rho_air, rho_s
! ----------------------------------------------------------------------------
sigma_r = ((rho_s/rho_air)-1)
sigma_w = (rho_s-rho_air)/rho_air
end subroutine
subroutine get_h_salt(h_salt, u_dyn0)
! ----------------------------------------------------------------------------
real, intent(out) :: h_salt
real, intent(in) :: u_dyn0
! ----------------------------------------------------------------------------
h_salt = 0.08436*u_dyn0**1.27
end subroutine
subroutine get_S_mean(S_mean, S_salt, h_salt, z, w_snow, u_dyn0)
!< @brief function for snow consentration
! ----------------------------------------------------------------------------
real, intent(out) :: S_mean !< snow consentration
real, intent(in) :: S_salt, h_salt, z, w_snow, u_dyn0
! ----------------------------------------------------------------------------
S_mean = S_salt * (z/h_salt)**(-w_snow/(kappa*u_dyn0))
end subroutine
subroutine get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt)
!< @brief function for snow consentration
! ----------------------------------------------------------------------------
real, intent(out) :: S_salt !< snow consentration
real, intent(in) :: u_dyn0,u_thsnow, g, h_salt
! ----------------------------------------------------------------------------
real qs
if (u_dyn0>u_thsnow) then
qs = (u_dyn0*u_dyn0-u_thsnow*u_thsnow)/(Csn*u_dyn0*g*h_salt)
S_salt=qs/(qs+rho_s/rho_air)
else
S_salt=0.0
endif
end subroutine
subroutine get_w_snow(w_snow, sigma_w, g, d_s, nu_air)
!< @brief function for smow velosity
! ----------------------------------------------------------------------------
real, intent(out) :: w_snow !<
real, intent(in) :: sigma_w, g, d_s, nu_air
! ----------------------------------------------------------------------------
w_snow = (sigma_w * g * d_s * d_s) / (18.0 * nu_air);
end subroutine
!-----------------------------------------------------------------------------
subroutine get_zeta(zeta, Rib, h, z0_m, z0_t)
real,intent(out) :: zeta
......
srcF/sfx_surface.f90
+ 13
− 13
View file @ 7fd94a0b
......@@ -20,8 +20,8 @@ module sfx_surface
#if defined(INCLUDE_CXX)
public :: set_c_struct_sfx_surface_param_values
#endif
public :: get_charnock_roughness
public :: get_thermal_roughness
!public :: get_charnock_roughness
!public :: get_thermal_roughness
public :: get_dynamic_roughness_all
public :: get_dynamic_roughness_definition
public :: get_thermal_roughness_all
......@@ -96,22 +96,22 @@ module sfx_surface
character(len = 16), parameter :: z0t_tag_user = 'zt_user'
integer, public, parameter :: ocean_z0m_id = z0m_fe !< ocean 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 :: snow_z0m_id = z0m_and !< snow covered 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
integer, public, parameter :: usersf_z0m_id = z0m_map_id !< user surface
integer, public, parameter :: ice_z0m_id = z0m_map_id !< ice surface
integer, public, parameter :: ocean_z0t_id = z0t_ca !< ocean surface
integer, public, parameter :: land_z0t_id = z0t_ot !< land surface
integer, public, parameter :: lake_z0t_id = z0t_ot !< 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 :: usersf_z0t_id = z0t_kl_water !< user surface
integer, public, parameter :: ice_z0t_id = z0t_ca !< 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_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_land !< user surface
integer, public, parameter :: ice_z0t_id = z0t_kl_land !< user surface
! --------------------------------------------------------------------------------
real, parameter, private :: kappa = 0.40 !< von Karman constant [n/d]
......
srcF/sfx_z0m_all_surface.f90
+ 3
− 3
View file @ 7fd94a0b
......@@ -84,7 +84,7 @@ module sfx_z0m_all_surface
! --- define dynamic velocity in neutral conditions
u_dyn0 = Uc / c
write(*,*) 'ch', z0_m, u_dyn0
end subroutine
! --------------------------------------------------------------------------------
......@@ -125,7 +125,7 @@ module sfx_z0m_all_surface
u_dyn0 = Uc / c
! write(*,*) z0_m, 'oewn'
write(*,*) 'ow', z0_m, u_dyn0
end subroutine
! --------------------------------------------------------------------------------
......@@ -204,7 +204,7 @@ 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(*,*) 'map', z0_m, u_dyn0
end subroutine
! --------------------------------------------------------------------------------
......
srcF/sfx_z0t_all_surface.f90
+ 2
− 2
View file @ 7fd94a0b
......@@ -49,7 +49,7 @@ module sfx_z0t_all_surface
B = min(B, B_max_land)
z0_t = z0_m / exp(B)
write(*,*) 'kl_land', z0_t, B
end subroutine
! --------------------------------------------------------------------------------
......@@ -75,7 +75,7 @@ module sfx_z0t_all_surface
B = min(B, B_max_ocean)
z0_t = z0_m / exp(B)
write(*,*) 'kl_water', z0_t, B
end subroutine
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment