#include "../includeF/sfx_def.fi" module sfx_esm !< @brief main Earth System Model surface flux module ! modules used ! -------------------------------------------------------------------------------- #ifdef SFX_CHECK_NAN use sfx_common #endif use sfx_data use sfx_surface use sfx_esm_param #if defined(INCLUDE_CXX) use iso_c_binding, only: C_LOC, C_PTR use C_FUNC #endif ! -------------------------------------------------------------------------------- ! directives list ! -------------------------------------------------------------------------------- implicit none private ! -------------------------------------------------------------------------------- ! public interface ! -------------------------------------------------------------------------------- public :: get_surface_fluxes public :: get_surface_fluxes_vec ! -------------------------------------------------------------------------------- ! -------------------------------------------------------------------------------- type, public :: numericsType integer :: maxiters_convection = 10 !< maximum (actual) number of iterations in convection integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness end type ! -------------------------------------------------------------------------------- contains ! -------------------------------------------------------------------------------- #if defined(INCLUDE_CXX) subroutine set_c_struct_sfx_esm_param_values(sfx_model_param) type (sfx_esm_param), intent(inout) :: sfx_model_param sfx_model_param%kappa = kappa sfx_model_param%Pr_t_0_inv = Pr_t_0_inv sfx_model_param%Pr_t_inf_inv = Pr_t_inf_inv sfx_model_param%alpha_m = alpha_m sfx_model_param%alpha_h = alpha_h sfx_model_param%alpha_h_fix = alpha_h_fix sfx_model_param%beta_m = beta_m sfx_model_param%beta_h = beta_h sfx_model_param%Rib_max = Rib_max end subroutine set_c_struct_sfx_esm_param_values #endif subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n) !< @brief surface flux calculation for array data !< @details contains C/C++ & CUDA interface ! ---------------------------------------------------------------------------- type (sfxDataVecType), intent(inout) :: sfx type (meteoDataVecType), intent(in) :: meteo type (numericsType), intent(in) :: numerics integer, intent(in) :: n ! ---------------------------------------------------------------------------- ! --- local variables type (meteoDataType) meteo_cell type (sfxDataType) sfx_cell integer i ! ---------------------------------------------------------------------------- #if defined(INCLUDE_CXX) type (meteoDataVecTypeC), target :: meteo_c !< meteorological data (input) type (sfxDataVecTypeC), target :: sfx_c !< surface flux data (output) type(C_PTR) :: meteo_c_ptr, sfx_c_ptr type (sfx_esm_param) :: model_param type (sfx_surface_param) :: surface_param type (sfx_esm_numericsTypeC) :: numerics_c type (sfx_phys_constants) :: phys_constants numerics_c%maxiters_convection = numerics%maxiters_convection numerics_c%maxiters_charnock = numerics%maxiters_charnock phys_constants%Pr_m = Pr_m; phys_constants%nu_air = nu_air; phys_constants%g = g; call set_c_struct_sfx_esm_param_values(model_param) call set_c_struct_sfx_surface_param_values(surface_param) call set_meteo_vec_c(meteo, meteo_c) call set_sfx_vec_c(sfx, sfx_c) meteo_c_ptr = C_LOC(meteo_c) sfx_c_ptr = C_LOC(sfx_c) call get_surface_fluxes_esm(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n) #else do i = 1, n #ifdef SFX_FORCE_DEPRECATED_ESM_CODE #else 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, meteo_cell, numerics) call push_sfx_data(sfx, sfx_cell, i) #endif end do #endif end subroutine get_surface_fluxes_vec ! -------------------------------------------------------------------------------- ! -------------------------------------------------------------------------------- subroutine get_surface_fluxes(sfx, 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 (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 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] real f_m_conv_lim !< stability function (momentum) value in free convection limit [n/d] real f_h_conv_lim !< stability function (heat) value in free convection limit [n/d] real psi_m, psi_h !< universal functions (momentum) & (heat) [n/d] 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] integer surface_type !< surface type = (ocean || land) real fval !< just a shortcut for partial calculations #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) 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 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 ! --- 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) ! --- define relative height [thermal] h0_t = h / z0_t ! --- define Ri-bulk Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2 ! --- 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 ! ---------------------------------------------------------------------------- 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) call get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B) fval = beta_m * zeta phi_m = 1.0 + fval phi_h = 1.0/Pr_t_0_inv + fval else if (Rib < Rib_conv_lim) then ! --- strong instability block call get_psi_convection(psi_m, psi_h, zeta, Rib, & zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, h0_m, h0_t, B, numerics%maxiters_convection) fval = (zeta_conv_lim / zeta)**(1.0/3.0) phi_m = fval / f_m_conv_lim phi_h = fval / (Pr_t_0_inv * f_h_conv_lim) else if (Rib > -0.001) then ! --- nearly neutral [-0.001, 0] block call get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B) zeta = 0.0 phi_m = 1.0 phi_h = 1.0 / Pr_t_0_inv else ! --- weak & semistrong instability block call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection) phi_m = (1.0 - alpha_m * zeta)**(-0.25) phi_h = 1.0 / (Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta)) end if ! ---------------------------------------------------------------------------- ! --- define transfer coeff. (momentum) & (heat) Cm = kappa / psi_m Ct = kappa / psi_h ! --- define eddy viscosity & inverse Prandtl number Km = kappa * Cm * U * h / phi_m 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) end subroutine get_surface_fluxes ! -------------------------------------------------------------------------------- ! convection universal functions shortcuts ! -------------------------------------------------------------------------------- function f_m_conv(zeta) ! ---------------------------------------------------------------------------- real :: f_m_conv real, intent(in) :: zeta ! ---------------------------------------------------------------------------- f_m_conv = (1.0 - alpha_m * zeta)**0.25 end function f_m_conv function f_h_conv(zeta) ! ---------------------------------------------------------------------------- real :: f_h_conv real, intent(in) :: zeta ! ---------------------------------------------------------------------------- f_h_conv = (1.0 - alpha_h * zeta)**0.5 end function f_h_conv ! -------------------------------------------------------------------------------- ! universal functions ! -------------------------------------------------------------------------------- subroutine get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B) !< @brief universal functions (momentum) & (heat): neutral case ! ---------------------------------------------------------------------------- real, intent(out) :: psi_m, psi_h !< universal functions real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h real, intent(in) :: B !< = log(z0_m / z0_h) ! ---------------------------------------------------------------------------- psi_m = log(h0_m) psi_h = log(h0_t) / Pr_t_0_inv !*: this looks redundant z0_t = z0_m in case |B| ~ 0 if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv end subroutine subroutine get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B) !< @brief universal functions (momentum) & (heat): stable case ! ---------------------------------------------------------------------------- real, intent(out) :: psi_m, psi_h !< universal functions [n/d] real, intent(out) :: zeta !< = z/L [n/d] real, intent(in) :: Rib !< bulk Richardson number [n/d] real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d] real, intent(in) :: B !< = log(z0_m / z0_h) [n/d] ! ---------------------------------------------------------------------------- ! --- local variables real :: Rib_coeff real :: psi0_m, psi0_h real :: phi real :: c ! ---------------------------------------------------------------------------- psi0_m = alog(h0_m) psi0_h = B / psi0_m Rib_coeff = beta_m * Rib c = (psi0_h + 1.0) / Pr_t_0_inv - 2.0 * Rib_coeff zeta = psi0_m * (sqrt(c**2 + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / & (2.0 * beta_m * (1.0 - Rib_coeff)) phi = beta_m * zeta psi_m = psi0_m + phi psi_h = (psi0_m + B) / Pr_t_0_inv + phi end subroutine subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters) !< @brief universal functions (momentum) & (heat): semi-strong convection case ! ---------------------------------------------------------------------------- real, intent(out) :: psi_m, psi_h !< universal functions [n/d] real, intent(out) :: zeta !< = z/L [n/d] real, intent(in) :: Rib !< bulk Richardson number [n/d] real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d] real, intent(in) :: B !< = log(z0_m / z0_h) [n/d] integer, intent(in) :: maxiters !< maximum number of iterations ! --- local variables real :: zeta0_m, zeta0_h real :: f0_m, f0_h real :: f_m, f_h integer :: i ! ---------------------------------------------------------------------------- psi_m = log(h0_m) psi_h = log(h0_t) if (abs(B) < 1.0e-10) psi_h = psi_m zeta = Rib * Pr_t_0_inv * psi_m**2 / psi_h do i = 1, maxiters + 1 zeta0_m = zeta / h0_m zeta0_h = zeta / h0_t if (abs(B) < 1.0e-10) zeta0_h = zeta0_m f_m = (1.0 - alpha_m * zeta)**0.25e0 f_h = sqrt(1.0 - alpha_h_fix * zeta) f0_m = (1.0 - alpha_m * zeta0_m)**0.25e0 f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h) f0_m = max(f0_m, 1.000001e0) f0_h = max(f0_h, 1.000001e0) psi_m = log((f_m - 1.0e0)*(f0_m + 1.0e0)/((f_m + 1.0e0)*(f0_m - 1.0e0))) + 2.0e0*(atan(f_m) - atan(f0_m)) psi_h = log((f_h - 1.0e0)*(f0_h + 1.0e0)/((f_h + 1.0e0)*(f0_h - 1.0e0))) / Pr_t_0_inv ! *: don't update zeta = z/L at last iteration if (i == maxiters + 1) exit zeta = Rib * psi_m**2 / psi_h end do end subroutine subroutine get_psi_convection(psi_m, psi_h, zeta, Rib, & zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, & h0_m, h0_t, B, maxiters) !< @brief universal functions (momentum) & (heat): fully convective case ! ---------------------------------------------------------------------------- real, intent(out) :: psi_m, psi_h !< universal functions [n/d] real, intent(out) :: zeta !< = z/L [n/d] real, intent(in) :: Rib !< bulk Richardson number [n/d] real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d] real, intent(in) :: B !< = log(z0_m / z0_h) [n/d] integer, intent(in) :: maxiters !< maximum number of iterations real, intent(in) :: zeta_conv_lim !< convective limit zeta real, intent(in) :: f_m_conv_lim, f_h_conv_lim !< universal function shortcuts in limiting case ! ---------------------------------------------------------------------------- ! --- local variables real :: zeta0_m, zeta0_h real :: f0_m, f0_h real :: p_m, p_h real :: a_m, a_h real :: c_lim, f integer :: i ! ---------------------------------------------------------------------------- p_m = 2.0 * atan(f_m_conv_lim) + log((f_m_conv_lim - 1.0) / (f_m_conv_lim + 1.0)) p_h = log((f_h_conv_lim - 1.0) / (f_h_conv_lim + 1.0)) zeta = zeta_conv_lim do i = 1, maxiters + 1 zeta0_m = zeta / h0_m zeta0_h = zeta / h0_t if (abs(B) < 1.0e-10) zeta0_h = zeta0_m f0_m = (1.0 - alpha_m * zeta0_m)**0.25 f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h) a_m = -2.0*atan(f0_m) + log((f0_m + 1.0)/(f0_m - 1.0)) a_h = log((f0_h + 1.0)/(f0_h - 1.0)) c_lim = (zeta_conv_lim / zeta)**(1.0 / 3.0) f = 3.0 * (1.0 - c_lim) psi_m = f / f_m_conv_lim + p_m + a_m psi_h = (f / f_h_conv_lim + p_h + a_h) / Pr_t_0_inv ! *: don't update zeta = z/L at last iteration if (i == maxiters + 1) exit zeta = Rib * psi_m**2 / psi_h end do end subroutine ! -------------------------------------------------------------------------------- ! convection limit definition ! -------------------------------------------------------------------------------- subroutine get_convection_lim(zeta_lim, Rib_lim, f_m_lim, f_h_lim, & h0_m, h0_t, B) ! ---------------------------------------------------------------------------- real, intent(out) :: zeta_lim !< limiting value of z/L real, intent(out) :: Rib_lim !< limiting value of Ri-bulk real, intent(out) :: f_m_lim, f_h_lim !< limiting values of universal functions shortcuts real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d] real, intent(in) :: B !< = log(z0_m / z0_h) [n/d] ! ---------------------------------------------------------------------------- ! --- local variables real :: psi_m, psi_h real :: f_m, f_h real :: c ! ---------------------------------------------------------------------------- ! --- define limiting value of zeta = z / L c = (Pr_t_inf_inv / Pr_t_0_inv)**4 zeta_lim = (2.0 * alpha_h - c * alpha_m - & sqrt((c * alpha_m)**2 + 4.0 * c * alpha_h * (alpha_h - alpha_m))) / (2.0 * alpha_h**2) f_m_lim = f_m_conv(zeta_lim) f_h_lim = f_h_conv(zeta_lim) ! --- universal functions f_m = zeta_lim / h0_m f_h = zeta_lim / h0_t if (abs(B) < 1.0e-10) f_h = f_m f_m = (1.0 - alpha_m * f_m)**0.25 f_h = sqrt(1.0 - alpha_h_fix * f_h) psi_m = 2.0 * (atan(f_m_lim) - atan(f_m)) + alog((f_m_lim - 1.0) * (f_m + 1.0)/((f_m_lim + 1.0) * (f_m - 1.0))) psi_h = alog((f_h_lim - 1.0) * (f_h + 1.0)/((f_h_lim + 1.0) * (f_h - 1.0))) / Pr_t_0_inv ! --- bulk Richardson number Rib_lim = zeta_lim * psi_h / (psi_m * psi_m) end subroutine ! -------------------------------------------------------------------------------- end module sfx_esm