add groundstations.py

components/groundstation.py

+#!/usr/bin/env python3
+
+import numpy as np
+import sys
+from argparse import ArgumentParser, FileType
+from pathlib import Path
+
+this_dir = Path(__file__).parent
+
+def lonlat(string):
+    pair = string.split(",")
+    return (float(pair[0]), float(pair[1]))
+
+parser = ArgumentParser()
+parser.add_argument("--descriptorfile", type=FileType('r'), default=this_dir.parent/"descriptorfiles"/"icon_1mom.py", help="") 
+parser.add_argument("--coords", type=lonlat, nargs="+", required=True, help="Coordinates of the groundstations")
+parser.add_argument("--frequency", type=float, nargs="+", default=94)
+parser.add_argument("--dt", type=str, default=None, help="Timestep")
+parser.add_argument("--source_component", type=str, default="atmo_output", help="component name of the source component")
+parser.add_argument("--source_grid", type=str, default="icon_atmos_grid", help="grid name of the source grid")
+args = parser.parse_args()
+
+from yac import YAC, UnstructuredGrid, Location, Field, TimeUnit, \
+    InterpolationStack, NNNReductionType, Reduction, Action
+import pyPamtra
+
+yac = YAC()
+
+comp = yac.def_comp("pamtra")
+
+coords = np.array(args.coords)
+ncoords = len(args.coords)
+grid = UnstructuredGrid("pamtra-grid", [3]*ncoords,
+                        [*coords[0,:], *(coords[0,:]+0.1), *coords[0,:]],
+                        [*coords[1,:], *coords[1,:], *(coords[1,:]+0.1)],
+                        zip(range(ncoords), range(ncoords, 2*ncoords), range(2*ncoords, 3*ncoords)))
+points = grid.def_points(Location.CELL,
+                         coords[0,:], coords[1,:])
+
+yac.sync_def()
+
+#print(yac.get_field_names(*source), file=sys.stderr)
+
+interp = InterpolationStack()
+interp.add_nnn(NNNReductionType.AVG, 1, 1.0)
+
+source = (args.source_component, args.source_grid)
+
+dt = args.dt or yac.get_field_timestep(*source, "temp")
+
+fields = {}
+for name in ["temp", "pres", "qv", "qi", "qc", "qs", "qr", "qg", "z_ifc", "u", "v", "w"]:
+    nlevel = yac.get_field_collection_size(*source, name)
+    fields[name] = Field.create(name, comp, points, nlevel,
+                                dt, TimeUnit.ISO_FORMAT)
+    yac.def_couple(*source, name,
+                   "pamtra", "pamtra-grid", name,
+                   dt, TimeUnit.ISO_FORMAT, Reduction.TIME_NONE,
+                   interp)
+
+yac.enddef()
+
+with open(args.descriptorfile, "r") as dfile:
+    exec(dfile.read()) # adds a variable "descriptorFile"
+
+def call_pamtra(temp, pres, qv, qi, qc, qs, qr, qg, z_ifc, u, v, w,
+                datetime):
+    pamtra = pyPamtra.pyPamtra()
+    for hydro in descriptorFile:
+        pamtra.df.addHydrometeor(hydro)
+    pamtra.nmlSet["passive"] = False
+    pamtra.nmlSet["active"] = True
+    pamtra.nmlSet["radar_mode"] = "simple"
+    rh = 100*pyPamtra.meteoSI.q2rh(qv, temp, pres)
+    pamtra.createProfile(hgt_lev=z_ifc.T[None, :, ::-1],
+                         press=pres.T[None, :, ::-1],
+                         temp=temp.T[None, :, ::-1],
+                         relhum=rh.T[None, :, ::-1],
+                         hydro_q=np.stack([qc.T, qi.T, qr.T, qs.T, qg.T], axis=-1)[:, ::-1, :],
+                         )
+    pamtra.runPamtra(args.frequency)
+    pamtra.writeResultsToNetCDF(f"output_{datetime}.nc", xarrayCompatibleOutput=True)
+
+while True:
+    datetime = fields["temp"].datetime
+    data = {}
+    for v, f in fields.items():
+        data[v], info = f.get()
+        print(v, data[v].shape, file=sys.stderr)
+    call_pamtra(**data, datetime=datetime)
+    if info == Action.GET_FOR_RESTART:
+        break
+

descriptorfiles/icon_1mom.py

+import numpy as np
+# Define descriptorFiles
+descriptorFile = np.array([ # TODO to be reviewed, coefficients for m-D and v-D are changin
+  
+   ('cwc_q', 1.0,  1, -99.0,   -99.0, -99.0,  -99.0, -99.0,  3,  1,   'mono',           -99.0, -99.0, -99.0, -99.0,  2.0e-5,  -99.0, 'mie-sphere', 'corPowerLaw_24388657.6_2.0', -99.0),
+   ('iwc_q', 0.2, -1, -99.0,   130.0,   3.0,  0.684,   2.0,  3,  1,   'mono_cosmo_ice', -99.0, -99.0, -99.0, -99.0,   -99.0,  -99.0, 'ssrg-rt3_0.18_0.89_2.06_0.08', 'corPowerLaw_30.606_0.5533', -99.0),
+   ('rwc_q', 1.0,  1, -99.0,   -99.0, -99.0,  -99.0, -99.0,  3,  100, 'exp',            -99.0, -99.0, 8.0e6, -99.0,  1.2e-4, 6.0e-3, 'mie-sphere', 'corPowerLaw_130.0_0.5', -99.0),
+   ('swc_q', 0.6, -1, -99.0,   0.038,   2.0, 0.3971,  1.88,  3,  100, 'exp_cosmo_snow', -99.0, -99.0, -99.0, -99.0, 5.1e-11, 1.0e-2, 'ssrg-rt3_0.25_1.00_1.66_0.04', 'corPowerLaw_4.9_0.25', -99.0),
+   ('gwc_q', 1.0, -1, -99.0,   169.6,   3.1,  -99.0, -99.0,  3,  100, 'exp',            -99.0, -99.0, 4.0e6, -99.0, 1.0e-10, 1.0e-2, 'mie-sphere', 'corPowerLaw_406.67_0.85', -99.0)
+],)
+
+#pam.df.addHydrometeor(df)