Skip to content

Commit

Permalink
+Eliminated h_neglect argument to remapping_core_h
Browse files Browse the repository at this point in the history 
  Eliminated the previously optional h_neglect and h_neglect_edge arguments to
remapping_core_h and remapping_core_w, and added optional h_neglect and
h_neglect_edge arguments to initialize_remapping for storage in the remapping
control structures.  The answer_date, h_neglect and h_neglect_edge arguments
to ALE_remap_scalar were also eliminated, as they are no longer used.
The new routine open_boundary_setup_vert was added to manage these changes.
A new verticalGrid_type argument was added to wave_speed_init.

  The h_neglect argument to 29 routines was made non-optional, because there is
no way to provide a consistent hard-coded default for a dimensional parameter.
The (mostly low-level) routines where this change to a non-optional h_neglect
was made include  build_reconstructions_1d, P1M_interpolation,
P3M_interpolation, P3M_boundary_extrapolation, PLM_reconstruction,
PLM_boundary_extrapolation, PPM_reconstruction, PPM_limiter_standard,
PPM_boundary_extrapolation, PQM_reconstruction, PQM_limiter,
PQM_boundary_extrapolation_v1, build_hycom1_column, build_rho_column,
bound_edge_values, edge_values_explicit_h4, edge_values_explicit_h4cw,
edge_values_implicit_h4, edge_slopes_implicit_h3, edge_slopes_implicit_h5,
edge_values_implicit_h6, regridding_set_ppolys, build_and_interpolate_grid,
remapByProjection, remapByDeltaZ, and integrateReconOnInterval.

  In two modules (MOM_open_boundary and MOM_wave_speed), two separate copies of
remapping_CS variables were needed to accommodate different negligible
thicknesses or units in different remapping calls.

  In those cases that also have an optional h_neglect_edge argument the default
is now set to h_neglect.  Improperly hard-coded dimensional default values for
h_neglect or h_neglect_edge were eliminated in 16 places as a result of this
change, but they were added back in 2 unit testing routines (one of these is
archaic and unused) that do not use exercise dimensional rescaling tests.

  Because the previously optional arguments were already present in all calls
from the dev/gfdl or main branches of the MOM6 code, and because the places
where arguments have been removed they are balanced by equivalent arguments
added to initialize_remapping, all answers are bitwise identical in the
regression tests, but this change in interfaces could impact other code that is
not in the main branch of MOM6.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Aug 7, 2024
1 parent f1ba822 commit 3d1b6cd
Show file tree
Hide file tree
Showing 27 changed files with 474 additions and 576 deletions.
7 changes: 5 additions & 2 deletions config_src/drivers/timing_tests/time_MOM_remapping.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,9 @@ program time_MOM_remapping
real, dimension(nschemes) :: tmin ! Shortest time for a call [s]
real, dimension(nschemes) :: tmax ! Longest time for a call [s]
real, dimension(:,:), allocatable :: u0, u1 ! Source/target values [arbitrary but same units as each other]
real, dimension(:,:), allocatable :: h0, h1 ! Source target thicknesses [0..1]
real, dimension(:,:), allocatable :: h0, h1 ! Source target thicknesses [0..1] [nondim]
real :: start, finish ! Times [s]
real :: h_neglect ! A negligible thickness [nondim]
real :: h0sum, h1sum ! Totals of h0 and h1 [nondim]
integer :: ij, k, isamp, iter, ischeme ! Indices and counters
integer :: seed_size ! Number of integers used by seed
Expand Down Expand Up @@ -50,6 +51,7 @@ program time_MOM_remapping
h0(:,ij) = h0(:,ij) / h0sum
h1(:,ij) = h1(:,ij) / h1sum
enddo
h_neglect = 1.0-30

! Loop over many samples of timing loop to collect statistics
tmean(:) = 0.
Expand All @@ -59,7 +61,8 @@ program time_MOM_remapping
do isamp = 1, nsamp
! Time reconstruction + remapping
do ischeme = 1, nschemes
call initialize_remapping(CS, remapping_scheme=trim(scheme_labels(ischeme)))
call initialize_remapping(CS, remapping_scheme=trim(scheme_labels(ischeme)), &
h_neglect=h_neglect, h_neglect_edge=h_neglect)
call cpu_time(start)
do iter = 1, nits ! Make many passes to reduce sampling error
do ij = 1, nij ! Calling nij times to make similar to cost in MOM_ALE()
Expand Down
105 changes: 27 additions & 78 deletions src/ALE/MOM_ALE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,6 +181,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
logical :: om4_remap_via_sub_cells
type(hybgen_regrid_CS), pointer :: hybgen_regridCS => NULL() ! Control structure for hybgen regridding
! for sharing parameters.
real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2]

if (associated(CS)) then
call MOM_error(WARNING, "ALE_init called with an associated "// &
Expand Down Expand Up @@ -248,20 +249,30 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
default=default_answer_date, do_not_log=.not.GV%Boussinesq)
if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701)

if (CS%answer_date >= 20190101) then
h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
elseif (GV%Boussinesq) then
h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10
else
h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10
endif

call initialize_remapping( CS%remapCS, 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)
answer_date=CS%answer_date, &
h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
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)
answer_date=CS%answer_date, &
h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)

call get_param(param_file, mdl, "PARTIAL_CELL_VELOCITY_REMAP", CS%partial_cell_vel_remap, &
"If true, use partial cell thicknesses at velocity points that are masked out "//&
Expand Down Expand Up @@ -345,7 +356,7 @@ 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 )
call remapping_set_param(CS%remapCS, om4_remap_via_sub_cells=om4_remap_via_sub_cells )

end subroutine ALE_set_OM4_remap_algorithm

Expand Down Expand Up @@ -591,8 +602,8 @@ subroutine ALE_offline_inputs(CS, G, GV, US, h, tv, Reg, uhtr, vhtr, Kd, debug,
endif
enddo ; enddo

call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%T, h_new, tv%T, answer_date=CS%answer_date)
call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%S, h_new, tv%S, answer_date=CS%answer_date)
call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%T, h_new, tv%T)
call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%S, h_new, tv%S)

if (debug) call MOM_tracer_chkinv("After ALE_offline_inputs", G, GV, h_new, Reg%Tr, Reg%ntr)

Expand Down Expand Up @@ -653,7 +664,6 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzInterface ! Interface height changes within
! an iteration [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzIntTotal ! Cumulative interface position changes [H ~> m or kg m-2]
real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2]

nz = GV%ke

Expand All @@ -680,14 +690,6 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d
if (present(dt)) &
call ALE_update_regrid_weights(dt, CS)

if (CS%answer_date >= 20190101) then
h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
elseif (GV%Boussinesq) then
h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10
else
h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10
endif

do itt = 1, n_itt

call do_group_pass(pass_T_S_h, G%domain)
Expand All @@ -704,10 +706,8 @@ subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, d

! remap from original grid onto new grid
do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), tv_local%S(i,j,:), &
h_neglect, h_neglect_edge)
call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), tv_local%T(i,j,:), &
h_neglect, h_neglect_edge)
call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), tv_local%S(i,j,:))
call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), tv_local%T(i,j,:))
enddo ; enddo

! starting grid for next iteration
Expand Down Expand Up @@ -763,22 +763,13 @@ subroutine ALE_remap_tracers(CS, G, GV, h_old, h_new, Reg, debug, dt, PCM_cell)
real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2]
real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2]
real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2]
logical :: show_call_tree
type(tracer_type), pointer :: Tr => NULL()
integer :: i, j, k, m, nz, ntr

show_call_tree = .false.
if (present(debug)) show_call_tree = debug

if (CS%answer_date >= 20190101) then
h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
elseif (GV%Boussinesq) then
h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10
else
h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10
endif

if (show_call_tree) call callTree_enter("ALE_remap_tracers(), MOM_ALE.F90")

nz = GV%ke
Expand All @@ -803,11 +794,9 @@ subroutine ALE_remap_tracers(CS, G, GV, h_old, h_new, Reg, debug, dt, PCM_cell)
h2(:) = h_new(i,j,:)
if (present(PCM_cell)) then
PCM(:) = PCM_cell(i,j,:)
call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column, &
h_neglect, h_neglect_edge, PCM_cell=PCM)
call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column, PCM_cell=PCM)
else
call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column, &
h_neglect, h_neglect_edge)
call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column)
endif

! Possibly underflow any very tiny tracer concentrations to 0. Note that this is not conservative!
Expand Down Expand Up @@ -1091,22 +1080,13 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
real :: v_tgt(GV%ke) ! A column of v-velocities on the target grid [L T-1 ~> m s-1]
real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2]
real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2]
real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2]
logical :: show_call_tree
integer :: i, j, k, nz

show_call_tree = .false.
if (present(debug)) show_call_tree = debug
if (show_call_tree) call callTree_enter("ALE_remap_velocities()")

if (CS%answer_date >= 20190101) then
h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
elseif (GV%Boussinesq) then
h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10
else
h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10
endif

nz = GV%ke

! --- Remap u profiles from the source vertical grid onto the new target grid.
Expand All @@ -1120,8 +1100,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
u_src(k) = u(I,j,k)
enddo

call remapping_core_h(CS%vel_remapCS, nz, h1, u_src, nz, h2, u_tgt, &
h_neglect, h_neglect_edge)
call remapping_core_h(CS%vel_remapCS, nz, h1, u_src, nz, h2, u_tgt)

if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) &
call mask_near_bottom_vel(u_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz)
Expand All @@ -1146,8 +1125,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
v_src(k) = v(i,J,k)
enddo

call remapping_core_h(CS%vel_remapCS, nz, h1, v_src, nz, h2, v_tgt, &
h_neglect, h_neglect_edge)
call remapping_core_h(CS%vel_remapCS, nz, h1, v_src, nz, h2, v_tgt)

if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) then
call mask_near_bottom_vel(v_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz)
Expand Down Expand Up @@ -1300,8 +1278,7 @@ end subroutine mask_near_bottom_vel
!> Remaps a single scalar between grids described by thicknesses h_src and h_dst.
!! h_dst must be dimensioned as a model array with GV%ke layers while h_src can
!! have an arbitrary number of layers specified by nk_src.
subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap, &
answers_2018, answer_date, h_neglect, h_neglect_edge)
subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap)
type(remapping_CS), intent(in) :: CS !< Remapping control structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
Expand All @@ -1319,44 +1296,16 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c
!! layers otherwise (default).
logical, optional, intent(in) :: old_remap !< If true, use the old "remapping_core_w"
!! method, otherwise use "remapping_core_h".
logical, optional, intent(in) :: answers_2018 !< If true, use the order of arithmetic
!! and expressions that recover the answers for
!! remapping from the end of 2018. Otherwise,
!! use more robust forms of the same expressions.
integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
!! for remapping
real, optional, intent(in) :: h_neglect !< A negligibly small thickness used in
!! remapping cell reconstructions, in the same
!! units as h_src, often [H ~> m or kg m-2]
real, optional, intent(in) :: h_neglect_edge !< A negligibly small thickness used in
!! remapping edge value calculations, in the same
!! units as h_src, often [H ~> m or kg m-2]
! Local variables
! Local variables
integer :: i, j, k, n_points
real :: dx(GV%ke+1) ! Change in interface position [H ~> m or kg m-2]
real :: h_neg, h_neg_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2]
logical :: ignore_vanished_layers, use_remapping_core_w, use_2018_remap
logical :: ignore_vanished_layers, use_remapping_core_w

ignore_vanished_layers = .false.
if (present(all_cells)) ignore_vanished_layers = .not. all_cells
use_remapping_core_w = .false.
if (present(old_remap)) use_remapping_core_w = old_remap
n_points = nk_src
use_2018_remap = .true. ; if (present(answers_2018)) use_2018_remap = answers_2018
if (present(answer_date)) use_2018_remap = (answer_date < 20190101)

if (present(h_neglect)) then
h_neg = h_neglect
h_neg_edge = h_neg ; if (present(h_neglect_edge)) h_neg_edge = h_neglect_edge
else
if (.not.use_2018_remap) then
h_neg = GV%H_subroundoff ; h_neg_edge = GV%H_subroundoff
elseif (GV%Boussinesq) then
h_neg = GV%m_to_H*1.0e-30 ; h_neg_edge = GV%m_to_H*1.0e-10
else
h_neg = GV%kg_m2_to_H*1.0e-30 ; h_neg_edge = GV%kg_m2_to_H*1.0e-10
endif
endif

!$OMP parallel do default(shared) firstprivate(n_points,dx)
do j = G%jsc,G%jec ; do i = G%isc,G%iec
Expand All @@ -1371,10 +1320,10 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c
if (use_remapping_core_w) then
call dzFromH1H2( n_points, h_src(i,j,1:n_points), GV%ke, h_dst(i,j,:), dx )
call remapping_core_w(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
GV%ke, dx, s_dst(i,j,:), h_neg, h_neg_edge)
GV%ke, dx, s_dst(i,j,:))
else
call remapping_core_h(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
GV%ke, h_dst(i,j,:), s_dst(i,j,:), h_neg, h_neg_edge)
GV%ke, h_dst(i,j,:), s_dst(i,j,:))
endif
else
s_dst(i,j,:) = 0.
Expand Down
Loading

0 comments on commit 3d1b6cd

Please sign in to comment.