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.