anchor.diff:          !> 10. Relax snow height to maximum (eg, 5 m)
           snow_to_ice   = dmax1(0.d0,now%hsnow-hsmax)
           now%hsnow   = now%hsnow - snow_to_ice
  -        now%smb_ice = snow_to_ice/par%tstic - now%melted_ice + refrozen_rain  ! Use to force ice sheet model
  +        now%smb_ice = snow_to_ice/par%tstic - now%melted_ice! + refrozen_rain  ! Use to force ice sheet model
           now%hice    = now%hice + now%smb_ice*par%tstic ! update new ice budget: remove or add ice
   
           !> 11. Total surface mass balance
  
positive.code:         end where
  
          now%subl = now%subl/rhow
  
      end subroutine mass_balance
  
  
negative.code:                 !f_rz = (1.0_dp+tanh(now%hsnow-par%rcrit))/2.0_dp
              !f_rz = dmin1(1.0_dp,now%hsnow/par%rcrit)
              !end where
              f_rz = 1.0_dp! - exp(-now%hsnow/par%rcrit)
              ! potential refreezing
              now%refr = qcold/(rhow*clm)
              refrozen_rain = dmin1(now%refr,now%rf)
  