Setting up grid for NorESM2-LM (vertical levels equal across variables)

[jdldeauna]

Hi! I'm trying to create a grid object for NorESM2-LM (grid label: gr) which has vertical levels equal across variables:

cat_url = "https://storage.googleapis.com/cmip6/pangeo-cmip6-noQC.json"
col = intake.open_esm_datastore(cat_url)
cat = col.search(table_id='Omon',experiment_id='historical',
             variable_id=['wmo','umo'], grid_label='gr', member_id='r1i1p1f1',source_id='NorESM2-LM')
noresm2_lm = cat.to_dataset_dict(
        zarr_kwargs={'consolidated':True, 'decode_times': True, 'use_cftime': True},
        preprocess=combined_preprocessing,
        aggregate=False)
wmo = noresm2_lm['CMIP.NCC.NorESM2-LM.historical.r1i1p1f1.Omon.wmo.gr.gs://cmip6/CMIP6/CMIP/NCC/NorESM2-LM/historical/r1i1p1f1/Omon/wmo/gr/v20190815/.nan.20190815.resolved.low.https://errata.es-doc.org/static/view.html?uid=5ebabff0-388f-07bc-b2cf-d44dcbb2940f']
umo = noresm2_lm[ 'CMIP.NCC.NorESM2-LM.historical.r1i1p1f1.Omon.umo.gr.gs://cmip6/CMIP6/CMIP/NCC/NorESM2-LM/historical/r1i1p1f1/Omon/umo/gr/v20190815/.nan.20190815.resolved.low.https://errata.es-doc.org/static/view.html?uid=5ebabff0-388f-07bc-b2cf-d44dcbb2940f']
diff = wmo.lev-umo.lev
diff.plot()

umo = umo.rename({'x':'x_g','lon':'lon_u', 'lat':'lat_u'})
ds_subset = xr.merge([umo, wmo], compat='override')
from xgcm import Grid
grid = Grid(ds_subset, coords={'X':{'center':'x', 'left':'x_g'}, 'Z':{'center':'lev'},},)
dw_dz = grid.diff(ds_subset.wmo, 'Z')

`KeyError: 'center'`

Is there a way to pre-format the grid through cmip6_pp before creating the grid object? cc: @jbusecke

[jbusecke]

I am looking into this right now. Quick comment: you are using the non-QC catalog. This is quite dangerous for 2 reasons:

1. If the dataset is in this catalog and not in the main one ("https://storage.googleapis.com/cmip6/pangeo-cmip6.json") there is probably something wrong with it.
2. Maybe more importantly, our retraction logic only applies filtering to the main catalog, so the data you work with might have been retracted!

I strongly recommend to only use "https://storage.googleapis.com/cmip6/pangeo-cmip6.json" for your work.

[jbusecke]

That said, your example still works with the newest catalog.

I slightly modified the example:

from cmip6_preprocessing.utils import google_cmip_col
from cmip6_preprocessing.preprocessing import combined_preprocessing
col =  google_cmip_col() # Just a little shortcut. This uses the most up to date catalog!
cat = col.search(table_id='Omon',experiment_id='historical',
             variable_id=['wmo','umo'], grid_label='gr', member_id='r1i1p1f1',source_id='NorESM2-LM')
noresm2_lm = cat.to_dataset_dict(
        zarr_kwargs={'consolidated':True, 'decode_times': True, 'use_cftime': True},
        preprocess=combined_preprocessing,
        aggregate=False)
wmo = noresm2_lm['CMIP.NCC.NorESM2-LM.historical.r1i1p1f1.Omon.wmo.gr.gs://cmip6/CMIP6/CMIP/NCC/NorESM2-LM/historical/r1i1p1f1/Omon/wmo/gr/v20190815/.nan.20190815']
umo = noresm2_lm[ 'CMIP.NCC.NorESM2-LM.historical.r1i1p1f1.Omon.umo.gr.gs://cmip6/CMIP6/CMIP/NCC/NorESM2-LM/historical/r1i1p1f1/Omon/umo/gr/v20190815/.nan.20190815']
diff = wmo.lev-umo.lev
diff.plot()

and it gives the same result!

[jbusecke]

Let me see if @AleksiNummelin knows how to proceeed here.

Hey Aleksi, to quickly summarize: @jdldeauna is trying to compute tracer budgets from the NorESM models, but for the regridded (grid_label='gr') output, the vertical transport is on the tracer points. Am I assuming rightfully, that both the horizontal/vertical transports have been interpolated onto the tracer point, and as such cannot be used to get exact budgets for the tracer cell (since we usually want the flux terms on the surrounding cell faces)?

[AleksiNummelin]

Hi Julius, I think we actually had a bit of a email exchange with @jdldeauna earlier this month around some of these issues, but I think this might be a bit distinct. Anyway, here's a bit of a re-iteration of some of the points from the email conversation. I wanted to include also @matsbn here just to double check as he knows the code better than I do (Mats, please correct if I'm talking nonsense here).

NorESM ocean model (https://github.com/NorESMhub/BLOM) works on a potential density coordinate so the 'lev' grid (z-coordinate) is arbitrary and indeed interpolated. wmo is in fact not a model variable, but diagnosed from (accumulated) horizontal transports for a given layer (code starting at https://github.com/NorESMhub/BLOM/blob/master/phy/mod_dia.F#L3863). Now, although the horizontal (umo,vmo) transports are defined at a given level, they are flux variables and depth=0 m is actually the flux between 0-2.5 m depth. Now wmo has the same lev coordinate as umo, vmo (but not lon,lat), but if I interpret the code correctly, it is actually the top interface flux i.e. wmo at depth=0 is really the flux trough depth=0 interface, whereas umo, vmo at depth=0 are really the fluxes trough the edges of a 2.5 m thick cell. Therefore, I think it should still be possible to close the budgets. I had a quick try and (sorry, just pseudo code here, but I think you get the idea)

conservation_error = - wmo.diff('lev').isel(lev=0) - (umo.isel(i=slice(1,nx),j=slice(0,-1))-umo.isel(i=slice(0,-1),j=slice(0,-1))).isel(lev=0) - (vmo.isel(j=slice(1,ny),i=slice(0,-1))-vmo.isel(j=slice(0,-1),i=slice(0,-1))).isel(lev=0)

is still considerable, but probably on the order I would expect and the global sum is tiny. Obviosuly, one should wrap around longitude properly and assume flux is 0 trough the bottom.

Hope this helps. I realize that this might not fit right into the cmip6_preprocessing workflow, but probably it would be easiest to redefine the umo,vmo depth axes.

[jbusecke]

Thanks for that thorough answer @AleksiNummelin.

Hope this helps. I realize that this might not fit right into the cmip6_preprocessing workflow, but probably it would be easiest to redefine the umo,vmo depth axes.

If I understand correctly I could try to add an additional layer of 0's to the bottom of wmo and then treat them as the 'outer' fluxes? Ill have a go at this in a bit and report back.

[jdldeauna]

Thank you for the discussion on this issue, and sorry for the late reply!

I strongly recommend to only use "https://storage.googleapis.com/cmip6/pangeo-cmip6.json" for your work.

I use the no-QC version when I need to use MPI-ESM1-2-HR salinity since it has an errata. It's for locations in the small bays of the Baltic, so I figured it would be fine for my purposes in the North Pacific. But I'll def try to use the no-errata version as much as I can.

If I understand correctly I could try to add an additional layer of 0's to the bottom of wmo and then treat them as the 'outer' fluxes?

Just to clarify, when would it count as a variable is shifted to the left vs when it should be treated as outer?