diff --git a/src/SIS_slow_thermo.F90 b/src/SIS_slow_thermo.F90
index 2b55c211..a4965e38 100644
--- a/src/SIS_slow_thermo.F90
+++ b/src/SIS_slow_thermo.F90
@@ -742,9 +742,24 @@ subroutine SIS2_thermodynamics(IST, dt_slow, CS, OSS, FIA, IOF, G, IG)
       qflx_res_ice(i,j) = -(LatHtFus*Rho_ice*Obs_h_ice(i,j)*Obs_cn_ice(i,j,2)-e2m_tot) / &
                            (86400.0*CS%ice_restore_timescale)
       if (qflx_res_ice(i,j) < 0.0) then
+        !There is less ice in model than Obs,
+        !so make some ice by increasing frazil heat
         FIA%frazil_left(i,j) = FIA%frazil_left(i,j) - qflx_res_ice(i,j)*dt_slow
+        !Note that ice should grow when frazil heat is positive
       elseif (qflx_res_ice(i,j) >  0.0) then
-        OSS%bheat(i,j) = OSS%bheat(i,j) + qflx_res_ice(i,j)
+        !There is more ice in model than Obs,
+        !so melt ice by increasing heat input to ice from ocean (bheat),
+        !        OSS%bheat(i,j) = OSS%bheat(i,j) + qflx_res_ice(i,j)
+        !Note that ice should melt when bheat increases.
+        !BUT, here it's too late for the bheat to have a negative feedback on the ice thickness
+        !since thickness is determined by the melting energies calculated in the fast ice
+        !module call ice_temp_SIS2() before this point.
+        !So, we should rather change the bottom melt energy directly here
+        !(as prescribed in ice_temp_SIS2) to have a restoring effect on the ice thickness
+        !later in the call ice_resize_SIS2() in this module.
+        do k=1,ncat
+          FIA%bmelt(i,j,k) = FIA%bmelt(i,j,k) + dt_slow*qflx_res_ice(i,j)
+        enddo
       endif
     enddo ; enddo
   endif
diff --git a/src/ice_model.F90 b/src/ice_model.F90
index 28cb6824..943015f2 100644
--- a/src/ice_model.F90
+++ b/src/ice_model.F90
@@ -1782,6 +1782,9 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow,
   logical :: redo_fast_update ! If true, recalculate the thermal updates from the fast
                               ! dynamics on the slowly evolving ice state, rather than
                               ! copying over the slow ice state to the fast ice state.
+  logical :: do_mask_restart  ! If true, apply the scaling and masks to mH_snow, mH_ice and part_size
+                              ! after a restart. However this may cause answers to diverge
+                              ! after a restart.Provide a switch to turn this option off.
   logical :: Verona
   logical :: Concurrent
   logical :: read_aux_restart
@@ -1996,6 +1999,9 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow,
                    "the sea ice that are handled by the fast ice PEs.", &
                    default="ice_model_fast.res.nc")
   endif
+  call get_param(param_file, mdl, "APPLY_MASKS_AFTER_RESTART", do_mask_restart, &
+                 "If true, applies masks to mH_ice,mH_snow and part_size after a restart.",&
+                  default=.true.)
 
   call get_param(param_file, mdl, "MASSLESS_ICE_ENTH", massless_ice_enth, &
                  "The ice enthalpy fill value for massless categories.", &
@@ -2472,10 +2478,21 @@ subroutine ice_model_init(Ice, Time_Init, Time, Time_step_fast, Time_step_slow,
       ! Deal with any ice masses or thicknesses over land, and rescale to
       ! account for differences between the current thickness units and whatever
       ! thickness units were in the input restart file.
-      do k=1,CatIce
-        sIST%mH_snow(:,:,k) = sIST%mH_snow(:,:,k) * H_rescale_snow * sG%mask2dT(:,:)
-        sIST%mH_ice(:,:,k) = sIST%mH_ice(:,:,k) * H_rescale_ice * sG%mask2dT(:,:)
-      enddo
+      ! However, in some model this causes answers to change after a restart
+      ! because these state variables are changed only after a restart and
+      ! not at each timestep. Provide a switch to turn this option off.
+      if(do_mask_restart) then
+        do k=1,CatIce
+          sIST%mH_snow(:,:,k) = sIST%mH_snow(:,:,k) * H_rescale_snow * sG%mask2dT(:,:)
+          sIST%mH_ice(:,:,k) = sIST%mH_ice(:,:,k) * H_rescale_ice * sG%mask2dT(:,:)
+          sIST%part_size(:,:,k) = sIST%part_size(:,:,k) * sG%mask2dT(:,:)
+        enddo
+        ! Since we masked out the part_size on land we should set
+        ! part_size(:,:,0) = 1. on land to satisfy the summation check
+        do j=jsc,jec ; do i=isc,iec
+          if (sG%mask2dT(i,j) < 0.5) sIST%part_size(i,j,0) = 1.
+        enddo ; enddo
+      endif
 
       if (sIG%ocean_part_min > 0.0) then ; do j=jsc,jec ; do i=isc,iec
         sIST%part_size(i,j,0) = max(sIST%part_size(i,j,0), sIG%ocean_part_min)