From e5528d01351f67c0caab3db3a28163d06407880f Mon Sep 17 00:00:00 2001 From: Graham McIntire Date: Wed, 15 Apr 2026 18:36:09 -0500 Subject: [PATCH] NarrClient.extract_profile_from_file/3 + fetch_profile_at/2 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit extract_profile_from_file shells out to cdo: cdo -outputtab,name,lev,value -remapnn,lon=X_lat=Y parses the text output (varNNN ↔ NCEP GRIB1 param number lookup), maps surface records to era5_profiles fields, builds the pressure profile from each (TMP, HGT, SPFH) triple, and merges in the derived fields from SoundingParams.derive/1. Returns the same attrs shape as Era5BatchClient.build_profile_attrs/5 so the existing bulk_insert_profiles path can ingest NARR rows unchanged. fetch_profile_at picks records from a fetched inventory, byte-range fetches each one into a per-record temp file, runs cdo -merge to build a composite, calls extract_profile_from_file, and cleans up via try/after. Drive-by fix: ThetaE.dewpoint_from_spfh/2 returns Kelvin (its existing test asserts that explicitly). Era5BatchClient was passing that Kelvin value straight into "dwpc" (Celsius), which would then crash SoundingParams.compute_refractivity_profile when it called Buck's saturation-vapor-pressure equation with t=294 instead of t=21. The bug never surfaced because era5_profiles is empty in production — no ERA5 backfill has ever succeeded. Fixed in both era5_batch_client and the new narr_client by subtracting 273.15 in the wrapper. Both wrappers now have a comment pointing at the K→C conversion. Test fixture test/fixtures/narr/narr_dfw_2010-06-15_12z.grb is the spike-captured composite (1.1 MB — Lambert conformal projection isn't sellonlatbox-subsettable so we can't shrink it further). --- .../weather/era5_batch_client.ex | 6 +- lib/microwaveprop/weather/narr_client.ex | 282 +++++++++++++++++- .../fixtures/narr/narr_dfw_2010-06-15_12z.grb | Bin 0 -> 1154646 bytes .../weather/narr_client_test.exs | 92 ++++++ 4 files changed, 376 insertions(+), 4 deletions(-) create mode 100644 test/fixtures/narr/narr_dfw_2010-06-15_12z.grb diff --git a/lib/microwaveprop/weather/era5_batch_client.ex b/lib/microwaveprop/weather/era5_batch_client.ex index e1f6853a..2bd81be4 100644 --- a/lib/microwaveprop/weather/era5_batch_client.ex +++ b/lib/microwaveprop/weather/era5_batch_client.ex @@ -337,9 +337,13 @@ defmodule Microwaveprop.Weather.Era5BatchClient do # ERA5 pressure-level moisture is specific humidity (kg/kg). Convert to # dewpoint in °C via the Magnus-Tetens inverse already used for HRRR native # profiles. `level_pa` is the pressure level itself in Pa. + # ThetaE.dewpoint_from_spfh/2 returns Kelvin. Profile `dwpc` entries + # are Celsius (the downstream SoundingParams.compute_refractivity_profile + # passes p["dwpc"] straight to the Buck saturation-vapor-pressure + # equation which expects °C). Convert here. defp dewpoint_from_spfh(nil, _level_pa), do: nil defp dewpoint_from_spfh(spfh, _level_pa) when spfh <= 0, do: nil - defp dewpoint_from_spfh(spfh, level_pa), do: ThetaE.dewpoint_from_spfh(spfh, level_pa) + defp dewpoint_from_spfh(spfh, level_pa), do: ThetaE.dewpoint_from_spfh(spfh, level_pa) - 273.15 # -- Bulk insert ----------------------------------------------------------- diff --git a/lib/microwaveprop/weather/narr_client.ex b/lib/microwaveprop/weather/narr_client.ex index a18c9be8..33822458 100644 --- a/lib/microwaveprop/weather/narr_client.ex +++ b/lib/microwaveprop/weather/narr_client.ex @@ -9,14 +9,44 @@ defmodule Microwaveprop.Weather.NarrClient do architecture — the filename is a historical artifact from an earlier MERRA-2 plan. - Currently provides URL/time helpers, the inventory parser, and the - inventory HTTP fetch. The byte-range `fetch_profile_at/2` pipeline - lands in a follow-up task of the same plan. + Provides URL/time helpers, the inventory parser, the inventory HTTP + fetch, the cdo-based composite-file profile extractor, and the + end-to-end byte-range `fetch_profile_at/2` pipeline. """ + alias Microwaveprop.Weather.SoundingParams + alias Microwaveprop.Weather.ThetaE + @base_url "https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/3-hourly" @valid_hours [0, 3, 6, 9, 12, 15, 18, 21] + # Surface records we need from the inventory. Tuples of {var, level} keyed + # exactly as `parse_inventory/2` returns them. Pressure-level records are + # added dynamically below from whatever (TMP|HGT|SPFH):"NN mb" tuples the + # parsed inventory contains. + @surface_records [ + {"HPBL", "sfc"}, + {"PRES", "sfc"}, + {"TMP", "2 m above gnd"}, + {"DPT", "2 m above gnd"}, + {"PWAT", "atmos col"} + ] + + # cdo emits NCEP GRIB1 parameters as `varNNN`. This table maps the cdo + # parameter token to the semantic NCEP short name, verified empirically + # in the spike (see plan doc cdo lookup table). + @cdo_var_to_name %{ + "var1" => "PRES", + "var7" => "HGT", + "var11" => "TMP", + "var17" => "DPT", + "var33" => "UGRD", + "var34" => "VGRD", + "var51" => "SPFH", + "var54" => "PWAT", + "var221" => "HPBL" + } + @doc """ Builds the NCEI HTTPS URL for the NARR analysis GRIB1 file at `valid_time`. @@ -110,6 +140,252 @@ defmodule Microwaveprop.Weather.NarrClient do end end + @doc """ + Extracts an `era5_profiles`-shaped attribute map from a composite NARR + GRIB1 file at the given lat/lon point. + + Pure-ish: shells out to `cdo -outputtab,name,lev,value -remapnn,lon=…_lat=…` + on the local file. No network. + + Returns `{:ok, attrs}` where `attrs` has the same shape + `Era5BatchClient.build_profile_attrs/5` produces, or `{:error, reason}` + on missing file, cdo failure, or empty parser output. + """ + @spec extract_profile_from_file(Path.t(), float(), float()) :: + {:ok, map()} | {:error, term()} + def extract_profile_from_file(grb_path, lat, lon) when is_binary(grb_path) and is_number(lat) and is_number(lon) do + if File.exists?(grb_path) do + run_cdo_outputtab(grb_path, lat, lon) + else + {:error, "NARR cdo extract: file not found: #{grb_path}"} + end + end + + @doc """ + Fetches all the NARR records we need for `valid_time` via byte-range + `Range:` GETs against the NCEI grib file, merges them into a composite + GRIB1 with `cdo -merge`, then extracts the profile at `{lat, lon}` + using `extract_profile_from_file/3`. + + Returns `{:ok, profile_attrs}` shaped like `Era5BatchClient.build_profile_attrs/5` + or `{:error, reason}`. Temp files are always cleaned up via `try/after`. + """ + @spec fetch_profile_at(DateTime.t(), {float(), float()}) :: + {:ok, map()} | {:error, term()} + def fetch_profile_at(%DateTime{} = valid_time, {lat, lon}) when is_number(lat) and is_number(lon) do + grb_url = url_for(valid_time) + + with {:ok, index} <- fetch_inventory(valid_time), + {:ok, record_keys} <- pick_records(index) do + do_fetch_profile(grb_url, index, record_keys, lat, lon) + end + end + + defp do_fetch_profile(grb_url, index, record_keys, lat, lon) do + tag = "narr_#{System.unique_integer([:positive])}" + tmp_dir = System.tmp_dir!() + merged_path = Path.join(tmp_dir, "#{tag}_merged.grb") + + record_paths = + Enum.with_index(record_keys, fn key, i -> + {key, Path.join(tmp_dir, "#{tag}_#{i}.grb")} + end) + + try do + with :ok <- fetch_records(grb_url, index, record_paths), + :ok <- merge_records(record_paths, merged_path) do + extract_profile_from_file(merged_path, lat, lon) + end + after + Enum.each(record_paths, fn {_key, path} -> File.rm(path) end) + File.rm(merged_path) + end + end + + # Builds the list of {var, level} keys we want to pull from the inventory: + # the fixed surface set plus every (TMP|HGT|SPFH) at "NN mb" present in + # the index. Returns `{:error, _}` if any of the required surface records + # are missing from the inventory. + defp pick_records(index) do + missing_surface = Enum.reject(@surface_records, &Map.has_key?(index, &1)) + + if missing_surface == [] do + pressure_keys = + index + |> Map.keys() + |> Enum.filter(fn + {var, level} when var in ["TMP", "HGT", "SPFH"] -> + String.ends_with?(level, " mb") + + _ -> + false + end) + + {:ok, @surface_records ++ pressure_keys} + else + {:error, "NARR inventory missing required surface records: #{inspect(missing_surface)}"} + end + end + + defp fetch_records(grb_url, index, record_paths) do + Enum.reduce_while(record_paths, :ok, fn {key, path}, _acc -> + {offset, length} = Map.fetch!(index, key) + + case fetch_byte_range(grb_url, offset, length, path) do + :ok -> {:cont, :ok} + {:error, reason} -> {:halt, {:error, reason}} + end + end) + end + + defp fetch_byte_range(url, offset, length, dest_path) do + range_header = "bytes=#{offset}-#{offset + length - 1}" + + case Req.get(url, [headers: [{"range", range_header}], receive_timeout: 60_000] ++ req_options()) do + {:ok, %{status: status, body: body}} when status in [200, 206] and is_binary(body) -> + File.write(dest_path, body) + + {:ok, %{status: status}} -> + {:error, "NARR byte-range HTTP #{status} for #{range_header}"} + + {:error, reason} -> + {:error, "NARR byte-range request failed: #{inspect(reason)}"} + end + end + + defp merge_records(record_paths, merged_path) do + inputs = Enum.map(record_paths, fn {_key, path} -> path end) + args = ["-merge"] ++ inputs ++ [merged_path] + + case System.cmd("cdo", args, stderr_to_stdout: true) do + {_output, 0} -> :ok + {output, code} -> {:error, "cdo -merge exited #{code}: #{output}"} + end + end + + defp run_cdo_outputtab(grb_path, lat, lon) do + args = [ + "-outputtab,name,lev,value", + "-remapnn,lon=#{lon}_lat=#{lat}", + grb_path + ] + + case System.cmd("cdo", args, stderr_to_stdout: false) do + {stdout, 0} -> + parse_cdo_outputtab(stdout) + + {output, code} -> + {:error, "cdo extract exited #{code}: #{output}"} + end + end + + defp parse_cdo_outputtab(stdout) do + raw = + stdout + |> String.split("\n", trim: true) + |> Enum.map(&String.trim/1) + |> Enum.filter(&String.starts_with?(&1, "var")) + |> Enum.reduce(%{}, &accumulate_cdo_row/2) + + if raw == %{} do + {:error, "NARR cdo extract: no parseable values in cdo output"} + else + {:ok, build_profile_attrs(raw)} + end + end + + defp accumulate_cdo_row(line, acc) do + case parse_cdo_row(line) do + {key, value} -> Map.put(acc, key, value) + :error -> acc + end + end + + defp parse_cdo_row(line) do + with [var_token, lev_str, value_str] <- String.split(line, ~r/\s+/, trim: true), + {:ok, name} <- Map.fetch(@cdo_var_to_name, var_token), + {lev_int, ""} <- Integer.parse(lev_str), + {value_float, ""} <- Float.parse(value_str) do + {{name, lev_int}, value_float} + else + _ -> :error + end + end + + defp build_profile_attrs(raw) do + profile = build_pressure_profile(raw) + + base = %{ + profile: profile, + surface_temp_c: kelvin_to_celsius(raw[{"TMP", 2}]), + surface_dewpoint_c: kelvin_to_celsius(raw[{"DPT", 2}]), + surface_pressure_mb: pa_to_mb(raw[{"PRES", 0}]), + hpbl_m: raw[{"HPBL", 0}], + pwat_mm: raw[{"PWAT", 0}] + } + + params = SoundingParams.derive(profile) + + if params do + Map.merge(base, %{ + surface_refractivity: params.surface_refractivity, + min_refractivity_gradient: params.min_refractivity_gradient, + ducting_detected: params.ducting_detected || false, + duct_characteristics: params.duct_characteristics + }) + else + Map.merge(base, %{ + surface_refractivity: nil, + min_refractivity_gradient: nil, + ducting_detected: false, + duct_characteristics: nil + }) + end + end + + defp build_pressure_profile(raw) do + raw + |> Map.keys() + |> Enum.filter(fn + {"TMP", lev} when lev > 100 -> true + _ -> false + end) + |> Enum.map(fn {"TMP", lev} -> lev end) + |> Enum.sort(:desc) + |> Enum.map(&pressure_level_entry(&1, raw)) + |> Enum.reject(&is_nil/1) + end + + defp pressure_level_entry(lev_pa, raw) do + with t_k when is_float(t_k) <- raw[{"TMP", lev_pa}], + h_m when is_float(h_m) <- raw[{"HGT", lev_pa}], + spfh when is_float(spfh) <- raw[{"SPFH", lev_pa}] do + %{ + "pres" => lev_pa / 100.0, + "tmpc" => kelvin_to_celsius(t_k), + "hght" => h_m, + "dwpc" => dewpoint_from_spfh(spfh, lev_pa) + } + else + _ -> nil + end + end + + # ThetaE.dewpoint_from_spfh/2 returns Kelvin (the function adds 273.15 + # at the end of the Magnus-Tetens inversion). Profile entries store + # `dwpc` in *Celsius* — the `c` is load-bearing — and downstream + # SoundingParams.compute_refractivity_profile passes p["dwpc"] straight + # to the Buck saturation-vapor-pressure function which expects °C. So + # convert here. + defp dewpoint_from_spfh(spfh, _lev_pa) when not is_number(spfh) or spfh <= 0, do: nil + defp dewpoint_from_spfh(spfh, lev_pa), do: ThetaE.dewpoint_from_spfh(spfh, lev_pa * 1.0) - 273.15 + + defp kelvin_to_celsius(nil), do: nil + defp kelvin_to_celsius(k) when is_number(k), do: k - 273.15 + + defp pa_to_mb(nil), do: nil + defp pa_to_mb(pa) when is_number(pa), do: pa / 100.0 + defp do_get_inventory(url) do case Req.get(url, [receive_timeout: 30_000] ++ req_options()) do {:ok, %{status: 200, body: body}} when is_binary(body) -> diff --git a/test/fixtures/narr/narr_dfw_2010-06-15_12z.grb b/test/fixtures/narr/narr_dfw_2010-06-15_12z.grb new file mode 100644 index 0000000000000000000000000000000000000000..d20eb656ad6abd1c54a97c1b7239b9434044faa0 GIT binary patch literal 1154646 zcmeFadsLHWx<0H<1;SL@+6siq83Lg#Q)vZ4tB$q;p)IYn0zo+ik`P)^P$5A;}vXX>=hl+(2MUH41(>_2~Dt*>j%{;l0L zD+}Z$ym_wYzV7S3?)!NsFqFS{{bS|&b?dgBF)tbU{HM}&8!XM9dU759h5p<66uo+` zd;F92QR}y^`^j9c-=TFm&5P@@eoz1V`MO`fNB{e;|9DPXK+T?4$`-#byKK|Ol^lR6!xZkZ|aqsUo{T9}+xQjo|;;yV=aTgs- z|74$=e(A~@7WeDjroV5F>6fmoVR7%DHvN6in11P^gX!1aH~rGbYgpXH-KM|q4bv}O z%i=DcHvN6qu(&G@roZi47Prv!xBYIn>DP{%e(B>i?A*01?xg9ru$IMLH2rNKuVHal zcANgb-s<8srl{qD5s*REl4A3K=-wrg43Uen+9$6DOqn|=#xS=_&w{J^AOFhqYuB>4e>MGW*Rr_(Si{?nnSSkx>6bqKkrwy+HQd7gE-o$BuS{RDe}Tm{ z{bnCO9y9&oe}SEA`prIm{3Gq$wcNu0E-tOsX?pqp|GzQQ-*PQG*YtP&pGFJ+o6rc` zABh(J?_Jz~0V(z)>D;x@!dfb{UzyG!Yf%foH~nqb!ntc%+(y%H;m2Cs-KO8dT6XRl z7We<+@^=l3yB1gMM_SyqxMGc_Z=3&Nu2{M0pKL9z*jg6X!Sq}BkrsC?X8Brn?iv<% zEoS+zO<%CJSdL$>;cb5;v;1G66S0QH{ca6A_eWaX#oeZZwuTozPR75VR7%T z;TG1JzPO8PxP{M6$M?H6+``vuz&G3UTUcDfEqreJEv&3FH+|cz;TEz@zlHDCa0~ZM z-{mW77$^RE4Y!bO`Yn97&h%Y=e+{>=&h%ZrxQ1KEUT1QwV{uQLx^KQ-XZkK*XZkHH zuHhCwH~kjy$xNY}m2Z9zF=)yy%+Q4U`myZ)^DTV;`F{<={~vDwAO64pTmJ{k&Dl@9 z@$^A6%W5};FycbHLKs%Sj^{@{X}Htk@wNF@xfZdZu$17pTP{6=nqO3bU3s5pui zaV=gFYdI>lL@8!x;SCIUgK3rcV0F7TIl4y8%Eg~t%7|5Qq%BFcYE_FwHKS5X@s;|+ zT8>;TR*Csl3CbB>ae|7&&Wfh{i=#V@sli{1u9-IG#?bA?E7jaG5iTyn9}?8!IzuY` zwNzov47WC3IU}tV<9%E3&RV)51-_cmUO94c#@KzhR;_PIRH^y3@hY(Ze@HB-64j~= zokuUKm#%$OJEOb&A%5)=?oW5+@I{Ue_s8KjB&x*n1_}OTadeGHpCq`@&I{Oo-i@{Q z{R?jF0IU!n--_7Gi+{+ugsul1g)sb1V^U;~Q_+VGbC&B4T79vjKihd|hIe+dKj)gpyZfCWn;py*@iGw|^Zc$r}r$0<210*YNM zRI1sdLy>!Qaxw+q7F2T-_|X+_eZ>UNbl?83BMh&xp!MO054Qe%~BR#}1)ODt5W z@nNtoY5>5>h^Acma7$sh!Yi3xrl3KBruH)BA@oV@>EvUv4^9rzY zSXoinIjJ!>Ln5CB~HsLq{z>u8;N>B>jgO)X;8jJyFa(rU>?j{X|Kcr6)U0jL@6<)qpfJ+1?8 zNwjlfus5o(eFSAO*e_{^sztnz_s&f%?R|H3aoV|ifg6{A;s(EW-gU`qD>iM&vy!`J zpZv$0)|#cqsoU~ys8U<}aMu8e_fMc+EpMHBxVxm4Pj`cHn|Y5c$FRc3N>JINem@2>e^<`{$R^f2702Dho3?PT$*HjeF6q3w!04fg1fHB3~GT7kzS6 z4c;Cspg_R3>sI5JLXvKm0T01RQHG^6q+x^b3h16|ffhRhZ*QP@p)NiwA8vz>BX6)Zcv}`+R-`=gO=gSM;FEP!qL^z zE~d|h9j2VY;s7Z8p>}QP{__b&&|Di{^5az@0E&gMTID}h@cta!$a7hZotgiZG}JzQ zA}La1mlb|t*7N1FlNQ}Cm$iEIyjB*k9Q7@m?3&?Ky7Fy<02Gsjk|P+>BA~Ez3d&Xr zA#GqLJ%L6L*oX!EvvP(@bV$j|4f*YYNFGTm9!V?y1B(Kvun(}2!aleFauWb^ES;M% zEt+TvNOQAds=JB-Xkr2S;|w1Vdy9-jhX&#U0!qxoz8BGQilxYxVr>IBj23%_2_{xgN}Ysv z`$WCF0Q{$6bsxDXCYn?0k&0$nt3ajHX}CQ^f}Mwd07{Fb>9(36L8>!6EDnd2y05&w zW8#DJQ;WAk&P{Y(eznl4Z-$46BA^)lL*B}FjGzyuPH8;Gp3jtL*)3#dzi+i8$B-6p z(Z9XI#jxS^au?+Uugs-rVrO{**DRp2O|_~N5Fe`Oa4a@t#VT6TfcF7e8n}do6g?wo zAR1vlG_n;^RvZxu@{u%>mM-)IKqf}k(52!HRJ4-(P@p9U3&*!&Y%_f-$rP}tnhydM z>4X4^7bdtctxLud%cGI*_@VT>2po}-qLu~ef}|M(a!bS7Gh#;`hJfV*9yfl-`{$;4 zA^Xoy^Ft82^v6XCqHq;buL$Y}?46AK!6x+}AQfZ)P-?zfj7U*aRA9w_0IEB+PQ$*1 zh~jl0z9?Z`ORnP(9!N?bUDU8Zj1u)#Bl6WkAcX+@){+=M@P&=Qw10EM2>ii6bN`kP*a~9HctU zONs(exnN*H3lcMBST!Z8RiKtKaGexfH?56E^aw~%7?4e;b=d9TKP;z!LQM;~fFFR7 zmIx>T!dhfVBZ{P0K#HjrNh99}QS1U^N~6P&dxS(^lm|&_AQd4lRngcM{Ads@m#Cf{ z@Gf$IPssU%r2v2|rEMms7X=_gbPJI9u!tOPHSU>}5exK-L8__-$UNQU!(g-Ph!5hm zSym5#;@(0;u`ea5INAmX2LnN2xB_-Gb{_>iJ;&-ag5hI8i++uMS8cc=uE93L{SnbZ zQIH0yO;9HId3&9)v);EPbDVP<&mWf+yizJ+AyKuH2%P2yzkhbxCU<*@utu>IW6?+f z1Xg1&k{3}1hh0XZ4XRl!MTk?P$7ZGcMtK1*2qsk|qDK$$-o6C2J}m|Hxp9MOw=+S4 zbcc*cS+@%MfCX-1;&)m{ETxSMJ3nRg=buY31RxrvktiA<=_q;OD;YlYqJo z6xjVoE+R*z)XgALaBHOWRUk$Ui_`<;IWFl*LA(v5*U$3%k1~OXm(Ihgc4G$72sWX} zPz>sjR-5ItSnc`ljGfK)bIHxtxz4A@7HTbWUHk^#JBMHbP|Qd0xNcSCV$-u-R(s+5 z3&NpJx9bXLztYab=f6FCLKqWx+po%%XTP^{if#6K*`)H(?(!+tmNzQfw9r=!03su0 zT|IaMoO2NTg3UUuN1X!xquYRRW%v&|jFJt00MWt70J_t(gBK76bTJTY*NG3{4T9fC zsFA!h7NAHHP}ZP=XnH9t5=@v%S^&R+4HAAHsW-TCO5y$Xc{c%Miklp=5i1V36{I2= zdZl1W07g|)dhg(Kt_IJbbFNzNch0#Q{b|RB;@1QOvI^<5oyoQOtT>h164=5`myA15 zwR{LlV|c9-{6|2!v-?45?4Hz%GYWhxHT1m-<5pUx5xO zc_m{tD1FjLz&`LRa8izn@}qrAE8_0Cv;J_j&Ldz?QM)wxpuA~FyrXoQ6h%8PH0(mF zMc8xstrk7obL*?ytHBO#SB^AwRJ&y=ckJ9BVP&@0sVG6g39M>UIKNrms@w2-MFQLG z&B`gMEuXS3acnMCxX4U+0SfYHMrV^2K)w{S%fwi0l2LrrM~F&X3_uAjM8y(?h7KwS z3LXlXh*87~Mxugvhd9oNh7M!O8`6&4v}m9JIgNS?pVF8V3CSKDQq`uj52$k0`y)R1 z$Zx2{!F4DBP&@`>)9}FZ2?p|V&iF8LdqA)7dJC~d8JiEZMA9&-5{N}r zVye7IzsQqOhD)nKeIyc)qrhu=gcU=567IE>IKUP zSBl?-{)0mxiimX1kF&VsSrYG|HHZyjK~~H~_bSV?PC@(6B;48h&U+u&D${mXCMcZW zs!9+bWhL-z->7s|IOC(rt3^a^iKzS`>5$67YsR!lWP70iNt1*+Q~Gl8;W+}b3AosQ zVir>(Q&FQ`B{$LxloH6saMO`_!N(p4J;c_iq2{QL2ler&j)zVp0tLSKLr9Ckq37DP zp8U!-Dd(+oZH53~gxETdoyZ7Wsur2A{NDjkK7o}9a@&2lc*~oWu5#y~Dp#H7+vlc? z#4pgD1Qc|Rl&TIr%3(;2@M_o-5EFXJQD83QnNE~YEJT!s+akKSO6#nuXeB%fNC-ZZ zpoV|@7?1%mO4tw`YH3+KggU;}fDY>bw+LL43`NW6Q=`3KL1UujhVdh-h6hiEUE)a%4KZ0jdsy2f-PqbOYNW z`3v6%|9J&hxyo%3E;>)(Vf>+@Ki8)70pS^Z7{Wt|qUXF-*~YWlTbUrR+*{Gc-HflK zRzX#5DsQALb^wtji9jjI869ziK}cVff#p>wp*Tiz91N(Bm~a8FCBViY_k#_xb{&tjw_s;Suu4N!F+KjLgSnVrsQ(%SKq{U(?%@aX_)kJPsc8E`+L+G^`J_Qph z@E^XyVk@Bu$)y8_iFDBd1QS(2_&59pd4gI-5Kz6PfIJejnG%$W{vTW+gq02@C_W&x zg-(o2!rvQm4k3`)^fYj+4y7-I2&m;EXXzF~jf-gHfd8;rv#^1QX^>c~$|o&4 z;n&&p_(S$EEI`>_nsf*bWOZ!uICqy%esf~2N6AEAv~^*E-uwB2`_`8?=EThHn9OwX zyS1ZaEbQ=k(UxyTf-^32E4wPDKz$Wc9?QF6D=lM%^v*$$bTGigM1nI4GTA6T@B%5s zD+6;$L0`Bv!i_-TGFTRn_RfN9pwD2ij3Zxwb}`^V1_MdUq4}o`KozPF$q+;f?qo`` zT;ve}76;*nT1XWNhBMgH9 z4VU19D3D>1dK2-ns)YzEVn-ttSkmDW^%B8wTAe%Z3M>n4abnGks`D_EAkikfM*5{<%X7u8iI**NTZ+10ENIbN z?XMIDB|mrO{$iTrnb_Hze=Btfx;X;rwuApxMaRI7vPr{+H!2f&Adl9e`Ca7<@wCA? zu2s{nAI93F1gErcJ{jh&adAW8@&@8%b;m zFl<()7^tR1MIwvPM{x`B^uS`%`@xhpT->HI@2NNn#6uRl74U;rFN1LNPLV0b8J`$Gfs!8g`(%)x-0Jyve>tR6QYc!r%n8NRiO6wkUMi1^zWH zbDg3Rj<#9^3b%T>WYVwUm9s4t8Jp8ulk;4SGhs9a)D zK+G^dg}7t-8e|=yOCbhKhui}SZ(Ur${c)|{ zK%`iL+$`qbK)H^>Jzh5ITM{pO@Jg9WV*jh9lYY1TN?jsme9I7Jua~tNHUKTbmVl}$ zg?%8Tr1$#>Y% %M9-?gXwMG&+^E?yPG^YGK+PVq_FP+vUS3@X6Isw7_rDw;U#q z0wPc_d>o>qHV$y3w4BjhLlMlo&A1vzv9lql#Qdtn%BiCMVC+-=a^aXVPokWC13aYWNV(=f~ z!6+UBxu#Izj9IgqTxX~-CI=V4wrbijl^6Mc@VL=B0-)?N^orOhGS&R7co-zPXlQAL zD2Q0z$U=fGdnpSsl!mg;Vr6$Zgpxm8uu-4V*2ICAN`#A@+e=#JZv6~WxBum{6N}yW zs6oR%rIX4DIMMnIyAhI>ZxUpCDp7ugUYM4`HiQRH1|2qtTA?O0sRy>J@iOIPEH*W- zuymbLeEe0bojO>xQ1wgwsr3?aFx-qu(de}`d<5$2M>)(#2AkzxKU|AVD6EJt>I8f9 zx4c!->i%Fuu47RF`81&nPzT_-LE=p7=tHw%LPH_PsKiXWOI{sYGO>8;mD2c#+piY2 zxcB2zrhoBG;JI%%zFxrATKSi?E`sXghlo%*Cj6my?L*H`Nj<4BNsX)(8z_frq1PZo zNsE#c=Z8k2{926)w37_b*=}TsLS!MiQa%z>76J-LKq#ew7SqU}@Y6HN3krn;?j(^< z20Oh0_zx+cx?23O{pS+&&HXk`1@mc{EYo^|fSS1uL<-leTY@HXQ+9=%S<>*wW@1~WvN-Unrx0Gkjd=Wcgxoix0gdX%Ytxtt|AQx85CZ(|SgJg&L z0#Hb3Fk8tkcY|F{7?E$Agu+-02Uny@B|AS7Y{B+N5iGTR1M9pDOSub#*GkPbm!jcU zN+#yiuEi*zxC&uVi$}3@N5czeTe_~iP(0D0_AH)IUcut@^?Q{}77e^w>LPmpeEUY} zQ?8(vq}XbglDA2Xs#PZBr+o;?4&xjy(}x$`S#XbFV}SaQ0}xS{Ay>7Gz_^xz>pCU) z9LZ7S@Fmr@^KAECoKz-1pY3Sl#V>K5yXFajcIR(rVuyk)aveRqSa}m&{oduEbpR!R z-3IFgTSZHTGvfxZv9r}&+7?Bx;uhujEs0=|YY@dLi7h%F?7$iT zt;xRm*UIA!W_v5)*M$~{U>G8F%fJyv!dzNKi%r>uOEGOoga60HiX%~B#dTztQ`d;B zY+Q$)D+jWX>Mhc-OwdsWkm`?BR|D!0b<}E~gfJrIf!6*7Hz}H;2?oz#$h-{*E8G=X z%E(8F_z2F6JufPpl~qm(XNzLzqLv`4*o#^`;aFdK0kMp)ghdUW#m-$