I noticed this error in my testing of MPAS-Dev/MPAS-Model#668.

Further work toward addressing MPAS-Dev/MPAS-Model#549

@mark-petersen and @dengwirda, I'm just looking for a review by inspection here and wanted to make sure you're aware of this change just in case there is code outside of MPAS-Tools and MPAS-Model that will be affected.

This will be part of the next mpas_tools release (0.0.13). I'm also working on fixing temp directory issues in #346 that will be part of that release.

As long as you've tested with a global case and see the expected radius, then this is working.

Finding that out is a bit involved. I have to make a new test version of mpas_tools and and compass to run with. Having just done that twice recently, I'm going to hold off until we get some other things taken care of. But I won't merge this in until I've verified that it works.

@xylar based on a quick review it looks to me like the flow is right here: earth_radius comes from constants => build_spherical_mesh => jigsaw_driver and is no longer hard-coded in jigsaw_to_netcdf.

Is there possibly a m <=> km issue though? In jigsaw_driver it looks like earth_radius / 1000. converts to km, which (I think) means the x, y, z read in jigsaw_to_netcdf are still in km?

If this seems right to you, perhaps x, y, z * 1000. could be done around ln 57 in jigsaw_to_netcdf? Or maybe the mesh could be generated in m direct in jigsaw_driver?

I can also do a full build and double-check I have the right understanding here...

@dengwirda, I agree, it would be cleaner if everything were in meters, rather than sometimes in km but I think that's separate from this PR. Here, the only thing that is getting modified is an attribute in the output NetCDF file and I'm sure both the value going in and the original hard-coded value were in meters, so I don't think there are any new problems introduced here.

I think we should do a separate PR if we want to convert the units that JIGSAW gets to meters instead of km.

I will double check, though, to see if the xCell, yCell and zCell in the NetCDF file are in km. I agree, if that's the case, that's a problem we need to fix here.

Update: I don't have an easy test of this branch but before the change here, xCell, etc. are definitely in meters:

xCell = 5754.987991377, -11939.736813286, 5700836.48866784, 
    4605609.42203388...

I agree, though, that I can't find anywhere where x, y and z from JIGSAW, which surely must be in km, get converted back to meters. Mysterious! @sbrus89, any insight? Were you involved in writing this code or was in Phil Wolfram?

On a side note, I need to fix the command-line tool. It needs a new flag for the Earth radius.

Oh, dear! I was looking at base_mesh.nc but If I look at mesh_triangles.nc instead, I see:

data:

 xCell = 5.75518671957636, -11.9401491099583, 5701.03334693617, 
    4605.76846050082, 1759.38392552772, -1763.68148114646, -4607.79250490292, ...

So definitely in km. It's getting converted to m via renormalizaiton that happens in the MPAS mesh converter. Hopefully, that's using the supplied sphere_radius and not yet another hard-coded value. Let me try to find out.

So, @dengwirda, you're completely right. We should be multiplying x, y and z by 1000 here and we might as well fix it in this PR.

@dengwirda, thanks for your comments. They should now be addressed.

@proteanplanet, would you be willing to review this as well. I believe this is the last of these darn hard-coded Earth radii on the MPAS-Tools and COMPASS side of things. But you may have a keener eye for something I missed.

It's getting converted to m via renormalizaiton that happens in the MPAS mesh converter. Hopefully, that's using the supplied sphere_radius and not yet another hard-coded value. Let me try to find out.

I verified that this is the case -- there's no hard-coded Earth radius in the mesh converter.

@xylar thanks for following up on all of this. To multiply out the km => m though I think we need to be careful about the sphere vs planar workflows, so as not to break anything for @matthewhoffman et al in LI.

I think the x, y, z * 1000. could be done inside the if on_sphere: block (line 57) in jigsaw_to_netcdf, since it looks like it's only the on_sphere path in jigsaw_driver that sees earth_radius / 1000.. Does this seem right?

Oh, shoot, you're right. I'll fix that.

@matthewhoffman isn't yet using this workflow for any of his test case as far as I'm aware but could certainly do so in the future. There may be some ocean test cases using JIGSAW in a planar configuration, not sure about that. In any case, it should work, so thanks for catching that. Let me know if it looks right now.

Looks good I think @xylar!

I ran the QU240 and QU240wISC test cases from COMPASS ocean. I can verify that:

• the sphere_radius attribute in both the base mesh and the culled mesh are the CIME value: sphere_radius = 6371220.
• the sphere_radius still gets replaced by the MPAS-Ocean value during initial_state: sphere_radius = 6371229. This will need to be dealt with separately and at a later date (see Earth Radius in COMPASS and MPAS-Ocean init mode MPAS-Model#549)
• xCell (and presumably yCell and zCell) in mesh_triangles.nc are now in m instead of km.

I will merge this PR and it will make it into the next release, which will happen shortly.

@dengwirda , @xylar , landice is using this jigsaw interface for mesh generation in COMPASS (we refer to the ocean dir within COMPASS to access these scripts, which seemed easier than asking the ocean to generalize this and put it in a common folder), so thanks for keeping an eye out for the planar workflow.

Ah, @matthewhoffman, I didn’t realize this. I probably need to make a PR so you can use the latest compass instead of jigsaw to MPAS. I’ll do that tomorrow if possible. All the more reason that the move to this repo was a good call.