Skip to content

Conversation

navidcy
Copy link
Collaborator

@navidcy navidcy commented Apr 23, 2024

Copy link

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

@navidcy
Copy link
Collaborator Author

navidcy commented Apr 23, 2024

We need to ensure the example runs.

Also we need to modify so that it uses the cosima-cookbook to load ACCESS-OM2 data; see COSIMA/regional-mom6#131

@navidcy
Copy link
Collaborator Author

navidcy commented Apr 25, 2024

@ashjbarnes I tried to run this and it errors on initial condition. See the output of Step 5 cell at https://app.reviewnb.com/COSIMA/cosima-recipes/pull/338/

@navidcy
Copy link
Collaborator Author

navidcy commented May 3, 2024

@ashjbarnes is this error related to the bug in COSIMA/regional-mom6#170? If so, it means that https://github.com/COSIMA/regional-mom6/releases/tag/v0.5.1 would be required.

cc @luweiyang

@navidcy
Copy link
Collaborator Author

navidcy commented May 6, 2024

I tried to see if regional-mom6 v0.5.2 healed the issues I was having. It didn't... :(

@luweiyang
Copy link
Collaborator

@ashjbarnes I tried to run the ACCESS-OM2 forced case and realised that the west and east unprocessed boundary files west_unprocessed.nc and east_unprocessed.nc are missing dimensions xu_ocean and xt_ocean (north_unprocessed.nc and south_unprocessed.nc have these two dimensions and can be processed by expt.rectangular_boundary )
xu_ocean = UNLIMITED ; // (0 currently)

@luweiyang
Copy link
Collaborator

luweiyang commented May 10, 2024

@ashjbarnes @navidcy

I've fixed the west and east boundary condition issues. I changed cell 6 in step 2 (changes in bold):
rmom6.longitude_slicer(om2_input, [longitude_extent[1]-0.2, longitude_extent[1]+0.2], ["xu_ocean", "xt_ocean"]).to_netcdf(tmp_dir + "/east_unprocessed.nc")

The notebook is running fine now and I've pushed the working version to the ncc/regional-mom6-v2 branch.

However, I couldn't run the experiment as I keep getting warning & error messages about initialization and extreme surface state. So somehow the date of ocean forcing and atmospheric forcing do not match? Since we use repeated year forcing, I didn't think that would be a problem... Any thoughts? And it turns out that the temperature unit is K (instead of C) in init_tracers.nc.

 Entering coupler_init at 20240510 154506.457
 Starting initializing ensemble_manager at 20240510 154506.458
 Finished initializing ensemble_manager at 20240510 154506.458
 Starting to initialize diag_manager at 20240510 154506.489

WARNING from PE     0: diag_manager_mod::diag_manager_init: DIAG_MANAGER_NML not found in input nml file.  Using defaults.

 Finished initializing diag_manager at 20240510 154506.493
 Starting to initialize tracer_manager at 20240510 154506.502
 Finished initializing tracer_manager at 20240510 154506.512
 Beginning to initialize component models at 20240510 154506.512
 Starting to initialize atmospheric model at 20240510 154506.534
 Finished initializing atmospheric model at 20240510 154506.581
 Starting to initialize land model at 20240510 154506.595
 Finished initializing land model at 20240510 154506.607
 Starting to initialize ice model at 20240510 154506.636
 Finished initializing ice model at 20240510 154506.731
 Starting to initialize ocean model at 20240510 154506.739
WARNING from PE    54: Extreme surface sfc_state detected: i=  57 j= 128 lon= 145.825 lat= -43.366 x= 145.825 y= -43.366 D= 8.3479E+01 SSH= 5.6639E-02 SST= 2.8748E+02 SSS= 3.4988E+01 U-= 2.0583E-02 U+= 2.0757E-02 V-=-6.1662E-03 V+= 4.9824E-03

WARNING from PE    54: There were more unreported extreme events!

FATAL from PE    39: There were a total of     29869 locations detected with extreme surface values!

@navidcy
Copy link
Collaborator Author

navidcy commented May 10, 2024

Great!

something else we’d want to do is to use the cookbook to load the ACCESS-OM2 output

@@ -0,0 +1,699 @@
{
Copy link
Collaborator Author

@navidcy navidcy May 10, 2024

Choose a reason for hiding this comment

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

why there are 3 commented out lines of code @luweiyang?

if we need them, then uncomment,

if we don't need them then delete?


Reply via ReviewNB

Copy link
Collaborator

Choose a reason for hiding this comment

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

I don't think we need them. L20 and L21 are the same as L6 and L7 for selecting data along the dims yu_ocean and yt_ocean. I will delete them.

Copy link
Collaborator Author

@navidcy navidcy May 10, 2024

Choose a reason for hiding this comment

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

Generally, commented code is bad practice...
Either we want the code or we delete it.
It's good to let go. :)

@@ -0,0 +1,699 @@
{
Copy link
Collaborator Author

@navidcy navidcy May 10, 2024

Choose a reason for hiding this comment

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

why was this cell commented out @luweiyang? is the dask client create issues?


Reply via ReviewNB

Copy link
Collaborator

Choose a reason for hiding this comment

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

No, I commented them out when I was debugging because at some stage I did get some client-related errors. But it's no longer an issue. I will uncomment them.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

OK, do we need to add this back or not?

@navidcy
Copy link
Collaborator Author

navidcy commented May 10, 2024

something else we’d want to do is to use the cookbook to load the ACCESS-OM2 output

this is related to loading from cookbook: COSIMA/regional-mom6#131

@luweiyang
Copy link
Collaborator

Yes, I checked the date as well and was confused.

As @angus-g pointed out, the data in this folder
/g/data/ik11/outputs/access-om2-01/01deg_jra55v13_ryf9091/output1077/
is for year 2170, it goes from 2170/ 4/ 1 0: 0: 0 to 2170/ 7/ 1 0: 0: 0 (you can see this info in access-om2.out). Looks like we're actually getting boundary conditions from April in year 2170 as our ocean forcing?

The info for the atmospheric forcing is provided in data_table:
"ATM", "u_bot", "uas_10m", "./INPUT/RYF.u_10.1990_1991.nc", "bicubic", 1.0
Maybe we should update the year in filename to 2003_2004 here? Or it wouldn't matter since it's repeated year forcing?

Does this mismatch in year explain the initialisation error?

@ashjbarnes
Copy link
Collaborator

ashjbarnes commented May 10, 2024

I've fixed the west and east boundary condition issues. I changed cell 6 in step 2 (changes in bold): rmom6.longitude_slicer(om2_input, [longitude_extent[1]-0.2, longitude_extent[1]+0.2], ["xu_ocean", "xt_ocean"]).to_netcdf(tmp_dir + "/east_unprocessed.nc")

The notebook is running fine now and I've pushed the working version to the ncc/regional-mom6-v2 branch.

However, I couldn't run the experiment as I keep getting warning & error messages about initialization and extreme surface state. So somehow the date of ocean forcing and atmospheric forcing do not match? Since we use repeated year forcing, I didn't think that would be a problem... Any thoughts? And it turns out that the temperature unit is K (instead of C) in init_tracers.nc.

Yes good find! It'll be crashing because of the temperatures being in Kelvin. This is strange - these lines of code check the minimum surface temperature value of the dataset and then assume that it's in Kelvin if this minimum is above 100. This used to work but apparently not any more! The code here hasn't been changed for a long time

This notebook worked fine a while ago but it was neglected whilst we made a lot of changes. Somewhere we've introduced a change that breaks it but I don't know where..

@navidcy
Copy link
Collaborator Author

navidcy commented May 10, 2024

This is good. It will force us to find whe K->C doesn’t work and implement a test to capture future similar failures.

Let’s not ruminate on the past. Sure it worked at some point. It doesn’t matter that much, unless it helps us find why it’s not working now. Let’s just make an example that works now. It doesn’t even have to resemble how it was in the past if we want to showcase something else.

@navidcy
Copy link
Collaborator Author

navidcy commented May 10, 2024

Yes, I checked the date as well and was confused.

As @angus-g pointed out, the data in this folder /g/data/ik11/outputs/access-om2-01/01deg_jra55v13_ryf9091/output1077/ is for year 2170, it goes from 2170/ 4/ 1 0: 0: 0 to 2170/ 7/ 1 0: 0: 0 (you can see this info in access-om2.out). Looks like we're actually getting boundary conditions from April in year 2170 as our ocean forcing?

The info for the atmospheric forcing is provided in data_table: "ATM", "u_bot", "uas_10m", "./INPUT/RYF.u_10.1990_1991.nc", "bicubic", 1.0 Maybe we should update the year in filename to 2003_2004 here? Or it wouldn't matter since it's repeated year forcing?

Does this mismatch in year explain the initialisation error?

If we are gonna load the RYF then there is no point referring to “year 2003” or anything. All years are the same! They are looped.

I do think it more intuitive if we force with IAF, some particular year that we load via the cookbook

@ashjbarnes
Copy link
Collaborator

This is good. It will force us to find whe K->C doesn’t work and implement a test to capture future similar failures.

Let’s not ruminate on the past. Sure it worked at some point. It doesn’t matter that much, unless it helps us find why it’s not working now. Let’s just make an example that works now. It doesn’t even have to resemble how it was in the past if we want to showcase something else.

Knowing that it worked in the past should help us troubleshoot though right? Nothing about the forcing dataset has changed, only our code. If our code used to correctly convert K to C and no longer does, that means something in the package (specifically something in the initial condition method) has changed that prevents it from doing this.

I've had a quick look in case there was a premature "fill nans with 0s" or something but that doesn't seem to be the case

@luweiyang
Copy link
Collaborator

Yes, same ocean forcing year and same atmos forcing did work in the past. I remember a while ago @ashjbarnes updated and recompiled MOM6, and with new exe file everything worked fine.

@navidcy
Copy link
Collaborator Author

navidcy commented May 11, 2024

I never doubted that it worked in the past.. But now it's not, so let's fix it.

There are seem to be at least 1 or more issues. One is that temperature is not converted to K->C despite there is a method that supposedly does that. Possibly:

  • the method is never called, or
  • the method does not do what we think it does

How does ERA5-forced example work? Is ERA5 in degrees C to begin with?

@navidcy
Copy link
Collaborator Author

navidcy commented May 11, 2024

I've had a quick look in case there was a premature "fill nans with 0s" or something but that doesn't seem to be the case

This statement is unclear to me. I don't get what you want to say here @ashjbarnes.
"Premature" implying that it shouldn't be there... and now it "doesn't seem to be the case" so that mean "now it's not there anyway" or does it mean that "the fill nans with 0s is not premature" anymore?

@luweiyang
Copy link
Collaborator

Well done adding the buffer to the initial conditions, @navidcy ! Where did you make changes (which lines)?

It was running! After 50 mins, it crashed because of the following error:

FATAL from PE 97: time_interp_external 2: time 731234 (20030105.000010 is after range of list 731230-731234(20030101.000000 - 20030105.000000),file=INPUT/forcing/forcing_obc_segment_002.nc,field=u_segment_002

In the notebook, I've include comments about the version of regional-mom6 (v0.6.0) that users should load and added cells to plot initial conditions. Once we fix the bug about the time, it will be ready to be included in the demos.

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 20, 2024

Well done adding the buffer to the initial conditions, @navidcy ! Where did you make changes (which lines)?

See 12b4f7c

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 20, 2024

We'd also like to replace

om2_input = xr.open_mfdataset(
    f"/g/data/ik11/outputs/access-om2-01/01deg_jra55v13_ryf9091/output1077/ocean/ocean_daily*.nc",
    parallel=True, compat='override', chunks='auto')[["u", "v", "salt", "temp", "eta_t"]].sel(
        yu_ocean = slice(latitude_extent[0] - buffer, latitude_extent[1] + buffer),
        yt_ocean = slice(latitude_extent[0] - buffer, latitude_extent[1] + buffer)).isel(time = slice(0, 5))

with using the cookbook to load this data. Could we do that?

@edoddridge
Copy link
Collaborator

We'd also like to replace

om2_input = xr.open_mfdataset(
    f"/g/data/ik11/outputs/access-om2-01/01deg_jra55v13_ryf9091/output1077/ocean/ocean_daily*.nc",
    parallel=True, compat='override', chunks='auto')[["u", "v", "salt", "temp", "eta_t"]].sel(
        yu_ocean = slice(latitude_extent[0] - buffer, latitude_extent[1] + buffer),
        yt_ocean = slice(latitude_extent[0] - buffer, latitude_extent[1] + buffer)).isel(time = slice(0, 5))

with using the cookbook to load this data. Could we do that?

You mean intake, right?

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 20, 2024

We'd also like to replace

om2_input = xr.open_mfdataset(
    f"/g/data/ik11/outputs/access-om2-01/01deg_jra55v13_ryf9091/output1077/ocean/ocean_daily*.nc",
    parallel=True, compat='override', chunks='auto')[["u", "v", "salt", "temp", "eta_t"]].sel(
        yu_ocean = slice(latitude_extent[0] - buffer, latitude_extent[1] + buffer),
        yt_ocean = slice(latitude_extent[0] - buffer, latitude_extent[1] + buffer)).isel(time = slice(0, 5))

with using the cookbook to load this data. Could we do that?

You mean intake, right?

Yes, let's do intake!

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 20, 2024

@anton-seaice or @dougiesquire, could you help us here? we'd like to replace:

om2_input = xr.open_mfdataset(
    f"/g/data/ik11/outputs/access-om2-01/01deg_jra55v13_ryf9091/output1077/ocean/ocean_daily*.nc",
    parallel=True, compat='override', chunks='auto')[["u", "v", "salt", "temp", "eta_t"]].sel(
        yu_ocean = slice(latitude_extent[0] - buffer, latitude_extent[1] + buffer),
        yt_ocean = slice(latitude_extent[0] - buffer, latitude_extent[1] + buffer)).isel(time = slice(0, 5))

with intake.

I tried

import intake
catalog = intake.cat.access_nri

ds = catalog["01deg_jra55v13_ryf9091"].search(
    variable=["u", "v", "salt", "temp", "eta_t"],
    frequency='1day').to_dask()

but got

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[13], line 6
      1 import intake
      2 catalog = intake.cat.access_nri
      4 ds = catalog["01deg_jra55v13_ryf9091"].search(
      5     variable=["u", "v", "salt", "temp", "eta_t"],
----> 6     frequency='1day').to_dask()

File [/g/data/hh5/public/apps/miniconda3/envs/analysis3-24.07/lib/python3.10/site-packages/intake_esm/core.py:811](http://localhost:8890/g/data/hh5/public/apps/miniconda3/envs/analysis3-24.07/lib/python3.10/site-packages/intake_esm/core.py#line=810), in esm_datastore.to_dask(self, **kwargs)
    796 """
    797 Convert result to an xarray dataset.
    798 
   (...)
    808 :py:class:`~xarray.Dataset`
    809 """
    810 if len(self) != 1:  # quick check to fail more quickly if there are many results
--> 811     raise ValueError(
    812         f'Expected exactly one dataset. Received {len(self)} datasets. Please refine your search or use `.to_dataset_dict()`.'
    813     )
    814 res = self.to_dataset_dict(**{**kwargs, 'progressbar': False})
    815 if len(res) != 1:  # extra check in case kwargs did modify something

ValueError: Expected exactly one dataset. Received 9 datasets. Please refine your search or use `.to_dataset_dict()`.

@dougiesquire
Copy link
Collaborator

dougiesquire commented Aug 22, 2024

Because each variable is stored in a different file, intake will return multiple datasets. So you want:

import intake
catalog = intake.cat.access_nri

esmds = catalog["01deg_jra55v13_ryf9091"]

ds_dict = esmds.search(
    path=".*output1077.*",
    variable=["u", "v", "salt", "temp", "eta_t"],
    frequency='1day',
).to_dataset_dict()

ds = xr.merge(ds_dict.values(), join="inner")

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 22, 2024

Because each variable is stored in a different file, intake will return multiple datasets. So you want:

import intake
catalog = intake.cat.access_nri

esmds = catalog["01deg_jra55v13_ryf9091"]

ds_dict = esmds.search(
    path=".*output1077.*",
    variable=["u", "v", "salt", "temp", "eta_t"],
    frequency='1day',
).to_dataset_dict()

ds = xr.merge(ds_dict.values())

can I even drop the path=".*output1077.*"?
that's actually not important... it also makes the example cumbersome, that is others might seem that output1077 is some sort of special output.

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 22, 2024

I tried the above and after 10min nothing was happening... I think the dask cluster dropped. :(

@dougiesquire
Copy link
Collaborator

These particular data are actually quite problematic.

  • For one, the domain actually changes depending on which output??? directory you are looking at. Some data are on a global domain, some are on a southern ocean domain. There is nothing in the filename or metadata that indicates the domain is different. I actually opened an issue in the cosima-cookbook about this issue last year - see here.
  • Also, there are multiple filename conventions used within the one dataset. For example, there are both ocean_daily_3d_u.nc and ocean_daily_3d_u_MM.nc files (where MM is the month).

If you're wanting a command that opens a block of consistent, global data, you could filter out the *_MM.nc filenames:

ds_dict = esmds.search(
    variable=["u", "v", "salt", "temp", "eta_t"],
    frequency='1day',
    filename=".*[a-z].nc"
).to_dataset_dict()

ds = xr.merge(ds_dict.values(), join="inner")

Note, ds is 40TB so you will probably have to think about chunking if you actually want to do any calculation with it. It's hard to say more without knowing the details of what you're doing.

@dougiesquire
Copy link
Collaborator

If you're only selecting 5 time points then it's definitely better not to open >10 years of daily data, which is what my previous suggested command does.

If you don't like the path=".*output1077.*", you could select the start_date you want (that might be less confusing to users). E.g.

ds_dict = esmds.search(
    variable=["u", "v", "salt", "temp", "eta_t"],
    frequency='1day',
    start_date="2170-04-01, 00:00:00",
).to_dataset_dict()

ds = xr.merge(ds_dict.values(), join="inner")

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 22, 2024

OK, I see your point. We just want few days of daily data for demo purposes. Nothing beyond that.

@dougiesquire thanks for elaborating! Do you have, perhaps, a different dataset in mind that does not have these issues? RYF 0.1 degree?

@dougiesquire
Copy link
Collaborator

What's wrong with specifying the start_date as above? Then you only open a month of data to select a few days from. You could add a comment that indicates the exact value of start_date doesn't matter - just stick to the period "2170-01-01, 00:00:00" >= start_date <= "2179-10-01, 00:00:00" to ensure that all variables are available on the global domain

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 22, 2024

What's wrong with specifying the start_date as above? Then you only open a month of data to select a few days from. You could add a comment that indicates the exact value of start_date doesn't matter - just stick to the period "2170-01-01, 00:00:00" >= start_date <= "2179-10-01, 00:00:00" to ensure that all variables are available on the global domain

Nothing. That sounds great! I'll do that ;)

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 23, 2024

Can I also use end_date?

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 23, 2024

Done!!!!

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 23, 2024

This is ready. Could somebody (have a look) and approve?

@anton-seaice
Copy link
Collaborator

Can I also use end_date?

Not in a useful way. These start_date and end_date are not ranges which are searched between, they are just an attribute of the netcdf file being referred to.

e.g.

esmds.search(
    variable=["temp"],
    frequency='1day',
    start_date="2170-04-01, 00:00:00",
)

returns the data assets from one file which starts on April 1

esmds.search(
    variable=["temp"],
    frequency='1day',
    start_date="2170-04-01, 00:00:00",
    end_date="2170-04-10, 00:00:00",
)

returns nothing because there is no file containing daily data which starts on the first of april AND finishes on the 9th of april. It doesn't search for all data in that range.

If you understand regex you select a range by using a regex of start dates (something like this, untested:)

esmds.search(
    variable=["temp"],
    frequency='1day',
    start_date="2170-04-0.*, 00:00:00",
)

@navidcy
Copy link
Collaborator Author

navidcy commented Aug 23, 2024

thanks @anton-seaice !

Copy link
Collaborator

@dougiesquire dougiesquire left a comment

Choose a reason for hiding this comment

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

I haven't looked in detail at the whole notebook, but the intake stuff looks good.

@navidcy navidcy merged commit 87911ff into main Aug 23, 2024
3 checks passed
@navidcy navidcy deleted the ncc/regional-mom6-v2 branch August 23, 2024 06:16
@navidcy
Copy link
Collaborator Author

navidcy commented Aug 23, 2024

I haven't looked in detail at the whole notebook, but the intake stuff looks good.

I’ve scrutinised it many times and it runs. I’m merging and we can open a new PR to tweak if needed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Development

Successfully merging this pull request may close these issues.

Convert demo that uses ACCESS-OM2 output to cosima-recipes example
8 participants