anchor.diff:          double precision, allocatable, dimension(:) :: smb_snow !< surface mass balance of snow [m/s]
           double precision, allocatable, dimension(:) :: smb_ice  !< surface mass balance of ice [m/s]
           double precision, allocatable, dimension(:) :: runoff   !< potential surface runoff [m/s]
  +        double precision, allocatable, dimension(:) :: qmr      !< heat flux from melting/refreezing [W/m2]
           ! Forcing variables
           double precision, allocatable, dimension(:) :: sf       !< snow fall [m/s]
           double precision, allocatable, dimension(:) :: rf       !< rain fall [m/s]
  
positive.code:             now%hsnow = 0.0_dp
          else where
              ! update snow height
              now%hsnow = dmax1(0.0_dp, now%hsnow + now%smb_snow*par%tsticsub)
          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%tsticsub - now%melted_ice + now%refr   ! Use to force ice sheet model
          now%hice    = now%hice + now%smb_ice*par%tsticsub ! 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%tsticsub
              elsewhere
                  now%smb = now%smb_snow + dmax1(0.0_dp,now%smb_ice - snow_to_ice/par%tsticsub)
              end where
          end if
  
  
negative.code: 
          where (now%mask >= 1)
              ! melt energy where temperature exceeds freezing (difference to heat at freezing)
              qmelt = dmax1(0.0_dp,(above)*par%ceff/par%tstic+now%qmr)
              ! watch the sign
              qcold = dmax1(0.0_dp,abs(below)*par%ceff/par%tstic)
          else where
              qmelt = 0.0_dp
              qcold = 0.0_dp
          end where
  
          ! 1) ablation: melt (m/s); potential melt resulting from available melt energy
          if (.not. bnd%melt) then
              ! potential melt
              now%melt = qmelt/(rhow*clm)
              ! separate potential melt into actual melt of snow and ice
              now%melted_snow = dmin1(now%melt,now%hsnow/par%tstic)
              now%melted_ice  = now%melt-now%melted_snow
  
              ! actual melt is sum of melted snow and ice (melted snow over land)
  