use crate::{
error::{AnalysisError, Result},
parcel::{convective_parcel, lowest_level_parcel, mixed_layer_parcel, Parcel},
parcel_profile::find_parcel_start_data,
};
use itertools::izip;
use metfor::{Celsius, CelsiusDiff, JpKg, Kelvin, Meters, Quantity};
use sounding_base::Sounding;
#[derive(Debug)]
pub struct PlumePotentialAnal {
plume_vals: Vec<(Celsius, JpKg)>,
}
impl PlumePotentialAnal {
pub fn analyze(snd: &Sounding) -> Result<PlumePotentialAnal> {
const MAX_HEATING: CelsiusDiff = CelsiusDiff(30.0);
const DT: CelsiusDiff = CelsiusDiff(0.1);
let mut plume_vals = Vec::with_capacity(MAX_HEATING.unpack() as usize);
for pcl in plume_parcels(snd, MAX_HEATING, DT)? {
if let Ok((cape, _)) = lift_parcel(pcl, snd) {
plume_vals.push((pcl.temperature, cape));
} else {
break;
}
}
Ok(PlumePotentialAnal { plume_vals })
}
pub fn analysis_values(&self) -> &[(Celsius, JpKg)] {
&self.plume_vals
}
}
pub fn convective_parcel_initiation_energetics(
snd: &Sounding,
) -> Result<(Celsius, CelsiusDiff, JpKg, JpKg)> {
let starting_parcel = convective_parcel(snd)?;
let mut no_cloud_pcl = starting_parcel;
let mut no_cloud_pcl_data = lift_parcel(starting_parcel, snd)?;
let mut cloud_pcl = starting_parcel;
let mut cloud_pcl_data = no_cloud_pcl_data;
let tgt_cloud_val = if no_cloud_pcl_data.1 {
false
} else if !cloud_pcl_data.1 {
true
} else {
unreachable!()
};
if tgt_cloud_val {
cloud_pcl.temperature += CelsiusDiff(1.0);
cloud_pcl_data = lift_parcel(cloud_pcl, snd)?;
} else {
no_cloud_pcl.temperature += CelsiusDiff(-1.0);
no_cloud_pcl_data = lift_parcel(no_cloud_pcl, snd)?;
}
while no_cloud_pcl_data.1 || !cloud_pcl_data.1 {
if tgt_cloud_val {
cloud_pcl.temperature += CelsiusDiff(1.0);
cloud_pcl_data = lift_parcel(cloud_pcl, snd)?;
} else {
no_cloud_pcl.temperature += CelsiusDiff(-1.0);
no_cloud_pcl_data = lift_parcel(no_cloud_pcl, snd)?;
}
}
let mut dt = cloud_pcl.temperature - no_cloud_pcl.temperature;
while dt > CelsiusDiff(0.1) {
let mid_t = no_cloud_pcl.temperature + CelsiusDiff(dt.unpack() / 2.0);
let test_pcl = Parcel {
temperature: mid_t,
..no_cloud_pcl
};
let test_pcl_data = lift_parcel(test_pcl, snd)?;
if test_pcl_data.1 {
cloud_pcl = test_pcl;
cloud_pcl_data = test_pcl_data;
} else {
no_cloud_pcl = test_pcl;
no_cloud_pcl_data = test_pcl_data;
}
dt = cloud_pcl.temperature - no_cloud_pcl.temperature;
}
let t0 = no_cloud_pcl.temperature;
let dt0 = lowest_level_parcel(snd)?.temperature - t0;
let e0 = no_cloud_pcl_data.0;
let d_e = cloud_pcl_data.0 - e0;
Ok((t0, dt0, e0, d_e))
}
fn plume_parcels(
snd: &Sounding,
plus_range: CelsiusDiff,
increment: CelsiusDiff,
) -> Result<impl Iterator<Item = Parcel>> {
let parcel = mixed_layer_parcel(snd)?;
let (row, mut parcel) = find_parcel_start_data(snd, &parcel)?;
if row.temperature.unwrap() > parcel.temperature {
parcel.temperature = row.temperature.unwrap();
}
let max_t = parcel.temperature + plus_range;
Ok(PlumeParcelIterator {
next_p: parcel,
max_t,
increment,
})
}
struct PlumeParcelIterator {
next_p: Parcel,
max_t: Celsius,
increment: CelsiusDiff,
}
impl Iterator for PlumeParcelIterator {
type Item = Parcel;
fn next(&mut self) -> Option<Self::Item> {
let next_t = self.next_p.temperature + self.increment;
if next_t > self.max_t {
None
} else {
self.next_p = Parcel {
temperature: next_t,
..self.next_p
};
Some(self.next_p)
}
}
}
fn lift_parcel(parcel: Parcel, snd: &Sounding) -> Result<(JpKg, bool)> {
let (lcl_pressure, _lcl_temperature) = metfor::pressure_and_temperature_at_lcl(
parcel.temperature,
parcel.dew_point,
parcel.pressure,
)
.ok_or(AnalysisError::MetForError)?;
let (_, parcel) = find_parcel_start_data(snd, &parcel)?;
let theta = parcel.theta();
let theta_e = parcel.theta_e()?;
let dry_mw = parcel.mixing_ratio()?;
let calc_parcel_t = |tgt_pres| {
if tgt_pres > lcl_pressure {
let t_k = metfor::temperature_from_theta(theta, tgt_pres);
metfor::virtual_temperature(
t_k,
metfor::dew_point_from_p_and_mw(tgt_pres, dry_mw)?,
tgt_pres,
)
.map(Kelvin::from)
} else {
metfor::temperature_from_theta_e_saturated_and_pressure(tgt_pres, theta_e)
.and_then(|t_c| metfor::virtual_temperature(t_c, t_c, tgt_pres))
.map(Kelvin::from)
}
};
let snd_pressure = snd.pressure_profile();
let hgt = snd.height_profile();
let env_t = snd.temperature_profile();
let env_dp = snd.dew_point_profile();
let (integrated_cape, cloud) = izip!(snd_pressure, hgt, env_t, env_dp)
.filter_map(|(p, h, env_t, env_dp)| {
if p.is_some() && h.is_some() && env_t.is_some() && env_dp.is_some() {
Some((p.unpack(), h.unpack(), env_t.unpack(), env_dp.unpack()))
} else {
None
}
})
.filter(move |(p, _, _, _)| *p < parcel.pressure)
.filter_map(|(p, h, env_t, env_dp)| {
calc_parcel_t(p).map(|pcl_t| (p, h, env_t, env_dp, pcl_t))
})
.filter_map(|(p, h, env_t, env_dp, pcl_t)| {
metfor::virtual_temperature(env_t, env_dp, p)
.map(|env_vt| (p, h, Kelvin::from(env_vt), pcl_t))
})
.take_while(|(_p, _h, env_t, pcl_t)| pcl_t > env_t)
.fold(
(
(0.0, false),
Meters(std::f64::MAX),
Kelvin(0.0),
Kelvin(0.0),
),
|acc, (p, h, env_t, pcl_t)| {
let ((mut cape, mut cloud), prev_h, prev_env_t, prev_pcl_t) = acc;
let dz = h - prev_h;
if dz <= Meters(0.0) {
((cape, cloud), h, env_t, pcl_t)
} else {
cape += ((pcl_t - env_t).unpack() / env_t.unpack()
+ (prev_pcl_t - prev_env_t).unpack() / prev_env_t.unpack())
* dz.unpack();
cloud = p < lcl_pressure;
((cape, cloud), h, env_t, pcl_t)
}
},
)
.0;
Ok((JpKg(integrated_cape / 2.0 * -metfor::g), cloud))
}