Advice on implementing sorting prior to regrid-remap (to allow diagnostic grids to be built for non-monotonic tracers)

[gmacgilchrist]

Intention/scope

I am planning to implement the capacity to build diagnostic grids for tracers that are non-monotonic in thee vertical dimension. The purpose is to be able to output diagnostics on a grid of any arbitrary tracer, e.g. temperature. I am proposing that this be used for diagnostic regrid-remap only, not for regrid-remap of the prognostic model fields and grid, for which the feature would be switched off.

I am implementing this initially for the rho coordinate type, and will subsequently extend it to anything in the tracer registry. I am only implementing this at the moment for diagnostics on the tracer points, not yet on the u, v, or w points.

Proposed procedure

To do this requires that the tracer, e.g. rho, has the option to be sorted prior to the regrid step. The regrid will then derive interface heights and thicknesses associated with the sorted column. This will be implemented with a simple sorting algorithm in build_*_column (the algorithm itself would be housed in MOM_regridding.F90. A boolean held in the coordinate control structure, e.g. needs_sorting=.true., could be used to turn this on and off as required.

To appropriately remap to this new grid requires that the model field is also sorted prior to its remap. The most efficient way to do this would be to create an array of the sorting indices and pass it to remapping_core_h. A key question is where and how to build this array. One option would be to create it as a 3d array at the same time as the 3d array of new grid thickness is created (i.e. in diag_remap_update for diagnostic grids and build_*_grid for prognostic grids), and carry this array in the remap control structure. This would mean additional 3d array(s) being carried around in memory.

I am posting this here to invite comments/thoughts/opinions before I start implementing this. Any thoughts welcome.

[gmacgilchrist]

Additional thoughts/intentions are here.

[gmacgilchrist]

Posted with apologies for the ramble: I am still trying to get my head around how best to do things in Fortran (and how people prefer to do things in MOM6)

I wonder if someone can advise me on the best strategy for having variables (including a logical to specify whether a given coordinate tracer needs to be sorted prior to regridding, and a 3D array of indices that would specify the sorting for that coordinate) available across different components of the regrid-remap algorithm.

More specifically, the logical needs to be available in coord_rho.F90 for calling in the build_rho_column subroutine to determine whether the column should be sorted, and whether to return the indices for the sorted column. Then, it needs to be available in MOM_diag_remap.F90 (for diagnostic coordinates) and MOM_regridding.F90 (for prognostics coordinates) to be called within diag_remap_update and build_rho_grid respectively, to ascertain whether a 3D array of indices should be filled. Finally, it needs to be available to remapping_core_h in MOM_remapping.F90 to determine whether the u0 variable should be sorted prior to remapping.

The 3D array needs to be available for filling in MOM_diag_remap.F90 (for diagnostic coordinates) and MOM_regridding.F90 (for prognostics coordinates) as noted above. It then needs to be available in MOM_ALE.F90 to be called by ALE_remap_tracers (prognostic coordinate) and MOM_diag_remap.F90 to be called by diag_remap_do_remap, where in each case it can be looped through and passed to remapping_core_h.

Perhaps it is not necessary to have both variables available in all of these locations. For example, after a 1d array of indices is created by build_rho_column, the presence of this array could be used to determine whether a 3d array should be built, the presence of which, in turn, would signify whether variables should be sorted prior to remapping.

Given the various places where these variables need to be accessed and amended, it seems to make sense to define them in a control structure. But which one? Both the regrid and the remap control structures are present in the diag_remap control structure and the ALE control structure. So perhaps holding the logical in the regrid control structure, and the array in the remap control structure would be the right approach?

[gmacgilchrist]

For the case of the diagnostic grids, I settled on placing ksort3d in the diag_remap_ctrl control structure - this seems to make sense as it is the same location that the destination thicknesses are stored. This means that it is possible to fill it in the diag_remap_update subroutine, and use it in the diag_remap_do_remap subroutine.

needs_sorting is held in the regridding_CS and passed to the rho_CS, which is consistent with the approach for the ref_pressure variable.

The whole case is a little trickier for the prognostic grid. Since this is clearly not a desirable feature for the prognostic grid, I wonder if it should just be made impossible to implement, i.e. needs_sorting is always set to false and ksort3d is never generated?

[gmacgilchrist]

I have a pretty hacky version of this "working", by which I mean that it seems to be visiting all of the places that it should be visiting, and passing arrays where it should be passing them. The modified version is here. The code slows down considerably, presumably because it is calling the sorting algorithm multiple times during each timestep (every time the diagnostic grid is updated). I am holding off submitting a PR because there is much that can be improved, but I will likely need further guidance on the best approach.

Unfortunately, none of the test cases that I can run on my laptop actually create a non-monotonic rho0 profile, so I can't yet see if the sorting is doing what I think it should. Indeed, creating a test case in which that will happen might be tricky given the model's desire to retain static stability. So, it may be easiest to press on and create a T or S diagnostic grid, and design a test case around this. Alternatively, I could run this on a global simulation, where rho0 is likely to be non-monotonic at depth.