#include <algorithm>
#include <cmath>
#include <limits>
#include <AMReX_GpuQualifiers.H>
#include <AMReX_REAL.H>
#include <ERF_Constants.H>
#include <ERF_MicrophysicsUtils.H>

Include dependency graph for ERF_KesslerUtils.H:

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

Go to the source code of this file.

Classes
struct	KesslerSourceTerms

struct	KesslerSaturationAdjustment

struct	KesslerFaceState

Functions
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int	kessler_face_donor_k (const int k, const int k_hi) noexcept

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	kessler_terminal_velocity (const amrex::Real rho, const amrex::Real qp) noexcept

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	kessler_precip_flux (const amrex::Real rho, const amrex::Real terminal_velocity, const amrex::Real qp) noexcept

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int	kessler_num_sedimentation_substeps (const amrex::Real current_reduced_value, const amrex::Real dt, const amrex::Real dzmin) noexcept

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	kessler_rain_accumulation_increment (const amrex::Real precip_mass_per_area) noexcept

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool	kessler_is_small_sedimentation_value (const amrex::Real value) noexcept

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE KesslerSaturationAdjustment	kessler_saturation_adjustment (const amrex::Real qv, const amrex::Real qc, const amrex::Real qsat, const amrex::Real dtqsat, const bool do_cond, const amrex::Real latent_over_cp) noexcept

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE KesslerFaceState	kessler_face_state (const int k, const int k_hi, const amrex::Real rho_km1, const amrex::Real rho_k, const amrex::Real qp_km1, const amrex::Real qp_k) noexcept

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	kessler_sedimentation_tendency (const amrex::Real fz_hi, const amrex::Real fz_lo, const amrex::Real rho, const amrex::Real dJinv, const amrex::Real coef) noexcept

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE KesslerSourceTerms	kessler_warm_rain_sources (const amrex::Real qv, const amrex::Real qc, const amrex::Real qp, const amrex::Real rho, const amrex::Real pressure_current_units, const amrex::Real qsat, const amrex::Real dtqsat, const amrex::Real dt, const bool do_cond, const amrex::Real latent_over_cp) noexcept

Variables
static constexpr amrex::Real	kessler_sedimentation_cfl_max = myhalf

static constexpr amrex::Real	kessler_water_depth_m_per_kg_m2 = amrex::Real(1.0) / rhor

static constexpr amrex::Real	kessler_millimeters_per_meter = amrex::Real(1000.0)

Function Documentation

◆ kessler_face_donor_k()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int kessler_face_donor_k	(	const int	k,
		const int	k_hi
	)

noexcept

 {
     return (k == k_hi + 1) ? k - 1 : k;
 }

Referenced by Kessler::AdvanceKessler().

Here is the caller graph for this function:

◆ kessler_face_state()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE KesslerFaceState kessler_face_state	(	const int	k,
		const int	k_hi,
		const amrex::Real	rho_km1,
		const amrex::Real	rho_k,
		const amrex::Real	qp_km1,
		const amrex::Real	qp_k
	)

noexcept

 {
     KesslerFaceState face_state{amrex::Real(0), amrex::Real(0)};
  
     // Sedimentation is one-way downward transport. Use a face-centered density
     // for the z-face flux, but upwind the transported rain mixing ratio from the
     // donor cell. A centered qp can create flux from downstream rain when the
     // donor cell is dry. The top boundary face uses the top cell below the
     // boundary.
     if (k == k_hi + 1) {
         face_state.rho = rho_km1;
         face_state.qp = qp_km1;
     } else if (k == 0) {
         face_state.rho = rho_k;
         face_state.qp = qp_k;
     } else {
         face_state.rho = amrex::Real(0.5) * (rho_km1 + rho_k);
         face_state.qp = qp_k;
     }
  
     // MUST MATCH: current AdvanceKessler nonnegative donor-qp clip before flux formation.
     face_state.qp = std::max(amrex::Real(0), face_state.qp);
     return face_state;
 }

Referenced by Kessler::AdvanceKessler().

Here is the caller graph for this function:

◆ kessler_is_small_sedimentation_value()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool kessler_is_small_sedimentation_value ( const amrex::Real value )

noexcept

 {
     // MUST MATCH: current AdvanceKessler sedimentation zero threshold literal.
     return std::fabs(value) < 1e-14;
 }

Referenced by Kessler::AdvanceKessler().

Here is the caller graph for this function:

◆ kessler_num_sedimentation_substeps()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int kessler_num_sedimentation_substeps	(	const amrex::Real	current_reduced_value,
		const amrex::Real	dt,
		const amrex::Real	dzmin
	)

noexcept

 {
     // MUST MATCH: current AdvanceKessler CFL_MAX precipitation substep formula.
     return static_cast<int>(std::ceil((current_reduced_value + std::numeric_limits<amrex::Real>::epsilon())
                                       * (dt / dzmin) / kessler_sedimentation_cfl_max));
 }

Referenced by Kessler::AdvanceKessler().

Here is the caller graph for this function:

◆ kessler_precip_flux()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real kessler_precip_flux	(	const amrex::Real	rho,
		const amrex::Real	terminal_velocity,
		const amrex::Real	qp
	)

noexcept

 {
     return rho * terminal_velocity * qp;
 }

Referenced by Kessler::AdvanceKessler().

Here is the caller graph for this function:

◆ kessler_rain_accumulation_increment()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real kessler_rain_accumulation_increment ( const amrex::Real precip_mass_per_area )

noexcept

 {
     // Convert precipitation mass per unit area [kg m^-2] to liquid-water depth
     // [mm]: divide by water density [kg m^-3], then multiply by mm per meter.
     return precip_mass_per_area
         * kessler_water_depth_m_per_kg_m2
         * kessler_millimeters_per_meter;
 }

Referenced by Kessler::AdvanceKessler().

Here is the caller graph for this function:

◆ kessler_saturation_adjustment()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE KesslerSaturationAdjustment kessler_saturation_adjustment	(	const amrex::Real	qv,
		const amrex::Real	qc,
		const amrex::Real	qsat,
		const amrex::Real	dtqsat,
		const bool	do_cond,
		const amrex::Real	latent_over_cp
	)

noexcept

 {
     KesslerSaturationAdjustment source_terms{amrex::Real(0), amrex::Real(0)};
  
     const amrex::Real fac = latent_over_cp * dtqsat;
  
     // MUST MATCH: current AdvanceKessler condensation branch behavior.
     if ((qv > qsat) && do_cond) {
         source_terms.dq_vapor_to_cloud = std::min(qv, (qv - qsat) / (amrex::Real(1) + fac));
     }
  
     // MUST MATCH: current AdvanceKessler cloud-evaporation branch behavior.
     if ((qv < qsat) && (qc > amrex::Real(0)) && do_cond) {
         source_terms.dq_cloud_to_vapor = std::min(qc, (qsat - qv) / (amrex::Real(1) + fac));
     }
  
     return source_terms;
 }

Referenced by Kessler::AdvanceKessler(), and kessler_warm_rain_sources().

Here is the caller graph for this function:

◆ kessler_sedimentation_tendency()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real kessler_sedimentation_tendency	(	const amrex::Real	fz_hi,
		const amrex::Real	fz_lo,
		const amrex::Real	rho,
		const amrex::Real	dJinv,
		const amrex::Real	coef
	)

noexcept

 {
     // MUST MATCH: current AdvanceKessler sedimentation tendency arithmetic order.
     return dJinv * (amrex::Real(1) / rho) * (fz_hi - fz_lo) * coef;
 }

Referenced by Kessler::AdvanceKessler().

Here is the caller graph for this function:

◆ kessler_terminal_velocity()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real kessler_terminal_velocity	(	const amrex::Real	rho,
		const amrex::Real	qp
	)

noexcept

 {
     // MUST MATCH: current AdvanceKessler terminal-velocity coefficients.
     return amrex::Real(36.34)
         * std::pow(rho * amrex::Real(0.001) * qp, amrex::Real(0.1346))
         * std::pow(rho / amrex::Real(1.16), -myhalf);
 }

Referenced by Kessler::AdvanceKessler().

Here is the caller graph for this function:

◆ kessler_warm_rain_sources()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE KesslerSourceTerms kessler_warm_rain_sources	(	const amrex::Real	qv,
		const amrex::Real	qc,
		const amrex::Real	qp,
		const amrex::Real	rho,
		const amrex::Real	pressure_current_units,
		const amrex::Real	qsat,
		const amrex::Real	dtqsat,
		const amrex::Real	dt,
		const bool	do_cond,
		const amrex::Real	latent_over_cp
	)

noexcept

 {
     // Units must match current AdvanceKessler: q* are mass fractions, rho is in
     // the production density units currently carried by Kessler, pressure is in
     // Kessler's mbar / hPa path for qsat helpers, and dt is in seconds.
     const KesslerSaturationAdjustment saturation_adjustment =
         kessler_saturation_adjustment(qv, qc, qsat, dtqsat, do_cond, latent_over_cp);
     KesslerSourceTerms source_terms{
         saturation_adjustment.dq_vapor_to_cloud,
         saturation_adjustment.dq_cloud_to_vapor,
         amrex::Real(0),
         amrex::Real(0)};
  
     const amrex::Real capacity = amrex::max(
         amrex::Real(0), (qsat - qv) / (amrex::Real(1) + latent_over_cp * dtqsat));
     const amrex::Real remaining_capacity = amrex::max(
         amrex::Real(0), capacity - saturation_adjustment.dq_cloud_to_vapor);
     // Cloud evaporation gets first claim on the linearized subsaturation
     // capacity. Rain evaporation can use only the remaining capacity.
     const amrex::Real qv_eff = qv
         - saturation_adjustment.dq_vapor_to_cloud
         + saturation_adjustment.dq_cloud_to_vapor;
  
     if ((qp > amrex::Real(0)) && (qsat > amrex::Real(0))
         && (pressure_current_units > amrex::Real(0))
         && (qv_eff < qsat) && (remaining_capacity > amrex::Real(0))) {
         // MUST MATCH: current AdvanceKessler rain-evaporation coefficients and exponents.
         const amrex::Real coeff = amrex::Real(1.6)
             + amrex::Real(124.9) * std::pow(amrex::Real(0.001) * rho * qp, amrex::Real(0.2046));
         const amrex::Real raw_rain_to_vapor = amrex::Real(1) / (amrex::Real(0.001) * rho)
             * (amrex::Real(1) - qv_eff / qsat)
             * coeff
             * std::pow(amrex::Real(0.001) * rho * qp, amrex::Real(0.525))
             / (amrex::Real(5.4e5) + amrex::Real(2.55e6) / (pressure_current_units * qsat))
             * dt;
         source_terms.dq_rain_to_vapor = amrex::min(qp, amrex::min(raw_rain_to_vapor, remaining_capacity));
     }
  
     if (qc > amrex::Real(0)) {
         const amrex::Real qcc = qc;
         // Cloud-to-rain conversion is capped by cloud water available after
         // saturation adjustment, not by the original qc.
         const amrex::Real available_cloud = amrex::max(
             amrex::Real(0), qc + saturation_adjustment.dq_vapor_to_cloud - saturation_adjustment.dq_cloud_to_vapor);
         amrex::Real auto_r = amrex::Real(0);
         if (qcc > qcw0) {
             // MUST MATCH: current AdvanceKessler autoconversion trigger.
             auto_r = alphaelq;
         }
  
         // MUST MATCH: current AdvanceKessler accretion coefficient and exponent.
         amrex::Real accrr = amrex::Real(0);
         accrr = amrex::Real(2.2) * std::pow(qp, amrex::Real(0.875));
         source_terms.dq_cloud_to_rain = dt * (accrr * qcc + auto_r * (qcc - qcw0));
         source_terms.dq_cloud_to_rain = amrex::min(source_terms.dq_cloud_to_rain, available_cloud);
     }
  
     return source_terms;
 }