#include "../includeF/sfx_def.fi" module sfx_sheba_noniterative ! modules used ! -------------------------------------------------------------------------------- #ifdef SFX_CHECK_NAN use sfx_common #endif use sfx_data use sfx_surface use sfx_sheba_noit_param #if defined(INCLUDE_CXX) use iso_c_binding, only: C_LOC, C_PTR, C_INT, C_FLOAT use C_FUNC #endif ! -------------------------------------------------------------------------------- ! directives list ! -------------------------------------------------------------------------------- implicit none private ! -------------------------------------------------------------------------------- ! public interface ! -------------------------------------------------------------------------------- public :: get_surface_fluxes public :: get_surface_fluxes_vec public :: get_psi_stable integer z0m_id integer z0t_id ! -------------------------------------------------------------------------------- ! -------------------------------------------------------------------------------- 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 ! -------------------------------------------------------------------------------- #if defined(INCLUDE_CXX) type, BIND(C), public :: sfx_sheba_noit_param_C real(C_FLOAT) :: kappa real(C_FLOAT) :: Pr_t_0_inv real(C_FLOAT) :: Pr_t_inf_inv real(C_FLOAT) :: alpha_m real(C_FLOAT) :: alpha_h real(C_FLOAT) :: alpha_h_fix real(C_FLOAT) :: Rib_max real(C_FLOAT) :: gamma real(C_FLOAT) :: zeta_a real(C_FLOAT) :: a_m real(C_FLOAT) :: b_m real(C_FLOAT) :: a_h real(C_FLOAT) :: b_h real(C_FLOAT) :: c_h end type type, BIND(C), public :: sfx_sheba_noit_numericsType_C integer(C_INT) :: maxiters_charnock integer(C_INT) :: maxiters_convection end type INTERFACE SUBROUTINE c_sheba_noit_compute_flux(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C, & name="c_sheba_noit_compute_flux") use sfx_data use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR Import :: sfx_sheba_noit_param_C, sfx_sheba_noit_numericsType_C implicit none integer(C_INT) :: grid_size type(C_PTR), value :: sfx type(C_PTR), value :: meteo type(sfx_sheba_noit_param_C) :: model_param type(sfx_surface_param) :: surface_param type(sfx_sheba_noit_numericsType_C) :: numerics type(sfx_phys_constants) :: constants END SUBROUTINE c_sheba_noit_compute_flux END INTERFACE #endif contains ! -------------------------------------------------------------------------------- #if defined(INCLUDE_CXX) subroutine set_c_struct_sfx_sheba_noit_param_values(sfx_model_param) type (sfx_sheba_noit_param_C), 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%Rib_max = Rib_max sfx_model_param%gamma = gamma sfx_model_param%zeta_a = zeta_a sfx_model_param%a_m = a_m sfx_model_param%b_m = b_m sfx_model_param%a_h = a_h sfx_model_param%b_h = b_h sfx_model_param%c_h = c_h end subroutine set_c_struct_sfx_sheba_noit_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_sheba_noit_param_C) :: model_param type (sfx_surface_param) :: surface_param type (sfx_sheba_noit_numericsType_C) :: 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_sheba_noit_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 c_sheba_noit_compute_flux(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n) #else 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), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i)) !write(*,*) 'get_flux', meteo%surface_type(i) call get_surface_fluxes(sfx_cell, meteo_cell, numerics) call push_sfx_data(sfx, sfx_cell, i) 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) real :: depth real :: lai integer :: surface_type ! ---------------------------------------------------------------------------- ! --- 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 Rib, Ri_sn !< bulk Richardson number, snow Richardson number real zeta !< = z/L [n/d] 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] 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 psi0_m, psi0_h !< universal functions (momentum) & (heat) [n/d] 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] 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 ! ---------------------------------------------------------------------------- #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_m1 = meteo%z0_m depth = meteo%depth lai = meteo%lai surface_type=meteo%surface_type !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) Re = u_dyn0 * z0_m / nu_air 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 ! --- restrict bulk Ri value ! *: note that this value is written to output ! Rib = min(Rib, Rib_max) if (surface_type == surface_snow) then write(*,*) 'sfx_snow1', Ri_sn, Rib, U, deltaS !write(*,*) 'sfx_snow2', deltaS, u_dyn0, S_mean 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) 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)) 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) if(Rib > 0)then 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 else Cm = kappa / psi_m Ct = kappa / psi_h end if ! --- 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 ! -------------------------------------------------------------------------------- !--------------------------------------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 real,intent(in) :: Rib, h, z0_m, z0_t real :: Ribl, C1, A1, A2, lne, lnet real :: psi_m, psi_h, psi0_m, psi0_h Ribl = (Rib * Pr_t_0_inv) * (1 - z0_t / h) / ((1 - z0_m / h)**2) call get_psi_stable(psi_m, psi_h, zeta_a, zeta_a) call get_psi_stable(psi0_m, psi0_h, zeta_a * z0_m / h, zeta_a * z0_t / h) lne = log(h/z0_m) lnet = log(h/z0_t) C1 = (lne**2)/lnet A1 = ((lne - psi_m + psi0_m)**(2*(gamma-1))) & & / ((zeta_a**(gamma-1))*((lnet-(psi_h-psi0_h)*Pr_t_0_inv)**(gamma-1))) A2 = ((lne - psi_m + psi0_m)**2) / (lnet-(psi_h-psi0_h)*Pr_t_0_inv) - C1 zeta = C1 * Ribl + A1 * A2 * (Ribl**gamma) end subroutine get_zeta ! 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_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 real :: q_m, q_h ! ---------------------------------------------------------------------------- q_m = ((1.0 - b_m) / b_m)**(1.0 / 3.0) x_m = (1.0 + zeta_m)**(1.0 / 3.0) psi_m = -3.0 * (a_m / b_m) * (x_m - 1.0) + 0.5 * (a_m / b_m) * q_m * (& 2.0 * log((x_m + q_m) / (1.0 + q_m)) - & log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + & 2.0 * sqrt(3.0) * (& atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - & atan((2.0 - q_m) / (sqrt(3.0) * q_m)))) q_h = sqrt(c_h * c_h - 4.0) x_h = zeta_h psi_h = -0.5 * b_h * log(1.0 + c_h * x_h + x_h * x_h) + & ((-a_h / q_h) + ((b_h * c_h) / (2.0 * q_h))) * (& log((2.0 * x_h + c_h - q_h) / (2.0 * x_h + c_h + q_h)) - & log((c_h - q_h) / (c_h + q_h))) end subroutine get_psi_stable 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_sheba_noniterative