diff --git a/components/elm/src/external_models/fates b/components/elm/src/external_models/fates
index e027a4030d2a..83863e9398e2 160000
--- a/components/elm/src/external_models/fates
+++ b/components/elm/src/external_models/fates
@@ -1 +1 @@
-Subproject commit e027a4030d2a0f09039fb337ad67ced7461dd4f0
+Subproject commit 83863e9398e2098b5ef04a8b8525733a881c1ca3
diff --git a/components/mpas-ocean/driver/mpaso_cpl_indices.F b/components/mpas-ocean/driver/mpaso_cpl_indices.F
index de86b439c2ea..e7daa4f9b57b 100644
--- a/components/mpas-ocean/driver/mpaso_cpl_indices.F
+++ b/components/mpas-ocean/driver/mpaso_cpl_indices.F
@@ -21,6 +21,8 @@ module mpaso_cpl_indices
   integer :: index_o2x_Fioo_q
   integer :: index_o2x_Foxo_q_li
   integer :: index_o2x_Fioo_frazil
+  integer :: index_o2x_Fioo_frazils
+  integer :: index_o2x_Fioo_frazilh
   integer :: index_o2x_Foxo_frazil_li
   integer :: index_o2x_Faoo_h2otemp
   integer :: index_o2x_Faoo_fco2_ocn
@@ -215,6 +217,8 @@ subroutine mpaso_cpl_indices_set( )
     index_o2x_Fioo_q        = mct_avect_indexra(o2x,'Fioo_q',perrWith='quiet')
     index_o2x_Foxo_q_li     = mct_avect_indexra(o2x,'Foxo_q_li',perrWith='quiet')
     index_o2x_Fioo_frazil   = mct_avect_indexra(o2x,'Fioo_frazil',perrWith='quiet')
+    index_o2x_Fioo_frazils  = mct_avect_indexra(o2x,'Fioo_frazils',perrWith='quiet')
+    index_o2x_Fioo_frazilh  = mct_avect_indexra(o2x,'Fioo_frazilh',perrWith='quiet')
     index_o2x_Foxo_frazil_li= mct_avect_indexra(o2x,'Foxo_frazil_li',perrWith='quiet')
     index_o2x_Faoo_h2otemp  = mct_avect_indexra(o2x,'Faoo_h2otemp',perrWith='quiet')
     index_o2x_Faoo_fco2_ocn = mct_avect_indexra(o2x,'Faoo_fco2_ocn',perrWith='quiet')
diff --git a/components/mpas-ocean/driver/ocn_comp_mct.F b/components/mpas-ocean/driver/ocn_comp_mct.F
index fc3c4edda820..9de9c00e0f7b 100644
--- a/components/mpas-ocean/driver/ocn_comp_mct.F
+++ b/components/mpas-ocean/driver/ocn_comp_mct.F
@@ -3166,6 +3166,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
    integer, dimension(:), pointer :: landIceMask
 
    real (kind=RKIND), dimension(:), pointer :: seaIceEnergy, accumulatedFrazilIceMass, frazilSurfacePressure, &
+                                               accumulatedFrazilIceEnergy, &
                                                avgTotalFreshWaterTemperatureFlux, &
                                                avgCO2_gas_flux, DMSFlux, surfaceUpwardCO2Flux, &
                                                avgOceanSurfaceDIC, &
@@ -3286,6 +3287,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
         call mpas_pool_get_array(forcingPool, 'seaIceEnergy', seaIceEnergy)
         call mpas_pool_get_array(forcingPool, 'frazilSurfacePressure', frazilSurfacePressure)
         call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMass, 1)
+        call mpas_pool_get_array(statePool, 'accumulatedFrazilIceEnergy', accumulatedFrazilIceEnergy, 1)
      end if
 
      ! Cryo fields
@@ -3405,11 +3407,17 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
 
              o2x_o(index_o2x_Fioo_q, n)  = seaIceEnergy(i) / ocn_cpl_dt
              o2x_o(index_o2x_Fioo_frazil, n) = accumulatedFrazilIceMass(i) / ocn_cpl_dt
+             o2x_o(index_o2x_Fioo_frazils, n) = accumulatedFrazilIceMass(i)/ 1000.0_RKIND * &
+                config_frazil_sea_ice_reference_salinity / ocn_cpl_dt
+             !defined positive to be consistent with Fioo_q when frazil is formed
+             o2x_o(index_o2x_Fioo_frazilh, n) = accumulatedFrazilIceEnergy(i)/ ocn_cpl_dt 
 
           else
 
              o2x_o(index_o2x_Fioo_q, n)  = 0.0_RKIND
              o2x_o(index_o2x_Fioo_frazil, n) = 0.0_RKIND
+             o2x_o(index_o2x_Fioo_frazils, n) = 0.0_RKIND
+             o2x_o(index_o2x_Fioo_frazilh, n) = 0.0_RKIND
              if (trim(config_land_ice_flux_mode) == 'standalone' .or. trim(config_land_ice_flux_mode) == 'data') then
                 o2x_o(index_o2x_Foxo_q_li, n) = accumulatedFrazilIceMass(i) * config_frazil_heat_of_fusion / ocn_cpl_dt
                 o2x_o(index_o2x_Foxo_frazil_li, n) = accumulatedFrazilIceMass(i) / ocn_cpl_dt
@@ -3420,6 +3428,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
           ! Reset SeaIce Energy and Accumulated Frazil Ice
           seaIceEnergy(i) = 0.0_RKIND
           accumulatedFrazilIceMass(i) = 0.0_RKIND
+          accumulatedFrazilIceEnergy(i) = 0.0_RKIND
           frazilSurfacePressure(i) = 0.0_RKIND
        end if
 
diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index 00f15a2961f9..31f950349fa8 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -2146,6 +2146,7 @@
 			<var_struct name="ecosysAuxiliary"/>
 			<var name="accumulatedFrazilIceMass"/>
 			<var name="accumulatedFrazilIceSalinity"/>
+			<var name="accumulatedFrazilIceEnergy"/>
 			<var name="accumulatedLandIceMass"/>
 			<var name="accumulatedLandIceHeat"/>
 			<var name="accumulatedLandIceFrazilMass"/>
@@ -2517,6 +2518,10 @@
 			 description="Mass per unit area of frazil ice produced. Reset to zero at each coupling interval"
 			 packages="frazilIce"
 		/>
+		<var name="accumulatedFrazilIceEnergy" type="real" dimensions="nCells Time" units="J m^-2"
+			 description="Accumulated frazil ice energy per unit area. Reset to zero at each coupling interval"
+			 packages="frazilIce"
+		/>
 		<var name="accumulatedFrazilIceSalinity" type="real" dimensions="nCells Time" units="kg m^-2"
 			 description="Salinity associated with accumulatedFrazilIceMass. Reset to zero at each coupling interval"
 			 packages="frazilIce"
diff --git a/components/mpas-ocean/src/analysis_members/Registry_conservation_check.xml b/components/mpas-ocean/src/analysis_members/Registry_conservation_check.xml
index f0096d237593..9795b52379f4 100644
--- a/components/mpas-ocean/src/analysis_members/Registry_conservation_check.xml
+++ b/components/mpas-ocean/src/analysis_members/Registry_conservation_check.xml
@@ -86,7 +86,7 @@
 			description="Iceberg heat flux from coupler, integrated in space and time. Positive into the ocean."
 		/>
 		<var name="accumulatedFrazilHeatFlux" type="real" dimensions="Time" units="W m^-2"
-			description="Heat flux from frazil, integrated in space and time, given to couipler. Positive into the ocean."
+			description="Heat flux from frazil, integrated in space and time, given to coupler. Positive into the ocean."
 		/>
 		<var name="accumulatedLandIceHeatFlux" type="real" dimensions="Time" units="W m^-2"
 			description="Land ice heat flux, integrated in space and time. Positive into the ocean."
diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F
index 95a9bca8f70c..8b881753d0fb 100644
--- a/components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F
+++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F
@@ -462,6 +462,8 @@ subroutine energy_conservation(domain, err)
            icebergTemperatureFlux, &
            accumulatedFrazilIceMassNew, &
            accumulatedFrazilIceMassOld, &
+           accumulatedFrazilIceEnergyNew, &
+           accumulatedFrazilIceEnergyOld, &
            accumulatedLandIceFrazilMassNew, &
            accumulatedLandIceFrazilMassOld
 
@@ -565,6 +567,8 @@ subroutine energy_conservation(domain, err)
             call mpas_pool_get_array(tracersSurfaceFluxPool, 'activeTracersSurfaceFluxSubglacialRunoff', activeTracersSurfaceFluxSubglacialRunoff)
             call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMassNew, 2)
             call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMassOld, 1)
+            call mpas_pool_get_array(statePool, 'accumulatedFrazilIceEnergy', accumulatedFrazilIceEnergyNew, 2)
+            call mpas_pool_get_array(statePool, 'accumulatedFrazilIceEnergy', accumulatedFrazilIceEnergyOld, 1)
 
             do iCell = 1, nCellsSolve
 
@@ -599,8 +603,10 @@ subroutine energy_conservation(domain, err)
                do iCell = 1, nCellsSolve
                   ! Frazil ice mass is negative. Negative coefficient makes heat
                   ! flux positive, because freezing ice releases heat.
-                  sumArray(10) = sumArray(10) - areaCell(iCell) * config_frazil_heat_of_fusion &
-                                 * (accumulatedFrazilIceMassNew(iCell) - accumulatedFrazilIceMassOld(iCell))/dt
+                  !sumArray(10) = sumArray(10) - areaCell(iCell) * config_frazil_heat_of_fusion &
+                  !               * (accumulatedFrazilIceMassNew(iCell) - accumulatedFrazilIceMassOld(iCell))/dt
+                  sumArray(10) = sumArray(10) - areaCell(iCell) &
+                                  * (accumulatedFrazilIceEnergyNew(iCell) - accumulatedFrazilIceEnergyOld(iCell))/dt                        
                enddo
             end if
 
@@ -749,7 +755,8 @@ subroutine energy_conservation(domain, err)
             s = 0.0_RKIND
             call mpas_log_write('HEAT FLUXES: explicit')
             call mpas_log_write('MPAS-Ocean name            W (raw sum)    coupler name    short name         W/m^2 (flux/A)')
-v=accumulatedFrazilHeatFlux          ; write(m,"('frazilMass * LHF /dt     ',es16.8,' o2x_Fioo_q      hfreeze         ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v
+!v=accumulatedFrazilHeatFlux          ; write(m,"('frazilMass * LHF /dt     ',es16.8,' o2x_Fioo_q      hfreeze         ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v
+v=accumulatedFrazilHeatFlux          ; write(m,"('frazilMass * LHF /dt     ',es16.8,' o2x_Fioo_frazilh hfreeze         ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v
 v=accumulatedSeaIceHeatFlux          ; write(m,"('seaIceHeatFlux           ',es16.8,' x2o_Fioi_melth  hmelt           ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v
 v=accumulatedShortWaveHeatFlux       ; write(m,"('shortWaveHeatFlux        ',es16.8,' x2o_Foxx_swnet  hnetsw          ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v
 v=accumulatedLongWaveHeatFluxDown    ; write(m,"('longWaveHeatFluxDown     ',es16.8,' x2o_Faxa_lwdn   hlwdn           ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v
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 fb92c4885404..d8cb7042b597 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_frazil_forcing.F
@@ -462,10 +462,13 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       real (kind=RKIND) :: newFrazilIceThickness, newThicknessWeightedSaltContent
       real (kind=RKIND) :: sumNewFrazilIceThickness, sumNewThicknessWeightedSaltContent
       real (kind=RKIND) :: meltedFrazilIceThickness, meltedThicknessWeightedSaltContent
+      real (kind=RKIND) :: sumNewFrazilEnergy, newFrazilEnergy, meltedFrazilEnergy
       real (kind=RKIND) :: oceanFreezingTemperature, frazilSalinity
 
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceMassNew
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceMassOld
+      real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceEnergyNew
+      real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceEnergyOld
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceSalinityNew
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedFrazilIceSalinityOld
       real (kind=RKIND), pointer, dimension(:)     :: accumulatedLandIceFrazilMassNew
@@ -515,6 +518,8 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1)
       call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMassNew, 2)
       call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMassOld, 1)
+      call mpas_pool_get_array(statePool, 'accumulatedFrazilIceEnergy', accumulatedFrazilIceEnergyNew, 2)
+      call mpas_pool_get_array(statePool, 'accumulatedFrazilIceEnergy', accumulatedFrazilIceEnergyOld, 1)
       call mpas_pool_get_array(statePool, 'accumulatedFrazilIceSalinity', accumulatedFrazilIceSalinityNew, 2)
       call mpas_pool_get_array(statePool, 'accumulatedFrazilIceSalinity', accumulatedFrazilIceSalinityOld, 1)
       call mpas_pool_get_array(statePool, 'accumulatedLandIceFrazilMass', accumulatedLandIceFrazilMassNew, 2)
@@ -552,7 +557,8 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
       !$omp         sumNewThicknessWeightedSaltContent, oceanFreezingTemperature, potential, &
       !$omp         freezingEnergy, meltingEnergy, frazilSalinity, newFrazilIceThickness, &
       !$omp         newThicknessWeightedSaltContent, meltedFrazilIceThickness, &
-      !$omp         meltedThicknessWeightedSaltContent)
+      !$omp         meltedThicknessWeightedSaltContent, newFrazilEnergy, meltedFrazilEnergy, &
+      !$omp         sumNewFrazilEnergy, newFrazilEnergy)
       do iCell=1,nCells
 
         underLandIce = .false.
@@ -604,6 +610,8 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
         ! initialize the sum of new frazil ice created
         sumNewFrazilIceThickness = 0.0_RKIND
         sumNewThicknessWeightedSaltContent = 0.0_RKIND
+        sumNewFrazilEnergy = 0.0_RKIND
+        newFrazilEnergy = 0.0_RKIND
 
         ! loop from maximum depth of frazil creation to surface
         do k = kBottomFrazil, minLevelCell(iCell), -1
@@ -632,31 +640,29 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
 
             ! new frazil ice formation measured in meters
             newFrazilIceThickness = freezingEnergy / (config_frazil_heat_of_fusion * config_frazil_ice_density)
-
             ! limit the frazil formed appropriately
             newFrazilIceThickness = min(newFrazilIceThickness, layerThickness(k,iCell) * config_frazil_fractional_thickness_limit)
 
             ! Determine new salt in frazil
             newThicknessWeightedSaltContent = newFrazilIceThickness * frazilSalinity
-
             ! compute tendency to thickness, temperature and salinity
             ! layerTendency is scaled so that mass of ice created == mass of ocean water removed
 
             ! layer thickness decreased due to creation of frazil
             ! note: -- this has to be density (not rho_sw) to keep buoyancy equal
             frazilLayerThicknessTendency(k,iCell) = - newFrazilIceThickness * config_frazil_ice_density / density(k,iCell) / dt
-
             ! salt is extracted with the frazil
             frazilSalinityTendency(k,iCell) = - newThicknessWeightedSaltContent / dt
 
             ! ocean fluid temperature is warmed due to creation of frazil
             frazilTemperatureTendency(k,iCell) =  + ( newFrazilIceThickness * config_frazil_heat_of_fusion &
                                                * config_frazil_ice_density ) / (config_specific_heat_sea_water * rho_sw) / dt
-
-            ! keep track of sum of frazil ice
+            ! Store the latent energy of the new frazil 
+            newFrazilEnergy = newFrazilIceThickness * config_frazil_ice_density * config_frazil_heat_of_fusion
+            !keep track of sum of frazil ice
             sumNewFrazilIceThickness = sumNewFrazilIceThickness + newFrazilIceThickness
             sumNewThicknessWeightedSaltContent = sumNewThicknessWeightedSaltContent + newThicknessWeightedSaltContent
-
+            sumNewFrazilEnergy = sumNewFrazilEnergy + newFrazilEnergy
           else
             ! ocean water is warm enough to melt frazil
 
@@ -688,18 +694,24 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
               frazilSalinityTendency(k,iCell) = + meltedThicknessWeightedSaltContent/ dt
 
               ! ocean fluid temperature is cooled due to melting of frazil
-              frazilTemperatureTendency(k,iCell) = - ( meltedFrazilIceThickness * config_frazil_heat_of_fusion &
+              frazilTemperatureTendency(k,iCell) =  - ( meltedFrazilIceThickness * config_frazil_heat_of_fusion &
                                                  * config_frazil_ice_density ) / (config_specific_heat_sea_water * rho_sw) &
                                                  / dt
 
               ! ocean fluid temperature is modified due to mixing of existing water with meltwater at the freezing point
+              ! more consistent with the layer tendency
+              !frazilTemperatureTendency(k,iCell) = frazilTemperatureTendency(k,iCell) + &
+              !     meltedFrazilIceThickness * oceanFreezingTemperature * config_frazil_ice_density / density(k,iCell) / dt
+              ! original formulation
               frazilTemperatureTendency(k,iCell) = frazilTemperatureTendency(k,iCell) + &
                    meltedFrazilIceThickness * oceanFreezingTemperature / dt
-
               ! keep track of new frazil ice
               sumNewThicknessWeightedSaltContent = max(0.0_RKIND,sumNewThicknessWeightedSaltContent - meltedThicknessWeightedSaltContent)
               sumNewFrazilIceThickness = max(0.0_RKIND,sumNewFrazilIceThickness - meltedFrazilIceThickness)
 
+              meltedFrazilEnergy = meltedFrazilIceThickness * config_frazil_ice_density * config_frazil_heat_of_fusion
+
+              sumNewFrazilEnergy = max(0.0_RKIND,sumNewFrazilEnergy - meltedFrazilEnergy)
             endif  ! if (sumNewFrazilIceThickness > 0.0_RKIND)
 
           endif  ! if (freezingEnergy < 0)
@@ -711,6 +723,9 @@ subroutine ocn_frazil_forcing_build_arrays(domain, meshPool, forcingPool, stateP
         accumulatedFrazilIceMassNew(iCell) = accumulatedFrazilIceMassOld(iCell) + sumNewFrazilIceThickness &
            * config_frazil_ice_density
 
+        !frazilIceEnergy > 0
+        accumulatedFrazilIceEnergyNew(iCell) = accumulatedFrazilIceEnergyOld(iCell) + sumNewFrazilEnergy
+
         accumulatedFrazilIceSalinityNew(iCell) = accumulatedFrazilIceSalinityOld(iCell) + &
            (1-usePorosity) * sumNewThicknessWeightedSaltContent
 
diff --git a/components/mpas-seaice/bld/namelist_files/namelist_definition_mpassi.xml b/components/mpas-seaice/bld/namelist_files/namelist_definition_mpassi.xml
index 0d8964e034f0..54a5153701f8 100644
--- a/components/mpas-seaice/bld/namelist_files/namelist_definition_mpassi.xml
+++ b/components/mpas-seaice/bld/namelist_files/namelist_definition_mpassi.xml
@@ -533,7 +533,7 @@ Default: Defined in namelist_defaults.xml
 	category="forcing" group="forcing">
 Type of frazil ice coupling.
 
-Valid values: 'external', 'internal', 'fresh_ice_correction'
+Valid values: 'external', 'internal', 'fresh_ice_correction', 'omega-fluxes'
 Default: Defined in namelist_defaults.xml
 </entry>
 
diff --git a/components/mpas-seaice/cime_config/buildnml b/components/mpas-seaice/cime_config/buildnml
index 4386f95711ad..d749b1424318 100755
--- a/components/mpas-seaice/cime_config/buildnml
+++ b/components/mpas-seaice/cime_config/buildnml
@@ -852,6 +852,9 @@ def buildnml(case, caseroot, compname):
             lines.append('    <var name="snowVolumeCell"/>')
             lines.append('    <var name="uVelocityGeo"/>')
             lines.append('    <var name="vVelocityGeo"/>')
+            lines.append('    <var name="frazilMassFlux"/>')
+            lines.append('    <var name="frazilSaltFlux"/>')
+            lines.append('    <var name="frazilEnthalpyFlux"/>')
             if ("WW3" in compset or "DWAV%FSD" in compset or "DWAV%ZEROS" in compset):
                 lines.append('    <var name="floeSizeDiameter"/> ')
                 lines.append('    <var name="waveSpectra"/> ')
@@ -913,6 +916,9 @@ def buildnml(case, caseroot, compname):
             lines.append('    <var name="airTemperature"/>')
             lines.append('    <var name="congelation"/>')
             lines.append('    <var name="frazilFormation"/>')
+            lines.append('    <var name="frazilMassFlux"/>')
+            lines.append('    <var name="frazilSaltFlux"/>')
+            lines.append('    <var name="frazilEnthalpyFlux"/>')
             lines.append('    <var name="snowiceFormation"/>')
             lines.append('    <var name="snowMelt"/>')
             lines.append('    <var name="surfaceIceMelt"/>')
diff --git a/components/mpas-seaice/driver/ice_comp_mct.F b/components/mpas-seaice/driver/ice_comp_mct.F
index 12e84a107bb1..2693e3e32a8c 100644
--- a/components/mpas-seaice/driver/ice_comp_mct.F
+++ b/components/mpas-seaice/driver/ice_comp_mct.F
@@ -2050,6 +2050,8 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
 !    o  swnet    -- sw: net
 !    o  q        -- ocn frazil heat flux(+) / melt potential(-)
 !    o  frazil   -- ocn frazil mass flux
+!    o  frazils  -- ocn frazil salt flux
+!    o  frazilh  -- ocn frazil enthalpy flux (+)
 !    o  bcphidry -- Black Carbon hydrophilic dry deposition flux
 !    o  bcphodry -- Black Carbon hydrophobic dry deposition flux
 !    o  bcphiwet -- Black Carbon hydrophilic wet deposition flux
@@ -2125,7 +2127,6 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
       i,n,iFreq
 
    real (kind=RKIND) :: &
-      frazilMassFlux, &
       frazilMassFluxRev
 
    type (block_type), pointer :: block_ptr
@@ -2154,7 +2155,8 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
    character(len=strKIND), pointer ::   &
       config_column_physics_type,       &
       config_thermodynamics_type,       &
-      config_ocean_surface_type
+      config_ocean_surface_type,        &
+      config_frazil_coupling_type
 
    type (field1DReal), pointer ::       &
       seaSurfaceTemperatureField,       &
@@ -2162,6 +2164,9 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
       seaFreezingTemperatureField,      &
       freezingMeltingPotentialField,    &
       frazilMassAdjustField,            &
+      frazilMassFluxField,              &
+      frazilSaltFluxField,              &
+      frazilEnthalpyFluxField,          &
       uOceanVelocityField,              &
       vOceanVelocityField,              &
       seaSurfaceTiltUField,             &
@@ -2208,6 +2213,9 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
       seaSurfaceSalinity,          &
       seaFreezingTemperature,      &
       freezingMeltingPotential,    &
+      frazilMassFlux,              &
+      frazilSaltFlux,              &
+      frazilEnthalpyFlux,          &
       frazilMassAdjust,            &
       uOceanVelocity,              &
       vOceanVelocity,              &
@@ -2280,6 +2288,7 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
       call mpas_pool_get_config(configs, "config_couple_biogeochemistry_fields", config_couple_biogeochemistry_fields)
       call mpas_pool_get_config(configs, "config_use_zaerosols", config_use_zaerosols)
       call mpas_pool_get_config(configs, "config_use_column_waves", config_use_column_waves)
+      call mpas_pool_get_config(configs, "config_frazil_coupling_type", config_frazil_coupling_type)
 
       call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool)
       call mpas_pool_get_subpool(block_ptr % structs, 'ocean_coupling', oceanCoupling)
@@ -2292,6 +2301,9 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
       call mpas_pool_get_array(oceanCoupling, 'seaSurfaceSalinity', seaSurfaceSalinity)
       call mpas_pool_get_array(oceanCoupling, 'seaFreezingTemperature', seaFreezingTemperature)
       call mpas_pool_get_array(oceanCoupling, 'freezingMeltingPotential', freezingMeltingPotential)
+      call mpas_pool_get_array(oceanCoupling, 'frazilMassFlux', frazilMassFlux)
+      call mpas_pool_get_array(oceanCoupling, 'frazilSaltFlux', frazilSaltFlux)
+      call mpas_pool_get_array(oceanCoupling, 'frazilEnthalpyFlux', frazilEnthalpyFlux)
       call mpas_pool_get_array(oceanCoupling, 'frazilMassAdjust', frazilMassAdjust)
       call mpas_pool_get_array(oceanCoupling, 'uOceanVelocity', uOceanVelocity)
       call mpas_pool_get_array(oceanCoupling, 'vOceanVelocity', vOceanVelocity)
@@ -2383,7 +2395,12 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
 
            freezingMeltingPotential(i) = x2i_i(index_x2i_Fioo_q, n)
 
-           frazilMassFlux              = x2i_i(index_x2i_Fioo_frazil, n)
+           frazilMassFlux(i)           = x2i_i(index_x2i_Fioo_frazil, n)
+
+           if (trim(config_frazil_coupling_type) == "omega-fluxes") then
+              frazilSaltFlux(i)           = x2i_i(index_x2i_Fioo_frazils, n)
+              frazilEnthalpyFlux(i)       = x2i_i(index_x2i_Fioo_frazilh, n)
+              frazilMassAdjust(i)         = 0.0_RKIND
 
            ! Now determine the sea ice mass associated with the frazil heat flux given when
            ! freezingMeltingPotential(i) is positive. This produces a revised mass flux, given
@@ -2392,14 +2409,23 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
            ! coupling step as a freshwater and salt flux. This step is required to balance mass
            ! and heat with the ocean.
 
-           call frazil_mass(freezingMeltingPotential(i), frazilMassFluxRev, seaSurfaceSalinity(i))
+           else
+              call frazil_mass(freezingMeltingPotential(i), frazilMassFluxRev, seaSurfaceSalinity(i))
 
-           frazilMassAdjust(i) = frazilMassFlux-frazilMassFluxRev
+              frazilMassAdjust(i) = frazilMassFlux(i)-frazilMassFluxRev
+              frazilSaltFlux(i) = 0.0_RKIND
+              frazilEnthalpyFlux(i) = 0.0_RKIND
+
+           endif
 
         else ! non-free surface (SOM)
 
            freezingMeltingPotential(i) = x2i_i(index_x2i_Fioo_q, n)
 
+           frazilSaltFlux(i)           = 0.0_RKIND
+           frazilEnthalpyFlux(i)       = 0.0_RKIND
+           frazilMassAdjust(i)         = 0.0_RKIND
+
         endif
 
         airLevelHeight(i)              = x2i_i(index_x2i_Sa_z, n)
@@ -2564,6 +2590,9 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
    call mpas_pool_get_field(oceanCoupling, 'seaFreezingTemperature', seaFreezingTemperatureField)
    call mpas_pool_get_field(oceanCoupling, 'freezingMeltingPotential', freezingMeltingPotentialField)
    call mpas_pool_get_field(oceanCoupling, 'frazilMassAdjust', frazilMassAdjustField)
+   call mpas_pool_get_field(oceanCoupling, 'frazilMassFlux', frazilMassFluxField)
+   call mpas_pool_get_field(oceanCoupling, 'frazilSaltFlux', frazilSaltFluxField)
+   call mpas_pool_get_field(oceanCoupling, 'frazilEnthalpyFlux', frazilEnthalpyFluxField)
    call mpas_pool_get_field(oceanCoupling, 'uOceanVelocity', uOceanVelocityField)
    call mpas_pool_get_field(oceanCoupling, 'vOceanVelocity', vOceanVelocityField)
    call mpas_pool_get_field(oceanCoupling, 'seaSurfaceTiltU', seaSurfaceTiltUField)
@@ -2627,6 +2656,9 @@ subroutine ice_import_mct(x2i_i, errorCode)!{{{
    call mpas_dmpar_exch_halo_field(seaFreezingTemperatureField)
    call mpas_dmpar_exch_halo_field(freezingMeltingPotentialField)
    call mpas_dmpar_exch_halo_field(frazilMassAdjustField)
+   call mpas_dmpar_exch_halo_field(frazilMassFluxField)
+   call mpas_dmpar_exch_halo_field(frazilSaltFluxField)
+   call mpas_dmpar_exch_halo_field(frazilEnthalpyFluxField)
    call mpas_dmpar_exch_halo_field(uOceanVelocityField)
    call mpas_dmpar_exch_halo_field(vOceanVelocityField)
    call mpas_dmpar_exch_halo_field(seaSurfaceTiltUField)
@@ -3025,12 +3057,15 @@ subroutine ice_export_mct(i2x_i, errorCode)   !{{{
             ! i/o fluxes computed by ice, as well as additional freshwater and salt calculated at the last
             ! coupling import and needed to grow sea ice from frazil passed from the ocean model in the
             ! field frazilMassAdjust.
-            i2x_i(index_i2x_Fioi_melth,n) = oceanHeatFlux(i)
-            i2x_i(index_i2x_Fioi_swpen,n) = oceanShortwaveFlux(i)
-            i2x_i(index_i2x_Fioi_meltw,n) = oceanFreshWaterFlux(i) + frazilMassAdjust(i)/ailohi
-            i2x_i(index_i2x_Fioi_salt ,n) = oceanSaltFlux(i) + seaiceReferenceSalinity*0.001_RKIND*frazilMassAdjust(i)/ailohi
-            i2x_i(index_i2x_Fioi_taux ,n) = tauxo
-            i2x_i(index_i2x_Fioi_tauy ,n) = tauyo
+            i2x_i(index_i2x_Fioi_melth  ,n) = oceanHeatFlux(i)
+            i2x_i(index_i2x_Fioi_swpen  ,n) = oceanShortwaveFlux(i)
+            i2x_i(index_i2x_Fioi_meltw  ,n) = oceanFreshWaterFlux(i) + frazilMassAdjust(i)/ailohi
+            i2x_i(index_i2x_Fioi_frazil ,n) = frazilMassAdjust(i)/ailohi
+            i2x_i(index_i2x_Fioi_salt   ,n) = oceanSaltFlux(i) + seaiceReferenceSalinity*0.001_RKIND*frazilMassAdjust(i)/ailohi
+            i2x_i(index_i2x_Fioi_frazils,n) = seaiceReferenceSalinity*0.001_RKIND*frazilMassAdjust(i)/ailohi
+            i2x_i(index_i2x_Fioi_frazilh,n) = 0.0_RKIND
+            i2x_i(index_i2x_Fioi_taux   ,n) = tauxo
+            i2x_i(index_i2x_Fioi_tauy   ,n) = tauyo
 
             ! export biogeochemistry fields, if configured
             if (config_use_column_biogeochemistry .and. config_couple_biogeochemistry_fields) then
diff --git a/components/mpas-seaice/driver/mpassi_cpl_indices.F b/components/mpas-seaice/driver/mpassi_cpl_indices.F
index 6e84fcba8533..c38387f31883 100644
--- a/components/mpas-seaice/driver/mpassi_cpl_indices.F
+++ b/components/mpas-seaice/driver/mpassi_cpl_indices.F
@@ -34,7 +34,10 @@ module mpassi_cpl_indices
   integer :: index_i2x_Fioi_swpen      ! sw: net penetrating ice
   integer :: index_i2x_Fioi_melth      ! heat  flux from melting ice (<0)
   integer :: index_i2x_Fioi_meltw      ! water flux from melting ice
-  integer :: index_i2x_Fioi_salt       ! salt  flux from meting  ice
+  integer :: index_i2x_Fioi_salt       ! salt  flux from melting  ice
+  integer :: index_i2x_Fioi_frazil     ! water flux from frazil ice
+  integer :: index_i2x_Fioi_frazils    ! salt  flux from frazil ice
+  integer :: index_i2x_Fioi_frazilh    ! heat  flux from frazil ice
   integer :: index_i2x_Fioi_taux       ! ice/ocn stress, zonal
   integer :: index_i2x_Fioi_tauy       ! ice/ocn stress, zonal
   integer :: index_i2x_Fioi_bergh      ! heat  flux from melting icebergs (<0)
@@ -87,7 +90,9 @@ module mpassi_cpl_indices
   integer :: index_x2i_Faxa_swvdf      ! sw: vis diffuse downward
   integer :: index_x2i_Faxa_swnet      ! sw: net
   integer :: index_x2i_Fioo_q          ! ocn freezing melting potential
-  integer :: index_x2i_Fioo_frazil     ! ocn frazil ice formation 
+  integer :: index_x2i_Fioo_frazil     ! ocn frazil ice mass flux
+  integer :: index_x2i_Fioo_frazils    ! ocn frazil ice salt flux
+  integer :: index_x2i_Fioo_frazilh    ! ocn frazil heat flux
   integer :: index_x2i_Faxa_bcphidry   ! flux: Black Carbon hydrophilic dry deposition
   integer :: index_x2i_Faxa_bcphodry   ! flux: Black Carbon hydrophobic dry deposition
   integer :: index_x2i_Faxa_bcphiwet   ! flux: Black Carbon hydrophilic wet deposition
@@ -176,6 +181,9 @@ subroutine mpassi_cpl_indices_set( )
     index_i2x_Fioi_melth    = mct_avect_indexra(i2x,'Fioi_melth')
     index_i2x_Fioi_meltw    = mct_avect_indexra(i2x,'Fioi_meltw')
     index_i2x_Fioi_salt     = mct_avect_indexra(i2x,'Fioi_salt')
+    index_i2x_Fioi_frazil   = mct_avect_indexra(i2x,'Fioi_frazil')
+    index_i2x_Fioi_frazils  = mct_avect_indexra(i2x,'Fioi_frazils')
+    index_i2x_Fioi_frazilh  = mct_avect_indexra(i2x,'Fioi_frazilh')
     index_i2x_Fioi_taux     = mct_avect_indexra(i2x,'Fioi_taux')
     index_i2x_Fioi_tauy     = mct_avect_indexra(i2x,'Fioi_tauy')
     index_i2x_Fioi_bergh    = mct_avect_indexra(i2x,'PFioi_bergh')
@@ -228,6 +236,8 @@ subroutine mpassi_cpl_indices_set( )
     index_x2i_Faxa_swvdf    = mct_avect_indexra(x2i,'Faxa_swvdf')
     index_x2i_Fioo_q        = mct_avect_indexra(x2i,'Fioo_q')
     index_x2i_Fioo_frazil   = mct_avect_indexra(x2i,'Fioo_frazil')
+    index_x2i_Fioo_frazils  = mct_avect_indexra(x2i,'Fioo_frazils')
+    index_x2i_Fioo_frazilh  = mct_avect_indexra(x2i,'Fioo_frazilh')
     index_x2i_Faxa_bcphidry = mct_avect_indexra(x2i,'Faxa_bcphidry')
     index_x2i_Faxa_bcphodry = mct_avect_indexra(x2i,'Faxa_bcphodry')
     index_x2i_Faxa_bcphiwet = mct_avect_indexra(x2i,'Faxa_bcphiwet')
diff --git a/components/mpas-seaice/src/Registry.xml b/components/mpas-seaice/src/Registry.xml
index de91fd55b530..46a8bd8bd8b5 100644
--- a/components/mpas-seaice/src/Registry.xml
+++ b/components/mpas-seaice/src/Registry.xml
@@ -566,7 +566,7 @@
 		/>
 		<nml_option name="config_frazil_coupling_type" type="character" default_value="external" units="unitless"
 			description="Type of frazil ice coupling."
-			possible_values="'external', 'internal', 'fresh_ice_correction'"
+			possible_values="'external', 'internal', 'fresh_ice_correction', 'omega-fluxes'"
 		/>
 		<nml_option name="config_include_pond_freshwater_feedback" type="logical" default_value="false" units="unitless"
 			description="Reduce the ocean fresh water flux by the pond fresh water flux for coupling."
@@ -2133,6 +2133,9 @@
 			<var name="oceanStressCellV"/>
 			<var name="seaSurfaceTemperature"/>
 			<var name="freezingMeltingPotential"/>
+			<var name="frazilMassFlux"/>
+			<var name="frazilSaltFlux"/>
+			<var name="frazilEnthalpyFlux"/>
 			<var name="airOceanDragCoefficientRatio"/>
 			<var name="skeletalAlgaeConc"/>
 			<var name="skeletalDOCConc"/>
@@ -2457,6 +2460,9 @@
 			<var name="uOceanVelocityVertexGeo"/>
 			<var name="vOceanVelocityVertexGeo"/>
 			<var name="freezingMeltingPotential"/>
+			<var name="frazilMassFlux"/>
+			<var name="frazilSaltFlux"/>
+			<var name="frazilEnthalpyFlux"/>
 			<var name="shortwaveScalingFactor"/>
 			<var name="airTemperature"/>
 			<var name="congelation"/>
@@ -4364,6 +4370,21 @@
 		     dimensions="nCells Time"
 		     units="kg m-2 s-1"
 		/>
+		<var name="frazilMassFlux"
+		     type="real"
+		     dimensions="nCells Time"
+		     units="kg m-2 s-1"
+		/>
+		<var name="frazilSaltFlux"
+		     type="real"
+		     dimensions="nCells Time"
+		     units="kg m-2 s-1"
+		/>
+		<var name="frazilEnthalpyFlux"
+		     type="real"
+		     dimensions="nCells Time"
+		     units="W m-2"
+		/>
 		<var name="uOceanVelocity"
 		     type="real"
 		     dimensions="nCells Time"
diff --git a/components/mpas-seaice/src/icepack b/components/mpas-seaice/src/icepack
index cf1d460fbb6a..30b06200adaa 160000
--- a/components/mpas-seaice/src/icepack
+++ b/components/mpas-seaice/src/icepack
@@ -1 +1 @@
-Subproject commit cf1d460fbb6a13d28988f798b33734f46d70a9dc
+Subproject commit 30b06200adaaa3bb5fb60605d840b23163766a51
diff --git a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F
index 0b9ed6d4cb7e..ca49db3905f4 100644
--- a/components/mpas-seaice/src/shared/mpas_seaice_icepack.F
+++ b/components/mpas-seaice/src/shared/mpas_seaice_icepack.F
@@ -2417,6 +2417,9 @@ subroutine column_itd_thermodynamics(domain, clock)
          lateralIceMeltRate, &
          lateralIceMelt, &
          freezingMeltingPotential, &
+         frazilMassFlux, &
+         frazilSaltFlux, &
+         frazilEnthalpyFlux, &
          frazilFormation, &
          rainfallRate, &
          pondFreshWaterFlux, &
@@ -2538,6 +2541,9 @@ subroutine column_itd_thermodynamics(domain, clock)
        call MPAS_pool_get_array(atmos_coupling, "rainfallRate", rainfallRate)
 
        call MPAS_pool_get_array(ocean_coupling, "freezingMeltingPotential", freezingMeltingPotential)
+       call MPAS_pool_get_array(ocean_coupling, "frazilMassFlux", frazilMassFlux)
+       call MPAS_pool_get_array(ocean_coupling, "frazilSaltFlux", frazilSaltFlux)
+       call MPAS_pool_get_array(ocean_coupling, "frazilEnthalpyFlux", frazilEnthalpyFlux)
        call MPAS_pool_get_array(ocean_coupling, "seaFreezingTemperature", seaFreezingTemperature)
        call MPAS_pool_get_array(ocean_coupling, "seaSurfaceSalinity", seaSurfaceSalinity)
 
@@ -2658,7 +2664,10 @@ subroutine column_itd_thermodynamics(domain, clock)
                d_afsd_latg=dFloeSizeLateralGrowth(:,iCell),    &
                d_afsd_newi=dFloeSizeNewIce(:,iCell),   &
                d_afsd_latm=dFloeSizeLateralMelt(:,iCell),   &
-               d_afsd_weld=dFloeSizeWeld(:,iCell))
+               d_afsd_weld=dFloeSizeWeld(:,iCell), &
+               frazilm=frazilMassFlux(iCell), &
+               frazils=frazilSaltFlux(iCell), &
+               frazilh=frazilEnthalpyFlux(iCell))
 
           do iBioTracers = 1, ciceTracerObject % nBioTracers
              oceanBioFluxes(iBioTracers,iCell) = oceanBioFluxes(iBioTracers,iCell) + oceanBioFluxesTemp(iBioTracers,iCell)
diff --git a/components/mpas-seaice/src/shared/mpas_seaice_initialize.F b/components/mpas-seaice/src/shared/mpas_seaice_initialize.F
index 2a4510bea396..db56e5f012d9 100644
--- a/components/mpas-seaice/src/shared/mpas_seaice_initialize.F
+++ b/components/mpas-seaice/src/shared/mpas_seaice_initialize.F
@@ -2780,9 +2780,10 @@ subroutine seaice_check_configs_coupled(domain)
            messageType=MPAS_LOG_CRIT)
     endif
 
-    if (trim(config_frazil_coupling_type) /= "external") then
+    if (trim(config_frazil_coupling_type) /= "external" .and. &
+      trim(config_frazil_coupling_type) /= "omega-fluxes") then
       call mpas_log_write(&
-           'seaice_check_configs_coupled: config_frazil_coupling_type must be external', &
+           'seaice_check_configs_coupled: config_frazil_coupling_type must be external or omega-fluxes', &
            messageType=MPAS_LOG_CRIT)
     endif
 
diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90
index 3c3872d48e50..9c5761ee9392 100644
--- a/driver-mct/main/seq_diag_mct.F90
+++ b/driver-mct/main/seq_diag_mct.F90
@@ -144,44 +144,50 @@ module seq_diag_mct
   integer(in),parameter :: f_hpolar    =11     ! heat : AIS imbalance
   integer(in),parameter :: f_hh2ot     =12     ! heat : water temperature
   integer(in),parameter :: f_hgsmb     =13     ! heat : Greenland ice sheet surface mass balance
-  integer(in),parameter :: f_wfrz      =14     ! water: freezing
-  integer(in),parameter :: f_wmelt     =15     ! water: melting
-  integer(in),parameter :: f_wrain     =16     ! water: precip, liquid
-  integer(in),parameter :: f_wsnow     =17     ! water: precip, frozen
-  integer(in),parameter :: f_wpolar    =18     ! water: AIS imbalance
-  integer(in),parameter :: f_wgsmb     =19     ! water: Greenland ice sheet surface mass balance
-  integer(in),parameter :: f_wevap     =20     ! water: evaporation
-  integer(in),parameter :: f_wroff     =21     ! water: runoff/flood
-  integer(in),parameter :: f_wioff     =22     ! water: frozen runoff
-  integer(in),parameter :: f_wirrig    =23     ! water: irrigation
-  integer(in),parameter :: f_wfrz_16O  =24     ! water: freezing
-  integer(in),parameter :: f_wmelt_16O =25     ! water: melting
-  integer(in),parameter :: f_wrain_16O =26     ! water: precip, liquid
-  integer(in),parameter :: f_wsnow_16O =27     ! water: precip, frozen
-  integer(in),parameter :: f_wevap_16O =28     ! water: evaporation
-  integer(in),parameter :: f_wroff_16O =29     ! water: runoff/flood
-  integer(in),parameter :: f_wioff_16O =30     ! water: frozen runoff
-  integer(in),parameter :: f_wfrz_18O  =31     ! water: freezing
-  integer(in),parameter :: f_wmelt_18O =32     ! water: melting
-  integer(in),parameter :: f_wrain_18O =33     ! water: precip, liquid
-  integer(in),parameter :: f_wsnow_18O =34     ! water: precip, frozen
-  integer(in),parameter :: f_wevap_18O =35     ! water: evaporation
-  integer(in),parameter :: f_wroff_18O =36     ! water: runoff/flood
-  integer(in),parameter :: f_wioff_18O =37     ! water: frozen runoff
-  integer(in),parameter :: f_wfrz_HDO  =38     ! water: freezing
-  integer(in),parameter :: f_wmelt_HDO =39     ! water: melting
-  integer(in),parameter :: f_wrain_HDO =40     ! water: precip, liquid
-  integer(in),parameter :: f_wsnow_HDO =41     ! water: precip, frozen
-  integer(in),parameter :: f_wevap_HDO =42     ! water: evaporation
-  integer(in),parameter :: f_wroff_HDO =43     ! water: runoff/flood
-  integer(in),parameter :: f_wioff_HDO =44     ! water: frozen runoff
-
-  integer(in),parameter :: f_size     = f_wioff_HDO   ! Total array size of all elements
+  integer(in),parameter :: f_hfrazo    =14     ! heat : frazil from ocean
+  integer(in),parameter :: f_hfrazi    =15     ! heat : frazil from ice
+  integer(in),parameter :: f_wfrz      =16     ! water: freezing (frazil ocean)
+  integer(in),parameter :: f_wfrzi     =17     ! water: frazil ice
+  integer(in),parameter :: f_wmelt     =18     ! water: melting
+  integer(in),parameter :: f_wrain     =19     ! water: precip, liquid
+  integer(in),parameter :: f_wsnow     =20     ! water: precip, frozen
+  integer(in),parameter :: f_wpolar    =21     ! water: AIS imbalance
+  integer(in),parameter :: f_wgsmb     =22     ! water: Greenland ice sheet surface mass balance
+  integer(in),parameter :: f_wevap     =23     ! water: evaporation
+  integer(in),parameter :: f_wroff     =24     ! water: runoff/flood
+  integer(in),parameter :: f_wioff     =25     ! water: frozen runoff
+  integer(in),parameter :: f_wirrig    =26     ! water: irrigation
+  integer(in),parameter :: f_wfrz_16O  =27     ! water: freezing
+  integer(in),parameter :: f_wmelt_16O =28     ! water: melting
+  integer(in),parameter :: f_wrain_16O =29     ! water: precip, liquid
+  integer(in),parameter :: f_wsnow_16O =30     ! water: precip, frozen
+  integer(in),parameter :: f_wevap_16O =31     ! water: evaporation
+  integer(in),parameter :: f_wroff_16O =32     ! water: runoff/flood
+  integer(in),parameter :: f_wioff_16O =33     ! water: frozen runoff
+  integer(in),parameter :: f_wfrz_18O  =34     ! water: freezing
+  integer(in),parameter :: f_wmelt_18O =35     ! water: melting
+  integer(in),parameter :: f_wrain_18O =36     ! water: precip, liquid
+  integer(in),parameter :: f_wsnow_18O =37     ! water: precip, frozen
+  integer(in),parameter :: f_wevap_18O =38     ! water: evaporation
+  integer(in),parameter :: f_wroff_18O =39     ! water: runoff/flood
+  integer(in),parameter :: f_wioff_18O =40     ! water: frozen runoff
+  integer(in),parameter :: f_wfrz_HDO  =41     ! water: freezing
+  integer(in),parameter :: f_wmelt_HDO =42     ! water: melting
+  integer(in),parameter :: f_wrain_HDO =43     ! water: precip, liquid
+  integer(in),parameter :: f_wsnow_HDO =44     ! water: precip, frozen
+  integer(in),parameter :: f_wevap_HDO =45     ! water: evaporation
+  integer(in),parameter :: f_wroff_HDO =46     ! water: runoff/flood
+  integer(in),parameter :: f_wioff_HDO =47     ! water: frozen runoff
+  integer(in),parameter :: f_salt      =48     ! salt: salinity flux from sea ice
+  integer(in),parameter :: f_sfrazo    =49     ! salt: salinity flux frazil from ocean
+  integer(in),parameter :: f_sfrazi    =50     ! salt: salinity flux frazil from ice
+
+  integer(in),parameter :: f_size     = f_sfrazi      ! Total array size of all elements
   integer(in),parameter :: f_a        = f_area        ! 1st index for area
   integer(in),parameter :: f_a_end    = f_area        ! last index for area
   integer(in),parameter :: f_h        = f_hfrz        ! 1st index for heat
   !integer(in),parameter :: f_h_end    = f_hh2ot       ! Last index for heat
-  integer(in),parameter :: f_h_end    = f_hgsmb       ! Last index for heat
+  integer(in),parameter :: f_h_end    = f_hfrazi      ! Last index for heat
   integer(in),parameter :: f_w        = f_wfrz        ! 1st index for water
   integer(in),parameter :: f_w_end    = f_wirrig      ! Last index for water
   integer(in),parameter :: f_16O      = f_wfrz_16O    ! 1st index for 16O water isotope
@@ -190,20 +196,24 @@ module seq_diag_mct
   integer(in),parameter :: f_16O_end  = f_wioff_16O   ! Last index for 16O water isotope
   integer(in),parameter :: f_18O_end  = f_wioff_18O   ! Last index for 18O water isotope
   integer(in),parameter :: f_HDO_end  = f_wioff_HDO   ! Last index for HDO water isotope
+  integer(in),parameter :: f_s        = f_salt        ! 1st index for salt
+  integer(in),parameter :: f_s_end    = f_sfrazi      ! last index for salt
 
   character(len=12),parameter :: fname(f_size) = &
 
        (/'        area','     hfreeze','       hmelt','      hnetsw','       hlwdn', &
        '       hlwup','     hlatvap','     hlatfus','      hiroff','        hsen', &
-       '      hpolar','    hh2otemp','       hgsmb','     wfreeze','       wmelt','       wrain', &
-       '       wsnow','      wpolar','       wgsmb','       wevap','     wrunoff','     wfrzrof', &
+       '      hpolar','    hh2otemp','       hgsmb','      hfrazo','      hfrazi', &
+       '     wfreeze','    wfreezei','       wmelt','       wrain','       wsnow', &
+       '       wpolar','      wgsmb','       wevap','     wrunoff','     wfrzrof', &
        '      wirrig',                                                             &
        ' wfreeze_16O','   wmelt_16O','   wrain_16O','   wsnow_16O',                &
        '   wevap_16O',' wrunoff_16O',' wfrzrof_16O',                               &
        ' wfreeze_18O','   wmelt_18O','   wrain_18O','   wsnow_18O',                &
        '   wevap_18O',' wrunoff_18O',' wfrzrof_18O',                               &
        ' wfreeze_HDO','   wmelt_HDO','   wrain_HDO','   wsnow_HDO',                &
-       '   wevap_HDO',' wrunoff_HDO',' wfrzrof_HDO'/)
+       '   wevap_HDO',' wrunoff_HDO',' wfrzrof_HDO',                               &
+       '        salt','      sfrazo','      sfrazi'/)
 
   !--- P for period ---
 
@@ -295,6 +305,8 @@ module seq_diag_mct
 
   integer :: index_o2x_Faoo_h2otemp
   integer :: index_o2x_Fioo_frazil
+  integer :: index_o2x_Fioo_frazils
+  integer :: index_o2x_Fioo_frazilh
   integer :: index_o2x_Foxo_frazil_li
   integer :: index_o2x_Fioo_q
   integer :: index_o2x_Foxo_q_li
@@ -325,10 +337,16 @@ module seq_diag_mct
   integer :: index_x2o_Fioi_bergh
   integer :: index_x2o_Fioi_bergw
   integer :: index_x2o_Fioi_salt
+  integer :: index_x2o_Fioi_frazil
+  integer :: index_x2o_Fioi_frazils
+  integer :: index_x2o_Fioi_frazilh
 
   integer :: index_i2x_Fioi_melth
   integer :: index_i2x_Fioi_meltw
   integer :: index_i2x_Fioi_salt
+  integer :: index_i2x_Fioi_frazil
+  integer :: index_i2x_Fioi_frazils
+  integer :: index_i2x_Fioi_frazilh
   integer :: index_i2x_Faii_swnet
   integer :: index_i2x_Fioi_swpen
   integer :: index_i2x_Faii_lwup
@@ -340,6 +358,8 @@ module seq_diag_mct
   integer :: index_x2i_Faxa_rain
   integer :: index_x2i_Faxa_snow
   integer :: index_x2i_Fioo_frazil
+  integer :: index_x2i_Fioo_frazils
+  integer :: index_x2i_Fioo_frazilh
   integer :: index_x2i_Fioo_q
   integer :: index_x2i_Fixx_rofi
 
@@ -1478,6 +1498,8 @@ subroutine seq_diag_ocn_mct( ocn, xao_o, frac_o, infodata, do_o2x, do_x2o, do_xa
     if (present(do_o2x)) then
        if (first_time) then
           index_o2x_Fioo_frazil  = mct_aVect_indexRA(o2x_o,'Fioo_frazil')
+          index_o2x_Fioo_frazils = mct_aVect_indexRA(o2x_o,'Fioo_frazils')
+          index_o2x_Fioo_frazilh = mct_aVect_indexRA(o2x_o,'Fioo_frazilh')
           index_o2x_Foxo_frazil_li= mct_aVect_indexRA(o2x_o,'Foxo_frazil_li',perrWith='quiet')
           index_o2x_Fioo_q       = mct_aVect_indexRA(o2x_o,'Fioo_q')
           index_o2x_Foxo_q_li    = mct_aVect_indexRA(o2x_o,'Foxo_q_li',perrWith='quiet')
@@ -1501,6 +1523,8 @@ subroutine seq_diag_ocn_mct( ocn, xao_o, frac_o, infodata, do_o2x, do_x2o, do_xa
           nf = f_area; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_o
           nf = f_wfrz;  budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - (ca_o+ca_i)*max(0.0_r8,o2x_o%rAttr(index_o2x_Fioo_frazil,n))
           nf = f_hfrz;  budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*max(0.0_r8,o2x_o%rAttr(index_o2x_Fioo_q,n))
+          nf = f_hfrazo;  budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*max(0.0_r8,o2x_o%rAttr(index_o2x_Fioo_frazilh,n))
+          nf = f_sfrazo;  budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - (ca_o+ca_i)*o2x_o%rAttr(index_o2x_Fioo_frazils,n)
           nf = f_hh2ot; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*o2x_o%rAttr(index_o2x_Faoo_h2otemp,n)
           if (flds_polar) then
              nf = f_wpolar;budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_c*o2x_o%rAttr(index_o2x_Foxo_frazil_li,n)
@@ -1560,8 +1584,11 @@ subroutine seq_diag_ocn_mct( ocn, xao_o, frac_o, infodata, do_o2x, do_x2o, do_xa
           index_x2o_Fioi_meltw  = mct_aVect_indexRA(x2o_o,'Fioi_meltw')
           index_x2o_Fioi_bergh  = mct_aVect_indexRA(x2o_o,'PFioi_bergh')
           index_x2o_Fioi_bergw  = mct_aVect_indexRA(x2o_o,'PFioi_bergw')
-          index_x2o_Fioi_salt   = mct_aVect_indexRA(x2o_o,'Fioi_salt')
           index_x2o_Foxx_swnet  = mct_aVect_indexRA(x2o_o,'Foxx_swnet')
+          index_x2o_Fioi_salt   = mct_aVect_indexRA(x2o_o,'Fioi_salt')
+          index_x2o_Fioi_frazils = mct_aVect_indexRA(x2o_o,'Fioi_frazils')
+          index_x2o_Fioi_frazilh = mct_aVect_indexRA(x2o_o,'Fioi_frazilh')
+          index_x2o_Fioi_frazil = mct_aVect_indexRA(x2o_o,'Fioi_frazil')
           index_x2o_Faxa_lwdn   = mct_aVect_indexRA(x2o_o,'Faxa_lwdn')
           index_x2o_Faxa_rain   = mct_aVect_indexRA(x2o_o,'Faxa_rain')
           index_x2o_Faxa_snow   = mct_aVect_indexRA(x2o_o,'Faxa_snow')
@@ -1624,7 +1651,10 @@ subroutine seq_diag_ocn_mct( ocn, xao_o, frac_o, infodata, do_o2x, do_x2o, do_xa
           nf = f_wpolar; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Fioi_bergw,n)
           nf = f_wroff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Foxx_rofl,n)
           nf = f_wioff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Foxx_rofi,n)
-
+          nf = f_salt  ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Fioi_salt,n)
+          nf = f_wfrzi ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Fioi_frazil,n)
+          nf = f_hfrazi; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Fioi_frazilh,n)
+          nf = f_sfrazi; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*x2o_o%rAttr(index_x2o_Fioi_frazils,n)
 
           if ( flds_wiso_ocn )then
              nf = f_wmelt_16O;
@@ -1747,6 +1777,9 @@ subroutine seq_diag_ice_mct( ice, frac_i, infodata, do_i2x, do_x2i)
        index_i2x_Faii_sen    = mct_aVect_indexRA(i2x_i,'Faii_sen')
        index_i2x_Faii_evap   = mct_aVect_indexRA(i2x_i,'Faii_evap')
        index_i2x_Fioi_salt   = mct_aVect_indexRA(i2x_i,'Fioi_salt')
+       index_i2x_Fioi_frazil = mct_aVect_indexRA(i2x_i,'Fioi_frazil')
+       index_i2x_Fioi_frazils= mct_aVect_indexRA(i2x_i,'Fioi_frazils')
+       index_i2x_Fioi_frazilh= mct_aVect_indexRA(i2x_i,'Fioi_frazilh')
 
        index_i2x_Fioi_meltw_16O   = mct_aVect_indexRA(i2x_i,'Fioi_meltw_16O',perrWith='quiet')
        if ( index_i2x_Fioi_meltw_16O /= 0 ) flds_wiso_ice = .true.
@@ -1777,6 +1810,10 @@ subroutine seq_diag_ice_mct( ice, frac_i, infodata, do_i2x, do_x2i)
           nf = f_hsen  ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_i*i2x_i%rAttr(index_i2x_Faii_sen,n)
           nf = f_wmelt ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_i*i2x_i%rAttr(index_i2x_Fioi_meltw,n)
           nf = f_wevap ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_i*i2x_i%rAttr(index_i2x_Faii_evap,n)
+          nf = f_salt  ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_i*i2x_i%rAttr(index_i2x_Fioi_salt,n)
+          nf = f_wfrzi ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_i*i2x_i%rAttr(index_i2x_Fioi_frazil,n)
+          nf = f_sfrazi; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_i*i2x_i%rAttr(index_i2x_Fioi_frazils,n)
+          nf = f_hfrazi; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_i*i2x_i%rAttr(index_i2x_Fioi_frazilh,n)
 
           if ( flds_wiso_ice )then
              nf = f_wmelt_16O;
@@ -1808,6 +1845,8 @@ subroutine seq_diag_ice_mct( ice, frac_i, infodata, do_i2x, do_x2i)
           index_x2i_Faxa_rain   = mct_aVect_indexRA(x2i_i,'Faxa_rain')
           index_x2i_Faxa_snow   = mct_aVect_indexRA(x2i_i,'Faxa_snow')
           index_x2i_Fioo_frazil = mct_aVect_indexRA(x2i_i,'Fioo_frazil')
+          index_x2i_Fioo_frazils= mct_aVect_indexRA(x2i_i,'Fioo_frazils')
+          index_x2i_Fioo_frazilh= mct_aVect_indexRA(x2i_i,'Fioo_frazilh')
           index_x2i_Fioo_q      = mct_aVect_indexRA(x2i_i,'Fioo_q')
           index_x2i_Fixx_rofi   = mct_aVect_indexRA(x2i_i,'Fixx_rofi')
 
@@ -1840,6 +1879,10 @@ subroutine seq_diag_ice_mct( ice, frac_i, infodata, do_i2x, do_x2i)
 
           nf = f_wfrz ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + &
                (ca_o+ca_i)*max(0.0_r8,x2i_i%rAttr(index_x2i_Fioo_frazil,n))
+          nf = f_sfrazo ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + &
+               (ca_o+ca_i)*max(0.0_r8,x2i_i%rAttr(index_x2i_Fioo_frazils,n))
+          nf = f_hfrazo ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - &
+               (ca_o+ca_i)*max(0.0_r8,x2i_i%rAttr(index_x2i_Fioo_frazilh,n))
           nf = f_hfrz ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - &
                (ca_o+ca_i)*max(0.0_r8,x2i_i%rAttr(index_x2i_Fioo_q,n))
 
@@ -2081,6 +2124,21 @@ SUBROUTINE seq_diag_print_mct(EClock, stop_alarm, do_bgc_budg, &
                            sum(dataGpr(iso0(is):isof(is),ico,ip))
                    end do
                 end if
+             
+                write(logunit,*) ' '
+                write(logunit,FAH) subname,trim(str)//' SALT BUDGET ((kg/s)/m^2*1e8): period = ',trim(pname(ip)),': date = ',cdate,sec
+                write(logunit,FA0) cname(ica),cname(icl),cname(icn),cname(ics),cname(ico),' *SUM*  '
+                do nf = f_s, f_s_end
+                   write(logunit,FA1)    fname(nf),dataGpr(nf,ica,ip),dataGpr(nf,icl,ip), &
+                        dataGpr(nf,icn,ip),dataGpr(nf,ics,ip),dataGpr(nf,ico,ip), &
+                        dataGpr(nf,ica,ip)+dataGpr(nf,icl,ip)+ &
+                        dataGpr(nf,icn,ip)+dataGpr(nf,ics,ip)+dataGpr(nf,ico,ip)
+                enddo
+                write(logunit,FA1)    '   *SUM*'   ,sum(dataGpr(f_s:f_s_end,ica,ip)),sum(dataGpr(f_s:f_s_end,icl,ip)), &
+                     sum(dataGpr(f_s:f_s_end,icn,ip)),sum(dataGpr(f_s:f_s_end,ics,ip)),sum(dataGpr(f_s:f_s_end,ico,ip)), &
+                     sum(dataGpr(f_s:f_s_end,ica,ip))+sum(dataGpr(f_s:f_s_end,icl,ip))+ &
+                     sum(dataGpr(f_s:f_s_end,icn,ip))+sum(dataGpr(f_s:f_s_end,ics,ip))+sum(dataGpr(f_s:f_s_end,ico,ip))
+                write(logunit,*) ' '
 
              enddo
           endif   ! plev
@@ -2187,6 +2245,22 @@ SUBROUTINE seq_diag_print_mct(EClock, stop_alarm, do_bgc_budg, &
                            -sum(dataGpr(iso0(is):isof(is),icas,ip))
                    end do
                 end if
+
+                write(logunit,*) ' '
+                write(logunit,FAH) subname,trim(str)//' SALT BUDGET (kg/m2s*1e8): period = ',trim(pname(ip)),': date = ',cdate,sec
+                write(logunit,FA0) cname(icar),cname(icxs),cname(icxr),cname(icas),' *SUM*  '
+                do nf = f_s, f_s_end
+                   write(logunit,FA1)    fname(nf),-dataGpr(nf,icar,ip),dataGpr(nf,icxs,ip), &
+                        dataGpr(nf,icxr,ip),-dataGpr(nf,icas,ip), &
+                        -dataGpr(nf,icar,ip)+dataGpr(nf,icxs,ip)+ &
+                        dataGpr(nf,icxr,ip)-dataGpr(nf,icas,ip)
+                enddo
+                write(logunit,FA1)    '   *SUM*'   ,sum(dataGpr(f_s:f_s_end,ica,ip)),sum(dataGpr(f_s:f_s_end,icl,ip)), &
+                     sum(dataGpr(f_s:f_s_end,icn,ip)),sum(dataGpr(f_s:f_s_end,ics,ip)),sum(dataGpr(f_s:f_s_end,ico,ip)), &
+                     sum(dataGpr(f_s:f_s_end,ica,ip))+sum(dataGpr(f_s:f_s_end,icl,ip))+ &
+                     sum(dataGpr(f_s:f_s_end,icn,ip))+sum(dataGpr(f_s:f_s_end,ics,ip))+sum(dataGpr(f_s:f_s_end,ico,ip))
+                write(logunit,*) ' '
+
              enddo
           endif   ! plev
 
@@ -2335,6 +2409,28 @@ SUBROUTINE seq_diag_print_mct(EClock, stop_alarm, do_bgc_budg, &
                 end do
              end if
 
+             write(logunit,*) ' '
+             write(logunit,FAH) subname,'NET SALT BUDGET (kg/m2s*1e8): period = ',trim(pname(ip)),': date = ',cdate,sec
+             write(logunit,FA0r) '     atm','     lnd','     rof','     ocn','  ice nh','  ice sh','     glc',' *SUM*  '
+             do nf = f_s, f_s_end
+                write(logunit,FA1r)   fname(nf),dataGpr(nf,c_atm_ar,ip)+dataGpr(nf,c_atm_as,ip), &
+                     dataGpr(nf,c_lnd_lr,ip)+dataGpr(nf,c_lnd_ls,ip), &
+                     dataGpr(nf,c_rof_rr,ip)+dataGpr(nf,c_rof_rs,ip), &
+                     dataGpr(nf,c_ocn_or,ip)+dataGpr(nf,c_ocn_os,ip), &
+                     dataGpr(nf,c_inh_ir,ip)+dataGpr(nf,c_inh_is,ip), &
+                     dataGpr(nf,c_ish_ir,ip)+dataGpr(nf,c_ish_is,ip), &
+                     dataGpr(nf,c_glc_gr,ip)+dataGpr(nf,c_glc_gs,ip), &
+                     dataGpr(nf,c_atm_ar,ip)+dataGpr(nf,c_atm_as,ip)+ &
+                     dataGpr(nf,c_lnd_lr,ip)+dataGpr(nf,c_lnd_ls,ip)+ &
+                     dataGpr(nf,c_rof_rr,ip)+dataGpr(nf,c_rof_rs,ip)+ &
+                     dataGpr(nf,c_ocn_or,ip)+dataGpr(nf,c_ocn_os,ip)+ &
+                     dataGpr(nf,c_inh_ir,ip)+dataGpr(nf,c_inh_is,ip)+ &
+                     dataGpr(nf,c_ish_ir,ip)+dataGpr(nf,c_ish_is,ip)+ &
+                     dataGpr(nf,c_glc_gr,ip)+dataGpr(nf,c_glc_gs,ip)
+             enddo
+             write(logunit,*) ' '
+
+
           endif
 
           ! ---- doprint ---- doprint ---- doprint ----
diff --git a/driver-mct/shr/seq_flds_mod.F90 b/driver-mct/shr/seq_flds_mod.F90
index f9960c1761ec..7eca8224f923 100644
--- a/driver-mct/shr/seq_flds_mod.F90
+++ b/driver-mct/shr/seq_flds_mod.F90
@@ -1602,6 +1602,24 @@ subroutine seq_flds_set(nmlfile, ID, infodata)
     attname  = 'Fioo_frazil'
     call metadata_set(attname, longname, stdname, units)
 
+    ! Ocean frazil salt flux
+    call seq_flds_add(o2x_fluxes,"Fioo_frazils")
+    call seq_flds_add(x2i_fluxes,"Fioo_frazils")
+    longname = 'Ocean frazil salt flux'
+    stdname  = 'ocean_frazil_salt_flux'
+    units    = 'kg m-2 s-1'
+    attname  = 'Fioo_frazils'
+    call metadata_set(attname, longname, stdname, units)
+
+    ! Ocean frazil head flux
+    call seq_flds_add(o2x_fluxes,"Fioo_frazilh")
+    call seq_flds_add(x2i_fluxes,"Fioo_frazilh")
+    longname = 'Ocean frazil heat flux'
+    stdname  = 'ocean_frazil_heat_flux'
+    units    = 'W m-2'
+    attname  = 'Fioo_frazilh'
+    call metadata_set(attname, longname, stdname, units)
+
     ! Heat flux from melting
     call seq_flds_add(i2x_fluxes,"Fioi_melth")
     call seq_flds_add(x2o_fluxes,"Fioi_melth")
@@ -1711,6 +1729,33 @@ subroutine seq_flds_set(nmlfile, ID, infodata)
     attname  = 'Fioi_salt'
     call metadata_set(attname, longname, stdname, units)
 
+    ! Frazil Mass flux
+    call seq_flds_add(i2x_fluxes,"Fioi_frazil")
+    call seq_flds_add(x2o_fluxes,"Fioi_frazil")
+    longname = 'Frazil mass flux'
+    stdname  = 'frazil_mass_flux_into_sea_water'
+    units    = 'kg m-2 s-1'
+    attname  = 'Fioi_frazil'
+    call metadata_set(attname, longname, stdname, units)
+
+    ! Frazil Salt flux
+    call seq_flds_add(i2x_fluxes,"Fioi_frazils")
+    call seq_flds_add(x2o_fluxes,"Fioi_frazils")
+    longname = 'Frazil salt flux'
+    stdname  = 'frazil_salt_flux_into_sea_water'
+    units    = 'kg m-2 s-1'
+    attname  = 'Fioi_frazils'
+    call metadata_set(attname, longname, stdname, units)
+
+    ! Frazil Heat flux
+    call seq_flds_add(i2x_fluxes,"Fioi_frazilh")
+    call seq_flds_add(x2o_fluxes,"Fioi_frazilh")
+    longname = 'Frazil heat flux'
+    stdname  = 'frazil_heat_flux_into_sea_water'
+    units    = 'W m-2'
+    attname  = 'Fioi_frazilh'
+    call metadata_set(attname, longname, stdname, units)
+
     ! Black Carbon hydrophilic deposition
     call seq_flds_add(i2x_fluxes,"Fioi_bcphi" )
     call seq_flds_add(x2o_fluxes,"Fioi_bcphi"   )
diff --git a/driver-moab/shr/seq_flds_mod.F90 b/driver-moab/shr/seq_flds_mod.F90
index 99d6580183e9..70bf634d672a 100644
--- a/driver-moab/shr/seq_flds_mod.F90
+++ b/driver-moab/shr/seq_flds_mod.F90
@@ -1617,6 +1617,24 @@ subroutine seq_flds_set(nmlfile, ID, infodata)
     attname  = 'Fioo_frazil'
     call metadata_set(attname, longname, stdname, units)
 
+    ! Ocean frazil salt production
+    call seq_flds_add(o2x_fluxes,"Fioo_frazils")
+    call seq_flds_add(x2i_fluxes,"Fioo_frazils")
+    longname = 'Ocean frazil salt production'
+    stdname  = 'ocean_frazil_salt_production'
+    units    = 'kg m-2 s-1'
+    attname  = 'Fioo_frazils'
+    call metadata_set(attname, longname, stdname, units)
+ 
+    ! Ocean frazil heat flux 
+    call seq_flds_add(o2x_fluxes,"Fioo_frazilh")
+    call seq_flds_add(x2i_fluxes,"Fioo_frazilh")
+    longname = 'Ocean frazil heat flux'
+    stdname  = 'ocean_frazil_heat_flux'
+    units    = 'W m-2'
+    attname  = 'Fioo_frazilh'
+    call metadata_set(attname, longname, stdname, units)
+
     ! Heat flux from melting
     call seq_flds_add(i2x_fluxes,"Fioi_melth")
     call seq_flds_add(x2o_fluxes,"Fioi_melth")
@@ -1725,6 +1743,33 @@ subroutine seq_flds_set(nmlfile, ID, infodata)
     attname  = 'Fioi_salt'
     call metadata_set(attname, longname, stdname, units)
 
+    ! Frazil mass flux
+    call seq_flds_add(i2x_fluxes,"Fioi_frazil")
+    call seq_flds_add(x2o_fluxes,"Fioi_frazil")
+    longname = 'Frazil mass flux'
+    stdname  = 'frazil_mass_flux_into_sea_water'
+    units    = 'kg m-2 s-1'
+    attname  = 'Fioi_frazil'
+    call metadata_set(attname, longname, stdname, units)
+
+    ! Frazil salt flux
+    call seq_flds_add(i2x_fluxes,"Fioi_frazils")
+    call seq_flds_add(x2o_fluxes,"Fioi_frazils")
+    longname = 'Frazil salt flux'
+    stdname  = 'frazil_salt_flux_into_sea_water'
+    units    = 'kg m-2 s-1'
+    attname  = 'Fioi_frazils'
+    call metadata_set(attname, longname, stdname, units)
+
+    ! Frazil heat flux
+    call seq_flds_add(i2x_fluxes,"Fioi_frazilh")
+    call seq_flds_add(x2o_fluxes,"Fioi_frazilh")
+    longname = 'Frazil heat flux'
+    stdname  = 'frazil_heat_flux_into_sea_water'
+    units    = 'W m-2'
+    attname  = 'Fioi_frazilh'
+    call metadata_set(attname, longname, stdname, units)
+
     ! Black Carbon hydrophilic deposition
     call seq_flds_add(i2x_fluxes,"Fioi_bcphi" )
     call seq_flds_add(x2o_fluxes,"Fioi_bcphi"   )