diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index a083402fde..0ae5fb1e92 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -178,6 +178,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) logical :: local_logical logical :: remap_boundary_extrap logical :: init_boundary_extrap + logical :: om4_remap_via_sub_cells type(hybgen_regrid_CS), pointer :: hybgen_regridCS => NULL() ! Control structure for hybgen regridding ! for sharing parameters. @@ -234,6 +235,11 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & default=99991231) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "This selects the remapping algorithm used in OM4 that does not use "//& + "the full reconstruction for the top- and lower-most sub-layers, but instead "//& + "assumes they are always vanished (untrue) and so just uses their edge values. "//& + "We recommend setting this option to false.", default=.true.) call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", CS%answer_date, & "The vintage of the expressions and order of arithmetic to use for remapping. "//& "Values below 20190101 result in the use of older, less accurate expressions "//& @@ -247,12 +253,14 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS) check_reconstruction=check_reconstruction, & check_remapping=check_remapping, & force_bounds_in_subcell=force_bounds_in_subcell, & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, & answer_date=CS%answer_date) call initialize_remapping( CS%vel_remapCS, vel_string, & boundary_extrapolation=init_boundary_extrap, & check_reconstruction=check_reconstruction, & check_remapping=check_remapping, & force_bounds_in_subcell=force_bounds_in_subcell, & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, & answer_date=CS%answer_date) call get_param(param_file, mdl, "PARTIAL_CELL_VELOCITY_REMAP", CS%partial_cell_vel_remap, & @@ -326,6 +334,21 @@ subroutine ALE_set_extrap_boundaries( param_file, CS) call remapping_set_param(CS%remapCS, boundary_extrapolation=remap_boundary_extrap) end subroutine ALE_set_extrap_boundaries +!> Sets the remapping algorithm to that of OM4 +!! +!! The remapping aglorithm used in OM4 made poor assumptions about the reconstructions +!! in the top/bottom layers, namely that they were always vanished and could be +!! represented solely by their upper/lower edge value respectively. +!! Passing .false. here uses the full reconstruction of those top and bottom layers +!! and properly sample those layers. +subroutine ALE_set_OM4_remap_algorithm( CS, om4_remap_via_sub_cells ) + type(ALE_CS), pointer :: CS !< Module control structure + logical, intent(in) :: om4_remap_via_sub_cells !< If true, use OM4 remapping algorithm + + call remapping_set_param(CS%remapCS, om4_remap_via_sub_cells =om4_remap_via_sub_cells ) + +end subroutine ALE_set_OM4_remap_algorithm + !> Initialize diagnostics for the ALE module. subroutine ALE_register_diags(Time, G, GV, US, diag, CS) type(time_type),target, intent(in) :: Time !< Time structure diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index ea7b1a3162..364a322f0f 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -5,7 +5,6 @@ module MOM_remapping ! Original module written by Laurent White, 2008.06.09 use MOM_error_handler, only : MOM_error, FATAL -use MOM_io, only : stdout, stderr use MOM_string_functions, only : uppercase use regrid_edge_values, only : edge_values_explicit_h4, edge_values_implicit_h4 use regrid_edge_values, only : edge_values_explicit_h4cw @@ -38,6 +37,9 @@ module MOM_remapping !> The vintage of the expressions to use for remapping. Values below 20190101 result !! in the use of older, less accurate expressions. integer :: answer_date = 99991231 + !> If true, use the OM4 version of the remapping algorithm that makes poor assumptions + !! about the reconstructions in top and bottom layers of the source grid + logical :: om4_remap_via_sub_cells = .false. end type !> Class to assist in unit tests @@ -111,13 +113,15 @@ 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, answers_2018, answer_date) + check_reconstruction, check_remapping, force_bounds_in_subcell, & + om4_remap_via_sub_cells, answers_2018, answer_date) 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) :: om4_remap_via_sub_cells !< If true, use OM4 remapping algorithm logical, optional, intent(in) :: answers_2018 !< If true use older, less accurate expressions. integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use @@ -136,6 +140,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(om4_remap_via_sub_cells)) then + CS%om4_remap_via_sub_cells = om4_remap_via_sub_cells + endif if (present(answers_2018)) then if (answers_2018) then CS%answer_date = 20181231 @@ -171,6 +178,10 @@ subroutine extract_member_remapping_CS(CS, remapping_scheme, degree, boundary_ex end subroutine extract_member_remapping_CS !> Remaps column of values u0 on grid h0 to grid h1 assuming the top edge is aligned. +!! +!! \todo Remove h_neglect argument by moving into remapping_CS +!! \todo Remove PCM_cell argument by adding new method in Recon1D class +!! \todo Inline the two versions of remap_via_sub_cells() in remapping_core_h() to eliminate remap_via_sub_cells() subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge, PCM_cell) type(remapping_CS), intent(in) :: CS !< Remapping control structure integer, intent(in) :: n0 !< Number of cells on source grid @@ -200,16 +211,26 @@ subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edg hNeglect_edge = 1.0e-10 ; if (present(h_neglect_edge)) hNeglect_edge = h_neglect_edge call build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, ppoly_r_E, ppoly_r_S, iMethod, & - hNeglect, hNeglect_edge, PCM_cell ) + hNeglect, hNeglect_edge, PCM_cell ) - if (CS%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, CS%degree, & - CS%boundary_extrapolation, ppoly_r_coefs, ppoly_r_E) + if (CS%om4_remap_via_sub_cells) then - call remap_via_sub_cells( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, n1, h1, iMethod, & - CS%force_bounds_in_subcell, u1, uh_err ) + if (CS%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, CS%degree, & + CS%boundary_extrapolation, ppoly_r_coefs, ppoly_r_E) - if (CS%check_remapping) call check_remapped_values(n0, h0, u0, ppoly_r_E, CS%degree, ppoly_r_coefs, & - n1, h1, u1, iMethod, uh_err, "remapping_core_h") + ! This calls the OM4 version of the remapmping algorithms + call remap_via_sub_cells_om4( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, n1, h1, iMethod, & + CS%force_bounds_in_subcell, u1, uh_err ) + + if (CS%check_remapping) call check_remapped_values(n0, h0, u0, ppoly_r_E, CS%degree, ppoly_r_coefs, & + n1, h1, u1, iMethod, uh_err, "remapping_core_h") + + else ! i.e. if (CS%om4_remap_via_sub_cells == .false.) + + call remap_via_sub_cells( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, n1, h1, iMethod, & + CS%force_bounds_in_subcell, u1, uh_err ) + + endif end subroutine remapping_core_h @@ -472,8 +493,11 @@ end subroutine check_reconstructions_1d !> Remaps column of n0 values u0 on grid h0 to grid h1 with n1 cells by calculating !! the n0+n1+1 sub-integrals of the intersection of h0 and h1, and the summing the !! appropriate integrals into the h1*u1 values. h0 and h1 must have the same units. -subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, method, & - force_bounds_in_subcell, u1, uh_err, ah_sub, aisub_src, aiss, aise ) +!! +!! This uses a buggy remap_via_sub_cells_om4() that produces the wrong value for +!! bottom layers that do not match the same total depth of the water column. +subroutine remap_via_sub_cells_om4(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, method, & + force_bounds_in_subcell, u1, uh_err) integer, intent(in) :: n0 !< Number of cells in source grid real, intent(in) :: h0(n0) !< Source grid widths (size n0) [H] real, intent(in) :: u0(n0) !< Source cell averages (size n0) [A] @@ -485,10 +509,6 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth logical, intent(in) :: force_bounds_in_subcell !< Force sub-cell values to be bounded real, intent(out) :: u1(n1) !< Target cell averages (size n1) [A] real, intent(out) :: uh_err !< Estimate of bound on error in sum of u*h [A H] - real, optional, intent(out) :: ah_sub(n0+n1+1) !< Overlapping sub-cell thicknesses, h_sub [H] - integer, optional, intent(out) :: aisub_src(n0+n1+1) !< i_sub_src - integer, optional, intent(out) :: aiss(n0) !< isrc_start - integer, optional, intent(out) :: aise(n0) !< isrc_ens ! Local variables real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell [H] real, dimension(n0+n1+1) :: uh_sub ! Integral of u*h over each sub-cell [A H] @@ -502,7 +522,6 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth integer, dimension(n1) :: itgt_end ! Index of last sub-cell within each target cell ! For error checking/debugging logical, parameter :: force_bounds_in_target = .true. ! To fix round-off issues - logical, parameter :: adjust_thickest_subcell = .true. ! To fix round-off conservation issues real :: u02_err ! Integrated reconstruction error estimates [H A] ! Calculate sub-layer thicknesses and indices connecting sub-layers to source and target grids @@ -512,7 +531,7 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth ! Loop over each sub-cell to calculate average/integral values within each sub-cell. ! Uses: h_sub, h0_eff, isub_src ! Sets: u_sub, uh_sub - call remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, & + call remap_src_to_sub_grid_om4(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, & h0_eff, isrc_start, isrc_end, isrc_max, isub_src, & method, force_bounds_in_subcell, u_sub, uh_sub, u02_err) @@ -525,10 +544,61 @@ subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, meth ! Include the error remapping from source to sub-cells in the estimate of total remapping error uh_err = uh_err + u02_err - if (present(ah_sub)) ah_sub(1:n0+n1+1) = h_sub(1:n0+n1+1) - if (present(aisub_src)) aisub_src(1:n0+n1+1) = isub_src(1:n0+n1+1) - if (present(aiss)) aiss(1:n0) = isrc_start(1:n0) - if (present(aise)) aise(1:n0) = isrc_end(1:n0) +end subroutine remap_via_sub_cells_om4 + +!> Remaps column of n0 values u0 on grid h0 to grid h1 with n1 cells by calculating +!! the n0+n1+1 sub-integrals of the intersection of h0 and h1, and the summing the +!! appropriate integrals into the h1*u1 values. h0 and h1 must have the same units. +!! +!! \todo Remove h0_eff which is not needed +!! \todo Undo force_bounds_in_target when switching to Recon1D class +subroutine remap_via_sub_cells( n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h1, method, & + force_bounds_in_subcell, u1, uh_err) + integer, intent(in) :: n0 !< Number of cells in source grid + real, intent(in) :: h0(n0) !< Source grid widths (size n0) [H] + real, intent(in) :: u0(n0) !< Source cell averages (size n0) [A] + real, intent(in) :: ppoly0_E(n0,2) !< 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(n1) !< Target grid widths (size n1) [H] + integer, intent(in) :: method !< Remapping scheme to use + logical, intent(in) :: force_bounds_in_subcell !< Force sub-cell values to be bounded + real, intent(out) :: u1(n1) !< Target cell averages (size n1) [A] + real, intent(out) :: uh_err !< Estimate of bound on error in sum of u*h [A H] + ! Local variables + real, dimension(n0+n1+1) :: h_sub ! Width of each each sub-cell [H] + real, dimension(n0+n1+1) :: uh_sub ! Integral of u*h over each sub-cell [A H] + real, dimension(n0+n1+1) :: u_sub ! Average of u over each sub-cell [A] + integer, dimension(n0+n1+1) :: isub_src ! Index of source cell for each sub-cell + integer, dimension(n0) :: isrc_start ! Index of first sub-cell within each source cell + integer, dimension(n0) :: isrc_end ! Index of last sub-cell within each source cell + integer, dimension(n0) :: isrc_max ! Index of thickest sub-cell within each source cell + real, dimension(n0) :: h0_eff ! Effective thickness of source cells [H] + integer, dimension(n1) :: itgt_start ! Index of first sub-cell within each target cell + integer, dimension(n1) :: itgt_end ! Index of last sub-cell within each target cell + ! For error checking/debugging + logical, parameter :: force_bounds_in_target = .true. ! To fix round-off issues + real :: u02_err ! Integrated reconstruction error estimates [H A] + + ! Calculate sub-layer thicknesses and indices connecting sub-layers to source and target grids + call intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, & + isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + + ! Loop over each sub-cell to calculate average/integral values within each sub-cell. + ! Uses: h_sub, h0_eff, isub_src + ! Sets: u_sub, uh_sub + call remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, & + isrc_start, isrc_end, isrc_max, isub_src, & + method, force_bounds_in_subcell, u_sub, uh_sub, u02_err) + + ! Loop over each target cell summing the integrals from sub-cells within the target cell. + ! Uses: itgt_start, itgt_end, h1, h_sub, uh_sub, u_sub + ! Sets: u1, uh_err + call remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, & + force_bounds_in_target, u1, uh_err) + + ! Include the error remapping from source to sub-cells in the estimate of total remapping error + uh_err = uh_err + u02_err end subroutine remap_via_sub_cells @@ -577,7 +647,6 @@ subroutine intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, & real :: dh ! The width of the sub-cell [H] real :: dh0_eff ! Running sum of source cell thickness [H] ! For error checking/debugging - logical, parameter :: adjust_thickest_subcell = .true. ! To fix round-off conservation issues integer :: i0_last_thick_cell logical :: src_has_volume !< True if h0 has not been consumed logical :: tgt_has_volume !< True if h1 has not been consumed @@ -716,7 +785,10 @@ subroutine intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, & end subroutine intersect_src_tgt_grids !> Remaps column of n0 values u0 on grid h0 to subgrid h_sub -subroutine remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, & +!! +!! This includes an error for the scenario where the source grid is much thicker than +!! the target grid and extrapolation is needed. +subroutine remap_src_to_sub_grid_om4(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, & h0_eff, isrc_start, isrc_end, isrc_max, isub_src, & method, force_bounds_in_subcell, u_sub, uh_sub, u02_err) integer, intent(in) :: n0 !< Number of cells in source grid @@ -747,7 +819,6 @@ subroutine remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, real :: dh0_eff ! Running sum of source cell thickness [H] real :: u0_min(n0), u0_max(n0) !< Min/max of u0 for each source cell [A] ! For error checking/debugging - logical, parameter :: force_bounds_in_target = .true. ! To fix round-off issues logical, parameter :: adjust_thickest_subcell = .true. ! To fix round-off conservation issues integer :: i0_last_thick_cell real :: u_orig ! The original value of the reconstruction in a cell [A] @@ -831,6 +902,147 @@ subroutine remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, enddo endif +end subroutine remap_src_to_sub_grid_om4 + +!> Remaps column of n0 values u0 on grid h0 to subgrid h_sub +subroutine remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, n1, h_sub, & + isrc_start, isrc_end, isrc_max, isub_src, & + method, force_bounds_in_subcell, u_sub, uh_sub, u02_err) + integer, intent(in) :: n0 !< Number of cells in source grid + real, intent(in) :: h0(n0) !< Source grid widths (size n0) [H] + real, intent(in) :: u0(n0) !< Source grid widths (size n0) [H] + real, intent(in) :: ppoly0_E(n0,2) !< 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) :: h_sub(n0+n1+1) !< Overlapping sub-cell thicknesses, h_sub [H] + integer, intent(in) :: isrc_start(n0) !< Index of first sub-cell within each source cell + integer, intent(in) :: isrc_end(n0) !< Index of last sub-cell within each source cell + integer, intent(in) :: isrc_max(n0) !< Index of thickest sub-cell within each source cell + integer, intent(in) :: isub_src(n0+n1+1) !< Index of source cell for each sub-cell + integer, intent(in) :: method !< Remapping scheme to use + logical, intent(in) :: force_bounds_in_subcell !< Force sub-cell values to be bounded + real, intent(out) :: u_sub(n0+n1+1) !< Sub-cell cell averages (size n1) [A] + real, intent(out) :: uh_sub(n0+n1+1) !< Sub-cell cell integrals (size n1) [A H] + real, intent(out) :: u02_err !< Integrated reconstruction error estimates [A H] + ! Local variables + integer :: i_sub ! Index of sub-cell + integer :: i0 ! Index into h0(1:n0), source column + integer :: i_max ! Used to record which sub-cell is the largest contribution of a source cell + real :: dh_max ! Used to record which sub-cell is the largest contribution of a source cell [H] + real :: xa, xb ! Non-dimensional position within a source cell (0..1) [nondim] + real :: dh ! The width of the sub-cell [H] + real :: duh ! The total amount of accumulated stuff (u*h) [A H] + real :: dh0_eff ! Running sum of source cell thickness [H] + real :: u0_min(n0), u0_max(n0) !< Min/max of u0 for each source cell [A] + ! For error checking/debugging + logical, parameter :: adjust_thickest_subcell = .true. ! To fix round-off conservation issues + integer :: i0_last_thick_cell + real :: u_orig ! The original value of the reconstruction in a cell [A] + + i0_last_thick_cell = 0 + do i0 = 1, n0 + u0_min(i0) = min(ppoly0_E(i0,1), ppoly0_E(i0,2)) + u0_max(i0) = max(ppoly0_E(i0,1), ppoly0_E(i0,2)) + if (h0(i0)>0.) i0_last_thick_cell = i0 + enddo + + ! Loop over each sub-cell to calculate average/integral values within each sub-cell. + ! Uses: h_sub, isub_src, h0_eff + ! Sets: u_sub, uh_sub + xa = 0. + dh0_eff = 0. + u02_err = 0. + do i_sub = 1, n0+n1 + + ! Sub-cell thickness from loop above + dh = h_sub(i_sub) + + ! Source cell + i0 = isub_src(i_sub) + + ! Evaluate average and integral for sub-cell i_sub. + ! Integral is over distance dh but expressed in terms of non-dimensional + ! positions with source cell from xa to xb (0 <= xa <= xb <= 1). + dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell + if (h0(i0)>0.) then + xb = dh0_eff / h0(i0) ! This expression yields xa <= xb <= 1.0 + xb = min(1., xb) ! This is only needed when the total target column is wider than the source column + u_sub(i_sub) = average_value_ppoly( n0, u0, ppoly0_E, ppoly0_coefs, method, i0, xa, xb) + else ! Vanished cell + xb = 1. + u_sub(i_sub) = u0(i0) + endif + if (force_bounds_in_subcell) then + ! These next two lines should not be needed but when using PQM we found roundoff + ! can lead to overshoots. These lines sweep issues under the rug which need to be + ! properly .. later. -AJA + u_orig = u_sub(i_sub) + u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) ) + u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) ) + u02_err = u02_err + dh*abs( u_sub(i_sub) - u_orig ) + endif + uh_sub(i_sub) = dh * u_sub(i_sub) + + if (isub_src(i_sub+1) /= i0) then + ! If the next sub-cell is in a different source cell, reset the position counters + dh0_eff = 0. + xa = 0. + else + xa = xb ! Next integral will start at end of last + endif + + enddo + i_sub = n0+n1+1 + ! Sub-cell thickness from loop above + dh = h_sub(i_sub) + + ! Source cell + i0 = isub_src(i_sub) + + ! Evaluate average and integral for sub-cell i_sub. + ! Integral is over distance dh but expressed in terms of non-dimensional + ! positions with source cell from xa to xb (0 <= xa <= xb <= 1). + dh0_eff = dh0_eff + dh ! Cumulative thickness within the source cell + if (h0(i0)>0.) then + xb = dh0_eff / h0(i0) ! This expression yields xa <= xb <= 1.0 + xb = min(1., xb) ! This is only needed when the total target column is wider than the source column + u_sub(i_sub) = average_value_ppoly( n0, u0, ppoly0_E, ppoly0_coefs, method, i0, xa, xb) + else ! Vanished cell + xb = 1. + u_sub(i_sub) = u0(i0) + endif + if (force_bounds_in_subcell) then + ! These next two lines should not be needed but when using PQM we found roundoff + ! can lead to overshoots. These lines sweep issues under the rug which need to be + ! properly .. later. -AJA + u_orig = u_sub(i_sub) + u_sub(i_sub) = max( u_sub(i_sub), u0_min(i0) ) + u_sub(i_sub) = min( u_sub(i_sub), u0_max(i0) ) + u02_err = u02_err + dh*abs( u_sub(i_sub) - u_orig ) + endif + uh_sub(i_sub) = dh * u_sub(i_sub) + + if (adjust_thickest_subcell) then + ! Loop over each source cell substituting the integral/average for the thickest sub-cell (within + ! the source cell) with the residual of the source cell integral minus the other sub-cell integrals + ! aka a genius algorithm for accurate conservation when remapping from Robert Hallberg (@Hallberg-NOAA). + ! Uses: i0_last_thick_cell, isrc_max, h_sub, isrc_start, isrc_end, uh_sub, u0, h0 + ! Updates: uh_sub + do i0 = 1, i0_last_thick_cell + i_max = isrc_max(i0) + dh_max = h_sub(i_max) + if (dh_max > 0.) then + ! duh will be the sum of sub-cell integrals within the source cell except for the thickest sub-cell. + duh = 0. + do i_sub = isrc_start(i0), isrc_end(i0) + if (i_sub /= i_max) duh = duh + uh_sub(i_sub) + enddo + uh_sub(i_max) = u0(i0)*h0(i0) - duh + u02_err = u02_err + max( abs(uh_sub(i_max)), abs(u0(i0)*h0(i0)), abs(duh) ) + endif + enddo + endif + end subroutine remap_src_to_sub_grid !> Remaps column of n0+n1+1 values usub on sub-grid h_sub to targets on grid h1 @@ -1305,7 +1517,8 @@ end subroutine dzFromH1H2 !> Constructor for remapping control structure subroutine initialize_remapping( CS, remapping_scheme, boundary_extrapolation, & - check_reconstruction, check_remapping, force_bounds_in_subcell, answers_2018, answer_date) + check_reconstruction, check_remapping, force_bounds_in_subcell, & + om4_remap_via_sub_cells, answers_2018, answer_date) ! Arguments type(remapping_CS), intent(inout) :: CS !< Remapping control structure character(len=*), intent(in) :: remapping_scheme !< Remapping scheme to use @@ -1313,13 +1526,15 @@ 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) :: om4_remap_via_sub_cells !< If true, use OM4 remapping algorithm logical, optional, intent(in) :: answers_2018 !< If true use older, less accurate expressions. integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use ! 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, answers_2018=answers_2018, answer_date=answer_date) + force_bounds_in_subcell=force_bounds_in_subcell, & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, answers_2018=answers_2018, answer_date=answer_date) end subroutine initialize_remapping @@ -1407,6 +1622,8 @@ logical function remapping_unit_tests(verbose) real :: err ! Errors in the remapped thicknesses [H] or values [A] real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H] type(testing) :: test ! Unit testing convenience functions + integer :: om4 + character(len=4) :: om4_tag call test%set( verbose=verbose ) ! Sets the verbosity flag in test @@ -1414,12 +1631,12 @@ logical function remapping_unit_tests(verbose) h_neglect = hNeglect_dflt h_neglect_edge = hNeglect_dflt ; if (answer_date < 20190101) h_neglect_edge = 1.0e-10 - if (verbose) write(stdout,*) ' ===== MOM_remapping: remapping_unit_tests =================' + if (verbose) write(test%stdout,*) ' ===== MOM_remapping: remapping_unit_tests =================' ! This line carries out tests on some older remapping schemes. call test%test( remapping_attic_unit_tests(verbose), 'attic remapping unit tests' ) - if (verbose) write(stdout,*) ' - - - - - 1st generation tests - - - - -' + if (verbose) write(test%stdout,*) ' - - - - - 1st generation tests - - - - -' call initialize_remapping(CS, 'PPM_H4', answer_date=answer_date) @@ -1466,11 +1683,12 @@ logical function remapping_unit_tests(verbose) call test%real_arr(3, u2, (/3.,-10.5,-12./), 'remap_via_sub_cells() 4') deallocate(h0, u0, h1, u1, h2, u2, ppoly0_E, ppoly0_S, ppoly0_coefs) + call end_remapping(CS) ! =============================================== ! This section tests the reconstruction functions ! =============================================== - if (verbose) write(stdout,*) ' - - - - - reconstruction tests - - - - -' + if (verbose) write(test%stdout,*) ' - - - - - reconstruction tests - - - - -' allocate( ppoly0_coefs(5,6), ppoly0_E(5,2), ppoly0_S(5,2), u2(2) ) @@ -1547,22 +1765,13 @@ logical function remapping_unit_tests(verbose) call test%real_arr(5, ppoly0_coefs(:,2), (/0.,6.,0.,0.,0./), 'Limits PPM: P1') call test%real_arr(5, ppoly0_coefs(:,3), (/0.,-3.,3.,0.,0./), 'Limits PPM: P2') + deallocate(ppoly0_E, ppoly0_S, ppoly0_coefs, u2) + ! ============================================================== - ! This section tests the components of remapping_via_sub_cells() + ! This section tests the components of remap_via_sub_cells() ! ============================================================== - call PLM_reconstruction(4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_E(1:4,:), & - ppoly0_coefs(1:4,:), h_neglect ) - call test%real_arr(4, ppoly0_E(1:4,1), (/5.,5.,3.,1./), 'PPM: left edges h=0110') - call test%real_arr(4, ppoly0_E(1:4,2), (/5.,3.,1.,1./), 'PPM: right edges h=0110') - call remap_via_sub_cells( 4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_E(1:4,:), & - ppoly0_coefs(1:4,:), & - 2, (/1.,1./), INTEGRATION_PLM, .false., u2, err ) - call test%real_arr(2, u2, (/4.,2./), 'PLM: remapped h=0110->h=11') - - deallocate(ppoly0_E, ppoly0_S, ppoly0_coefs, u2) - - if (verbose) write(stdout,*) ' - - - - - interpolation tests - - - - -' + if (verbose) write(test%stdout,*) ' - - - - - remapping algororithm tests - - - - -' ! Test 1: n0=2, n1=3 Maps uniform grids with one extra target layer and no implicitly-vanished interior sub-layers ! h_src = | 3 | 3 | @@ -1579,9 +1788,9 @@ logical function remapping_unit_tests(verbose) 3, (/2., 2., 2./), & ! n1, h1 h_sub, h0_eff, & isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) - if (verbose) write(stdout,*) "intersect_src_tgt_grids test 1: n0=2, n1=3" - if (verbose) write(stdout,*) " h_src = | 3 | 3 |" - if (verbose) write(stdout,*) " h_tgt = | 2 | 2 | 2 |" + if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 1: n0=2, n1=3" + if (verbose) write(test%stdout,*) " h_src = | 3 | 3 |" + if (verbose) write(test%stdout,*) " h_tgt = | 2 | 2 | 2 |" call test%real_arr(6, h_sub, (/0.,2.,1.,1.,2.,0./), 'h_sub') call test%real_arr(2, h0_eff, (/3.,3./), 'h0_eff') call test%int_arr(2, isrc_start, (/1,4/), 'isrc_start') @@ -1607,9 +1816,9 @@ logical function remapping_unit_tests(verbose) 2, (/3., 3./), & ! n1, h1 h_sub, h0_eff, & isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) - if (verbose) write(stdout,*) "intersect_src_tgt_grids test 2: n0=3, n1=2" - if (verbose) write(stdout,*) " h_src = | 2 | 2 | 2 |" - if (verbose) write(stdout,*) " h_tgt = | 3 | 3 |" + if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 2: n0=3, n1=2" + if (verbose) write(test%stdout,*) " h_src = | 2 | 2 | 2 |" + if (verbose) write(test%stdout,*) " h_tgt = | 3 | 3 |" call test%real_arr(6, h_sub, (/0.,2.,1.,1.,2.,0./), 'h_sub') call test%real_arr(3, h0_eff, (/2.,2.,2./), 'h0_eff') call test%int_arr(3, isrc_start, (/1,3,5/), 'isrc_start') @@ -1639,9 +1848,9 @@ logical function remapping_unit_tests(verbose) ! isub_src |1| 1 |2| 2 | 2 |2| call intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, & isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) - if (verbose) write(stdout,*) "intersect_src_tgt_grids test 3: n0=2, n1=3" - if (verbose) write(stdout,*) " h_src = | 2 | 4 |" - if (verbose) write(stdout,*) " h_tgt = | 2 | 2 | 2 |" + if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 3: n0=2, n1=3" + if (verbose) write(test%stdout,*) " h_src = | 2 | 4 |" + if (verbose) write(test%stdout,*) " h_tgt = | 2 | 2 | 2 |" call test%real_arr(6, h_sub, (/0.,2.,0.,2.,2.,0./), 'h_sub') call test%real_arr(2, h0_eff, (/2.,4./), 'h0_eff') call test%int_arr(2, isrc_start, (/1,3/), 'isrc_start') @@ -1659,14 +1868,18 @@ logical function remapping_unit_tests(verbose) call PLM_reconstruction(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect ) call PLM_boundary_extrapolation(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect) allocate(u_sub(n0+n1+1), uh_sub(n0+n1+1)) - call remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, & + call remap_src_to_sub_grid_om4(n0, h0, u0, ppoly0_E, ppoly0_coefs, & n1, h_sub, h0_eff, isrc_start, isrc_end, isrc_max, isub_src, & INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) + call test%real_arr(6, u_sub, (/1.,2.,3.,4.,6.,7./), 'u_sub om4') + call remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, & + n1, h_sub, isrc_start, isrc_end, isrc_max, isub_src, & + INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) + call test%real_arr(6, u_sub, (/1.,2.,3.,4.,6.,7./), 'u_sub') ! h_sub = |0|<- 2 ->|0|<- 2 ->|<- 2 ->|0| ! u_sub = |1| 2 |3| 4 | 6 |7| ! h_tgt = |<- 2 ->|<- 2 ->|<- 2 ->| ! u_tgt = | 2 | 4 | 6 | - call test%real_arr(6, u_sub, (/1.,2.,3.,4.,6.,7./), 'u_sub') call remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, & .false., u1, u02_err) call test%real_arr(3, u1, (/2.,4.,6./), 'u1') @@ -1695,9 +1908,9 @@ logical function remapping_unit_tests(verbose) ! isub_src |1| 1 |2| 2 | 2 | 2 | call intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, & isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) - if (verbose) write(stdout,*) "intersect_src_tgt_grids test 4: n0=2, n1=3" - if (verbose) write(stdout,*) " h_src = | 2 | 4 |" - if (verbose) write(stdout,*) " h_tgt = | 2 | 2 | 1 |" + if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 4: n0=2, n1=3" + if (verbose) write(test%stdout,*) " h_src = | 2 | 4 |" + if (verbose) write(test%stdout,*) " h_tgt = | 2 | 2 | 1 |" call test%real_arr(6, h_sub, (/0.,2.,0.,2.,1.,1./), 'h_sub') call test%real_arr(2, h0_eff, (/2.,3./), 'h0_eff') call test%int_arr(2, isrc_start, (/1,3/), 'isrc_start') @@ -1707,18 +1920,18 @@ logical function remapping_unit_tests(verbose) call test%int_arr(3, itgt_end, (/3,4,5/), 'itgt_end') call test%int_arr(6, isub_src, (/1,1,2,2,2,2/), 'isub_src') allocate(ppoly0_coefs(n0,2), ppoly0_E(n0,2), ppoly0_S(n0,2)) -! ! h_src = |<- 2 ->|<- 4 ->| -! ! h_sub = |0|<- 2 ->|0|<- 2 ->|<-1->|<-1->| -! ! u_src = | 2 | 5 | -! ! edge = |1 3|3 7| -! ! u_sub = |1| 2 |3| 4 | 5.5 | 6.5 | -! call PLM_reconstruction(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect ) -! call PLM_boundary_extrapolation(n0, h0, u0, poly0_E, ppoly0_coefs, h_neglect) + ! h_src = |<- 2 ->|<- 4 ->| + ! h_sub = |0|<- 2 ->|0|<- 2 ->|<-1->|<-1->| + ! u_src = | 2 | 5 | + ! edge = |1 3|3 7| + ! u_sub = |1| 2 |3| 4 | 5.5 | 6.5 | + call PLM_reconstruction(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect ) + call PLM_boundary_extrapolation(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect) allocate(u_sub(n0+n1+1), uh_sub(n0+n1+1)) -! call remap_src_to_sub_grid(2, (/2.,4./), (/2.,5./), ppoly0_E, ppoly0_coefs, & -! 3, h_sub, h0_eff, isrc_start, isrc_end, isrc_max, isub_src, & -! INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) -! call test%real_arr(6, u_sub, (/1.,2.,3.,4.,5.5,6.5/), 'u_sub') + call remap_src_to_sub_grid(2, (/2.,4./), (/2.,5./), ppoly0_E, ppoly0_coefs, & + 3, h_sub, isrc_start, isrc_end, isrc_max, isub_src, & + INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) + call test%real_arr(6, u_sub, (/1.,2.,3.,4.,5.5,6.5/), 'u_sub') deallocate( ppoly0_coefs, ppoly0_E, ppoly0_S, u_sub, uh_sub, h0, u0, h1, u1) deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) @@ -1741,9 +1954,9 @@ logical function remapping_unit_tests(verbose) ! isub_src |1| 1 |2| 2 | 3 | 3 | call intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, & isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) - if (verbose) write(stdout,*) "intersect_src_tgt_grids test 5: n0=3, n1=2" - if (verbose) write(stdout,*) " h_src = | 2 | 2 | 1 |" - if (verbose) write(stdout,*) " h_tgt = | 2 | 4 |" + if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 5: n0=3, n1=2" + if (verbose) write(test%stdout,*) " h_src = | 2 | 2 | 1 |" + if (verbose) write(test%stdout,*) " h_tgt = | 2 | 4 |" call test%real_arr(6, h_sub, (/0.,2.,0.,2.,1.,1./), 'h_sub') call test%real_arr(3, h0_eff, (/2.,2.,1./), 'h0_eff') call test%int_arr(3, isrc_start, (/1,3,5/), 'isrc_start') @@ -1761,14 +1974,18 @@ logical function remapping_unit_tests(verbose) call PLM_reconstruction(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect ) call PLM_boundary_extrapolation(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect) allocate(u_sub(n0+n1+1), uh_sub(n0+n1+1)) - call remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, & + call remap_src_to_sub_grid_om4(n0, h0, u0, ppoly0_E, ppoly0_coefs, & n1, h_sub, h0_eff, isrc_start, isrc_end, isrc_max, isub_src, & INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) + call test%real_arr(6, u_sub, (/1.,2.,3.,4.,5.5,6./), 'u_sub om4') + call remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, & + n1, h_sub, isrc_start, isrc_end, isrc_max, isub_src, & + INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) + call test%real_arr(6, u_sub, (/1.,2.,3.,4.,5.5,6./), 'u_sub') ! h_sub = |0|<- 2 ->|0|<- 2 ->|< 1 >|< 1 >| ! u_sub = |1| 2 |3| 4 | 5.5 | 6 | ! h_tgt = |<- 2 ->|<- 4 ->| ! u_tgt = | 2 | 4 7/8 | - call test%real_arr(6, u_sub, (/1.,2.,3.,4.,5.5,6./), 'u_sub') call remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, & .false., u1, u02_err) call test%real_arr(2, u1, (/2.,4.875/), 'u1') @@ -1797,9 +2014,9 @@ logical function remapping_unit_tests(verbose) ! isub_src |1| 1 |1| 1 |2|3|3| 3 |3| call intersect_src_tgt_grids( n0, h0, n1, h1, h_sub, h0_eff, & isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) - if (verbose) write(stdout,*) "intersect_src_tgt_grids test 6: n0=3, n1=5" - if (verbose) write(stdout,*) " h_src = | 2 |0| 2 |" - if (verbose) write(stdout,*) " h_tgt = | 1 |0| 1 |0| 2 |" + if (verbose) write(test%stdout,*) "intersect_src_tgt_grids test 6: n0=3, n1=5" + if (verbose) write(test%stdout,*) " h_src = | 2 |0| 2 |" + if (verbose) write(test%stdout,*) " h_tgt = | 1 |0| 1 |0| 2 |" call test%real_arr(9, h_sub, (/0.,1.,0.,1.,0.,0.,0.,2.,0./), 'h_sub') call test%real_arr(3, h0_eff, (/2.,0.,2./), 'h0_eff') call test%int_arr(3, isrc_start, (/1,5,6/), 'isrc_start') @@ -1817,14 +2034,18 @@ logical function remapping_unit_tests(verbose) call PLM_reconstruction(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect ) call PLM_boundary_extrapolation(n0, h0, u0, ppoly0_E, ppoly0_coefs, h_neglect) allocate(u_sub(n0+n1+1), uh_sub(n0+n1+1)) - call remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, & + call remap_src_to_sub_grid_om4(n0, h0, u0, ppoly0_E, ppoly0_coefs, & n1, h_sub, h0_eff, isrc_start, isrc_end, isrc_max, isub_src, & INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) + call test%real_arr(9, u_sub, (/1.,1.5,2.,2.5,3.,3.,3.,4.,5./), 'u_sub om4') + call remap_src_to_sub_grid(n0, h0, u0, ppoly0_E, ppoly0_coefs, & + n1, h_sub, isrc_start, isrc_end, isrc_max, isub_src, & + INTEGRATION_PLM, .false., u_sub, uh_sub, u02_err) + call test%real_arr(9, u_sub, (/1.,1.5,2.,2.5,3.,3.,3.,4.,5./), 'u_sub') ! h_sub = |0|< 1 ->|0|< 1 >|0|0|0|<- 2 ->|0| ! u_sub = |1| 1.5 |2| 2.5 |3|3|3| 4 |5| ! h_tgt = |<- 1 ->|0|<- 1 ->|0|<- 2 ->| ! u_tgt = | 1.5 |2| 2.5 |3| 4 | - call test%real_arr(9, u_sub, (/1.,1.5,2.,2.5,3.,3.,3.,4.,5./), 'u_sub') call remap_sub_to_tgt_grid(n0, n1, h1, h_sub, u_sub, uh_sub, itgt_start, itgt_end, & .false., u1, u02_err) call test%real_arr(5, u1, (/1.5,2.,2.5,3.,4./), 'u1') @@ -1834,10 +2055,129 @@ logical function remapping_unit_tests(verbose) deallocate( ppoly0_coefs, ppoly0_E, ppoly0_S, u_sub, uh_sub, h0, u0, h1, u1) deallocate( h_sub, h0_eff, isrc_start, isrc_end, isrc_max, itgt_start, itgt_end, isub_src ) + ! ============================================================ + ! This section tests remap_via_sub_cells() + ! ============================================================ + if (verbose) write(test%stdout,*) '- - - - - - - - - - remap_via_sub_cells() tests - - - - - - - - -' + + allocate(ppoly0_E(4,2), ppoly0_S(4,2), ppoly0_coefs(4,3), u2(2)) + + ! These tests has vanished top and bottom layers with a linear profile + call PLM_reconstruction(4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_E, ppoly0_coefs, 0.) + ! Reconstruction has piecewise constant 5 and 1 for top and bottom layers respectively + ! but which are vanished so give no wieght to integrals. + ! Interface values are 5, 5, 3, 1, 1. + call test%real_arr(4, ppoly0_E(1:4,1), (/5.,5.,3.,1./), 'PPM: left edges h=0110') + call test%real_arr(4, ppoly0_E(1:4,2), (/5.,3.,1.,1./), 'PPM: right edges h=0110') + + ! Remapping to just the two interior layers yields the same values as u_src(2:3) + call remap_via_sub_cells_om4(4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_E, ppoly0_coefs, & + 2, (/1.,1./), INTEGRATION_PLM, .false., u2, err) + call test%real_arr(2, u2, (/4.,2./), 'PLM: remapped h=0110->h=11 om4') + call remap_via_sub_cells(4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_E, ppoly0_coefs, & + 2, (/1.,1./), INTEGRATION_PLM, .false., u2, err) + call test%real_arr(2, u2, (/4.,2./), 'PLM: remapped h=0110->h=11') + + ! Remapping to two layers that are deeper. For the bottom layer of thickness 4, + ! the first 1/4 has average 2, the remaining 3/4 has the bottom edge value or 1 + ! yield ing and average or 1.25 + call remap_via_sub_cells_om4(4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_E, ppoly0_coefs, & + 2, (/1.,4./), INTEGRATION_PLM, .false., u2, err) + call test%real_arr(2, u2, (/4.,1.25/), 'PLM: remapped h=0110->h=14 om4') + call remap_via_sub_cells(4, (/0.,1.,1.,0./), (/5.,4.,2.,1./), ppoly0_E, ppoly0_coefs, & + 2, (/1.,4./), INTEGRATION_PLM, .false., u2, err) + call test%real_arr(2, u2, (/4.,1.25/), 'PLM: remapped h=0110->h=14') + + ! Remapping to two layers with lowest layer not reach the bottom. + ! Here, the bottom layer samples top half of source yeilding 2.5. + ! Note: OM4 used the value as if the target layer was the same thickness as source. + call remap_via_sub_cells_om4(4, (/0.,4.,4.,0./), (/5.,4.,2.,1./), ppoly0_E, ppoly0_coefs, & + 2, (/4.,2./), INTEGRATION_PLM, .false., u2, err) + call test%real_arr(2, u2, (/4.,2./), 'PLM: remapped h=0440->h=42 om4 (with known bug)') + call remap_via_sub_cells(4, (/0.,4.,4.,0./), (/5.,4.,2.,1./), ppoly0_E, ppoly0_coefs, & + 2, (/4.,2./), INTEGRATION_PLM, .false., u2, err) + call test%real_arr(2, u2, (/4.,2.5/), 'PLM: remapped h=0440->h=42') + + ! Remapping to two layers with no layers sampling the bottom source layer + ! The first layer samples the top half of u1, yielding 4.5 + ! The second layer samples the next quarter of u1, yielding 3.75 + call remap_via_sub_cells_om4(4, (/0.,5.,5.,0./), (/5.,4.,2.,1./), ppoly0_E, ppoly0_coefs, & + 2, (/2.,2./), INTEGRATION_PLM, .false., u2, err) + call test%real_arr(2, u2, (/4.5,3.5/), 'PLM: remapped h=0880->h=21 om4 (with known bug)') + call remap_via_sub_cells(4, (/0.,4.,4.,0./), (/5.,4.,2.,1./), ppoly0_E, ppoly0_coefs, & + 2, (/2.,1./), INTEGRATION_PLM, .false., u2, err) + call test%real_arr(2, u2, (/4.5,3.75/), 'PLM: remapped h=0440->h=21') + + deallocate(ppoly0_E, ppoly0_S, ppoly0_coefs, u2) + + ! ============================================================ + ! This section tests remapping_core_h() + ! ============================================================ + if (verbose) write(test%stdout,*) '- - - - - - - - - - remapping_core_h() tests - - - - - - - - -' + + ! Profile 0: 8 layers, 1x top/2x bottom vanished, and the rest with thickness 1.0, total depth 5, u(z) = 1 + z + n0 = 8 + allocate( h0(n0), u0(n0) ) + h0 = (/0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0/) + u0 = (/1.0, 1.5, 2.5, 3.5, 4.5, 5.5, 6.0, 6.0/) + allocate( u1(8) ) + + call initialize_remapping(CS, 'PLM', answer_date=99990101) + + do om4 = 0, 1 + if ( om4 == 0 ) then + CS%om4_remap_via_sub_cells = .false. + om4_tag(:) = ' ' + else + CS%om4_remap_via_sub_cells = .true. + om4_tag(:) = ' om4' + endif + + ! Unchanged grid + call remapping_core_h( CS, n0, h0, u0, 8, [0.,1.,1.,1.,1.,1.,0.,0.], u1, 1.e-17, 1.e-2) + call test%real_arr(8, u1, (/1.0,1.5,2.5,3.5,4.5,5.5,6.0,6.0/), 'PLM: remapped h=01111100->h=01111100'//om4_tag) + + ! Removing vanished layers (unchanged values for non-vanished layers, layer centers 0.5, 1.5, 2.5, 3.5, 4.5) + call remapping_core_h( CS, n0, h0, u0, 5, [1.,1.,1.,1.,1.], u1, 1.e-17, 1.e-2) + call test%real_arr(5, u1, (/1.5,2.5,3.5,4.5,5.5/), 'PLM: remapped h=01111100->h=11111'//om4_tag) + + ! Remapping to variable thickness layers (layer centers 0.25, 1.0, 2.25, 4.0) + call remapping_core_h( CS, n0, h0, u0, 4, [0.5,1.,1.5,2.], u1, 1.e-17, 1.e-2) + call test%real_arr(4, u1, (/1.25,2.,3.25,5./), 'PLM: remapped h=01111100->h=h1t2'//om4_tag) + + ! Remapping to variable thickness + vanished layers (layer centers 0.25, 1.0, 1.5, 2.25, 4.0) + call remapping_core_h( CS, n0, h0, u0, 6, [0.5,1.,0.,1.5,2.,0.], u1, 1.e-17, 1.e-2) + call test%real_arr(6, u1, (/1.25,2.,2.5,3.25,5.,6./), 'PLM: remapped h=01111100->h=h10t20'//om4_tag) + + ! Remapping to deeper water column (layer centers 0.75, 2.25, 3., 5., 8.) + call remapping_core_h( CS, n0, h0, u0, 5, [1.5,1.5,0.,4.,2.], u1, 1.e-17, 1.e-2) + call test%real_arr(5, u1, (/1.75,3.25,4.,5.5,6./), 'PLM: remapped h=01111100->h=tt02'//om4_tag) + + ! Remapping to slightly shorter water column (layer centers 0.5, 1.5, 2.5,, 3.5, 4.25) + call remapping_core_h( CS, n0, h0, u0, 5, [1.,1.,1.,1.,0.5], u1, 1.e-17, 1.e-2) + if ( om4 == 0 ) then + call test%real_arr(5, u1, (/1.5,2.5,3.5,4.5,5.25/), 'PLM: remapped h=01111100->h=1111h') + else + call test%real_arr(5, u1, (/1.5,2.5,3.5,4.5,5.5/), 'PLM: remapped h=01111100->h=1111h om4 (known bug)') + endif + + ! Remapping to much shorter water column (layer centers 0.25, 0.5, 1.) + call remapping_core_h( CS, n0, h0, u0, 3, [0.5,0.,1.], u1, 1.e-17, 1.e-2) + if ( om4 == 0 ) then + call test%real_arr(3, u1, (/1.25,1.5,2./), 'PLM: remapped h=01111100->h=h01') + else + call test%real_arr(3, u1, (/1.25,1.5,1.875/), 'PLM: remapped h=01111100->h=h01 om4 (known bug)') + endif + + enddo ! om4 + + call end_remapping(CS) + deallocate( h0, u0, u1 ) + ! ============================================================ ! This section tests interpolation and reintegration functions ! ============================================================ - if (verbose) write(stdout,*) '- - - - - - - - - - interpolation tests - - - - - - - - -' + if (verbose) write(test%stdout,*) '- - - - - - - - - - interpolation tests - - - - - - - - -' call test_interp(test, 'Identity: 3 layer', & 3, (/1.,2.,3./), (/1.,2.,3.,4./), & @@ -1875,7 +2215,7 @@ logical function remapping_unit_tests(verbose) 3, (/1.,2.,4./), (/1.,2.,4.,8./), & 4, (/0.,2.,6.,0./), (/0.,1.,3.,8.,0./) ) - if (verbose) write(stdout,*) ' - - - - - reintegration tests - - - - -' + if (verbose) write(test%stdout,*) ' - - - - - reintegration tests - - - - -' call test_reintegrate(test, 'Identity: 3 layer', & 3, (/1.,2.,3./), (/-5.,2.,1./), & diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 8394735cb9..3674e63c31 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -441,6 +441,7 @@ subroutine open_boundary_config(G, US, param_file, OBC) real :: Lscale_in, Lscale_out ! parameters controlling tracer values at the boundaries [L ~> m] integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: check_reconstruction, check_remapping, force_bounds_in_subcell + logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm character(len=64) :: remappingScheme ! This include declares and sets the variable "version". # include "version_variable.h" @@ -695,10 +696,15 @@ subroutine open_boundary_config(G, US, param_file, OBC) "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date) + call get_param(param_file, mdl, "OBC_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. "//& + "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& + "We recommend setting this option to false.", default=.true.) allocate(OBC%remap_CS) call initialize_remapping(OBC%remap_CS, remappingScheme, boundary_extrapolation = .false., & check_reconstruction=check_reconstruction, check_remapping=check_remapping, & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, & force_bounds_in_subcell=force_bounds_in_subcell, answer_date=OBC%remap_answer_date) endif ! OBC%number_of_segments > 0 diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index fd8057c38f..e5996826e2 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1580,6 +1580,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag logical :: better_speed_est ! If true, use a more robust estimate of the first ! mode wave speed as the starting point for iterations. logical :: split ! True if using the barotropic-baroclinic split algorithm + logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for calculating the EBT structure ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_diagnostics" ! This module's name. @@ -1617,6 +1618,10 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, & "If true, use a more robust estimate of the first mode wave speed as the "//& "starting point for iterations.", default=.true.) + call get_param(param_file, mdl, "INTWAVE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//& + "See REMAPPING_USE_OM4_SUBCELLS for details. "//& + "We recommend setting this option to false.", default=.true.) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & default=99991231) @@ -1858,7 +1863,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag (CS%id_cg_ebt>0) .or. (CS%id_Rd_ebt>0) .or. (CS%id_p_ebt>0)) then call wave_speed_init(CS%wave_speed, remap_answer_date=remap_answer_date, & better_speed_est=better_speed_est, min_speed=wave_speed_min, & - wave_speed_tol=wave_speed_tol) + wave_speed_tol=wave_speed_tol, om4_remap_via_sub_cells=om4_remap_via_sub_cells) endif CS%id_mass_wt = register_diag_field('ocean_model', 'mass_wt', diag%axesT1, Time, & diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 5caf47a51c..8ee271f315 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -1611,7 +1611,8 @@ end subroutine tridiag_det !> Initialize control structure for MOM_wave_speed subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth, remap_answers_2018, & - remap_answer_date, better_speed_est, min_speed, wave_speed_tol, c1_thresh) + remap_answer_date, better_speed_est, om4_remap_via_sub_cells, & + min_speed, wave_speed_tol, c1_thresh) type(wave_speed_CS), intent(inout) :: CS !< Wave speed control struct logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent !! barotropic mode instead of the first baroclinic mode. @@ -1630,6 +1631,8 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de !! forms of the same remapping expressions. logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first !! mode speed as the starting point for iterations. + logical, optional, intent(in) :: om4_remap_via_sub_cells !< Use the OM4-era ramap_via_sub_cells + !! for calculating the EBT structure real, optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed !! below which 0 is returned [L T-1 ~> m s-1]. real, optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the @@ -1656,6 +1659,7 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de ! The remap_answers_2018 argument here is irrelevant, because remapping is hard-coded to use PLM. call initialize_remapping(CS%remapping_CS, 'PLM', boundary_extrapolation=.false., & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, & answer_date=CS%remap_answer_date) end subroutine wave_speed_init diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index b3194af3d8..541e349d29 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -3251,6 +3251,7 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) ! answers from 2018, while higher values use more robust ! forms of the same remapping expressions. integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. + logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for diagnostics character(len=8) :: this_pe character(len=240) :: doc_file, doc_file_dflt, doc_path character(len=240), allocatable :: diag_coords(:) @@ -3282,6 +3283,10 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & default=99991231) + call get_param(param_file, mdl, "DIAG_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "If true, use the OM4 remapping-via-subcells algorithm for diagnostics. "//& + "See REMAPPING_USE_OM4_SUBCELLS for details. "//& + "We recommend setting this option to false.", default=.true.) call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, & "The vintage of the expressions and order of arithmetic to use for remapping. "//& "Values below 20190101 result in the use of older, less accurate expressions "//& @@ -3311,7 +3316,7 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) allocate(diag_cs%diag_remap_cs(diag_cs%num_diag_coords)) ! Initialize each diagnostic vertical coordinate do i=1, diag_cs%num_diag_coords - call diag_remap_init(diag_cs%diag_remap_cs(i), diag_coords(i), answer_date=remap_answer_date, GV=GV) + call diag_remap_init(diag_cs%diag_remap_cs(i), diag_coords(i), om4_remap_via_sub_cells, remap_answer_date, GV) enddo deallocate(diag_coords) endif diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index bbefa3808b..e8e6a756e9 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -92,6 +92,7 @@ module MOM_diag_remap !! vertical extents in [Z ~> m] for remapping extensive variables integer :: interface_axes_id = 0 !< Vertical axes id for remapping at interfaces integer :: layer_axes_id = 0 !< Vertical axes id for remapping on layers + logical :: om4_remap_via_sub_cells !< Use the OM4-era ramap_via_sub_cells integer :: answer_date !< The vintage of the order of arithmetic and expressions !! to use for remapping. Values below 20190101 recover !! the answers from 2018, while higher values use more @@ -102,10 +103,11 @@ module MOM_diag_remap contains !> Initialize a diagnostic remapping type with the given vertical coordinate. -subroutine diag_remap_init(remap_cs, coord_tuple, answer_date, GV) +subroutine diag_remap_init(remap_cs, coord_tuple, om4_remap_via_sub_cells, answer_date, GV) type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure character(len=*), intent(in) :: coord_tuple !< A string in form of !! MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME + logical, intent(in) :: om4_remap_via_sub_cells !< Use the OM4-era ramap_via_sub_cells integer, intent(in) :: answer_date !< The vintage of the order of arithmetic and expressions !! to use for remapping. Values below 20190101 recover !! the answers from 2018, while higher values use more @@ -127,6 +129,7 @@ subroutine diag_remap_init(remap_cs, coord_tuple, answer_date, GV) remap_cs%configured = .false. remap_cs%initialized = .false. remap_cs%used = .false. + remap_cs%om4_remap_via_sub_cells = om4_remap_via_sub_cells remap_cs%answer_date = answer_date remap_cs%nz = 0 @@ -309,6 +312,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe if (.not. remap_cs%initialized) then ! Initialize remapping and regridding on the first call call initialize_remapping(remap_cs%remap_cs, 'PPM_IH4', boundary_extrapolation=.false., & + om4_remap_via_sub_cells=remap_cs%om4_remap_via_sub_cells, & answer_date=remap_cs%answer_date) remap_cs%initialized = .true. endif diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index c18752c83d..3e71a98f55 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2504,6 +2504,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just character(len=64) :: remappingScheme real :: tempAvg ! Spatially averaged temperatures on a layer [C ~> degC] real :: saltAvg ! Spatially averaged salinities on a layer [S ~> ppt] + logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm (only used if useALEremapping) logical :: do_conv_adj, ignore integer :: nPoints integer :: id_clock_routine, id_clock_ALE @@ -2583,6 +2584,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=just_read.or.(.not.GV%Boussinesq)) + call get_param(PF, mdl, "Z_INIT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "If true, use the OM4 remapping-via-subcells algorithm for initialization. "//& + "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& + "We recommend setting this option to false.", default=.true.) if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701) endif call get_param(PF, mdl, "HOR_REGRID_ANSWER_DATE", hor_regrid_answer_date, & @@ -2753,7 +2758,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ! Build the target grid (and set the model thickness to it) call ALE_initRegridding( GV, US, G%max_depth, PF, mdl, regridCS ) ! sets regridCS - call initialize_remapping( remapCS, remappingScheme, boundary_extrapolation=.false., answer_date=remap_answer_date ) + call initialize_remapping( remapCS, remappingScheme, boundary_extrapolation=.false., & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=remap_answer_date ) ! Now remap from source grid to target grid, first setting reconstruction parameters if (remap_general) then diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index cac8a5cd6c..bafa5d8c36 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -96,6 +96,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ ! remapping cell reconstructions [Z ~> m] real :: dz_neglect_edge ! A negligibly small vertical layer extent used in ! remapping edge value calculations [Z ~> m] + logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm integer :: nPoints ! The number of valid input data points in a column integer :: id_clock_routine, id_clock_ALE integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. @@ -137,6 +138,10 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=.not.GV%Boussinesq) + call get_param(PF, mdl, "Z_INIT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "If true, use the OM4 remapping-via-subcells algorithm for initialization. "//& + "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& + "We recommend setting this option to false.", default=.true.) if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701) endif call get_param(PF, mdl, "HOR_REGRID_ANSWER_DATE", hor_regrid_answer_date, & @@ -174,7 +179,8 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ allocate( dzSrc(isd:ied,jsd:jed,kd) ) allocate( hSrc(isd:ied,jsd:jed,kd) ) ! Set parameters for reconstructions - call initialize_remapping( remapCS, remapScheme, boundary_extrapolation=.false., answer_date=remap_answer_date ) + call initialize_remapping( remapCS, remapScheme, boundary_extrapolation=.false., & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=remap_answer_date ) ! Next we initialize the regridding package so that it knows about the target grid do j = js, je ; do i = is, ie diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 8453ceb497..9275555afc 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -183,6 +183,7 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS) character(len=80) :: remap_scheme character(len=80) :: bias_correction_file, inc_file integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. + logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm if (associated(CS)) call MOM_error(FATAL, 'Calling oda_init with associated control structure') allocate(CS) @@ -320,8 +321,10 @@ subroutine init_oda(Time, G, GV, US, diag_CS, CS) call get_param(PF, 'oda_driver', "REGRIDDING_COORDINATE_MODE", coord_mode, & "Coordinate mode for vertical regridding.", & default="ZSTAR", fail_if_missing=.false.) + call get_param(PF, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call initialize_regridding(CS%regridCS, CS%GV, CS%US, dG%max_depth,PF,'oda_driver',coord_mode,'','') - call initialize_remapping(CS%remapCS,remap_scheme) + call initialize_remapping(CS%remapCS, remap_scheme, om4_remap_via_sub_cells=om4_remap_via_sub_cells) call set_regrid_params(CS%regridCS, min_thickness=0.) isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index d54d34506e..94d09554c2 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -143,6 +143,8 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h, nz_data, re real :: dt, dt_therm ! Model timesteps [T ~> s] character(len=256) :: mesg character(len=64) :: remapScheme + logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm + if (.not.associated(CS)) then call MOM_error(WARNING, "initialize_oda_incupd called without an associated "// & "control structure.") @@ -195,6 +197,8 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h, nz_data, re "When defined, the incoming oda_incupd data are "//& "assumed to be on the model horizontal grid " , & default=.true.) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) CS%nz = GV%ke @@ -236,7 +240,7 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h, nz_data, re ! Call the constructor for remapping control structure !### Revisit this hard-coded answer_date. call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=bndExtrapolation, & - answer_date=20190101) + om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=20190101) end subroutine initialize_oda_incupd diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 5b9ce4934c..225781ce0c 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -2550,6 +2550,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) real, dimension(:,:), allocatable :: ridge_temp ! array for temporary storage of flags ! of cells with double-reflecting ridges [nondim] logical :: use_int_tides, use_temperature + logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for calculating the EBT structure real :: IGW_c1_thresh ! A threshold first mode internal wave speed below which all higher ! mode speeds are not calculated but simply assigned a speed of 0 [L T-1 ~> m s-1]. real :: kappa_h2_factor ! A roughness scaling factor [nondim] @@ -2726,6 +2727,10 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) "mode speeds are not calculated but are simply reported as 0. This must be "//& "non-negative for the wave_speeds routine to be used.", & units="m s-1", default=0.01, scale=US%m_s_to_L_T) + call get_param(param_file, mdl, "INTWAVE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//& + "See REMAPPING_USE_OM4_SUBCELLS for details. "//& + "We recommend setting this option to false.", default=.true.) call get_param(param_file, mdl, "UNIFORM_TEST_CG", CS%uniform_test_cg, & "If positive, a uniform group velocity of internal tide for test case", & @@ -3107,7 +3112,8 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) enddo ! Initialize the module that calculates the wave speeds. - call wave_speed_init(CS%wave_speed, c1_thresh=IGW_c1_thresh) + call wave_speed_init(CS%wave_speed, c1_thresh=IGW_c1_thresh, & + om4_remap_via_sub_cells=om4_remap_via_sub_cells) end subroutine internal_tides_init diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index cd7e235274..7db0dc4e90 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1190,6 +1190,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) ! mode wave speed as the starting point for iterations. real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale ! temperature variance [nondim] + logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for calculating the EBT structure ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name. @@ -1593,10 +1594,14 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, & "If true, use a more robust estimate of the first mode wave speed as the "//& "starting point for iterations.", default=.true.) + call get_param(param_file, mdl, "EBT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//& + "See REMAPPING_USE_OM4_SUBCELLS for details. "//& + "We recommend setting this option to false.", default=.true.) call wave_speed_init(CS%wave_speed, use_ebt_mode=CS%Resoln_use_ebt, & mono_N2_depth=N2_filter_depth, remap_answer_date=remap_answer_date, & better_speed_est=better_speed_est, min_speed=wave_speed_min, & - wave_speed_tol=wave_speed_tol) + om4_remap_via_sub_cells=om4_remap_via_sub_cells, wave_speed_tol=wave_speed_tol) endif ! Leith parameters diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index f1739485d6..af33806a90 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -190,6 +190,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, logical :: data_h_to_Z logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. + logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v if (associated(CS)) then @@ -234,6 +235,10 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=.not.GV%Boussinesq) if (.not.GV%Boussinesq) CS%remap_answer_date = max(CS%remap_answer_date, 20230701) + call get_param(param_file, mdl, "SPONGE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "If true, use the OM4 remapping-via-subcells algorithm for ALE sponge. "//& + "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& + "We recommend setting this option to false.", default=.true.) call get_param(param_file, mdl, "HOR_REGRID_ANSWER_DATE", CS%hor_regrid_answer_date, & "The vintage of the order of arithmetic for horizontal regridding. "//& @@ -284,6 +289,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, ! Call the constructor for remapping control structure call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=bndExtrapolation, & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, & answer_date=CS%remap_answer_date) call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, & @@ -468,6 +474,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, US, param_file, CS, I logical :: use_sponge logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. + logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm integer :: i, j, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v if (associated(CS)) then @@ -510,6 +517,10 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, US, param_file, CS, I "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date) + call get_param(param_file, mdl, "SPONGE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "If true, use the OM4 remapping-via-subcells algorithm for ALE sponge. "//& + "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& + "We recommend setting this option to false.", default=.true.) call get_param(param_file, mdl, "HOR_REGRID_ANSWER_DATE", CS%hor_regrid_answer_date, & "The vintage of the order of arithmetic for horizontal regridding. "//& "Dates before 20190101 give the same answers as the code did in late 2018, "//& @@ -549,6 +560,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, US, param_file, CS, I ! Call the constructor for remapping control structure call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=bndExtrapolation, & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, & answer_date=CS%remap_answer_date) call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, & "The total number of columns where sponges are applied at h points.", like_default=.true.) diff --git a/src/tracer/MOM_hor_bnd_diffusion.F90 b/src/tracer/MOM_hor_bnd_diffusion.F90 index 13e91e8973..b6714148ea 100644 --- a/src/tracer/MOM_hor_bnd_diffusion.F90 +++ b/src/tracer/MOM_hor_bnd_diffusion.F90 @@ -89,6 +89,7 @@ logical function hor_bnd_diffusion_init(Time, G, GV, US, param_file, diag, diaba ! local variables character(len=80) :: string ! Temporary strings logical :: boundary_extrap ! controls if boundary extrapolation is used in the HBD code + logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for HBD logical :: debug !< If true, write verbose checksums for debugging purposes if (ASSOCIATED(CS)) then @@ -142,10 +143,15 @@ logical function hor_bnd_diffusion_init(Time, G, GV, US, param_file, diag, diaba "for vertical remapping for all variables. "//& "It can be one of the following schemes: "//& trim(remappingSchemesDoc), default=remappingDefaultScheme) + call get_param(param_file, mdl, "HBD_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "If true, use the OM4 remapping-via-subcells algorithm for horizontal boundary diffusion. "//& + "See REMAPPING_USE_OM4_SUBCELLS for details. "//& + "We recommend setting this option to false.", default=.true.) ! GMM, TODO: add HBD params to control optional arguments in initialize_remapping. - call initialize_remapping( CS%remap_CS, string, boundary_extrapolation = boundary_extrap ,& - check_reconstruction=.false., check_remapping=.false.) + call initialize_remapping( CS%remap_CS, string, boundary_extrapolation=boundary_extrap, & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, & + check_reconstruction=.false., check_remapping=.false.) call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) call get_param(param_file, mdl, "DEBUG", debug, default=.false., do_not_log=.true.) call get_param(param_file, mdl, "HBD_DEBUG", CS%debug, & @@ -849,8 +855,9 @@ logical function near_boundary_unit_tests( verbose ) allocate(CS) ! fill required fields in CS CS%linear=.false. - call initialize_remapping( CS%remap_CS, 'PLM', boundary_extrapolation=.true. ,& - check_reconstruction=.true., check_remapping=.true.) + call initialize_remapping( CS%remap_CS, 'PLM', boundary_extrapolation=.true., & + om4_remap_via_sub_cells=.true., & ! ### see fail below when using fixed remapping alg. + check_reconstruction=.true., check_remapping=.true.) call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) CS%H_subroundoff = 1.0E-20 CS%debug=.false. @@ -1040,6 +1047,7 @@ logical function near_boundary_unit_tests( verbose ) call hbd_grid_test(SURFACE, hbl_L, hbl_R, h_L, h_R, CS) call fluxes_layer_method(SURFACE, nk, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, & khtr_u, F_layer, 1., 1., CS%hbd_u_kmax(1,1), CS%hbd_grd_u(1,1,:), CS) + ! ### This test fails when om4_remap_via_sub_cells=.false. near_boundary_unit_tests = near_boundary_unit_tests .or. & test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.0,-4.0/) ) diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 402a008244..6b60769144 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -152,6 +152,7 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, logical :: debug ! If true, write verbose checksums for debugging purposes. logical :: boundary_extrap ! Indicate whether high-order boundary !! extrapolation should be used within boundary cells. + logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm if (associated(CS)) then call MOM_error(FATAL, "neutral_diffusion_init called with associated control structure.") @@ -232,8 +233,13 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=.not.GV%Boussinesq) + call get_param(param_file, mdl, "NDIFF_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + "If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. "//& + "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& + "We recommend setting this option to false.", default=.true.) if (.not.GV%Boussinesq) CS%remap_answer_date = max(CS%remap_answer_date, 20230701) call initialize_remapping( CS%remap_CS, string, boundary_extrapolation=boundary_extrap, & + om4_remap_via_sub_cells=om4_remap_via_sub_cells, & answer_date=CS%remap_answer_date ) call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg) call get_param(param_file, mdl, "NEUTRAL_POS_METHOD", CS%neutral_pos_method, &