Regrid/remap for arbitrary, non-monotonic tracers

[gmacgilchrist]

This is a discussion on the approach to implement regrid/remap capabilities for any arbitrary tracer in the model. Presently, the regrid/remp procedure requires (assumes) that the coordinate to which the grid is regridded is montonically increasing. This forgoes the possibility of regridding to, for example, temperature.

The implementation of this feature will be carried out in four stages:

1. implement the sorting potential for the rho coordinate
2. copy this procedure and implement for temperature
3. generalize this to work for any tracer
4. integrate to work seamlessly with the diagnostics module

[gmacgilchrist]

Implementing sorting for rho coordinate

(for comment from @ashao; there are aspects of this that I'm not sure of, and sure that this approach can be improved, but thought it best to post anyway as a starting point; I am mostly unsure how to handle the different ways that the regrid-remap is done in the diagnostic and prognostic implementations)

Preliminary note

There are some differences between the algorithmic approach in the prognostic model remapping and the diagnostic grid remapping. Both approaches build the new grid one column at a time with build_rho_column, but the diagnostic approach has its own implementation for filling out a 3d array with the new thicknesses (see diag_remap_update in MOM_diag_remap). The prognostic model remapping uses build_rho_grid from MOM_regridding.F90 to build a 3d array of interface differences between old and new model grids (this is turned back to a new thickness array via calc_h_new_by_dz). This difference in approach means that it is not completely clear where to build and store the 3d array of indices for sorting variables prior to the remapping. One could feasibly do this within build_rho_grid or within diag_remap_update, but in each case it would mean it was not available for the other procedure. Below, I'm offering one approach that is focused on building it into the diagnostic procedure. Another approach would be to reorganize the diagnostic approach to more closely align with the prognostic approach (i.e. implement build_rho_grid in diag_remap_update).

Pseudo-algorithm

1. write a subroutine for sorting of 1d column, implemented in MOM_regridding.F90, returns sorted tracer column and sorted index (ksort). @ashao, I think you had previously noted that the sorted tracer was not required to be returned here, but I think that when the new grid is built, it needs to be done for the sorted tracer, thus making sense to return that array here...?
2. sorting called within build_rho_column (perhaps based on boolean, e.g. needs_sorting, present in rho_CS) on densities column before z_interface evaluated. Return z_interface as normal, return ksort 1d array ( in rho_CS? Or as new 1d array?)
3. (Procedure splits because of different approaches) In diag_remap_update, get ksort for each column and fill a 3d array, ksort3d.
4. Via diag_update_remap_grids in MOM_diag_mediator.F90, pass ksort3d to diag_remap_cs (just as is done for h_target).
5. In remapping_core_h, as called in diag_remap_do_remap, use ksort3d (now present in diag_remap_cs) to sort h0 and u0 prior to remapping.

[gmacgilchrist]

Collection of thoughts/issues

Updated 4/18/23

Primary

• Where to specify a needs_sorting logical so that it is accessible everywhere? (don't wish to do any sorting if it is not necessary or desired) Seems to make sense to place in the coordinate CS, but how to make it available everywhere it is needed (getting needs_sorting is PRIVATE ... errors when I try to call it in diag_remap_update for example). I guess it needs to show up in all of the coordinate, remapping, and regridding control structures? What's the best way to pass it between all of these...
• Where to place scalar_sort_k_1d? Currently in regrid_interp.F90, which doesn't feel like the right place, but had some issues calling it from MOM_regridding.F90 into coord_rho.F90.
• Specifying ksort3d in the remap CS feels cumbersome. What are the alternatives to this?

Secondary

• How to make logical allocatable at run time (i.e. how to put option into MOM_input)?

[gmacgilchrist]

I navigated the needs_sorting challenge in a rather hacky way : by writing a check_if_sorting_needed function in the MOM_regridding module. I can call this from another module and retrieve the logical value. This seems like a pretty hacky approach.

[gmacgilchrist]

I've managed to implement the capacity for sorting within remapping_core_h, but still unsure where and how to store the ksort array and its 3d counterpart.

[gmacgilchrist]

Latest perspective on procedure

[ Presently, the code is compiling with all of the following components added, but is failing when attempting to run a test case. I may start over, adding these elements piece by piece until I can work out where the error is. ]

1. Add needs_sorting logical to regrid_cs (MOM_regridding) and rho_cs (coord_rho); add procedure for reading this variable from MOM_input (MOM_regridding) and posting to regrid_cs via set_regrid_params and rho_cs via init_coord_rho.
2. Write subroutine for bubble sort algorithm (note : I'm still not completely sure where to put this subroutine).
3. Implement calling of sorting algorithm within build_rho_column, based on needs_sorting parameter.
4. Add ksort3d array to diag_remap_ctrl control structure; copy structure for h.
5. Implement filling of ksort3d via diag_remap_update (MOM_diag_remap).
6. Write subroutine to apply sorting from ksort, place in MOM_remapping.
7. Implement calling of sorting via ksort within remapping_core_h (MOM_remapping).
8. Implement looping of ksort3d to pass ksort column to remapping_core_h in diag_remap_do_remap.

I think that these stages are currently implemented correctly, but that there are issues with how and where arrays are allocated, the order in which that is happening, and whether or not they are allocated at any given moment.

Current error

In the current configuration, I am getting an error Array bound mismatch for dimension 1 of array 'ksort' at the point where I am trying to index a 1d ksort from ksort3d. This is within diag_remap_do_remap just prior to the call to remapping_core_h. Importantly, it seems that this is being called prior to any call to diag_remap_update, such that ksort3d is actually not going to be filled at that moment (indeed, writing the size of remap_cs%ksort3d returns 0). I need to understand the order of these calls better. Perhaps I could put an if statement to discern whether ksort3d has yet been allocated, and if not, don't try to use it.

[gmacgilchrist]

Quick note before I finish for the week : applying a logical check on whether ksort3d is allocated allows the code to run, but it seems that diag_remap_update is never called, so ksort3d never gets allocated.

[hdrake]

Hi @gmacgilchrist, do you have any updates on this? Is there a PR I can take a look at?

[gmacgilchrist]

https://github.com/gmacgilchrist/MOM6/tree/gmac_diagmods

This is the branch in which I implement a sorting algorithm and build a diagnostic grid based on the sorted tracer (here it is done for RHO). It does not yet include the capacity to build a grid based on any tracer in the model, e.g. temperature.

[gmacgilchrist]

Also, see here for more details: NOAA-GFDL#351

[gmacgilchrist]

A working version of this approach, for temperature only, is here. However, in completing this, I recognised that a column-sorting approach introduces fundamental errors. Specifically, the division of a sorted source column into the coordinate of the tracer is not the same as a division of the original source column into the coordinate of that same tracer.

Consider, for example, a tracer that has a local mid-depth maximum, such that there are two locations where a given value of tracer intersects the column. Each of the cells that contain that value should be partitioned appropriately when transforming the column to a coordinate of that tracer. However, in the sorted column, there is only one intersection of this tracer value, based on whichever cell has the greater tracer value at its center. The other cell, therefore, is not partitioned at all, and included entirely within the layer of the new coordinate for which its center tracer value corresponds.

The core point is that a histogram of water column quantities in tracer space is not the same as a histogram of the same column when it is first sorted based on the tracer value. This seems rather obvious when written explicitly.

However, the regrid-remap architecture in MOM6 relies intrinsically on the coordinate of the target grid being monotonically increasing. It is only possible to assure this by applying a sorting step. Thus, I am currently of the opinion that the regrid-remap approach cannot be used to conservatively bin diagnostics into tracer space when the tracer is nonmonotonic. I could be missing some smarter way of doing this that does not actually sort the column during the regrid step, but I can't see this at the moment.

Absent this smarter way of doing things, I am attempting to implement a histogramming approach. It will copy the basic structure of diagnostic remapping, but rather than finding a new grid and remapping based on this grids, it will determine weights for each source grid cell to be included in each target grid cell.

[hdrake]

@gmacgilchrist, can you point me towards any examples (e.g. MOM6-examples configurations) where you apply the column-sorting to temperature-coordinate diagnostics? I found this MOM6-examples/gmac_fds_thermo branch of yours, but it looks like that only considers the rho0 case.

I'm not sure I understand your point about the histograms. It seems to me that sorting should not affect the histogram.

Also, if you have a branch of the new histogram-based approach you're trying out, I am happy to take a look at that too.

[gmacgilchrist]

@hdrake MOM6-examples folder for remapping to a scalar coordinate here and for invoking the histogram approach (just adding a single line in the MOM_override2 file) here.

[gmacgilchrist]

FYI, the relevant MOM6 branch for doing the histogramming (rather than a sort-remap approach) is here.

[hdrake]

Hi @gmacgilchrist, I finally tried checking out your implementation of these temperature-remapping diagnostics to run them myself. However, I am getting compiler errors, I think because I am using incompatible commits.

What I tried to do is to checkout these two branches locally:
https://github.com/gmacgilchrist/MOM6/tree/gmac_addhistogramandscalarcoord
https://github.com/gmacgilchrist/MOM6-examples/tree/gmac_fds_thermo

The problem seems to be that the submodules linked in MOM6-examples/gmac_fds_thermo are not consistent with MOM6/gmac_addhistogramandscalarcoord. For instance, the coord_scalar.F90 source code file that you introduce in MOM6 is not actually present in the linked commit of the MOM6 submodule.

Relatedly, your branch https://github.com/gmacgilchrist/MOM6/tree/gmac_addhistogramandscalarcoord has a confusing commit history, because your specific changes that enable the temperature remapping seem to be mixed in with thousands of other commits in this branch's commit history. Is there any way you can refactor this branch's history to isolate just the changes you've made?

It would be helpful if you could let me know precisely how you've configured your compile-time environment so I can replicate these results and potentially contribute to this development!

[hdrake]

Maybe it would help to make https://github.com/gmacgilchrist/MOM6/tree/gmac_addhistogramandscalarcoord into a draft PR so the changes you've made can be more easily isolated and commented on?

[gmacgilchrist]

@hdrake Apologies for the confusion - submodules continue to confound me.

I worked on the scalar coordinate implementation outside of the MOM6-examples repository, so they're not intrinsically linked. Perhaps that was an error - I'm not actually sure how to work on independent source code changes from within the MOM6-examples repo in a way that retains them when the parent repo is updated. But anyway, a way forward in the meantime would be to compile a version of the scalar coordinate MOM6 code outside of the MOM-examples repo, and point to that executable when executing from the run folder within MOM6-examples. That was how I was working while testing the code. As I say, that may not be the preferred approach but it was easier to me than wrangling with submodules.

I will work on cleaning up the commit history of the scalar coordinate code - I'm not sure what happened here, something in the way that I rebased the repository meant that it merged the commit history too.

[hdrake]

I've cleaned up the commit history in this branch/PR: NOAA-GFDL#803

I don't intend for this to be the actual PR that gets merged, but creating it helped me sort out the commit history for myself.

[gmacgilchrist]

Description of implementation of histogram approach invoked in mom-ocean/MOM6@main...gmacgilchrist:MOM6:gmac_addhistogramandscalarcoord.

Histogram procedure

Core functions

The core functions for deriving histogram weights can be found in regrid_interp.F90:

• get_interface_indices searches within each grid cell to find if a target interface if it is present. If it is not present, the subroutine records whether the grid cell values are less than or greater than the target interface value.
• get_histogram_weights returns a 2D array that maps the grid onto the target grid. It is implemented by looping over the target grid values, detecting their presence (or absence) in each native grid cell (using get_interface_indices), and determining the fraction of each native grid cell that falls within each target grid layer.
• build_histogram_weights follows the build_and_interpolate_grid procedure to set the polynomials and retrieve the 2D array from get_histogram_weights.

Build (and update) weights

Building the histogram weights for a given coordinate is set up in coord_*.F90 and implemented in MOM_diag_remap:

• The set up is most easily seen in diff for coord_rho.F90, with all aspects carried over to coord_scalar.F90 (see below):
  • Into build_rho_column, the capability to return the 2D array of histogram weights (called histogram_weights) is included (see Core functions above). It is only invoked if histogram_extensive_diags (see below) is True for this coordinate.
• The building of the histogram weights happens in MOM_diag_remap in the subroutine diag_remap_update. The 2D array of weights is derived in the call to build_rho_column (or build_scalar_column) and used to fill out a 4D array of weights.

Implement histogram of extensive diagnostics

Implementing a histogram approach for extensive diagnostic variables is initiated in MOM_diag_mediator.F90, and uses the vertically_histogram_field subroutine (from MOM_diag_remap.F90), which is itself implementing the core function histogram_column (from MOM_remapping.F90):

• The underlying function to histogram an extensive variable is histogram_column (in MOM_remapping.F90), which maps native to target grids using the weights array, simply multiplying the native grid by the weights and summing in the vertical dimension.
• This function is invoked in MOM_diag_remap.F90 by the subroutine vertically_histogram_field, which takes the 4D array of histogram weights as input, loops through i and j and histograms the column. This is wrapped in a diagnostic-specific subroutine vertically_histogram_diag_field.
• The actual implementation of a histogram approach and posting of the diagnostic is done in the post_data_3d subroutine in MOM_diag_mediator, wherein the above function (vertically_histogram_diag_field) is called for the specific diagnostic to be posted.

Infrastructure to implement histogram approach

Various hooks and tags are included to allow for a histogram approach to be invoked:

• histogram_extensive_diags is added as a logical into coordinate control structure, and invocation implemented in set_*_params (where * could be either rho or scalar).
• In MOM_regridding.F90, the procedure for reading the HISTOGRAM_EXTENSIVE_DIAGS parameter from the namelist is implemented. This parameter can then be set as a regrid parameter (which can itself be set for a diagnostic coordinate control structure) using set_regrid_params.
• Also in MOM_regridding.F90, a public subroutine (check_if_histogram_extensive_diags) is included to check whether the histogram_extensive_diags is True (this is to get past the fact that some control structure variables are not public, but the public subroutine allows us to check).

Scalar coordinate

Adding a "scalar" coordinate is done by replicating all of the relevant code for the RHO coordinate, primarily in two locations (coord_*.F90 and MOM_regridding.F90) with some small modifications.

1. Adds coord_scalar.F90 as an exact replica of coord_rho.F90. The file is kept identical for simplicity, with only names changed from rho to scalar. See above for changes implemented to coord_rho.F90 that are carried over.
2. In MOM_regridding.F90, all aspects of the RHO coordinate are copied for a SCALAR coordinate. At present, units and conversions are set as thought SCALAR were temperature.

[gmacgilchrist]

@hdrake In refreshing my memory of the histogram approach, I put together the above, which may be helpful for interpreting the mess of file changes in that commit.

[hdrake]

Thanks @gmacgilchrist, these notes provide some helpful structure for understanding your code changes.