diff --git a/srcF/sfx_most_snow.f90 b/srcF/sfx_most_snow.f90 new file mode 100644 index 0000000000000000000000000000000000000000..5ff26f84a2fc617ab6334218a9db0354dd5744b8 --- /dev/null +++ b/srcF/sfx_most_snow.f90 @@ -0,0 +1,431 @@ +#include "../includeF/sfx_def.fi" + +module sfx_most_snow + !< @brief MOST surface flux module + + ! modules used + ! -------------------------------------------------------------------------------- +#ifdef SFX_CHECK_NAN + use sfx_common +#endif + use sfx_data + use sfx_surface + use sfx_most_snow_param + ! -------------------------------------------------------------------------------- + + ! directives list + ! -------------------------------------------------------------------------------- + implicit none + private + ! -------------------------------------------------------------------------------- + + ! public interface + ! -------------------------------------------------------------------------------- + public :: get_surface_fluxes + public :: get_surface_fluxes_vec + ! -------------------------------------------------------------------------------- + + ! -------------------------------------------------------------------------------- + type, public :: numericsType + integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness + integer :: maxiters_snow = 10 !< maximum (actual) number of iterations in snow roughness + end type + ! -------------------------------------------------------------------------------- + +contains + + ! -------------------------------------------------------------------------------- + subroutine get_surface_fluxes_vec(sfx, sfx2, meteo, numerics, n) + !< @brief surface flux calculation for array data + !< @details contains C/C++ & CUDA interface + ! ---------------------------------------------------------------------------- + type (sfxDataVecType), intent(inout) :: sfx + type (sfxDataAddVecType), intent(inout) :: sfx2 + + type (meteoDataVecType), intent(in) :: meteo + type (numericsType), intent(in) :: numerics + integer, intent(in) :: n + ! ---------------------------------------------------------------------------- + + ! --- local variables + type (meteoDataType) meteo_cell + type (sfxDataType) sfx_cell + type (sfxDataAddType) sfx2_cell + integer i + ! ---------------------------------------------------------------------------- + + do i = 1, n + meteo_cell = meteoDataType(& + h = meteo%h(i), & + U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), & + z0_m = meteo%z0_m(i)) + + call get_surface_fluxes(sfx_cell, sfx2_cell, meteo_cell, numerics) + + call push_sfx_data(sfx, sfx_cell, i) + call push_sfx_data_add(sfx2, sfx2_cell, i) + end do + + end subroutine get_surface_fluxes_vec + ! -------------------------------------------------------------------------------- + + ! -------------------------------------------------------------------------------- + subroutine get_surface_fluxes(sfx, sfx2, meteo, numerics) + !< @brief surface flux calculation for single cell + !< @details contains C/C++ interface + ! ---------------------------------------------------------------------------- +#ifdef SFX_CHECK_NAN + use ieee_arithmetic +#endif + + type (sfxDataType), intent(out) :: sfx + type (sfxDataAddType), intent(out) :: sfx2 + type (meteoDataType), intent(in) :: meteo + type (numericsType), intent(in) :: numerics + ! ---------------------------------------------------------------------------- + + ! --- meteo derived datatype name shadowing + ! ---------------------------------------------------------------------------- + real :: h !< constant flux layer height [m] + real :: U !< abs(wind speed) at 'h' [m/s] + real :: dT !< difference between potential temperature at 'h' and at surface [K] + 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) + ! ---------------------------------------------------------------------------- + + ! --- local variables + ! ---------------------------------------------------------------------------- + real z0_t !< thermal roughness [m] + real B !< = ln(z0_m / z0_t) [n/d] + 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 zeta !< = z/L [n/d] + real Rib !< bulk Richardson number + + real Udyn, Tdyn, Qdyn !< dynamic scales + + real phi_m, phi_h !< stability functions (momentum) & (heat) [n/d] + + real Km !< eddy viscosity coeff. at h [m^2/s] + real Pr_t_inv !< invese Prandt number [n/d] + + real Cm, Ct !< transfer coeff. for (momentum) & (heat) [n/d] + + ! --- local additional variables + ! ---------------------------------------------------------------------------- + !real phi_m, phi_m + real hfx, mfx + real S_mean, Lsnow + real z0_s, h_salt + + + integer surface_type !< surface type = (ocean || land) + + +#ifdef SFX_CHECK_NAN + real NaN +#endif + ! ---------------------------------------------------------------------------- + +#ifdef SFX_CHECK_NAN + ! --- checking if arguments are finite + if (.not.(is_finite(meteo%U).and.is_finite(meteo%Tsemi).and.is_finite(meteo%dT).and.is_finite(meteo%dQ) & + .and.is_finite(meteo%z0_m).and.is_finite(meteo%h))) then + + !NaN = ieee_value(0.0, ieee_quiet_nan) ! setting NaN + sfx = sfxDataType(zeta = NaN, Rib = NaN, & + Re = NaN, B = NaN, z0_m = NaN, z0_t = NaN, & + Rib_conv_lim = NaN, & + Cm = NaN, Ct = NaN, Km = NaN, Pr_t_inv = NaN) + + sfx2 = sfxDataAddType(phi_m = NaN, phi_h = NaN, & + hfx = NaN, mfx = NaN, Udyn = NaN, S_mean = NaN, & + Lsnow = NaN, & + z0_s = NaN, h_salt = NaN) + + + return + end if +#endif + + ! --- shadowing names for clarity + U = meteo%U + Tsemi = meteo%Tsemi + dT = meteo%dT + dQ = meteo%dQ + h = meteo%h + z0_m = meteo%z0_m + + ! --- define surface type + if (z0_m < 0.0) then + surface_type = surface_ocean + else if (z0_m == 0.0) then + surface_type = surface_snow + 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_snow) then + ! --- define surface roughness [momentum] & dynamic velocity in neutral conditions + call get_snow_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_snow) + ! --- 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 + + ! --- 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, u_dyn0) + + ! --- define relative height [thermal] + h0_t = h / z0_t + + ! --- define Ri-bulk + Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2 + + + ! --- get the fluxes + ! ---------------------------------------------------------------------------- + call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, & + Lsnow, S_mean, h_salt, & + U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10) + ! ---------------------------------------------------------------------------- + + call get_phi(phi_m, phi_h, zeta) + ! ---------------------------------------------------------------------------- + + ! --- define transfer coeff. (momentum) & (heat) + Cm = 0.0 + if (U > 0.0) then + Cm = Udyn / U + end if + Ct = 0.0 + if (abs(dT) > 0.0) then + Ct = Tdyn / dT + end if + + ! --- define eddy viscosity & inverse Prandtl number + Km = kappa * Cm * U * h / phi_m + Pr_t_inv = phi_m / phi_h + ! --- define heat flux and momentum flux + hfx=-rho_air*U*dT*Cm*Ct + mfx=-rho_air*Cm*Cm*U*U + h_salt=h_s + + ! --- setting output + sfx = sfxDataType(zeta = zeta, Rib = Rib, & + Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, & + Rib_conv_lim = 0.0, & + Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv) + ! --- setting additional output + sfx2 = sfxDataAddType(phi_m = phi_m, phi_h = phi_h, & + hfx = hfx, mfx = mfx, Udyn = Udyn, S_mean = S_mean, & + Lsnow = Lsnow, & + z0_s = z0_m, h_salt = h_s) + + end subroutine get_surface_fluxes + ! -------------------------------------------------------------------------------- + + !< @brief get dynamic scales + ! -------------------------------------------------------------------------------- + subroutine get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, & + Lsnow, S_mean, h_salt, & + U, Tsemi, dT, dQ, z, z0_m, z0_t, beta, maxiters) + ! ---------------------------------------------------------------------------- + real, intent(out) :: Udyn, Tdyn, Qdyn !< dynamic scales + real, intent(out) :: zeta !< = z/L + real, intent(out) :: Lsnow, S_mean + real, intent(out) :: h_salt + + real, intent(in) :: U !< abs(wind speed) at z + real, intent(in) :: Tsemi !< semi-sum of temperature at z and at surface + real, intent(in) :: dT, dQ !< temperature & humidity difference between z and at surface + real, intent(in) :: z !< constant flux layer height + real, intent(in) :: z0_m, z0_t !< roughness parameters + real, intent(in) :: beta !< buoyancy parameter + + integer, intent(in) :: maxiters !< maximum number of iterations + ! ---------------------------------------------------------------------------- + + ! --- local variables + real, parameter :: gamma = 0.61 + + real :: psi_m, psi_h + real :: psi0_m, psi0_h + real :: Linv + real :: w_snow, sigma_m + integer :: i + ! ---------------------------------------------------------------------------- + + + Udyn = kappa * U / log(z / z0_m) + Tdyn = kappa * dT * Pr_t_0_inv / log(z / z0_t) + Qdyn = kappa * dQ * Pr_t_0_inv / log(z / z0_t) + zeta = 0.0 + + ! --- no wind + if (Udyn < 1e-5) return + + Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn) + zeta = z * Linv + + ! --- near neutral case + if (Linv < 1e-5) return + + do i = 1, maxiters + + call get_psi(psi_m, psi_h, zeta) + call get_psi_mh(psi0_m, psi0_h, z0_m * Linv, z0_t * Linv) + + Udyn = kappa * U / (log(z / z0_m) - (psi_m - psi0_m)) + Tdyn = kappa * dT * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h)) + 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 + + if (Udyn>u_thsnow) then + + 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 + !pause + !stop + + endif + + + + end do + + end subroutine get_dynamic_scales + ! -------------------------------------------------------------------------------- + + ! stability functions + ! -------------------------------------------------------------------------------- + subroutine get_phi(phi_m, phi_h, zeta) + !< @brief stability functions (momentum) & (heat): neutral case + ! ---------------------------------------------------------------------------- + real, intent(out) :: phi_m, phi_h !< stability functions + + real, intent(in) :: zeta !< = z/L + ! ---------------------------------------------------------------------------- + + + if (zeta >= 0.0) then + phi_m = 1.0 + beta_m * zeta + phi_h = 1.0 + beta_h * zeta + else + phi_m = (1.0 - alpha_m * zeta)**(-0.25) + phi_h = (1.0 - alpha_h * zeta)**(-0.5) + end if + + end subroutine + ! -------------------------------------------------------------------------------- + subroutine get_S_mean(S_mean, S_salt, h_s, z) + !< @brief function for snow consentration + ! ---------------------------------------------------------------------------- + real, intent(out) :: S_mean !< snow consentration + real, intent(in) :: S_salt, h_s, z + ! ---------------------------------------------------------------------------- + S_mean = S_salt * h_s/z + end subroutine + ! -------------- + subroutine get_sigma_m(sigma_m, rho_air, rho_s) + !< @brief function for + ! ---------------------------------------------------------------------------- + real, intent(out) :: sigma_m !< s + real, intent(in) :: rho_air, rho_s + ! ---------------------------------------------------------------------------- + sigma_m = (rho_s - rho_air)/rho_air + end subroutine + subroutine get_w_snow(w_snow, sigma_m, g, d_s, nu_air) + !< @brief function for smow velosity + ! ---------------------------------------------------------------------------- + real, intent(out) :: w_snow !< + real, intent(in) :: sigma_m, g, d_s, nu_air + ! ---------------------------------------------------------------------------- + w_snow = sigma_m * g * d_s * d_s / (18.0 * nu_air); + + end subroutine + ! universal functions + ! -------------------------------------------------------------------------------- + subroutine get_psi(psi_m, psi_h, zeta) + !< @brief universal functions (momentum) & (heat): neutral case + ! ---------------------------------------------------------------------------- + real, intent(out) :: psi_m, psi_h !< universal functions + + real, intent(in) :: zeta !< = z/L + ! ---------------------------------------------------------------------------- + + ! --- local variables + real :: x_m, x_h + ! ---------------------------------------------------------------------------- + + + if (zeta >= 0.0) then + psi_m = -beta_m * zeta; + psi_h = -beta_h * zeta; + else + x_m = (1.0 - alpha_m * zeta)**(0.25) + x_h = (1.0 - alpha_h * zeta)**(0.25) + + psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m) + psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h)) + end if + + end subroutine + + + subroutine get_psi_mh(psi_m, psi_h, zeta_m, zeta_h) + !< @brief universal functions (momentum) & (heat): neutral case + ! ---------------------------------------------------------------------------- + real, intent(out) :: psi_m, psi_h !< universal functions + + real, intent(in) :: zeta_m, zeta_h !< = z/L + ! ---------------------------------------------------------------------------- + + ! --- local variables + real :: x_m, x_h + ! ---------------------------------------------------------------------------- + + + if (zeta_m >= 0.0) then + psi_m = -beta_m * zeta_m + else + x_m = (1.0 - alpha_m * zeta_m)**(0.25) + psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m) + end if + + if (zeta_h >= 0.0) then + psi_h = -beta_h * zeta_h; + else + x_h = (1.0 - alpha_h * zeta_h)**(0.25) + psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h)) + end if + + end subroutine + ! -------------------------------------------------------------------------------- + +end module sfx_most_snow \ No newline at end of file diff --git a/srcF/sfx_most_snow_param.f90 b/srcF/sfx_most_snow_param.f90 new file mode 100644 index 0000000000000000000000000000000000000000..c72f10ced323731ace8d36dd3d1c168ff2e24aa8 --- /dev/null +++ b/srcF/sfx_most_snow_param.f90 @@ -0,0 +1,38 @@ +module sfx_most_snow_param + !< @brief MOST BD71 surface flux model parameters + !< @details all in SI units + + ! modules used + ! -------------------------------------------------------------------------------- + use sfx_phys_const + ! -------------------------------------------------------------------------------- + + ! directives list + ! -------------------------------------------------------------------------------- + implicit none + ! -------------------------------------------------------------------------------- + + + !< von Karman constant [n/d] + real, parameter :: kappa = 0.40 + !< inverse Prandtl (turbulent) number in neutral conditions [n/d] + real, parameter :: Pr_t_0_inv = 1.15 + + + !< stability function coeff. (unstable) + real, parameter :: alpha_m = 16.0 + real, parameter :: alpha_h = 16.0 + + !< stability function coeff. (stable) + real, parameter :: beta_m = 4.7 + real, parameter :: beta_h = beta_m + + !< snow parameters + real, parameter :: rho_s =900 + real, parameter :: d_s=0.0000886 + real, parameter :: h_s=0.07 + real, parameter :: u_thsnow=0.25 + real, parameter :: S_salt=0.0004 + real, parameter :: rho_air=1.2 + +end module sfx_most_snow_param