diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90
index 324808e374..8e988ccce8 100644
--- a/src/framework/MOM_horizontal_regridding.F90
+++ b/src/framework/MOM_horizontal_regridding.F90
@@ -591,7 +591,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr
 
     ! Horizontally homogenize data to produce perfectly "flat" initial conditions
     if (PRESENT(homogenize)) then ; if (homogenize) then
-      call homogenize_field(tr_out, mask_out, G, scale, answer_date)
+      call homogenize_field(tr_out, G, tmp_scale=I_scale, weights=mask_out, answer_date=answer_date)
     endif ; endif
 
     ! tr_out contains input z-space data on the model grid with missing values
@@ -908,7 +908,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, &
 
       ! Horizontally homogenize data to produce perfectly "flat" initial conditions
       if (PRESENT(homogenize)) then ; if (homogenize) then
-        call homogenize_field(tr_out, mask_out, G, scale, answer_date)
+        call homogenize_field(tr_out, G, tmp_scale=I_scale, weights=mask_out, answer_date=answer_date)
       endif ; endif
 
       ! tr_out contains input z-space data on the model grid with missing values
@@ -950,14 +950,15 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(field, Time, G, tr_z, mask_z, &
 end subroutine horiz_interp_and_extrap_tracer_fms_id
 
 !> Replace all values of a 2-d field with the weighted average over the valid points.
-subroutine homogenize_field(field, weight, G, scale, answer_date, wt_unscale)
+subroutine homogenize_field(field, G, tmp_scale, weights, answer_date, wt_unscale)
   type(ocean_grid_type),            intent(inout) :: G      !< Ocean grid type
   real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: field  !< The tracer on the model grid in arbitrary units [A ~> a]
-  real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: weight !< The weights for the tracer in arbitrary units that
+  real,                   optional, intent(in)    :: tmp_scale !< A temporary rescaling factor for the
+                                                            !! variable that is reversed in the
+                                                            !! return value [a A-1 ~> 1]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                          optional, intent(in)    :: weights !< The weights for the tracer in arbitrary units that
                                                             !! typically differ from those used by field [B ~> b]
-  real,                             intent(in)    :: scale  !< A rescaling factor that has been used for the
-                                                            !! variable and has to be undone before the
-                                                            !! reproducing sums [A a-1 ~> 1]
   integer,                optional, intent(in)    :: answer_date !< The vintage of the expressions in the code.
                                                             !! Dates before 20230101 use non-reproducing sums
                                                             !! in their averages, while later versions use
@@ -971,12 +972,11 @@ subroutine homogenize_field(field, weight, G, scale, answer_date, wt_unscale)
   ! In the following comments, [A] and [B] are used to indicate the arbitrary, possibly rescaled
   ! units of the input field and the weighting array, while [a] and [b] indicate the corresponding
   ! unscaled (e.g., mks) units that can be used with the reproducing sums
-  real, dimension(SZI_(G),SZJ_(G)) :: field_for_Sums  ! The field times the weights with the scaling undone [a b]
-  real, dimension(SZI_(G),SZJ_(G)) :: wts_for_Sums    ! A copy of the wieghts with the scaling undone [b]
+  real, dimension(G%isc:G%iec, G%jsc:G%jec) :: field_for_Sums  ! The field times the weights [A B ~> a b]
+  real, dimension(G%isc:G%iec, G%jsc:G%jec) :: weight ! A copy of weights, if it is present, or the
+                      ! tracer-point grid mask if it weights is absent [B ~> b]
   real :: var_unscale ! The reciprocal of the scaling factor for the field and weights [a b A-1 B-1 ~> 1]
-  real :: wt_descale  ! A factor that undoes any dimensional scaling of the weights so that they
-                      ! can be used with reproducing sums [b B-1 ~> 1]
-  real :: wt_sum      ! The sum of the weights, in [b] (reproducing) or [B ~> b] (non-reproducing)
+  real :: wt_sum      ! The sum of the weights, in [B ~> b]
   real :: varsum      ! The weighted sum of field being averaged [A B ~> a b]
   real :: varAvg      ! The average of the field [A ~> a]
   logical :: use_repro_sums  ! If true, use reproducing sums.
@@ -988,23 +988,27 @@ subroutine homogenize_field(field, weight, G, scale, answer_date, wt_unscale)
 
   use_repro_sums = .false. ; if (present(answer_date)) use_repro_sums = (answer_date >= 20230101)
 
-  if (scale == 0.0) then
-    ! This seems like an unlikely case to ever be used, but dealing with it is better than having NaNs arise?
-    varAvg = 0.0
-  elseif (use_repro_sums) then
-    wt_descale = 1.0 ; if (present(wt_unscale)) wt_descale = wt_unscale
-    var_unscale = wt_descale / scale
+  if (present(weights)) then
+    do j=js,je ; do i=is,ie
+      weight(i,j) = weights(i,j)
+    enddo ; enddo
+  else
+    do j=js,je ; do i=is,ie
+      weight(i,j) = G%mask2dT(i,j)
+    enddo ; enddo
+  endif
+
+  if (use_repro_sums) then
+    var_unscale = 1.0 ; if (present(tmp_scale)) var_unscale = tmp_scale
+    if (present(wt_unscale)) var_unscale = wt_unscale * var_unscale
 
-    field_for_Sums(:,:) = 0.0
-    wts_for_Sums(:,:) = 0.0
     do j=js,je ; do i=is,ie
-      wts_for_Sums(i,j) = wt_descale * weight(i,j)
-      field_for_Sums(i,j) = var_unscale * (field(i,j) * weight(i,j))
+      field_for_Sums(i,j) = field(i,j) * weight(i,j)
     enddo ; enddo
 
-    wt_sum = reproducing_sum(wts_for_Sums)
+    wt_sum = reproducing_sum(weight, unscale=wt_unscale)
     if (abs(wt_sum) > 0.0) &
-      varAvg = reproducing_sum(field_for_Sums) * (scale / wt_sum)
+      varAvg = reproducing_sum(field_for_Sums, unscale=var_unscale) * (1.0 / wt_sum)
 
   else  ! Do the averages with order-dependent sums to reproduce older answers.
     wt_sum = 0 ; varsum = 0.
@@ -1021,8 +1025,12 @@ subroutine homogenize_field(field, weight, G, scale, answer_date, wt_unscale)
       call sum_across_PEs(varsum)
       varAvg = varsum / wt_sum
     endif
+
   endif
 
+  ! This seems like an unlikely case to ever be used, but it is needed to recreate previous behavior.
+  if (present(tmp_scale)) then ; if (tmp_scale == 0.0) varAvg = 0.0 ; endif
+
   field(:,:) = varAvg
 
 end subroutine homogenize_field
diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90
index 41b407d6a1..d908ec23a0 100644
--- a/src/initialization/MOM_state_initialization.F90
+++ b/src/initialization/MOM_state_initialization.F90
@@ -2928,8 +2928,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
     if (homogenize) then
       ! Horizontally homogenize data to produce perfectly "flat" initial conditions
       do k=1,nz
-        call homogenize_field(tv%T(:,:,k), G%mask2dT, G, scale=US%degC_to_C, answer_date=hor_regrid_answer_date)
-        call homogenize_field(tv%S(:,:,k), G%mask2dT, G, scale=US%ppt_to_S, answer_date=hor_regrid_answer_date)
+        call homogenize_field(tv%T(:,:,k), G, tmp_scale=US%C_to_degC, answer_date=hor_regrid_answer_date)
+        call homogenize_field(tv%S(:,:,k), G, tmp_scale=US%S_to_ppt, answer_date=hor_regrid_answer_date)
       enddo
     endif