MOM Supergrid 100km standalone

[anton-seaice]

PR checklist

• Short (1 sentence) summary of your PR:
  This adds a standalone cice test case to cover MOM supergrid, and tripolar grids

Closes #977

• Developer(s):
  @anton-seaice @claude @apcraig
• Suggest PR reviewers from list in the column to the right.
  @apcraig
• Please copy the PR test results link or provide a summary of testing completed below.
  ENTER INFORMATION HERE
• How much do the PR code changes differ from the unmodified code?
  • bit for bit
  • different at roundoff level
  • more substantial
• Does this PR create or have dependencies on Icepack or any other models?
  • Yes
  • No
• Does this PR update the Icepack submodule? If so, the Icepack submodule must point to a hash on Icepack's main branch.
  • Yes
  • No
• Does this PR add any new test cases?
  • Yes
  • No
• Is the documentation being updated? ("Documentation" includes information on the wiki or in the .rst files from doc/source/, which are used to create the online technical docs at https://readthedocs.org/projects/cice-consortium-cice/. A test build of the technical docs will be performed as part of the PR testing.)
  • Yes
  • No, does the documentation need to be updated at a later time?
    • Yes
    • No
• Please document the changes in detail, including why the changes are made. This will become part of the PR commit log.

This change improves test coverage for repro testing, by including the MOM supergrid and tripolar grid from a netcdf file. New grid is named tx1m, generated through FRENCTools. See the metadata in the grid and forcing files for how to regenerate files. A python script is added to generate CICE standalone forcing files from JRA55do data. The resultant forcing data for 2005 will also be added to the CICE consortium zenodo archive.

[anton-seaice]

@apcraig - if the gist of this looks ok, ill get the test data to you.

The script to make the test data is here. It might standalone ok, so could be added to this repo as well ?

I added a new test suite to cover smoke and restart cases, and both b and c grid. If that looks ok, I guess it should be incorporated into base_suite ?

[apcraig]

Thanks @anton-seaice, it's great to get mom grid testing in the standalone model! I ran some tests with the new grid on derecho. A few comments about the current implementation. If you want me to work on some of these issues and PR changes to your branch, I'm happy to do that. But also happy for you to take the lead.

• With debug on, tested with a few compilers, I had to make the following changes to get the model to run. The issue is that on non-master tasks, G_T is allocated (1,1) but never initialized. With debug, these are initialized to nan and are trapped when multiplied by deg_to_rad. Also, there are specific array references to G_T like (1:nx_global,1:ny_global) in the scatter call in the mom grid implementation, but when the arrays are allocated (1,1), an error is produced on the gnu compiler on non-master tasks even though it doesn't impact the results. I used the following mods, but other fixes could be implemented.

diff --git a/cicecore/cicedyn/infrastructure/ice_grid.F90 b/cicecore/cicedyn/infrastructure/ice_grid.F90
index 6ae922be..18f68dc4 100644
--- a/cicecore/cicedyn/infrastructure/ice_grid.F90
+++ b/cicecore/cicedyn/infrastructure/ice_grid.F90
@@ -1721,7 +1721,7 @@ subroutine mom_grid
       ! lat, lon, angle
       !-----------------------------------------------------------------
 
-      if (my_task == master_task) then
-         allocate( &
-            work_mom(nx_global*2+1, ny_global*2+1), &
-            work_gE(nx_global+1,ny_global+1)      , & !include left and top
-            work_gN(nx_global+1,ny_global+1)      , & !include right and top
-            G_ULAT(nx_global+1,ny_global+1)       , & !include left and bottom
-            G_TLAT(nx_global+1,ny_global+1)       , & !include top and right
-            G_TLON(nx_global+1,ny_global+1)       , & !include top and right
-            G_ULON(nx_global+1,ny_global+1)       , & !include left and bottom
-            stat = ierr &
-         )
-      else
-         allocate(work_mom(1,1), work_gE(1,1), work_gN(1,1), &
-            G_ULAT(1,1), G_TLAT(1,1), G_TLON(1,1), G_ULON(1,1), &
-            stat=ierr)
-      endif
+      allocate( &
+         work_mom(nx_global*2+1, ny_global*2+1), &
+         work_gE(nx_global+1,ny_global+1)      , & !include left and top
+         work_gN(nx_global+1,ny_global+1)      , & !include right and top
+         G_ULAT(nx_global+1,ny_global+1)       , & !include left and bottom
+         G_TLAT(nx_global+1,ny_global+1)       , & !include top and right
+         G_TLON(nx_global+1,ny_global+1)       , & !include top and right
+         G_ULON(nx_global+1,ny_global+1)       , & !include left and bottom
+         stat = ierr &
+      )
       if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory', file=__FILE__, line=__LINE__)
 
       ! populate all LAT fields
@@ -2012,10 +2012,12 @@ subroutine mom_corners_scatter(G_U, G_T, G_E, G_N, U, T, E, N )
       deg_to_rad = pi/c180
 
       ! convert to rad
-      G_T = G_T * deg_to_rad
-      G_U = G_U * deg_to_rad
-      G_N = G_N * deg_to_rad
-      G_E = G_E * deg_to_rad
+      if (my_task == master_task) then
+         G_T = G_T * deg_to_rad
+         G_U = G_U * deg_to_rad
+         G_N = G_N * deg_to_rad
+         G_E = G_E * deg_to_rad
+      endif
 
       ! distribute to processors
       ! subset G_T to active cells by dropping right/top halo

• The name of the grid in this PR right now is access_100km. I think we should change the name of the grid. In particular, we should have a short but somewhat identifiable and extendable name, we probably want to avoid associating the grid with any particular model or group, we prefer lowercase, and we want to avoid grid names (and run options) with underscores or other special characters in them. Is the gridname in ACCESS actually access_100km? Maybe the name should be associated with mom, but that's the grid format, and this grid could also exist as a pop grid. Maybe this should be called tx1a, mx1, or something like it?

• I think the set_nml."$grid" file should have some of the settings that are required to run. For instance, in the test suite, this grid needs to have "jra55do,icnone,bctripole" all set. I think these should be set by default for the grid in the set_nml."$grid" file. See set_nml.[gx3,gx1,tx1] to see how the other grid settings are doing this. I would even explicitly set the timestep and stuff like year_init because those are kind of needed for the forcing data we have (even if the defaults work). I know, many of the defaults are NOT set in these files, but I think it's a good idea to more-or-less follow the established pattern.

• Should we have more decomp options for the grid in cice_decomp.csh?

• We will want to add some tests to the test suites. We should have some basic smoke and restart tests in the base suite, without and without debug. We can add some gridc/cd tests in the gridsys_suite and maybe some decomp tests in the decomp_suite. We might even want one or two unit tests to use this grid, not sure it's necessary, but might be good to do. For instance, grep for "tx1" in configuration/scripts/tests/*.ts to see what we're doing for the tx1 grid. This grid might have similar tests, but lets worry about the tests after we clean up the current implementation and decide on a gridname. I think we'll want to get rid of the access_100km_cgrid test suite unless it's useful. In that case, happy to keep it around.

• I love the idea of including the python forcing data generation script if it would be helpful. There are some scripts (configuration/tools/jra55_datasets) and documentation, https://cice-consortium-cice.readthedocs.io/en/main/developer_guide/dg_forcing.html already. We would probably add your script there and would want to add some info in the user guide if that makes sense.

• I dropped the new dataset into an area I can point the model to. We need to formalize the grid names, setup the input files as needed, probably do a 1 year test run (happy to do that) to make sure the full year of forcing data is stable, put the inputdata on zenodo, and migrate the new dataset onto various test machines.

[anton-seaice]

• With debug on, tested with a few compilers, I had to make the following changes to get the model to run.

Oh thanks - ok replicated and applied these.

• The name of the grid in this PR right now is access_100km. I think we should change the name of the grid. In particular, we should have a short but somewhat identifiable and extendable name, we probably want to avoid associating the grid with any particular model or group, we prefer lowercase, and we want to avoid grid names (and run options) with underscores or other special characters in them. Is the gridname in ACCESS actually access_100km? Maybe the name should be associated with mom, but that's the grid format, and this grid could also exist as a pop grid. Maybe this should be called tx1a, mx1, or something like it?

ACCESS don't really have names for grids as such - just names for configurations. e.g. this is the grid from access-om3 100km configurations.

The current cice tests grids are from POP. I think the names have POP specific meanings ? So the idea here was just to make the grid clearly different. I do see the issue with the underscore now though.

I guess there might be MOM tripole, MOM regional and MOM offset pole grids in the future, so some name to reflect those might make more sense. I think the name depends a bit on how likely we are to add configurations for each modelling centre (e.g. #971 )

• I think the set_nml."$grid" file should have some of the settings that are required to run. For instance, in the test suite, this grid needs to have "jra55do,icnone,bctripole" all set. I think these should be set by default for the grid in the set_nml."$grid" file. See set_nml.[gx3,gx1,tx1] to see how the other grid settings are doing this. I would even explicitly set the timestep and stuff like year_init because those are kind of needed for the forcing data we have (even if the defaults work). I know, many of the defaults are NOT set in these files, but I think it's a good idea to more-or-less follow the established pattern.

I'll add the bctripole settings. In theory, there's nothing specific about this grid to the initial conditions or forcing. e.g. an idealised forcing would work too, as would other initial conditions, so it make some sense for the other settings to be customisable ?

• Should we have more decomp options for the grid in cice_decomp.csh?

I don't see this as critical, the code paths for decomposition are fairly independent to the grid setup

• We will want to add some tests to the test suites. ... I think we'll want to get rid of the access_100km_cgrid test suite unless it's useful. In that case, happy to keep it around.

Yes agreed - it was just an example suite to get it started.

• I love the idea of including the python forcing data generation script if it would be helpful. There are some scripts (configuration/tools/jra55_datasets) and documentation, https://cice-consortium-cice.readthedocs.io/en/main/developer_guide/dg_forcing.html already. We would probably add your script there and would want to add some info in the user guide if that makes sense.

Ok, ill work on this. I didn't look at the existing script too long because the variable names didn't make sense to me. I guess they are from jra55 instead of jra55do ?

• I dropped the new dataset into an area I can point the model to. We need to formalize the grid names, setup the input files as needed, probably do a 1 year test run (happy to do that) to make sure the full year of forcing data is stable, put the inputdata on zenodo, and migrate the new dataset onto various test machines.

I have done 1 year runs already, and they worked for me

[apcraig]

Thanks Anton. How about if we call the grid tx1a? @eclare108213 @dabail10 any thoughts? This is a tripole, 100km grid (about 1 degree), in mom netcdf format, similar to tx1 in some ways but not tx1. Our grid prefixes now are gx (dispaced pole), tx (tripole), and gbox (box grids) and these are pop formats so far, whether binary or netcdf, but we have not really indicated whether a grid is binary, netcdf, pop, etc in the grid name which I think is right. In that case, this grid should be "tx1" but given a unique name. tx1a would be a new tripole 1deg grid. It happens to be mom netcdf format, but I don't think that should be reflected in the grid name, although it will be defined in the set_nml.$grid file.

Also, I think we should set the default forcing, timestep, and some other things in the set.nml.$grid like the other grids. I understand that other forcing and settings can work, but the settings in the grid file establishes the baseline/preferred test configuration and the default settings for that grid. Those can be changed with other options or manual changes to the ice_in file. Fundamentally, I think we want to avoid having to know that other options SHOULD be set anytime a particular grid is chosen. That gets complicated and confusing.

[dabail10]

We still use the tx for our tripole grids. Our workhorse is now 2/3 degree and we call it tx2_3v3 with a version number attached.

[anton-seaice]

Thanks Anton. How about if we call the grid tx1a? @eclare108213 @dabail10 any thoughts? This is a tripole, 100km grid (about 1 degree), in mom netcdf format, similar to tx1 in some ways but not tx1. Our grid prefixes now are gx (dispaced pole), tx (tripole), and gbox (box grids) and these are pop formats so far, whether binary or netcdf, but we have not really indicated whether a grid is binary, netcdf, pop, etc in the grid name which I think is right. In that case, this grid should be "tx1" but given a unique name. tx1a would be a new tripole 1deg grid. It happens to be mom netcdf format, but I don't think that should be reflected in the grid name, although it will be defined in the set_nml.$grid file.

Could tx1a get mixed up as a variant of tx1 ? What about tm1 as a name ? Tripole + mom + 1 degree.

Also, I think we should set the default forcing, timestep, and some other things in the set.nml.$grid like the other grids. I understand that other forcing and settings can work, but the settings in the grid file establishes the baseline/preferred test configuration and the default settings for that grid. Those can be changed with other options or manual changes to the ice_in file. Fundamentally, I think we want to avoid having to know that other options SHOULD be set anytime a particular grid is chosen. That gets complicated and confusing.

Sorry, I think i get it now. I guess the default values in ice_in work in other grids but not this one.

[eclare108213]

Could tx1a get mixed up as a variant of tx1 ? What about tm1 as a name ? Tripole + mom + 1 degree.

In the pop grids, I think the 'x' means 'by' as in '2 by 2' but only referring to the second part (~longitude). Grid names have become extremely complicated (and long) with version numbers and all kinds of other information included in the name itself. I'd prefer to keep it simple for our purposes, and I also prefer not to name anything that might be used more generally after a specific institution. But perhaps that is appropriate in this case. Is this a general-purpose mom grid that access is using, or is it a mom grid that access developed for its own use? Is 100 km a global average? I think we need the t for tripole and either 100km or x1 for resolution, then maybe something to distinguish this grid from the pop grids. How about something like t100km? tx1m or tmx1 would be okay too. Is the 'a' is for 'access' in tx1a? Is this a "supergrid"? If so, maybe tx1super?

[anton-seaice]

Is this a general-purpose mom grid that access is using, or is it a mom grid that access developed for its own use?

"Developed" is a stretch, but yes we made the grid (using MOM tools) based on some goals relevant to us. (e.g. more Southern Ocean cells)

Is 100 km a global average?

it's a fairly approximate name. Its 1 degree in longitudes (ignoring the tripole), which is very roughly 100km at the equator.

I think we need the t for tripole and either 100km or x1 for resolution, then maybe something to distinguish this grid from the pop grids. How about something like t100km? tx1m or tmx1 would be okay too. Is the 'a' is for 'access' in tx1a? Is this a "supergrid"? If so, maybe tx1super?

Yeah, its a supergrid, but we've used the name mom in the code (rather that supergrid). tx1m or tmx1 are both good to me ?

[anton-seaice]

I think the 'x' means 'by' as in '2 by 2' but only referring to the second part (~latitude).

ah. The 1 degree in this grid is longitude, latitudes vary a fair bit :

The grid is defined using the conventional tripolar definition(Murray, 1996)1. The grid is Mercator (i.e. the meridional spacing scales as the cosine of latitude) between 65° N and 75° S; south of 75° S, the meridional grid spacing is held at the same value as at 75° S. The meridional spacing is refined between 10° S and 10° N to increase equatorial resolution.

[apcraig]

Could tx1a get mixed up as a variant of tx1 ? What about tm1 as a name ? Tripole + mom + 1 degree.

In the pop grids, I think the 'x' means 'by' as in '2 by 2' but only referring to the second part (~latitude). Grid names have become extremely complicated (and long) with version numbers and all kinds of other information included in the name itself. I'd prefer to keep it simple for our purposes, and I also prefer not to name anything that might be used more generally after a specific institution. But perhaps that is appropriate in this case. Is this a general-purpose mom grid that access is using, or is it a mom grid that access developed for its own use? Is 100 km a global average? I think we need the t for tripole and either 100km or x1 for resolution, then maybe something to distinguish this grid from the pop grids. How about something like t100km? tx1m or tmx1 would be okay too. Is the 'a' is for 'access' in tx1a? Is this a "supergrid"? If so, maybe tx1super?

My feeling is the grid name should not be associated with the grid format, whether that be pop or mom or supergrid, binary or netcdf. Nor should it necessarily be associated with the modeling center although it could be if it's unique enough. The grid specifically just defines some 2d lon/lat. The type of grid, displaced pole, tripole, rectangular, etc are part of our naming standards, so that's useful. The "a" in tx1a just means version "a" although it could also mean access (or anton) if anyone wanted to think of it that way. Dave's example of tx2_3v3 is going to be a little problematic for standalone CICE, I'm trying to avoid special characters. Maybe CICE would call that grid tx23v3 or txp7v3? We'll deal with that if it ever comes up. Like many things in software, conventions get established that are then sometimes difficult to extend later. I think we want unique gridnames that contain some info about the grid and are unique and short. This new grid is a "tx1" grid and we just need to uniquely identify it. I would suggest something like "tx1a", "tx1m", "tx1v1" "tx1v2", "tx1a1", "tx1mv1" or similar where the a, m, v1, v2, a1, mv1 suffix can be interpreted anyway you like :)

[anton-seaice]

Ok - ill go for tx1m

[anton-seaice]

I think I've addressed all the comments except for updating the test suites. Could you do that @apcraig ?

[apcraig]

I have started testing again and will propose some changes to the test suite shortly. Thanks @anton-seaice

[apcraig]

I just submitted a PR with updates for the access_100km_standalone branch. Let me know if there are any concerns. Will retest once we converge on the updates.

[anton-seaice]

Suggested change

	``gx3``, ``gx1``, ``tx1``, ``tx1m`` are associate with grid specific settings
	``gx3``, ``gx1``, ``tx1``, ``tx1m`` are associated with grid specific settings

[apcraig]

@anton-seaice, feel free to continue to make changes, especially if you find fixed in documentation. I will start a larger test suite, do you agree we're ready to start "final" testing?

[anton-seaice]

Ready to go @apcraig :

• should years 2006-2009 be added to the forcing dataset?
• someone will need to put the forcing data on zenodo

[eclare108213]

Looks like this PR is nearly ready. @anton-seaice you should be able to submit the forcing dataset to the CICE Consortium community on zenodo, and I can accept it. Is the data in our standard format? In addition to zenodo, links and other info about the forcing will need to be added to the CICE forcing wiki page. Are you able to edit that one?

[anton-seaice]

Looks like this PR is nearly ready. @anton-seaice you should be able to submit the forcing dataset to the CICE Consortium community on zenodo, and I can accept it. Is the data in our standard format?

Sounds good. The format follows the other forcing data as much as possible (e.g. variable names, folder structure etc).

I've only made data for 2005, should I include 2006-2009 as well?

In addition to zenodo, links and other info about the forcing will need to be added to the CICE forcing wiki page. Are you able to edit that one?

I don't have access but can write the update

[apcraig]

This is for testing the mom grid capability, I think one year of forcing data is more than adequate. I (or Elizabeth) can update the wiki, thanks Anton for putting some words together.

[apcraig]

I'm running a test suite, will approve and merge if that goes smoothly. Thanks!