Source

MOZART-global to AQMEII-NA / README.md

Full commit

table of contents

  1. description
    1. problem
    2. goal
    3. plan
  2. implementation

description

problem

I have a 3D global N2O inventory (i.e., atmospheric N2O concentrations for each cell in a 3D grid), which is the input to this project (which is a small part of a larger project involving CMAQ). The input is unprojected, with grid dimensions=(longitude, latitude, hybrid sigma-pressure). I need to convert that input into a form (my output) that can be consumed by a particular CMAQ run. For that, my output must be

Note also that the input and output differ in both the vertical and horizontal extents: the layer heights, layer count, and top height all differ.

goal

3D-regrid from the global geographical input to the regional LCC output (as conservatively as possible). Unfortunately this particular 3D-regridding is not currently (early 2013) well-supported by existing tools; particularly

  • GRASS, NCL, and R interpolations do not support 3D directly
  • ESMF regridding tools do not support spherical->Cartesian directly
  • no tools support vertical HσP directly

plan

constraints

Plan and implementation constraints include:

nodal grids vs center grids

I have found only one tool that currently supports 3D regridding in any form, ESMF_RegridWeightGen. The ESMF regridding tools consider LCC to be a Cartesian grid, so I will need to convert my input to Cartesian form (i.e., (x,y,z) where the dimensions are all distances, not angles). But my input and output grids are also not grids in the same sense as ESMF uses. The input and output grids attribute values (concentrations) to the spaces interior to sets of grid nodes: implicitly, they attribute values to centers of hexahedronal gridcells. But the ESMF regridding tools operate on grids of nodes, i.e., the vertices or corners of hexahedrons. So, from both the input and output, I need to extract the grid of nodes.

two-pass regridding

The ESMF regridding tools work by

  1. generating weights that map (size of) the input grid to the output grid
  2. applying those weights to the values of the input grid

... but, currently, only for spherical grids. For Cartesian grids, weight application is not currently supported, so I will need to generate the weights, then apply them "manually."

need for surface elevation data

details

Use ESMF, NCL, and R (driven by bash scripts) to

  1. Find surface elevation data. (Went with ETOPO1.)
  2. "Cartesianize" the input grid.

    1. (R) Plot global, spherical/lon-lat MOZART-outputted N2O concentrations for eyeballing against future assimilations.
    2. (R) Crop raw input (global, spherical/lon-lat data on N2O concentrations, surface pressure (PS), and surface elevation (zS)) to the horizontal extents of the (output) AQMEII-NA domain (and plot for verification). Note that
      • only data variables={N2O, PS, zS} have horizontal coordinates
      • data variable=PS has only horizontal coordinates (not vertical)
      • HσP data variables={hyai, hyam, hybi, hybm} have only vertical coordinates (and no units)
      • data variable=P0 (reference pressure) has neither coordinates nor units
    3. (R) Regrid cropped surface elevation grid to match cropped {N2O, PS} grids.
    4. (NCL) Compute layer-interface elevations, producing new datavar=double z_int_geo(ilev, lat, lon) in 2008N2O_restart_region_z.nc using NCL::pres_hybrid_ccm and NCL::stdatmus_p2tdz. Replaced missing values output by NCL::stdatmus_p2tdz with geographically-corresponding values from ETOPO1:

      • if ETOPO1 value >= 0: replace missing value with ETOPO1 value, and add latter to all above-surface layer interfaces
      • if ETOPO1 value < 0 (mostly oceanic): replace missing value with 0
    5. (NCL) Compute layer-midpoint elevations, producing new datavar=double z_mid_geo(lev, lat, lon) in 2008N2O_restart_region_z.nc by taking the mean of the adjacent layer interfaces (and checking that they are vertically monotonic).

    6. (NCL, R) Create purely Cartesian (x,y,z) datavars for consumption by ESMF_RegridWeightGen. Caveats:

      1. Datavar for input (i.e., regrid from) completed (in this step=2), datavar for output (i.e., regrid to) in process (following step=3). Latter will require some code refactoring (or just plain dumb copy/mod).
  3. (in process) Cartesianize the output grid.

    1. (NCL) Copy one datavar (sufficient to define dimensions) from original output to pruned output.
    2. (R) Convert horizontal grid nodes from LCC ([hσp][hybrid sigma-pressure @ wikipedia],col,row) in pruned output to geographical (hσp,lat,lon) in geographicalized output.
    3. (planned) (NCL) Convert vertical layer-interface (i.e., grid node) values from HσP to elevations using NCL::pres_hybrid_ccm and NCL::stdatmus_p2tdz, replaced missing values output by NCL::stdatmus_p2tdz with geographically-corresponding values from ETOPO1:

      • if ETOPO1 value >= 0: replace missing value with ETOPO1 value, and add latter to all above-surface layer interfaces
      • if ETOPO1 value < 0 (mostly oceanic): replace missing value with 0
    4. (planned) Compute layer-midpoint elevations by taking the mean of the adjacent layer-interface elevations (and checking that they are vertically monotonic).

    5. (planned) (NCL, R) Create purely Cartesian (x,y,z) datavars for consumption by ESMF_RegridWeightGen.
    6. (planned) Use ESMF_RegridWeightGen to generate regrid weights from the Cartesianized input and output grids. Caveats:
      1. There is currently no code for generating the nodeCoords and elementConn datavars for the ESMF Unstructured Grid File Format. However I do have (on unfortunately private list esmf_support@list.woc.noaa.gov) pseudocode from Bob Oehmke for this.
    7. (planned) Use regrid weights to calculate and attribute values (concentrations) to output gridcells.

implementation

To run the examples:

  1. Clone this repo.
  2. cd to its working directory.
  3. Process MOZART-side data:

    1. Visualize the global data:
      1. Edit ./levelplot_netCDF_global.sh to match your local configuration (e.g., set your PDF viewer).
      2. Run ./levelplot_netCDF_global.sh to visualize the raw data.
      3. Check your resulting visualization against the reference (download).
    2. Create and visualize regional subsets of the global data:
      1. Edit ./restrict_lon_lat_to_region__MOZART.sh to match your local configuration.
      2. Run ./restrict_lon_lat_to_region__MOZART.sh to:
        1. crop the global raw input data to the region (and save separately).
        2. visualize the regional data.
      3. Check your resulting visualizations against provided references (downloads):
    3. Regrid the surface elevation data to match the grid of the other regional inputs:
      1. Edit ./regrid_surface_elevation__MOZART.sh to match your local configuration.
      2. Run ./regrid_surface_elevation__MOZART.sh to:
        1. regrid the surface elevation data (and save separately).
        2. visualize the regridded data.
      3. Check your resulting visualizations against the reference (download).
    4. Compute layer-interface and layer-midpoint elevations:

      1. Edit ./hsp_to_gh_amsl__MOZART.sh to match your local configuration.
      2. Run ./hsp_to_gh_amsl__MOZART.sh.
      3. Check the values of the resulting new elevation datavars (e.g.,

        double z_int_geo(ilev, lat, lon) # layer-interface elevations on geographical coordinates
        double z_mid_geo(lev, lat, lon)  # layer-midpoint elevations on geographical coordinates
        

      ) in the resulting netCDF output against the reference (download). (Yes, this needs visualizations ...) 1. Check the values of z_int_geo written (as ASCII CSV, for data interchange with R) to the reference (download). (Yes, this also needs visualizations ...)

    5. Compute fully-Cartesian layer-interface elevations:

      1. Edit ./cartesianize_elevations.sh to match your local configuration.
      2. Run ./cartesianize_elevations.sh. Tested on both tlrPanP5 and terrae.
      3. TODO: provide visualization to check the values of the resulting elevation datavar==nodeCoords
  4. Process CMAQ-side data:

<!-- Above see markdown for nesting code, and continuation of a (parent) list item after a nested sublist. Note that BB doesn't render markdown as well as Markdown.pl! E.g., this comment is visible in BB, but not in Markdown.pl-generated HTML :-( -->