diff --git a/components/mpas-ocean/cime_config/buildnml b/components/mpas-ocean/cime_config/buildnml
index 7dca58023e0d..e276df0a3f56 100755
--- a/components/mpas-ocean/cime_config/buildnml
+++ b/components/mpas-ocean/cime_config/buildnml
@@ -597,6 +597,8 @@ def buildnml(case, caseroot, compname):
             lines.append('    <var name="maxLevelEdgeBot"/>')
             lines.append('    <var name="columnIntegratedSpeed"/>')
             lines.append('    <var name="landIceFreshwaterFlux"/>')
+            lines.append('    <var name="landIceFreshwaterFluxTotal"/>')
+            lines.append('    <var name="frazilIceFreshwaterFlux"/>')
             lines.append('    <var name="landIceHeatFlux"/>')
             lines.append('    <var name="heatFluxToLandIce"/>')
             if ocn_bgc in ['eco_only', 'eco_and_dms', 'eco_and_macromolecules', 'eco_and_dms_and_macromolecules']:
@@ -873,7 +875,7 @@ def buildnml(case, caseroot, compname):
             lines.append('    <var name="tThreshMLD"/>')
             lines.append('    <var name="columnIntegratedSpeed"/>')
             lines.append('    <var name="barotropicSpeed"/>')
-            lines.append('    <var name="landIceFreshwaterFlux"/>')
+            lines.append('    <var name="landIceFreshwaterFluxTotal"/>')
             lines.append('    <var name="pressureAdjustedSSH"/>')
             lines.append('    <var name="atmosphericPressure"/>')
             if ocn_bgc in ['eco_only', 'eco_and_dms', 'eco_and_macromolecules', 'eco_and_dms_and_macromolecules']:
@@ -1297,6 +1299,8 @@ def buildnml(case, caseroot, compname):
             lines.append('    <var name="boundaryLayerDepth"/>')
             lines.append('    <var name="columnIntegratedSpeed"/>')
             lines.append('    <var name="landIceFreshwaterFlux"/>')
+            lines.append('    <var name="landIceFreshwaterFluxTotal"/>')
+            lines.append('    <var name="frazilIceFreshwaterFlux"/>')
             lines.append('    <var name="landIceHeatFlux"/>')
             lines.append('    <var name="heatFluxToLandIce"/>')
             lines.append('    <var name="landIceFrictionVelocity"/>')
@@ -1475,6 +1479,8 @@ def buildnml(case, caseroot, compname):
             lines.append('    <var name="snowFlux"/>')
             lines.append('    <var name="boundaryLayerDepth"/>')
             lines.append('    <var name="landIceFreshwaterFlux"/>')
+            lines.append('    <var name="landIceFreshwaterFluxTotal"/>')
+            lines.append('    <var name="frazilIceFreshwaterFlux"/>')
             lines.append('    <var name="landIceHeatFlux"/>')
             lines.append('    <var name="heatFluxToLandIce"/>')
             lines.append('    <var name="surfaceBuoyancyForcing"/>')
@@ -1537,6 +1543,8 @@ def buildnml(case, caseroot, compname):
             lines.append('    <var name="snowFlux"/>')
             lines.append('    <var name="boundaryLayerDepth"/>')
             lines.append('    <var name="landIceFreshwaterFlux"/>')
+            lines.append('    <var name="landIceFreshwaterFluxTotal"/>')
+            lines.append('    <var name="frazilIceFreshwaterFlux"/>')
             lines.append('    <var name="landIceHeatFlux"/>')
             lines.append('    <var name="heatFluxToLandIce"/>')
             lines.append('    <var name="surfaceBuoyancyForcing"/>')
diff --git a/components/mpas-ocean/driver/ocn_comp_mct.F b/components/mpas-ocean/driver/ocn_comp_mct.F
index 0097ac3f302a..be40de0622a9 100644
--- a/components/mpas-ocean/driver/ocn_comp_mct.F
+++ b/components/mpas-ocean/driver/ocn_comp_mct.F
@@ -753,13 +753,12 @@ end subroutine xml_stream_get_attributes
         call mpas_pool_get_subpool(block_ptr % structs, 'scratch', scratchPool)
 
         call ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, forcingPool, ierr, 1)
+        call ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, statePool, ierr)
         call mpas_timer_start("land_ice_build_arrays", .false.)
         call ocn_surface_land_ice_fluxes_build_arrays(meshPool, forcingPool, scratchPool, &
                                                       statePool, err=ierr)
         call mpas_timer_stop("land_ice_build_arrays")
 
-        call ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, statePool, ierr)
-
         block_ptr => block_ptr % next
     end do
 
@@ -1095,7 +1094,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{
             call mpas_pool_get_subpool(block_ptr % structs, 'scratch', scratchPool)
 
             call ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, forcingPool, ierr, 1)
-
+            call ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, statePool, ierr)
             call mpas_timer_start("land_ice_build_arrays", .false.)
             call ocn_surface_land_ice_fluxes_build_arrays(meshPool, &
                                                           forcingPool, scratchPool, statePool, ierr)
@@ -1103,8 +1102,6 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{
             call ocn_surface_land_ice_fluxes_accumulate_fluxes(meshPool, forcingPool, &
                                                                statePool, dt, ierr)
 
-            call ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, statePool, ierr)
-
             call ocn_eddy_compute_mixed_layer_depth(statePool, forcingPool)
             if (config_use_GM .or. config_submesoscale_enable) then
                 call ocn_eddy_compute_buoyancy_gradient()
diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index bfa8d3cc3f99..905343f3e408 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -3760,7 +3760,11 @@
 			/>
 		</var_array>
 		<var name="landIceFreshwaterFlux" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
-			description="Flux of mass through the ocean surface. Positive into ocean."
+			description="Flux of mass through the ocean surface due to interface melt/freeze. Positive into ocean."
+			packages="landIceFluxesPKG;dataLandIceFluxesPKG"
+		/>
+		<var name="landIceFreshwaterFluxTotal" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
+			description="Total flux of mass through the ocean surface that includes interface melt freeze and column averaged frazil freeze rate. Positive into ocean."
 			packages="landIceFluxesPKG;dataLandIceFluxesPKG"
 		/>
 		<var name="landIceHeatFlux"  type="real" dimensions="nCells Time" units="W m^-2"
@@ -3859,6 +3863,10 @@
 			 description="surface pressure forcing due to weight of frazil ice"
 			 packages="frazilIce"
 		/>
+		<var name="frazilIceFreshwaterFlux" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
+			description="Flux of frazil ice mass through the ocean surface. Negative is defined as mass leaving the ocean."
+			packages="frazilIce"
+		/>
 
 		<!-- Fields for tidal potential forcing -->
 		<var name="tidalPotentialEta" type="real" dimensions="nCells Time" units="m"
diff --git a/components/mpas-ocean/src/analysis_members/Registry_global_stats.xml b/components/mpas-ocean/src/analysis_members/Registry_global_stats.xml
index 1f76e701fd29..480ad08ef3f3 100644
--- a/components/mpas-ocean/src/analysis_members/Registry_global_stats.xml
+++ b/components/mpas-ocean/src/analysis_members/Registry_global_stats.xml
@@ -148,6 +148,9 @@
 			<var name="landIceFreshwaterFluxMin" array_group="mins" units="kg m^-2 s^-1"
 				 description="Minimum global value of landIceFreshwaterFlux in ocean cells."
 			/>
+			<var name="landIceFreshwaterFluxTotalMin" array_group="mins" units="kg m^-2 s^-1"
+				 description="Minimum global value of landIceFreshwaterFluxTotal in ocean cells."
+			/>
 			<var name="accumulatedLandIceMassMin" array_group="mins" units="kg m^-2"
 				 description="Minimum global value of accumulatedLandIceMass in ocean cells."
 			/>
@@ -157,6 +160,9 @@
 			<var name="accumulatedLandIceFrazilMassMin" array_group="mins" units="kg m^-2"
 				 description="Minimum global value of accumulatedLandIceFrazilMass in ocean cells."
 			/>
+			<var name="frazilIceFreshwaterFluxMin" array_group="mins" units="kg m^-2"
+				 description="Minimum global value of frazilIceFreshwaterFlux in ocean cells."
+			/>
 		</var_array>
 
 		<var_array name="maxGlobalStats" type="real" dimensions="Time">
@@ -247,6 +253,9 @@
 			<var name="landIceFreshwaterFluxMax" array_group="maxes" units="kg m^-2 s^-1"
 				 description="Maximum global value of landIceFreshwaterFlux in ocean cells."
 			/>
+			<var name="landIceFreshwaterFluxTotalMax" array_group="maxes" units="kg m^-2 s^-1"
+				 description="Maximum global value of landIceFreshwaterFluxTotal in ocean cells."
+			/>
 			<var name="accumulatedLandIceMassMax" array_group="maxes" units="kg m^-2"
 				 description="Maximum global value of accumulatedLandIceMass in ocean cells."
 			/>
@@ -256,6 +265,9 @@
 			<var name="accumulatedLandIceFrazilMassMax" array_group="maxes" units="kg m^-2"
 				 description="Maximum global value of accumulatedLandIceFrazilMass in ocean cells."
 			/>
+			<var name="frazilIceFreshwaterFluxMax" array_group="maxes" units="kg m^-2"
+				 description="Maximum global value of frazilIceFreshwaterFlux in ocean cells."
+			/>
 		</var_array>
 
 		<var_array name="sumGlobalStats" type="real" dimensions="Time">
@@ -346,6 +358,9 @@
 			<var name="landIceFreshwaterFluxSum" array_group="sums" units="kg s^-1"
 				 description="Accumulated global value of landIceFreshwaterFlux in ocean cells."
 			/>
+			<var name="landIceFreshwaterFluxTotalSum" array_group="sums" units="kg s^-1"
+				 description="Accumulated global value of landIceFreshwaterFluxTotal in ocean cells."
+			/>
 			<var name="accumulatedLandIceMassSum" array_group="sums" units="kg"
 				 description="Accumulated global value of accumulatedLandIceMass in ocean cells."
 			/>
@@ -355,6 +370,9 @@
 			<var name="accumulatedLandIceFrazilMassSum" array_group="sums" units="kg"
 				 description="Accumulated global value of accumulatedLandIceFrazilMass in ocean cells."
 			/>
+			<var name="frazilIceFreshwaterFluxSum" array_group="sums" units="kg"
+				 description="Accumulated global value of frazilIceFreshwaterFlux in ocean cells."
+			/>
 		</var_array>
 
 		<var_array name="rmsGlobalStats" type="real" dimensions="Time">
@@ -445,6 +463,9 @@
 			<var name="landIceFreshwaterFluxRms" array_group="rms" units="kg m^-2 s^-1"
 				 description="Global root mean square value of landIceFreshwaterFlux in ocean cells."
 			/>
+			<var name="landIceFreshwaterFluxTotalRms" array_group="rms" units="kg m^-2 s^-1"
+				 description="Global root mean square value of landIceFreshwaterFluxTotal in ocean cells."
+			/>
 			<var name="accumulatedLandIceMassRms" array_group="rms" units="kg m^-2"
 				 description="Global root mean square value of accumulatedLandIceMass in ocean cells."
 			/>
@@ -454,6 +475,9 @@
 			<var name="accumulatedLandIceFrazilMassRms" array_group="rms" units="kg m^-2"
 				 description="Global root mean square value of accumulatedLandIceFrazilMass in ocean cells."
 			/>
+			<var name="frazilIceFreshwaterFluxRms" array_group="rms" units="kg m^-2"
+				 description="Global root mean square value of frazilIceFreshwaterFlux in ocean cells."
+			/>
 		</var_array>
 
 		<var_array name="avgGlobalStats" type="real" dimensions="Time">
@@ -544,6 +568,9 @@
 			<var name="landIceFreshwaterFluxAvg" array_group="avg" units="kg m^-2 s^-1"
 				 description="Average value of landIceFreshwaterFlux in ocean cells."
 			/>
+			<var name="landIceFreshwaterFluxTotalAvg" array_group="avg" units="kg m^-2 s^-1"
+				 description="Average value of landIceFreshwaterFluxTotal in ocean cells."
+			/>
 			<var name="accumulatedLandIceMassAvg" array_group="avg" units="kg m^-2"
 				 description="Average value of accumulatedLandIceMass in ocean cells."
 			/>
@@ -553,6 +580,9 @@
 			<var name="accumulatedLandIceFrazilMassAvg" array_group="avg" units="kg m^-2"
 				 description="Average value of accumulatedLandIceFrazilMass in ocean cells."
 			/>
+			<var name="frazilIceFreshwaterFluxAvg" array_group="avg" units="kg m^-2"
+				 description="Average value of frazilIceFreshwaterFlux in ocean cells."
+			/>
 		</var_array>
 
 		<var_array name="vertSumMinGlobalStats" type="real" dimensions="Time">
@@ -643,6 +673,9 @@
 			<var name="landIceFreshwaterFluxMinVertSum" array_group="vertSumMin" units="kg m^-2 s^-1"
 				 description="Minimum vertical sum of landIceFreshwaterFlux in ocean cells."
 			/>
+			<var name="landIceFreshwaterFluxTotalMinVertSum" array_group="vertSumMin" units="kg m^-2 s^-1"
+				 description="Minimum vertical sum of landIceFreshwaterFluxTotal in ocean cells."
+			/>
 			<var name="accumulatedLandIceMassMinVertSum" array_group="vertSumMin" units="kg m^-2"
 				 description="Minimum vertical sum of accumulatedLandIceMass in ocean cells."
 			/>
@@ -652,6 +685,9 @@
 			<var name="accumulatedLandIceFrazilMassMinVertSum" array_group="vertSumMin" units="kg m^-2"
 				 description="Minimum vertical sum of accumulatedLandIceFrazilMass in ocean cells."
 			/>
+			<var name="frazilIceFreshwaterFluxMinVertSum" array_group="vertSumMin" units="kg m^-2"
+				 description="Minimum vertical sum of frazilIceFreshwaterFlux in ocean cells."
+			/>
 		</var_array>
 
 		<var_array name="vertSumMaxGlobalStats" type="real" dimensions="Time">
@@ -742,6 +778,9 @@
 			<var name="landIceFreshwaterFluxMaxVertSum" array_group="vertSumMax" units="kg m^-2 s^-1m"
 				 description="Maximum vertical sum of landIceFreshwaterFlux in ocean cells."
 			/>
+			<var name="landIceFreshwaterFluxTotalMaxVertSum" array_group="vertSumMax" units="kg m^-2 s^-1m"
+				 description="Maximum vertical sum of landIceFreshwaterFluxTotal in ocean cells."
+			/>
 			<var name="accumulatedLandIceMassMaxVertSum" array_group="vertSumMax" units="kg m^-2"
 				 description="Maximum vertical sum of accumulatedLandIceMass in ocean cells."
 			/>
@@ -751,6 +790,9 @@
 			<var name="accumulatedLandIceFrazilMassMaxVertSum" array_group="vertSumMax" units="kg m^-2"
 				 description="Maximum vertical sum of accumulatedLandIceFrazilMass in ocean cells."
 			/>
+			<var name="frazilIceFreshwaterFluxMaxVertSum" array_group="vertSumMax" units="kg m^-2"
+				 description="Maximum vertical sum of frazilIceFreshwaterFlux in ocean cells."
+			/>
 		</var_array>
 		<var name="totalVolumeChange" type="real" dimensions="Time" units="m"
 			 description="Total volume change of the ocean relative to previous timestep"
diff --git a/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_max.xml b/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_max.xml
index 3e79990bfd8f..c2a815e38e0b 100644
--- a/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_max.xml
+++ b/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_max.xml
@@ -128,6 +128,8 @@
 			<var name="snowFlux"/>
 			<var name="boundaryLayerDepth"/>
 			<var name="landIceFreshwaterFlux"/>
+			<var name="landIceFreshwaterFluxTotal"/>
+			<var name="frazilIceFreshwaterFlux"/>
 			<var name="landIceHeatFlux"/>
 			<var name="heatFluxToLandIce"/>
 			<var name="surfaceBuoyancyForcing"/>
diff --git a/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_mean.xml b/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_mean.xml
index 58ac7e95af4c..2858e01961d2 100644
--- a/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_mean.xml
+++ b/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_mean.xml
@@ -187,6 +187,8 @@
 			<var name="boundaryLayerDepth"/>
 			<var name="columnIntegratedSpeed"/>
 			<var name="landIceFreshwaterFlux"/>
+			<var name="landIceFreshwaterFluxTotal"/>
+			<var name="frazilIceFreshwaterFlux"/>
 			<var name="landIceHeatFlux"/>
 			<var name="heatFluxToLandIce"/>
 			<var name="landIceFrictionVelocity"/>
diff --git a/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_min.xml b/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_min.xml
index 5ee051585d71..30121568ad46 100644
--- a/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_min.xml
+++ b/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_min.xml
@@ -126,6 +126,8 @@
 			<var name="snowFlux"/>
 			<var name="boundaryLayerDepth"/>
 			<var name="landIceFreshwaterFlux"/>
+			<var name="landIceFreshwaterFluxTotal"/>
+			<var name="frazilIceFreshwaterFlux"/>
 			<var name="landIceHeatFlux"/>
 			<var name="heatFluxToLandIce"/>
 			<var name="surfaceBuoyancyForcing"/>
diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_global_stats.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_global_stats.F
index c21d1d4b5b50..203665999a86 100644
--- a/components/mpas-ocean/src/analysis_members/mpas_ocn_global_stats.F
+++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_global_stats.F
@@ -254,9 +254,9 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
 
       real (kind=RKIND), dimension(:), pointer :: evaporationFlux, snowFlux
       real (kind=RKIND), dimension(:), pointer :: seaIceFreshWaterFlux, icebergFreshWaterFlux, riverRunoffFlux, iceRunoffFlux
-      real (kind=RKIND), dimension(:), pointer :: rainFlux, landIceFreshwaterFlux
+      real (kind=RKIND), dimension(:), pointer :: rainFlux, landIceFreshwaterFlux, landIceFreshwaterFluxTotal
       real (kind=RKIND), dimension(:), pointer :: accumulatedLandIceMass, accumulatedLandIceHeat, &
-          accumulatedLandIceFrazilMass
+          accumulatedLandIceFrazilMass, frazilIceFreshwaterFlux
 
       real (kind=RKIND), dimension(:,:), pointer :: frazilLayerThicknessTendency
 
@@ -364,6 +364,8 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
          call mpas_pool_get_array(forcingPool, 'rainFlux', rainFlux)
          call mpas_pool_get_array(forcingPool, 'frazilLayerThicknessTendency', frazilLayerThicknessTendency)
          call mpas_pool_get_array(forcingPool, 'landIceFreshwaterFlux', landIceFreshwaterFlux)
+         call mpas_pool_get_array(forcingPool, 'landIceFreshwaterFluxTotal', landIceFreshwaterFluxTotal)
+         call mpas_pool_get_array(forcingPool, 'frazilIceFreshwaterFlux', frazilIceFreshwaterFlux)
          call mpas_pool_get_array(forcingPool, 'landIceFloatingMask', landIceFloatingMask)
 
          call mpas_pool_get_array(tracersSurfaceFluxPool, 'activeTracersSurfaceFlux', activeTracersSurfaceFlux)
@@ -377,6 +379,7 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
          call mpas_pool_get_array(statePool, 'accumulatedLandIceHeat', accumulatedLandIceHeat, 1)
          call mpas_pool_get_array(statePool, 'accumulatedLandIceFrazilMass', accumulatedLandIceFrazilMass, 1)
 
+
          allocate(areaEdge(1:nEdgesSolve))
          areaEdge = dcEdge(1:nEdgesSolve)*dvEdge(1:nEdgesSolve)
 
@@ -958,6 +961,21 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
             maxes_tmp(variableIndex) = 0.0_RKIND
          end if
 
+         ! landIceFreshwaterFluxTotal
+         variableIndex = variableIndex + 1
+         if ( associated(landIceFreshwaterFluxTotal) ) then
+            call ocn_compute_field_area_weighted_local_stats_surface(dminfo, nCellsSolve, &
+                         landIceFloatingArea, &
+                         landIceFreshwaterFluxTotal(1:nCellsSolve), sums_tmp(variableIndex), &
+                         sumSquares_tmp(variableIndex), mins_tmp(variableIndex), &
+                         maxes_tmp(variableIndex))
+         else
+            sums_tmp(variableIndex) = 0.0_RKIND
+            sumSquares_tmp(variableIndex) = 0.0_RKIND
+            mins_tmp(variableIndex) = 0.0_RKIND
+            maxes_tmp(variableIndex) = 0.0_RKIND
+         end if
+
          sums(variableIndex) = sums(variableIndex) + sums_tmp(variableIndex)
          sumSquares(variableIndex) = sumSquares(variableIndex) + sumSquares_tmp(variableIndex)
          mins(variableIndex) = min(mins(variableIndex), mins_tmp(variableIndex))
@@ -1024,6 +1042,21 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
             maxes_tmp(variableIndex) = 0.0_RKIND
          end if
 
+         ! frazilIceFreshwaterFlux
+         variableIndex = variableIndex + 1
+         if ( associated(frazilIceFreshwaterFlux) ) then
+            call ocn_compute_field_area_weighted_local_stats_surface(dminfo, nCellsSolve, &
+                         areaCell(1:nCellsSolve), &
+                         frazilIceFreshwaterFlux(1:nCellsSolve), sums_tmp(variableIndex), &
+                         sumSquares_tmp(variableIndex), mins_tmp(variableIndex), &
+                         maxes_tmp(variableIndex))
+         else
+            sums_tmp(variableIndex) = 0.0_RKIND
+            sumSquares_tmp(variableIndex) = 0.0_RKIND
+            mins_tmp(variableIndex) = 0.0_RKIND
+            maxes_tmp(variableIndex) = 0.0_RKIND
+         end if
+
          sums(variableIndex) = sums(variableIndex) + sums_tmp(variableIndex)
          sumSquares(variableIndex) = sumSquares(variableIndex) + sumSquares_tmp(variableIndex)
          mins(variableIndex) = min(mins(variableIndex), mins_tmp(variableIndex))
@@ -1299,6 +1332,16 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
          rms(variableIndex) = 0.0_RKIND
       end if
 
+      ! landIceFreshwaterFluxTotal
+      variableIndex = variableIndex + 1
+      if (associated(landIceFreshwaterFluxTotal) .and. landIceFloatingAreaSum > 0.0_RKIND) then
+         averages(variableIndex) = sums(variableIndex)/(landIceFloatingAreaSum)
+         rms(variableIndex) = sqrt(sumSquares(variableIndex)/(landIceFloatingAreaSum))
+      else
+         averages(variableIndex) = 0.0_RKIND
+         rms(variableIndex) = 0.0_RKIND
+      end if
+
       ! continue accumulating fresh water inputs
       netFreshwaterInput = netFreshwaterInput + sums(variableIndex) * dt/rho_sw
 
@@ -1332,6 +1375,16 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
          rms(variableIndex) = 0.0_RKIND
       end if
 
+      ! frazilIceFreshwaterFlux
+      variableIndex = variableIndex + 1
+      if (frazilIcePkgActive) then
+         averages(variableIndex) = sums(variableIndex)/(areaCellGlobal)
+         rms(variableIndex) = sqrt(sumSquares(variableIndex)/(areaCellGlobal))
+      else
+         averages(variableIndex) = 0.0_RKIND
+         rms(variableIndex) = 0.0_RKIND
+      end if
+
       ! calculate fresh water conservation check quantities
       absoluteFreshWaterConservation = totalVolumeChange - netFreshwaterInput
       if (abs(totalVolumeChange) < 1e-12_RKIND) then
diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F
index af69a63e46cf..daf0417a2091 100644
--- a/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F
+++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F
@@ -700,6 +700,7 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
            call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool)
            call mpas_pool_get_subpool(block_ptr % structs, 'scratch', scratchPool)
            call ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, forcingPool, ierr, 1)
+           call ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, statePool, err)
            call mpas_timer_start("land_ice_build_arrays")
            call ocn_surface_land_ice_fluxes_build_arrays(meshPool, &
                                                          forcingPool, scratchPool, statePool, err)
@@ -707,7 +708,6 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
            call ocn_surface_land_ice_fluxes_accumulate_fluxes(meshPool, forcingPool, &
                                                               statePool, dt, err)
 
-           call ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, statePool, err)
            call ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, statePool, err)
 
            ! Compute normalGMBolusVelocity, relativeSlope and RediDiffVertCoef if respective flags are turned on
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F b/components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F
index 16c223986840..0f06a8d39c16 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F
@@ -380,6 +380,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceSalinityOld
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedLandIceFrazilMassNew
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedLandIceFrazilMassOld
+      real (kind=RKIND), pointer, dimension(:)     :: frazilIceFreshwaterFlux
       real (kind=RKIND), pointer, dimension(:)     :: ssh
       real (kind=RKIND), pointer, dimension(:,:)   :: layerThickness
       real (kind=RKIND), pointer, dimension(:,:,:) :: activeTracers
@@ -436,6 +437,7 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       call mpas_pool_get_array(forcingPool, 'frazilSurfacePressure', frazilSurfacePressure)
       call mpas_pool_get_array(forcingPool, 'landIceMask', landIceMask)
       call mpas_pool_get_array(forcingPool, 'landIceFloatingMask', landIceFloatingMask)
+      call mpas_pool_get_array(forcingPool, 'frazilIceFreshwaterFlux', frazilIceFreshwaterFlux)
 
       ! get time step in units of seconds
       timeStep = mpas_get_clock_timestep(domain % clock, ierr=err)
@@ -607,6 +609,8 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
 
         enddo   ! do k=kBottom,minLevelCell,-1
 
+        ! frazilIceFreshwaterFlux forms part of the landIceFreshwaterFluxTotal computed in surface land ice fluxes
+        frazilIceFreshwaterFlux(iCell) = - (sumNewFrazilIceThickness * config_frazil_ice_density) / dt
         ! accumulate frazil mass to column total
         ! note: the accumulatedFrazilIceMass (at both time levels) is reset to zero after being sent to the coupler
         accumulatedFrazilIceMassNew(iCell) = accumulatedFrazilIceMassOld(iCell) + sumNewFrazilIceThickness &
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F b/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F
index 69c60c79ca51..c249f202757d 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F
@@ -468,6 +468,8 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, &
                                                   landIceFreshwaterFlux, &
                                                   landIceHeatFlux, heatFluxToLandIce
 
+      real (kind=RKIND), dimension(:), pointer :: frazilIceFreshwaterFlux, landIceFreshwaterFluxTotal
+
       integer, dimension(:), pointer :: landIceFloatingMask
 
       real (kind=RKIND), dimension(:,:), pointer :: &
@@ -505,6 +507,9 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, &
       call mpas_pool_get_array(forcingPool, 'landIceHeatFlux', landIceHeatFlux)
       call mpas_pool_get_array(forcingPool, 'heatFluxToLandIce', heatFluxToLandIce)
 
+      call mpas_pool_get_array(forcingPool, 'frazilIceFreshwaterFlux', frazilIceFreshwaterFlux)
+      call mpas_pool_get_array(forcingPool, 'landIceFreshwaterFluxTotal', landIceFreshwaterFluxTotal)
+
       call mpas_pool_get_array(forcingPool, 'landIceInterfaceTracers', landIceInterfaceTracers)
       call mpas_pool_get_dimension(forcingPool, &
                                 'index_landIceInterfaceTemperature', &
@@ -650,6 +655,17 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, &
 
       endif ! jenkinsOn or hollandJenkinsOn
 
+      ! Add frazil and interface melt/freeze to get total fresh water flux
+      if ( associated(frazilIceFreshwaterFlux) ) then
+          do iCell = 1, nCells
+                landIceFreshwaterFluxTotal(iCell) = landIceFreshwaterFlux(iCell) + landIceFloatingMask(iCell)*frazilIceFreshwaterFlux(iCell)
+          end do
+      else
+          do iCell = 1, nCells
+                landIceFreshwaterFluxTotal(iCell) = landIceFreshwaterFlux(iCell)
+          end do
+      end if
+
       if(useHollandJenkinsAdvDiff) then
          deallocate(freezeInterfaceSalinity, &
                     freezeInterfaceTemperature, &