Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Met_AREAM2 grid cell area #338

Closed
cbutenhoff opened this issue Aug 15, 2023 · 8 comments · Fixed by geoschem/geos-chem#1931
Closed

Met_AREAM2 grid cell area #338

cbutenhoff opened this issue Aug 15, 2023 · 8 comments · Fixed by geoschem/geos-chem#1931
Assignees
Milestone

Comments

@cbutenhoff
Copy link

Name and Institution (Required)

Name: Chris Butenhoff
Institution: Portland State University

Confirm you have reviewed the following documentation

Description of your issue or question

We are running GCHP at C24 and outputting diagnostics, including Met_AREAM2, which contains the grid cell areas, on a 2x2.5 lat-lon grid. When I sum over all the grid cell areas I would expect the total to be about Earth's surface area, 5.2e14 m^2 but instead I get 2.17e15 m^2, about 4 times larger.

Perhaps I'm not correctly understanding Met_AREAM2, or is there an issue with using a non-native grid for the output, or in outputting at 2x2.5 when the native resolution is closer to 4x5?

Thank you.

@sdeastham
Copy link
Contributor

Thanks @cbutenhoff ! Using a non-native grid for output is fine. Met_AREAM2 is an unusual quantity because it is extensive - it varies depending on the grid cell size. However, the regridding routines are area-conserving, and assume that the quantities being regridded are intensive. For example, to regrid mixing ratio (an intensive quantity - it is independent of grid cell size) from 1x1 to 2x2, you would take the area-weighted average of each set of four 1x1 grid cells to get the mean value within a single 2x2 grid cell. That is fundamentally what GCHP (really MAPL's HISTORY component) is doing. This generally works fine because almost all the quantities output from GEOS-Chem are intensive. However this gives you the wrong answer when the quantity in question is extensive, because what you really need is the sum, not the average. The key examples of this are grid cell area (Met_AREAM2) and dry air mass (Met_AD).

To cut a long story short: to get the grid cell area for a 2x2.5 grid I would rely on calculating it externally. The algorithm for a rectilinear grid is pretty simple - there's an actually an implementation in the GCPy grid module, here: https://github.com/geoschem/gcpy/blob/84d11cb1769e1fff938d1256cb2b90d84b75edfd/gcpy/grid.py#L903-L948

@cbutenhoff
Copy link
Author

Thanks much @sdeastham for your detailed explanation and pointing me to the GCPy routine. I will give that a try. Maybe this issue with MET_AREAM2 doesn't pop up much, but it might be worthwhile to make a note of this in the HISTORY.rc file.

@cbutenhoff
Copy link
Author

One other question.

When GCHP or MAPL regrids say to a 2x2.5 grid for output, are the regridded values calculated on the boundaries of grid cells or are they grid cell averages? I mention this because I noticed the size of the regridded 2x2.5 'lat' variable is 92 and the latitudes are specified on the grid cell boundaries, i.e. they run from -90.0, ..., +90.0. This would imply there are 91 cells along the latitude dimension. Other variables such as 'SpecConc' that hold the mixing ratios also have a latitude dimension of size 92 which must mean they too are calculated on the cell boundary, rather than grid cell averages, in which case they should have a latitude dimension of size 91.

@sdeastham
Copy link
Contributor

Whether you are using conservative or bilinear regridding in the HISTORY specification, the values should correspond to the grid cell average (or center, in the case of bilinear). There is one nuance though: the standard 2x2.5 grid cell centers do indeed start and end at +/-90 degrees because the grid is "pole-centered". That's what the "PC" means in the example of a 1x1 grid given in the standard HISTORY.rc file. To help debug, could you post the grid definition you're using under GRID_LABELS in HISTORY.rc for the 2x2.5 grid? I'll also mention that @lizziel may well have better insights on this as I rarely use the in-built regridding function.

@cbutenhoff
Copy link
Author

cbutenhoff commented Aug 16, 2023 via email

@sdeastham
Copy link
Contributor

Ah, got it! The pole definition is indeed pretty odd. If you just change GLOBAL2x25.JM_WORLD to be 91 instead of 92, you should find that the grid has the right dimensions (and corresponds to the 2x2.5 grid you expect, just with the rather funky definition of the "centers" for the northernmost and southernmost latitude bands).

@lizziel
Copy link
Contributor

lizziel commented Aug 17, 2023

Thanks for this discussion. I am going to update the docs with a caveat about outputting to non-native grids. This issue has not been brought up before and is an important one to highlight. Let's leave the issue open until I follow-up with docs updates so I don't forget.

@lizziel
Copy link
Contributor

lizziel commented Sep 1, 2023

Comments on this added in geoschem/geos-chem#1931. I also updated ReadTheDocs with a warning about this on the HISTORY.rc page.

@lizziel lizziel closed this as completed Sep 1, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants