diff --git a/src/ALE/MOM_hybgen_regrid.F90 b/src/ALE/MOM_hybgen_regrid.F90 index 524f9b8ff2..491693549f 100644 --- a/src/ALE/MOM_hybgen_regrid.F90 +++ b/src/ALE/MOM_hybgen_regrid.F90 @@ -41,8 +41,10 @@ module MOM_hybgen_regrid dp0k, & !< minimum deep z-layer separation [H ~> m or kg m-2] ds0k !< minimum shallow z-layer separation [H ~> m or kg m-2] - real :: coord_scale = 1.0 !< A scaling factor to restores the depth coordinates to values in m - real :: Rho_coord_scale = 1.0 !< A scaling factor to restores the denesity coordinates to values in kg m-3 + real :: coord_scale = 1.0 !< A scaling factor to restores the depth coordinates to + !! values in m [m H-1 ~> 1 or m3 kg-1] + real :: Rho_coord_scale = 1.0 !< A scaling factor to restores the denesity coordinates to + !! values in kg m-3 [kg m-3 R-1 ~> 1] real :: dpns !< depth to start terrain following [H ~> m or kg m-2] real :: dsns !< depth to stop terrain following [H ~> m or kg m-2] @@ -68,7 +70,7 @@ module MOM_hybgen_regrid !! the bottom that certain adjustments can be made in the Hybgen regridding !! code [H ~> m or kg m-2]. In Hycom, this is set to onem (nominally 1 m). real :: h_thin !< A layer thickness below which a layer is considered to be too thin for - !! certain adjustments to be made in the Hybgen regridding code. + !! certain adjustments to be made in the Hybgen regridding code [H ~> m or kg m-2]. !! In Hycom, this is set to onemm (nominally 0.001 m). real :: rho_eps !< A small nonzero density that is used to prevent division by zero @@ -284,7 +286,7 @@ subroutine get_hybgen_regrid_params(CS, nk, ref_pressure, hybiso, nsigma, dp00i, real, optional, intent(out) :: ref_pressure !< Reference pressure for density calculations [R L2 T-2 ~> Pa] real, optional, intent(out) :: hybiso !< Hybgen uses PCM if layer is within hybiso of target density [R ~> kg m-3] integer, optional, intent(out) :: nsigma !< Number of sigma levels used by HYBGEN - real, optional, intent(out) :: dp00i !< Deep isopycnal spacing minimum thickness (m) + real, optional, intent(out) :: dp00i !< Deep isopycnal spacing minimum thickness [H ~> m or kg m-2] real, optional, intent(out) :: qhybrlx !< Fractional relaxation amount per timestep, 0 < qyhbrlx <= 1 [nondim] real, optional, intent(out) :: dp0k(:) !< minimum deep z-layer separation [H ~> m or kg m-2] real, optional, intent(out) :: ds0k(:) !< minimum shallow z-layer separation [H ~> m or kg m-2] @@ -687,8 +689,8 @@ real function cushn(delp, dp0) ! These are derivative nondimensional parameters. ! real, parameter :: cusha = qqmn**2 * (qqmx-1.0) / (qqmx-qqmn)**2 ! real, parameter :: I_qqmn = 1.0 / qqmn - real, parameter :: qq_scale = (qqmx-1.0) / (qqmx-qqmn)**2 - real, parameter :: I_qqmx = 1.0 / qqmx + real, parameter :: qq_scale = (qqmx-1.0) / (qqmx-qqmn)**2 ! A scaling factor based on qqmn and qqmx [nondim] + real, parameter :: I_qqmx = 1.0 / qqmx ! The inverse of qqmx [nondim] ! --- if delp >= qqmx*dp0 >> dp0, cushn returns delp. ! --- if delp <= qqmn*dp0 << -dp0, cushn returns dp0. diff --git a/src/ALE/MOM_hybgen_remap.F90 b/src/ALE/MOM_hybgen_remap.F90 index 5ab3e162db..f97b0e9c62 100644 --- a/src/ALE/MOM_hybgen_remap.F90 +++ b/src/ALE/MOM_hybgen_remap.F90 @@ -126,7 +126,7 @@ subroutine hybgen_ppm_coefs(s, h_src, edges, nk, ns, thin, PCM_lay) real :: da ! Difference between the unlimited scalar edge value estimates [A] real :: a6 ! Scalar field differences that are proportional to the curvature [A] real :: slk, srk ! Differences between adjacent cell averages of scalars [A] - real :: sck ! Scalar differences across a cell. + real :: sck ! Scalar differences across a cell [A] real :: as(nk) ! Scalar field difference across each cell [A] real :: al(nk), ar(nk) ! Scalar field at the left and right edges of a cell [A] real :: h112(nk+1), h122(nk+1) ! Combinations of thicknesses [H ~> m or kg m-2] diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index eeb4590a08..abcd821790 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -1723,7 +1723,7 @@ logical function test_interp(verbose, msg, nsrc, h_src, u_src, ndest, h_dest, u_ ! Local variables real, dimension(ndest+1) :: u_dest ! Interpolated value at destination cell interfaces [A] integer :: k - real :: error + real :: error ! The difference between the evaluated and expected solutions [A] ! Interpolate from src to dest call interpolate_column(nsrc, h_src, u_src, ndest, h_dest, u_dest, .true.) @@ -1760,7 +1760,7 @@ logical function test_reintegrate(verbose, msg, nsrc, h_src, uh_src, ndest, h_de ! Local variables real, dimension(ndest) :: uh_dest ! Reintegrated value on destination cells [A H] integer :: k - real :: error + real :: error ! The difference between the evaluated and expected solutions [A H] ! Interpolate from src to dest call reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, uh_dest) diff --git a/src/ALE/P1M_functions.F90 b/src/ALE/P1M_functions.F90 index b17b35c85c..7889966135 100644 --- a/src/ALE/P1M_functions.F90 +++ b/src/ALE/P1M_functions.F90 @@ -36,7 +36,7 @@ subroutine P1M_interpolation( N, h, u, edge_values, ppoly_coef, h_neglect, answe ! Local variables integer :: k ! loop index - real :: u0_l, u0_r ! edge values (left and right) + real :: u0_l, u0_r ! edge values (left and right) [A] ! Bound edge values (routine found in 'edge_values.F90') call bound_edge_values( N, h, u, edge_values, h_neglect, answer_date=answer_date ) @@ -74,10 +74,10 @@ subroutine P1M_boundary_extrapolation( N, h, u, edge_values, ppoly_coef ) real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly [A] ! Local variables - real :: u0, u1 ! cell averages - real :: h0, h1 ! corresponding cell widths - real :: slope ! retained PLM slope - real :: u0_l, u0_r ! edge values + real :: u0, u1 ! cell averages [A] + real :: h0, h1 ! corresponding cell widths [H] + real :: slope ! retained PLM slope [A] + real :: u0_l, u0_r ! edge values [A] ! ----------------------------------------- ! Left edge value in the left boundary cell diff --git a/src/ALE/PCM_functions.F90 b/src/ALE/PCM_functions.F90 index 4f64e4a96d..f5899339e4 100644 --- a/src/ALE/PCM_functions.F90 +++ b/src/ALE/PCM_functions.F90 @@ -17,11 +17,11 @@ module PCM_functions !! defining 'grid' and 'ppoly'. No consistency check is performed. subroutine PCM_reconstruction( N, u, edge_values, ppoly_coef ) integer, intent(in) :: N !< Number of cells - real, dimension(:), intent(in) :: u !< cell averages + real, dimension(:), intent(in) :: u !< cell averages in arbitrary units [A] real, dimension(:,:), intent(inout) :: edge_values !< Edge value of polynomial, - !! with the same units as u. + !! with the same units as u [A]. real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, - !! with the same units as u. + !! with the same units as u [A]. ! Local variables integer :: k diff --git a/src/ALE/PLM_functions.F90 b/src/ALE/PLM_functions.F90 index bc7f100a04..c0c4516fe2 100644 --- a/src/ALE/PLM_functions.F90 +++ b/src/ALE/PLM_functions.F90 @@ -16,20 +16,21 @@ module PLM_functions contains -!> Returns a limited PLM slope following White and Adcroft, 2008. [units of u] +!> Returns a limited PLM slope following White and Adcroft, 2008, in the same arbitrary +!! units [A] as the input values. !! Note that this is not the same as the Colella and Woodward method. real elemental pure function PLM_slope_wa(h_l, h_c, h_r, h_neglect, u_l, u_c, u_r) - real, intent(in) :: h_l !< Thickness of left cell [units of grid thickness] - real, intent(in) :: h_c !< Thickness of center cell [units of grid thickness] - real, intent(in) :: h_r !< Thickness of right cell [units of grid thickness] - real, intent(in) :: h_neglect !< A negligible thickness [units of grid thickness] - real, intent(in) :: u_l !< Value of left cell [units of u] - real, intent(in) :: u_c !< Value of center cell [units of u] - real, intent(in) :: u_r !< Value of right cell [units of u] + real, intent(in) :: h_l !< Thickness of left cell in arbitrary grid thickness units [H] + real, intent(in) :: h_c !< Thickness of center cell in arbitrary grid thickness units [H] + real, intent(in) :: h_r !< Thickness of right cell in arbitrary grid thickness units [H] + real, intent(in) :: h_neglect !< A negligible thickness [H] + real, intent(in) :: u_l !< Value of left cell in arbitrary units [A] + real, intent(in) :: u_c !< Value of center cell in arbitrary units [A] + real, intent(in) :: u_r !< Value of right cell in arbitrary units [A] ! Local variables real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as - ! differences across the cell [units of u] - real :: u_min, u_max ! Minimum and maximum value across cell [units of u] + ! differences across the cell [A] + real :: u_min, u_max ! Minimum and maximum value across cell [A] ! Side differences sigma_r = u_r - u_c @@ -63,20 +64,21 @@ real elemental pure function PLM_slope_wa(h_l, h_c, h_r, h_neglect, u_l, u_c, u_ end function PLM_slope_wa -!> Returns a limited PLM slope following Colella and Woodward 1984. +!> Returns a limited PLM slope following Colella and Woodward 1984, in the same +!! arbitrary units as the input values [A]. real elemental pure function PLM_slope_cw(h_l, h_c, h_r, h_neglect, u_l, u_c, u_r) - real, intent(in) :: h_l !< Thickness of left cell [units of grid thickness] - real, intent(in) :: h_c !< Thickness of center cell [units of grid thickness] - real, intent(in) :: h_r !< Thickness of right cell [units of grid thickness] - real, intent(in) :: h_neglect !< A negligible thickness [units of grid thickness] - real, intent(in) :: u_l !< Value of left cell [units of u] - real, intent(in) :: u_c !< Value of center cell [units of u] - real, intent(in) :: u_r !< Value of right cell [units of u] + real, intent(in) :: h_l !< Thickness of left cell in arbitrary grid thickness units [H] + real, intent(in) :: h_c !< Thickness of center cell in arbitrary grid thickness units [H] + real, intent(in) :: h_r !< Thickness of right cell in arbitrary grid thickness units [H] + real, intent(in) :: h_neglect !< A negligible thickness [H] + real, intent(in) :: u_l !< Value of left cell in arbitrary units [A] + real, intent(in) :: u_c !< Value of center cell in arbitrary units [A] + real, intent(in) :: u_r !< Value of right cell in arbitrary units [A] ! Local variables real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as - ! differences across the cell [units of u] - real :: u_min, u_max ! Minimum and maximum value across cell [units of u] - real :: h_cn ! Thickness of center cell [units of grid thickness] + ! differences across the cell [A] + real :: u_min, u_max ! Minimum and maximum value across cell [A] + real :: h_cn ! Thickness of center cell [H] h_cn = h_c + h_neglect @@ -117,18 +119,19 @@ real elemental pure function PLM_slope_cw(h_l, h_c, h_r, h_neglect, u_l, u_c, u_ end function PLM_slope_cw -!> Returns a limited PLM slope following Colella and Woodward 1984. +!> Returns a limited PLM slope following Colella and Woodward 1984, in the same +!! arbitrary units as the input values [A]. real elemental pure function PLM_monotonized_slope(u_l, u_c, u_r, s_l, s_c, s_r) - real, intent(in) :: u_l !< Value of left cell [units of u] - real, intent(in) :: u_c !< Value of center cell [units of u] - real, intent(in) :: u_r !< Value of right cell [units of u] - real, intent(in) :: s_l !< PLM slope of left cell [units of u] - real, intent(in) :: s_c !< PLM slope of center cell [units of u] - real, intent(in) :: s_r !< PLM slope of right cell [units of u] + real, intent(in) :: u_l !< Value of left cell in arbitrary units [A] + real, intent(in) :: u_c !< Value of center cell in arbitrary units [A] + real, intent(in) :: u_r !< Value of right cell in arbitrary units [A] + real, intent(in) :: s_l !< PLM slope of left cell [A] + real, intent(in) :: s_c !< PLM slope of center cell [A] + real, intent(in) :: s_r !< PLM slope of right cell [A] ! Local variables - real :: e_r, e_l, edge ! Right, left and temporary edge values [units of u] - real :: almost_two ! The number 2, almost. - real :: slp ! Magnitude of PLM central slope [units of u] + real :: e_r, e_l, edge ! Right, left and temporary edge values [A] + real :: almost_two ! The number 2, almost [nondim] + real :: slp ! Magnitude of PLM central slope [A] almost_two = 2. * ( 1. - epsilon(s_c) ) @@ -155,17 +158,18 @@ real elemental pure function PLM_monotonized_slope(u_l, u_c, u_r, s_l, s_c, s_r) end function PLM_monotonized_slope -!> Returns a PLM slope using h2 extrapolation from a cell to the left. +!> Returns a PLM slope using h2 extrapolation from a cell to the left, in the same +!! arbitrary units as the input values [A]. !! Use the negative to extrapolate from the cell to the right. real elemental pure function PLM_extrapolate_slope(h_l, h_c, h_neglect, u_l, u_c) - real, intent(in) :: h_l !< Thickness of left cell [units of grid thickness] - real, intent(in) :: h_c !< Thickness of center cell [units of grid thickness] - real, intent(in) :: h_neglect !< A negligible thickness [units of grid thickness] - real, intent(in) :: u_l !< Value of left cell [units of u] - real, intent(in) :: u_c !< Value of center cell [units of u] + real, intent(in) :: h_l !< Thickness of left cell in arbitrary grid thickness units [H] + real, intent(in) :: h_c !< Thickness of center cell in arbitrary grid thickness units [H] + real, intent(in) :: h_neglect !< A negligible thickness [H] + real, intent(in) :: u_l !< Value of left cell in arbitrary units [A] + real, intent(in) :: u_c !< Value of center cell in arbitrary units [A] ! Local variables - real :: left_edge ! Left edge value [units of u] - real :: hl, hc ! Left and central cell thicknesses [units of grid thickness] + real :: left_edge ! Left edge value [A] + real :: hl, hc ! Left and central cell thicknesses [H] ! Avoid division by zero for vanished cells hl = h_l + h_neglect @@ -185,24 +189,26 @@ end function PLM_extrapolate_slope !! defining 'grid' and 'ppoly'. No consistency check is performed here. subroutine PLM_reconstruction( N, h, u, edge_values, 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(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A] real, dimension(:,:), intent(inout) :: edge_values !< edge values of piecewise polynomials, - !! with the same units as u. + !! with the same units as u [A]. real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly - !! with the same units as u. + !! with the same units as u [A]. real, optional, intent(in) :: h_neglect !< A negligibly small width for !! the purpose of cell reconstructions - !! in the same units as h + !! in the same units as h [H] ! Local variables - integer :: k ! loop index - real :: u_l, u_r ! left and right cell averages - real :: slope ! retained PLM slope - real :: e_r, edge - real :: almost_one - real, dimension(N) :: slp, mslp - real :: hNeglect + integer :: k ! loop index + real :: u_l, u_r ! left and right cell averages [A] + real :: slope ! retained PLM slope for a normalized cell width [A] + real :: e_r ! The edge value in the neighboring cell [A] + real :: edge ! The projected edge value in the cell [A] + real :: almost_one ! A value that is slightly smaller than 1 [nondim] + real, dimension(N) :: slp ! The first guess at the normalized tracer slopes [A] + real, dimension(N) :: mslp ! The monotonized normalized tracer slopes [A] + real :: hNeglect ! A negligibly small width used in cell reconstructions [H] hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect @@ -265,18 +271,18 @@ end subroutine PLM_reconstruction !! defining 'grid' and 'ppoly'. No consistency check is performed here. subroutine PLM_boundary_extrapolation( N, h, u, edge_values, 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(in) :: h !< cell widths (size N) [H] + real, dimension(:), intent(in) :: u !< cell averages (size N) in arbitrary units [A] real, dimension(:,:), intent(inout) :: edge_values !< edge values of piecewise polynomials, - !! with the same units as u. + !! with the same units as u [A]. real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly - !! with the same units as u. + !! with the same units as u [A]. real, optional, intent(in) :: h_neglect !< A negligibly small width for !! the purpose of cell reconstructions - !! in the same units as h + !! in the same units as h [H] ! Local variables - real :: slope ! retained PLM slope - real :: hNeglect + real :: slope ! retained PLM slope for a normalized cell width [A] + real :: hNeglect ! A negligibly small width used in cell reconstructions [H] hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect diff --git a/src/ALE/PPM_functions.F90 b/src/ALE/PPM_functions.F90 index 805a70d502..ef6841f635 100644 --- a/src/ALE/PPM_functions.F90 +++ b/src/ALE/PPM_functions.F90 @@ -28,7 +28,7 @@ module PPM_functions subroutine PPM_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect, answer_date) integer, intent(in) :: N !< Number of cells real, dimension(N), intent(in) :: h !< Cell widths [H] - real, dimension(N), intent(in) :: u !< Cell averages [A] + real, dimension(N), intent(in) :: u !< Cell averages in arbitrary coordinates [A] real, dimension(N,2), intent(inout) :: edge_values !< Edge values [A] real, dimension(N,3), intent(inout) :: ppoly_coef !< Polynomial coefficients, mainly [A] real, optional, intent(in) :: h_neglect !< A negligibly small width [H] @@ -36,7 +36,7 @@ subroutine PPM_reconstruction( N, h, u, edge_values, ppoly_coef, h_neglect, answ ! Local variables integer :: k ! Loop index - real :: edge_l, edge_r ! Edge values (left and right) + real :: edge_l, edge_r ! Edge values (left and right) [A] ! PPM limiter call PPM_limiter_standard( N, h, u, edge_values, h_neglect, answer_date=answer_date ) @@ -69,9 +69,9 @@ subroutine PPM_limiter_standard( N, h, u, edge_values, h_neglect, answer_date ) ! Local variables integer :: k ! Loop index - real :: u_l, u_c, u_r ! Cell averages (left, center and right) - real :: edge_l, edge_r ! Edge values (left and right) - real :: expr1, expr2 + real :: u_l, u_c, u_r ! Cell averages (left, center and right) [A] + real :: edge_l, edge_r ! Edge values (left and right) [A] + real :: expr1, expr2 ! Temporary expressions [A2] ! Bound edge values call bound_edge_values( N, h, u, edge_values, h_neglect, answer_date=answer_date ) @@ -135,8 +135,8 @@ subroutine PPM_monotonicity( N, u, edge_values ) real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A] ! Local variables - integer :: k ! Loop index - real :: a6,da ! scalar temporaries + integer :: k ! Loop index + real :: a6, da ! Normalized scalar curvature and slope [A] ! Loop on interior cells to impose monotonicity ! Eq. 1.10 of (Colella & Woodward, JCP 84) @@ -195,14 +195,16 @@ subroutine PPM_boundary_extrapolation( N, h, u, edge_values, ppoly_coef, h_negle ! Local variables integer :: i0, i1 - real :: u0, u1 - real :: h0, h1 - real :: a, b, c - real :: u0_l, u0_r - real :: u1_l, u1_r - real :: slope - real :: exp1, exp2 - real :: hNeglect + real :: u0, u1 ! Average concentrations in the two neighboring cells [A] + real :: h0, h1 ! Thicknesses of the two neighboring cells [H] + real :: a, b, c ! An edge value, normalized slope and normalized curvature + ! of a reconstructed distribution [A] + real :: u0_l, u0_r ! Edge values of a neighboring cell [A] + real :: u1_l, u1_r ! Neighboring cell slopes renormalized by the thickness of + ! the cell being worked on [A] + real :: slope ! The normalized slope [A] + real :: exp1, exp2 ! Temporary expressions [A2] + real :: hNeglect ! A negligibly small width used in cell reconstructions [H] hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect diff --git a/src/ALE/PQM_functions.F90 b/src/ALE/PQM_functions.F90 index 1159652858..ef42fb9f01 100644 --- a/src/ALE/PQM_functions.F90 +++ b/src/ALE/PQM_functions.F90 @@ -30,10 +30,10 @@ subroutine PQM_reconstruction( N, h, u, edge_values, edge_slopes, ppoly_coef, h_ ! Local variables integer :: k ! loop index - real :: h_c ! cell width + real :: h_c ! cell width [H] real :: u0_l, u0_r ! edge values (left and right) [A] real :: u1_l, u1_r ! edge slopes (left and right) [A H-1] - real :: a, b, c, d, e ! parabola coefficients + real :: a, b, c, d, e ! quartic fit coefficients [A] ! PQM limiter call PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answer_date=answer_date ) @@ -90,14 +90,15 @@ subroutine PQM_limiter( N, h, u, edge_values, edge_slopes, h_neglect, answer_dat 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 - real :: slope ! retained PLM slope - real :: a, b, c, d, e - real :: alpha1, alpha2, alpha3 - real :: rho, sqrt_rho - real :: gradient1, gradient2 - real :: x1, x2 - real :: hNeglect + 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 :: a, b, c, d, e ! quartic fit coefficients [A] + real :: alpha1, alpha2, alpha3 ! Normalized second derivative coefficients [A] + real :: rho ! A temporary expression [A2] + real :: sqrt_rho ! The square root of rho [A] + real :: gradient1, gradient2 ! Normalized gradients [A] + real :: x1, x2 ! Fractional inflection point positions in a cell [nondim] + real :: hNeglect ! A negligibly small width for the purpose of cell reconstructions [H] hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect @@ -359,13 +360,13 @@ subroutine PQM_boundary_extrapolation( N, h, u, edge_values, ppoly_coef ) real, dimension(:,:), intent(inout) :: ppoly_coef !< Coefficients of polynomial, mainly [A] ! Local variables integer :: i0, i1 - real :: u0, u1 - real :: h0, h1 - real :: a, b, c, d, e - real :: u0_l, u0_r - real :: u1_l, u1_r - real :: slope - real :: exp1, exp2 + real :: u0, u1 ! Successive cell averages [A] + real :: h0, h1 ! Successive cell thicknesses [H] + real :: a, b, c, d, e ! quartic fit coefficients [A] + real :: u0_l, u0_r ! Edge values [A] + real :: u1_l, u1_r ! Edge slopes [A H-1] + real :: slope ! The integrated slope across the cell [A] + real :: exp1, exp2 ! Two temporary expressions [A2] ! ----- Left boundary ----- i0 = 1 @@ -511,19 +512,21 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, edge_values, edge_slopes, ppo integer :: i0, i1 integer :: inflexion_l integer :: inflexion_r - real :: u0, u1, um - real :: h0, h1 - real :: a, b, c, d, e - real :: ar, br, beta - real :: u0_l, u0_r - real :: u1_l, u1_r - real :: u_plm - real :: slope - real :: alpha1, alpha2, alpha3 - real :: rho, sqrt_rho - real :: gradient1, gradient2 - real :: x1, x2 - real :: hNeglect + real :: u0, u1, um ! Successive cell averages [A] + real :: h0, h1 ! Successive cell thicknesses [H] + real :: a, b, c, d, e ! quartic fit coefficients [A] + real :: ar, br ! Temporary variables in [A] + real :: beta ! A rational function coefficient [nondim] + real :: u0_l, u0_r ! Edge values [A] + real :: u1_l, u1_r ! Edge slopes [A H-1] + real :: u_plm ! The integrated piecewise linear method slope [A] + real :: slope ! The integrated slope across the cell [A] + real :: alpha1, alpha2, alpha3 ! Normalized second derivative coefficients [A] + real :: rho ! A temporary expression [A2] + real :: sqrt_rho ! The square root of rho [A] + real :: gradient1, gradient2 ! Normalized gradients [A] + real :: x1, x2 ! Fractional inflection point positions in a cell [nondim] + real :: hNeglect ! A negligibly small width for the purpose of cell reconstructions [H] hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect diff --git a/src/ALE/coord_zlike.F90 b/src/ALE/coord_zlike.F90 index f2ed7f0035..7f284217b2 100644 --- a/src/ALE/coord_zlike.F90 +++ b/src/ALE/coord_zlike.F90 @@ -67,13 +67,15 @@ subroutine build_zstar_column(CS, depth, total_thickness, zInterface, & !! output units), units may be [Z ~> m] or [H ~> m or kg m-2] real, intent(in) :: total_thickness !< Column thickness (positive definite in the same !! units as depth) [Z ~> m] or [H ~> m or kg m-2] - real, dimension(CS%nk+1), intent(inout) :: zInterface !< Absolute positions of interfaces + real, dimension(CS%nk+1), intent(inout) :: zInterface !< Absolute positions of interfaces (in the same + !! units as depth) [Z ~> m] or [H ~> m or kg m-2] real, optional, intent(in) :: z_rigid_top !< The height of a rigid top (positive upward in the same !! units as depth) [Z ~> m] or [H ~> m or kg m-2] - real, optional, intent(in) :: eta_orig !< The actual original height of the top in the same + real, optional, intent(in) :: eta_orig !< The actual original height of the top (in the same !! units as depth) [Z ~> m] or [H ~> m or kg m-2] real, optional, intent(in) :: zScale !< Scaling factor from the target coordinate resolution - !! in Z to desired units for zInterface, perhaps Z_to_H + !! in Z to desired units for zInterface, perhaps Z_to_H, + !! often [nondim] or [H Z-1 ~> 1 or kg m-3] ! Local variables real :: eta ! Free surface height [Z ~> m] or [H ~> m or kg m-2] real :: stretching ! A stretching factor for the coordinate [nondim] diff --git a/src/ALE/polynomial_functions.F90 b/src/ALE/polynomial_functions.F90 index e5c90fe31d..b01e097b83 100644 --- a/src/ALE/polynomial_functions.F90 +++ b/src/ALE/polynomial_functions.F90 @@ -9,7 +9,7 @@ module polynomial_functions contains -!> Pointwise evaluation of a polynomial at x +!> Pointwise evaluation of a polynomial in arbitrary units [A] at x !! !! The polynomial is defined by the coefficients contained in the !! array of the same name, as follows: C(1) + C(2)x + C(3)x^2 + C(4)x^3 + ... @@ -17,12 +17,14 @@ module polynomial_functions !! The number of coefficients is given by ncoef and x !! is the coordinate where the polynomial is to be evaluated. real function evaluation_polynomial( coeff, ncoef, x ) - real, dimension(:), intent(in) :: coeff !< The coefficients of the polynomial + real, dimension(:), intent(in) :: coeff !< The coefficients of the polynomial, in units that + !! vary with the index k as [A H^(k-1)] integer, intent(in) :: ncoef !< The number of polynomial coefficients real, intent(in) :: x !< The position at which to evaluate the polynomial + !! in arbitrary thickness units [H] ! Local variables integer :: k - real :: f ! value of polynomial at x + real :: f ! value of polynomial at x in arbitrary units [A] f = 0.0 do k = 1,ncoef @@ -33,7 +35,8 @@ real function evaluation_polynomial( coeff, ncoef, x ) end function evaluation_polynomial -!> Calculates the first derivative of a polynomial evaluated at a point x +!> Calculates the first derivative of a polynomial evaluated in arbitrary units of [A H-1] +!! at a point x !! !! The polynomial is defined by the coefficients contained in the !! array of the same name, as follows: C(1) + C(2)x + C(3)x^2 + C(4)x^3 + ... @@ -41,12 +44,14 @@ end function evaluation_polynomial !! The number of coefficients is given by ncoef and x !! is the coordinate where the polynomial's derivative is to be evaluated. real function first_derivative_polynomial( coeff, ncoef, x ) - real, dimension(:), intent(in) :: coeff !< The coefficients of the polynomial + real, dimension(:), intent(in) :: coeff !< The coefficients of the polynomial, in units that + !! vary with the index k as [A H^(k-1)] integer, intent(in) :: ncoef !< The number of polynomial coefficients real, intent(in) :: x !< The position at which to evaluate the derivative + !! in arbitrary thickness units [H] ! Local variables integer :: k - real :: f ! value of polynomial at x + real :: f ! value of the derivative at x in [A H-1] f = 0.0 do k = 2,ncoef @@ -57,17 +62,20 @@ real function first_derivative_polynomial( coeff, ncoef, x ) end function first_derivative_polynomial -!> Exact integration of polynomial of degree npoly +!> Exact integration of polynomial of degree npoly in arbitrary units of [A H] !! !! The array of coefficients (Coeff) must be of size npoly+1. real function integration_polynomial( xi0, xi1, Coeff, npoly ) - real, intent(in) :: xi0 !< The lower bound of the integral - real, intent(in) :: xi1 !< The lower bound of the integral - real, dimension(:), intent(in) :: Coeff !< The coefficients of the polynomial + real, intent(in) :: xi0 !< The lower bound of the integral in arbitrary + !! thickness units [H] + real, intent(in) :: xi1 !< The upper bound of the integral in arbitrary + !! thickness units [H] + real, dimension(:), intent(in) :: Coeff !< The coefficients of the polynomial, in units that + !! vary with the index k as [A H^(k-1)] integer, intent(in) :: npoly !< The degree of the polynomial ! Local variables - integer :: k - real :: integral + integer :: k + real :: integral ! The integral of the polynomial over the specified range in [A H] integral = 0.0 diff --git a/src/ALE/regrid_edge_values.F90 b/src/ALE/regrid_edge_values.F90 index 9b574348af..0814c6a907 100644 --- a/src/ALE/regrid_edge_values.F90 +++ b/src/ALE/regrid_edge_values.F90 @@ -27,7 +27,7 @@ module regrid_edge_values !! thickness for sum(h) in edge value inversions real, parameter :: hNeglect_dflt = 1.e-30 !< The default value for cut-off minimum !! thickness for sum(h) in other calculations -real, parameter :: hMinFrac = 1.e-5 !< A minimum fraction for min(h)/sum(h) +real, parameter :: hMinFrac = 1.e-5 !< A minimum fraction for min(h)/sum(h) [nondim] contains @@ -119,7 +119,7 @@ subroutine average_discontinuous_edge_values( N, edge_val ) !! second index is for the two edges of each cell. ! Local variables integer :: k ! loop index - real :: u0_avg ! avg value at given edge + real :: u0_avg ! avg value at given edge [A] ! Loop on interior edges do k = 1,N-1 @@ -231,19 +231,24 @@ subroutine edge_values_explicit_h4( N, h, u, edge_val, h_neglect, answer_date ) ! Local variables real :: h0, h1, h2, h3 ! temporary thicknesses [H] real :: h_min ! A minimal cell width [H] - real :: f1, f2, f3 ! auxiliary variables with various units + real :: f1 ! An auxiliary variable [H] + real :: f2 ! An auxiliary variable [A H] + real :: f3 ! An auxiliary variable [H-1] real :: et1, et2, et3 ! terms the expression for edge values [A H] real :: I_h12 ! The inverse of the sum of the two central thicknesses [H-1] real :: I_h012, I_h123 ! Inverses of sums of three successive thicknesses [H-1] real :: I_den_et2, I_den_et3 ! Inverses of denominators in edge value terms [H-2] - real, dimension(5) :: x ! Coordinate system with 0 at edges [H] - real, dimension(4) :: dz ! A temporary array of limited layer thicknesses [H] - real, dimension(4) :: u_tmp ! A temporary array of cell average properties [A] - real, parameter :: C1_12 = 1.0 / 12.0 - real :: dx ! Difference of successive values of x [H] - real, dimension(4,4) :: A ! values near the boundaries - real, dimension(4) :: B, C - real :: hNeglect ! A negligible thickness in the same units as h. + real, dimension(5) :: x ! Coordinate system with 0 at edges [H] + real, dimension(4) :: dz ! A temporary array of limited layer thicknesses [H] + real, dimension(4) :: u_tmp ! A temporary array of cell average properties [A] + real, parameter :: C1_12 = 1.0 / 12.0 ! A rational constant [nondim] + real :: dx ! Difference of successive values of x [H] + real, dimension(4,4) :: A ! Differences in successive positions raised to various powers, + ! in units that vary with the second (j) index as [H^j] + real, dimension(4) :: B ! The right hand side of the system to solve for C [A H] + real, dimension(4) :: C ! The coefficients of a fit polynomial in units that vary + ! with the index (j) as [A H^(j-1)] + real :: hNeglect ! A negligible thickness in the same units as h [H]. integer :: i, j logical :: use_2018_answers ! If true use older, less accurate expressions. @@ -383,11 +388,11 @@ subroutine edge_values_explicit_h4cw( N, h, u, edge_val, h_neglect ) ! Local variables real :: dp(N) ! Input grid layer thicknesses, but with a minimum thickness [H ~> m or kg m-2] - real :: hNeglect ! A negligible thickness in the same units as h + real :: hNeglect ! A negligible thickness in the same units as h [H] real :: da ! Difference between the unlimited scalar edge value estimates [A] real :: a6 ! Scalar field differences that are proportional to the curvature [A] real :: slk, srk ! Differences between adjacent cell averages of scalars [A] - real :: sck ! Scalar differences across a cell. + real :: sck ! Scalar differences across a cell [A] real :: au(N) ! Scalar field difference across each cell [A] real :: al(N), ar(N) ! Scalar field at the left and right edges of a cell [A] real :: h112(N+1), h122(N+1) ! Combinations of thicknesses [H ~> m or kg m-2] @@ -496,22 +501,26 @@ subroutine edge_values_implicit_h4( N, h, u, edge_val, h_neglect, answer_date ) integer :: i, j ! loop indexes real :: h0, h1 ! cell widths [H] real :: h_min ! A minimal cell width [H] - real :: h0_2, h1_2, h0h1 - real :: h0ph1_2, h0ph1_4 + real :: h0_2, h1_2, h0h1 ! Squares or products of thicknesses [H2] + real :: h0ph1_2 ! The square of a sum of thicknesses [H2] + real :: h0ph1_4 ! The fourth power of a sum of thicknesses [H4] real :: alpha, beta ! stencil coefficients [nondim] real :: I_h2, abmix ! stencil coefficients [nondim] - real :: a, b + real :: a, b ! Combinations of stencil coefficients [nondim] real, dimension(5) :: x ! Coordinate system with 0 at edges [H] - real, parameter :: C1_12 = 1.0 / 12.0 - real, parameter :: C1_3 = 1.0 / 3.0 + real, parameter :: C1_12 = 1.0 / 12.0 ! A rational constant [nondim] + real, parameter :: C1_3 = 1.0 / 3.0 ! A rational constant [nondim] real, dimension(4) :: dz ! A temporary array of limited layer thicknesses [H] real, dimension(4) :: u_tmp ! A temporary array of cell average properties [A] real :: dx ! Differences and averages of successive values of x [H] - real, dimension(4,4) :: Asys ! boundary conditions - real, dimension(4) :: Bsys, Csys + real, dimension(4,4) :: Asys ! Differences in successive positions raised to various powers, + ! in units that vary with the second (j) index as [H^j] + real, dimension(4) :: Bsys ! The right hand side of the system to solve for C [A H] + real, dimension(4) :: Csys ! The coefficients of a fit polynomial in units that vary + ! with the index (j) as [A H^(j-1)] real, dimension(N+1) :: tri_l, & ! tridiagonal system (lower diagonal) [nondim] tri_d, & ! tridiagonal system (middle diagonal) [nondim] - tri_c, & ! tridiagonal system central value, with tri_d = tri_c+tri_l+tri_u + tri_c, & ! tridiagonal system central value [nondim], with tri_d = tri_c+tri_l+tri_u tri_u, & ! tridiagonal system (upper diagonal) [nondim] tri_b, & ! tridiagonal system (right hand side) [A] tri_x ! tridiagonal system (solution vector) [A] @@ -667,7 +676,7 @@ subroutine end_value_h4(dz, u, Csys) real :: I_denom ! The inverse of the denominator some expressions [H-3] real :: I_denB3 ! The inverse of the product of three sums of thicknesses [H-3] real :: min_frac = 1.0e-6 ! The square of min_frac should be much larger than roundoff [nondim] - real, parameter :: C1_3 = 1.0 / 3.0 + real, parameter :: C1_3 = 1.0 / 3.0 ! A rational parameter [nondim] ! These are only used for code verification ! real, dimension(4) :: Atest ! The coefficients of an expression that is being tested. @@ -810,17 +819,21 @@ subroutine edge_slopes_implicit_h3( N, h, u, edge_slopes, h_neglect, answer_date real :: I_h ! Inverses of thicknesses [H-1] real :: alpha, beta ! stencil coefficients [nondim] real :: a, b ! weights of cells [H-1] - real, parameter :: C1_12 = 1.0 / 12.0 + real, parameter :: C1_12 = 1.0 / 12.0 ! A rational parameter [nondim] real, dimension(4) :: dz ! A temporary array of limited layer thicknesses [H] real, dimension(4) :: u_tmp ! A temporary array of cell average properties [A] - real, dimension(5) :: x ! Coordinate system with 0 at edges [H] - real :: dx ! 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(5) :: x ! Coordinate system with 0 at edges [H] + real :: dx ! Differences and averages of successive values of x [H] + real, dimension(4,4) :: Asys ! Differences in successive positions raised to various powers, + ! in units that vary with the second (j) index as [H^j] + real, dimension(4) :: Bsys ! The right hand side of the system to solve for C [A H] + real, dimension(4) :: Csys ! The coefficients of a fit polynomial in units that vary with the + ! index (j) as [A H^(j-1)] + real, dimension(3) :: Dsys ! The coefficients of the first derivative of the fit polynomial + ! in units that vary with the index (j) as [A H^(j-2)] real, dimension(N+1) :: tri_l, & ! tridiagonal system (lower diagonal) [nondim] tri_d, & ! tridiagonal system (middle diagonal) [nondim] - tri_c, & ! tridiagonal system central value, with tri_d = tri_c+tri_l+tri_u + tri_c, & ! tridiagonal system central value [nondim], with tri_d = tri_c+tri_l+tri_u tri_u, & ! tridiagonal system (upper diagonal) [nondim] tri_b, & ! tridiagonal system (right hand side) [A H-1] tri_x ! tridiagonal system (solution vector) [A H-1] @@ -1009,23 +1022,27 @@ subroutine edge_slopes_implicit_h5( N, h, u, edge_slopes, h_neglect, answer_date real :: h01, h01_2 ! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n] real :: h23, h23_2 ! Summed thicknesses to various powers [H^n ~> m^n or kg^n m-2n] real :: hNeglect ! A negligible thickness [H]. - real :: h1_2, h2_2 ! the coefficients of the - real :: h1_3, h2_3 ! tridiagonal system - real :: h1_4, h2_4 ! ... - real :: h1_5, h2_5 ! ... - real :: alpha, beta ! stencil coefficients - 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(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 :: h_Min_Frac = 1.0e-4 + real :: h1_2, h2_2 ! Squares of thicknesses [H2] + real :: h1_3, h2_3 ! Cubes of thicknesses [H3] + real :: h1_4, h2_4 ! Fourth powers of thicknesses [H4] + real :: h1_5, h2_5 ! Fifth powers of thicknesses [H5] + real :: alpha, beta ! stencil coefficients [nondim] + real, dimension(7) :: x ! Coordinate system with 0 at edges in the same units as h [H] + real, parameter :: C1_12 = 1.0 / 12.0 ! A rational parameter [nondim] + real, parameter :: C5_6 = 5.0 / 6.0 ! A rational parameter [nondim] + real :: dx, xavg ! Differences and averages of successive values of x [same units as h] + real, dimension(6,6) :: Asys ! The matrix that is being inverted for a solution, + ! in units that might vary with the second (j) index as [H^j] + real, dimension(6) :: Bsys ! The right hand side of the system to solve for C in various + ! units that sometimes vary with the intex (j) as [H^(j-1)] or [H^j] + ! or might be [A] + real, dimension(6) :: Csys ! The solution to a matrix equation usually [nondim] in this routine. + 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 (rhs) [A H-1] + tri_x ! trid. system (unknowns vector) [A H-1] + real :: h_Min_Frac = 1.0e-4 ! A minimum fractional thickness [nondim] integer :: i, k ! loop indexes hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect @@ -1249,18 +1266,24 @@ subroutine edge_values_implicit_h6( N, h, u, edge_val, h_neglect, answer_date ) real :: hNeglect ! A negligible thickness [H]. real :: h1_2, h2_2, h1_3, h2_3 ! Cell widths raised to the 2nd and 3rd powers [H2] or [H3] real :: h1_4, h2_4, h1_5, h2_5 ! Cell widths raised to the 4th and 5th powers [H4] or [H5] - real :: alpha, beta ! stencil coefficients - 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) - tri_d, & ! trid. system (middle diagonal) - tri_u, & ! trid. system (upper diagonal) - tri_b, & ! trid. system (unknowns vector) - tri_x ! trid. system (rhs) + real :: alpha, beta ! stencil coefficients [nondim] + real, dimension(7) :: x ! Coordinate system with 0 at edges in the same units as h [H] + real, parameter :: C1_12 = 1.0 / 12.0 ! A rational parameter [nondim] + real, parameter :: C5_6 = 5.0 / 6.0 ! A rational parameter [nondim] + real :: dx, xavg ! Differences and averages of successive values of x [H] + real, dimension(6,6) :: Asys ! The matrix that is being inverted for a solution, + ! in units that might vary with the second (j) index as [H^j] + real, dimension(6) :: Bsys ! The right hand side of the system to solve for C in various + ! units that sometimes vary with the intex (j) as [H^(j-1)] or [H^j] + ! or might be [A] + real, dimension(6) :: Csys ! The solution to a matrix equation, which might be [nondim] or the + ! coefficients of a fit polynomial in units that vary with the + ! index (j) as [A H^(j-1)] + 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 (rhs) [A] + tri_x ! trid. system (unknowns vector) [A] integer :: i, k ! loop indexes hNeglect = hNeglect_edge_dflt ; if (present(h_neglect)) hNeglect = h_neglect @@ -1432,16 +1455,16 @@ end subroutine edge_values_implicit_h6 !> Test that A*C = R to within a tolerance, issuing a fatal error with an explanatory message if they do not. subroutine test_line(msg, N, A, C, R, mag, tol) - real, intent(in) :: mag !< The magnitude of leading order terms in this line - integer, intent(in) :: N !< The number of points in the system - real, dimension(4), intent(in) :: A !< One of the two vectors being multiplied - real, dimension(4), intent(in) :: C !< One of the two vectors being multiplied - real, intent(in) :: R !< The expected solution of the equation character(len=*), intent(in) :: msg !< An identifying message for this test - real, optional, intent(in) :: tol !< The fractional tolerance for the two solutions - - real :: sum, sum_mag - real :: tolerance + integer, intent(in) :: N !< The number of points in the system + real, dimension(4), intent(in) :: A !< One of the two vectors being multiplied in arbitrary units [A] + real, dimension(4), intent(in) :: C !< One of the two vectors being multiplied in arbitrary units [B] + real, intent(in) :: R !< The expected solution of the equation [A B] + real, intent(in) :: mag !< The magnitude of leading order terms in this line [A B] + real, optional, intent(in) :: tol !< The fractional tolerance for the sums [nondim] + + real :: sum, sum_mag ! The sum of the products and their magnitude in arbitrary units [A B] + real :: tolerance ! The fractional tolerance for the sums [nondim] character(len=128) :: mesg2 integer :: i diff --git a/src/ALE/regrid_solvers.F90 b/src/ALE/regrid_solvers.F90 index 0655d31062..6e5b3a0cb0 100644 --- a/src/ALE/regrid_solvers.F90 +++ b/src/ALE/regrid_solvers.F90 @@ -18,15 +18,19 @@ module regrid_solvers !! The matrix A must be square, with the first index varing down the column. subroutine solve_linear_system( A, R, X, N, 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] + real, dimension(N,N), intent(inout) :: A !< The matrix being inverted in arbitrary units [A] on + !! input, but internally modified to become nondimensional + !! during the solver. + real, dimension(N), intent(inout) :: R !< system right-hand side in arbitrary units [A B] on + !! input, but internally modified to have units of [B] + !! during the solver + real, dimension(N), intent(inout) :: X !< solution vector in arbitrary units [B] 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. - real :: pivot, I_pivot ! The pivot value and its reciprocal [nondim] - real :: swap_a, swap_b + real, parameter :: eps = 0.0 ! Minimum pivot magnitude allowed [A] + real :: factor ! The factor that eliminates the leading nonzero element in a row [A-1] + real :: pivot, I_pivot ! The pivot value and its reciprocal, in [A] and [A-1] + real :: swap_a, swap_b ! Swap space in various units [various] logical :: found_pivot ! If true, a pivot has been found logical :: old_answers ! If true, use expressions that give the original (2008 through 2018) MOM6 answers integer :: i, j, k @@ -110,15 +114,18 @@ end subroutine solve_linear_system !! The matrix A must be square, with the first index varing along the row. subroutine linear_solver( N, A, R, X ) 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] + real, dimension(N,N), intent(inout) :: A !< The matrix being inverted in arbitrary units [A] on + !! input, but internally modified to become nondimensional + !! during the solver. + real, dimension(N), intent(inout) :: R !< system right-hand side in [A B] on input, but internally + !! modified to have units of [B] during the solver + real, dimension(N), intent(inout) :: X !< solution vector [B] ! Local variables - real, parameter :: eps = 0.0 ! Minimum pivot magnitude allowed - real :: factor ! The factor that eliminates the leading nonzero element in a row. - real :: I_pivot ! The reciprocal of the pivot value [inverse of the input units of a row of A] - real :: swap + real, parameter :: eps = 0.0 ! Minimum pivot magnitude allowed [A] + real :: factor ! The factor that eliminates the leading nonzero element in a row [A-1]. + real :: I_pivot ! The reciprocal of the pivot value [A-1] + real :: swap ! Swap space used in various units [various] integer :: i, j, k ! Loop on rows to transform the problem into multiplication by an upper-right matrix. @@ -175,16 +182,17 @@ end subroutine linear_solver !! (A is made up of lower, middle and upper diagonals) subroutine solve_tridiagonal_system( Al, Ad, Au, R, X, N, 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 - real, dimension(N), intent(in) :: Au !< Matrix upper diagonal - real, dimension(N), intent(in) :: R !< system right-hand side - real, dimension(N), intent(out) :: X !< solution vector + real, dimension(N), intent(in) :: Ad !< Matrix center diagonal in arbitrary units [A] + real, dimension(N), intent(in) :: Al !< Matrix lower diagonal [A] + real, dimension(N), intent(in) :: Au !< Matrix upper diagonal [A] + real, dimension(N), intent(in) :: R !< system right-hand side in arbitrary units [A B] + real, dimension(N), intent(out) :: X !< solution vector in arbitrary units [B] 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 - real :: I_pivot ! The inverse of the most recent pivot + real, dimension(N) :: pivot ! The pivot value [A] + real, dimension(N) :: Al_piv ! The lower diagonal divided by the pivot value [nondim] + real, dimension(N) :: c1 ! Au / pivot for the backward sweep [nondim] + real :: I_pivot ! The inverse of the most recent pivot [A-1] integer :: k ! Loop index logical :: old_answers ! If true, use expressions that give the original (2008 through 2018) MOM6 answers @@ -237,16 +245,16 @@ end subroutine solve_tridiagonal_system !! roundoff compared with (Al+Au), the answers are prone to inaccuracy. subroutine solve_diag_dominant_tridiag( Al, Ac, Au, R, X, N ) integer, intent(in) :: N !< The size of the system - real, dimension(N), intent(in) :: Ac !< Matrix center diagonal offset from Al + Au - real, dimension(N), intent(in) :: Al !< Matrix lower diagonal - real, dimension(N), intent(in) :: Au !< Matrix upper diagonal - real, dimension(N), intent(in) :: R !< system right-hand side - real, dimension(N), intent(out) :: X !< solution vector + real, dimension(N), intent(in) :: Ac !< Matrix center diagonal offset from Al + Au in arbitrary units [A] + real, dimension(N), intent(in) :: Al !< Matrix lower diagonal [A] + real, dimension(N), intent(in) :: Au !< Matrix upper diagonal [A] + real, dimension(N), intent(in) :: R !< system right-hand side in arbitrary units [A B] + real, dimension(N), intent(out) :: X !< solution vector in arbitrary units [B] ! Local variables - real, dimension(N) :: c1 ! Au / pivot for the backward sweep - real :: d1 ! The next value of 1.0 - c1 - real :: I_pivot ! The inverse of the most recent pivot - real :: denom_t1 ! The first term in the denominator of the inverse of the pivot. + real, dimension(N) :: c1 ! Au / pivot for the backward sweep [nondim] + real :: d1 ! The next value of 1.0 - c1 [nondim] + real :: I_pivot ! The inverse of the most recent pivot [A-1] + real :: denom_t1 ! The first term in the denominator of the inverse of the pivot [A] integer :: k ! Loop index ! Factorization and forward sweep, in a form that will never give a division by a diff --git a/src/ALE/remapping_attic.F90 b/src/ALE/remapping_attic.F90 index 534428aaed..be20a27466 100644 --- a/src/ALE/remapping_attic.F90 +++ b/src/ALE/remapping_attic.F90 @@ -46,11 +46,13 @@ module remapping_attic function isPosSumErrSignificant(n1, sum1, n2, sum2) integer, intent(in) :: n1 !< Number of values in sum1 integer, intent(in) :: n2 !< Number of values in sum2 - real, intent(in) :: sum1 !< Sum of n1 values [A] + real, intent(in) :: sum1 !< Sum of n1 values in arbitrary units [A] real, intent(in) :: sum2 !< Sum of n2 values [A] logical :: isPosSumErrSignificant !< True if difference in sums is large ! Local variables - real :: sumErr, allowedErr, eps + real :: sumErr ! The absolutde difference in the sums [A] + real :: allowedErr ! The tolerance for the integrated reconstruction [A] + real :: eps ! A tiny fractional error [nondim] if (sum1<0.) call MOM_error(FATAL,'isPosSumErrSignificant: sum1<0 is not allowed!') if (sum2<0.) call MOM_error(FATAL,'isPosSumErrSignificant: sum2<0 is not allowed!') @@ -73,22 +75,22 @@ end function isPosSumErrSignificant subroutine remapByProjection( n0, h0, u0, ppoly0_E, ppoly0_coefs, & n1, h1, method, u1, h_neglect ) integer, intent(in) :: n0 !< Number of cells in source grid - real, intent(in) :: h0(:) !< Source grid widths (size n0) - real, intent(in) :: u0(:) !< Source cell averages (size n0) - real, intent(in) :: ppoly0_E(:,:) !< Edge value of polynomial - real, intent(in) :: ppoly0_coefs(:,:) !< Coefficients of polynomial + real, intent(in) :: h0(:) !< Source grid widths (size n0) in thickness units [H] + real, intent(in) :: u0(:) !< Source cell averages (size n0) in arbitrary units [A] + real, intent(in) :: ppoly0_E(:,:) !< Edge value of polynomial [A] + real, intent(in) :: ppoly0_coefs(:,:) !< Coefficients of polynomial [A] integer, intent(in) :: n1 !< Number of cells in target grid - real, intent(in) :: h1(:) !< Target grid widths (size n1) + real, intent(in) :: h1(:) !< Target grid widths (size n1) [H] integer, intent(in) :: method !< Remapping scheme to use - real, intent(out) :: u1(:) !< Target cell averages (size n1) + real, intent(out) :: u1(:) !< Target cell averages (size n1) [A] real, optional, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions - !! in the same units as h. + !! in the same units as h [H]. ! Local variables integer :: iTarget - real :: xL, xR ! coordinates of target cell edges + real :: xL, xR ! coordinates of target cell edges [H] integer :: jStart ! Used by integrateReconOnInterval() - real :: xStart ! Used by integrateReconOnInterval() + real :: xStart ! Used by integrateReconOnInterval() [H] ! Loop on cells in target grid (grid1). For each target cell, we need to find ! in which source cells the target cell edges lie. The associated indexes are @@ -120,27 +122,32 @@ end subroutine remapByProjection subroutine remapByDeltaZ( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, dx1, & method, u1, h1, h_neglect ) integer, intent(in) :: n0 !< Number of cells in source grid - real, dimension(:), intent(in) :: h0 !< Source grid sizes (size n0) - real, dimension(:), intent(in) :: u0 !< Source cell averages (size n0) - real, dimension(:,:), intent(in) :: ppoly0_E !< Edge value of polynomial - real, dimension(:,:), intent(in) :: ppoly0_coefs !< Coefficients of polynomial + real, dimension(:), intent(in) :: h0 !< Source grid sizes (size n0) in thickness units [H] + real, dimension(:), intent(in) :: u0 !< Source cell averages (size n0) in arbitrary units [A] + real, dimension(:,:), intent(in) :: ppoly0_E !< Edge value of polynomial [A] + real, dimension(:,:), intent(in) :: ppoly0_coefs !< Coefficients of polynomial [A] integer, intent(in) :: n1 !< Number of cells in target grid - real, dimension(:), intent(in) :: dx1 !< Target grid edge positions (size n1+1) + real, dimension(:), intent(in) :: dx1 !< Target grid edge positions (size n1+1) [H] integer, intent(in) :: method !< Remapping scheme to use - real, dimension(:), intent(out) :: u1 !< Target cell averages (size n1) + real, dimension(:), intent(out) :: u1 !< Target cell averages (size n1) [A] real, dimension(:), & - optional, intent(out) :: h1 !< Target grid widths (size n1) + optional, intent(out) :: h1 !< Target grid widths (size n1) [H] real, optional, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions - !! in the same units as h. + !! in the same units as h [H]. ! Local variables integer :: iTarget - real :: xL, xR ! coordinates of target cell edges - real :: xOld, hOld, uOld - real :: xNew, hNew, h_err - real :: uhNew, hFlux, uAve, fluxL, fluxR + real :: xL, xR ! Coordinates of target cell edges [H] + real :: xOld, xNew ! Edge positions on the old and new grids [H] + real :: hOld, hNew ! Cell thicknesses on the old and new grids [H] + real :: uOld ! A source cell average of u [A] + real :: h_err ! An estimate of the error in the reconstructed thicknesses [H] + real :: uhNew ! Cell integrated u on the new grid [A H] + real :: hFlux ! Width of the remapped volume [H] + real :: uAve ! Target cell average of u [A] + real :: fluxL, fluxR ! Fluxes of u through the two cell faces [A H] integer :: jStart ! Used by integrateReconOnInterval() - real :: xStart ! Used by integrateReconOnInterval() + real :: xStart ! Used by integrateReconOnInterval() [H] ! Loop on cells in target grid. For each cell, iTarget, the left flux is ! the right flux of the cell to the left, iTarget-1. @@ -198,36 +205,35 @@ end subroutine remapByDeltaZ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, & xL, xR, hC, uAve, jStart, xStart, h_neglect ) integer, intent(in) :: n0 !< Number of cells in source grid - real, dimension(:), intent(in) :: h0 !< Source grid sizes (size n0) - real, dimension(:), intent(in) :: u0 !< Source cell averages - real, dimension(:,:), intent(in) :: ppoly0_E !< Edge value of polynomial - real, dimension(:,:), intent(in) :: ppoly0_coefs !< Coefficients of polynomial + real, dimension(:), intent(in) :: h0 !< Source grid sizes (size n0) in thickness units [H] + real, dimension(:), intent(in) :: u0 !< Source cell averages in arbitrary units [A] + real, dimension(:,:), intent(in) :: ppoly0_E !< Edge value of polynomial [A] + real, dimension(:,:), intent(in) :: ppoly0_coefs !< Coefficients of polynomial [A] integer, intent(in) :: method !< Remapping scheme to use - real, intent(in) :: xL !< Left edges of target cell - real, intent(in) :: xR !< Right edges of target cell - real, intent(in) :: hC !< Cell width hC = xR - xL - real, intent(out) :: uAve !< Average value on target cell + real, intent(in) :: xL !< Left edges of target cell [H] + real, intent(in) :: xR !< Right edges of target cell [H] + real, intent(in) :: hC !< Cell width hC = xR - xL [H] + real, intent(out) :: uAve !< Average value on target cell [A] integer, intent(inout) :: jStart !< The index of the cell to start searching from !< On exit, contains index of last cell used - real, intent(inout) :: xStart !< The left edge position of cell jStart + real, intent(inout) :: xStart !< The left edge position of cell jStart [H] !< On first entry should be 0. real, optional, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions - !! in the same units as h + !! in the same units as h [H] ! Local variables integer :: j, k - integer :: jL, jR ! indexes of source cells containing target - ! cell edges - real :: q ! complete integration - real :: xi0, xi1 ! interval of integration (local -- normalized - ! -- coordinates) - real :: x0jLl, x0jLr ! Left/right position of cell jL - real :: x0jRl, x0jRr ! Left/right position of cell jR + integer :: jL, jR ! indexes of source cells containing target cell edges + real :: q ! complete integration [A H] + real :: xi0, xi1 ! interval of integration (local -- normalized -- coordinates) [nondim] + real :: x0jLl, x0jLr ! Left/right position of cell jL [H] + real :: x0jRl, x0jRr ! Left/right position of cell jR [H] real :: hAct ! The distance actually used in the integration - ! (notionally xR - xL) which differs due to roundoff. - real :: x0_2, x1_2, x02px12, x0px1 ! Used in evaluation of integrated polynomials - real :: hNeglect ! A negligible thickness in the same units as h - real, parameter :: r_3 = 1.0/3.0 ! Used in evaluation of integrated polynomials + ! (notionally xR - xL) which differs due to roundoff [H]. + real :: x0_2, x1_2 ! Squares of normalized positions used to evaluate polynomials [nondim] + real :: x0px1, x02px12 ! Sums of normalized positions and their squares [nondim] + real :: hNeglect ! A negligible thickness in the same units as h [H] + real, parameter :: r_3 = 1.0/3.0 ! Used in evaluation of integrated polynomials [nondim] hNeglect = hNeglect_dflt ; if (present(h_neglect)) hNeglect = h_neglect @@ -281,7 +287,7 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefs, method, if ( h0(jL) == 0.0 ) then uAve = 0.5 * ( ppoly0_E(jL,1) + ppoly0_E(jL,2) ) else - !### WHY IS THIS NOT WRITTEN AS xi0 = ( xL - x0jLl ) / h0(jL) ---AJA + ! WHY IS THIS NOT WRITTEN AS xi0 = ( xL - x0jLl ) / h0(jL) ---AJA xi0 = xL / ( h0(jL) + hNeglect ) - x0jLl / ( h0(jL) + hNeglect ) select case ( method ) @@ -540,17 +546,24 @@ logical function remapping_attic_unit_tests(verbose) logical, intent(in) :: verbose !< If true, write results to stdout ! Local variables integer, parameter :: n0 = 4, n1 = 3, n2 = 6 - real :: h0(n0), x0(n0+1), u0(n0) - real :: h1(n1), x1(n1+1), u1(n1), hn1(n1), dx1(n1+1) - real :: h2(n2), x2(n2+1), u2(n2), hn2(n2), dx2(n2+1) + real :: h0(n0), x0(n0+1) ! Test cell widths and edge coordinates [H] + real :: u0(n0) ! Test values for remapping in arbitrary units [A] + real :: h1(n1), x1(n1+1) ! Test cell widths and edge coordinates [H] + real :: u1(n1) ! Test values for remapping [A] + real :: h2(n2), x2(n2+1) ! Test cell widths and edge coordinates [H] + real :: u2(n2) ! Test values for remapping [A] + real :: hn1(n1), hn2(n2) ! Updated grid thicknesses [H] + real :: dx1(n1+1), dx2(n2+1) ! Differences in interface positions [H] data u0 /9., 3., -3., -9./ ! Linear profile, 4 at surface to -4 at bottom data h0 /4*0.75/ ! 4 uniform layers with total depth of 3 data h1 /3*1./ ! 3 uniform layers with total depth of 3 data h2 /6*0.5/ ! 6 uniform layers with total depth of 3 - real, allocatable, dimension(:,:) :: ppoly0_E, ppoly0_S, ppoly0_coefs + real, allocatable, dimension(:,:) :: ppoly0_E, ppoly0_S ! Polynomial edge values [A] + real, allocatable, dimension(:,:) :: ppoly0_coefs ! Polynomial reconstruction coefficients [A] integer :: answer_date ! The vintage of the expressions to test integer :: i, degree - real :: err, h_neglect, h_neglect_edge + real :: err ! Difference between a remapped value and its expected value [A] + real :: h_neglect, h_neglect_edge ! Negligible thicknesses used in remapping [H] logical :: thisTest, v v = verbose