Class updates: add wind, piecewise lapse rate, and time-dependent surface fluxes + various QOL improvements (#144)

* Implement basic zoom; TODO: fix angle in skewT

* Add logarithmic zoom for skew-T diagram; TODO: T calculations should use real axis extent rather than fixes base/top pressure

* Fix skew-T lines responding to original extent instead of actual; now they don't tilt anymore

* Add panning effect, but it is stroboscopic and doesn't work for skewT yet

* combine side-effects for both axes in a single callback

* Remove animationframe

* Use produce to update both scales in a single call

* Don't update panstart; this fixes the jittering

* Also work in log space

* Make consistent for x-direction

* zoom towards cursor

* Add reset plot button

* Thinner lines

* Higher resolution plot

* Round time to steps of 10 minutes

* Use runClass from package, skipping BMI altogher + rich metadata for output vars; however, app/model hangs on second or third run

* Fix hanging issue: wrap only once...

* Make sure initial state is included in output

* Reset pan/zoom when variable changes

* formatting

* Use output metadata in plot labels and variable pickers

* ditch BMI

* Fix xlabel in timeseries plot; fix axes extent for non-time on x-axis

* Remove config.JSON everywhere; only keep the embedded version

* Change units to use superscript

* Drop _0 for initial fields

* Make wq number array in form/config; TODO: interpolate value in model

* Make gammaq and gammatheta arrays; TODO: fix profiles in app, use in model

* Add free-troposphere anchor points for theta and q; TODO: use in model and app; update reference/default/test configs

* Add wind variables to config. TODO: use in model / display profiles in app?

* Implement hourly interpolation of surface fluxes + add them to output + move output spec to separate module

* Implement piecewise lapse rate in model

* formatting

* Fix display of piecewise lapse rate in web app

* Add basic wind logistics; rethink optional variable groups; use getValue rather than direct getters to access nested object groups; TODO in config only allow wind if mixed layer is on; TODO complete wind implementation

* Complete wind implementation in model; include it in timeseries plot for validation

* Add varnavas reference config as preset

* Add u and v to profile plots

* Make sure thermodynamic plot works with piecewise lapse rate

* Determine axis extents of profile plots posterior to generating the profiles, so lapse part is fully included

* Fix tests, set some values back to their previous default

* Make wind from observations also visible in profile plots

* Added Varnavas config with observations to test folder + increased max height in preset

* Make labels of du and dv unique

* Combine tooltip and label + for number[] input show extra description about comma seperator

* Fix types in tests

* fix type issue

---------

Co-authored-by: sverhoeven <[email protected]>

Author: Peter9192

apps/class-solid/src/components/Analysis.tsx

@@ -1,5 +1,5 @@
 import type { Config } from "@classmodel/class/config";
-import { type ClassOutput, outputVariables } from "@classmodel/class/runner";
+import { type ClassOutput, outputVariables } from "@classmodel/class/output";
 import * as d3 from "d3";
 import { saveAs } from "file-saver";
 import { toBlob } from "html-to-image";
@@ -226,33 +226,16 @@ export function VerticalProfilePlot({
   const variableOptions = {
     "Potential temperature [K]": "theta",
     "Specific humidity [kg/kg]": "q",
+    "u-wind component [m/s]": "u",
+    "v-wind component [m/s]": "v",
   };
 
   const classVariable = () =>
     variableOptions[analysis.variable as keyof typeof variableOptions];
 
   const observations = () =>
     flatObservations().map((o) => observationsForProfile(o, classVariable()));
-  const obsAllX = () =>
-    observations().flatMap((obs) => obs.data.map((d) => d.x));
-  const obsAllY = () =>
-    observations().flatMap((obs) => obs.data.map((d) => d.y));
-
-  const allValues = () => [
-    ...flatExperiments().flatMap((e) =>
-      e.output ? e.output[classVariable()] : [],
-    ),
-    ...obsAllX(),
-  ];
-  const allHeights = () => [
-    ...flatExperiments().flatMap((e) => (e.output ? e.output.h : [])),
-    ...obsAllY(),
-  ];
 
-  // TODO: better to include jump at top in extent calculation rather than adding random margin.
-  const xLim = () => getNiceAxisLimits(allValues(), 1);
-  const yLim = () =>
-    [0, getNiceAxisLimits(allHeights(), 0)[1]] as [number, number];
   const profileData = () =>
     flatExperiments().map((e) => {
       const { config, output, ...formatting } = e;
@@ -271,6 +254,19 @@ export function VerticalProfilePlot({
       };
     });
 
+  const allX = () => [
+    ...profileData().flatMap((p) => p.data.map((d) => d.x)),
+    ...observations().flatMap((obs) => obs.data.map((d) => d.x)),
+  ];
+  const allY = () => [
+    ...profileData().flatMap((p) => p.data.map((d) => d.y)),
+    ...observations().flatMap((obs) => obs.data.map((d) => d.y)),
+  ];
+
+  // TODO: better to include jump at top in extent calculation rather than adding random margin.
+  const xLim = () => getNiceAxisLimits(allX(), 1);
+  const yLim = () => [0, getNiceAxisLimits(allY(), 0)[1]] as [number, number];
+
   function chartData() {
     return [...profileData(), ...observations()];
   }

apps/class-solid/src/lib/download.ts

@@ -1,4 +1,4 @@
-import type { ClassOutput } from "@classmodel/class/runner";
+import type { ClassOutput } from "@classmodel/class/output";
 import { BlobReader, BlobWriter, ZipWriter } from "@zip.js/zip.js";
 import { toPartial } from "./encode";
 import type { ExperimentConfig } from "./experiment_config";

apps/class-solid/src/lib/presets.ts

@@ -8,13 +8,15 @@ import { overwriteDefaultsInJsonSchema } from "@classmodel/form/utils";
 import type { DefinedError, JSONSchemaType, ValidateFunction } from "ajv";
 // TODO replace with preset of a forest fire
 import deathValley from "./presets/death-valley.json";
+import varnavas from "./presets/varnavas.json";
 
 const presetConfigs = [
   {
     name: "Default",
     description: "The classic default configuration",
   },
   deathValley,
+  varnavas,
 ] as const;
 
 export interface Preset {

apps/class-solid/src/lib/presets/death-valley.json

@@ -1,19 +1,19 @@
 {
   "name": "Death Valley",
   "description": "Preset with Death Valley conditions",
-  "theta_0": 323,
-  "h_0": 200,
-  "dtheta_0": 1,
-  "q_0": 0.008,
-  "dq_0": -0.001,
+  "theta": 323,
+  "h": 200,
+  "dtheta": 1,
+  "q": 0.008,
+  "dq": -0.001,
   "dt": 60,
   "runtime": 43200,
   "wtheta": [0.1],
   "advtheta": 0,
-  "gammatheta": 0.006,
-  "wq": 0.0001,
+  "gammatheta": [0.006],
+  "wq": [0.0001],
   "advq": 0,
-  "gammaq": 0,
+  "gammaq": [0],
   "divU": 0,
   "beta": 0.2
 }

apps/class-solid/src/lib/presets/varnavas.json

@@ -0,0 +1,42 @@
+{
+  "name": "Varnavas",
+  "description": "EWED Reference configuration",
+  "h": 665.4086303710938,
+  "theta": 301.971435546875,
+  "dtheta": 0.889312744140625,
+  "gammatheta": [0.0015336532378569245, 0.006792282219976187],
+  "z_theta": [1667.025390625, 3720],
+  "q": 0.010794480331242085,
+  "dq": -0.001073240302503109,
+  "gammaq": [-2.382889988439274e-6, -2.5458646177867195e-6],
+  "z_q": [1360.52392578125, 3720],
+  "divU": 1.0981982995872386e-5,
+  "u": -1.6769592761993408,
+  "du": -1.1726601123809814,
+  "gamma_u": [-0.001995581667870283, -0.002064943080767989],
+  "z_u": [1573.5614013671875, 3720],
+  "v": -10.087028503417969,
+  "dv": -1.5727958679199219,
+  "gamma_v": [-0.0005029001622460783, 0.004509935155510902],
+  "z_v": [3037.957763671875, 3720],
+  "wtheta": [
+    0.33137200988670706, 0.28088430353706567, 0.21148639643745845,
+    0.1354453158978963, 0.0558820378352221, 0.0001602229911650289, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0
+  ],
+  "wq": [
+    1.9756698748096824e-5, 2.111953654093668e-5, 2.2018328309059143e-5,
+    2.2002774130669422e-5, 1.9648385205073282e-5, 1.6331905499100685e-5,
+    1.3947577826911584e-5, 1.2534562301880214e-5, 1.1651739441731479e-5,
+    1.1526963135111146e-5, 1.3352756468520965e-5, 1.3273539479996543e-5
+  ],
+  "sw_wind": true,
+  "advtheta": 0,
+  "advq": 0,
+  "advu": 0,
+  "advv": 0,
+  "beta": 0.2,
+  "ustar": 0.3,
+  "dt": 60,
+  "runtime": 43200
+}

apps/class-solid/src/lib/profiles.ts

@@ -1,5 +1,6 @@
 import type { Config } from "@classmodel/class/config";
-import type { ClassOutput } from "@classmodel/class/runner";
+import type { ClassOutput } from "@classmodel/class/output";
+import { findInsertIndex } from "@classmodel/class/utils";
 import type { Point } from "~/components/plots/Line";
 import type { Observation } from "./experiment_config";
 
@@ -15,31 +16,126 @@ export function getVerticalProfiles(
     return [];
   }
 
-  // Extract height profile
-  const height = output.h.slice(t)[0];
-  const dh = 1600; // how much free troposphere to display?
-  const hProfile = [0, height, height, height + dh];
   if (variable === "theta") {
-    // Extract potential temperature profile
-    const theta = output.theta.slice(t)[0];
+    let z = output.h.slice(t)[0];
+    let theta = output.theta.slice(t)[0];
     const dtheta = output.dtheta.slice(t)[0];
     const gammatheta = config.gammatheta;
-    const thetaProfile = [
-      theta,
-      theta,
-      theta + dtheta,
-      theta + dtheta + dh * gammatheta,
+    const z_theta = config.z_theta;
+    const maxHeight = z_theta.slice(-1)[0];
+
+    // Mixed layer
+    const profile = [
+      { x: theta, y: 0 },
+      { x: theta, y: z },
     ];
-    return hProfile.map((h, i) => ({ x: thetaProfile[i], y: h }));
+
+    // Inversion
+    theta += dtheta;
+    profile.push({ x: theta, y: z });
+
+    // Free troposphere
+    while (z < maxHeight) {
+      const idx = findInsertIndex(z_theta, z);
+      const lapse_rate = gammatheta[idx] ?? 0;
+      const dz = z_theta[idx] - z;
+      z += dz;
+      theta += lapse_rate * dz;
+      profile.push({ x: theta, y: z });
+    }
+    return profile;
   }
+
   if (variable === "q") {
-    // Extract humidity profile
-    const q = output.q.slice(t)[0];
+    let z = output.h.slice(t)[0];
+    let q = output.q.slice(t)[0];
     const dq = output.dq.slice(t)[0];
     const gammaq = config.gammaq;
-    const qProfile = [q, q, q + dq, q + dq + dh * gammaq];
-    return hProfile.map((h, i) => ({ x: qProfile[i], y: h }));
+    const z_q = config.z_q;
+    const maxHeight = z_q.slice(-1)[0];
+
+    // Mixed layer
+    const profile = [
+      { x: q, y: 0 },
+      { x: q, y: z },
+    ];
+
+    // Inversion
+    q += dq;
+    profile.push({ x: q, y: z });
+
+    // Free troposphere
+    while (z < maxHeight) {
+      const idx = findInsertIndex(z_q, z);
+      const lapse_rate = gammaq[idx] ?? 0;
+      const dz = z_q[idx] - z;
+      z += dz;
+      q += lapse_rate * dz;
+      profile.push({ x: q, y: z });
+    }
+    return profile;
+  }
+
+  if (config.sw_wind && variable === "u") {
+    let z = output.h.slice(t)[0];
+    let u = output.u.slice(t)[0];
+    const du = output.du.slice(t)[0];
+    const gammau = config.gamma_u;
+    const z_u = config.z_u;
+    const maxHeight = z_u.slice(-1)[0];
+
+    // Mixed layer
+    const profile = [
+      { x: u, y: 0 },
+      { x: u, y: z },
+    ];
+
+    // Inversion
+    u += du;
+    profile.push({ x: u, y: z });
+
+    // Free troposphere
+    while (z < maxHeight) {
+      const idx = findInsertIndex(z_u, z);
+      const lapse_rate = gammau[idx] ?? 0;
+      const dz = z_u[idx] - z;
+      z += dz;
+      u += lapse_rate * dz;
+      profile.push({ x: u, y: z });
+    }
+    return profile;
+  }
+
+  if (config.sw_wind && variable === "v") {
+    let z = output.h.slice(t)[0];
+    let v = output.v.slice(t)[0];
+    const dv = output.dv.slice(t)[0];
+    const gammav = config.gamma_v;
+    const z_v = config.z_v;
+    const maxHeight = z_v.slice(-1)[0];
+
+    // Mixed layer
+    const profile = [
+      { x: v, y: 0 },
+      { x: v, y: z },
+    ];
+
+    // Inversion
+    v += dv;
+    profile.push({ x: v, y: z });
+
+    // Free troposphere
+    while (z < maxHeight) {
+      const idx = findInsertIndex(z_v, z);
+      const lapse_rate = gammav[idx] ?? 0;
+      const dz = z_v[idx] - z;
+      z += dz;
+      v += lapse_rate * dz;
+      profile.push({ x: v, y: z });
+    }
+    return profile;
   }
+
   return [];
 }
 
@@ -142,6 +238,8 @@ export function getThermodynamicProfiles(
   const h = output.h.slice(t)[0];
   const gammaTheta = config.gammatheta;
   const gammaq = config.gammaq;
+  const z_theta = config.z_theta;
+  const z_q = config.z_q;
 
   const nz = 25;
   let dz = h / nz;
@@ -169,11 +267,18 @@ export function getThermodynamicProfiles(
   soundingData.push({ p, T, Td });
 
   // Free troposphere
+  let z = zArrayMixedLayer.slice(-1)[0];
   dz = 200;
   while (p > 100) {
-    theta += dz * gammaTheta;
-    q += dz * gammaq;
+    // Note: idx can exceed length of anchor points, then lapse becomes undefined and profile stops
+    const idx_th = findInsertIndex(z_theta, z);
+    const lapse_theta = gammaTheta[idx_th];
+    const idx_q = findInsertIndex(z_q, z);
+    const lapse_q = gammaq[idx_q];
+    theta += dz * lapse_theta;
+    q += dz * lapse_q;
     p += pressureDiff(T, q, p, dz);
+    z += dz;
     T = thetaToT(theta, p);
     Td = dewpoint(q, p);
 
@@ -194,10 +299,23 @@ export function observationsForProfile(obs: Observation, variable = "theta") {
       const p = obs.pressure[i];
       const theta = tToTheta(T, p);
       const q = calculateSpecificHumidity(T, p, rh);
-      if (variable === "theta") {
-        return { y: h, x: theta };
+      const { u, v } = windSpeedDirectionToUV(
+        obs.windSpeed[i],
+        obs.windDirection[i],
+      );
+
+      switch (variable) {
+        case "theta":
+          return { y: h, x: theta };
+        case "q":
+          return { y: h, x: q };
+        case "u":
+          return { y: h, x: u };
+        case "v":
+          return { y: h, x: v };
+        default:
+          throw new Error(`Unknown variable '${variable}'`);
       }
-      return { y: h, x: q };
     }),
   };
 }
@@ -217,3 +335,15 @@ export function observationsForSounding(obs: Observation) {
     linestyle: "none",
   };
 }
+
+function windSpeedDirectionToUV(
+  speed: number,
+  directionDeg: number,
+): { u: number; v: number } {
+  const directionRad = (directionDeg * Math.PI) / 180;
+
+  const u = -speed * Math.sin(directionRad); // zonal (east-west)
+  const v = -speed * Math.cos(directionRad); // meridional (north-south)
+
+  return { u, v };
+}

apps/class-solid/src/lib/runner.ts

@@ -1,5 +1,6 @@
 import type { Config } from "@classmodel/class/config";
-import type { ClassOutput, runClass } from "@classmodel/class/runner";
+import type { ClassOutput } from "@classmodel/class/output";
+import type { runClass } from "@classmodel/class/runner";
 import { wrap } from "comlink";
 
 const worker = new Worker(new URL("./worker.ts", import.meta.url), {

apps/class-solid/src/lib/store.ts

@@ -2,7 +2,7 @@ import { createUniqueId } from "solid-js";
 import { createStore, produce, unwrap } from "solid-js/store";
 
 import type { Config } from "@classmodel/class/config";
-import type { ClassOutput } from "@classmodel/class/runner";
+import type { ClassOutput } from "@classmodel/class/output";
 
 import {
   mergeConfigurations,