anchor.diff:          double precision    :: tmin      !< minimum temperature for which albedo decline becomes effective ("slater") [K]
           double precision    :: tmax      !< maximum temperature for which albedo decline becomes effective ("slater") [K]
           double precision    :: tstic     !< time step [s]
  +        double precision    :: tsticsub  !< sub-time step [s]
           double precision    :: shf_enh   !< sensible heat flux enhancement factor for positive u*(tsurf-t2m) [no unit]
           double precision    :: lhf_enh   !< latent heat flux enhancement factor over land (mask==1) [no unit]
           double precision    :: tau_a     !< dry albedo decline for "isba" albedo scheme [1/day]
  
positive.code:         end if
          
          ! 5) surface mass balance
          now%smb_snow = now%sf - now%subl/rhow - melt_snow
  
          where (now%mask == 0)
              now%hsnow = 0.0_dp
          else where
              ! update snow height
              now%hsnow = dmax1(0.0_dp, now%hsnow + now%smb_snow*par%tstic)
          end where
          
          ! 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 - melt_ice + now%refr   ! Use to force ice sheet model
          now%hice    = now%hice + now%smb_ice*par%tstic ! update new ice budget: remove or add ice
  
          if (.not. bnd%smb) then
              where (now%mask == 2)
                  now%smb = now%smb_snow + now%smb_ice - snow_to_ice/par%tstic
              elsewhere
                  now%smb = now%smb_snow + dmax1(0.0_dp,now%smb_ice - snow_to_ice/par%tstic)
              end where
          end if
  
          ! reset temperature according to melt/refreezing
          if (.not. bnd%tsurf) then
              where (now%mask == 2)         ! ice
                  now%tsurf = tsurf_new - (melt_snow + melt_ice - now%refr)*par%tstic*rhow*clm/par%ceff
              elsewhere                               ! land and sea
                  now%tsurf = tsurf_new - (melt_snow - now%refr)*par%tstic*rhow*clm/par%ceff
              end where
          end if
  
  
negative.code:         deallocate(now%alb)
          deallocate(now%alb_snow)
          deallocate(now%melt)
          deallocate(now%melted_ice)
          deallocate(now%melted_snow)
          deallocate(now%refr)
          deallocate(now%smb)
          deallocate(now%smb_snow)
  