diff --git a/components/groundstation.py b/components/groundstation.py new file mode 100644 index 0000000000000000000000000000000000000000..d0ca43878f45415dbb2edeec695263b013263131 --- /dev/null +++ b/components/groundstation.py @@ -0,0 +1,93 @@ +#!/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 + diff --git a/descriptorfiles/icon_1mom.py b/descriptorfiles/icon_1mom.py new file mode 100644 index 0000000000000000000000000000000000000000..c385361663977f3d93d5b37a9204bca65ab43153 --- /dev/null +++ b/descriptorfiles/icon_1mom.py @@ -0,0 +1,12 @@ +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)