Fix the ice restoring issues

- This update addresses two independent issues with the ice restoring mechanism:
  1. SM4p model blows up with too much ice when it tries to do ice restoring to climatology.
     This only happens when SIS1 is replaced by SIS2.
     Upon investigation I found that the algorithm in SIS_slow_thermo.F90 has an update ordering issue.
     At the point when qflx_res_ice is calculated it updates the heat input to ice through bheat.
     But this does not to have an immediate direct effect on the ice thickness
     and its effect is delayed till the next timestep through other fluxes (flux_sh_ocn_top).
     This leads to a runaway increasing of the ice thickness (wrong feedback).
     To get around this we can increase bmelt instead of bheat which will decrease the ice thickness
     subsequently in the same timestep , in the same module through call ice_resize_SIS2().
     Standalone tests shows this has the correct restoring effect on the ice.

  2. SM4p has a restart issue, the answers of straight runs and restarted runs do not match (1x8days != 2x4days)
     I tracked this to a bit of code in ice_model.F90 that applies land mask to mH_ice.
     We subsequently tried applying lad mask to part_size (PR !95). It fixed a NaN issue but did not
     fix the restart issue.
     Applying masks to state variables such as mH_ice is leading to the restart issue because these
     state variables are being modified only after restarts and not in straight (no-restart) runs.
     I provided a switch to avoid this masking if the user needs to. With it SM4p reproduces.

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

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,16 +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(:,:)
-        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
+      ! 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)