xmitgcm icon indicating copy to clipboard operation
xmitgcm copied to clipboard

Use with CS grid?

Open brian-rose opened this issue 7 years ago • 12 comments

Does xmitgcm support reading output on the Cubed Sphere? I suspect the answer is no since I didn't see any thing in the docs about this.

We are starting up some new work with a coupled CS grid MITgcm setup, and trying to figure out the best python-based analysis workflow.

brian-rose avatar Aug 23 '18 16:08 brian-rose

Right now it only supports rectangular, LLC and ASTE grids.

But it should be pretty easy to extend it to work with CS, which is actually simpler than those other two. All the hard parts of this are done. @raphaeldussin is currently working on refactoring the internals (see #87, #89, #90).

Do you, or a someone in your group, have time to devote to this? If so, it can be done quickly.

rabernat avatar Aug 23 '18 17:08 rabernat

I might be able to put a bit of time into this, if you can hold my hand :-) Would this be a valuable contribution? I'm not sure how widespread the use of CS is these days.

Weighing this against using mnc output and existing (non-xarray) python scripts for our project. Or [shudder] going back to matlab.

brian-rose avatar Aug 23 '18 17:08 brian-rose

The new internals allow for an arbitrary numbers of faces and padding if necessary. I have looked at the documentation here cube sphere and I think it might work out of the box. Do you have a test dataset ? or is it close enough from one of the test cases (MITgcm/verification) ?

raphaeldussin avatar Aug 23 '18 18:08 raphaeldussin

My case is close to MITgcm/verification/cpl_aim+ocn

brian-rose avatar Aug 23 '18 20:08 brian-rose

I tested the ocean part (global_ocean.cs32x15) with my latest code and it almost worked. The only issue is that the llc option (that allows the faces) assumes nx = ny for each face. Once we got all those PR in, it should take a few lines to get the cube sphere working.

raphaeldussin avatar Aug 24 '18 18:08 raphaeldussin

Because the numpy / xarray / netCDF data models do not permit "ragged arrays", it is necessary for all "faces" in an xarray representation of these grids to have the same shape. This is not how MITgcm represents its "tiles" internally. For example, for the LLC grids, there are only five logical tiles, with the following shapes and orientations:

tile 0   tile 1   tile 2
(270,90) (270,90) (90, 90)

######   ######   ######
######   ######   ######
######   ######   ######
######   ######
######   ######
######   ######
######   ######
######   ######
######   ######

tile 3 (90, 270)
##############################
##############################
##############################

tile 4 (90, 270)
##############################
##############################
##############################

In order to fit this into a single xarray dataset, I chose to break up the tiles into their lowest common denominator. These are called faces. Because the smallest tile in the LLC grid is square (the arctic tile, 90x90), this constrains the geometry of the other faces. So there are 13 90x90 "faces" in an xmitgcm llc dataset:

image

My understanding is that the only difference between the LLC grid and the CS grid is that the CS grid contains a tile for the south pole, upping the number of faces to 14.

The only issue is that the llc option (that allows the faces) assumes nx = ny for each face

I can't see how to get around this. What is the original shape of the CS tiles?

rabernat avatar Aug 24 '18 18:08 rabernat

I think the CS grid is simpler. It's just 6 square tiles, i.e. nx = ny for each tile already. E.g. the horizontal grid for MITgcm/verification/global_ocean.cs32x15 is 6x32x32

brian-rose avatar Aug 24 '18 19:08 brian-rose

Good news! That maps much more cleanly to the xarray data model. So "assuming nx = ny" for each face should definitely not cause a problem..

rabernat avatar Aug 24 '18 19:08 rabernat

I got confused with the number of tiles in the test case (12 tiles of 16x32)... I can make xmitgcm read the data but it doesn't look right. Is it possible that they write into mds files in a different order ? (face, z, y, x ) instead of the llc (z, face, y, x) ?

raphaeldussin avatar Aug 24 '18 20:08 raphaeldussin

The answer might be buried in here: https://github.com/MITgcm/MITgcm/blob/master/pkg/exch2/w2_set_cs6_facets.F

rabernat avatar Aug 24 '18 20:08 rabernat

Hi @brian-rose ,

sorry for the delay! I have been working on this today. The cube sphere output goes like this: (z,y,face,x). Although it is not straightforward to make it work with dask (lazy evaluation of data) in the current implementation of xmitgcm, I have found a way to make it work eagerly. We're in the process of refactoring the low level functions but I can send a new PR with CS support once we're done.

raphaeldussin avatar Sep 11 '18 16:09 raphaeldussin

Great! For our purposes we don't need lazy evaluation since the grids are small. But let me know I can help with anything.

brian-rose avatar Sep 11 '18 19:09 brian-rose