Regridding Cases and Current Source Code

Note all these codes can handle quantities on either center or edges on the sigma levels but since the actual vertical coordinate is defined by the edge pressure on the sigma levels you need PLE for any hybrid sigma coordinate.

Model Levels to Model Levels

We have at least 3 versions of the same code to do this, but I think Bill has said in the past the on in GMAO_Shared/GEOS_Shared is probably the one we should look at. It's what MKIAU uses for example to do the vertical remapping if needed.

SJ Code in GMAO_Shared/Hermes

SJ Code in GMAO_Shared/GEOS_Shared

SJ Code in FV3

Model Levels to Fixed Pressure (Heights)

Code used by History to do model to pressure level regridding.

Fixed Pressure (Heights) to Model Levels

Module Jules wrote to perform pressure level to GEOS model level interpolation.

File Metadata Proposal or What We are Looking For

For this I will attempt to follow what the CF convention says we should do and be putting in the files to correctly identify the vertical information.

First we must find the coordinate variable (i.e. a variable with a name matching a dimension) that corresponds to the vertical dimension. According to the CF conventions, such a variable will be either have units of pressure (height) and/or the presents of the positive attribute. Therefore the first step is to identify this variable if it exists. Once we have found it there are two possibilities. Note that for the following it is worth it to review what CF says about dimensional vs dimensionless variables

Pressure(Height) Coordinate

The first option is if the vertical coordinate presents a dimensional pressure (height). In that case the units will clearly be a units of vertical height in units of length or units of pressure as defined by UDUNITS. Since this has no dependence on the horizontal coordinate this is a single 1D array to read which can be used for each column to regrid to model levels. An example level variable is given below:

float lev(lev) ;
  lev:units = "hPa" ;
  lev:long_name = "air _pressure" ;
  lev:positive = "down" ;
  lev:standard_name = "air_pressure" ;

GEOS Model levels

The second option is the vertical coordinate represents a dimensionless unit which define a parametric vertical coordinate.. As specified in the CF convention dimensional variables may have no units attribute. In this case the presence of the positive attribute will be the give away the this is the vertical coordinate. Note that the CF convention says the user may optionally give a unit of level, layer, or sigma_level for backwards compliance with COARDS but this is not part of the CF standard so we should not rely on this.

Once we have identified the vertical coordinate as a dimensionless vertical coordinate we need to look for further information. CF defines the GEOS model levels as a parametric coordinate. In other words the actual dimensional vertical coordinate cannot be store as as 1D array but depends on the horizontal location. The standard_name vertical coordinate is used to find the definition of the associated dimensional vertical coordinate. The formula_terms attribute defines what other variables in the file are needed to compute the dimensional coordinate. GEOS uses hybrid sigma and the correct standard_name and formula_terms can be found in the link but are repeated here:

float lev(lev) ;
  lev:long_name = "???" ;
  lev:positive = "down" ;
  lev:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
  lev:formula_terms = ""a: var1 b: var2 ps: var3 p0: var4"" ;
  lev:computed_standard_name = "air_pressure" ;

where

p(t,k,j,i) = a(k)*p0 + b(k)*ps(t,j,i)

Note that we fold p0 into a(k) so NO need need for the file to actually define p0.

Given all this information we can determine the input file is on GEOS model levels and these can be computed via the formula with information in the file. So the file needs PS, ak, and bk which we are definitely not putting in files now. Note that once could also have PLE in the file but I would think for space reasons this would be undesireable.

Problems

Note that our current History output nor the files prepared and ingested by ExtData in many cases do not follow this standard. For example we some input data for ExtData that is on model levels with units of pressure.

Also the current AK/BK values, rather than in say a simple ASCII file we keep or YAML file, are hard coded in source code in GMAO_Shared so not in an easy place to get at them for people creating new datasets.

Nor do I think dynamics exports them if we wanted to write them say in History output, but this is easily fixable.

What Do We Currently Do/Support

History

Note that our output not CF compliant in cases. Output on GEOS model levels:

double lev(lev) ;
       lev:long_name = "vertical level" ;
       lev:units = "layer" ;
       lev:positive = "down" ;
       lev:coordinate = "eta" ;
       lev:standard_name = "model_layers" ;

and output on fixed pressure levels:

double lev(lev) ;
       lev:long_name = "vertical level" ;
       lev:units = "hPa" ;
       lev:positive = "down" ;
       lev:coordinate = "PLE" ;
       lev:standard_name = "PLE_level" ;

ExtData

Currently ExtData uses the following algorithm to find the coordinate variable for the vertical. This is not following what CF says to do.

1. It first checks the there a variable named lev or height.
2. If this is not found, check the units of each variable until it finds on that matches hPa, mb, sigma_level, or millibar.

What Should Be in the Files for Center/Edge Variables.

Now to the real world. What we do know files will have. If we LM model levels, a variable on a file on level centers, and enough information to vertically regrid this then:

1. The file will have to have a dimension of size LM with an associated coordinate variable
2. The file will have to have a dimension of size LM+1 with an associated coordinate variable.
3. Reading the earlier section it seems the coordinate variable of size LM+1 will have standard name of atmosphere_hybrid_sigma_pressure_coordinate with the formula term attribute to describe how to compute PLE from other variables in the file (PS, AK, BK, PTOP). But there is this section which throws this into question.

So the questions

1. What should the names of these variables be? CF is silent on this.
2. What should the standard name of be coordinate variable of size LM so we can say ah ah, this is the centers of the parametric coordinate heights, my first thought is that I see a CF standard name that is model_level_number so that's this should be, but then...
3. There is this section, that is perhaps what we should be doing for any of cell boundaries including the vertical, but in practice we would not want to do this for any multidimensional coordinate...

Proposals for vertical metadata

Here are two options that pretty much follow CF. I'm putting a sample ncdump of a hypothetical 72 level file for each. I've included data variables on both the center and edge to reflect the most complex use case. I've omitted the horizontal coordinates.

Proposal 1

This is one proposal, you could say, we have not violated anything that CF says, in fact this what I got the from first reading. One could say the lev represents an alternative coordinate.

dimensions:
        lon = 1440 ;
        lat = 721 ;
        lev = 72 ;
        eta = 73 ;
        time = UNLIMITED ; // (1 currently)

float lev(lev) ;
  lev:long_name = "model_level_number" ;
  lev:standard_name = "model_level_number" ;

float eta(eta) ;
  eta:long_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
  eta:positive = "down" ;
  eta:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
  eta:formula_terms = ""ak: var1 bk: var2 PS: var3 p0: var4"" ;
  eta:computed_standard_name = "air_pressure" ;


float O3(time, lev, lat, lon) ;
  O3:long_name = "ozone_mass_mixing_ratio" ;
  O3:units = "kg kg-1" ;
  O3:standard_name = "ozone_mass_mixing_ratio" ;

float ZLE(time, eta, lat, lon) ;
  ZLE:long_name = "height" ;
  ZLE:units = "m" ;
  ZLE:standard_name = "height" ;

float PS(time, lat, lon) ;
  PS:standard_name = "surface_air_pressure"
  PS:units = "Pa"

float ak(edge) ;
float bk(edge) ;

Or if we follow 7.1 of CF Standard

This section says every coordinate variable should have a bounds attribute (which of course we should be doing for horizontal coordinates supposedly), which contains the vertices. A literally reading says we should be doing something like this. But then the problem, you still need an LM+1 size dimension and coordinate variable, you certainly wouldn't want to store 3D data variables using the extra dimension for the bounds if you have variables on the level edges.. What the heck do you call this (it is `XXX`` in the example below). You still need to figure out, ok, this is the edges, that are defined by bounds of some other coordinate variable.

The standard is just lacking here. Also their recommendations for storing bounds of say 2D grids while very general results in wasteful store. They acknowledge this, but then say, it doesn't matter, but it does in the real world once we are not talking about 1D arrays. But here is what they say we should, notice all the question marks. for the edge coordinate variable.

Also, the actual ak and bk that we think of would be stored in the AK_bnds and BK_bnds arrays with duplication as they say in 7.1. Then the arrays AK and BK would have to be defined as AK(i)=(AK_bnds(i,2)+AK_bnds(i,1))/2 and BK(i)=(BK_bnds(i,2)+BK_bnds(i,1))/2 which are NOT what the average GEOS user would think of as AK and BK

dimensions:
        lon = 1440 ;
        lat = 721 ;
        eta_centers = 72 ;
        eta_edges = 73 ;
        time = UNLIMITED ; // (1 currently)

float eta_centers(eta) ;
  eta_centers:long_name = "eta full levels" ;
  eta_centers:positive = "down" ;
  eta_centers:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ;
  eta_centers:formula_terms = "ak: AK bk: BK ps: PS p0: P04" ;
  eta_centers:computed_standard_name = "air_pressure" ;
  eta_centers:bounds = "eta_bnds"

float eta_bnds(eta, 2) ;
  eta_bnds:formula_terms = "ak: AK_bnds bk: BK_bnds ps: PS p0: P0" ;

float eta_edges(eta_edges) ;
  eta_edges:long_name = ??? ;
  eta_edges:positive = ??? ;
  eta_edges:standard_name = ??? ;

float O3(time, eta, lat, lon) ;
  O3:long_name = "ozone_mass_mixing_ratio" ;
  O3:units = "kg kg-1" ;
  O3:standard_name = "ozone_mass_mixing_ratio" ;

float ZLE(time, XXX, lat, lon) ;
  ZLE:long_name = "height" ;
  ZLE:units = "m" ;
  ZLE:standard_name = "height" ;

float PS(time, lat, lon) ;
  PS:standard_name = "surface_air_pressure"
  PS:units = "Pa"

float AK_centers(eta) ;
float bk_centers(eta) ;

float AK_bnds(eta,2) ;
float BK_bnds(eta,2) ;

Fixed Pressure Levels with Bounds

Of course what if you have fixed pressure levels and want to define edges or bounds. Well same problem, you could follow 7.1 and use the bounds attribute.

float lev(lev) ;
  lev:units = "hPa" ;
  lev:long_name = "air _pressure" ;
  lev:positive = "down" ;
  lev:standard_name = "air_pressure" ;
  lev:bounds = "lev_bnds"

float lev_bnds(lev, 2)