+(*)Refactor MOM_diag_remap for global indices

  Refactored MOM_diag_remap to work with global indices and to move logical
branches outside of do loops.  This was done by adding internal routines that
set the loop indices consistently with the NE c-grid convention used throughout
the MOM6 code and converting the optionally associated mask pointer into an
optional argument.  It also simplifies the logic of many of the expressions
within the remapping code.  There is also a new element, Z_based_coord, in the
diag_remap_ctrl type, to indicate whether the remapping is working in thickness
or height units, but for now it is always set to false.

  The function set_h_neglect or set_dz_neglect is used to set the negligible
thicknesses used for remapping in diag_remap_update and diag_remap_do_remap,
depending on whether the remapping is being done in thickness or vertical height
coordinates.  Diag_remap_init has a new vertical_grid_type argument, and
diag_remap_do_remap has a new unit_scale_type argument.

  For REMAPPING_ANSWER_DATES later than 20240201, diag_remap_updated does an
explicit sum to determine the total water column thickness, rather than using
sum function, which is indeterminate of the order of the sums.  For some
compilers, this could change the vertical grids used for remapping diagnostics
at roundoff, but no such change was detected for any of the compilers used with
the MOM6 regression test suite.

  All answers and diagnostics in cases that worked before are bitwise identical,
but there are new arguments to two publicly visible interfaces.

src/framework/MOM_diag_mediator.F90

@@ -1583,7 +1583,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h)
       if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)
       allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz))
       call diag_remap_do_remap(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), &
-              diag_cs%G, diag_cs%GV, h_diag, staggered_in_x, staggered_in_y, &
+              diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, staggered_in_x, staggered_in_y, &
               diag%axes%mask3d, field, remapped_field)
       if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)
       if (associated(diag%axes%mask3d)) then
@@ -3207,7 +3207,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)
+      call diag_remap_init(diag_cs%diag_remap_cs(i), diag_coords(i), answer_date=remap_answer_date, GV=GV)
     enddo
     deallocate(diag_coords)
   endif