anchor.diff:          allocate(now%alb(npts))
           allocate(now%alb_snow(npts))
           allocate(now%melt(npts))
  +        allocate(now%melted_snow(npts))
  +        allocate(now%melted_ice(npts))
           allocate(now%refr(npts))
           allocate(now%smb(npts))
           allocate(now%smb_snow(npts))
  
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:                 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%tsticsub,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)
  