+Add answer_date optional arguments

  Added optional answer_date arguments to various remapping routines.  These are
vintage-encoding integers intended to replace the logical answers_2018
arguments, and allow for multiple generations of improved algorithms rather than
just two choices, without requiring added interface changes.  However, this
change is backward compatible, and these two arguments are both offered for
now.  All answers are bitwise identical, but there are new optional arguments
to numerous publicly visible routines.

src/ALE/P1M_functions.F90

@@ -24,7 +24,7 @@ module P1M_functions
 !!
 !! It is assumed that the size of the array 'u' is equal to the number of cells
 !! defining 'grid' and 'ppoly'. No consistency check is performed here.
-subroutine P1M_interpolation( N, h, u, edge_values, ppoly_coef, h_neglect, answers_2018 )
+subroutine P1M_interpolation( N, h, u, edge_values, ppoly_coef, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(:),   intent(in)    :: h !< cell widths (size N) [H]
   real, dimension(:),   intent(in)    :: u !< cell average properties (size N) [A]
@@ -33,13 +33,15 @@ subroutine P1M_interpolation( N, h, u, edge_values, ppoly_coef, h_neglect, answe
                                            !! piecewise polynomial coefficients, mainly [A]
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
 
   ! Local variables
   integer   :: k            ! loop index
   real      :: u0_l, u0_r   ! edge values (left and right)
 
   ! Bound edge values (routine found in 'edge_values.F90')
-  call bound_edge_values( N, h, u, edge_values, h_neglect, answers_2018 )
+  call bound_edge_values( N, h, u, edge_values, h_neglect, &
+                          answers_2018=answers_2018, answer_date=answer_date )
 
   ! Systematically average discontinuous edge values (routine found in
   ! 'edge_values.F90')

src/ALE/P3M_functions.F90

@@ -25,7 +25,7 @@ module P3M_functions
 !!
 !! It is assumed that the size of the array 'u' is equal to the number of cells
 !! defining 'grid' and 'ppoly'. No consistency check is performed here.
-subroutine P3M_interpolation( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, answers_2018 )
+subroutine P3M_interpolation( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(:),   intent(in)    :: h !< cell widths (size N) [H]
   real, dimension(:),   intent(in)    :: u !< cell averages (size N) in arbitrary units [A]
@@ -35,13 +35,15 @@ subroutine P3M_interpolation( N, h, u, edge_values, ppoly_S, ppoly_coef, h_negle
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width for the
                                           !! purpose of cell reconstructions [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
 
   ! Call the limiter for p3m, which takes care of everything from
   ! computing the coefficients of the cubic to monotonizing it.
   ! This routine could be called directly instead of having to call
   ! 'P3M_interpolation' first but we do that to provide an homogeneous
   ! interface.
-  call P3M_limiter( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, answers_2018 )
+  call P3M_limiter( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, &
+                    answers_2018=answers_2018, answer_date=answer_date )
 
 end subroutine P3M_interpolation
 
@@ -58,7 +60,7 @@ end subroutine P3M_interpolation
 !!    c. If not, monotonize cubic curve and rebuild it
 !!
 !! Step 3 of the monotonization process leaves all edge values unchanged.
-subroutine P3M_limiter( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, answers_2018 )
+subroutine P3M_limiter( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(:),   intent(in)    :: h !< cell widths (size N) [H]
   real, dimension(:),   intent(in)    :: u !< cell averages (size N) in arbitrary units [A]
@@ -68,6 +70,7 @@ subroutine P3M_limiter( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, an
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width for
                                            !! the purpose of cell reconstructions [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
 
   ! Local variables
   integer :: k            ! loop index
@@ -86,7 +89,7 @@ subroutine P3M_limiter( N, h, u, edge_values, ppoly_S, ppoly_coef, h_neglect, an
   eps = 1e-10
 
   ! 1. Bound edge values (boundary cells are assumed to be local extrema)
-  call bound_edge_values( N, h, u, edge_values, hNeglect, answers_2018 )
+  call bound_edge_values( N, h, u, edge_values, hNeglect, answers_2018=answers_2018, answer_date=answer_date )
 
   ! 2. Systematically average discontinuous edge values
   call average_discontinuous_edge_values( N, edge_values )

src/ALE/PPM_functions.F90

@@ -25,21 +25,22 @@ module PPM_functions
 contains
 
 !> Builds quadratic polynomials coefficients from cell mean and edge values.
-subroutine PPM_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect, answers_2018)
+subroutine PPM_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect, answers_2018, answer_date)
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(N),   intent(in)    :: h !< Cell widths [H]
   real, dimension(N),   intent(in)    :: u !< Cell averages [A]
   real, dimension(N,2), intent(inout) :: edge_values !< Edge values [A]
   real, dimension(N,3), intent(inout) :: ppoly_coef !< Polynomial coefficients, mainly [A]
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
 
   ! Local variables
   integer   :: k              ! Loop index
   real      :: edge_l, edge_r ! Edge values (left and right)
 
   ! PPM limiter
-  call PPM_limiter_standard( N, h, u, edge_values, h_neglect, answers_2018 )
+  call PPM_limiter_standard( N, h, u, edge_values, h_neglect, answers_2018=answers_2018, answer_date=answer_date )
 
   ! Loop over all cells
   do k = 1,N
@@ -59,13 +60,14 @@ end subroutine PPM_reconstruction
 !> Adjusts edge values using the standard PPM limiter (Colella & Woodward, JCP 1984)
 !! after first checking that the edge values are bounded by neighbors cell averages
 !! and that the edge values are monotonic between cell averages.
-subroutine PPM_limiter_standard( N, h, u, edge_values, h_neglect, answers_2018 )
+subroutine PPM_limiter_standard( N, h, u, edge_values, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(:),   intent(in)    :: h !< cell widths (size N) [H]
   real, dimension(:),   intent(in)    :: u !< cell average properties (size N) [A]
   real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A]
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
 
   ! Local variables
   integer   :: k              ! Loop index
@@ -74,7 +76,7 @@ subroutine PPM_limiter_standard( N, h, u, edge_values, h_neglect, answers_2018 )
   real      :: expr1, expr2
 
   ! Bound edge values
-  call bound_edge_values( N, h, u, edge_values, h_neglect, answers_2018 )
+  call bound_edge_values( N, h, u, edge_values, h_neglect, answers_2018=answers_2018, answer_date=answer_date )
 
   ! Make discontinuous edge values monotonic
   call check_discontinuous_edge_values( N, u, edge_values )

src/ALE/PQM_functions.F90

@@ -17,7 +17,7 @@ module PQM_functions
 !!
 !! It is assumed that the dimension of 'u' is equal to the number of cells
 !! defining 'grid' and 'ppoly'. No consistency check is performed.
-subroutine PQM_reconstruction( N, h, u, edge_values, edge_slopes, ppoly_coef, h_neglect, answers_2018 )
+subroutine PQM_reconstruction( N, h, u, edge_values, edge_slopes, ppoly_coef, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(:),   intent(in)    :: h !< cell widths (size N) [H]
   real, dimension(:),   intent(in)    :: u !< cell averages (size N) [A]
@@ -27,6 +27,7 @@ subroutine PQM_reconstruction( N, h, u, edge_values, edge_slopes, ppoly_coef, h_
   real,       optional, intent(in)    :: h_neglect  !< A negligibly small width for
                                            !! the purpose of cell reconstructions [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
 
   ! Local variables
   integer   :: k                ! loop index
@@ -36,7 +37,7 @@ subroutine PQM_reconstruction( N, h, u, edge_values, edge_slopes, ppoly_coef, h_
   real      :: a, b, c, d, e    ! parabola coefficients
 
   ! PQM limiter
-  call PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answers_2018 )
+  call PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answers_2018=answers_2018, answer_date=answer_date )
 
   ! Loop on cells to construct the cubic within each cell
   do k = 1,N
@@ -72,7 +73,7 @@ end subroutine PQM_reconstruction
 !!
 !! It is assumed that the dimension of 'u' is equal to the number of cells
 !! defining 'grid' and 'ppoly'. No consistency check is performed.
-subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answers_2018 )
+subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(:),   intent(in)    :: h !< cell widths (size N) [H]
   real, dimension(:),   intent(in)    :: u !< cell average properties (size N) [A]
@@ -81,6 +82,7 @@ subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answers_20
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width for
                                            !! the purpose of cell reconstructions [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
 
   ! Local variables
   integer :: k            ! loop index
@@ -102,7 +104,7 @@ subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answers_20
   hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect
 
   ! Bound edge values
-  call bound_edge_values( N, h, u, edge_values, hNeglect, answers_2018 )
+  call bound_edge_values( N, h, u, edge_values, hNeglect, answers_2018=answers_2018, answer_date=answer_date )
 
   ! Make discontinuous edge values monotonic (thru averaging)
   call check_discontinuous_edge_values( N, u, edge_values )

src/ALE/regrid_edge_values.F90

@@ -41,14 +41,16 @@ module regrid_edge_values
 !! Both boundary edge values are set equal to the boundary cell averages.
 !! Any extrapolation scheme is applied after this routine has been called.
 !! Therefore, boundary cells are treated as if they were local extrama.
-subroutine bound_edge_values( N, h, u, edge_val, h_neglect, answers_2018 )
+subroutine bound_edge_values( N, h, u, edge_val, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(N),   intent(in)    :: h !< cell widths [H]
   real, dimension(N),   intent(in)    :: u !< cell average properties in arbitrary units [A]
   real, dimension(N,2), intent(inout) :: edge_val !< Potentially modified edge values [A]; the
                                            !! second index is for the two edges of each cell.
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
+
   ! Local variables
   real    :: sigma_l, sigma_c, sigma_r    ! left, center and right van Leer slopes [A H-1] or [A]
   real    :: slope_x_h     ! retained PLM slope times  half grid step [A]
@@ -57,6 +59,7 @@ subroutine bound_edge_values( N, h, u, edge_val, h_neglect, answers_2018 )
   integer :: k, km1, kp1   ! Loop index and the values to either side.
 
   use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018
+  if (present(answer_date)) use_2018_answers = (answer_date < 20190101)
   if (use_2018_answers) then
     hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect
   endif
@@ -218,14 +221,15 @@ end subroutine edge_values_explicit_h2
 !! available interpolant.
 !!
 !! For this fourth-order scheme, at least four cells must exist.
-subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answers_2018 )
+subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(N),   intent(in)    :: h !< cell widths [H]
   real, dimension(N),   intent(in)    :: u !< cell average properties in arbitrary units [A]
   real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the second index
                                            !! is for the two edges of each cell.
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
 
   ! Local variables
   real :: h0, h1, h2, h3        ! temporary thicknesses [H]
@@ -247,6 +251,7 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answers_2018 )
   logical   :: use_2018_answers  ! If true use older, less acccurate expressions.
 
   use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018
+  if (present(answer_date)) use_2018_answers = (answer_date < 20190101)
   if (use_2018_answers) then
     hNeglect = hNeglect_edge_dflt ; if (present(h_neglect)) hNeglect = h_neglect
   else
@@ -382,14 +387,15 @@ end subroutine edge_values_explicit_h4
 !!
 !! There are N+1 unknowns and we are able to write N-1 equations. The
 !! boundary conditions close the system.
-subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answers_2018 )
+subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(N),   intent(in)    :: h !< cell widths [H]
   real, dimension(N),   intent(in)    :: u !< cell average properties in arbitrary units [A]
   real, dimension(N,2), intent(inout) :: edge_val !< Returned edge values [A]; the second index
                                            !! is for the two edges of each cell.
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
 
   ! Local variables
   integer               :: i, j                 ! loop indexes
@@ -418,6 +424,7 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answers_2018 )
   logical   :: use_2018_answers  ! If true use older, less acccurate expressions.
 
   use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018
+  if (present(answer_date)) use_2018_answers = (answer_date < 20190101)
   if (use_2018_answers) then
     hNeglect = hNeglect_edge_dflt ; if (present(h_neglect)) hNeglect = h_neglect
   else
@@ -690,14 +697,16 @@ end subroutine end_value_h4
 !!
 !! There are N+1 unknowns and we are able to write N-1 equations. The
 !! boundary conditions close the system.
-subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answers_2018 )
+subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(N),   intent(in)    :: h !< cell widths [H]
   real, dimension(N),   intent(in)    :: u !< cell average properties in arbitrary units [A]
   real, dimension(N,2), intent(inout) :: edge_slopes !< Returned edge slopes [A H-1]; the
                                            !! second index is for the two edges of each cell.
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
+
   ! Local variables
   integer               :: i, j                 ! loop indexes
   real                  :: h0, h1               ! cell widths [H or nondim]
@@ -729,6 +738,7 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answers_201
   hNeglect = hNeglect_dflt ; if (present(h_neglect))  hNeglect = h_neglect
   hNeglect3 = hNeglect**3
   use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018
+  if (present(answer_date)) use_2018_answers = (answer_date < 20190101)
 
   ! Loop on cells (except last one)
   do i = 1,N-1
@@ -859,14 +869,16 @@ end subroutine edge_slopes_implicit_h3
 
 !------------------------------------------------------------------------------
 !> Compute ih5 edge slopes (implicit fifth order accurate)
-subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answers_2018 )
+subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(N),   intent(in)    :: h !< cell widths [H]
   real, dimension(N),   intent(in)    :: u !< cell average properties in arbitrary units [A]
   real, dimension(N,2), intent(inout) :: edge_slopes !< Returned edge slopes [A H-1]; the
                                            !! second index is for the two edges of each cell.
   real, optional,       intent(in)    :: h_neglect !< A negligibly small width [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
+
 ! -----------------------------------------------------------------------------
 ! Fifth-order implicit estimates of edge slopes are based on a four-cell,
 ! three-edge stencil. A tridiagonal system is set up and is based on
@@ -1129,14 +1141,15 @@ end subroutine edge_slopes_implicit_h5
 !!          become computationally expensive if regridding is carried out
 !!          often. Figuring out closed-form expressions for these coefficients
 !!          on nonuniform meshes turned out to be intractable.
-subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answers_2018 )
+subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answers_2018, answer_date )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(N),   intent(in)    :: h !< cell widths [H]
   real, dimension(N),   intent(in)    :: u !< cell average properties (size N) in arbitrary units [A]
   real, dimension(N,2), intent(inout) :: edge_val  !< Returned edge values [A]; the second index
                                            !! is for the two edges of each cell.
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width [H]
   logical,    optional, intent(in)    :: answers_2018 !< If true use older, less acccurate expressions.
+  integer,    optional, intent(in)    :: answer_date  !< The vintage of the expressions to use
 
   ! Local variables
   real :: h0, h1, h2, h3       ! cell widths [H]