diff --git a/srcF/sfx_sheba_noniterative.f90 b/srcF/sfx_sheba_noniterative.f90 index 5c6462c5dc11fb032d2d4df9965576c6164e68b9..f76028ff54ff52b29b0c7fe93cfa48e0aa6e71b2 100644 --- a/srcF/sfx_sheba_noniterative.f90 +++ b/srcF/sfx_sheba_noniterative.f90 @@ -230,6 +230,7 @@ contains real sigma_w, sigma_r real S_mean, S_salt real z0_m1 + integer i #ifdef SFX_CHECK_NAN real NaN #endif @@ -288,36 +289,33 @@ contains h0_m, h0_t, B) ! --- 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) + !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 + !deltaS=S_salt-S_mean - Ri_sn = (g * sigma_r * deltaS * h) / U**2 - + !Ri_sn = (g * sigma_r * deltaS * h) / U**2 + Ri_sn=0.0 - ! --- get the fluxes ! ---------------------------------------------------------------------------- - !write(*,*) 'surface_type', surface_type if (Rib > 0.0) then ! --- stable stratification block ! --- restrict bulk Ri value ! *: note that this value is written to output ! Rib = min(Rib, Rib_max) - do i = 1, maxiters + do i = 1, numerics%maxiters_convection if (surface_type == surface_snow) then - write(*,*) 'sfx_snow1', Ri_sn, Rib, U, deltaS, i - !write(*,*) 'sfx_snow2', deltaS, u_dyn0, S_mean + !write(*,*) 'RIsnow', Rib, Ri_sn 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) call get_psi_stable(psi0_m, psi0_h, zeta * z0_m / h, zeta * z0_t / h) @@ -327,16 +325,20 @@ contains 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)) + + 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 - + Ri_sn = (g * sigma_r * deltaS * h) / U**2 end do - else if (Rib < Rib_conv_lim) then + !write(*,*) 'sfx_snow1', Ri_sn, Rib, Udyn, S_salt, S_mean + else if (Rib < Rib_conv_lim) then ! --- strong instability block call get_psi_convection(psi_m, psi_h, zeta, Rib, & @@ -384,10 +386,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 = 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) + 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) end subroutine get_surface_fluxes ! -------------------------------------------------------------------------------- @@ -418,6 +425,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 end subroutine subroutine get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt) @@ -430,6 +438,7 @@ contains 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) + !write(*,*) 'S_salt', S_salt else S_salt=0.0 endif