Refactor of vertical reconstruction with six new schemes

This refactor is motivated by i) adding a set of new reconstructions schemes
that are ready-to-use-now for remapping, but ii) will be an essential part of a
new grid generator that will use some of the new member methods available only
as a result of the refactor. The "refactor" actually duplicates many existing
schemes, leaving the original code in place to provide a period of transition
just in case there's some major reason to revert to the old style of code.

What changed
============
- added 16 modules name Recon1d_*.F90, which extend a base Recon1d_type.F90
- MOM_remapping.F90
  - rewrote remapping_unit_tests to be tidier
  - added tests for
    1) invariance when remapping to/from the same grid;
    2) preservation of a uniform value (many of the old schemes fail);
    3) internal consistency such as monotonic ordering of values at
       edges, no internal extrema, etc., but only for the new Recon1d_* schemes.
- added more schemes to config_src/drivesr/timing_tests/time_MOM_remapping.F90
- added documentation via _Vertical_Reconstruction.dox with references

Six new schemes are added: C_PLM_CW, C_PLM_CWK, C_MPLM_CWK, C_EMPLM_CWK,
C_PPM_CWK, and C_EPPM_CWK.

Ten existing schemes were ported to the new class, along with their
non-monotonic behaviors or ill-formed expressions.  Others can be ported later
if needed, although one result of the refactor was the discovery that none of
the existing schemes worked as desired.

The new form uses a class (Recon1d) to define a uniform API for all the
reconstruction schemes. The methods of the class include:
- init() that creates work arrays and stores parameters such as `h_neglect`
- reconstruct() that causes the reconstruction parameters to be calculated
- average() that returns the average between two points in a reconstruction
  (used in the remap_src_to_sub_grid() functions)
- unit_tests() method
- check_reconstruction() checks that a reconstruction has the
  properties that it should (e.g. reconstruction values are monotonic)
  Theses properties are defined by each scheme, although many are shared
  between schemes.
- f() that evaluates the reconstruction at a point (this will be used in the
  new grid-generator).
- dfdx() that evaluates the derivative of the reconstruction at a point (this
  will be used in the new grid-generator)

Testing
=======

New tests have been added that check uniform states are preserved, states are
preserved under unchanging grid, and that the reconstructions are internally
consistent according to their own definition. This last is only available for
the new class-based schemes via the check_reconstruction() method which can be
customized to each reconstruction, which allows us to document which properties
(e.g. monotonicity) are satisfied and which are not.  These new tests use
seeded random numbers which has proven useful for debugging/development, but
are not necessarily useful in the CI. Nevertheless, I've made it part of
remapping_unit_tests() so there is a template for debugging/testing new schemes
when we add them. The same sequence of random numbers is used for each scheme,
but the numbers are compiler dependent so the documented "fails" record the
compiler version in comments.

The helper type was previously broken out of MOM_remapping to make testing convenient.

Performance
===========
- the new class-based schemes are all faster than their counterparts in the
  original code style;
- example numbers can be found in the build-test-perfmon job  of this commit
  (expand the "Display timing results");
- e.g. timings for PPM_HYBGEN and C_PPM_HYBGEN show a marginal (~10%) speed
  up due to code style, and C_PPM_CWK has a more substantial speed up (~20%)
  due to a lighter algorithm;
- I imagine we will see further speed ups if we worked on arrays rather than
  columns, which is an extension that could be implemented later.

Documentation
=============
- the new module APIs are documented;
- a new page, `_Vertical_Reconstruction.dox`, details the vertical
  reconstruction schemes, tabulates the values of `REMPAPPING_SCHEME`,
  and summarize some properties.

Minor (unrelated)
=================
- removed `-j` for nested make in `.testing/Makefile`

Author: adcroft

.testing/Makefile

@@ -267,16 +267,16 @@ $(BUILD)/timing/Makefile: MOM_ACFLAGS += --with-driver=timing_tests
 # Build executables
 .NOTPARALLEL:$(foreach e,$(UNIT_EXECS),$(BUILD)/unit/$(e))
 $(BUILD)/unit/test_%: $(BUILD)/unit/Makefile FORCE
-	cd $(@D) && $(TIME) $(MAKE) $(@F) -j
+	cd $(@D) && $(TIME) $(MAKE) $(@F)
 $(BUILD)/unit/Makefile: $(foreach e,$(UNIT_EXECS),../config_src/drivers/unit_tests/$(e).F90)
 
 .NOTPARALLEL:$(foreach e,$(TIMING_EXECS),$(BUILD)/timing/$(e))
 $(BUILD)/timing/time_%: $(BUILD)/timing/Makefile FORCE
-	cd $(@D) && $(TIME) $(MAKE) $(@F) -j
+	cd $(@D) && $(TIME) $(MAKE) $(@F)
 $(BUILD)/timing/Makefile: $(foreach e,$(TIMING_EXECS),../config_src/drivers/timing_tests/$(e).F90)
 
 $(BUILD)/%/MOM6: $(BUILD)/%/Makefile FORCE
-	cd $(@D) && $(TIME) $(MAKE) $(@F) -j
+	cd $(@D) && $(TIME) $(MAKE) $(@F)
 
 # Target codebase should use its own build system
 $(BUILD)/target/MOM6: $(BUILD)/target FORCE | $(TARGET_CODEBASE)

config_src/drivers/timing_tests/time_MOM_remapping.F90

@@ -9,8 +9,30 @@ program time_MOM_remapping
 implicit none
 
 type(remapping_CS) :: CS
-integer, parameter :: nk=75, nij=20*20, nits=10, nsamp=100, nschemes = 2
-character(len=10) :: scheme_labels(nschemes)
+integer, parameter :: nk=75, nij=20*20, nits=10, nsamp=100, nschemes = 22
+character(len=16) :: scheme_labels(nschemes) = [ character(len=16) :: &
+      'PCM', &
+      'C_PCM', &
+      'PLM', &
+      'C_MPLM_WA', &
+      'C_EMPLM_WA', &
+      'C_PLM_HYBGEN', &
+      'C_PLM_CW', &
+      'C_PLM_CWK', &
+      'C_MPLM_WA_POLY', &
+      'C_EMPLM_WA_POLY', &
+      'C_MPLM_CWK', &
+      'PPM_H4', &
+      'PPM_IH4', &
+      'PQM_IH4IH3', &
+      'PPM_CW', &
+      'PPM_HYBGEN', &
+      'C_PPM_H4_2018', &
+      'C_PPM_H4_2019', &
+      'C_PPM_HYBGEN', &
+      'C_PPM_CW', &
+      'C_PPM_CWK', &
+      'C_EPPM_CWK' ]
 real, dimension(nschemes) :: timings ! Time for nits of nij calls for each scheme [s]
 real, dimension(nschemes) :: tmean ! Mean time for a call [s]
 real, dimension(nschemes) :: tstd ! Standard deviation of time for a call [s]
@@ -31,9 +53,6 @@ program time_MOM_remapping
 seed(:) = 102030405
 call random_seed(put=seed)
 
-scheme_labels(1) = 'PCM'
-scheme_labels(2) = 'PLM'
-
 ! Set up some test data (note: using k,i indexing rather than i,k)
 allocate( u0(nk,nij), h0(nk,nij), u1(nk,nij), h1(nk,nij) )
 call random_number(u0) ! In range 0-1
@@ -61,8 +80,8 @@ program time_MOM_remapping
 do isamp = 1, nsamp
   ! Time reconstruction + remapping
   do ischeme = 1, nschemes
-    call initialize_remapping(CS, remapping_scheme=trim(scheme_labels(ischeme)), &
-                 h_neglect=h_neglect, h_neglect_edge=h_neglect)
+    call initialize_remapping(CS, remapping_scheme=trim(scheme_labels(ischeme)), nk=nk, &
+                              h_neglect=h_neglect, h_neglect_edge=h_neglect)
     call cpu_time(start)
     do iter = 1, nits ! Make many passes to reduce sampling error
       do ij = 1, nij ! Calling nij times to make similar to cost in MOM_ALE()

config_src/drivers/unit_tests/test_MOM_remapping.F90

@@ -2,6 +2,18 @@ program test_MOM_remapping
 
 use MOM_remapping, only : remapping_unit_tests
 
-if (remapping_unit_tests(.true.)) stop 1
+integer :: n !< Number of arguments, or tests
+character(len=12) :: cmd_ln_arg !< Command line argument (if any)
+
+n = command_argument_count()
+
+if (n==1) then
+  call get_command_argument(1, cmd_ln_arg)
+  read(cmd_ln_arg,*) n
+else
+  n = 3000 ! Fallback value if no argument provided
+endif
+
+if (remapping_unit_tests(.true., num_comp_samp=n)) stop 1
 
 end program test_MOM_remapping

docs/discrete_space.rst

@@ -17,4 +17,5 @@ algorithm.
     api/generated/pages/Discrete_Coriolis
     api/generated/pages/Discrete_PG
     api/generated/pages/Energetic_Consistency
+    api/generated/pages/Vertical_Reconstruction
     api/generated/pages/Discrete_OBC

docs/zotero.bib

@@ -2946,3 +2946,17 @@ @article{Young1994
 	pages={1812--1826},
 	year={1994}
 }
+
+@article{van_leer_1977,
+	title = {Towards the ultimate conservative difference scheme. {IV}. {A} new approach to numerical convection},
+	volume = {23},
+	issn = {0021-9991},
+	doi = {10.1016/0021-9991(77)90095-X},
+	number = {3},
+	journal = {Journal of Computational Physics},
+	author = {Van Leer, Bram},
+	month = mar,
+	year = {1977},
+	pages = {276--299},
+}
+

src/ALE/MOM_ALE.F90

@@ -257,15 +257,15 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
     h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10
   endif
 
-  call initialize_remapping( CS%remapCS, string, &
+  call initialize_remapping( CS%remapCS, string, nk=GV%ke, &
                              boundary_extrapolation=init_boundary_extrap, &
                              check_reconstruction=check_reconstruction, &
                              check_remapping=check_remapping, &
                              force_bounds_in_subcell=force_bounds_in_subcell, &
                              om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
                              answer_date=CS%answer_date, &
                              h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
-  call initialize_remapping( CS%vel_remapCS, vel_string, &
+  call initialize_remapping( CS%vel_remapCS, vel_string, nk=GV%ke, &
                              boundary_extrapolation=init_boundary_extrap, &
                              check_reconstruction=check_reconstruction, &
                              check_remapping=check_remapping, &