Earth Radius in COMPASS and MPAS-Ocean init mode

[xylar]

A (sometimes a bit heated) email discussion has been going on about inconsistent Earth radii between:

• COMPASS when calling Jigsaw: 6371000
• MPAS framework constant a: 6371229.0
• CIME, EAM and ELM: 6.37122e6

[xylar]

@mark-petersen summarized as follows (with some formatting changes from me):

Here are the steps in the mesh generation

1. Base mesh (jigsaw): x,y,z coordinates have r=6371000 from triangle_jigsaw_to_netcdf. This value does not matter in the end.
2. MpasMeshConverter.x computes all mesh variables, including dcEdge and dvEdge
3. Cell culler – x,y,z remain the same
4. MPAS-Ocean init mode: Expands x,y,z to have r=6371229, expands dcEdge and dvEdge proportionally.

More details below. The gist is that step 1 could use a unit sphere. It doesn’t matter, because step 4 expands everything to the MPAS radius anyway.

1. triangle_jigsaw_to_netcdf.py

triangle_jigsaw_to_netcdf.py has a hardwired value of 6371000m that determines the (x,y,z) coordinates of the cell centers. Those coordinates are written out by jigsaw, along with the cellsOnCells array, i.e. neighbor relations.

2.MpasMeshConverter.x

Those coordinates are then passed to MpasMeshConverter.x which calculates dcEdge and dvEdge, as well as edge and vertex coordinates and all other index relations in the MPAS mesh. The distance calculation is here:
https://github.com/MPAS-Dev/MPAS-Tools/blob/master/mesh_tools/mesh_conversion_tools/mpas_mesh_converter.cpp#L1212

dvEdge.at(iEdge) = vert_loc2.sphereDistance(vert_loc1);
dcEdge.at(iEdge) = cell_loc2.sphereDistance(cell_loc1);

which is here:
https://github.com/MPAS-Dev/MPAS-Tools/blob/master/mesh_tools/mesh_conversion_tools/pnt.h#L253

arg = sqrt( powf(sin(0.5*((*this).getLat()-p.getLat())),2) +
cos(p.getLat())*cos((*this).getLat())*powf(sin(0.5*((*this).getLon()-p.getLon())),2));
return 2.0*asin(arg);

3. cell culler

The cell culler does not change radius

4. MPAS-Ocean init mode

mpas init mode does it here:
If we set config_expand_sphere = .true.
https://github.com/MPAS-Dev/MPAS-Model/blob/ocean/develop/src/core_ocean/mode_init/mpas_ocn_init_spherical_utils.F#L71

71    subroutine ocn_init_expand_sphere(domain, stream_manager, newRadius, err)!{{{
...
192            xCell(iCell) = (oldX / norm) * newRadius
193            yCell(iCell) = (oldY / norm) * newRadius
194            zCell(iCell) = (oldZ / norm) * newRadius
...
240            dvEdge(iEdge) = (dvEdge(iEdge) / oldRadius) * newRadius
241            dcEdge(iEdge) = (dcEdge(iEdge) / oldRadius) * newRadius

[xylar]

@proteanplanet suggests that the sphere radius for both (all) cases should be taken from CIME to be consistent with EAM and ELM:
https://github.com/E3SM-Project/E3SM/blob/master/cime/src/share/util/shr_const_mod.F90#L20
This value is: SHR_CONST_REARTH = 6.37122e6

[xylar]

@mark-petersen wrote:

More details below. The gist is that step 1 could use a unit sphere. It doesn’t matter, because step 4 expands everything to the MPAS radius anyway.

This was discussed further and it does not seem to be true for at least 2 reasons:

1. As both @dengwirda and @xylar pointed out, JIGSAW uses real lengths to determine the target resolution of the mesh it produces (so resolution would have to be normalized if we wanted to use a unit sphere)
2. There is an interest in using these meshes in situations where MPAS-Ocean init mode is not used (e.g. future land and river components).

[xylar]

cc @dengwirda, @proteanplanet, @sbrus89, @pwolfram, @mark-petersen, @jonbob

[xylar]

I like @proteanplanet idea of using CIME constants for our Earth radius.  We will just need to figure out the logistics of doing that:

• Which branch or commit of E3SM is used?
• How is the appropriate branch/commit downloaded by the tool (e.g. compute nodes may not be able to access the internet)?
• How to robustly parse out the constant we need from the Fortran file in python?
• How are the CIME constants made available to MPAS-Ocean during init mode?
• How to change MPAS-Ocean init mode to use the CIME constant instead of a from the MPAS framework?
• (Who will do all this work?!)

[dengwirda]

Quick additional background ---

This is important for mixed model workflows (i.e. ICoM) in which the MPAS mesh files created by MPAS_mesh_converter.cpp are handed off to various other components (ELM, MOSART, HexWatershed, etc) as well as to MPAS-O.

If MPAS-O re-normalises its mesh data using its own radius value there could be a mismatch in the key mesh variables (dvEdge, dcEdge, etc) used within each dycore.

Rather than requiring mesh re-normalisation on a per model basis, it'd be nice to have consistency throughout.

[xylar]

Makes sense. Eliminating, rather than replacing, @mark-petersen's step 4.

[proteanplanet]

Thanks @xylar for opening up this issue. Just to add to the comment on using CIME constants, it would be ideal if we were grabbing all of our physical constants, wherever possible, from this location. Like you, I'm not sure of the logistics. I'll give it some thought.

[klindsay28]

FYI, we're dealing with the basically the same problem in CESM/CIME/MOM.
In MOM, the radius of the Earth is read from a text-based parameter file, while in CIME it is in Fortran source code. We would like to use the CIME value in MOM when running with CIME.

To achieve this, we're leaning towards adding a set parameter method to MOM's parameter infrastructure, and the CIME cap for MOM would call that, passing in values from CIME's shr_const_mod.F90. This applies to other constants besides Earth's radius. A feature of this approach that we like is not introducing multiple copies of the same values that have the potential to get out of sync.

This approach may or may not work for CIME/MPAS-O.

[xylar]

@klindsay28, thanks, that's some really helpful insight! I was thinking along the same lines that it would be ideal if CIME could inform MPAS-Ocean of parameters. (This would not apply to the Earth radius, which is needed at initialization but to a lot of other constants.)