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

ocean depth fix for cnrm_esm2_1, gfdl_esm4, ipsl_cm6a_lr datasets + mcm_ua_1_0 #1098

Merged
merged 14 commits into from
May 10, 2021

Conversation

tomaslovato
Copy link
Contributor

@tomaslovato tomaslovato commented May 5, 2021

Based on PR #1095 create the ocean depth mip-level fixes for other problematic CMIP6 datastest, namely cnrm_esm2_1, gfdl_esm4, ipsl_cm6a_lr

Closes #1069 (and also #527 and #1100)

Link to documentation:


Checklist

It is the responsibility of the author to make sure the pull request is ready to review. The icons indicate whether the item will be subject to the 🛠 Technical or 🧪 Scientific review.


To help with the number pull requests:

@tomaslovato tomaslovato added the fix for dataset Related to dataset-specific fix files label May 5, 2021
@tomaslovato tomaslovato requested a review from jvegreg as a code owner May 5, 2021 13:10
@tomaslovato tomaslovato self-assigned this May 5, 2021
@tomaslovato tomaslovato requested review from schlunma and removed request for jvegreg May 5, 2021 13:10
@tomaslovato
Copy link
Contributor Author

@schlunma could you please give me a feedback on test_cnrm_esm2_1.py before I continue with the other datasets (I feel like something is missing)... Thanks!

@schlunma
Copy link
Contributor

schlunma commented May 5, 2021

Looks good already! The only thing that's missing is a test of the actual fix, so something like this:

@pytest.fixture
def thetao_cubes():
"""Cubes to test fixes for ``thetao``."""
time_coord = iris.coords.DimCoord(
[0.0004, 1.09776], var_name='time', standard_name='time',
units='days since 1850-01-01 00:00:00')
lat_coord = iris.coords.DimCoord(
[0.0, 1.0], var_name='lat', standard_name='latitude', units='degrees')
lon_coord = iris.coords.DimCoord(
[0.0, 1.0], var_name='lon', standard_name='longitude', units='degrees')
lev_coord = iris.coords.DimCoord(
[500.0, 1000.0], bounds=[[2.5, 7.5], [7.5, 12.5]],
var_name='lev', standard_name=None, units='cm',
attributes={'positive': 'up'})
coord_specs = [
(time_coord, 0),
(lev_coord, 1),
(lat_coord, 2),
(lon_coord, 3),
]
thetao_cube = iris.cube.Cube(
np.ones((2, 2, 2, 2)),
var_name='thetao',
dim_coords_and_dims=coord_specs,
)
return iris.cube.CubeList([thetao_cube])

and this

def test_thetao_fix_metadata(thetao_cubes):
"""Test ``fix_metadata`` for ``thetao``."""
vardef = get_var_info('CMIP6', 'Omon', 'thetao')
fix = Omon(vardef)
out_cubes = fix.fix_metadata(thetao_cubes)
assert out_cubes is thetao_cubes
assert len(out_cubes) == 1
out_cube = out_cubes[0]
# Check metadata of depth coordinate
depth_coord = out_cube.coord('depth')
assert depth_coord.standard_name == 'depth'
assert depth_coord.var_name == 'lev'
assert depth_coord.long_name == 'ocean depth coordinate'
assert depth_coord.units == 'm'
assert depth_coord.attributes == {'positive': 'down'}

I think you can just copy-pasting these lines should work!

@tomaslovato tomaslovato added this to the v2.3.0 milestone May 5, 2021
@tomaslovato
Copy link
Contributor Author

@schlunma would it make this PR too long if I include also changes for MCM-Ua-1-0 dataset (#527)?

@schlunma
Copy link
Contributor

schlunma commented May 6, 2021

No, that should be fine! Just let me know when it's ready for review 👍

@tomaslovato tomaslovato changed the title ocean depth fix for cnrm_esm2_1, gfdl_esm4, ipsl_cm6a_lr datasets ocean depth fix for cnrm_esm2_1, gfdl_esm4, ipsl_cm6a_lr datasets + mcm_ua_1_0 (#527) May 6, 2021
@tomaslovato tomaslovato changed the title ocean depth fix for cnrm_esm2_1, gfdl_esm4, ipsl_cm6a_lr datasets + mcm_ua_1_0 (#527) ocean depth fix for cnrm_esm2_1, gfdl_esm4, ipsl_cm6a_lr datasets + mcm_ua_1_0 May 6, 2021
@remi-kazeroni remi-kazeroni self-requested a review May 6, 2021 09:10
@remi-kazeroni remi-kazeroni added the AR6 Necessary changes for IPCC AR6 label May 6, 2021
@tomaslovato
Copy link
Contributor Author

tomaslovato commented May 6, 2021

@schlunma @remi-kazeroni the code is ready for review 🚀 .
I already made a test with a recipe to check that it works for ocean variables thetao and no3.
Note that for mcm_ua_1_0 I decided to integrate needed changes in the original code under the Allvars fix.

@schlunma
Copy link
Contributor

schlunma commented May 6, 2021

Great stuff! Having a look at it right now!

@remi-kazeroni
Copy link
Contributor

@schlunma @remi-kazeroni the code is ready for review 🚀 .
I already made a test with a recipe to check that it works for ocean variables thetao and no3.
Note that for mcm_ua_1_0 I decided to integrate needed changes in the original code under the Allvars fix.

Thanks a lot! I'll have a look later today, in particular for the MCM_UA_1_0 fixes.

Copy link
Contributor

@schlunma schlunma left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice work @tomaslovato, looks really good already!

I have some suggestions to the fix of MCM and its test.

@tomaslovato
Copy link
Contributor Author

@schlunma I commited all the review changes in d9988fb

@tomaslovato
Copy link
Contributor Author

@schlunma I got a bit confused about test_mcm_ua_1_0.py (I'm not skilled enough on this) as it is by now it only goes through Omon without applying changes from AllVars to roll the cube and coordinates... any suggestion?

@schlunma
Copy link
Contributor

schlunma commented May 6, 2021

No worries! There are two options here:

  1. Split the test function into two tests (maybe test_omon_thetao_fix_metadata and test_allvars_thetao_fix_metadata), use fix = Omon(...) and fix = AllVars(...), respectively, and only inlcude the asserts that belong to the respective fix.
  2. Apply both fixes in the current test function:
    fix_omon = Omon(vardef)
    fix_allvars = AllVars(vardef)
    out_cubes = fix_omon.fix_metadata(thetao_cubes)
    out_cubes = fix_allvars.fix_metadata(out_cubes)

Option 2 is probably easier to implement, option 1 clearly seperates the two individual fixes. Optimally we would have both tests (i.e. option 1 and 2), but I'm fine with whatever you like to implement.

@tomaslovato
Copy link
Contributor Author

    fix_omon = Omon(vardef)
    fix_allvars = AllVars(vardef)
    out_cubes = fix_omon.fix_metadata(thetao_cubes)
    out_cubes = fix_allvars.fix_metadata(out_cubes)

I think that for the time being option 2 will be ok, as it applies the fix to thetao including both Omon and AllVars as it should do when calling it from a recipe ...

@schlunma
Copy link
Contributor

schlunma commented May 6, 2021

Tests are passing now, thank you really much for the fixes @tomaslovato! 🎉

Would you mind if I push 1-2 commits to this branch directly? I just realized that a test for Fgco2 of GFDL is missing (which is also missing for CESM2 and which I wanted to add for a long time) and that I forgot to add some documentation about the mip-level fixes.

I will approve afterwards and merge when @remi-kazeroni had the chance to test the fixes. Cheers!

@tomaslovato
Copy link
Contributor Author

@schlunma thank you for the nice interaction! About you addition just go for it !
I've a little request ... could you have a look to #1097, that is only one line code change leftover from the former objective of this issue.

Copy link
Contributor

@schlunma schlunma left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@remi-kazeroni Please merge when you tested successfully 👍

@tomaslovato
Copy link
Contributor Author

While working on the time coordinate issue for climate_statistic I found a problem with the dataset fix for ipsl_cm6a_lr.

The cell_area variable was oddly added to almost all ocean variables of this dataset, so a specific fix was setup about 2 years ago

try:
cell_area = cubes.extract_cube('cell_area')
except ConstraintMismatchError:
return cubes
cell_area = AuxCoord(
cell_area.data,
standard_name=cell_area.standard_name,
long_name=cell_area.long_name,
var_name=cell_area.var_name,
units=cell_area.units,
)
new_list = CubeList()
for cube in cubes:
if cube.name() == 'cell_area':
continue
cube.add_aux_coord(cell_area, cube.coord_dims('latitude'))
cube.coord('latitude').var_name = 'lat'
cube.coord('longitude').var_name = 'lon'
new_list.append(cube)
return CubeList(new_list)

After the introduction of cube cell_measures trough _ancillary_vars.py, the ingestion of areacello from this dataset fails with an empty cubelist error.
Looking at the dump of the areacello file (here below), You will see that there are two variables (cell_area, areacello) both having as standard_name cell_area which is used by iris to assign cube name.
The general fix control

if cube.name() == 'cell_area':
continue

is responsible for trashing both variables and produce an empty cubelist.

For me the simplest solution is to modify the if statement as in the following
if cube.name() == 'cell_area' and cube.var_name != 'areacello':

but probably we should also think to remove the creation of the auxcoord 'cell_area' (likely kept in there at the time for some potential use) which is a potentially misleading duplicate of the cube cell_measure property.

areacello ipsl_cm6a_lr
netcdf areacello_Ofx_IPSL-CM6A-LR_historical_r1i1p1f1_gn {
dimensions:
	axis_nbounds = 2 ;
	x = 362 ;
	y = 332 ;
	nvertex = 4 ;
	time = UNLIMITED ; // (0 currently)
variables:
	float nav_lat(y, x) ;
		nav_lat:standard_name = "latitude" ;
		nav_lat:long_name = "Latitude" ;
		nav_lat:units = "degrees_north" ;
		nav_lat:bounds = "bounds_nav_lat" ;
	float nav_lon(y, x) ;
		nav_lon:standard_name = "longitude" ;
		nav_lon:long_name = "Longitude" ;
		nav_lon:units = "degrees_east" ;
		nav_lon:bounds = "bounds_nav_lon" ;
	float bounds_nav_lon(y, x, nvertex) ;
	float bounds_nav_lat(y, x, nvertex) ;
	float area(y, x) ;
		area:standard_name = "cell_area" ;
		area:units = "m2" ;
		area:coordinates = "nav_lon nav_lat" ;
	float areacello(y, x) ;
		areacello:standard_name = "cell_area" ;
		areacello:long_name = "Grid-Cell Area" ;
		areacello:units = "m2" ;
		areacello:online_operation = "once" ;
		areacello:cell_methods = "area: sum" ;
		areacello:cell_measures = "area: area" ;
		areacello:_FillValue = 1.e+20f ;
		areacello:missing_value = 1.e+20f ;
		areacello:coordinates = "nav_lat nav_lon" ;
		areacello:description = "Cell areas for any grid used to report ocean variables and variables which are requested as used on the model ocean grid (e.g. hfsso, which is a downward heat flux from the atmosphere interpolated onto the ocean grid). These cell areas should be defined to enable exact calculation of global integrals (e.g., of vertical fluxes of energy at the surface and top of the atmosphere)." ;
		areacello:history = "none" ;

// global attributes:
		:Conventions = "CF-1.7 CMIP-6.2" ;
		:creation_date = "2018-07-11T07:26:30Z" ;
		:tracking_id = "hdl:21.14100/718e5096-9a8f-4794-a867-0ac583371b62" ;
		:description = "CMIP6 historical" ;
		:title = "IPSL-CM6A-LR model output prepared for CMIP6 / CMIP historical" ;
		:activity_id = "CMIP" ;
...

@schlunma
Copy link
Contributor

schlunma commented May 7, 2021

Good find @tomaslovato! I just pushed a tiny fix for this, please let me know if it works for you! As you said, we do not need the cell_area coordinate since this is included via cell_measures.

@tomaslovato
Copy link
Contributor Author

@schlunma I re-run the recipe with the modified code and it works great!

@schlunma
Copy link
Contributor

schlunma commented May 7, 2021

Great! Ready to merge from my side

@zklaus
Copy link

zklaus commented May 10, 2021

@remi-kazeroni, could you please have a quick look?

@remi-kazeroni
Copy link
Contributor

@remi-kazeroni, could you please have a quick look?

Yes, I'm on it. I'm done opening issues and PRs for the Core. It's time to start closing them 🍺

Copy link
Contributor

@remi-kazeroni remi-kazeroni left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Everything looks and works fine! Thanks a lot for your work @tomaslovato 👍 This can be merged

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
AR6 Necessary changes for IPCC AR6 fix for dataset Related to dataset-specific fix files
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Apply ocean depth fixes for cnrm_esm2_1, gfdl_esm4, ipsl_cm6a_lr
4 participants