What is the official gx1v7 → r1x1 conservative regridding procedure for CESM2 OMIP2?

[Zhi-Jian]

I am trying to regrid CESM2 ocean output from the native gx1v7 (384×320) grid to the standard regular r1x1 (180×360) grid.
For comparison with the official NCAR products, I used the following two files as a test case:

• tos_Omon_CESM2_omip2_r1i1p1f1_gn_024501-030512.nc
• tos_Omon_CESM2_omip2_r1i1p1f1_gr_024501-030512.nc

In the gr file, it says :grid = "ocean data regridded from native gx1v7 displaced pole grid (384x320 latxlon) to 180x360 latxlon using conservative regridding", so I attempted to reproduce this regridding using NCO with conservative remapping (the ESMF Offline Regridding Weight Generator has been installed internally in the NCO):

ncremap -v tos --alg_typ=conserve \
-d tos_Omon_CESM2_omip2_r1i1p1f1_gr_024501-030512.nc \
-i tos_Omon_CESM2_omip2_r1i1p1f1_gn_024501-030512.nc \
-o tos_regridded_180x360.nc

I then compared my regridded result with the official NCAR “gr” file. While the large-scale patterns are very similar, I found noticeable differences near coastlines, where the official gr product appears to be much better behaved (see attached figure).

Therefore, I would like to ask:

1. What exact procedure and tools were used at NCAR to convert CESM2 ocean output from gx1v7 (gn) to r1x1 (gr) for OMIP2?
2. Were pre-generated mapping/weight files used for this conversion?
   (This question was raised in a previous thread but does not seem to have received a definitive answer; see posts #11–14 in:
   CESM raw output variable corresponding to cland)
3. If possible, would it be feasible to share the official gx1v7 → r1x1 conservative remapping weight file used in production?

I have also done some searching. For example, the repository
- Revision 70792: /trunk/ncl
contains several pre-generated regridding weight files, but I could not find a conservative mapping file for gx1v7 → r1x1 (it has gx1v6 → r1x1 using bilinear remapping).

Overall, I have not yet been able to reproduce the quality of the official CMIP6 gn → gr product, which is remarkably good, especially in marginal seas and coastal regions.
Any guidance or pointers would be greatly appreciated.

 

[sacks]

Thanks for all of the details in your question. I'll reach out to see if I can find someone who knows this answer definitively, but due to the recent departure of some staff, I may have a hard time finding someone who can answer this definitively.

In the meantime, I can suggest a few things based on what I know in general about how we do regridding of data in general for CESM.

An important piece in doing this remapping with many tools is providing a mask on the source grid (and sometimes the destination grid), and describing how to do renormalization. I don't have much experience with ncremap, but from a quick skim of the documentation, the following options seem like they could be important: --msk_src (and maybe --msk_dst), --preserve (or maybe -r) (note: I see, "NCO version 5.1.5, released in March, 2023, fixed a longstanding problem with the implmentation of this option, which had essentially been broken since its inception. The option finally works as documented.")

If you have a build of ESMF, you could also simplify the problem by directly using the command-line tools ESMF_RegridWeightGen and ESMF_Regrid (see Regridding). These have fewer options than ncremap but I would guess that these tools support any option that was used for the CMIP6 workflow.

 

[mlevy]

I then compared my regridded result with the official NCAR “gr” file. While the large-scale patterns are very similar, I found noticeable differences near coastlines, where the official gr product appears to be much better behaved

The values near the coastline look different because the CESM ocean community works from the assumption that any active cell on the ocean grid is entirely ocean, there is no concept of partial ocean cells. However, the ocean mask on the gx1 grid does not exactly match the ocean mask on the uniform 1° grid... if you map a field of all 1s conservatively, you will have values between 0 and 1 along the coast. These values are the fraction of the cell on the destination grid that receive values from ocean cells on the source grid (the rest of the cell is mapped from the masked out land values).

An alternate way to do the mapping is to divide the values on the destination grid by the ocean fraction, in which case a field of 1s on the source grid would also be a field of 1s on the destination grid... but when computing global integrals you then need to remember to compute sum(values*area*fraction) instead of just sum(values*area), and as mentioned that is not the way oceanographers typically compute integrals.

What exact procedure and tools were used at NCAR to convert CESM2 ocean output from gx1v7 (gn) to r1x1 (gr) for OMIP2?

I hesitate to call this the "exact" procedure, but my understanding is that we came up with a pseudo-conservative alternative that allows users to continue to compute integrals with sum(values*area) and also maintains a field of constant values on the destination grid. This is done by modifying the existing conservative grid in the following manner:

1. Compute the ocean area on the source grid (this is just sum(area*mask), where mask is 1 for ocean cells and 0 for land cells); call this area_src
2. Compute the ocean fraction on the destination grid by mapping a field of 1s (as mentioned, this will be 1 in the open ocean and have values between 0 and 1 along the coast
3. Define a function f(x) = sum(dst_area where ocn_frac >= x)
4. It should be clear that f(0) = sum(dst_area) >= area_src = sum(dst_area*ocn_frac) >= sum(dst_area in open ocean) = f(1)
5. If f(x) were a continuous function, we'd be guaranteed the existence of a value of x_0 such that f(x_0) = area_src... however, f(x) is not continuous. Stlll, we can find an x that minimizes abs(f(x) - area_src)

We then modify our conservative map to remove destination cells where ocn_frac < x and ensure the weights sum to 1 where ocn_frac >=x - the result is that mapping a field of 1s results in a field of 1s on the destination grid, and the global integral on the destination grid will be as close to area_src as possible given that restriction. You should be able to test this out with the results you have so far... the following should be true:

1. the sum of the product of a field in tos_Omon_CESM2_omip2_r1i1p1f1_gn_024501-030512.nc times the area should match the sum of the same product in your remapped output (there might numerical differences on the order of round-off; 1e-6 for single precision fields and 1e-12 for double-precision fields)
2. the sum of that product in tos_Omon_CESM2_omip2_r1i1p1f1_gr_024501-030512.nc should be close to zero but several orders of magnitude larger than what is observed in [1]

I'll spend some time this coming week trying to find the exact mapping file used, though it might be the case that the modification to the map was done in memory without ever being written to disk. If that's the case, I should be able to create a map for you and make it available.