diff --git a/srcF/sfx_sheba_noniterative.f90 b/srcF/sfx_sheba_noniterative.f90
index a2dbcbcb616f5b695623b8db5603b3262424cfec..8d87893e615e1a12c284e0b8ae3ed618c1f0f26e 100644
--- a/srcF/sfx_sheba_noniterative.f90
+++ b/srcF/sfx_sheba_noniterative.f90
@@ -228,7 +228,7 @@ contains
real w_snow, deltaS, h_salt
real sigma_w, sigma_r
- real S_mean, S_salt
+ real S_mean, S_salt, Linv
real z0_m1
integer i
#ifdef SFX_CHECK_NAN
@@ -315,29 +315,50 @@ contains
!endif
call get_zeta(zeta, Rib, h, z0_m, z0_t)
+ !write(*,*) 'get_zeta', zeta, h
call get_psi_stable(psi_m, psi_h, zeta, zeta)
call get_psi_stable(psi0_m, psi0_h, zeta * z0_m / h, zeta * z0_t / h)
phi_m = 1.0 + (a_m * zeta * (1.0 + zeta)**(1.0 / 3.0)) / (1.0 + b_m * zeta)
phi_h = 1.0 + (a_h * zeta + b_h * zeta * zeta) / (1.0 + c_h * zeta + zeta * zeta)
-
+ S_mean=0.0
Udyn = kappa * U / (log(h / z0_m) - (psi_m - psi0_m))
Tdyn = kappa * dT * Pr_t_0_inv / (log(h / z0_t) - (psi_h - psi0_h))
+ !write(*,*) 'sfx_before_snow', Udyn, zeta, S_mean
+
+ if (surface_type==3.or.surface_type==6.) then
+ if (Udyn>u_thsnow) then
-
- ! 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, Udyn)
- ! call get_S_salt(S_salt, Udyn, u_thsnow, g, h_salt)
- ! call get_S_mean(S_mean, S_salt, h_salt, h, w_snow, Udyn)
+ 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, Udyn)
+ call get_S_salt(S_salt, Udyn, u_thsnow, g, h_salt)
+ call get_S_mean(S_mean, S_salt, h_salt, h, w_snow, Udyn)
- ! deltaS=S_salt-S_mean
-
+ Linv=Linv*((1-S_mean)/(1+sigma_w*S_mean))+(g*w_snow*sigma_w*S_mean/(Udyn**3.0))/(1+sigma_w*S_mean)
+ zeta = h * Linv
+ !write(*,*) 'sfx_snow1', Udyn, Linv, zeta, S_mean
+
+
+
+ call get_psi_stable(psi_m, psi_h, zeta, zeta)
+ call get_psi_stable(psi0_m, psi0_h, zeta * z0_m / h, zeta * z0_t / h)
+
+ phi_m = 1.0 + (a_m * zeta * (1.0 + zeta)**(1.0 / 3.0)) / (1.0 + b_m * zeta)
+ phi_h = 1.0 + (a_h * zeta + b_h * zeta * zeta) / (1.0 + c_h * zeta + zeta * zeta)
+
+ Udyn = kappa * U / (log(h / z0_m) - (psi_m - psi0_m))
+ Tdyn = kappa * dT * Pr_t_0_inv / (log(h / z0_t) - (psi_h - psi0_h))
+
+ write(*,*) 'sfx_snow2', Udyn, zeta, S_mean, Linv
+ endif
+ endif
+
! Ri_sn = (g * sigma_r * deltaS * h) / U**2
!end do
- !write(*,*) 'sfx_snow', Ri_sn, Rib, Udyn, surface_type
+
else if (Rib < Rib_conv_lim) then
! --- strong instability block
@@ -386,15 +407,15 @@ contains
Pr_t_inv = phi_m / phi_h
! --- setting output
- sfx = sfxDataType(zeta = zeta, Rib = Rib, &
- Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
- Rib_conv_lim = Rib_conv_lim, &
- Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
-
!sfx = sfxDataType(zeta = zeta, Rib = Rib, &
- ! Re = Ri_sn, B = B, z0_m = z0_m, z0_t = z0_t, &
- ! Rib_conv_lim = S_mean, &
- ! Cm = Cm, Ct = Ct, Km = S_salt, Pr_t_inv = Udyn)
+ ! Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
+ ! Rib_conv_lim = Rib_conv_lim, &
+ ! Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
+ !write(*,*) 'Smean_0ut', S_mean
+ sfx = sfxDataType(zeta = zeta, Rib = Rib, &
+ Re = Linv, B = B, z0_m = z0_m, z0_t = z0_t, &
+ Rib_conv_lim = S_mean, &
+ Cm = Cm, Ct = Ct, Km = S_salt, Pr_t_inv = Udyn)
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
@@ -425,7 +446,7 @@ contains
real, intent(in) :: S_salt, h_salt, z, w_snow, u_dyn0
! ----------------------------------------------------------------------------
S_mean = S_salt * (z/h_salt)**(-w_snow/(kappa*u_dyn0))
- !write(*,*) 'Smean_Udyn', S_mean, u_dyn0
+ !write(*,*) 'Smean_func', S_mean, u_dyn0
end subroutine
subroutine get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt)
diff --git a/srcF/sfx_z0m_all_surface.f90 b/srcF/sfx_z0m_all_surface.f90
index 249ce699569617e050b827cfec42e57894756620..183403cfaff2ff31ff8133af9c7c90f0a6a6df86 100644
--- a/srcF/sfx_z0m_all_surface.f90
+++ b/srcF/sfx_z0m_all_surface.f90
@@ -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
+ !write(*,*) 'ch', 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
+ !write(*,*) 'map', z0_m, u_dyn0
end subroutine
! --------------------------------------------------------------------------------
diff --git a/srcF/sfx_z0t_all_surface.f90 b/srcF/sfx_z0t_all_surface.f90
index c6581f3b82f83bed460de75bfed918af5a7be352..88d2e0ae66fe93c297150815fc6137824556b027 100644
--- a/srcF/sfx_z0t_all_surface.f90
+++ b/srcF/sfx_z0t_all_surface.f90
@@ -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
+ !write(*,*) 'kl_land', z0_t, B
end subroutine
! --------------------------------------------------------------------------------