prop/lib/microwaveprop/weather/gefs_client.ex

defmodule Microwaveprop.Weather.GefsClient do
  @moduledoc """
  Client for the NCEP Global Ensemble Forecast System (GEFS) served
  anonymously from NOMADS at
  `https://nomads.ncep.noaa.gov/pub/data/nccf/com/gens/prod/`.

  GEFS is the 0.5° / 31-member physics-based global ensemble. It runs
  four times per day (00/06/12/18Z) with 3-hourly output to +240 h and
  6-hourly output to +384 h. The `pgrb2ap5/` product family contains
  the subset of surface and pressure-level variables the propagation
  scorer needs for extended-horizon (Day 2-7) outlooks, once HRRR's
  18 h horizon runs out.

  This module exposes the URL, time, and message-inventory helpers.
  Byte-range idx parsing and HTTP range fetch logic are shared with
  `Microwaveprop.Weather.HrrrClient` since the GRIB2/idx format is
  identical; the GEFS-specific pieces are the NOMADS path shape, the
  3-then-6-hourly forecast-hour list, and the Magnus dewpoint
  derivation (GEFS pgrb2a publishes RH rather than Td).
  """

  @base_url "https://nomads.ncep.noaa.gov/pub/data/nccf/com/gens/prod"

  @run_hours [0, 6, 12, 18]

  # Variables the propagation scorer reads. All present in pgrb2a f000-f384
  # per an April 2026 .idx probe. DPT is absent from GEFS output — derive
  # it from TMP + RH via Magnus (see `dewpoint_from_rh/2`). HPBL is also
  # absent; the scorer tolerates a nil HPBL without penalising the score.
  @surface_messages [
    %{var: "TMP", level: "2 m above ground"},
    %{var: "RH", level: "2 m above ground"},
    %{var: "PRES", level: "surface"},
    %{var: "PWAT", level: "entire atmosphere (considered as a single layer)"},
    %{var: "UGRD", level: "10 m above ground"},
    %{var: "VGRD", level: "10 m above ground"},
    %{var: "TCDC", level: "entire atmosphere"},
    %{var: "APCP", level: "surface"}
  ]

  @doc """
  Rounds a wall-clock `DateTime` down to the nearest GEFS run cycle
  (00/06/12/18Z) on the same UTC day.
  """
  @spec nearest_run(DateTime.t()) :: DateTime.t()
  def nearest_run(%DateTime{} = dt) do
    run_hour = dt.hour |> div(6) |> Kernel.*(6)

    dt
    |> DateTime.truncate(:second)
    |> Map.put(:hour, run_hour)
    |> Map.put(:minute, 0)
    |> Map.put(:second, 0)
  end

  @doc """
  Builds the NOMADS URL for the GEFS ensemble-mean (`geavg`) pgrb2a 0.5°
  GRIB2 file at the given `run_date`, `run_hour`, and `forecast_hour`.
  """
  @spec ensmean_url(Date.t(), non_neg_integer(), non_neg_integer()) :: String.t()
  def ensmean_url(%Date{} = run_date, run_hour, forecast_hour) when run_hour in @run_hours and forecast_hour >= 0 do
    date_str = Calendar.strftime(run_date, "%Y%m%d")
    hh = pad2(run_hour)
    fff = pad3(forecast_hour)

    "#{@base_url}/gefs.#{date_str}/#{hh}/atmos/pgrb2ap5/geavg.t#{hh}z.pgrb2a.0p50.f#{fff}"
  end

  @doc "Appends `.idx` to a GRIB2 URL to locate its wgrib2 byte-range index."
  @spec idx_url(String.t()) :: String.t()
  def idx_url(grib_url) when is_binary(grib_url), do: grib_url <> ".idx"

  @doc """
  Returns the list of GEFS forecast hours published per run: 3-hourly
  from f000 through f240, then 6-hourly from f246 through f384.
  """
  @spec forecast_hours() :: [non_neg_integer()]
  def forecast_hours do
    short = Enum.to_list(0..240//3)
    long = Enum.to_list(246..384//6)
    short ++ long
  end

  @doc """
  GRIB2 message descriptors the propagation scorer needs from the
  GEFS pgrb2a product. Values match the `var` and `level` strings that
  appear in the file's `.idx` sidecar exactly.
  """
  @spec surface_messages() :: [%{var: String.t(), level: String.t()}]
  def surface_messages, do: @surface_messages

  @doc """
  Derives dewpoint in °C from air temperature (°C) and relative humidity
  (%) using the Magnus formula with the Bolton (1980) coefficients.

  GEFS pgrb2a publishes 2 m RH but not 2 m DPT, so downstream code that
  consumed HRRR's `surface_dewpoint_c` must call this helper to fill in
  the equivalent value. RH is clamped to a small positive floor to keep
  the logarithm finite at 0 %.
  """
  @spec dewpoint_from_rh(float(), float()) :: float()
  def dewpoint_from_rh(temp_c, rh_pct) when is_number(temp_c) and is_number(rh_pct) do
    a = 17.625
    b = 243.04
    rh = max(rh_pct, 0.01) / 100.0
    gamma = :math.log(rh) + a * temp_c / (b + temp_c)
    b * gamma / (a - gamma)
  end

  defp pad2(n), do: n |> Integer.to_string() |> String.pad_leading(2, "0")
  defp pad3(n), do: n |> Integer.to_string() |> String.pad_leading(3, "0")
end