Creating warm start forcing data

We want to add or replace soil data in the forcing data files. We show an example here, where we take the soil data from the output file of a remo run and add this to the forcing file of another run for use as initial lateral boundary conditions.

[15]:

import xarray as xr

import pyremo as pr

afile = xr.open_dataset(
    "/work/ch0636/remo/forcing-data/EUR-11/ERA5/vc_27lev/1979/a056530a1979010100.nc"
)
tfile = pr.remo_ds.open_remo_dataset(
    "/work/ch0636/remo/forcing-data/era5_warm/a056524a1979010100",
    update_meta=True,
    parse_time=True,
)

The output dataset (tfile) contains the spin up soil data. This could also come from another point in time of the same run or another run. The soil is defined by a number of codes or variable names:

[18]:

# soil_codes = [54,55,56,206,207,208,209,232,170,183,84,140,194,141]
warm_codes = [54, 206, 207, 208, 209, 232, 170, 183, 140, 194, 141]
warm_vars = [pr.codes.get_dict(c)["variable"] for c in soil_codes]
warm_vars

[18]:

['TSL', 'TSN', 'TD3', 'TD4', 'TD5', 'GLAC', 'TD', 'TDCL', 'WS', 'WL', 'SN']

Now we copy and replace the data from the output file (tfile) in the forcing file (afile).

[19]:

for var in soil_vars:
    afile[var] = tfile[var]

We should check the data with a simple plot, we will look at the temperatures here.

[20]:

def get_temps(ds):
    temps = ["TSW", "TSI", "TSL", "TSN", "TD3", "TD4", "TD5", "TD", "TDCL"]
    dim = xr.DataArray(data=temps, dims="var", name="var")
    tda = xr.concat([ds[var] - 273.5 for var in temps], dim=dim)
    tda.name = "temperature"
    return tda


get_temps(afile).plot(col="var", col_wrap=3)

[20]:

<xarray.plot.facetgrid.FacetGrid at 0x2ba23714bf70>

Adding static fields to the forcing file

The initial forcing file also has to contain static fields (e.g., like the land sea mask or the orography). Remo will read these fields in from the first initial forcing file.

[21]:

from pyremo import data as rd

surflib = rd.surflib("EUR-11", update_meta=True)
surflib

We have to take care if the coordinate values in the afile and surface library are not exactly the same. To make sure, we set join="override" to simply use the coordinate values from the afile.

[22]:

afile = xr.merge([afile, surflib], join="override")
afile

Now, we can write back the afile to Netcdf. We have to take care of correct encoding of missing values here.

[24]:

fillvars = ["TSW", "SEAICE", "TSI"]

for var in afile:
    # we have to set the encoding here explicitly, otherwise xarray.to_netcdf will
    # encode missing values by NaN, which will crash REMO...
    if var in fillvars:
        afile[var].encoding["_FillValue"] = 1.0e20
    else:
        afile[var].encoding["_FillValue"] = None

afile.to_netcdf(
    "/work/ch0636/remo/forcing-data/EUR-11/ERA5/vc_27lev/1979/a056530a1979010100_warm_soil.nc"
)