3D cmor example

[larsbuntemeyer]

Here is an example of how to prepare a dataset for cmorization that has a depth dimension including bounds which is created from single layer temperature fields:

def tsl(TD3, TD4, TD5, TD, TDCL):
    """Create 3D dataset for cmorization from 5 single layer fields"""
    # layer thicknesses
    delta = np.array([0.065, 0.254, 0.913, 2.902, 5.7])
    # layer depths as mids of cumulative sums of thicknesses
    # depth dimension name is a key from the cmor coordinates table, e.g. 'sdepth' (not the out name 'depth')
    # https://github.com/WCRP-CORDEX/cordex-cmip6-cmor-tables/blobf1d5e98e4ef4d5538f4169403abfe61d4637884f/Tables/CORDEX-CMIP6_coordinate.json#L647-L669
    sdepth = xr.DataArray(np.cumsum(delta) - 0.5 * delta, dims="sdepth")
    # layer bounds
    bounds = np.column_stack((np.cumsum(delta) - delta, np.cumsum(delta)))
    # set appropriate CF metadata for cf_xarray to understand
    sdepth.attrs["units"] = "m"
    sdepth.attrs["axis"] = "Z"
    sdepth.attrs["bounds"] = "depth_bnds"
    # create dataset
    # concat fields along the sdepth coordinate
    res = xr.concat([TD3, TD4, TD5, TD, TDCL], dim=sdepth)
    res.name = "tsl"
    res.attrs["units"] = "K"
    res = res.to_dataset()
    # assign bounds coordinate
    res = res.assign_coords(depth_bnds=xr.DataArray(bounds, dims=("sdepth", "bnds")))

    return res.transpose("time", "sdepth", "rlat", "rlon", ...)