From 8fd122990afab650b2aa6cfd14a5f24256953cd7 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 29 Jul 2022 12:02:19 -0400 Subject: [PATCH] +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 | 6 ++- src/ALE/P3M_functions.F90 | 11 ++++-- src/ALE/PPM_functions.F90 | 10 +++-- src/ALE/PQM_functions.F90 | 10 +++-- src/ALE/regrid_edge_values.F90 | 25 ++++++++++--- src/ALE/regrid_solvers.F90 | 8 +++- src/framework/MOM_horizontal_regridding.F90 | 41 ++++++++++++++------- 7 files changed, 76 insertions(+), 35 deletions(-) diff --git a/src/ALE/P1M_functions.F90 b/src/ALE/P1M_functions.F90 index d99c611229..281971cca4 100644 --- a/src/ALE/P1M_functions.F90 +++ b/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') diff --git a/src/ALE/P3M_functions.F90 b/src/ALE/P3M_functions.F90 index e3a9f75a3c..4d39542337 100644 --- a/src/ALE/P3M_functions.F90 +++ b/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 ) diff --git a/src/ALE/PPM_functions.F90 b/src/ALE/PPM_functions.F90 index bbf93b4a81..16441565ac 100644 --- a/src/ALE/PPM_functions.F90 +++ b/src/ALE/PPM_functions.F90 @@ -25,7 +25,7 @@ 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] @@ -33,13 +33,14 @@ subroutine PPM_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect, answ 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 ) diff --git a/src/ALE/PQM_functions.F90 b/src/ALE/PQM_functions.F90 index 630ecb34fc..d3809a5d1c 100644 --- a/src/ALE/PQM_functions.F90 +++ b/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 ) diff --git a/src/ALE/regrid_edge_values.F90 b/src/ALE/regrid_edge_values.F90 index a972fc3444..08425fc92d 100644 --- a/src/ALE/regrid_edge_values.F90 +++ b/src/ALE/regrid_edge_values.F90 @@ -41,7 +41,7 @@ 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] @@ -49,6 +49,8 @@ subroutine bound_edge_values( N, h, u, edge_val, h_neglect, answers_2018 ) !! 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,7 +221,7 @@ 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] @@ -226,6 +229,7 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answers_2018 ) !! 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,7 +387,7 @@ 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] @@ -390,6 +395,7 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answers_2018 ) !! 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,7 +697,7 @@ 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] @@ -698,6 +705,8 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answers_201 !! 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,7 +869,7 @@ 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] @@ -867,6 +877,8 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answers_201 !! 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,7 +1141,7 @@ 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] @@ -1137,6 +1149,7 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answers_2018 ) !! 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] diff --git a/src/ALE/regrid_solvers.F90 b/src/ALE/regrid_solvers.F90 index b7cc3b5402..022946a29d 100644 --- a/src/ALE/regrid_solvers.F90 +++ b/src/ALE/regrid_solvers.F90 @@ -16,12 +16,13 @@ module regrid_solvers !! This routine uses Gauss's algorithm to transform the system's original !! matrix into an upper triangular matrix. Back substitution yields the answer. !! The matrix A must be square, with the first index varing down the column. -subroutine solve_linear_system( A, R, X, N, answers_2018 ) +subroutine solve_linear_system( A, R, X, N, answers_2018, answer_date ) integer, intent(in) :: N !< The size of the system real, dimension(N,N), intent(inout) :: A !< The matrix being inverted [nondim] real, dimension(N), intent(inout) :: R !< system right-hand side [A] real, dimension(N), intent(inout) :: X !< solution vector [A] logical, optional, intent(in) :: answers_2018 !< If true or absent use older, less efficient expressions. + integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use ! Local variables real, parameter :: eps = 0.0 ! Minimum pivot magnitude allowed real :: factor ! The factor that eliminates the leading nonzero element in a row. @@ -32,6 +33,7 @@ subroutine solve_linear_system( A, R, X, N, answers_2018 ) integer :: i, j, k old_answers = .true. ; if (present(answers_2018)) old_answers = answers_2018 + if (present(answer_date)) old_answers = (answer_date < 20190101) ! Loop on rows to transform the problem into multiplication by an upper-right matrix. do i = 1,N-1 @@ -173,7 +175,7 @@ end subroutine linear_solver !! !! This routine uses Thomas's algorithm to solve the tridiagonal system AX = R. !! (A is made up of lower, middle and upper diagonals) -subroutine solve_tridiagonal_system( Al, Ad, Au, R, X, N, answers_2018 ) +subroutine solve_tridiagonal_system( Al, Ad, Au, R, X, N, answers_2018, answer_date ) integer, intent(in) :: N !< The size of the system real, dimension(N), intent(in) :: Ad !< Matrix center diagonal real, dimension(N), intent(in) :: Al !< Matrix lower diagonal @@ -181,6 +183,7 @@ subroutine solve_tridiagonal_system( Al, Ad, Au, R, X, N, answers_2018 ) real, dimension(N), intent(in) :: R !< system right-hand side real, dimension(N), intent(out) :: X !< solution vector 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, dimension(N) :: pivot, Al_piv real, dimension(N) :: c1 ! Au / pivot for the backward sweep @@ -189,6 +192,7 @@ subroutine solve_tridiagonal_system( Al, Ad, Au, R, X, N, answers_2018 ) logical :: old_answers ! If true, use expressions that give the original (2008 through 2018) MOM6 answers old_answers = .true. ; if (present(answers_2018)) old_answers = answers_2018 + if (present(answer_date)) old_answers = (answer_date < 20190101) if (old_answers) then ! This version gives the same answers as the original (2008 through 2018) MOM6 code diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 05e3e393b6..bbb5ae0e15 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -81,7 +81,7 @@ end subroutine myStats !> Use ICE-9 algorithm to populate points (fill=1) with valid data (good=1). If no information !! is available, use a previous guess (prev). Optionally (smooth) blend the filled points to !! achieve a more desirable result. -subroutine fill_miss_2d(aout, good, fill, prev, G, acrit, num_pass, relc, debug, answers_2018) +subroutine fill_miss_2d(aout, good, fill, prev, G, acrit, num_pass, relc, debug, answer_date) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. real, dimension(SZI_(G),SZJ_(G)), & intent(inout) :: aout !< The array with missing values to fill [A] @@ -98,8 +98,9 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, acrit, num_pass, relc, debug, integer, optional, intent(in) :: num_pass !< The maximum number of iterations real, optional, intent(in) :: relc !< A relaxation coefficient for Laplacian [nondim] logical, optional, intent(in) :: debug !< If true, write verbose debugging messages. - logical, optional, intent(in) :: answers_2018 !< If true, use expressions that give the same - !! answers as the code did in late 2018. Otherwise + integer, optional, intent(in) :: answer_date !< The vintage of the expressions in the code. + !! Dates before 20190101 give the same answers + !! as the code did in late 2018, while later versions !! add parentheses for rotational symmetry. real, dimension(SZI_(G),SZJ_(G)) :: a_filled ! The aout with missing values filled in [A] @@ -135,7 +136,7 @@ subroutine fill_miss_2d(aout, good, fill, prev, G, acrit, num_pass, relc, debug, relax_coeff = relc_default if (PRESENT(relc)) relax_coeff = relc - ans_2018 = .true. ; if (PRESENT(answers_2018)) ans_2018 = answers_2018 + ans_2018 = .true. ; if (PRESENT(answer_date)) ans_2018 = (answer_date < 20190101) fill_pts(:,:) = fill(:,:) @@ -251,7 +252,7 @@ end subroutine fill_miss_2d !> Extrapolate and interpolate from a file record subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, recnum, G, tr_z, mask_z, & z_in, z_edges_in, missing_value, reentrant_x, tripolar_n, & - homogenize, m_to_Z, answers_2018, ongrid, tr_iter_tol) + homogenize, m_to_Z, answers_2018, ongrid, tr_iter_tol, answer_date) character(len=*), intent(in) :: filename !< Path to file containing tracer to be !! interpolated. @@ -287,6 +288,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, real, optional, intent(in) :: tr_iter_tol !< The tolerance for changes in tracer concentrations !! between smoothing iterations that determines when to !! stop iterating [CU ~> conc] + integer, optional, intent(in) :: answer_date !< The vintage of the expressions in the code. + !! Dates before 20190101 give the same answers + !! as the code did in late 2018, while later versions + !! add parentheses for rotational symmetry. ! Local variables real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its @@ -313,6 +318,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, real :: missing_val_in ! The missing value in the input field [conc] real :: roundoff ! The magnitude of roundoff, usually ~2e-16 [nondim] real :: add_offset, scale_factor ! File-specific conversion factors. + integer :: ans_date ! The vintage of the expressions and order of arithmetic to use logical :: found_attr logical :: add_np logical :: is_ongrid @@ -356,6 +362,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, PI_180 = atan(1.0)/45. + ans_date = 20181231 + if (present(answers_2018)) then ; if (.not.answers_2018) ans_date = 20190101 ; endif + if (present(answer_date)) ans_date = answer_date + ! Open NetCDF file and if present, extract data and spatial coordinate information ! The convention adopted here requires that the data be written in (i,j,k) ordering. @@ -565,8 +575,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, conversion, good2(:,:) = good(:,:) fill2(:,:) = fill(:,:) - call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, dtr_iter_stop, & - answers_2018=answers_2018) + call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, dtr_iter_stop, answer_date=ans_date) if (debug) then call myStats(tr_outf, missing_value, is, ie, js, je, k, 'field from fill_miss_2d()', scale=I_scale) endif @@ -589,7 +598,8 @@ end subroutine horiz_interp_and_extrap_tracer_record !> Extrapolate and interpolate using a FMS time interpolation handle subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, tr_z, mask_z, & z_in, z_edges_in, missing_value, reentrant_x, tripolar_n, & - homogenize, spongeOngrid, m_to_Z, answers_2018, tr_iter_tol) + homogenize, spongeOngrid, m_to_Z, & + answers_2018, tr_iter_tol, answer_date) integer, intent(in) :: fms_id !< A unique id used by the FMS time interpolator type(time_type), intent(in) :: Time !< A FMS time type @@ -621,6 +631,10 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t real, optional, intent(in) :: tr_iter_tol !< The tolerance for changes in tracer concentrations !! between smoothing iterations that determines when to !! stop iterating [CU ~> conc] + integer, optional, intent(in) :: answer_date !< The vintage of the expressions in the code. + !! Dates before 20190101 give the same answers + !! as the code did in late 2018, while later versions + !! add parentheses for rotational symmetry. ! Local variables real, dimension(:,:), allocatable :: tr_in !< A 2-d array for holding input data on its @@ -658,7 +672,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t integer, dimension(4) :: fld_sz logical :: debug=.false. logical :: is_ongrid - logical :: ans_2018 + integer :: ans_date ! The vintage of the expressions and order of arithmetic to use real :: I_scale ! The inverse of the conversion factor for diagnostic output [conc CU-1 ~> 1] real :: dtr_iter_stop ! The tolerance for changes in tracer concentrations between smoothing ! iterations that determines when to stop iterating [CU ~> conc] @@ -692,7 +706,9 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t PI_180 = atan(1.0)/45. - ans_2018 = .true.;if (present(answers_2018)) ans_2018 = answers_2018 + ans_date = 20181231 + if (present(answers_2018)) then ; if (.not.answers_2018) ans_date = 20190101 ; endif + if (present(answer_date)) ans_date = answer_date ! Open NetCDF file and if present, extract data and spatial coordinate information ! The convention adopted here requires that the data be written in (i,j,k) ordering. @@ -872,8 +888,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t good2(:,:) = good(:,:) fill2(:,:) = fill(:,:) - call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, dtr_iter_stop, & - answers_2018=answers_2018) + call fill_miss_2d(tr_outf, good2, fill2, tr_prev, G, dtr_iter_stop, answer_date=ans_date) ! if (debug) then ! call hchksum(tr_outf, 'field from fill_miss_2d ', G%HI, scale=I_scale) @@ -895,7 +910,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id, Time, conversion, G, t do j=js,je do i=is,ie tr_z(i,j,k) = data_in(i,j,k) * conversion - if (.not. ans_2018) mask_z(i,j,k) = 1. + if (ans_date >= 20190101) mask_z(i,j,k) = 1. if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0. enddo enddo