The Pycmor Cookbook#

A showcase of some more complicated use cases.

If you’d like to contribute with your own recipe, or ask for a recipe, please open a documentation issue on our GitHub repository.

Default Unit conversions#

In this example, we show how to do “default unit conversions”. pycmor can handle some standard cases, assuming your NetCDF files are sensible. We show two examples:

mmolC/m2/d –> kg m-2 s-2, as in fgco2 for CMIP6_Omon.json:

"fgco2": {
    "frequency": "mon",
    "modeling_realm": "ocnBgchem",
    "standard_name": "surface_downward_mass_flux_of_carbon_dioxide_expressed_as_carbon",
    "units": "kg m-2 s-1",
    "cell_methods": "area: mean where sea time: mean",
    "cell_measures": "area: areacello",
    "long_name": "Surface Downward Mass Flux of Carbon as CO2 [kgC m-2 s-1]",
    "comment": "Gas exchange flux of CO2 (positive into ocean)",
    "dimensions": "longitude latitude time depth0m",
    "out_name": "fgco2",
    "type": "real",
    "positive": "down",
    "valid_min": "",
    "valid_max": "",
    "ok_min_mean_abs": "",
    "ok_max_mean_abs": ""
},

µatm –> Pa, as in spco2 for CMIP6_Omon.json:

"spco2": {
    "frequency": "mon",
    "modeling_realm": "ocnBgchem",
    "standard_name": "surface_partial_pressure_of_carbon_dioxide_in_sea_water",
    "units": "Pa",
    "cell_methods": "area: mean where sea time: mean",
    "cell_measures": "area: areacello",
    "long_name": "Surface Aqueous Partial Pressure of CO2",
    "comment": "The surface called 'surface' means the lower boundary of the atmosphere. The partial pressure of a dissolved gas in sea water is the partial pressure in air with which it would be in equilibrium. The partial pressure of a gaseous constituent of air is the pressure which it alone would exert with unchanged temperature and number of moles per unit volume. The chemical formula for carbon dioxide is CO2.",
    "dimensions": "longitude latitude time depth0m",
    "out_name": "spco2",
    "type": "real",
    "positive": "",
    "valid_min": "",
    "valid_max": "",
    "ok_min_mean_abs": "",
    "ok_max_mean_abs": ""
},

You can download test data with the provided script:

$ ./download-example-ata.sh

This extracts ten years of example data for the two variables in questions.

In our configuration file, we don’t need to specify anything extra, since the default pipeline can handle these cases.

Dealing with wrong input units#

In this example we show how to deal with incorrect units set on the input data. The motivating example comes from the variables dissic and talk for oceanic biogeochemistry:

cmip6-cmor-tables/Tables/CMIP6_Omon.json#

    "dissic": {
        "standard_name": "mole_concentration_of_dissolved_inorganic_carbon_in_sea_water",
        "units": "mol m-3",
    },
    "talk": {
        "standard_name": "sea_water_alkalinity_expressed_as_mole_equivalent",
        "units": "mol m-3",
    },

These two variables need DIC and Talk, which are saved as bgc02 and bgc03 by REcoM. Unfortunately, the units attribute in the NetCDF files is wrong, e.g.:

Excerpt of ncdump -h bgc02_3d.nc#

 float bgc02(time, nodes_3d) ;
         bgc02:description = "bgc tracer 02" ;
         bgc02:units = "" ;
         bgc02:grid_type = "unstructured" ;
         bgc02:_FillValue = 1.e+30f ;

The actual input unit in both cases is \(\frac{mmol}{m^{3}}\), i.e. use mult_factor = 1/1e3 to convert from mmol to mol.

Happily enough, pycmor provides an easy way to rectify this without needing to edit the source files with something like ncatted. You can specify the correct units in the model_unit setting of the rule. Here is an example of how to do that:

general:
  name: "AWI-ESM-1-1-lr PI Control"
  description: "CMOR configuration for AWIESM 1.1 LR"
  maintainer: "pgierz"
  email: "pgierz@awi.de"
  cmor_version: "CMIP6"
  mip: "CMIP"
  CMIP_Tables_Dir: "/work/ab0246/a270077/SciComp/Projects/pycmor/cmip6-cmor-tables/Tables"
  CV_Dir: "/work/ab0246/a270077/SciComp/Projects/pycmor/cmip6-cmor-tables/CMIP6_CVs/"
pycmor:
  # parallel: True
  warn_on_no_rule: False
  use_flox: True
  dask_cluster: "slurm"
  dask_cluster_scaling_mode: fixed
  fixed_jobs: 1
  # minimum_jobs: 8
  # maximum_jobs: 30
  # You can add your own path to the dimensionless mapping table
  # If nothing is specified here, it will use the built-in one.
rules:
  - name: Dissolved Inorganic Carbon in Seawater
    description: "dissic from REcoM, showing missing units in NetCDF"
    inputs:
      - path: "/work/ab0246/a270077/SciComp/Projects/pymor/examples/03-incorrect-units-in-source-files/model_runs/piControl_LUtrans1850/outdata/recom/"
        pattern: bgc02_fesom_.*.nc
    grid_file: /pool/data/AWICM/FESOM1/MESHES/core/griddes.nc
    mesh_path: /pool/data/AWICM/FESOM1/MESHES/core
    cmor_variable: dissic
    model_variable: "bgc02"
    model_unit: "mmol m-3"
    output_directory: .
    variant_label: r1i1p1f1
    experiment_id: piControl
    source_id: AWI-CM-1-1-HR
    model_component: ocnBgchem
    grid_label: gn
  - name: Seawater Alkalinity
    description: "talk from REcoM, showing missing units in NetCDF"
    inputs:
      - path: "/work/ab0246/a270077/SciComp/Projects/pymor/examples/03-incorrect-units-in-source-files/model_runs/piControl_LUtrans1850/outdata/recom/"
        pattern: bgc03_fesom_.*.nc
    grid_file: /pool/data/AWICM/FESOM1/MESHES/core/griddes.nc
    mesh_path: /pool/data/AWICM/FESOM1/MESHES/core
    cmor_variable: talk
    model_variable: "bgc03"
    model_unit: "mmol m-3"
    output_directory: .
    variant_label: r1i1p1f1
    experiment_id: piControl
    source_id: AWI-CM-1-1-HR
    model_component: ocnBgchem
    grid_label: gn
distributed:
  worker:
    memory:
      target: 0.6 # Target 60% of worker memory usage
      spill: 0.7 # Spill to disk when 70% of memory is used
      pause: 0.8 # Pause workers if memory usage exceeds 80%
      terminate: 0.95 # Terminate workers at 95% memory usage
    resources:
      CPU: 4 # Assign 4 CPUs per worker
    death-timeout: 600 # Worker timeout if no heartbeat (seconds): Keep workers alive for 5 minutes
# SLURM-specific settings for launching workers
jobqueue:
  slurm:
    name: pycmor-worker
    queue: compute # SLURM queue/partition to submit jobs
    account: ab0246 # SLURM project/account name
    cores: 4 # Number of cores per worker
    memory: 128GB # Memory per worker
    walltime: '00:30:00' # Maximum walltime per job
    interface: ib0 # Network interface for communication
    job-extra-directives: # Additional SLURM job options
      - '--exclusive' # Run on exclusive nodes
      - '--nodes=1'
    # Worker template
    worker-extra:
      - "--nthreads"
      - 4
      - "--memory-limit"
      - "128GB"
      - "--lifetime"
      - "25m"
      - "--lifetime-stagger"
      - "4m"
    # How to launch workers and scheduler
    job-cpu: 128
    job-mem: 256GB
    # worker-command: dask-worker
    processes: 4 # Limited by memory per worker!
    # scheduler-command: dask-scheduler

Multivariable Inputs and vertical integration#

Here we show how use two input variables and perform a vertical integration. The motivating example is from the CMIP6_Omon.json table:

   "intpp": {
       "standard_name": "net_primary_mole_productivity_of_biomass_expressed_as_carbon_by_phytoplankton",
       "units": "mol m-2 s-1",
       "long_name": "Primary Organic Carbon Production by All Types of Phytoplankton",
       "comment": "Vertically integrated total primary (organic carbon) production by phytoplankton.  This should equal the sum of intpdiat+intpphymisc, but those individual components may be unavailable in some models.",
}

In the AWI-ESM-1-REcoM case two phytoplanktion groups contribute to primary production:

nanophytoplankton (variable diags3d01)
diatoms (variable diags3d02)

Hence, the workflow to get implemented would be:

\[\text{NPP} = \int_{\text{bottom}}^{\text{top}} (\text{NPP}_{\text{nanopythoplankton}} + \text{NPP}_{\text{diatoms}})\]

In pseudo-code:

pp = diags3d01 + diags3d03
pp.sum(dim="depth")
mmol_mol = 1e-3
per_day_to_per_sec = 1/86400
conversion_factor_mmolC_m2_d_to_molC_m2_s = 1/1e3/86400  # About 1.157407e-08

We implement custom steps for the addition and vertical integration.

Preparing the example

You should run this example on DKRZ Levante. Start by downloading the sample data:

$ ./download-example-data.sh

This will extract 10 years of example data for the two variables in question.

Next, we define some custom functions in intpp_recom.py:

def add_pp_components(data, rule):
    """Adds together Net Primary Production components"""
    # This is a USER FUNCTION and therefore you need some human brainpower and knowledge
    # of your model in order to know what you want to do.

    data["pp"] = data["diags3d01"] + data["diags3d02"]
    # This return is mandatory for the pipeline to usefully continue
    # We return our new variable to be further processed:
    return data["pp"]  # Return type: DataArray!


def set_pp_units(data, rule):
    """Corrects missing units on PP"""
    # Again, a user function. Insider knowledge of REcoM is needed
    data["pp"].attrs["units"] = "molC m-2 day-1"
    # Actually, data is in mmol, so we do a "by hand" conversion here to get the mol correct
    data["pp"] *= 1e-6
    return data


def vertical_integration(data, rule):
    data = data.sum(dim="depth")
    return data


def manual_breakpoint(data, rule):
    breakpoint()
    return data

Working with Dimensionless Units#

Problem#

You need to work with variables that have ambiguous dimensionless units in CMIP6, such as:

Pure dimensionless quantities (unit: “1”)
Percentage values (unit: “%”)
Ratios (e.g., “kg kg-1”)
Salinity values (unit: “0.001”)
Parts-per-million concentrations (unit: “1e-06”)

And you need these units to be properly converted or recognized by Pycmor. In essence, we need to tell Pycmor more about the physical meaning of the units so that it can convert from one unit to another arbitrarily. An example is a mass ratio of 0.001. If mass is not specified that is ambiguous because it could also be, for example, a volume ratio. The ratios of mass and volume are different depending on the density. To help Pycmor to convert to 0.001 (mass ratio), we need to tell Pycmor that it is a mass ratio and we do this by indicating that 0.001 means g/kg or 10-3 kg kg-1.

Solution#

First, check if the variable exists in the dimensionless mappings file. The file is typically located at:

<your_pycmor_installation>/src/pycmor/data/dimensionless_mappings.yaml

For example: /Users/username/Codes/pycmor/src/pycmor/data/dimensionless_mappings.yaml

Open this file and search for your variable name (e.g., “sisali”) to see if it already exists.

If your variable exists but has an empty mapping (or if it’s missing entirely), you need to add the appropriate mapping:

# For salinity variables (with unit "0.001")
sisali:  # sea_ice_salinity
  "0.001": g/kg

# For pure dimensionless variables (with unit "1")
abs550aer:  # atmosphere_absorption_optical_thickness_due_to_aerosol_particles
  "1":

# For variables with mole fraction units
co2:  # mole_fraction_of_carbon_dioxide_in_air
  "1e-06":

If you have added a new mapping, you can now use it in your regular Pycmor workflow. The cmorize function will automatically use the dimensionless mapping to interpret and convert the units correctly.
To contribute your dimensionless mappings back to the Pycmor repository:
1. Fork the Pycmor repository on GitHub: esm-tools/pycmor
2. Clone your fork and create a branch for your changes
3. Update the dimensionless_mappings.yaml file with your additions/corrections
4. Commit your changes with a descriptive message explaining the mappings you’ve added
5. Push your changes and create a pull request to the main repository
Your contributions help improve Pycmor for the entire climate science community!

The Pycmor Cookbook

Contents