Question regarding MITgcm code comment on "coefficient 4 changed"

[navidcy]

There is a note here:

http://wwwcvs.mitgcm.org/viewvc/MITgcm/MITgcm_contrib/high_res_cube/matlab-grid-generator/map_xy2xyz.m?view=markup#l73
regarding coefficient 4 changed -0.00895884801823 to -0.01895884801823.

Does anyone have any idea why that was done there?

@christophernhill @jm-c @ali-ramadhan

[ali-ramadhan]

I did notice that and might have asked about it but I don't remember why it was done.

With that change for coefficient $A_4$ I think the $B$ coefficients for $Z(w)$ in map_xy2xyz.m are actually wrong since the $B$ coefficients are obtained by inverting the Taylor series for $W(z)$ with the $A$ coefficients. The coefficients in this repo should be correct since we have a test that checks that they are correct and they match table B1 of Rančić et al. (1996).

[navidcy]

Sure. But if there is a typo in Rancic paper then the tests would still be passing. :)

[navidcy]

Seems like GFDL's FRE-NCtools also use -0.0189....

https://github.com/NOAA-GFDL/FRE-NCtools/blob/83acb799f47dfa27b433a131e9f7c1310767cc59/tools/make_hgrid/create_conformal_cubic_grid.c#L272

But possibly this is just copied from the MITgcm scripts (as the comment above mentions)

[glwagner]

Nice locating the FV3 code!

[navidcy]

I reached out to Alistair to see if he remembers why that extra digit was added in coefficient $A_4$. He came back with a plot of the two versions of the coefficients. If one had to guess from the plot you'd say that Alistair's version is the correct.

using CubedSphere, GLMakie

A_Rancic = deepcopy(CubedSphere.A_Rancic)

A_Alistair = deepcopy(CubedSphere.A_Rancic)
A_Alistair[5] = -0.01895883606818

k = 3:11

fig = Figure(fontsize=28)
ax = Axis(fig[1, 1], xlabel = "k", ylabel = "log10(- Aₖ)")
scatter!(ax, k .- 1, log10.(-A_Rancic[k]), label="Rancic", markersize=20, marker=:diamond)
scatter!(ax, k .- 1, log10.(-A_Alistair[k]), label="Alistair", markersize=16)
axislegend(ax)
display(fig)

However, I followed Rancic et al (1996) steps to compute the coefficients. I created a script and when I follow Rancic iterative method, after 30 iterations I converge to coefficients exactly as reported in the 1996 paper. (See table in #18).

I'm trying to reconcile the figure above with the fact that iterative method gives me exactly the coeffs that Rancic reported back in 1996. (I'm drafting some notes to explain the method by Rancic and will post here soon.)

[navidcy]

Here are some notes regarding how I computed the conformal map coefficients. Seems like my method agrees with the coefficients reported by Rancic. However, the plot above with the $A_4$ corrected by Alistair seems too pretty to be wrong...

CubedSpheranigans.pdf

[kburns]

I think it's pretty likely we can write an entirely independent and much more efficient remapping using rational functions. What are the absolute tolerances you want to actually compute the remapping to? And does the rest of the code just require the forward and backwards maps as functions, or also the local scaling and local orientation? Or are those computed upstream via finite differences?

[kburns]

Here's an independent algorithm using least squares fitting and rational functions ("lightning" rational approximations) to compute the forward and inverse maps: https://github.com/kburns/conformal_cubed_sphere.

The resulting maps match those from this repository using the original Rancic coefficients to full precision. The algorithm is much simpler than Rancic, but in the end it looks like you still need about as many coefficients to store the map to full precision (although we may be able to improve on this using log-lightning approximations). So it's not clear it would offer any performance benefits to switch implementations, but at least it confirms the original coefficients. Modifying $A_4$ to -0.0189... results in an O(1e-5) error in the maps, so this should probably be reverted where implemented.