ERF_MOSTUtils.H File Reference

#include <AMReX_Algorithm.H>
#include <AMReX_BLassert.H>
#include <AMReX_Math.H>
#include "ERF_Constants.H"

Include dependency graph for ERF_MOSTUtils.H:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	most_data
struct	similarity_funs
struct	most_flux_scale
struct	most_aitken

Functions
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	most_bracket_floor (amrex::Real C)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	most_bracket (amrex::Real C, amrex::Real psi)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	most_rib_resid (const similarity_funs &sfuns, amrex::Real C, amrex::Real Rib, amrex::Real zeta)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool	most_rib_zeta_is_root (const similarity_funs &sfuns, amrex::Real C, amrex::Real Rib, amrex::Real zeta, amrex::Real delta)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	most_rib_zeta_bracketed (const similarity_funs &sfuns, amrex::Real C, amrex::Real Rib, amrex::Real zeta0, amrex::Real tol)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	most_clamp_zeta (amrex::Real zeta, amrex::Real zeta_fold, bool &at_fold)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE most_flux_scale	most_fold_flux_scale (bool at_fold, amrex::Real zeta_fold, amrex::Real ustar, amrex::Real thv, amrex::Real tflux_s, amrex::Real tflux_q, bool spec_qflux, amrex::Real kappa, amrex::Real gravity, amrex::Real zref)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	most_flux_outputs (const most_data &mdata, bool spec_qflux, bool at_fold, amrex::Real zeta_fold, amrex::Real C, amrex::Real psi_h, amrex::Real ustar, amrex::Real zref, amrex::Real tm, amrex::Real tvm, amrex::Real qvm, amrex::Real qflux, amrex::Real &olen, amrex::Real &t_surf, amrex::Real &t_star, amrex::Real &q_surf, amrex::Real &q_star)
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real	air_viscosity (amrex::Real T_degK)
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real	Charnock_roughness (amrex::Real ustar, amrex::Real Cnk_a, amrex::Real nu, bool visc)
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real	Mod_Charnock_roughness (amrex::Real ustar, amrex::Real Cnk_b)
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real	COARE3_roughness (amrex::Real zref, amrex::Real umm, amrex::Real ustar, amrex::Real nu, bool visc)
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real	Donelan_roughness (amrex::Real ustar, amrex::Real nu, bool visc)
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real	WaveCoupled_roughness (amrex::Real Hwave, amrex::Real Lwave, amrex::Real eps, amrex::Real ustar, amrex::Real nu, bool visc)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	bound_inverted_roughness (amrex::Real z0, amrex::Real zref) noexcept
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real	Compute_roughness (amrex::Real zref, amrex::Real Olen, amrex::Real umm, amrex::Real ustar)

Variables
constexpr amrex::Real	most_zeta_lim = amrex::Real(100.0)
constexpr amrex::Real	most_min_ustar = amrex::Real(0.001)
constexpr amrex::Real	inverted_roughness_lo = amrex::Real(1.0e-7)

Function Documentation

air_viscosity()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real air_viscosity

(

amrex::Real

T_degK

)

Empirical kinematic viscosity [m2/s] formula from Andreas (1989) CRREL Rep. 89-11, valid between -173 and 277 deg C.

Parameters

[in]

T_degK

air temperature in Kelvin

{
    amrex::Real TC = T_degK - amrex::Real(273.15);
    return amrex::Real(1.326e-5)*(one + amrex::Real(6.542e-3)*TC + amrex::Real(8.301e-6)*TC*TC - amrex::Real(4.84e-9)*TC*TC*TC);
}

Referenced by adiabatic_charnock::iterate_flux(), adiabatic_donelan::iterate_flux(), adiabatic_wave_coupled::iterate_flux(), surface_flux_charnock::iterate_flux(), surface_flux_donelan::iterate_flux(), surface_temp_charnock::iterate_flux(), surface_temp_donelan::iterate_flux(), surface_flux_wave_coupled::iterate_flux(), and surface_temp_wave_coupled::iterate_flux().

Here is the caller graph for this function:

bound_inverted_roughness()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real bound_inverted_roughness ( amrex::Real  z0, amrex::Real  zref  )

noexcept

Bounds on the roughness recovered by Compute_roughness. The inversion is unbounded above: in weak wind under strong stability (u* from an LSM or bulk flux large relative to kappa*U, psi_m << 0) it returns z0 >> zref. The upper bound is above every tabulated land-cover roughness in NoahmpTable.TBL (Z0MVT <= 1.10 m) and below the 2 m diagnostic height; z0 is also kept below half of zref so that ln(zref/z0) stays positive.

{
    const amrex::Real z0_lo = inverted_roughness_lo;
    const amrex::Real z0_hi = amrex::min(amrex::Real(1.5), amrex::Real(0.5) * zref);
    return amrex::min(amrex::max(z0, z0_lo), z0_hi);
}

Referenced by Compute_roughness().

Here is the caller graph for this function:

Charnock_roughness()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Charnock_roughness	(	amrex::Real	ustar,
		amrex::Real	Cnk_a,
		amrex::Real	nu,
		bool	visc
	)

Charnock roughness with optional aerodynamically-smooth contribution

Smooth-flow limit follows z0 = 0.11 nu / u* (Smith 1988; ECMWF; WRF)

Parameters

[in]	ustar	friction velocity
[in]	cnk_a	charnock constant
[in]	nu	intrinsic fluid viscosity
[in]	visc	flag to apply smooth flow limit

{
    amrex::Real z0 = (Cnk_a / CONST_GRAV) * ustar * ustar;
    if (visc) {
        z0 += amrex::Real(0.11) * nu / std::max(ustar, amrex::Real(0.01));
    }
    // From Davis et al 2008 MWR
    return std::min(std::max(z0, amrex::Real(1.27e-7)), amrex::Real(2.85e-3));
}

Referenced by adiabatic_charnock::iterate_flux(), surface_flux_charnock::iterate_flux(), and surface_temp_charnock::iterate_flux().

Here is the caller graph for this function:

COARE3_roughness()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real COARE3_roughness	(	amrex::Real	zref,
		amrex::Real	umm,
		amrex::Real	ustar,
		amrex::Real	nu,
		bool	visc
	)

COARE3.0, following WRF phys/module_sf_mynn.F

Smooth-flow limit follows z0 = 0.11 nu / u* (Smith 1988; ECMWF; WRF)

Parameters

[in]	zref	reference height for the wind speed
[in]	umm	horizontal wind speed magnitude at zref
[in]	ustar	friction velocity
[in]	nu	intrinsic fluid viscosity
[in]	visc	flag to apply smooth flow limit

{
    amrex::Real ws10m = umm*std::log(amrex::Real(10.)/amrex::Real(1e-4)) / std::log(zref/amrex::Real(1e-4));
    amrex::Real Czc = amrex::Real(0.011) + amrex::Real(0.007)*std::min(std::max((ws10m-amrex::Real(10.))/amrex::Real(8.), amrex::Real(0.0)), amrex::Real(1.0));
    amrex::Real z0 = (Czc / CONST_GRAV) * ustar * ustar;
    if (visc) {
        z0 += amrex::Real(0.11) * nu / std::max(ustar, amrex::Real(0.01));
    }
    // From Davis et al 2008 MWR
    return std::min(std::max(z0, amrex::Real(1.27e-7)), amrex::Real(2.85e-3));
}

Referenced by adiabatic_charnock::iterate_flux(), surface_flux_charnock::iterate_flux(), and surface_temp_charnock::iterate_flux().

Here is the caller graph for this function:

Compute_roughness()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_roughness	(	amrex::Real	zref,
		amrex::Real	Olen,
		amrex::Real	umm,
		amrex::Real	ustar
	)

Invert log-law to get roughness

Parameters

[in]	zref	reference height for the wind speed
[in]	Olen	Obukhov length
[in]	umm	horizontal wind speed magnitude at zref
[in]	ustar	friction velocity

{
    similarity_funs sfuns;
    amrex::Real zeta  = zref / Olen;
    amrex::Real psi_m = sfuns.calc_psi_m2(zeta);
    // z0 = zref / exp(arg) is below the lower bound once arg > log(zref / z0_lo),
    // so the bound is returned without evaluating exp(arg), which overflows when
    // u* sits at its floor of machine epsilon (zero stress: bulk_coeff with
    // Cd = 0, or an LSM that returns none). At or below the threshold exp(arg)
    // is finite and the inversion runs as before, so no other input changes.
    const amrex::Real arg = KAPPA*umm/ustar + psi_m;
    if (arg > std::log(zref / inverted_roughness_lo)) {
        return bound_inverted_roughness(inverted_roughness_lo, zref);
    }
    amrex::Real z0    = zref / std::exp(arg);
    return bound_inverted_roughness(z0, zref);
}

Referenced by SurfaceLayer::compute_sfc_params_from_lsm_fluxes(), and SurfaceLayer::compute_SurfaceLayer_bcs().

Here is the call graph for this function:

Here is the caller graph for this function:

Donelan_roughness()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Donelan_roughness	(	amrex::Real	ustar,
		amrex::Real	nu,
		bool	visc
	)

Updated Donelan 2004 as found in AHW

Matches implementations in WRF phys/physics_mmm/sf_sfclayrev.F90 and phys/module_sf_mynn.F

Smooth flow limit is intrinsic in zn2.

Parameters

[in]	ustar	friction velocity
[in]	nu	intrinsic fluid viscosity
[in]	visc	flag to apply smooth flow limit

{
    constexpr amrex::Real ozo{amrex::Real(1.59e-5)};
    amrex::Real zw  = std::min(std::pow(ustar/amrex::Real(1.06), amrex::Real(0.3)),amrex::Real(1.0));
    amrex::Real zn1 = amrex::Real(0.011)*ustar*ustar/CONST_GRAV + ozo;
    amrex::Real zn2 = amrex::Real(10.)*std::exp(-amrex::Real(9.5)*std::pow(ustar, -amrex::Real(.3333)));
    if (visc) {
        zn2 += amrex::Real(0.11) * nu / std::max(ustar, amrex::Real(0.01));
    }
    amrex::Real z0 = (amrex::Real(1.0)-zw) * zn1 + zw * zn2;
    // From Davis et al 2008 MWR
    return std::min(std::max(z0, amrex::Real(1.27e-7)), amrex::Real(2.85e-3));
}

Referenced by adiabatic_donelan::iterate_flux(), surface_flux_donelan::iterate_flux(), and surface_temp_donelan::iterate_flux().

Here is the caller graph for this function:

Mod_Charnock_roughness()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Mod_Charnock_roughness	(	amrex::Real	ustar,
		amrex::Real	Cnk_b
	)

Modified Charnock roughness

Smooth-flow limit is intrinsic.

Parameters

[in]	ustar	friction velocity
[in]	cnk_b	modified charnock constant

{
    amrex::Real z0 = std::exp( (amrex::Real(2.7)*ustar - amrex::Real(1.8)/Cnk_b) / (ustar + amrex::Real(0.17)/Cnk_b) );
    return std::min(std::max(z0, amrex::Real(1.27e-7)), amrex::Real(1.0e-2));
}

Referenced by adiabatic_mod_charnock::iterate_flux(), surface_flux_mod_charnock::iterate_flux(), and surface_temp_mod_charnock::iterate_flux().

Here is the caller graph for this function:

most_bracket()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real most_bracket	(	amrex::Real	C,
		amrex::Real	psi
	)

Log-law bracket ln(zref/z0) - psi of the MOST flux kernels, held at or above most_bracket_floor(C).

In strongly unstable air psi exceeds C = ln(zref/z0) and the bracket changes sign, which turns u* negative and puts the diagnosed surface temperature on the wrong side of the air for a prescribed flux. With the floor the bracket is at least 1/2, so u* stays positive and finite for every C.

Parameters

[in]	C	ln(zref/z0)
[in]	psi	stability correction (momentum or heat)

{
    return amrex::max(C - psi, most_bracket_floor(C));
}

Referenced by adiabatic::iterate_flux(), adiabatic_eb::iterate_flux(), adiabatic_mod_charnock::iterate_flux(), adiabatic_charnock::iterate_flux(), adiabatic_donelan::iterate_flux(), adiabatic_wave_coupled::iterate_flux(), surface_flux_charnock::iterate_flux(), surface_flux_mod_charnock::iterate_flux(), surface_flux_donelan::iterate_flux(), surface_temp_charnock::iterate_flux(), surface_temp_mod_charnock::iterate_flux(), surface_temp_donelan::iterate_flux(), surface_flux_wave_coupled::iterate_flux(), surface_temp_wave_coupled::iterate_flux(), surface_flux::iterate_flux(), surface_temp::iterate_flux(), surface_temp_eb::iterate_flux(), surface_flux_eb::iterate_flux(), most_flux_outputs(), and most_rib_resid().

Here is the call graph for this function:

Here is the caller graph for this function:

most_bracket_floor()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real most_bracket_floor

(

amrex::Real

C

)

Lower bound of the log-law bracket ln(zref/z0) - psi of the MOST flux kernels: C/2 held between 1/2 and 1.

The bound is 1 for zref >= e^2 z0 and C/2 down to zref = e z0, so it does not exceed the neutral bracket C wherever zref >= e^(1/2) z0. The lower limit of 1/2 keeps the bracket, and so u*, positive and finite when z0 approaches or exceeds zref (C <= 0).

The bound is set by the under-relaxed u* iteration, not by the physics. At an unstable fixed point the relaxed map has slope (1 - alpha) - alpha s / bracket, with s = -3 (z/L) dpsi_m/d(z/L) between 1.3 (z/L = -0.5) and 2.4 (z/L = -50) for Businger-Dyer. With alpha = 0.5 the iteration converges only for bracket > s/3 (0.4 to 0.8), so a bound of 1 keeps the switch between the bounded and free branches inside the convergent range, while 0.1 C (0.3 at C = 3) does not and the iteration oscillates between two values. If alpha changes, the bound needs revisiting; see PR #4866.

Parameters

[in]

C

ln(zref/z0)

{
    return amrex::max(amrex::min(amrex::Real(0.5) * C, amrex::Real(1.0)), amrex::Real(0.5));
}

Referenced by similarity_funs::calc_zeta_fold(), and most_bracket().

Here is the caller graph for this function:

most_clamp_zeta()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real most_clamp_zeta	(	amrex::Real	zeta,
		amrex::Real	zeta_fold,
		bool &	at_fold
	)

z/L of the prescribed-flux kernels, bounded below by -most_zeta_lim and above by the fold (see similarity_funs::calc_zeta_fold), where the prescribed stable flux stops having a solution.

Parameters

[in]	zeta	z/L from the current u* and flux
[in]	zeta_fold	z/L at the fold
[out]	at_fold	whether the fold bound was applied

{
    at_fold = (zeta >= zeta_fold);
    return amrex::max(amrex::min(zeta, zeta_fold), -most_zeta_lim);
}

Referenced by surface_flux_charnock::iterate_flux(), surface_flux_mod_charnock::iterate_flux(), surface_flux_donelan::iterate_flux(), surface_flux_wave_coupled::iterate_flux(), surface_flux::iterate_flux(), and surface_flux_eb::iterate_flux().

Here is the caller graph for this function:

most_flux_outputs()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void most_flux_outputs	(	const most_data &	mdata,
		bool	spec_qflux,
		bool	at_fold,
		amrex::Real	zeta_fold,
		amrex::Real	C,
		amrex::Real	psi_h,
		amrex::Real	ustar,
		amrex::Real	zref,
		amrex::Real	tm,
		amrex::Real	tvm,
		amrex::Real	qvm,
		amrex::Real	qflux,
		amrex::Real &	olen,
		amrex::Real &	t_surf,
		amrex::Real &	t_star,
		amrex::Real &	q_surf,
		amrex::Real &	q_star
	)

Surface state of a prescribed-flux kernel once its u* iteration has converged.

Sets t_surf, t_star and q_star (and q_surf when the moisture flux is prescribed). At the fold the prescribed fluxes are reduced by most_fold_flux_scale and the Obukhov length is set to the fold value. u* is floored at most_min_ustar here, so every kernel applies the same floor.

Parameters

[in]	mdata	MOST data (prescribed fluxes, kappa, gravity)
[in]	spec_qflux	whether the moisture flux is prescribed
[in]	at_fold	whether z/L was held at the fold on the last pass
[in]	zeta_fold	z/L at the fold
[in]	C	ln(zref/z0)
[in]	psi_h	heat stability correction of the last pass
[in]	ustar	friction velocity (floored at most_min_ustar)
[in]	zref	reference height
[in]	tm	potential temperature at zref
[in]	tvm	virtual potential temperature at zref
[in]	qvm	water vapor at zref
[in]	qflux	moisture flux of the last pass
[in,out]	olen	Obukhov length
[out]	t_surf	surface potential temperature
[out]	t_star	temperature scale
[in,out]	q_surf	surface water vapor (read when the moisture flux is diagnosed)
[out]	q_star	moisture scale

{
    ustar = amrex::max(ustar, amrex::Real(most_min_ustar));
    const most_flux_scale scale = most_fold_flux_scale(at_fold, zeta_fold, ustar, tvm,
                                                       mdata.surf_temp_flux * (one + epsv * qvm),
                                                       qflux * epsv * tm, spec_qflux,
                                                       mdata.kappa, mdata.gravity, zref);
    if (at_fold) { olen = zref / zeta_fold; }
    const amrex::Real CmPsih = most_bracket(C, psi_h);
    const amrex::Real hflux  = scale.temp * mdata.surf_temp_flux;
    t_surf = hflux * CmPsih / (ustar * mdata.kappa) + tm;
    t_star = -hflux / ustar;
    if (spec_qflux) {
        const amrex::Real mflux = scale.moist * mdata.surf_moist_flux;
        q_surf = mflux * CmPsih / (ustar * mdata.kappa) + qvm;
        q_star = -mflux / ustar;
    } else {
        q_star = mdata.kappa * (qvm - q_surf) / CmPsih;
    }
}

Referenced by surface_flux_charnock::iterate_flux(), surface_flux_mod_charnock::iterate_flux(), surface_flux_donelan::iterate_flux(), surface_flux_wave_coupled::iterate_flux(), surface_flux::iterate_flux(), and surface_flux_eb::iterate_flux().

Here is the call graph for this function:

Here is the caller graph for this function:

most_fold_flux_scale()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE most_flux_scale most_fold_flux_scale	(	bool	at_fold,
		amrex::Real	zeta_fold,
		amrex::Real	ustar,
		amrex::Real	thv,
		amrex::Real	tflux_s,
		amrex::Real	tflux_q,
		bool	spec_qflux,
		amrex::Real	kappa,
		amrex::Real	gravity,
		amrex::Real	zref
	)

Part of a prescribed stable flux the wind can carry at the fold.

The virtual heat flux is tflux_s + tflux_q, with tflux_s = surf_temp_flux (1 + epsv qv) and tflux_q = qflux epsv theta. When z/L is held at the fold (see similarity_funs::calc_zeta_fold), the wind is below the minimum the prescribed flux needs and carries only -u*^3 thv zeta_fold / (kappa g zref). Only the prescribed components that are downward drive the fold, so only they are reduced, by a common factor, until the net virtual flux equals that value. An upward component, or a moisture flux diagnosed from q_surf, is left unchanged; if no prescribed component is downward, nothing is reduced. The factor is not taken below 0, so when a diagnosed condensation flux alone is stronger than the fold value, the downward prescribed components are removed and the net flux remains stronger than the fold value. Returns 1 for both factors away from the fold.

Parameters

[in]	at_fold	whether z/L was held at the fold
[in]	zeta_fold	z/L at the fold
[in]	ustar	friction velocity
[in]	thv	virtual potential temperature at zref
[in]	tflux_s	sensible part of the virtual heat flux
[in]	tflux_q	moisture part of the virtual heat flux
[in]	spec_qflux	whether the moisture flux is prescribed
[in]	kappa	von Karman constant
[in]	gravity	gravitational acceleration
[in]	zref	reference height

{
    most_flux_scale scale;
    const amrex::Real tflux = tflux_s + tflux_q;
    if (!at_fold || tflux >= zero) { return scale; }
    const amrex::Real tflux_fold = -ustar * ustar * ustar * thv * zeta_fold / (kappa * gravity * zref);
    if (tflux_fold <= tflux) { return scale; }
    const bool down_s = (tflux_s < zero);
    const bool down_q = (spec_qflux && tflux_q < zero);
    const amrex::Real down = (down_s ? tflux_s : zero) + (down_q ? tflux_q : zero);
    if (down >= zero) { return scale; }
    const amrex::Real r = amrex::max((down + tflux_fold - tflux) / down, amrex::Real(0.0));
    if (down_s) { scale.temp  = r; }
    if (down_q) { scale.moist = r; }
    return scale;
}

Referenced by most_flux_outputs().

Here is the caller graph for this function:

most_rib_resid()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real most_rib_resid	(	const similarity_funs &	sfuns,
		amrex::Real	C,
		amrex::Real	Rib,
		amrex::Real	zeta
	)

Residual zeta - Rib (C - psi_m)^2 / (C - psi_h) of the prescribed-temperature kernels' z/L equation, with the Jimenez functions and the brackets of most_bracket.

Parameters

[in]	sfuns	similarity functions
[in]	C	ln(zref/z0)
[in]	Rib	bulk Richardson number
[in]	zeta	z/L

{
    const amrex::Real CmPsim = most_bracket(C, sfuns.calc_psi_m2(zeta));
    const amrex::Real CmPsih = most_bracket(C, sfuns.calc_psi_h2(zeta));
    return zeta - Rib * CmPsim * CmPsim / CmPsih;
}

Referenced by most_rib_zeta_bracketed(), and most_rib_zeta_is_root().

Here is the call graph for this function:

Here is the caller graph for this function:

most_rib_zeta_bracketed()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real most_rib_zeta_bracketed	(	const similarity_funs &	sfuns,
		amrex::Real	C,
		amrex::Real	Rib,
		amrex::Real	zeta0,
		amrex::Real	tol
	)

Solve zeta = Rib (C - psi_m)^2 / (C - psi_h) of the prescribed-temperature kernels by bracketing and bisection, with the Jimenez functions and the brackets of most_bracket.

The relaxed fixed-point iteration stalls where zeta - Rib F(zeta) nearly touches zero without crossing it, which the Jimenez stable functions produce when z0 is comparable to or above zref (C below about 1): each pass then moves zeta by about the tolerance. From zeta0 this steps, with doubling steps, in the direction the iteration moves until the residual changes sign, and bisects that interval, so it returns the root the iteration was heading to, or the next one past a near-touch. The first step is limited so that two nearby roots are not stepped over together. The search stays within +-most_zeta_lim, and returns that bound when no root lies inside it, as the prescribed-flux kernels bound zeta.

Parameters

[in]	sfuns	similarity functions
[in]	C	ln(zref/z0)
[in]	Rib	bulk Richardson number
[in]	zeta0	z/L to start from
[in]	tol	width of the final interval

{
    auto resid = [&] (amrex::Real zeta) { return most_rib_resid(sfuns, C, Rib, zeta); };
 
    amrex::Real a  = zeta0;
    amrex::Real ra = resid(a);
    if (ra == zero) { return a; }
 
    // A non-finite residual (a NaN or infinite C, Rib or zeta0) cannot be bracketed: with it the
    // step below is NaN and the search never reaches zlim. Return it, so the bad input shows up
    // in the caller's results rather than as a hang.
    if (!amrex::Math::isfinite(ra)) { return ra; }
 
    // Step against the residual, within +-most_zeta_lim, until it changes sign
    const amrex::Real dir  = (ra > zero) ? amrex::Real(-1.0) : amrex::Real(1.0);
    const amrex::Real zlim = dir * most_zeta_lim;
    auto step = [&] (amrex::Real z, amrex::Real dz) {
        return (dir > zero) ? amrex::min(z + dz, zlim) : amrex::max(z - dz, zlim);
    };
    amrex::Real h  = amrex::max(amrex::min(std::abs(ra),
                                           amrex::Real(0.1) * amrex::max(amrex::Real(one), std::abs(a))),
                                tol);
    amrex::Real b  = step(a, h);
    amrex::Real rb = resid(b);
    // With finite values the doubling step reaches zlim in well under 200 steps (from tol to
    // 2 most_zeta_lim); the cap only guarantees that this loop ends
    int nstep = 0;
    while ((rb > zero) == (ra > zero) && b != zlim && nstep < 200) {
        a  = b;
        ra = rb;
        h *= two;
        b  = step(a, h);
        rb = resid(b);
        ++nstep;
    }
    // No root within the bound: hold zeta there
    if ((rb > zero) == (ra > zero)) { return zlim; }
 
    // Bisect [a, b]
    int n = 0;
    while (std::abs(b - a) > tol && n < 60) {
        const amrex::Real m  = myhalf * (a + b);
        const amrex::Real rm = resid(m);
        if ((rm > zero) == (ra > zero)) {
            a  = m;
            ra = rm;
        } else {
            b = m;
        }
        ++n;
    }
    return myhalf * (a + b);
}

Referenced by surface_temp::iterate_flux(), and surface_temp_eb::iterate_flux().

Here is the call graph for this function:

Here is the caller graph for this function:

most_rib_zeta_is_root()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool most_rib_zeta_is_root	(	const similarity_funs &	sfuns,
		amrex::Real	C,
		amrex::Real	Rib,
		amrex::Real	zeta,
		amrex::Real	delta
	)

Whether a root of most_rib_resid lies within delta of zeta.

The relaxed fixed-point iteration stops on a small step, which a near-touch of the residual also gives where no root is near; the residual must change sign across [zeta - delta, zeta + delta] for zeta to be taken as a solution.

Parameters

[in]	sfuns	similarity functions
[in]	C	ln(zref/z0)
[in]	Rib	bulk Richardson number
[in]	zeta	z/L the relaxed iteration stopped at
[in]	delta	half-width of the interval checked

{
    const amrex::Real rlo = most_rib_resid(sfuns, C, Rib, zeta - delta);
    const amrex::Real rhi = most_rib_resid(sfuns, C, Rib, zeta + delta);
    return (rlo > zero) != (rhi > zero) || rlo == zero || rhi == zero;
}

Referenced by surface_temp::iterate_flux(), and surface_temp_eb::iterate_flux().

Here is the call graph for this function:

Here is the caller graph for this function:

WaveCoupled_roughness()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real WaveCoupled_roughness	(	amrex::Real	Hwave,
		amrex::Real	Lwave,
		amrex::Real	eps,
		amrex::Real	ustar,
		amrex::Real	nu,
		bool	visc
	)

Wave coupled roughness

Smooth-flow limit follows z0 = 0.11 nu / u* (Smith 1988; ECMWF; WRF)

Parameters

[in]	Hwave	wave height
[in]	Lwave	wave length
[in]	eps	avoid divide by zero
[in]	ustar	friction velocity
[in]	nu	intrinsic fluid viscosity
[in]	visc	flag to apply smooth flow limit

{
 
    amrex::Real z0 = amrex::Real(1200.0) * Hwave * std::pow( Hwave/(Lwave+eps), amrex::Real(4.5) );
    if (visc) {
        z0 += amrex::Real(0.11) * nu / std::max(ustar, amrex::Real(0.01));
    }
    return std::min(std::max(z0, amrex::Real(1.27e-7)), amrex::Real(1.0e-2));
}

Referenced by adiabatic_wave_coupled::iterate_flux(), surface_flux_wave_coupled::iterate_flux(), and surface_temp_wave_coupled::iterate_flux().

Here is the caller graph for this function:

Variable Documentation

inverted_roughness_lo

constexpr amrex::Real inverted_roughness_lo = amrex::Real(1.0e-7)

constexpr

Lower bound of the roughness recovered by Compute_roughness; shared with the guard in Compute_roughness so the two cannot drift apart.

Referenced by bound_inverted_roughness(), and Compute_roughness().

most_min_ustar

constexpr amrex::Real most_min_ustar = amrex::Real(0.001)

inlineconstexpr

Minimum friction velocity [m/s] of the prescribed-flux kernels' outputs (see most_flux_outputs).

Referenced by most_flux_outputs().

most_zeta_lim

constexpr amrex::Real most_zeta_lim = amrex::Real(100.0)

inlineconstexpr

Bound on |z/L| of the MOST kernels (see most_rib_zeta_bracketed and most_clamp_zeta).

Referenced by most_clamp_zeta(), and most_rib_zeta_bracketed().