From e072bc7cce861b38df7ef1d2bde825d2d010f694 Mon Sep 17 00:00:00 2001 From: Jess <20195932+wrongkindofdoctor@users.noreply.github.com> Date: Fri, 6 Dec 2019 11:35:56 -0500 Subject: [PATCH] Merge in latest dev/gfdl updates (#40) * (*)Fixed dimensional inconsistency in P3M_functions Corrected dimensionally inconsistent expressions in P3M_functions.F90, notably in P3M_limiter and monotonize_cubic and a complete rewrite and simplification of is_cubic_monotonic. Also added comments documenting the units of all real variables in this module, and changed the code to use logical variables in place of integer "booleans", including in the return value from is_cubic_monotonic. These changes will change (fix) the answers when remapping variables with small numerical values, but no answers change in the MOM6-examples test cases. * +Added REMAPPING_2018 runtime option Added a new runtime option, REMAPPING_2018, which if set to false triggers the use of new, more accurate expressions in various parts of the ALE remapping code. By default, the older expressions are used, and all answers are bitwise identical, but there are new optional arguments to various routines related to remapping to trigger the use of new mathematically equivalent expressions. By default all answers are bitwise identical, but there are new and reordered entries in the MOM6_parameter_doc files. * Corrected the formatting of a doxygen comment * Added conversion factors to forcing diagnostics Added conversion factors to 4 mass-flux diagnostics and comments to 4 others on why no conversion factors are needed. All answers are bitwise identical. * Added correct scaling factors to chksum calls Added scale arguments to 5 chksum calls and grouped another two chksum calls while also adding the right scaling argument. All answers are bitwise identical. * +Unscales area before taking global sum Undoes the dimensional scaling of the cell areas before taking their global sum, so that the reproducing sum does not overflow when there is dimensional rescaling. All answers are bitwise identical when there is no rescaling, but this eliminates a source of inadvertent overflows or underflows in the global sums, and there is a new optional argument to compute_global_grid_integrals. * (*)Correct dimensionally inconsistent advective CFL Corrects the dimensionally inconsistent expressions for the CFL number in the tracer advection code, in which a negligible thickness had been added to the cell volume to avoid division by zero. This change does not alter the solutions in the MOM6-examples test cases, but now it permits dimensional rescaling of lengths over a much larger range, and it could change answers if the minimum layer thicknesses are small enough. * Unscale sea level before averaging Unscale interface heights before taking a global average via a reproducing sum in non-Boussinesq mode global diagnostics to permit dimensional consistency testing over a larger range. All answers are bitwise identical. * +Added an optional tmp_scale arg to global_i_mean Added an optional tmp_scale argument to global_i_mean and global_j_mean to specify an internal rescaling of variables being averaged before the reproducing sum. All answers are bitwise identical, but there are new optional arguments to two public interfaces. * Expand consistency testing with i-mean sponges Use tmp_scale when taking the i-mean interface heights for i-mean sponges, to give a greatly expanded range of dimensional consistency testing. All answers are bitwise identical. --- src/ALE/MOM_ALE.F90 | 43 ++- src/ALE/MOM_remapping.F90 | 37 ++- src/ALE/P3M_functions.F90 | 261 +++++++----------- src/ALE/regrid_edge_slopes.F90 | 142 ++++++---- src/ALE/regrid_edge_values.F90 | 229 +++++++++------ src/ALE/regrid_interp.F90 | 60 ++-- src/core/MOM.F90 | 2 +- src/core/MOM_dynamics_split_RK2.F90 | 6 +- src/core/MOM_forcing_type.F90 | 28 +- src/diagnostics/MOM_sum_output.F90 | 5 +- src/framework/MOM_spatial_means.F90 | 32 ++- .../MOM_fixed_initialization.F90 | 2 +- .../MOM_shared_initialization.F90 | 10 +- .../MOM_state_initialization.F90 | 7 +- src/parameterizations/lateral/MOM_MEKE.F90 | 2 +- .../vertical/MOM_set_diffusivity.F90 | 4 +- src/parameterizations/vertical/MOM_sponge.F90 | 2 +- src/tracer/MOM_tracer_advect.F90 | 18 +- 18 files changed, 494 insertions(+), 396 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index d7917f8cad..97232b22ca 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -63,9 +63,9 @@ module MOM_ALE !> ALE control structure type, public :: ALE_CS ; private - logical :: remap_uv_using_old_alg !< If true, uses the old "remapping via a delta z" - !! method. If False, uses the new method that - !! remaps between grids described by h. + logical :: remap_uv_using_old_alg !< If true, uses the old "remapping via a delta z" + !! method. If False, uses the new method that + !! remaps between grids described by h. real :: regrid_time_scale !< The time-scale used in blending between the current (old) grid !! and the target (new) grid [T ~> s] @@ -73,9 +73,13 @@ module MOM_ALE type(regridding_CS) :: regridCS !< Regridding parameters and work arrays type(remapping_CS) :: remapCS !< Remapping parameters and work arrays - integer :: nk !< Used only for queries, not directly by this module + integer :: nk !< Used only for queries, not directly by this module - logical :: remap_after_initialization !< Indicates whether to regrid/remap after initializing the state. + logical :: remap_after_initialization !< Indicates whether to regrid/remap after initializing the state. + + logical :: answers_2018 !< If true, use the order of arithmetic and expressions for remapping + !! that recover the answers from the end of 2018. Otherwise, use more + !! robust and accurate forms of mathematically equivalent expressions. logical :: show_call_tree !< For debugging @@ -145,6 +149,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) character(len=40) :: mdl = "MOM_ALE" ! This module's name. character(len=80) :: string ! Temporary strings real :: filter_shallow_depth, filter_deep_depth + logical :: default_2018_answers logical :: check_reconstruction logical :: check_remapping logical :: force_bounds_in_subcell @@ -192,11 +197,19 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) call get_param(param_file, mdl, "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, & "If true, values at the interfaces of boundary cells are "//& "extrapolated instead of piecewise constant", default=.false.) + call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & + "This sets the default value for the various _2018_ANSWERS parameters.", & + default=.true.) + call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", CS%answers_2018, & + "If true, use the order of arithmetic and expressions that recover the "//& + "answers from the end of 2018. Otherwise, use updated and more robust "//& + "forms of the same expressions.", default=default_2018_answers) call initialize_remapping( CS%remapCS, string, & boundary_extrapolation=remap_boundary_extrap, & check_reconstruction=check_reconstruction, & check_remapping=check_remapping, & - force_bounds_in_subcell=force_bounds_in_subcell) + force_bounds_in_subcell=force_bounds_in_subcell, & + answers_2018=CS%answers_2018) call get_param(param_file, mdl, "REMAP_AFTER_INITIALIZATION", CS%remap_after_initialization, & "If true, applies regridding and remapping immediately after "//& @@ -220,7 +233,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) "REGRID_FILTER_SHALLOW_DEPTH the filter weights adopt a cubic profile.", & units="m", default=0., scale=GV%m_to_H) call set_regrid_params(CS%regridCS, depth_of_time_filter_shallow=filter_shallow_depth, & - depth_of_time_filter_deep=filter_deep_depth) + depth_of_time_filter_deep=filter_deep_depth) call get_param(param_file, mdl, "REGRID_USE_OLD_DIRECTION", local_logical, & "If true, the regridding ntegrates upwards from the bottom for "//& "interface positions, much as the main model does. If false "//& @@ -1089,13 +1102,13 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext ! Local variables integer :: i, j, k - real :: hTmp(GV%ke) - real :: tmp(GV%ke) + real :: hTmp(GV%ke) ! A 1-d copy of h [H ~> m or kg m-2] + real :: tmp(GV%ke) ! A 1-d copy of a column of temperature [degC] or salinity [ppt] real, dimension(CS%nk,2) :: & - ppol_E !Edge value of polynomial + ppol_E ! Edge value of polynomial in [degC] or [ppt] real, dimension(CS%nk,3) :: & - ppol_coefs !Coefficients of polynomial - real :: h_neglect, h_neglect_edge + ppol_coefs ! Coefficients of polynomial, all in [degC] or [ppt] + real :: h_neglect, h_neglect_edge ! Tiny thicknesses [H ~> m or kg m-2] !### Try replacing both of these with GV%H_subroundoff if (GV%Boussinesq) then @@ -1116,7 +1129,8 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext ppol_E(:,:) = 0.0 ppol_coefs(:,:) = 0.0 !### Try to replace the following value of h_neglect with GV%H_subroundoff. - call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=h_neglect_edge ) + call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=h_neglect_edge, & + answers_2018=CS%answers_2018 ) call PPM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect ) if (bdry_extrap) & call PPM_boundary_extrapolation( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect ) @@ -1131,7 +1145,8 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext ppol_coefs(:,:) = 0.0 tmp(:) = tv%T(i,j,:) !### Try to replace the following value of h_neglect with GV%H_subroundoff. - call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=1.0e-10*GV%m_to_H ) + call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=1.0e-10*GV%m_to_H, & + answers_2018=CS%answers_2018 ) call PPM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect ) if (bdry_extrap) & call PPM_boundary_extrapolation(GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect ) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index f399aa2c0f..d7f8343993 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -33,6 +33,8 @@ module MOM_remapping logical :: check_remapping = .false. !> If true, the intermediate values used in remapping are forced to be bounded. logical :: force_bounds_in_subcell = .false. + !> If true use older, less acccurate expressions. + logical :: answers_2018 = .true. end type ! The following routines are visible to the outside world @@ -84,13 +86,14 @@ module MOM_remapping !> Set parameters within remapping object subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation, & - check_reconstruction, check_remapping, force_bounds_in_subcell) + check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018) type(remapping_CS), intent(inout) :: CS !< Remapping control structure character(len=*), optional, intent(in) :: remapping_scheme !< Remapping scheme to use logical, optional, intent(in) :: boundary_extrapolation !< Indicate to extrapolate in boundary cells logical, optional, intent(in) :: check_reconstruction !< Indicate to check reconstructions logical, optional, intent(in) :: check_remapping !< Indicate to check results of remapping logical, optional, intent(in) :: force_bounds_in_subcell !< Force subcells values to be bounded + logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions. if (present(remapping_scheme)) then call setReconstructionType( remapping_scheme, CS ) @@ -107,6 +110,9 @@ subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation, & if (present(force_bounds_in_subcell)) then CS%force_bounds_in_subcell = force_bounds_in_subcell endif + if (present(answers_2018)) then + CS%answers_2018 = answers_2018 + endif end subroutine remapping_set_param subroutine extract_member_remapping_CS(CS, remapping_scheme, degree, boundary_extrapolation, check_reconstruction, & @@ -392,22 +398,22 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & endif iMethod = INTEGRATION_PLM case ( REMAPPING_PPM_H4 ) - call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge ) + call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 ) call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) if ( CS%boundary_extrapolation ) then call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) endif iMethod = INTEGRATION_PPM case ( REMAPPING_PPM_IH4 ) - call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge ) + call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 ) call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) if ( CS%boundary_extrapolation ) then call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) endif iMethod = INTEGRATION_PPM case ( REMAPPING_PQM_IH4IH3 ) - call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge ) - call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_S, h_neglect ) + call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 ) + call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_S, h_neglect, answers_2018=CS%answers_2018 ) call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefs, h_neglect ) if ( CS%boundary_extrapolation ) then call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, & @@ -415,8 +421,8 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & endif iMethod = INTEGRATION_PQM case ( REMAPPING_PQM_IH6IH5 ) - call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E, h_neglect_edge ) - call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_S, h_neglect ) + call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 ) + call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_S, h_neglect, answers_2018=CS%answers_2018 ) call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefs, h_neglect ) if ( CS%boundary_extrapolation ) then call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, & @@ -1537,7 +1543,7 @@ end subroutine dzFromH1H2 !> Constructor for remapping control structure subroutine initialize_remapping( CS, remapping_scheme, boundary_extrapolation, & - check_reconstruction, check_remapping, force_bounds_in_subcell) + check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018) ! Arguments type(remapping_CS), intent(inout) :: CS !< Remapping control structure character(len=*), intent(in) :: remapping_scheme !< Remapping scheme to use @@ -1545,11 +1551,12 @@ subroutine initialize_remapping( CS, remapping_scheme, boundary_extrapolation, & logical, optional, intent(in) :: check_reconstruction !< Indicate to check reconstructions logical, optional, intent(in) :: check_remapping !< Indicate to check results of remapping logical, optional, intent(in) :: force_bounds_in_subcell !< Force subcells values to be bounded + logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions. - ! Note that remapping_scheme is mandatory fir initialize_remapping() + ! Note that remapping_scheme is mandatory for initialize_remapping() call remapping_set_param(CS, remapping_scheme=remapping_scheme, boundary_extrapolation=boundary_extrapolation, & check_reconstruction=check_reconstruction, check_remapping=check_remapping, & - force_bounds_in_subcell=force_bounds_in_subcell) + force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018) end subroutine initialize_remapping @@ -1615,6 +1622,7 @@ logical function remapping_unit_tests(verbose) data h2 /6*0.5/ ! 6 uniform layers with total depth of 3 type(remapping_CS) :: CS !< Remapping control structure real, allocatable, dimension(:,:) :: ppoly0_E, ppoly0_S, ppoly0_coefs + logical :: answers_2018 ! If true use older, less acccurate expressions. integer :: i real :: err, h_neglect, h_neglect_edge logical :: thisTest, v @@ -1622,6 +1630,7 @@ logical function remapping_unit_tests(verbose) v = verbose h_neglect = hNeglect_dflt h_neglect_edge = 1.0e-10 + answers_2018 = .true. write(*,*) '==== MOM_remapping: remapping_unit_tests =================' remapping_unit_tests = .false. ! Normally return false @@ -1643,7 +1652,7 @@ logical function remapping_unit_tests(verbose) remapping_unit_tests = remapping_unit_tests .or. thisTest thisTest = .false. - call initialize_remapping(CS, 'PPM_H4') + call initialize_remapping(CS, 'PPM_H4', answers_2018=answers_2018) if (verbose) write(*,*) 'h0 (test data)' if (verbose) call dumpGrid(n0,h0,x0,u0) @@ -1667,7 +1676,7 @@ logical function remapping_unit_tests(verbose) ppoly0_S(:,:) = 0.0 ppoly0_coefs(:,:) = 0.0 - call edge_values_explicit_h4( n0, h0, u0, ppoly0_E, h_neglect=1e-10 ) + call edge_values_explicit_h4( n0, h0, u0, ppoly0_E, h_neglect=1e-10, answers_2018=answers_2018 ) call PPM_reconstruction( n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect ) call PPM_boundary_extrapolation( n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect ) u1(:) = 0. @@ -1798,7 +1807,7 @@ logical function remapping_unit_tests(verbose) test_answer(v, 3, ppoly0_coefs(:,2), (/0.,4.,0./), 'Non-uniform line PLM: P1') call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,3.,5.,7.,9./), ppoly0_E, & - h_neglect=1e-10 ) + h_neglect=1e-10, answers_2018=answers_2018 ) ! The next two tests currently fail due to roundoff. thisTest = test_answer(v, 5, ppoly0_E(:,1), (/0.,2.,4.,6.,8./), 'Line H4: left edges') thisTest = test_answer(v, 5, ppoly0_E(:,2), (/2.,4.,6.,8.,10./), 'Line H4: right edges') @@ -1814,7 +1823,7 @@ logical function remapping_unit_tests(verbose) test_answer(v, 5, ppoly0_coefs(:,3), (/0.,0.,0.,0.,0./), 'Line PPM: P2') call edge_values_explicit_h4( 5, (/1.,1.,1.,1.,1./), (/1.,1.,7.,19.,37./), ppoly0_E, & - h_neglect=1e-10 ) + h_neglect=1e-10, answers_2018=answers_2018 ) ! The next two tests currently fail due to roundoff. thisTest = test_answer(v, 5, ppoly0_E(:,1), (/3.,0.,3.,12.,27./), 'Parabola H4: left edges') thisTest = test_answer(v, 5, ppoly0_E(:,2), (/0.,3.,12.,27.,48./), 'Parabola H4: right edges') diff --git a/src/ALE/P3M_functions.F90 b/src/ALE/P3M_functions.F90 index 1964cd25dd..da3fe5bb6b 100644 --- a/src/ALE/P3M_functions.F90 +++ b/src/ALE/P3M_functions.F90 @@ -25,20 +25,15 @@ 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, ppoly_E, ppoly_S, ppoly_coef, & - h_neglect ) +subroutine P3M_interpolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: h !< cell widths (size N) - real, dimension(:), intent(in) :: u !< cell averages (size N) - real, dimension(:,:), intent(inout) :: ppoly_E !< Edge value of polynomial, - !! with the same units as u. - real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial, - !! in the units of u over the units of h. - real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly - !! with the same units as u. + real, dimension(:), intent(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A] + real, dimension(:,:), intent(inout) :: ppoly_E !< Edge value of polynomial [A] + real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial [A H-1]. + real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A] real, optional, intent(in) :: h_neglect !< A negligibly small width for the - !! purpose of cell reconstructions - !! in the same units as h. + !! purpose of cell reconstructions [H] ! Call the limiter for p3m, which takes care of everything from ! computing the coefficients of the cubic to monotonizing it. @@ -64,28 +59,24 @@ end subroutine P3M_interpolation !! Step 3 of the monotonization process leaves all edge values unchanged. subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: h !< cell widths (size N) - real, dimension(:), intent(in) :: u !< cell averages (size N) - real, dimension(:,:), intent(inout) :: ppoly_E !< Edge value of polynomial, - !! with the same units as u. - real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial, - !! in the units of u over the units of h. - real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial + real, dimension(:), intent(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A] + real, dimension(:,:), intent(inout) :: ppoly_E !< Edge value of polynomial [A] + real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial [A H-1] + real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A] real, optional, intent(in) :: h_neglect !< A negligibly small width for - !! the purpose of cell reconstructions - !! in the same units as h. + !! the purpose of cell reconstructions [H] ! Local variables integer :: k ! loop index - integer :: monotonic ! boolean indicating whether the cubic is monotonic - real :: u0_l, u0_r ! edge values - real :: u1_l, u1_r ! edge slopes - real :: u_l, u_c, u_r ! left, center and right cell averages - real :: h_l, h_c, h_r ! left, center and right cell widths - real :: sigma_l, sigma_c, sigma_r ! left, center and right - ! van Leer slopes - real :: slope ! retained PLM slope + logical :: monotonic ! boolean indicating whether the cubic is monotonic + real :: u0_l, u0_r ! edge values [A] + real :: u1_l, u1_r ! edge slopes [A H-1] + real :: u_l, u_c, u_r ! left, center and right cell averages [A] + real :: h_l, h_c, h_r ! left, center and right cell widths [H] + real :: sigma_l, sigma_c, sigma_r ! left, center and right van Leer slopes [A H-1] + real :: slope ! retained PLM slope [A H-1] real :: eps - real :: hNeglect + real :: hNeglect ! A negligibly small thickness [H] hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect @@ -142,16 +133,9 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) slope = 0.0 endif - ! If the slopes are close to zero in machine precision and in absolute - ! value, we set the slope to zero. This prevents asymmetric representation - ! near extrema. These expressions are both nondimensional. - if ( abs(u1_l*h_c) < eps ) then - u1_l = 0.0 - endif - - if ( abs(u1_r*h_c) < eps ) then - u1_r = 0.0 - endif + ! If the slopes are small, set them to zero to prevent asymmetric representation near extrema. + if ( abs(u1_l*h_c) < epsilon(u_c)*abs(u_c) ) u1_l = 0.0 + if ( abs(u1_r*h_c) < epsilon(u_c)*abs(u_c) ) u1_r = 0.0 ! The edge slopes are limited from above by the respective ! one-sided slopes @@ -172,7 +156,7 @@ subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coef, h_neglect ) ! If cubic is not monotonic, monotonize it by modifiying the ! edge slopes, store the new edge slopes and recompute the ! cubic coefficients - if ( monotonic == 0 ) then + if ( .not.monotonic ) then call monotonize_cubic( h_c, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ) endif @@ -204,30 +188,25 @@ end subroutine P3M_limiter subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, & h_neglect, h_neglect_edge ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: h !< cell widths (size N) - real, dimension(:), intent(in) :: u !< cell averages (size N) - real, dimension(:,:), intent(inout) :: ppoly_E !< Edge value of polynomial, - !! with the same units as u. - real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial, - !! in the units of u over the units of h. - real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly - !! with the same units as u. + real, dimension(:), intent(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A] + real, dimension(:,:), intent(inout) :: ppoly_E !< Edge value of polynomial [A] + real, dimension(:,:), intent(inout) :: ppoly_S !< Edge slope of polynomial [A H-1] + real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A] real, optional, intent(in) :: h_neglect !< A negligibly small width for the - !! purpose of cell reconstructions - !! in the same units as h. + !! purpose of cell reconstructions [H] real, optional, intent(in) :: h_neglect_edge !< A negligibly small width - !! for the purpose of finding edge values - !! in the same units as h. + !! for the purpose of finding edge values [H] ! Local variables integer :: i0, i1 - integer :: monotonic - real :: u0, u1 - real :: h0, h1 - real :: b, c, d - real :: u0_l, u0_r - real :: u1_l, u1_r - real :: slope - real :: hNeglect, hNeglect_edge + logical :: monotonic ! boolean indicating whether the cubic is monotonic + real :: u0, u1 ! Values of u in two adjacent cells [A] + real :: h0, h1 ! Values of h in two adjacent cells, plus a smal increment [H] + real :: b, c, d ! Temporary variables [A] + real :: u0_l, u0_r ! Left and right edge values [A] + real :: u1_l, u1_r ! Left and right edge slopes [A H-1] + real :: slope ! The cell center slope [A H-1] + real :: hNeglect, hNeglect_edge ! Negligibly small thickness [H] hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect hNeglect_edge = hNeglect_edge_dflt @@ -281,7 +260,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, & call build_cubic_interpolant( h, i0, ppoly_E, ppoly_S, ppoly_coef ) monotonic = is_cubic_monotonic( ppoly_coef, i0 ) - if ( monotonic == 0 ) then + if ( .not.monotonic ) then call monotonize_cubic( h0, u0_l, u0_r, 0.0, slope, slope, u1_l, u1_r ) ! Rebuild cubic after monotonization @@ -340,7 +319,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coef, & call build_cubic_interpolant( h, i1, ppoly_E, ppoly_S, ppoly_coef ) monotonic = is_cubic_monotonic( ppoly_coef, i1 ) - if ( monotonic == 0 ) then + if ( .not.monotonic ) then call monotonize_cubic( h1, u0_l, u0_r, slope, 0.0, slope, u1_l, u1_r ) ! Rebuild cubic after monotonization @@ -360,19 +339,17 @@ end subroutine P3M_boundary_extrapolation !! NOTE: edge values and slopes MUST have been properly calculated prior to !! calling this routine. subroutine build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coef ) - real, dimension(:), intent(in) :: h !< cell widths (size N) + real, dimension(:), intent(in) :: h !< cell widths (size N) [H] integer, intent(in) :: k !< The index of the cell to work on - real, dimension(:,:), intent(in) :: ppoly_E !< Edge value of polynomial, - !! with the same units as u. - real, dimension(:,:), intent(in) :: ppoly_S !< Edge slope of polynomial, - !! in the units of u over the units of h. - real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly - !! with the same units as u. + real, dimension(:,:), intent(in) :: ppoly_E !< Edge value of polynomial in arbitrary units [A] + real, dimension(:,:), intent(in) :: ppoly_S !< Edge slope of polynomial [A H-1] + real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial [A] + ! Local variables - real :: u0_l, u0_r ! edge values - real :: u1_l, u1_r ! edge slopes - real :: h_c ! cell width - real :: a0, a1, a2, a3 ! cubic coefficients + real :: u0_l, u0_r ! edge values [A] + real :: u1_l, u1_r ! edge slopes times the cell width [A] + real :: h_c ! cell width [H] + real :: a0, a1, a2, a3 ! cubic coefficients [A] h_c = h(k) @@ -400,63 +377,30 @@ end subroutine build_cubic_interpolant !! This function checks whether the cubic curve in cell k is monotonic. !! If so, returns 1. Otherwise, returns 0. !! -!! The cubic is monotonic if the first derivative is single-signed in [0,1]. +!! The cubic is monotonic if the first derivative is single-signed in (0,1). !! Hence, we check whether the roots (if any) lie inside this interval. If there !! is no root or if both roots lie outside this interval, the cubic is monotonic. -integer function is_cubic_monotonic( ppoly_coef, k ) - real, dimension(:,:), intent(in) :: ppoly_coef !< Coefficients of cubic polynomial +logical function is_cubic_monotonic( ppoly_coef, k ) + real, dimension(:,:), intent(in) :: ppoly_coef !< Coefficients of cubic polynomial in arbitary units [A] integer, intent(in) :: k !< The index of the cell to work on ! Local variables - integer :: monotonic ! boolean indicating if monotonic or not - real :: a0, a1, a2, a3 ! cubic coefficients - real :: a, b, c ! coefficients of first derivative - real :: xi_0, xi_1 ! roots of first derivative (if any !) - real :: rho - real :: eps - - ! Define the radius of the ball around 0 and 1 in which all values are assumed - ! to be equal to 0 or 1, respectively - eps = 1e-14 - - a0 = ppoly_coef(k,1) - a1 = ppoly_coef(k,2) - a2 = ppoly_coef(k,3) - a3 = ppoly_coef(k,4) - - a = a1 - b = 2.0 * a2 - c = 3.0 * a3 - - xi_0 = -1.0 - xi_1 = -1.0 - - rho = b*b - 4.0*a*c - - if ( rho >= 0.0 ) then - if ( abs(c) > 1e-15 ) then - xi_0 = 0.5 * ( -b - sqrt( rho ) ) / c - xi_1 = 0.5 * ( -b + sqrt( rho ) ) / c - elseif ( abs(b) > 1e-15 ) then - xi_0 = - a / b - xi_1 = - a / b - endif - - ! If one of the roots of the first derivative lies in (0,1), - ! the cubic is not monotonic. - if ( ( (xi_0 > eps) .AND. (xi_0 < 1.0-eps) ) .OR. & - ( (xi_1 > eps) .AND. (xi_1 < 1.0-eps) ) ) then - monotonic = 0 - else - monotonic = 1 - endif - - else ! there are no real roots --> cubic is monotonic - monotonic = 1 + real :: a, b, c ! Coefficients of the first derivative of the cubic [A] + + a = ppoly_coef(k,2) + b = 2.0 * ppoly_coef(k,3) + c = 3.0 * ppoly_coef(k,4) + + ! Look for real roots of the quadratic derivative equation, c*x**2 + b*x + a = 0, in (0, 1) + if (b*b - 4.0*a*c <= 0.0) then ! The cubic is monotonic everywhere. + is_cubic_monotonic = .true. + elseif (a * (a + (b + c)) < 0.0) then ! The derivative changes sign between the endpoints of (0, 1) + is_cubic_monotonic = .false. + elseif (b * (b + 2.0*c) < 0.0) then ! The second derivative changes sign inside of (0, 1) + is_cubic_monotonic = .false. + else + is_cubic_monotonic = .true. endif - ! Set the return value - is_cubic_monotonic = monotonic - end function is_cubic_monotonic !> Monotonize a cubic curve by modifying the edge slopes. @@ -487,30 +431,27 @@ end function is_cubic_monotonic !! edge or onto the right edge. subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ) - real, intent(in) :: h !< cell width - real, intent(in) :: u0_l !< left edge value - real, intent(in) :: u0_r !< right edge value - real, intent(in) :: sigma_l !< left 2nd-order slopes - real, intent(in) :: sigma_r !< right 2nd-order slopes - real, intent(in) :: slope !< limited PLM slope - real, intent(inout) :: u1_l !< left edge slopes - real, intent(inout) :: u1_r !< right edge slopes + real, intent(in) :: h !< cell width [H] + real, intent(in) :: u0_l !< left edge value in arbitrary units [A] + real, intent(in) :: u0_r !< right edge value [A] + real, intent(in) :: sigma_l !< left 2nd-order slopes [A H-1] + real, intent(in) :: sigma_r !< right 2nd-order slopes [A H-1] + real, intent(in) :: slope !< limited PLM slope [A H-1] + real, intent(inout) :: u1_l !< left edge slopes [A H-1] + real, intent(inout) :: u1_r !< right edge slopes [A H-1] ! Local variables - integer :: found_ip - integer :: inflexion_l ! bool telling if inflex. pt must be on left - integer :: inflexion_r ! bool telling if inflex. pt must be on right - real :: eps - real :: a1, a2, a3 - real :: u1_l_tmp ! trial left edge slope - real :: u1_r_tmp ! trial right edge slope - real :: xi_ip ! location of inflexion point - real :: slope_ip ! slope at inflexion point - - eps = 1e-14 - - found_ip = 0 - inflexion_l = 0 - inflexion_r = 0 + logical :: found_ip + logical :: inflexion_l ! bool telling if inflex. pt must be on left + logical :: inflexion_r ! bool telling if inflex. pt must be on right + real :: a1, a2, a3 ! Temporary slopes times the cell width [A] + real :: u1_l_tmp ! trial left edge slope [A H-1] + real :: u1_r_tmp ! trial right edge slope [A H-1] + real :: xi_ip ! location of inflexion point in cell coordinates (0,1) [nondim] + real :: slope_ip ! slope at inflexion point times cell width [A] + + found_ip = .false. + inflexion_l = .false. + inflexion_r = .false. ! If the edge slopes are inconsistent w.r.t. the limited PLM slope, ! set them to zero @@ -537,7 +478,7 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ! If the inflexion point lies in [0,1], change boolean value if ( (xi_ip >= 0.0) .AND. (xi_ip <= 1.0) ) then - found_ip = 1 + found_ip = .true. endif endif @@ -546,25 +487,25 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r ! decide on which side we want to collapse the inflexion point. ! If the inflexion point lies on one of the edges, the cubic is ! guaranteed to be monotonic - if ( found_ip == 1 ) then + if ( found_ip ) then slope_ip = a1 + 2.0*a2*xi_ip + 3.0*a3*xi_ip*xi_ip ! Check whether slope is consistent if ( slope_ip*slope < 0.0 ) then if ( abs(sigma_l) < abs(sigma_r) ) then - inflexion_l = 1 + inflexion_l = .true. else - inflexion_r = 1 + inflexion_r = .true. endif endif endif ! found_ip ! At this point, if the cubic is not monotonic, we know where the ! inflexion point should lie. When the cubic is monotonic, both - ! 'inflexion_l' and 'inflexion_r' are set to 0 and nothing is to be done. + ! 'inflexion_l' and 'inflexion_r' are false and nothing is to be done. ! Move inflexion point on the left - if ( inflexion_l == 1 ) then + if ( inflexion_l ) then u1_l_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_r u1_r_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_l @@ -594,7 +535,7 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r endif ! end treating case with inflexion point on the left ! Move inflexion point on the right - if ( inflexion_r == 1 ) then + if ( inflexion_r ) then u1_l_tmp = 3.0*(u0_r-u0_l)/h - 2.0*u1_r u1_r_tmp = 1.5*(u0_r-u0_l)/h - 0.5*u1_l @@ -623,13 +564,9 @@ subroutine monotonize_cubic( h, u0_l, u0_r, sigma_l, sigma_r, slope, u1_l, u1_r endif ! end treating case with inflexion point on the right - if ( abs(u1_l*h) < eps ) then - u1_l = 0.0 - endif - - if ( abs(u1_r*h) < eps ) then - u1_r = 0.0 - endif + ! Zero out negligibly small slopes. + if ( abs(u1_l*h) < epsilon(u0_l) * (abs(u0_l) + abs(u0_r)) ) u1_l = 0.0 + if ( abs(u1_r*h) < epsilon(u0_l) * (abs(u0_l) + abs(u0_r)) ) u1_r = 0.0 end subroutine monotonize_cubic diff --git a/src/ALE/regrid_edge_slopes.F90 b/src/ALE/regrid_edge_slopes.F90 index c22a524683..8d5c055907 100644 --- a/src/ALE/regrid_edge_slopes.F90 +++ b/src/ALE/regrid_edge_slopes.F90 @@ -46,35 +46,39 @@ module regrid_edge_slopes !! !! 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 ) +subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answers_2018 ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: h !< cell widths (size N) - real, dimension(:), intent(in) :: u !< cell average properties (size N) - real, dimension(:,:), intent(inout) :: edge_slopes !< Returned edge slopes, with the - !! same units as u divided by the units of h. + real, dimension(:), intent(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A] + real, dimension(:,:), intent(inout) :: edge_slopes !< Returned edge slopes [A H-1] real, optional, intent(in) :: h_neglect !< A negligibly small width + logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions. ! Local variables integer :: i, j ! loop indexes - real :: h0, h1 ! cell widths - real :: h0_2, h1_2, h0h1 - real :: h0_3, h1_3 - real :: d - real :: alpha, beta ! stencil coefficients - real :: a, b - real, dimension(5) :: x ! system used to enforce - real, dimension(4,4) :: Asys ! boundary conditions + real :: h0, h1 ! cell widths [H] + real :: h0_2, h1_2, h0h1 ! products of cell widths [H2] + real :: h0_3, h1_3 ! products of three cell widths [H3] + real :: d ! A demporary variable [H3] + real :: alpha, beta ! stencil coefficients [nondim] + real :: a, b ! weights of cells [H-1] + real, parameter :: C1_12 = 1.0 / 12.0 + real, dimension(5) :: x ! Coordinate system with 0 at edges [H] + real :: dx, xavg ! Differences and averages of successive values of x [H] + real, dimension(4,4) :: Asys ! matrix used to find boundary conditions real, dimension(4) :: Bsys, Csys real, dimension(3) :: Dsys - real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal) - tri_d, & ! trid. system (middle diagonal) - tri_u, & ! trid. system (upper diagonal) - tri_b, & ! trid. system (unknowns vector) - tri_x ! trid. system (rhs) - real :: hNeglect ! A negligible thicness in the same units as h. - real :: hNeglect3 ! hNeglect^3 in the same units as h^3. + real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal) [nondim] + tri_d, & ! trid. system (middle diagonal) [nondim] + tri_u, & ! trid. system (upper diagonal) [nondim] + tri_b, & ! trid. system (unknowns vector) [A H-1] + tri_x ! trid. system (rhs) [A H-1] + real :: hNeglect ! A negligible thickness [H]. + real :: hNeglect3 ! hNeglect^3 [H3]. + logical :: use_2018_answers ! If true use older, less acccurate expressions. 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 ! Loop on cells (except last one) do i = 1,N-1 @@ -113,12 +117,18 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect ) enddo do i = 1,4 - - do j = 1,4 - Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - enddo - - Bsys(i) = u(i) * ( h(i) ) + dx = h(i) + if (use_2018_answers) then + do j = 1,4 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo + else ! Use expressions with less sensitivity to roundoff + xavg = 0.5 * (x(i+1) + x(i)) + Asys(i,1) = dx + Asys(i,2) = dx * xavg + Asys(i,3) = dx * (xavg**2 + C1_12*dx**2) + Asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2) + endif + + Bsys(i) = u(i) * dx enddo @@ -139,12 +149,17 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect ) enddo do i = 1,4 - - do j = 1,4 - Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - enddo - - Bsys(i) = u(N-4+i) * ( h(N-4+i) ) + dx = h(N-4+i) + if (use_2018_answers) then + do j = 1,4 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo + else ! Use expressions with less sensitivity to roundoff + xavg = 0.5 * (x(i+1) + x(i)) + Asys(i,1) = dx + Asys(i,2) = dx * xavg + Asys(i,3) = dx * (xavg**2 + C1_12*dx**2) + Asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2) + endif + Bsys(i) = u(N-4+i) * dx enddo @@ -173,14 +188,13 @@ end subroutine edge_slopes_implicit_h3 !------------------------------------------------------------------------------ !> Compute ih5 edge values (implicit fifth order accurate) -subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect ) +subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answers_2018 ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: h !< cell widths (size N) - real, dimension(:), intent(in) :: u !< cell average properties (size N) - real, dimension(:,:), intent(inout) :: edge_slopes !< Returned edge slopes, with the - !! same units as u divided by the units of h. - real, optional, intent(in) :: h_neglect !< A negligibly small width - !! in the same units as h. + real, dimension(:), intent(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A] + real, dimension(:,:), intent(inout) :: edge_slopes !< Returned edge slopes [A H-1] + real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions. ! ----------------------------------------------------------------------------- ! Fifth-order implicit estimates of edge values are based on a four-cell, ! three-edge stencil. A tridiagonal system is set up and is based on @@ -232,8 +246,11 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect ) real :: h2ph3_3, h2ph3_4 ! ... real :: alpha, beta ! stencil coefficients real :: a, b, c, d ! " - real, dimension(7) :: x ! system used to enforce - real, dimension(6,6) :: Asys ! boundary conditions + real, dimension(7) :: x ! Coordinate system with 0 at edges [same units as h] + real, parameter :: C1_12 = 1.0 / 12.0 + real, parameter :: C5_6 = 5.0 / 6.0 + real :: dx, xavg ! Differences and averages of successive values of x [same units as h] + real, dimension(6,6) :: Asys ! matrix used to find boundary conditions real, dimension(6) :: Bsys, Csys ! ... real, dimension(5) :: Dsys ! derivative real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal) @@ -241,9 +258,11 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect ) tri_u, & ! trid. system (upper diagonal) tri_b, & ! trid. system (unknowns vector) tri_x ! trid. system (rhs) - real :: hNeglect ! A negligible thicness in the same units as h. + real :: hNeglect ! A negligible thickness in the same units as h. + logical :: use_2018_answers ! If true use older, less acccurate expressions. hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect + use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018 ! Loop on cells (except last one) do k = 2,N-2 @@ -473,11 +492,20 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect ) do i = 1,6 - do j = 1,6 - Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - enddo - - Bsys(i) = u(i) * h(i) + dx = h(i) + if (use_2018_answers) then + do j = 1,6 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo + else ! Use expressions with less sensitivity to roundoff + xavg = 0.5 * (x(i+1) + x(i)) + Asys(i,1) = dx + Asys(i,2) = dx * xavg + Asys(i,3) = dx * (xavg**2 + C1_12*dx**2) + Asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2) + Asys(i,5) = dx * (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4) + Asys(i,6) = dx * xavg * (xavg**4 + C5_6*xavg**2*dx**2 + 0.0625*dx**4) + endif + + Bsys(i) = u(i) * dx enddo @@ -612,13 +640,19 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect ) enddo do i = 1,6 - - do j = 1,6 - Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - enddo - - Bsys(i) = u(N-6+i) * h(N-6+i) - + dx = h(N-6+i) + if (use_2018_answers) then + do j = 1,6 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo + else ! Use expressions with less sensitivity to roundoff + xavg = 0.5 * (x(i+1) + x(i)) + Asys(i,1) = dx + Asys(i,2) = dx * xavg + Asys(i,3) = dx * (xavg**2 + C1_12*dx**2) + Asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2) + Asys(i,5) = dx * (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4) + Asys(i,6) = dx * xavg * (xavg**4 + C5_6*xavg**2*dx**2 + 0.0625*dx**4) + endif + Bsys(i) = u(N-6+i) * dx enddo call solve_linear_system( Asys, Bsys, Csys, 6 ) diff --git a/src/ALE/regrid_edge_values.F90 b/src/ALE/regrid_edge_values.F90 index d27d69153c..f82e42e0e6 100644 --- a/src/ALE/regrid_edge_values.F90 +++ b/src/ALE/regrid_edge_values.F90 @@ -46,23 +46,20 @@ module regrid_edge_values !! Therefore, boundary cells are treated as if they were local extrama. subroutine bound_edge_values( N, h, u, edge_val, h_neglect ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: h !< cell widths (size N) - real, dimension(:), intent(in) :: u !< cell average properties (size N) - real, dimension(:,:), intent(inout) :: edge_val !< Potentially modified edge values, - !! with the same units as u. - real, optional, intent(in) :: h_neglect !< A negligibly small width - !! in the same units as h. + real, dimension(:), intent(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A] + real, dimension(:,:), intent(inout) :: edge_val !< Potentially modified edge values [A] + real, optional, intent(in) :: h_neglect !< A negligibly small width [H] ! Local variables integer :: k ! loop index integer :: k0, k1, k2 - real :: h_l, h_c, h_r - real :: u_l, u_c, u_r - real :: u0_l, u0_r + real :: h_l, h_c, h_r ! Layer thicknesses [H] + real :: u_l, u_c, u_r ! Cell average properties [A] + real :: u0_l, u0_r ! Edge values of properties [A] real :: sigma_l, sigma_c, sigma_r ! left, center and right - ! van Leer slopes - real :: slope ! retained PLM slope - - real :: hNeglect ! A negligible thicness in the same units as h. + ! van Leer slopes [A H-1] + real :: slope ! retained PLM slope [A H-1] + real :: hNeglect ! A negligible thickness [H]. hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect @@ -175,15 +172,15 @@ end subroutine average_discontinuous_edge_values !! If so and if they are not monotonic, replace each edge value by their average. subroutine check_discontinuous_edge_values( N, u, edge_val ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: u !< cell averages (size N) - real, dimension(:,:), intent(inout) :: edge_val !< Cell edge values with the same units as u. + real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A] + real, dimension(:,:), intent(inout) :: edge_val !< Cell edge values [A]. ! Local variables integer :: k ! loop index - real :: u0_minus ! left value at given edge - real :: u0_plus ! right value at given edge - real :: um_minus ! left cell average - real :: um_plus ! right cell average - real :: u0_avg ! avg value at given edge + real :: u0_minus ! left value at given edge [A] + real :: u0_plus ! right value at given edge [A] + real :: um_minus ! left cell average [A] + real :: um_plus ! right cell average [A] + real :: u0_avg ! avg value at given edge [A] ! Loop on interior cells do k = 1,N-1 @@ -227,16 +224,15 @@ end subroutine check_discontinuous_edge_values !! Boundary edge values are set to be equal to the boundary cell averages. subroutine edge_values_explicit_h2( N, h, u, edge_val, h_neglect ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: h !< cell widths (size N) - real, dimension(:), intent(in) :: u !< cell average properties (size N) - real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values, with the - !! same units as u; the second index size is 2. - real, optional, intent(in) :: h_neglect !< A negligibly small width + real, dimension(:), intent(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A] + real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values [A]; the second index size is 2. + real, optional, intent(in) :: h_neglect !< A negligibly small width [H] ! Local variables integer :: k ! loop index - real :: h0, h1 ! cell widths - real :: u0, u1 ! cell averages - real :: hNeglect ! A negligible thicness in the same units as h. + real :: h0, h1 ! cell widths [H] + real :: u0, u1 ! cell averages [A] + real :: hNeglect ! A negligible thickness [H] hNeglect = hNeglect_edge_dflt ; if (present(h_neglect)) hNeglect = h_neglect @@ -289,24 +285,29 @@ 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 ) +subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answers_2018 ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: h !< cell widths (size N) - real, dimension(:), intent(in) :: u !< cell average properties (size N) - real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values, with the - !! same units as u; the second index size is 2. - real, optional, intent(in) :: h_neglect !< A negligibly small width + real, dimension(:), intent(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A] + real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values [A]; the second index size is 2. + real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions. + ! Local variables integer :: i, j - real :: u0, u1, u2, u3 - real :: h0, h1, h2, h3 - real :: f1, f2, f3 ! auxiliary variables + real :: u0, u1, u2, u3 ! temporary properties [A] + real :: h0, h1, h2, h3 ! temporary thicknesses [H] + real :: f1, f2, f3 ! auxiliary variables with various units real :: e ! edge value - real, dimension(5) :: x ! used to compute edge + real, dimension(5) :: x ! Coordinate system with 0 at edges [H] + real, parameter :: C1_12 = 1.0 / 12.0 + real :: dx, xavg ! Differences and averages of successive values of x [same units as h] real, dimension(4,4) :: A ! values near the boundaries real, dimension(4) :: B, C - real :: hNeglect ! A negligible thicness in the same units as h. + real :: hNeglect ! A negligible thickness in the same units as h. + logical :: use_2018_answers ! If true use older, less acccurate expressions. + use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018 hNeglect = hNeglect_edge_dflt ; if (present(h_neglect)) hNeglect = h_neglect ! Loop on interior cells @@ -372,12 +373,18 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect ) enddo do i = 1,4 + dx = max(f1, h(i) ) + if (use_2018_answers) then + do j = 1,4 ; A(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) ; enddo + else ! Use expressions with less sensitivity to roundoff + xavg = 0.5 * (x(i+1) + x(i)) + A(i,1) = dx + A(i,2) = dx * xavg + A(i,3) = dx * (xavg**2 + C1_12*dx**2) + A(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2) + endif - do j = 1,4 - A(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) - enddo - - B(i) = u(i) * max(f1, h(i) ) + B(i) = u(i) * dx enddo @@ -410,12 +417,18 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect ) enddo do i = 1,4 + dx = max(f1, h(N-4+i) ) + if (use_2018_answers) then + do j = 1,4 ; A(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) ; enddo + else ! Use expressions with less sensitivity to roundoff + xavg = 0.5 * (x(i+1) + x(i)) + A(i,1) = dx + A(i,2) = dx * xavg + A(i,3) = dx * (xavg**2 + C1_12*dx**2) + A(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2) + endif - do j = 1,4 - A(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / real(j) - enddo - - B(i) = u(N-4+i) * max(f1, h(N-4+i) ) + B(i) = u(N-4+i) * dx enddo @@ -475,21 +488,24 @@ 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 ) +subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answers_2018 ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: h !< cell widths (size N) - real, dimension(:), intent(in) :: u !< cell average properties (size N) - real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values, with the - !! same units as u; the second index size is 2. - real, optional, intent(in) :: h_neglect !< A negligibly small width + real, dimension(:), intent(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A] + real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values [A]; the second index size is 2. + real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions. + ! Local variables integer :: i, j ! loop indexes - real :: h0, h1 ! cell widths + real :: h0, h1 ! cell widths [H] real :: h0_2, h1_2, h0h1 real :: d2, d4 real :: alpha, beta ! stencil coefficients real :: a, b - real, dimension(5) :: x ! system used to enforce + real, dimension(5) :: x ! Coordinate system with 0 at edges [H] + real, parameter :: C1_12 = 1.0 / 12.0 + real :: dx, xavg ! Differences and averages of successive values of x [H] real, dimension(4,4) :: Asys ! boundary conditions real, dimension(4) :: Bsys, Csys real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal) @@ -497,8 +513,10 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect ) tri_u, & ! trid. system (upper diagonal) tri_b, & ! trid. system (unknowns vector) tri_x ! trid. system (rhs) - real :: hNeglect ! A negligible thicness in the same units as h. + real :: hNeglect ! A negligible thickness [H] + logical :: use_2018_answers ! If true use older, less acccurate expressions. + use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018 hNeglect = hNeglect_edge_dflt ; if (present(h_neglect)) hNeglect = h_neglect ! Loop on cells (except last one) @@ -543,12 +561,18 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect ) enddo do i = 1,4 + dx = max(h0, h(i) ) + if (use_2018_answers) then + do j = 1,4 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo + else ! Use expressions with less sensitivity to roundoff + xavg = 0.5 * (x(i+1) + x(i)) + Asys(i,1) = dx + Asys(i,2) = dx * xavg + Asys(i,3) = dx * (xavg**2 + C1_12*dx**2) + Asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2) + endif - do j = 1,4 - Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - enddo - - Bsys(i) = u(i) * max( h0, h(i) ) + Bsys(i) = u(i) * dx enddo @@ -566,12 +590,17 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect ) enddo do i = 1,4 - - do j = 1,4 - Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - enddo - - Bsys(i) = u(N-4+i) * max( h0, h(N-4+i) ) + dx = max(h0, h(N-4+i) ) + if (use_2018_answers) then + do j = 1,4 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo + else ! Use expressions with less sensitivity to roundoff + xavg = 0.5 * (x(i+1) + x(i)) + Asys(i,1) = dx + Asys(i,2) = dx * xavg + Asys(i,3) = dx * (xavg**2 + C1_12*dx**2) + Asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2) + endif + Bsys(i) = u(N-4+i) * dx enddo @@ -628,16 +657,17 @@ end subroutine edge_values_implicit_h4 !! 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 ) +subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answers_2018 ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: h !< cell widths (size N) - real, dimension(:), intent(in) :: u !< cell average properties (size N) - real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values, with the - !! same units as u; the second index size is 2. - real, optional, intent(in) :: h_neglect !< A negligibly small width + real, dimension(:), intent(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell average properties (size N) in arbitrary units [A] + real, dimension(:,:), intent(inout) :: edge_val !< Returned edge values [A]; the second index size is 2. + real, optional, intent(in) :: h_neglect !< A negligibly small width [H] + logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions. + ! Local variables integer :: i, j, k ! loop indexes - real :: h0, h1, h2, h3 ! cell widths + real :: h0, h1, h2, h3 ! cell widths [H] real :: g, g_2, g_3 ! the following are real :: g_4, g_5, g_6 ! auxiliary variables real :: d2, d3, d4, d5, d6 ! to set up the systems @@ -654,7 +684,10 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect ) real :: h0ph1_5, h2ph3_5 ! ... real :: alpha, beta ! stencil coefficients real :: a, b, c, d ! " - real, dimension(7) :: x ! system used to enforce + real, dimension(7) :: x ! Coordinate system with 0 at edges [same units as h] + real, parameter :: C1_12 = 1.0 / 12.0 + real, parameter :: C5_6 = 5.0 / 6.0 + real :: dx, xavg ! Differences and averages of successive values of x [same units as h] real, dimension(6,6) :: Asys ! boundary conditions real, dimension(6) :: Bsys, Csys ! ... real, dimension(N+1) :: tri_l, & ! trid. system (lower diagonal) @@ -662,8 +695,10 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect ) tri_u, & ! trid. system (upper diagonal) tri_b, & ! trid. system (unknowns vector) tri_x ! trid. system (rhs) - real :: hNeglect ! A negligible thicness in the same units as h. + real :: hNeglect ! A negligible thickness [H]. + logical :: use_2018_answers ! If true use older, less acccurate expressions. + use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018 hNeglect = hNeglect_edge_dflt ; if (present(h_neglect)) hNeglect = h_neglect ! Loop on cells (except last one) @@ -913,12 +948,19 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect ) enddo do i = 1,6 - - do j = 1,6 - Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - enddo - - Bsys(i) = u(i) * max( g, h(i) ) + dx = max( g, h(i) ) + if (use_2018_answers) then + do j = 1,6 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo + else ! Use expressions with less sensitivity to roundoff + xavg = 0.5 * (x(i+1) + x(i)) + Asys(i,1) = dx + Asys(i,2) = dx * xavg + Asys(i,3) = dx * (xavg**2 + C1_12*dx**2) + Asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2) + Asys(i,5) = dx * (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4) + Asys(i,6) = dx * xavg * (xavg**4 + C5_6*xavg**2*dx**2 + 0.0625*dx**4) + endif + Bsys(i) = u(i) * dx enddo @@ -1058,12 +1100,19 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect ) enddo do i = 1,6 - - do j = 1,6 - Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j - enddo - - Bsys(i) = u(N-6+i) * max( g, h(N-6+i) ) + dx = max( g, h(N-6+i) ) + if (use_2018_answers) then + do j = 1,6 ; Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j ; enddo + else ! Use expressions with less sensitivity to roundoff + xavg = 0.5 * (x(i+1) + x(i)) + Asys(i,1) = dx + Asys(i,2) = dx * xavg + Asys(i,3) = dx * (xavg**2 + C1_12*dx**2) + Asys(i,4) = dx * xavg * (xavg**2 + 0.25*dx**2) + Asys(i,5) = dx * (xavg**4 + 0.5*xavg**2*dx**2 + 0.0125*dx**4) + Asys(i,6) = dx * xavg * (xavg**4 + C5_6*xavg**2*dx**2 + 0.0625*dx**4) + endif + Bsys(i) = u(N-6+i) * dx enddo diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index d2c384c15e..ace311cc21 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -30,6 +30,9 @@ module regrid_interp !> Indicate whether high-order boundary extrapolation should be used within !! boundary cells logical :: boundary_extrapolation + + !> If true use older, less acccurate expressions. + logical :: answers_2018 = .true. end type interp_CS_type public regridding_set_ppolys, interpolate_grid @@ -112,7 +115,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_P1M_H4 ) degree = DEGREE_1 if ( n0 >= 4 ) then - call edge_values_explicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge ) + call edge_values_explicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answers_2018=CS%answers_2018 ) else call edge_values_explicit_h2( n0, h0, densities, ppoly0_E, h_neglect_edge ) endif @@ -124,7 +127,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_P1M_IH4 ) degree = DEGREE_1 if ( n0 >= 4 ) then - call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge ) + call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answers_2018=CS%answers_2018 ) else call edge_values_explicit_h2( n0, h0, densities, ppoly0_E, h_neglect_edge ) endif @@ -143,7 +146,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_PPM_H4 ) if ( n0 >= 4 ) then degree = DEGREE_2 - call edge_values_explicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge ) + call edge_values_explicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answers_2018=CS%answers_2018 ) call PPM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, & @@ -161,7 +164,7 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_PPM_IH4 ) if ( n0 >= 4 ) then degree = DEGREE_2 - call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge ) + call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answers_2018=CS%answers_2018 ) call PPM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefs, h_neglect ) if (extrapolate) then call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, & @@ -179,8 +182,8 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_P3M_IH4IH3 ) if ( n0 >= 4 ) then degree = DEGREE_3 - call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge ) - call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_S, h_neglect ) + call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answers_2018=CS%answers_2018 ) + call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_S, h_neglect, answers_2018=CS%answers_2018 ) call P3M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect ) if (extrapolate) then @@ -199,8 +202,8 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_P3M_IH6IH5 ) if ( n0 >= 6 ) then degree = DEGREE_3 - call edge_values_implicit_h6( n0, h0, densities, ppoly0_E, h_neglect_edge ) - call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_S, h_neglect ) + call edge_values_implicit_h6( n0, h0, densities, ppoly0_E, h_neglect_edge, answers_2018=CS%answers_2018 ) + call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_S, h_neglect, answers_2018=CS%answers_2018 ) call P3M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect ) if (extrapolate) then @@ -219,8 +222,8 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_PQM_IH4IH3 ) if ( n0 >= 4 ) then degree = DEGREE_4 - call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge ) - call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_S, h_neglect ) + call edge_values_implicit_h4( n0, h0, densities, ppoly0_E, h_neglect_edge, answers_2018=CS%answers_2018 ) + call edge_slopes_implicit_h3( n0, h0, densities, ppoly0_S, h_neglect, answers_2018=CS%answers_2018 ) call PQM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect ) if (extrapolate) then @@ -239,8 +242,8 @@ subroutine regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, & case ( INTERPOLATION_PQM_IH6IH5 ) if ( n0 >= 6 ) then degree = DEGREE_4 - call edge_values_implicit_h6( n0, h0, densities, ppoly0_E, h_neglect_edge ) - call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_S, h_neglect ) + call edge_values_implicit_h6( n0, h0, densities, ppoly0_E, h_neglect_edge, answers_2018=CS%answers_2018 ) + call edge_slopes_implicit_h5( n0, h0, densities, ppoly0_S, h_neglect, answers_2018=CS%answers_2018 ) call PQM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_S, & ppoly0_coefs, h_neglect ) if (extrapolate) then @@ -264,7 +267,7 @@ end subroutine regridding_set_ppolys !! 'ppoly0' (possibly discontinuous), the coordinates of the new grid 'grid1' !! are determined by finding the corresponding target interface densities. subroutine interpolate_grid( n0, h0, x0, ppoly0_E, ppoly0_coefs, & - target_values, degree, n1, h1, x1 ) + target_values, degree, n1, h1, x1, answers_2018 ) integer, intent(in) :: n0 !< Number of points on source grid real, dimension(:), intent(in) :: h0 !< Thicknesses of source grid cells real, dimension(:), intent(in) :: x0 !< Source interface positions @@ -275,7 +278,10 @@ subroutine interpolate_grid( n0, h0, x0, ppoly0_E, ppoly0_coefs, & integer, intent(in) :: n1 !< Number of points on target grid real, dimension(:), intent(inout) :: h1 !< Thicknesses of target grid cells real, dimension(:), intent(inout) :: x1 !< Target interface positions + logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions. + ! Local variables + logical :: use_2018_answers ! If true use older, less acccurate expressions. integer :: k ! loop index real :: t ! current interface target density @@ -287,7 +293,8 @@ subroutine interpolate_grid( n0, h0, x0, ppoly0_E, ppoly0_coefs, & ! Find coordinates for interior target values do k = 2,n1 t = target_values(k) - x1(k) = get_polynomial_coordinate ( n0, h0, x0, ppoly0_E, ppoly0_coefs, t, degree ) + x1(k) = get_polynomial_coordinate ( n0, h0, x0, ppoly0_E, ppoly0_coefs, t, degree, & + answers_2018=answers_2018 ) h1(k-1) = x1(k) - x1(k-1) enddo h1(n1) = x1(n1+1) - x1(n1) @@ -320,7 +327,7 @@ subroutine build_and_interpolate_grid(CS, densities, n0, h0, x0, target_values, call regridding_set_ppolys(CS, densities, n0, h0, ppoly0_E, ppoly0_S, ppoly0_C, & degree, h_neglect, h_neglect_edge) call interpolate_grid(n0, h0, x0, ppoly0_E, ppoly0_C, target_values, degree, & - n1, h1, x1) + n1, h1, x1, answers_2018=CS%answers_2018) end subroutine build_and_interpolate_grid !> Given a target value, find corresponding coordinate for given polynomial @@ -340,7 +347,7 @@ end subroutine build_and_interpolate_grid !! It is assumed that the number of cells defining 'grid' and 'ppoly' are the !! same. function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, & - target_value, degree ) result ( x_tgt ) + target_value, degree, answers_2018 ) result ( x_tgt ) ! Arguments integer, intent(in) :: N !< Number of grid cells real, dimension(:), intent(in) :: h !< Grid cell thicknesses @@ -349,6 +356,7 @@ function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, & real, dimension(:,:), intent(in) :: ppoly_coefs !< Coefficients of interpolating polynomials real, intent(in) :: target_value !< Target value to find position for integer, intent(in) :: degree !< Degree of the interpolating polynomials + logical, optional, intent(in) :: answers_2018 !< If true use older, less acccurate expressions. real :: x_tgt !< The position of x_g at which target_value is found. ! Local variables integer :: i, k ! loop indices @@ -363,9 +371,11 @@ function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, & real :: eps ! offset used to get away from ! boundaries real :: grad ! gradient during N-R iterations + logical :: use_2018_answers ! If true use older, less acccurate expressions. eps = NR_OFFSET k_found = -1 + use_2018_answers = .true. ; if (present(answers_2018)) use_2018_answers = answers_2018 ! If the target value is outside the range of all values, we ! force the target coordinate to be equal to the lowest or @@ -441,10 +451,14 @@ function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, & exit endif - numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + & - a(5)*xi0*xi0*xi0*xi0 - target_value - - denominator = a(2) + 2*a(3)*xi0 + 3*a(4)*xi0*xi0 + 4*a(5)*xi0*xi0*xi0 + if (use_2018_answers) then + numerator = a(1) + a(2)*xi0 + a(3)*xi0*xi0 + a(4)*xi0*xi0*xi0 + & + a(5)*xi0*xi0*xi0*xi0 - target_value + denominator = a(2) + 2*a(3)*xi0 + 3*a(4)*xi0*xi0 + 4*a(5)*xi0*xi0*xi0 + else ! These expressions are mathematicaly equivalent but more accurate. + numerator = (a(1) - target_value) + xi0*(a(2) + xi0*(a(3) + xi0*(a(4) + a(5)*xi0))) + denominator = a(2) + xi0*(2.*a(3) + xi0*(3.*a(4) + 4.*a(5)*xi0)) + endif delta = -numerator / denominator @@ -463,7 +477,11 @@ function get_polynomial_coordinate( N, h, x_g, ppoly_E, ppoly_coefs, & if ( xi0 > 1.0 ) then xi0 = 1.0 - grad = a(2) + 2*a(3) + 3*a(4) + 4*a(5) + if (use_2018_answers) then + grad = a(2) + 2*a(3) + 3*a(4) + 4*a(5) + else ! These expressions are mathematicaly equivalent but more accurate. + grad = a(2) + (2.*a(3) + (3.*a(4) + 4.*a(5))) + endif if ( grad == 0.0 ) xi0 = xi0 - eps endif diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index ad9e235b27..4b16730fee 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1180,7 +1180,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call cpu_clock_begin(id_clock_thermo) if (.not.CS%adiabatic) then if (CS%debug) then - call uvchksum("Pre-diabatic [uv]", u, v, G%HI, haloshift=2) + call uvchksum("Pre-diabatic [uv]", u, v, G%HI, haloshift=2, scale=US%L_T_to_m_s) call hchksum(h,"Pre-diabatic h", G%HI, haloshift=1, scale=GV%H_to_m) call uvchksum("Pre-diabatic [uv]h", CS%uhtr, CS%vhtr, G%HI, & haloshift=0, scale=GV%H_to_m*US%L_to_m**2) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index c479550847..8c016b11b0 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -480,7 +480,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & call disable_averaging(CS%diag) if (CS%debug) then - call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym) + call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s) endif call vertvisc_coef(up, vp, h, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp) @@ -670,7 +670,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (CS%debug) then call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym) - call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym) + call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) call hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_m) ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US) call check_redundant("Predictor up ", up, vp, G) @@ -860,7 +860,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (CS%id_v_BT_accel > 0) call post_data(CS%id_v_BT_accel, CS%v_accel_bt, CS%diag) if (CS%debug) then - call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym, vel_scale=1.0) + call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym) call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) call hchksum(h_av, "Corrector avg h", G%HI, haloshift=1, scale=GV%H_to_m) ! call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV, US) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 05f2cac00a..9794070f20 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -1284,20 +1284,22 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, ! surface mass flux maps handles%id_prcme = register_diag_field('ocean_model', 'PRCmE', diag%axesT1, Time, & - 'Net surface water flux (precip+melt+lrunoff+ice calving-evap)', 'kg m-2 s-1',& + 'Net surface water flux (precip+melt+lrunoff+ice calving-evap)', 'kg m-2 s-1', & standard_name='water_flux_into_sea_water', cmor_field_name='wfo', & cmor_standard_name='water_flux_into_sea_water',cmor_long_name='Water Flux Into Sea Water') + ! This diagnostic is rescaled to MKS units when combined. - handles%id_evap = register_diag_field('ocean_model', 'evap', diag%axesT1, Time, & - 'Evaporation/condensation at ocean surface (evaporation is negative)', 'kg m-2 s-1',& - standard_name='water_evaporation_flux', cmor_field_name='evs', & - cmor_standard_name='water_evaporation_flux', & + handles%id_evap = register_diag_field('ocean_model', 'evap', diag%axesT1, Time, & + 'Evaporation/condensation at ocean surface (evaporation is negative)', & + 'kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, & + standard_name='water_evaporation_flux', cmor_field_name='evs', & + cmor_standard_name='water_evaporation_flux', & cmor_long_name='Water Evaporation Flux Where Ice Free Ocean over Sea') ! smg: seaice_melt field requires updates to the sea ice model handles%id_seaice_melt = register_diag_field('ocean_model', 'seaice_melt', & diag%axesT1, Time, 'water flux to ocean from snow/sea ice melting(> 0) or formation(< 0)', & - 'kg m-2 s-1', & + 'kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, & standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics', & cmor_field_name='fsitherm', & cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics',& @@ -1305,6 +1307,7 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, handles%id_precip = register_diag_field('ocean_model', 'precip', diag%axesT1, Time, & 'Liquid + frozen precipitation into ocean', 'kg m-2 s-1') + ! This diagnostic is rescaled to MKS units when combined. handles%id_fprec = register_diag_field('ocean_model', 'fprec', diag%axesT1, Time, & 'Frozen precipitation into ocean', & @@ -1324,32 +1327,39 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T) handles%id_frunoff = register_diag_field('ocean_model', 'frunoff', diag%axesT1, Time, & - 'Frozen runoff (calving) and iceberg melt into ocean', 'kg m-2 s-1', & + 'Frozen runoff (calving) and iceberg melt into ocean', & + units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, & standard_name='water_flux_into_sea_water_from_icebergs', & cmor_field_name='ficeberg', & cmor_standard_name='water_flux_into_sea_water_from_icebergs', & cmor_long_name='Water Flux into Seawater from Icebergs') handles%id_lrunoff = register_diag_field('ocean_model', 'lrunoff', diag%axesT1, Time, & - 'Liquid runoff (rivers) into ocean', 'kg m-2 s-1', & + 'Liquid runoff (rivers) into ocean', & + units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, & standard_name='water_flux_into_sea_water_from_rivers', cmor_field_name='friver', & cmor_standard_name='water_flux_into_sea_water_from_rivers', & cmor_long_name='Water Flux into Sea Water From Rivers') handles%id_net_massout = register_diag_field('ocean_model', 'net_massout', diag%axesT1, Time, & 'Net mass leaving the ocean due to evaporation, seaice formation', 'kg m-2 s-1') + ! This diagnostic is rescaled to MKS units when combined. handles%id_net_massin = register_diag_field('ocean_model', 'net_massin', diag%axesT1, Time, & 'Net mass entering ocean due to precip, runoff, ice melt', 'kg m-2 s-1') + ! This diagnostic is rescaled to MKS units when combined. handles%id_massout_flux = register_diag_field('ocean_model', 'massout_flux', diag%axesT1, Time, & 'Net mass flux of freshwater out of the ocean (used in the boundary flux calculation)', & 'kg m-2', conversion=diag%GV%H_to_kg_m2) + ! This diagnostic is calculated in MKS units. handles%id_massin_flux = register_diag_field('ocean_model', 'massin_flux', diag%axesT1, Time, & 'Net mass flux of freshwater into the ocean (used in boundary flux calculation)', 'kg m-2') + ! This diagnostic is calculated in MKS units. + !========================================================================= - ! area integrated surface mass transport + ! area integrated surface mass transport, all are rescaled to MKS units before area integration. handles%id_total_prcme = register_scalar_field('ocean_model', 'total_PRCmE', Time, diag, & long_name='Area integrated net surface water flux (precip+melt+liq runoff+ice calving-evap)',& diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index f99b6d7f7c..668f185152 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -534,9 +534,10 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ call find_eta(h, tv, G, GV, US, eta) do k=1,nz ; do j=js,je ; do i=is,ie - tmp1(i,j,k) = (eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j) + tmp1(i,j,k) = US%Z_to_m*(eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j) enddo ; enddo ; enddo - vol_tot = US%Z_to_m*reproducing_sum(tmp1, sums=vol_lay) + vol_tot = reproducing_sum(tmp1, sums=vol_lay) + do k=1,nz ; vol_lay(k) = US%m_to_Z * vol_lay(k) ; enddo else do k=1,nz ; do j=js,je ; do i=is,ie tmp1(i,j,k) = H_to_kg_m2 * h(i,j,k) * areaTm(i,j) diff --git a/src/framework/MOM_spatial_means.F90 b/src/framework/MOM_spatial_means.F90 index 829afb851f..85d5ce452b 100644 --- a/src/framework/MOM_spatial_means.F90 +++ b/src/framework/MOM_spatial_means.F90 @@ -43,7 +43,7 @@ function global_area_mean(var, G, scale) do j=js,je ; do i=is,ie tmpForSumming(i,j) = var(i,j) * (scalefac * G%areaT(i,j) * G%mask2dT(i,j)) enddo ; enddo - global_area_mean = reproducing_sum(tmpForSumming) * (G%US%m_to_L**2 * G%IareaT_global) + global_area_mean = reproducing_sum(tmpForSumming) * G%IareaT_global end function global_area_mean @@ -182,17 +182,20 @@ end function global_mass_integral !> Determine the global mean of a field along rows of constant i, returning it !! in a 1-d array using the local indexing. This uses reproducing sums. -subroutine global_i_mean(array, i_mean, G, mask, scale) +subroutine global_i_mean(array, i_mean, G, mask, scale, tmp_scale) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G)), intent(in) :: array !< The variable being averaged real, dimension(SZJ_(G)), intent(out) :: i_mean !< Global mean of array along its i-axis real, dimension(SZI_(G),SZJ_(G)), & - optional, intent(in) :: mask !< An array used for weighting the i-mean - real, optional, intent(in) :: scale !< A rescaling factor for the variable + optional, intent(in) :: mask !< An array used for weighting the i-mean + real, optional, intent(in) :: scale !< A rescaling factor for the output variable + real, optional, intent(in) :: tmp_scale !< A rescaling factor for the internal + !! calculations that is removed from the output ! Local variables type(EFP_type), allocatable, dimension(:) :: asum, mask_sum real :: scalefac ! A scaling factor for the variable. + real :: unscale ! A factor for undoing any internal rescaling before output. real :: mask_sum_r integer :: is, ie, js, je, idg_off, jdg_off integer :: i, j @@ -201,6 +204,10 @@ subroutine global_i_mean(array, i_mean, G, mask, scale) idg_off = G%idg_offset ; jdg_off = G%jdg_offset scalefac = 1.0 ; if (present(scale)) scalefac = scale + unscale = 1.0 + if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then + scalefac = scalefac * tmp_scale ; unscale = 1.0 / tmp_scale + endif ; endif call reset_EFP_overflow_error() allocate(asum(G%jsg:G%jeg)) @@ -253,24 +260,29 @@ subroutine global_i_mean(array, i_mean, G, mask, scale) enddo endif + if (unscale /= 1.0) then ; do j=js,je ; i_mean(j) = unscale*i_mean(j) ; enddo ; endif + deallocate(asum) end subroutine global_i_mean !> Determine the global mean of a field along rows of constant j, returning it !! in a 1-d array using the local indexing. This uses reproducing sums. -subroutine global_j_mean(array, j_mean, G, mask, scale) +subroutine global_j_mean(array, j_mean, G, mask, scale, tmp_scale) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure real, dimension(SZI_(G),SZJ_(G)), intent(in) :: array !< The variable being averaged real, dimension(SZI_(G)), intent(out) :: j_mean !< Global mean of array along its j-axis real, dimension(SZI_(G),SZJ_(G)), & - optional, intent(in) :: mask !< An array used for weighting the j-mean - real, optional, intent(in) :: scale !< A rescaling factor for the variable + optional, intent(in) :: mask !< An array used for weighting the j-mean + real, optional, intent(in) :: scale !< A rescaling factor for the output variable + real, optional, intent(in) :: tmp_scale !< A rescaling factor for the internal + !! calculations that is removed from the output ! Local variables type(EFP_type), allocatable, dimension(:) :: asum, mask_sum real :: mask_sum_r real :: scalefac ! A scaling factor for the variable. + real :: unscale ! A factor for undoing any internal rescaling before output. integer :: is, ie, js, je, idg_off, jdg_off integer :: i, j @@ -278,6 +290,10 @@ subroutine global_j_mean(array, j_mean, G, mask, scale) idg_off = G%idg_offset ; jdg_off = G%jdg_offset scalefac = 1.0 ; if (present(scale)) scalefac = scale + unscale = 1.0 + if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then + scalefac = scalefac * tmp_scale ; unscale = 1.0 / tmp_scale + endif ; endif call reset_EFP_overflow_error() allocate(asum(G%isg:G%ieg)) @@ -330,6 +346,8 @@ subroutine global_j_mean(array, j_mean, G, mask, scale) enddo endif + if (unscale /= 1.0) then ; do i=is,ie ; j_mean(i) = unscale*j_mean(i) ; enddo ; endif + deallocate(asum) end subroutine global_j_mean diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 8ed9a0a4c7..0ddca45c51 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -159,7 +159,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir) call initialize_grid_rotation_angle(G, PF) ! Compute global integrals of grid values for later use in scalar diagnostics ! - call compute_global_grid_integrals(G) + call compute_global_grid_integrals(G, US=US) ! Write out all of the grid data used by this run. if (write_geom) call write_ocean_geometry_file(G, PF, output_dir, US=US) diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index 3d0fe6f1ed..3338f1fedb 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -1145,17 +1145,21 @@ end subroutine set_velocity_depth_min ! ----------------------------------------------------------------------------- !> Pre-compute global integrals of grid quantities (like masked ocean area) for !! later use in reporting diagnostics -subroutine compute_global_grid_integrals(G) - type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid +subroutine compute_global_grid_integrals(G, US) + type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming + real :: area_scale ! A scaling factor for area into MKS units integer :: i,j + area_scale = 1.0 ; if (present(US)) area_scale = US%L_to_m**2 + tmpForSumming(:,:) = 0. G%areaT_global = 0.0 ; G%IareaT_global = 0.0 do j=G%jsc,G%jec ; do i=G%isc,G%iec - tmpForSumming(i,j) = G%areaT(i,j) * G%mask2dT(i,j) + tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j) enddo ; enddo G%areaT_global = reproducing_sum(tmpForSumming) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 03310d70f3..ff08912191 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1889,16 +1889,19 @@ end subroutine set_velocity_depth_max !> Subroutine to pre-compute global integrals of grid quantities for !! later use in reporting diagnostics -subroutine compute_global_grid_integrals(G) +subroutine compute_global_grid_integrals(G, US) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming + real :: area_scale integer :: i,j + area_scale = US%L_to_m**2 tmpForSumming(:,:) = 0. G%areaT_global = 0.0 ; G%IareaT_global = 0.0 do j=G%jsc,G%jec ; do i=G%isc,G%iec - tmpForSumming(i,j) = G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j) + tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j) enddo ; enddo G%areaT_global = reproducing_sum(tmpForSumming) G%IareaT_global = 1. / (G%areaT_global) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 38bf24ee60..a2257369a8 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -188,7 +188,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h call hchksum(MEKE%GM_src, 'MEKE GM_src', G%HI, scale=US%R_to_kg_m3*US%Z_to_m*US%L_to_m**2*US%s_to_T**3) if (associated(MEKE%MEKE)) call hchksum(MEKE%MEKE, 'MEKE MEKE', G%HI, scale=US%L_T_to_m_s**2) call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI, scale=US%s_to_T) - call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=1, scale=GV%H_to_m) + call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=1, scale=GV%H_to_m*US%L_to_m**2) endif sdt = dt*CS%MEKE_dtScale ! Scaled dt to use for time-stepping diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index b4c100dc5d..eb1afb6bb8 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -3,6 +3,7 @@ module MOM_set_diffusivity ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_checksums, only : hchksum_pair use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : diag_ctrl, time_type @@ -344,8 +345,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%useKappaShear) then if (CS%debug) then - call hchksum(u_h, "before calc_KS u_h",G%HI) - call hchksum(v_h, "before calc_KS v_h",G%HI) + call hchksum_pair("before calc_KS [uv]_h", u_h, v_h, G%HI, scale=US%L_T_to_m_s) endif call cpu_clock_begin(id_clock_kappaShear) if (CS%Vertex_shear) then diff --git a/src/parameterizations/vertical/MOM_sponge.F90 b/src/parameterizations/vertical/MOM_sponge.F90 index dd0887845c..6016dbb98b 100644 --- a/src/parameterizations/vertical/MOM_sponge.F90 +++ b/src/parameterizations/vertical/MOM_sponge.F90 @@ -420,7 +420,7 @@ subroutine apply_sponge(h, dt, G, GV, US, ea, eb, CS, Rcv_ml) eta_anom(i,j) = e_D(i,j,k) - CS%Ref_eta_im(j,k) if (CS%Ref_eta_im(j,K) < -G%bathyT(i,j)) eta_anom(i,j) = 0.0 enddo ; enddo - call global_i_mean(eta_anom(:,:), eta_mean_anom(:,K), G) + call global_i_mean(eta_anom(:,:), eta_mean_anom(:,K), G, tmp_scale=US%Z_to_m) enddo if (CS%fldno > 0) allocate(fld_mean_anom(G%isd:G%ied,nz,CS%fldno)) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index e050933dc2..e425629c77 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -485,8 +485,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & else uhh(I) = uhr(I,j,k) endif - !ts2(I) = 0.5*(1.0 + uhh(I)/(hprev(i+1,j,k)+h_neglect)) - CFL(I) = - uhh(I)/(hprev(i+1,j,k)+h_neglect) ! CFL is positive + !ts2(I) = 0.5*(1.0 + uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j))) + CFL(I) = - uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j)) ! CFL is positive else hup = hprev(i,j,k) - G%areaT(i,j)*min_h hlos = MAX(0.0,-uhr(I-1,j,k)) @@ -497,8 +497,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & else uhh(I) = uhr(I,j,k) endif - !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect)) - CFL(I) = uhh(I)/(hprev(i,j,k)+h_neglect) ! CFL is positive + !ts2(I) = 0.5*(1.0 - uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j))) + CFL(I) = uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive endif enddo @@ -573,7 +573,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! Original implementation of PLM !flux_x(I,m) = uhh(I)*(Tr(m)%t(i+1,j,k) - slope_x(i+1,m)*ts2(I)) endif - !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect)) + !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect*G%areaT(i,j))) enddo ; enddo endif ! usePPM @@ -856,8 +856,8 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & else vhh(i,J) = vhr(i,J,k) endif - !ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k)+h_neglect)) - CFL(i) = - vhh(i,J) / (hprev(i,j+1,k)+h_neglect) ! CFL is positive + !ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1)) + CFL(i) = - vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1)) ! CFL is positive else hup = hprev(i,j,k) - G%areaT(i,j)*min_h hlos = MAX(0.0,-vhr(i,J-1,k)) @@ -868,8 +868,8 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & else vhh(i,J) = vhr(i,J,k) endif - !ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k)+h_neglect)) - CFL(i) = vhh(i,J) / (hprev(i,j,k)+h_neglect) ! CFL is positive + !ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j))) + CFL(i) = vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive endif enddo