anchor.diff:          call longwave_upward(now%tsurf,sigm,now%lwu)
   
           ! residual energy from melting and refreezing
  -        qmr = 0.0_dp
  +        now%qmr = 0.0_dp
           where (now%mask == 2)         ! ice
  -            qmr = (now%melted_snow + now%melted_ice - now%refr)*rhow*clm
  +            now%qmr = (now%melted_snow + now%melted_ice - now%refr)*rhow*clm
           elsewhere                               ! land and sea
  -            qmr = (now%melted_snow - now%refr)*rhow*clm
  +            now%qmr = (now%melted_snow - now%refr)*rhow*clm
           end where
   
           ! surface energy balance of incoming and outgoing surface fluxes (W/m^2)
  -        qsb  = (1.0_dp-now%alb)*now%swd + now%lwd - now%lwu - now%shf - now%lhf - qmr
  +        qsb  = (1.0_dp-now%alb)*now%swd + now%lwd - now%lwu - now%shf - now%lhf - now%qmr
   
           ! update surface temperature according to surface energy balance
           if (.not. bnd%tsurf) then
               now%tsurf = now%tsurf + qsb*par%tsticsub/par%ceff 
  +            ! store residual energy for subfreezing tsurf over ice in qmr for later use in mass balance
  +            where (now%mask == 2 .and. now%tsurf > t0)
  +                now%qmr = (now%tsurf-t0)*par%ceff/par%tstic
  +                now%tsurf = t0
  +            else where
  +                now%qmr = 0.0_dp
  +            end where
           end if
           now%t2m   = now%t2m + (now%shf+now%lhf)*par%tsticsub/par%ceff  
   
  
positive.code: 
          end if
          
          !now%t2m = t2m_new 
          ! add remainder of melted snow/ice and refreezing to surface temperature
          !where (now%mask == 2)         ! ice
          !    now%tsurf = now%tsurf  - (now%melted_snow + now%melted_ice - now%refr)*par%tsticsub*rhow*clm/par%ceff
          !elsewhere                     ! land and sea
          !    now%tsurf = now%tsurf - (now%melted_snow - now%refr)*par%tsticsub*rhow*clm/par%ceff
          !end where
  
          ! write prognostic and diagnostic output
          now%subl = now%subl/rhow
  
  
negative.code:                 call albedo_denby(now%alb_snow,now%melt,par%alb_smax,par%alb_smin,par%mcrit)
              end if
              if (trim(par%alb_scheme) .eq. "isba") then
                  call albedo_isba(now%alb_snow,now%sf,now%melt,par%tstic,par%tstic,par%tau_a,par%tau_f,&
                                   par%w_crit,par%mcrit,par%alb_smin,par%alb_smax)
              end if
              where (now%mask == 2)
  