diff --git a/config_src/mct_driver/MOM_surface_forcing.F90 b/config_src/mct_driver/MOM_surface_forcing.F90
index 173ffd1e3d..fc9e7b7eeb 100644
--- a/config_src/mct_driver/MOM_surface_forcing.F90
+++ b/config_src/mct_driver/MOM_surface_forcing.F90
@@ -152,32 +152,34 @@ module MOM_surface_forcing
 ! the elements, units, and conventions that exactly conform to the use for
 ! MOM-based coupled models.
 type, public :: ice_ocean_boundary_type
-  real, pointer, dimension(:,:) :: latent_flux     =>NULL() !< latent flux (W/m2)
-  real, pointer, dimension(:,:) :: rofl_flux       =>NULL() !< liquid runoff (W/m2)
-  real, pointer, dimension(:,:) :: rofi_flux       =>NULL() !< ice runoff (W/m2)
-  real, pointer, dimension(:,:) :: u_flux          =>NULL() !< i-direction wind stress (Pa)
-  real, pointer, dimension(:,:) :: v_flux          =>NULL() !< j-direction wind stress (Pa)
-  real, pointer, dimension(:,:) :: t_flux          =>NULL() !< sensible heat flux (W/m2)
-  real, pointer, dimension(:,:) :: q_flux          =>NULL() !< specific humidity flux (kg/m2/s)
-  real, pointer, dimension(:,:) :: salt_flux       =>NULL() !< salt flux (kg/m2/s)
-  real, pointer, dimension(:,:) :: lw_flux         =>NULL() !< long wave radiation (W/m2)
-  real, pointer, dimension(:,:) :: sw_flux_vis_dir =>NULL() !< direct visible sw radiation (W/m2)
-  real, pointer, dimension(:,:) :: sw_flux_vis_dif =>NULL() !< diffuse visible sw radiation (W/m2)
-  real, pointer, dimension(:,:) :: sw_flux_nir_dir =>NULL() !< direct Near InfraRed sw radiation (W/m2)
-  real, pointer, dimension(:,:) :: sw_flux_nir_dif =>NULL() !< diffuse Near InfraRed sw radiation (W/m2)
-  real, pointer, dimension(:,:) :: lprec           =>NULL() !< mass flux of liquid precip (kg/m2/s)
-  real, pointer, dimension(:,:) :: fprec           =>NULL() !< mass flux of frozen precip (kg/m2/s)
-  real, pointer, dimension(:,:) :: runoff          =>NULL() !< mass flux of liquid runoff (kg/m2/s)
-  real, pointer, dimension(:,:) :: calving         =>NULL() !< mass flux of frozen runoff (kg/m2/s)
-  real, pointer, dimension(:,:) :: ustar_berg      =>NULL() !< frictional velocity beneath icebergs [m s-1]
-  real, pointer, dimension(:,:) :: area_berg       =>NULL() !< area covered by icebergs(m2/m2)
-  real, pointer, dimension(:,:) :: mass_berg       =>NULL() !< mass of icebergs(kg/m2)
-  real, pointer, dimension(:,:) :: runoff_hflx     =>NULL() !< heat content of liquid runoff (W/m2)
-  real, pointer, dimension(:,:) :: calving_hflx    =>NULL() !< heat content of frozen runoff (W/m2)
-  real, pointer, dimension(:,:) :: p               =>NULL() !< pressure of overlying ice and atmosphere
-                                                            !< on ocean surface (Pa)
-  real, pointer, dimension(:,:) :: mi              =>NULL() !< mass of ice (kg/m2)
-  real, pointer, dimension(:,:) :: ice_rigidity    =>NULL() !< rigidity of the sea ice, sea-ice and
+  real, pointer, dimension(:,:) :: latent_flux      =>NULL() !< latent flux (W/m2)
+  real, pointer, dimension(:,:) :: rofl_flux        =>NULL() !< liquid runoff (W/m2)
+  real, pointer, dimension(:,:) :: rofi_flux        =>NULL() !< ice runoff (W/m2)
+  real, pointer, dimension(:,:) :: u_flux           =>NULL() !< i-direction wind stress (Pa)
+  real, pointer, dimension(:,:) :: v_flux           =>NULL() !< j-direction wind stress (Pa)
+  real, pointer, dimension(:,:) :: t_flux           =>NULL() !< sensible heat flux (W/m2)
+  real, pointer, dimension(:,:) :: seaice_melt_heat =>NULL() !< sea ice and snow melt heat flux (W/m2)
+  real, pointer, dimension(:,:) :: seaice_melt      =>NULL() !< water flux due to sea ice and snow melting (kg/m2/s)
+  real, pointer, dimension(:,:) :: q_flux           =>NULL() !< specific humidity flux (kg/m2/s)
+  real, pointer, dimension(:,:) :: salt_flux        =>NULL() !< salt flux (kg/m2/s)
+  real, pointer, dimension(:,:) :: lw_flux          =>NULL() !< long wave radiation (W/m2)
+  real, pointer, dimension(:,:) :: sw_flux_vis_dir  =>NULL() !< direct visible sw radiation (W/m2)
+  real, pointer, dimension(:,:) :: sw_flux_vis_dif  =>NULL() !< diffuse visible sw radiation (W/m2)
+  real, pointer, dimension(:,:) :: sw_flux_nir_dir  =>NULL() !< direct Near InfraRed sw radiation (W/m2)
+  real, pointer, dimension(:,:) :: sw_flux_nir_dif  =>NULL() !< diffuse Near InfraRed sw radiation (W/m2)
+  real, pointer, dimension(:,:) :: lprec            =>NULL() !< mass flux of liquid precip (kg/m2/s)
+  real, pointer, dimension(:,:) :: fprec            =>NULL() !< mass flux of frozen precip (kg/m2/s)
+  real, pointer, dimension(:,:) :: runoff           =>NULL() !< mass flux of liquid runoff (kg/m2/s)
+  real, pointer, dimension(:,:) :: calving          =>NULL() !< mass flux of frozen runoff (kg/m2/s)
+  real, pointer, dimension(:,:) :: ustar_berg       =>NULL() !< frictional velocity beneath icebergs [m s-1]
+  real, pointer, dimension(:,:) :: area_berg        =>NULL() !< area covered by icebergs(m2/m2)
+  real, pointer, dimension(:,:) :: mass_berg        =>NULL() !< mass of icebergs(kg/m2)
+  real, pointer, dimension(:,:) :: runoff_hflx      =>NULL() !< heat content of liquid runoff (W/m2)
+  real, pointer, dimension(:,:) :: calving_hflx     =>NULL() !< heat content of frozen runoff (W/m2)
+  real, pointer, dimension(:,:) :: p                =>NULL() !< pressure of overlying ice and atmosphere
+                                                             !< on ocean surface (Pa)
+  real, pointer, dimension(:,:) :: mi               =>NULL() !< mass of ice (kg/m2)
+  real, pointer, dimension(:,:) :: ice_rigidity     =>NULL() !< rigidity of the sea ice, sea-ice and
                                                             !! ice-shelves, expressed as a coefficient
                                                             !! for divergence damping, as determined
                                                             !! outside of the ocean model in (m3/s)
@@ -252,6 +254,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
   real :: delta_sst           ! temporary storage for sst diff from restoring value
 
   real :: C_p                 ! heat capacity of seawater ( J/(K kg) )
+  real :: sign_for_net_FW_bug ! Should be +1. but an old bug can be recovered by using -1.
 
   call cpu_clock_begin(id_clock_forcing)
 
@@ -444,6 +447,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
     if (associated(fluxes%sens)) &
       fluxes%sens(i,j) = G%mask2dT(i,j) * IOB%t_flux(i-i0,j-j0)
 
+    ! sea ice and snow melt heat flux (W/m2)
+    if (associated(fluxes%seaice_melt_heat)) &
+      fluxes%seaice_melt_heat(i,j) = G%mask2dT(i,j) * IOB%seaice_melt_heat(i-i0,j-j0)
+
+    ! water flux due to sea ice and snow melt (kg/m2/s)
+    if (associated(fluxes%seaice_melt)) &
+      fluxes%seaice_melt(i,j) = G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0)
+
     ! latent heat flux (W/m^2)
     if (associated(fluxes%latent)) &
       fluxes%latent(i,j) = G%mask2dT(i,j) * IOB%latent_flux(i-i0,j-j0)
@@ -463,17 +474,19 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
     ! salt flux
     ! more salt restoring logic
     if (associated(fluxes%salt_flux)) &
-      fluxes%salt_flux(i,j) = G%mask2dT(i,j)*(IOB%salt_flux(i-i0,j-j0) + fluxes%salt_flux(i,j))
+      fluxes%salt_flux(i,j) = G%mask2dT(i,j)*(fluxes%salt_flux(i,j) - IOB%salt_flux(i-i0,j-j0))
 
     if (associated(fluxes%salt_flux_in)) &
-      fluxes%salt_flux_in(i,j) = G%mask2dT(i,j)*IOB%salt_flux(i-i0,j-j0)
+      fluxes%salt_flux_in(i,j) = G%mask2dT(i,j)*(-IOB%salt_flux(i-i0,j-j0))
 
   enddo; enddo
 
   ! adjust the NET fresh-water flux to zero, if flagged
   if (CS%adjust_net_fresh_water_to_zero) then
+    sign_for_net_FW_bug = 1.
+    if (CS%use_net_FW_adjustment_sign_bug) sign_for_net_FW_bug = -1.
     do j=js,je ; do i=is,ie
-      net_FW(i,j) = (((fluxes%lprec(i,j)   + fluxes%fprec(i,j)) + &
+      net_FW(i,j) = (((fluxes%lprec(i,j)   + fluxes%fprec(i,j) + fluxes%seaice_melt(i,j)) + &
           (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + &
           (fluxes%evap(i,j)    + fluxes%vprec(i,j)) ) * G%areaT(i,j)
       !   The following contribution appears to be calculating the volume flux of sea-ice
@@ -482,9 +495,12 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
       !   Bob thinks this is trying ensure the net fresh-water of the ocean + sea-ice system
       ! is constant.
       !   To do this correctly we will need a sea-ice melt field added to IOB. -AJA
+      ! GMM: as stated above, the following is wrong. CIME deals with volume/mass and
+      ! heat from sea ice/snow via seaice_melt and seaice_melt_heat, respectively.
       if (associated(fluxes%salt_flux) .and. (CS%ice_salt_concentration>0.0)) &
           net_FW(i,j) = net_FW(i,j) + G%areaT(i,j) * &
           (fluxes%salt_flux(i,j) / CS%ice_salt_concentration)
+
       net_FW2(i,j) = net_FW(i,j)/G%areaT(i,j)
     enddo; enddo
 
@@ -599,13 +615,14 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
         forces%p_surf(i,j) = forces%p_surf_full(i,j)
       endif
 
-      if (CS%use_limited_P_SSH) then
-        forces%p_surf_SSH => forces%p_surf
-      else
-        forces%p_surf_SSH => forces%p_surf_full
-      endif
-    end if
-  end do; end do
+    endif
+  enddo; enddo
+
+  if (CS%use_limited_P_SSH) then
+    forces%p_surf_SSH => forces%p_surf
+  else
+    forces%p_surf_SSH => forces%p_surf_full
+  endif
 
   ! GMM, CIME uses AGRID. All the BGRID_NE code can be cleaned later
   wind_stagger = AGRID
@@ -774,6 +791,8 @@ subroutine IOB_allocate(IOB, isc, iec, jsc, jec)
              IOB% u_flux (isc:iec,jsc:jec),          &
              IOB% v_flux (isc:iec,jsc:jec),          &
              IOB% t_flux (isc:iec,jsc:jec),          &
+             IOB% seaice_melt_heat (isc:iec,jsc:jec),&
+             IOB% seaice_melt (isc:iec,jsc:jec),           &
              IOB% q_flux (isc:iec,jsc:jec),          &
              IOB% salt_flux (isc:iec,jsc:jec),       &
              IOB% lw_flux (isc:iec,jsc:jec),         &
@@ -783,7 +802,6 @@ subroutine IOB_allocate(IOB, isc, iec, jsc, jec)
              IOB% sw_flux_nir_dif (isc:iec,jsc:jec), &
              IOB% lprec (isc:iec,jsc:jec),           &
              IOB% fprec (isc:iec,jsc:jec),           &
-             IOB% runoff (isc:iec,jsc:jec),          &
              IOB% ustar_berg (isc:iec,jsc:jec),      &
              IOB% area_berg (isc:iec,jsc:jec),       &
              IOB% mass_berg (isc:iec,jsc:jec),       &
@@ -793,30 +811,31 @@ subroutine IOB_allocate(IOB, isc, iec, jsc, jec)
              IOB% mi (isc:iec,jsc:jec),              &
              IOB% p (isc:iec,jsc:jec))
 
-  IOB%latent_flux     = 0.0
-  IOB%rofl_flux       = 0.0
-  IOB%rofi_flux       = 0.0
-  IOB%u_flux          = 0.0
-  IOB%v_flux          = 0.0
-  IOB%t_flux          = 0.0
-  IOB%q_flux          = 0.0
-  IOB%salt_flux       = 0.0
-  IOB%lw_flux         = 0.0
-  IOB%sw_flux_vis_dir = 0.0
-  IOB%sw_flux_vis_dif = 0.0
-  IOB%sw_flux_nir_dir = 0.0
-  IOB%sw_flux_nir_dif = 0.0
-  IOB%lprec           = 0.0
-  IOB%fprec           = 0.0
-  IOB%runoff          = 0.0
-  IOB%ustar_berg      = 0.0
-  IOB%area_berg       = 0.0
-  IOB%mass_berg       = 0.0
-  IOB%calving         = 0.0
-  IOB%runoff_hflx     = 0.0
-  IOB%calving_hflx    = 0.0
-  IOB%mi              = 0.0
-  IOB%p               = 0.0
+  IOB%latent_flux      = 0.0
+  IOB%rofl_flux        = 0.0
+  IOB%rofi_flux        = 0.0
+  IOB%u_flux           = 0.0
+  IOB%v_flux           = 0.0
+  IOB%t_flux           = 0.0
+  IOB%seaice_melt_heat = 0.0
+  IOB%seaice_melt      = 0.0
+  IOB%q_flux           = 0.0
+  IOB%salt_flux        = 0.0
+  IOB%lw_flux          = 0.0
+  IOB%sw_flux_vis_dir  = 0.0
+  IOB%sw_flux_vis_dif  = 0.0
+  IOB%sw_flux_nir_dir  = 0.0
+  IOB%sw_flux_nir_dif  = 0.0
+  IOB%lprec            = 0.0
+  IOB%fprec            = 0.0
+  IOB%ustar_berg       = 0.0
+  IOB%area_berg        = 0.0
+  IOB%mass_berg        = 0.0
+  IOB%calving          = 0.0
+  IOB%runoff_hflx      = 0.0
+  IOB%calving_hflx     = 0.0
+  IOB%mi               = 0.0
+  IOB%p                = 0.0
 
 end subroutine IOB_allocate
 
@@ -1048,6 +1067,11 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
                  CS%adjust_net_fresh_water_to_zero, &
                  "If true, adjusts the net fresh-water forcing seen \n"//&
                  "by the ocean (including restoring) to zero.", default=.false.)
+  if (CS%adjust_net_fresh_water_to_zero) &
+    call get_param(param_file, mdl, "USE_NET_FW_ADJUSTMENT_SIGN_BUG", &
+                 CS%use_net_FW_adjustment_sign_bug, &
+                   "If true, use the wrong sign for the adjustment to\n"//&
+                   "the net fresh-water.", default=.false.)
   call get_param(param_file, mdl, "ADJUST_NET_FRESH_WATER_BY_SCALING", &
                  CS%adjust_net_fresh_water_by_scaling, &
                  "If true, adjustments to net fresh water to achieve zero net are\n"//&
@@ -1312,27 +1336,30 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
   outunit = stdout()
 
   write(outunit,*) "BEGIN CHECKSUM(ice_ocean_boundary_type):: ", id, timestep
-  write(outunit,100) 'iobt%u_flux         ', mpp_chksum( iobt%u_flux         )
-  write(outunit,100) 'iobt%v_flux         ', mpp_chksum( iobt%v_flux         )
-  write(outunit,100) 'iobt%t_flux         ', mpp_chksum( iobt%t_flux         )
-  write(outunit,100) 'iobt%q_flux         ', mpp_chksum( iobt%q_flux         )
-  write(outunit,100) 'iobt%salt_flux      ', mpp_chksum( iobt%salt_flux      )
-  write(outunit,100) 'iobt%lw_flux        ', mpp_chksum( iobt%lw_flux        )
-  write(outunit,100) 'iobt%sw_flux_vis_dir', mpp_chksum( iobt%sw_flux_vis_dir)
-  write(outunit,100) 'iobt%sw_flux_vis_dif', mpp_chksum( iobt%sw_flux_vis_dif)
-  write(outunit,100) 'iobt%sw_flux_nir_dir', mpp_chksum( iobt%sw_flux_nir_dir)
-  write(outunit,100) 'iobt%sw_flux_nir_dif', mpp_chksum( iobt%sw_flux_nir_dif)
-  write(outunit,100) 'iobt%lprec          ', mpp_chksum( iobt%lprec          )
-  write(outunit,100) 'iobt%fprec          ', mpp_chksum( iobt%fprec          )
-  write(outunit,100) 'iobt%runoff         ', mpp_chksum( iobt%runoff         )
-  write(outunit,100) 'iobt%calving        ', mpp_chksum( iobt%calving        )
-  write(outunit,100) 'iobt%p              ', mpp_chksum( iobt%p              )
+  write(outunit,100) 'iobt%u_flux          ', mpp_chksum( iobt%u_flux          )
+  write(outunit,100) 'iobt%v_flux          ', mpp_chksum( iobt%v_flux          )
+  write(outunit,100) 'iobt%t_flux          ', mpp_chksum( iobt%t_flux          )
+  write(outunit,100) 'iobt%seaice_melt_heat', mpp_chksum( iobt%seaice_melt_heat)
+  write(outunit,100) 'iobt%seaice_melt     ', mpp_chksum( iobt%seaice_melt     )
+  write(outunit,100) 'iobt%q_flux          ', mpp_chksum( iobt%q_flux          )
+  write(outunit,100) 'iobt%rofl_flux       ', mpp_chksum( iobt%rofl_flux       )
+  write(outunit,100) 'iobt%rofi_flux       ', mpp_chksum( iobt%rofi_flux       )
+  write(outunit,100) 'iobt%salt_flux       ', mpp_chksum( iobt%salt_flux       )
+  write(outunit,100) 'iobt%lw_flux         ', mpp_chksum( iobt%lw_flux         )
+  write(outunit,100) 'iobt%sw_flux_vis_dir ', mpp_chksum( iobt%sw_flux_vis_dir )
+  write(outunit,100) 'iobt%sw_flux_vis_dif ', mpp_chksum( iobt%sw_flux_vis_dif )
+  write(outunit,100) 'iobt%sw_flux_nir_dir ', mpp_chksum( iobt%sw_flux_nir_dir )
+  write(outunit,100) 'iobt%sw_flux_nir_dif ', mpp_chksum( iobt%sw_flux_nir_dif )
+  write(outunit,100) 'iobt%lprec           ', mpp_chksum( iobt%lprec           )
+  write(outunit,100) 'iobt%fprec           ', mpp_chksum( iobt%fprec           )
+  write(outunit,100) 'iobt%calving         ', mpp_chksum( iobt%calving         )
+  write(outunit,100) 'iobt%p               ', mpp_chksum( iobt%p               )
   if (associated(iobt%ustar_berg)) &
-    write(outunit,100) 'iobt%ustar_berg     ', mpp_chksum( iobt%ustar_berg )
+    write(outunit,100) 'iobt%ustar_berg    ', mpp_chksum( iobt%ustar_berg      )
   if (associated(iobt%area_berg)) &
-    write(outunit,100) 'iobt%area_berg      ', mpp_chksum( iobt%area_berg  )
+    write(outunit,100) 'iobt%area_berg     ', mpp_chksum( iobt%area_berg       )
   if (associated(iobt%mass_berg)) &
-    write(outunit,100) 'iobt%mass_berg      ', mpp_chksum( iobt%mass_berg  )
+    write(outunit,100) 'iobt%mass_berg     ', mpp_chksum( iobt%mass_berg       )
 100 FORMAT("   CHECKSUM::",A20," = ",Z20)
 
   call coupler_type_write_chksums(iobt%fluxes, outunit, 'iobt%')
diff --git a/config_src/mct_driver/ocn_cap_methods.F90 b/config_src/mct_driver/ocn_cap_methods.F90
index a04e3af4aa..a8ac5310c7 100644
--- a/config_src/mct_driver/ocn_cap_methods.F90
+++ b/config_src/mct_driver/ocn_cap_methods.F90
@@ -33,7 +33,7 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
   real(kind=8), optional        , intent(in)    :: c1, c2, c3, c4     !< Coeffs. used in the shortwave decomposition
 
   ! Local variables
-  integer         :: i, j, ig, jg, isc, iec, jsc, jec  ! Grid indices
+  integer         :: i, j, isc, iec, jsc, jec  ! Grid indices
   integer         :: k
   integer         :: day, secs, rc
   type(ESMF_time) :: currTime
@@ -44,16 +44,17 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
 
   k = 0
   do j = jsc, jec
-    jg = j + grid%jsc - jsc
     do i = isc, iec
-      ig = i + grid%jsc - isc
       k = k + 1 ! Increment position within gindex
 
+      ! rotate taux and tauy from true zonal/meridional to local coordinates
       ! taux
-      ice_ocean_boundary%u_flux(i,j) = x2o(ind%x2o_Foxx_taux,k)
+      ice_ocean_boundary%u_flux(i,j) = GRID%cos_rot(i,j) * x2o(ind%x2o_Foxx_taux,k) &
+                                      + GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_tauy,k)
 
       ! tauy
-      ice_ocean_boundary%v_flux(i,j) = x2o(ind%x2o_Foxx_tauy,k)
+      ice_ocean_boundary%v_flux(i,j) = GRID%cos_rot(i,j) * x2o(ind%x2o_Foxx_tauy,k) &
+                                      - GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_taux,k)
 
       ! liquid precipitation (rain)
       ice_ocean_boundary%lprec(i,j) = x2o(ind%x2o_Faxa_rain,k)
@@ -65,25 +66,31 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
       ice_ocean_boundary%lw_flux(i,j) = (x2o(ind%x2o_Faxa_lwdn,k) + x2o(ind%x2o_Foxx_lwup,k))
 
       ! specific humitidy flux
-      ice_ocean_boundary%q_flux(i,j) = x2o(ind%x2o_Foxx_evap,k) !???TODO: should this be a minus sign
+      ice_ocean_boundary%q_flux(i,j) = x2o(ind%x2o_Foxx_evap,k)
 
       ! sensible heat flux (W/m2)
-      ice_ocean_boundary%t_flux(i,j) = x2o(ind%x2o_Foxx_sen,k)  !???TODO: should this be a minus sign
+      ice_ocean_boundary%t_flux(i,j) = x2o(ind%x2o_Foxx_sen,k)
 
       ! latent heat flux (W/m^2)
-      ice_ocean_boundary%latent_flux(i,j) = x2o(ind%x2o_Foxx_lat,k) !???TODO: should this be a minus sign
+      ice_ocean_boundary%latent_flux(i,j) = x2o(ind%x2o_Foxx_lat,k)
+
+      ! snow&ice melt heat flux  (W/m^2)
+      ice_ocean_boundary%seaice_melt_heat(i,j) = x2o(ind%x2o_Fioi_melth,k)
+
+      ! water flux from snow&ice melt (kg/m2/s)
+      ice_ocean_boundary%seaice_melt(i,j) = x2o(ind%x2o_Fioi_meltw,k)
 
       ! liquid runoff
-      ice_ocean_boundary%rofl_flux(i,j) = x2o(ind%x2o_Foxx_rofl,k) * GRID%mask2dT(ig,jg)
+      ice_ocean_boundary%rofl_flux(i,j) = x2o(ind%x2o_Foxx_rofl,k) * GRID%mask2dT(i,j)
 
       ! ice runoff
-      ice_ocean_boundary%rofi_flux(i,j) = x2o(ind%x2o_Foxx_rofi,k) * GRID%mask2dT(ig,jg)
+      ice_ocean_boundary%rofi_flux(i,j) = x2o(ind%x2o_Foxx_rofi,k) * GRID%mask2dT(i,j)
 
       ! surface pressure
-      ice_ocean_boundary%p(i,j) = x2o(ind%x2o_Sa_pslv,k) * GRID%mask2dT(ig,jg)
+      ice_ocean_boundary%p(i,j) = x2o(ind%x2o_Sa_pslv,k) * GRID%mask2dT(i,j)
 
-      ! salt flux
-      ice_ocean_boundary%salt_flux(i,j) = x2o(ind%x2o_Fioi_salt,k) * GRID%mask2dT(ig,jg)
+      ! salt flux (minus sign needed here -GMM)
+      ice_ocean_boundary%salt_flux(i,j) = -x2o(ind%x2o_Fioi_salt,k) * GRID%mask2dT(i,j)
 
       ! 1) visible, direct shortwave  (W/m2)
       ! 2) visible, diffuse shortwave (W/m2)
@@ -91,15 +98,15 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
       ! 4) near-IR, diffuse shortwave (W/m2)
       if (present(c1) .and. present(c2) .and. present(c3) .and. present(c4)) then
         ! Use runtime coefficients to decompose net short-wave heat flux into 4 components
-        ice_ocean_boundary%sw_flux_vis_dir(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c1 * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_vis_dif(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c2 * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_nir_dir(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c3 * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_nir_dif(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c4 * GRID%mask2dT(ig,jg)
+        ice_ocean_boundary%sw_flux_vis_dir(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c1 * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_vis_dif(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c2 * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_nir_dir(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c3 * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_nir_dif(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c4 * GRID%mask2dT(i,j)
       else
-        ice_ocean_boundary%sw_flux_vis_dir(i,j) = x2o(ind%x2o_Faxa_swvdr,k) * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_vis_dif(i,j) = x2o(ind%x2o_Faxa_swvdf,k) * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_nir_dir(i,j) = x2o(ind%x2o_Faxa_swndr,k) * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_nir_dif(i,j) = x2o(ind%x2o_Faxa_swndf,k) * GRID%mask2dT(ig,jg)
+        ice_ocean_boundary%sw_flux_vis_dir(i,j) = x2o(ind%x2o_Faxa_swvdr,k) * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_vis_dif(i,j) = x2o(ind%x2o_Faxa_swvdf,k) * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_nir_dir(i,j) = x2o(ind%x2o_Faxa_swndr,k) * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_nir_dif(i,j) = x2o(ind%x2o_Faxa_swndf,k) * GRID%mask2dT(i,j)
       endif
     enddo
   enddo
@@ -110,12 +117,16 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
 
     do j = GRID%jsc, GRID%jec
       do i = GRID%isc, GRID%iec
-        write(logunit,F01)'import: day, secs, j, i, u_flux          = ',day,secs,j,i,ice_ocean_boundary%u_flux(i,j)
-        write(logunit,F01)'import: day, secs, j, i, v_flux          = ',day,secs,j,i,ice_ocean_boundary%v_flux(i,j)
-        write(logunit,F01)'import: day, secs, j, i, lprec           = ',day,secs,j,i,ice_ocean_boundary%lprec(i,j)
-        write(logunit,F01)'import: day, secs, j, i, lwrad           = ',day,secs,j,i,ice_ocean_boundary%lw_flux(i,j)
-        write(logunit,F01)'import: day, secs, j, i, q_flux          = ',day,secs,j,i,ice_ocean_boundary%q_flux(i,j)
-        write(logunit,F01)'import: day, secs, j, i, t_flux          = ',day,secs,j,i,ice_ocean_boundary%t_flux(i,j)
+        write(logunit,F01)'import: day, secs, j, i, u_flux           = ',day,secs,j,i,ice_ocean_boundary%u_flux(i,j)
+        write(logunit,F01)'import: day, secs, j, i, v_flux           = ',day,secs,j,i,ice_ocean_boundary%v_flux(i,j)
+        write(logunit,F01)'import: day, secs, j, i, lprec            = ',day,secs,j,i,ice_ocean_boundary%lprec(i,j)
+        write(logunit,F01)'import: day, secs, j, i, lwrad            = ',day,secs,j,i,ice_ocean_boundary%lw_flux(i,j)
+        write(logunit,F01)'import: day, secs, j, i, q_flux           = ',day,secs,j,i,ice_ocean_boundary%q_flux(i,j)
+        write(logunit,F01)'import: day, secs, j, i, t_flux           = ',day,secs,j,i,ice_ocean_boundary%t_flux(i,j)
+        write(logunit,F01)'import: day, secs, j, i, seaice_melt_heat = ',&
+                          day,secs,j,i,ice_ocean_boundary%seaice_melt_heat(i,j)
+        write(logunit,F01)'import: day, secs, j, i, seaice_melt      = ',&
+                          day,secs,j,i,ice_ocean_boundary%seaice_melt(i,j)
         write(logunit,F01)'import: day, secs, j, i, latent_flux     = ',&
                           day,secs,j,i,ice_ocean_boundary%latent_flux(i,j)
         write(logunit,F01)'import: day, secs, j, i, runoff          = ',&
@@ -152,6 +163,8 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
 
   ! Local variables
   real, dimension(grid%isd:grid%ied,grid%jsd:grid%jed) :: ssh !< Local copy of sea_lev with updated halo
+  real, dimension(grid%isd:grid%ied,grid%jsd:grid%jed) :: sshx!< Zonal SSH gradient, local coordinate.
+  real, dimension(grid%isd:grid%ied,grid%jsd:grid%jed) :: sshy!< Meridional SSH gradient, local coordinate.
   integer :: i, j, n, ig, jg  !< Grid indices
   real    :: slp_L, slp_R, slp_C, slope, u_min, u_max
   real :: I_time_int  !< The inverse of coupling time interval [s-1].
@@ -174,8 +187,13 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
       ! surface temperature in Kelvin
       o2x(ind%o2x_So_t, n) = ocn_public%t_surf(ig,jg) * grid%mask2dT(i,j)
       o2x(ind%o2x_So_s, n) = ocn_public%s_surf(ig,jg) * grid%mask2dT(i,j)
-      o2x(ind%o2x_So_u, n) = ocn_public%u_surf(ig,jg) * grid%mask2dT(i,j)
-      o2x(ind%o2x_So_v, n) = ocn_public%v_surf(ig,jg) * grid%mask2dT(i,j)
+      ! rotate ocn current from local tripolar grid to true zonal/meridional (inverse transformation)
+      o2x(ind%o2x_So_u, n) = (grid%cos_rot(i,j) * ocn_public%u_surf(ig,jg) - &
+                              grid%sin_rot(i,j) * ocn_public%v_surf(ig,jg)) * grid%mask2dT(i,j)
+      o2x(ind%o2x_So_v, n) = (grid%cos_rot(i,j) * ocn_public%v_surf(ig,jg) + &
+                              grid%sin_rot(i,j) * ocn_public%u_surf(ig,jg)) * grid%mask2dT(i,j)
+
+      ! boundary layer depth (m)
       o2x(ind%o2x_So_bldepth, n) = ocn_public%OBLD(ig,jg) * grid%mask2dT(i,j)
       ! ocean melt and freeze potential (o2x_Fioo_q), W m-2
       if (ocn_public%frazil(ig,jg) > 0.0) then
@@ -197,9 +215,7 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
   call pass_var(ssh, grid%domain)
 
   ! d/dx ssh
-  n = 0
   do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
-    n = n+1
     ! This is a simple second-order difference
     ! o2x(ind%o2x_So_dhdx, n) = 0.5 * (ssh(i+1,j) - ssh(i-1,j)) * grid%IdxT(i,j) * grid%mask2dT(i,j)
     ! This is a PLM slope which might be less prone to the A-grid null mode
@@ -219,14 +235,12 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
       ! larger extreme values.
       slope = 0.0
     endif
-    o2x(ind%o2x_So_dhdx, n) = slope * grid%IdxT(i,j) * grid%mask2dT(i,j)
-    if (grid%mask2dT(i,j)==0.) o2x(ind%o2x_So_dhdx, n) = 0.0
+    sshx(i,j) = slope * grid%IdxT(i,j) * grid%mask2dT(i,j)
+    if (grid%mask2dT(i,j)==0.) sshx(i,j) = 0.0
   enddo; enddo
 
   ! d/dy ssh
-  n = 0
   do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
-    n = n+1
     ! This is a simple second-order difference
     ! o2x(ind%o2x_So_dhdy, n) = 0.5 * (ssh(i,j+1) - ssh(i,j-1)) * grid%IdyT(i,j) * grid%mask2dT(i,j)
     ! This is a PLM slope which might be less prone to the A-grid null mode
@@ -237,7 +251,6 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
     if (grid%mask2dCv(i,J+1)==0.) slp_R = 0.
 
     slp_C = 0.5 * (slp_L + slp_R)
-    !write(6,*)'slp_L, slp_R,i,j,slp_L*slp_R', slp_L, slp_R,i,j,slp_L*slp_R
     if ((slp_L * slp_R) > 0.0) then
       ! This limits the slope so that the edge values are bounded by the
       ! two cell averages spanning the edge.
@@ -249,8 +262,16 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
       ! larger extreme values.
       slope = 0.0
     endif
-    o2x(ind%o2x_So_dhdy, n) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j)
-    if (grid%mask2dT(i,j)==0.) o2x(ind%o2x_So_dhdy, n) = 0.0
+    sshy(i,j) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j)
+    if (grid%mask2dT(i,j)==0.) sshy(i,j) = 0.0
+  enddo; enddo
+
+  ! rotate ssh gradients from local coordinates to true zonal/meridional (inverse transformation)
+  n = 0
+  do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
+    n = n+1
+    o2x(ind%o2x_So_dhdx, n) = grid%cos_rot(i,j) * sshx(i,j) - grid%sin_rot(i,j) * sshy(i,j)
+    o2x(ind%o2x_So_dhdy, n) = grid%cos_rot(i,j) * sshy(i,j) + grid%sin_rot(i,j) * sshx(i,j)
   enddo; enddo
 
 end subroutine ocn_export
diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90
index c2e8423e4b..5ce89fc9f7 100644
--- a/config_src/mct_driver/ocn_comp_mct.F90
+++ b/config_src/mct_driver/ocn_comp_mct.F90
@@ -773,8 +773,6 @@ end subroutine ocean_model_init_sfc
 !! mi, mass of ice (kg/m2)
 !!
 !! Variables in the coupler that are **NOT** used in MOM6 (i.e., no corresponding field in fluxes):
-!! x2o_Fioi_melth, heat flux from snow & ice melt (W/m2)
-!! x2o_Fioi_meltw, snow melt flux (kg/m2/s)
 !! x2o_Si_ifrac, fractional ice wrt ocean
 !! x2o_So_duu10n, 10m wind speed squared (m^2/s^2)
 !! x2o_Sa_co2prog, bottom atm level prognostic CO2
diff --git a/config_src/mct_driver/ocn_cpl_indices.F90 b/config_src/mct_driver/ocn_cpl_indices.F90
index 645b358ec1..a701083c0c 100644
--- a/config_src/mct_driver/ocn_cpl_indices.F90
+++ b/config_src/mct_driver/ocn_cpl_indices.F90
@@ -41,7 +41,7 @@ module ocn_cpl_indices
     integer :: x2o_Faxa_swndr    !< near-IR, direct shortwave (W/m2)
     integer :: x2o_Faxa_swndf    !< near-IR, direct shortwave (W/m2)
     integer :: x2o_Fioi_melth    !< Heat flux from snow & ice melt (W/m2)
-    integer :: x2o_Fioi_meltw    !< Snow melt flux (kg/m2/s)
+    integer :: x2o_Fioi_meltw    !< Water flux from sea ice and snow melt (kg/m2/s)
     integer :: x2o_Fioi_bcpho    !< Black Carbon hydrophobic release from sea ice component
     integer :: x2o_Fioi_bcphi    !< Black Carbon hydrophilic release from sea ice component
     integer :: x2o_Fioi_flxdst   !< Dust release from sea ice component
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index 4cc413926d..465cdf2c28 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -69,9 +69,10 @@ module MOM_forcing_type
 
   ! turbulent heat fluxes into the ocean [W m-2]
   real, pointer, dimension(:,:) :: &
-    latent     => NULL(), & !< latent [W m-2]   (typically < 0)
-    sens       => NULL(), & !< sensible [W m-2] (typically negative)
-    heat_added => NULL()    !< additional heat flux from SST restoring or flux adjustments [W m-2]
+    latent           => NULL(), & !< latent [W m-2]   (typically < 0)
+    sens             => NULL(), & !< sensible [W m-2] (typically negative)
+    seaice_melt_heat => NULL(), & !< sea ice and snow melt or formation [W m-2] (typically negative)
+    heat_added       => NULL()    !< additional heat flux from SST restoring or flux adjustments [W m-2]
 
   ! components of latent heat fluxes used for diagnostic purposes
   real, pointer, dimension(:,:) :: &
@@ -87,7 +88,7 @@ module MOM_forcing_type
     vprec       => NULL(), & !< virtual liquid precip associated w/ SSS restoring [kg m-2 s-1]
     lrunoff     => NULL(), & !< liquid river runoff entering ocean [kg m-2 s-1]
     frunoff     => NULL(), & !< frozen river runoff (calving) entering ocean [kg m-2 s-1]
-    seaice_melt => NULL(), & !< seaice melt (positive) or formation (negative) [kg m-2 s-1]
+    seaice_melt => NULL(), & !< snow/seaice melt (positive) or formation (negative) [kg m-2 s-1]
     netMassIn   => NULL(), & !< Sum of water mass flux out of the ocean [kg m-2 s-1]
     netMassOut  => NULL(), & !< Net water mass flux into of the ocean [kg m-2 s-1]
     netSalt     => NULL()    !< Net salt entering the ocean [kgSalt m-2 s-1]
@@ -96,11 +97,11 @@ module MOM_forcing_type
   real, pointer, dimension(:,:) :: &
     heat_content_cond    => NULL(), & !< heat content associated with condensating water [W m-2]
     heat_content_lprec   => NULL(), & !< heat content associated with liquid >0 precip   [W m-2] (diagnostic)
+    heat_content_icemelt => NULL(), & !< heat content associated with snow/seaice melt/formation [W/m^2]
     heat_content_fprec   => NULL(), & !< heat content associated with frozen precip      [W m-2]
     heat_content_vprec   => NULL(), & !< heat content associated with virtual >0 precip  [W m-2]
     heat_content_lrunoff => NULL(), & !< heat content associated with liquid runoff      [W m-2]
     heat_content_frunoff => NULL(), & !< heat content associated with frozen runoff      [W m-2]
-    heat_content_icemelt => NULL(), & !< heat content associated with liquid sea ice     [W m-2]
     heat_content_massout => NULL(), & !< heat content associated with mass leaving ocean [W m-2]
     heat_content_massin  => NULL()    !< heat content associated with mass entering ocean [W m-2]
 
@@ -269,6 +270,7 @@ module MOM_forcing_type
   integer :: id_heat_content_vprec  = -1, id_heat_content_massout  = -1
   integer :: id_heat_added          = -1, id_heat_content_massin   = -1
   integer :: id_hfrainds            = -1, id_hfrunoffds            = -1
+  integer :: id_seaice_melt_heat    = -1, id_heat_content_icemelt  = -1
 
   ! global area integrated heat flux diagnostic handles
   integer :: id_total_net_heat_coupler    = -1, id_total_net_heat_surface      = -1
@@ -281,6 +283,7 @@ module MOM_forcing_type
   integer :: id_total_heat_content_cond   = -1, id_total_heat_content_surfwater= -1
   integer :: id_total_heat_content_vprec  = -1, id_total_heat_content_massout  = -1
   integer :: id_total_heat_added          = -1, id_total_heat_content_massin   = -1
+  integer :: id_total_seaice_melt_heat    = -1, id_total_heat_content_icemelt  = -1
 
   ! global area averaged heat flux diagnostic handles
   integer :: id_net_heat_coupler_ga = -1, id_net_heat_surface_ga = -1
@@ -501,20 +504,22 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
     endif
 
     ! net volume/mass of liquid and solid passing through surface boundary fluxes
-    netMassInOut(i) = dt * (scale * ((((( fluxes%lprec(i,j)      &
-                                        + fluxes%fprec(i,j)   )  &
-                                        + fluxes%evap(i,j)    )  &
-                                        + fluxes%lrunoff(i,j) )  &
-                                        + fluxes%vprec(i,j)   )  &
-                                        + fluxes%frunoff(i,j) ) )
+    netMassInOut(i) = dt * (scale * (((((( fluxes%lprec(i,j)        &
+                                        + fluxes%fprec(i,j)      )  &
+                                        + fluxes%evap(i,j)       )  &
+                                        + fluxes%lrunoff(i,j)    )  &
+                                        + fluxes%vprec(i,j)      )  &
+                                        + fluxes%seaice_melt(i,j))  &
+                                        + fluxes%frunoff(i,j)   ))
 
     if (do_NMIOr) then  ! Repeat the above code w/ dt=1s for legacy reasons
-      netMassInOut_rate(i) = (scale * ((((( fluxes%lprec(i,j)      &
-                                        + fluxes%fprec(i,j)   )  &
-                                        + fluxes%evap(i,j)    )  &
-                                        + fluxes%lrunoff(i,j) )  &
-                                        + fluxes%vprec(i,j)   )  &
-                                        + fluxes%frunoff(i,j) ) )
+      netMassInOut_rate(i) = (scale * (((((( fluxes%lprec(i,j)      &
+                                        + fluxes%fprec(i,j)      )  &
+                                        + fluxes%evap(i,j)       )  &
+                                        + fluxes%lrunoff(i,j)    )  &
+                                        + fluxes%vprec(i,j)      )  &
+                                        + fluxes%seaice_melt(i,j))  &
+                                        + fluxes%frunoff(i,j)   ))
     endif
 
     ! smg:
@@ -545,6 +550,11 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
       netMassOut(i) = netMassOut(i) + fluxes%lprec(i,j)
     endif
 
+    ! seaice_melt < 0 means sea ice formation taking water from the ocean.
+    if (fluxes%seaice_melt(i,j) < 0.0) then
+      netMassOut(i) = netMassOut(i) + fluxes%seaice_melt(i,j)
+    endif
+
     ! vprec < 0 means virtual evaporation arising from surface salinity restoring,
     ! in which case heat_content_vprec is computed in MOM_diabatic_driver.F90.
     if (fluxes%vprec(i,j) < 0.0) then
@@ -559,11 +569,24 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
 
     ! surface heat fluxes from radiation and turbulent fluxes (K * H)
     ! (H=m for Bouss, H=kg/m2 for non-Bouss)
-    net_heat(i) = scale * dt * J_m2_to_H * &
-                  ( fluxes%sw(i,j) +  ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j)) )
-    !Repeats above code w/ dt=1. for legacy reason
-    if (do_NHR)  net_heat_rate(i) = scale * J_m2_to_H * &
-         ( fluxes%sw(i,j) +  ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j)) )
+
+    ! CIME provides heat flux from snow&ice melt (seaice_melt_heat), so this is added below
+    if (associated(fluxes%seaice_melt_heat)) then
+      net_heat(i) = scale * dt * J_m2_to_H * &
+                    ( fluxes%sw(i,j) +  ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j) + &
+                      fluxes%seaice_melt_heat(i,j)) )
+      !Repeats above code w/ dt=1. for legacy reason
+      if (do_NHR)  net_heat_rate(i) = scale * J_m2_to_H * &
+           ( fluxes%sw(i,j) +  ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j) + &
+             fluxes%seaice_melt_heat(i,j)))
+    else
+      net_heat(i) = scale * dt * J_m2_to_H * &
+                    ( fluxes%sw(i,j) +  ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j)) )
+      !Repeats above code w/ dt=1. for legacy reason
+      if (do_NHR)  net_heat_rate(i) = scale * J_m2_to_H * &
+           ( fluxes%sw(i,j) +  ((fluxes%lw(i,j) + fluxes%latent(i,j)) + fluxes%sens(i,j)) )
+    endif
+
     ! Add heat flux from surface damping (restoring) (K * H) or flux adjustments.
     if (associated(fluxes%heat_added)) then
        net_heat(i) = net_heat(i) + (scale * (dt * J_m2_to_H)) * fluxes%heat_added(i,j)
@@ -719,6 +742,15 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
         endif
       endif
 
+      ! Following lprec and fprec, water flux due to sea ice melt (seaice_melt) enters at SST - GMM
+      if (associated(fluxes%heat_content_icemelt)) then
+        if (fluxes%seaice_melt(i,j) > 0.0) then
+          fluxes%heat_content_icemelt(i,j) = fluxes%C_p*fluxes%seaice_melt(i,j)*T(i,1)
+        else
+          fluxes%heat_content_icemelt(i,j) = 0.0
+        endif
+      endif
+
       ! virtual precip associated with salinity restoring
       ! vprec > 0 means add water to ocean, assumed to be at SST
       ! vprec < 0 means remove water from ocean; set heat_content_vprec in MOM_diabatic_driver.F90
@@ -1012,6 +1044,8 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift)
     call hchksum(fluxes%vprec, mesg//" fluxes%vprec",G%HI,haloshift=hshift)
   if (associated(fluxes%seaice_melt)) &
     call hchksum(fluxes%seaice_melt, mesg//" fluxes%seaice_melt",G%HI,haloshift=hshift)
+  if (associated(fluxes%seaice_melt_heat)) &
+    call hchksum(fluxes%seaice_melt_heat, mesg//" fluxes%seaice_melt_heat",G%HI,haloshift=hshift)
   if (associated(fluxes%p_surf)) &
     call hchksum(fluxes%p_surf, mesg//" fluxes%p_surf",G%HI,haloshift=hshift)
   if (associated(fluxes%salt_flux)) &
@@ -1032,6 +1066,8 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift)
     call hchksum(fluxes%heat_content_lprec, mesg//" fluxes%heat_content_lprec",G%HI,haloshift=hshift)
   if (associated(fluxes%heat_content_fprec)) &
     call hchksum(fluxes%heat_content_fprec, mesg//" fluxes%heat_content_fprec",G%HI,haloshift=hshift)
+  if (associated(fluxes%heat_content_icemelt)) &
+    call hchksum(fluxes%heat_content_icemelt, mesg//" fluxes%heat_content_icemelt",G%HI,haloshift=hshift)
   if (associated(fluxes%heat_content_cond)) &
     call hchksum(fluxes%heat_content_cond, mesg//" fluxes%heat_content_cond",G%HI,haloshift=hshift)
   if (associated(fluxes%heat_content_massout)) &
@@ -1123,6 +1159,7 @@ subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg)
   call locMsg(fluxes%fprec,'fprec')
   call locMsg(fluxes%vprec,'vprec')
   call locMsg(fluxes%seaice_melt,'seaice_melt')
+  call locMsg(fluxes%seaice_melt_heat,'seaice_melt_heat')
   call locMsg(fluxes%p_surf,'p_surf')
   call locMsg(fluxes%salt_flux,'salt_flux')
   call locMsg(fluxes%TKE_tidal,'TKE_tidal')
@@ -1133,6 +1170,7 @@ subroutine forcing_SinglePointPrint(fluxes, G, i, j, mesg)
   call locMsg(fluxes%heat_content_frunoff,'heat_content_frunoff')
   call locMsg(fluxes%heat_content_lprec,'heat_content_lprec')
   call locMsg(fluxes%heat_content_fprec,'heat_content_fprec')
+  call locMsg(fluxes%heat_content_icemelt,'heat_content_icemelt')
   call locMsg(fluxes%heat_content_vprec,'heat_content_vprec')
   call locMsg(fluxes%heat_content_cond,'heat_content_cond')
   call locMsg(fluxes%heat_content_cond,'heat_content_massout')
@@ -1231,13 +1269,13 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
        cmor_long_name='Water Evaporation Flux Where Ice Free Ocean over Sea')
 
   ! smg: seaice_melt field requires updates to the sea ice model
-  !handles%id_seaice_melt = register_diag_field('ocean_model', 'seaice_melt',       &
-  !   diag%axesT1, Time, 'water flux to ocean from sea ice melt(> 0) or form(< 0)', &
-  !   'kg m-2 s-1',                                                                 &
-  !    standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics',     &
-  !    cmor_field_name='fsitherm',                                                  &
-  !    cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics',&
-  !    cmor_long_name='water flux to ocean from sea ice melt(> 0) or form(< 0)')
+  handles%id_seaice_melt = register_diag_field('ocean_model', 'seaice_melt',       &
+     diag%axesT1, Time, 'water flux to ocean from snow/sea ice melting(> 0) or formation(< 0)', &
+     'kg m-2 s-1',                                                                 &
+      standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics',     &
+      cmor_field_name='fsitherm',                                                  &
+      cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics',&
+      cmor_long_name='water flux to ocean from sea ice melt(> 0) or form(< 0)')
 
   handles%id_precip = register_diag_field('ocean_model', 'precip', diag%axesT1, Time, &
         'Liquid + frozen precipitation into ocean', 'kg m-2 s-1')
@@ -1299,12 +1337,12 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
       cmor_long_name='Evaporation Where Ice Free Ocean over Sea Area Integrated')
 
   ! seaice_melt field requires updates to the sea ice model
-  !handles%id_total_seaice_melt = register_scalar_field('ocean_model', 'total_seaice_melt', Time, diag, &
-  !    long_name='Area integrated sea ice melt (>0) or form (<0)', units='kg s-1',                      &
-  !    standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics_area_integrated',         &
-  !    cmor_field_name='total_fsitherm',                                                                &
-  !    cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics_area_integrated',    &
-  !    cmor_long_name='Water Melt/Form from Sea Ice Area Integrated')
+  handles%id_total_seaice_melt = register_scalar_field('ocean_model', 'total_icemelt', Time, diag, &
+      long_name='Area integrated sea ice melt (>0) or form (<0)', units='kg s-1',                      &
+      standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics_area_integrated',         &
+      cmor_field_name='total_fsitherm',                                                                &
+      cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics_area_integrated',    &
+      cmor_long_name='Water Melt/Form from Sea Ice Area Integrated')
 
   handles%id_total_precip = register_scalar_field('ocean_model', 'total_precip', Time, diag, &
       long_name='Area integrated liquid+frozen precip into ocean', units='kg s-1')
@@ -1404,6 +1442,10 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
         diag%axesT1,Time,'Heat content (relative to 0degC) of frozen prec entering ocean',&
         'W m-2')
 
+  handles%id_heat_content_icemelt = register_diag_field('ocean_model', 'heat_content_icemelt',&
+        diag%axesT1,Time,'Heat content (relative to 0degC) of water flux due to sea ice melting/freezing',&
+        'W m-2')
+
   handles%id_heat_content_vprec = register_diag_field('ocean_model', 'heat_content_vprec',   &
         diag%axesT1,Time,'Heat content (relative to 0degC) of virtual precip entering ocean',&
         'W m-2')
@@ -1434,14 +1476,16 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
         'W m-2')
 
   handles%id_net_heat_coupler = register_diag_field('ocean_model', 'net_heat_coupler',          &
-        diag%axesT1,Time,'Surface ocean heat flux from SW+LW+latent+sensible (via the coupler)',&
+        diag%axesT1,Time,'Surface ocean heat flux from SW+LW+latent+sensible+seaice_melt_heat (via the coupler)',&
         'W m-2')
 
-  handles%id_net_heat_surface = register_diag_field('ocean_model', 'net_heat_surface',diag%axesT1,  &
-        Time,'Surface ocean heat flux from SW+LW+lat+sens+mass transfer+frazil+restore or flux adjustments', 'W m-2',&
+  handles%id_net_heat_surface = register_diag_field('ocean_model', 'net_heat_surface',diag%axesT1, Time,  &
+        'Surface ocean heat flux from SW+LW+lat+sens+mass transfer+frazil+restore+seaice_melt_heat or '// &
+        'flux adjustments',&
+        'W m-2',&
         standard_name='surface_downward_heat_flux_in_sea_water', cmor_field_name='hfds',            &
         cmor_standard_name='surface_downward_heat_flux_in_sea_water',           &
-        cmor_long_name='Surface ocean heat flux from SW+LW+latent+sensible+masstransfer+frazil')
+        cmor_long_name='Surface ocean heat flux from SW+LW+latent+sensible+masstransfer+frazil+seaice_melt_heat')
 
   handles%id_sw = register_diag_field('ocean_model', 'SW', diag%axesT1, Time,  &
         'Shortwave radiation flux into ocean', 'W m-2',                        &
@@ -1494,6 +1538,13 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
         cmor_standard_name='surface_downward_sensible_heat_flux',                    &
         cmor_long_name='Surface Downward Sensible Heat Flux')
 
+  handles%id_seaice_melt_heat = register_diag_field('ocean_model', 'seaice_melt_heat', diag%axesT1, Time,&
+        'Heat flux into ocean due to snow and sea ice melt/freeze', 'W m-2',      &
+        standard_name='snow_ice_melt_heat_flux',                         &
+  !GMM TODO cmor_field_name='hfsso',                                                     &
+        cmor_standard_name='snow_ice_melt_heat_flux',                    &
+        cmor_long_name='Heat flux into ocean from snow and sea ice melt')
+
   handles%id_heat_added = register_diag_field('ocean_model', 'heat_added', diag%axesT1, Time, &
         'Flux Adjustment or restoring surface heat flux into ocean', 'W m-2')
 
@@ -1533,6 +1584,11 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
       long_name='Area integrated heat content (relative to 0C) of frozen precip',&
       units='W')
 
+  handles%id_total_heat_content_icemelt = register_scalar_field('ocean_model',     &
+      'total_heat_content_icemelt', Time, diag,long_name=                          &
+      'Area integrated heat content (relative to 0C) of water flux due sea ice melting/freezing', &
+      units='W')
+
   handles%id_total_heat_content_vprec = register_scalar_field('ocean_model',      &
       'total_heat_content_vprec', Time, diag,                                     &
       long_name='Area integrated heat content (relative to 0C) of virtual precip',&
@@ -1564,7 +1620,7 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
 
   handles%id_total_net_heat_coupler = register_scalar_field('ocean_model',                       &
       'total_net_heat_coupler', Time, diag,                                                      &
-      long_name='Area integrated surface heat flux from SW+LW+latent+sensible (via the coupler)',&
+      long_name='Area integrated surface heat flux from SW+LW+latent+sensible+seaice_melt_heat (via the coupler)',&
       units='W')
 
   handles%id_total_net_heat_surface = register_scalar_field('ocean_model',                      &
@@ -1645,18 +1701,22 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
       long_name='Area integrated surface heat flux from restoring and/or flux adjustment',   &
       units='W')
 
+  handles%id_total_seaice_melt_heat = register_scalar_field('ocean_model',&
+      'total_seaice_melt_heat', Time, diag,                               &
+      long_name='Area integrated surface heat flux from snow and sea ice melt',   &
+      units='W')
 
   !===============================================================
   ! area averaged surface heat fluxes
 
   handles%id_net_heat_coupler_ga = register_scalar_field('ocean_model',                       &
       'net_heat_coupler_ga', Time, diag,                                                      &
-      long_name='Area averaged surface heat flux from SW+LW+latent+sensible (via the coupler)',&
+      long_name='Area averaged surface heat flux from SW+LW+latent+sensible+seaice_melt_heat (via the coupler)',&
       units='W m-2')
 
   handles%id_net_heat_surface_ga = register_scalar_field('ocean_model',                       &
-      'net_heat_surface_ga', Time, diag,                                                      &
-      long_name='Area averaged surface heat flux from SW+LW+lat+sens+mass+frazil+restore or flux adjustments', &
+      'net_heat_surface_ga', Time, diag, long_name=                                                     &
+      'Area averaged surface heat flux from SW+LW+lat+sens+mass+frazil+restore+seaice_melt_heat or flux adjustments', &
       units='W m-2',                                                                          &
       cmor_field_name='ave_hfds',                                                             &
       cmor_standard_name='surface_downward_heat_flux_in_sea_water_area_averaged',             &
@@ -1855,8 +1915,7 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, dt, G, wt2, forces)
     fluxes%vprec(i,j) = wt1*fluxes%vprec(i,j) + wt2*flux_tmp%vprec(i,j)
     fluxes%lrunoff(i,j) = wt1*fluxes%lrunoff(i,j) + wt2*flux_tmp%lrunoff(i,j)
     fluxes%frunoff(i,j) = wt1*fluxes%frunoff(i,j) + wt2*flux_tmp%frunoff(i,j)
- ! ### ADD LATER fluxes%seaice_melt(i,j) = wt1*fluxes%seaice_melt(i,j) + wt2*flux_tmp%seaice_melt(i,j)
-
+    fluxes%seaice_melt(i,j) = wt1*fluxes%seaice_melt(i,j) + wt2*flux_tmp%seaice_melt(i,j)
     fluxes%sw(i,j) = wt1*fluxes%sw(i,j) + wt2*flux_tmp%sw(i,j)
     fluxes%sw_vis_dir(i,j) = wt1*fluxes%sw_vis_dir(i,j) + wt2*flux_tmp%sw_vis_dir(i,j)
     fluxes%sw_vis_dif(i,j) = wt1*fluxes%sw_vis_dif(i,j) + wt2*flux_tmp%sw_vis_dif(i,j)
@@ -1889,6 +1948,11 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, dt, G, wt2, forces)
       fluxes%heat_content_fprec(i,j) = wt1*fluxes%heat_content_fprec(i,j) + wt2*flux_tmp%heat_content_fprec(i,j)
     enddo ; enddo
   endif
+  if (associated(fluxes%heat_content_icemelt) .and. associated(flux_tmp%heat_content_icemelt)) then
+    do j=js,je ; do i=is,ie
+      fluxes%heat_content_icemelt(i,j) = wt1*fluxes%heat_content_icemelt(i,j) + wt2*flux_tmp%heat_content_icemelt(i,j)
+    enddo ; enddo
+  endif
   if (associated(fluxes%heat_content_vprec) .and. associated(flux_tmp%heat_content_vprec)) then
     do j=js,je ; do i=is,ie
       fluxes%heat_content_vprec(i,j) = wt1*fluxes%heat_content_vprec(i,j) + wt2*flux_tmp%heat_content_vprec(i,j)
@@ -2017,7 +2081,7 @@ subroutine set_derived_forcing_fields(forces, fluxes, G, US, Rho0)
 end subroutine set_derived_forcing_fields
 
 
-!> This subroutine calculates determines the net mass source to the ocean from
+!> This subroutine determines the net mass source to the ocean from
 !! a (thermodynamic) forcing type and stores it in a mech_forcing type.
 subroutine set_net_mass_forcing(fluxes, forces, G)
   type(forcing),           intent(in)    :: fluxes   !< A structure containing thermodynamic forcing fields
@@ -2059,6 +2123,9 @@ subroutine get_net_mass_forcing(fluxes, G, net_mass_src)
   if (associated(fluxes%evap)) then ; do j=js,je ; do i=is,ie
     net_mass_src(i,j) = net_mass_src(i,j) + fluxes%evap(i,j)
   enddo ; enddo ; endif
+  if (associated(fluxes%seaice_melt)) then ; do j=js,je ; do i=is,ie
+    net_mass_src(i,j) = net_mass_src(i,j) + fluxes%seaice_melt(i,j)
+  enddo ; enddo ; endif
 
 end subroutine get_net_mass_forcing
 
@@ -2156,6 +2223,7 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
         if (associated(fluxes%lrunoff))     res(i,j) = res(i,j)+fluxes%lrunoff(i,j)
         if (associated(fluxes%frunoff))     res(i,j) = res(i,j)+fluxes%frunoff(i,j)
         if (associated(fluxes%vprec))       res(i,j) = res(i,j)+fluxes%vprec(i,j)
+        if (associated(fluxes%seaice_melt)) res(i,j) = res(i,j)+fluxes%seaice_melt(i,j)
       enddo ; enddo
       if (handles%id_prcme > 0) call post_data(handles%id_prcme, res, diag)
       if (handles%id_total_prcme > 0) then
@@ -2171,9 +2239,10 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
     if (handles%id_net_massout > 0 .or. handles%id_total_net_massout > 0) then
       do j=js,je ; do i=is,ie
         res(i,j) = 0.0
-        if (fluxes%lprec(i,j) < 0.0) res(i,j) = res(i,j) + fluxes%lprec(i,j)
-        if (fluxes%vprec(i,j) < 0.0) res(i,j) = res(i,j) + fluxes%vprec(i,j)
-        if (fluxes%evap(i,j)  < 0.0) res(i,j) = res(i,j) + fluxes%evap(i,j)
+        if (fluxes%lprec(i,j)       < 0.0) res(i,j) = res(i,j) + fluxes%lprec(i,j)
+        if (fluxes%vprec(i,j)       < 0.0) res(i,j) = res(i,j) + fluxes%vprec(i,j)
+        if (fluxes%evap(i,j)        < 0.0) res(i,j) = res(i,j) + fluxes%evap(i,j)
+        if (fluxes%seaice_melt(i,j) < 0.0) res(i,j) = res(i,j) + fluxes%seaice_melt(i,j)
       enddo ; enddo
       if (handles%id_net_massout > 0) call post_data(handles%id_net_massout, res, diag)
       if (handles%id_total_net_massout > 0) then
@@ -2187,10 +2256,11 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
     if (handles%id_net_massin > 0 .or. handles%id_total_net_massin > 0) then
       do j=js,je ; do i=is,ie
         res(i,j) = fluxes%fprec(i,j) + fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)
-        if (fluxes%lprec(i,j) > 0.0) res(i,j) = res(i,j) + fluxes%lprec(i,j)
-        if (fluxes%vprec(i,j) > 0.0) res(i,j) = res(i,j) + fluxes%vprec(i,j)
+        if (fluxes%lprec(i,j)       > 0.0) res(i,j) = res(i,j) + fluxes%lprec(i,j)
+        if (fluxes%vprec(i,j)       > 0.0) res(i,j) = res(i,j) + fluxes%vprec(i,j)
         ! fluxes%cond is not needed because it is derived from %evap > 0
-        if (fluxes%evap(i,j)  > 0.0) res(i,j) = res(i,j) + fluxes%evap(i,j)
+        if (fluxes%evap(i,j)        > 0.0) res(i,j) = res(i,j) + fluxes%evap(i,j)
+        if (fluxes%seaice_melt(i,j) > 0.0) res(i,j) = res(i,j) + fluxes%seaice_melt(i,j)
       enddo ; enddo
       if (handles%id_net_massin > 0) call post_data(handles%id_net_massin, res, diag)
       if (handles%id_total_net_massin > 0) then
@@ -2279,6 +2349,14 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
       endif
     endif
 
+    if (associated(fluxes%seaice_melt)) then
+      if (handles%id_seaice_melt > 0) call post_data(handles%id_seaice_melt, fluxes%seaice_melt, diag)
+      if (handles%id_total_seaice_melt > 0) then
+        total_transport = global_area_integral(fluxes%seaice_melt,G)
+        call post_data(handles%id_total_seaice_melt, total_transport, diag)
+      endif
+    endif
+
     ! post diagnostics for boundary heat fluxes ====================================
 
     if ((handles%id_heat_content_lrunoff > 0) .and. associated(fluxes%heat_content_lrunoff))  &
@@ -2309,6 +2387,13 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
       call post_data(handles%id_total_heat_content_fprec, total_transport, diag)
     endif
 
+    if ((handles%id_heat_content_icemelt > 0) .and. associated(fluxes%heat_content_icemelt))      &
+      call post_data(handles%id_heat_content_icemelt, fluxes%heat_content_icemelt, diag)
+    if ((handles%id_total_heat_content_icemelt > 0) .and. associated(fluxes%heat_content_icemelt)) then
+      total_transport = global_area_integral(fluxes%heat_content_icemelt,G)
+      call post_data(handles%id_total_heat_content_icemelt, total_transport, diag)
+    endif
+
     if ((handles%id_heat_content_vprec > 0) .and. associated(fluxes%heat_content_vprec))      &
       call post_data(handles%id_heat_content_vprec, fluxes%heat_content_vprec, diag)
     if ((handles%id_total_heat_content_vprec > 0) .and. associated(fluxes%heat_content_vprec)) then
@@ -2341,10 +2426,11 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
         handles%id_net_heat_coupler_ga > 0. ) then
       do j=js,je ; do i=is,ie
       res(i,j) = 0.0
-      if (associated(fluxes%LW))         res(i,j) = res(i,j) + fluxes%LW(i,j)
-      if (associated(fluxes%latent))     res(i,j) = res(i,j) + fluxes%latent(i,j)
-      if (associated(fluxes%sens))       res(i,j) = res(i,j) + fluxes%sens(i,j)
-      if (associated(fluxes%SW))         res(i,j) = res(i,j) + fluxes%SW(i,j)
+      if (associated(fluxes%LW))               res(i,j) = res(i,j) + fluxes%LW(i,j)
+      if (associated(fluxes%latent))           res(i,j) = res(i,j) + fluxes%latent(i,j)
+      if (associated(fluxes%sens))             res(i,j) = res(i,j) + fluxes%sens(i,j)
+      if (associated(fluxes%SW))               res(i,j) = res(i,j) + fluxes%SW(i,j)
+      if (associated(fluxes%seaice_melt_heat)) res(i,j) = res(i,j) + fluxes%seaice_melt_heat(i,j)
       enddo ; enddo
       if (handles%id_net_heat_coupler > 0) call post_data(handles%id_net_heat_coupler, res, diag)
       if (handles%id_total_net_heat_coupler > 0) then
@@ -2365,18 +2451,20 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
         if (associated(fluxes%latent))               res(i,j) = res(i,j) + fluxes%latent(i,j)
         if (associated(fluxes%sens))                 res(i,j) = res(i,j) + fluxes%sens(i,j)
         if (associated(fluxes%SW))                   res(i,j) = res(i,j) + fluxes%SW(i,j)
+        if (associated(fluxes%seaice_melt_heat))     res(i,j) = res(i,j) + fluxes%seaice_melt_heat(i,j)
         if (associated(sfc_state%frazil))            res(i,j) = res(i,j) + sfc_state%frazil(i,j) * I_dt
-      ! if (associated(sfc_state%TempXpme)) then
-      !    res(i,j) = res(i,j) + sfc_state%TempXpme(i,j) * fluxes%C_p * I_dt
-      ! else
-        if (associated(fluxes%heat_content_lrunoff)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j)
-        if (associated(fluxes%heat_content_frunoff)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff(i,j)
-        if (associated(fluxes%heat_content_lprec))   res(i,j) = res(i,j) + fluxes%heat_content_lprec(i,j)
-        if (associated(fluxes%heat_content_fprec))   res(i,j) = res(i,j) + fluxes%heat_content_fprec(i,j)
-        if (associated(fluxes%heat_content_vprec))   res(i,j) = res(i,j) + fluxes%heat_content_vprec(i,j)
-        if (associated(fluxes%heat_content_cond))    res(i,j) = res(i,j) + fluxes%heat_content_cond(i,j)
-        if (associated(fluxes%heat_content_massout)) res(i,j) = res(i,j) + fluxes%heat_content_massout(i,j)
-      ! endif
+        !if (associated(sfc_state%TempXpme)) then
+        !  res(i,j) = res(i,j) + sfc_state%TempXpme(i,j) * fluxes%C_p * I_dt
+        !else
+          if (associated(fluxes%heat_content_lrunoff)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j)
+          if (associated(fluxes%heat_content_frunoff)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff(i,j)
+          if (associated(fluxes%heat_content_lprec))   res(i,j) = res(i,j) + fluxes%heat_content_lprec(i,j)
+          if (associated(fluxes%heat_content_fprec))   res(i,j) = res(i,j) + fluxes%heat_content_fprec(i,j)
+          if (associated(fluxes%heat_content_icemelt)) res(i,j) = res(i,j) + fluxes%heat_content_icemelt(i,j)
+          if (associated(fluxes%heat_content_vprec))   res(i,j) = res(i,j) + fluxes%heat_content_vprec(i,j)
+          if (associated(fluxes%heat_content_cond))    res(i,j) = res(i,j) + fluxes%heat_content_cond(i,j)
+          if (associated(fluxes%heat_content_massout)) res(i,j) = res(i,j) + fluxes%heat_content_massout(i,j)
+        !endif
         if (associated(fluxes%heat_added))         res(i,j) = res(i,j) + fluxes%heat_added(i,j)
       enddo ; enddo
       if (handles%id_net_heat_surface > 0) call post_data(handles%id_net_heat_surface, res, diag)
@@ -2400,6 +2488,7 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
           if (associated(fluxes%heat_content_lrunoff)) res(i,j) = res(i,j) + fluxes%heat_content_lrunoff(i,j)
           if (associated(fluxes%heat_content_frunoff)) res(i,j) = res(i,j) + fluxes%heat_content_frunoff(i,j)
           if (associated(fluxes%heat_content_lprec))   res(i,j) = res(i,j) + fluxes%heat_content_lprec(i,j)
+          if (associated(fluxes%heat_content_icemelt)) res(i,j) = res(i,j) + fluxes%heat_content_icemelt(i,j)
           if (associated(fluxes%heat_content_fprec))   res(i,j) = res(i,j) + fluxes%heat_content_fprec(i,j)
           if (associated(fluxes%heat_content_vprec))   res(i,j) = res(i,j) + fluxes%heat_content_vprec(i,j)
           if (associated(fluxes%heat_content_cond))    res(i,j) = res(i,j) + fluxes%heat_content_cond(i,j)
@@ -2531,6 +2620,16 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
     if ((handles%id_sens > 0) .and. associated(fluxes%sens)) then
       call post_data(handles%id_sens, fluxes%sens, diag)
     endif
+
+    if ((handles%id_seaice_melt_heat > 0) .and. associated(fluxes%seaice_melt_heat)) then
+      call post_data(handles%id_seaice_melt_heat, fluxes%seaice_melt_heat, diag)
+    endif
+
+    if ((handles%id_total_seaice_melt_heat > 0) .and. associated(fluxes%seaice_melt_heat)) then
+      total_transport = global_area_integral(fluxes%seaice_melt_heat,G)
+      call post_data(handles%id_total_seaice_melt_heat, total_transport, diag)
+    endif
+
     if ((handles%id_total_sens > 0) .and. associated(fluxes%sens)) then
       total_transport = global_area_integral(fluxes%sens,G)
       call post_data(handles%id_total_sens, total_transport, diag)
@@ -2648,7 +2747,7 @@ subroutine allocate_forcing_type(G, fluxes, water, heat, ustar, press, shelf, ic
   call myAlloc(fluxes%netMassOut,isd,ied,jsd,jed, water)
   call myAlloc(fluxes%netMassIn,isd,ied,jsd,jed, water)
   call myAlloc(fluxes%netSalt,isd,ied,jsd,jed, water)
-
+  call myAlloc(fluxes%seaice_melt_heat,isd,ied,jsd,jed, heat)
   call myAlloc(fluxes%sw,isd,ied,jsd,jed, heat)
   call myAlloc(fluxes%lw,isd,ied,jsd,jed, heat)
   call myAlloc(fluxes%latent,isd,ied,jsd,jed, heat)
@@ -2661,6 +2760,7 @@ subroutine allocate_forcing_type(G, fluxes, water, heat, ustar, press, shelf, ic
 
   if (present(heat) .and. present(water)) then ; if (heat .and. water) then
     call myAlloc(fluxes%heat_content_cond,isd,ied,jsd,jed, .true.)
+    call myAlloc(fluxes%heat_content_icemelt,isd,ied,jsd,jed, .true.)
     call myAlloc(fluxes%heat_content_lprec,isd,ied,jsd,jed, .true.)
     call myAlloc(fluxes%heat_content_fprec,isd,ied,jsd,jed, .true.)
     call myAlloc(fluxes%heat_content_vprec,isd,ied,jsd,jed, .true.)
@@ -2743,6 +2843,7 @@ subroutine deallocate_forcing_type(fluxes)
   if (associated(fluxes%ustar_gustless))       deallocate(fluxes%ustar_gustless)
   if (associated(fluxes%buoy))                 deallocate(fluxes%buoy)
   if (associated(fluxes%sw))                   deallocate(fluxes%sw)
+  if (associated(fluxes%seaice_melt_heat))     deallocate(fluxes%seaice_melt_heat)
   if (associated(fluxes%sw_vis_dir))           deallocate(fluxes%sw_vis_dir)
   if (associated(fluxes%sw_vis_dif))           deallocate(fluxes%sw_vis_dif)
   if (associated(fluxes%sw_nir_dir))           deallocate(fluxes%sw_nir_dir)
@@ -2756,6 +2857,7 @@ subroutine deallocate_forcing_type(fluxes)
   if (associated(fluxes%heat_added))           deallocate(fluxes%heat_added)
   if (associated(fluxes%heat_content_lrunoff)) deallocate(fluxes%heat_content_lrunoff)
   if (associated(fluxes%heat_content_frunoff)) deallocate(fluxes%heat_content_frunoff)
+  if (associated(fluxes%heat_content_icemelt)) deallocate(fluxes%heat_content_icemelt)
   if (associated(fluxes%heat_content_lprec))   deallocate(fluxes%heat_content_lprec)
   if (associated(fluxes%heat_content_fprec))   deallocate(fluxes%heat_content_fprec)
   if (associated(fluxes%heat_content_cond))    deallocate(fluxes%heat_content_cond)
@@ -2905,7 +3007,7 @@ end subroutine deallocate_mech_forcing
 !!  * non-penetrative = non-downwelling shortwave; portion of SW
 !!                      totally absorbed in the k=1 cell.
 !!                      The non-penetrative SW is combined with
-!!                      LW+LAT+SENS in net_heat inside routine
+!!                      LW+LAT+SENS+seaice_melt_heat in net_heat inside routine
 !!                      extractFluxes1d. Notably, for many cases,
 !!                      non-penetrative SW = 0.
 !!  * penetrative     = that portion of shortwave penetrating below
diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90
index e6dcbf11bc..cfc74b47fc 100644
--- a/src/diagnostics/MOM_sum_output.F90
+++ b/src/diagnostics/MOM_sum_output.F90
@@ -956,6 +956,10 @@ subroutine accumulate_net_input(fluxes, sfc_state, dt, G, CS)
     endif
   endif
 
+  if (associated(fluxes%seaice_melt)) then ; do j=js,je ; do i=is,ie
+    FW_in(i,j) = FW_in(i,j) + dt * G%areaT(i,j) * fluxes%seaice_melt(i,j)
+  enddo ; enddo ; endif
+
   salt_in(:,:) = 0.0 ; heat_in(:,:) = 0.0
   if (CS%use_temperature) then
 
@@ -964,6 +968,10 @@ subroutine accumulate_net_input(fluxes, sfc_state, dt, G, CS)
              (fluxes%lw(i,j) + (fluxes%latent(i,j) + fluxes%sens(i,j))))
     enddo ; enddo ; endif
 
+    if (associated(fluxes%seaice_melt_heat)) then ; do j=js,je ; do i=is,ie
+       heat_in(i,j) = heat_in(i,j) + dt*G%areaT(i,j) * fluxes%seaice_melt_heat(i,j)
+    enddo ; enddo ; endif
+
     ! smg: new code
     ! include heat content from water transport across ocean surface
 !    if (associated(fluxes%heat_content_lprec)) then ; do j=js,je ; do i=is,ie
diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90
index 7cb566ddba..e8b4500bbc 100644
--- a/src/parameterizations/vertical/MOM_diabatic_aux.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90
@@ -1355,7 +1355,10 @@ subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgori
       call get_param(param_file, mdl, "RIVERMIX_DEPTH", CS%rivermix_depth, &
                  "The depth to which rivers are mixed if DO_RIVERMIX is \n"//&
                  "defined.", units="m", default=0.0, scale=US%m_to_Z)
-  else ; CS%do_rivermix = .false. ; CS%rivermix_depth = 0.0 ; endif
+  else
+    CS%do_rivermix = .false. ; CS%rivermix_depth = 0.0 ; CS%ignore_fluxes_over_land = .false.
+  endif
+
   if (GV%nkml == 0) then
     call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", CS%use_river_heat_content, &
                    "If true, use the fluxes%runoff_Hflx field to set the \n"//&