anchor.diff:              now%refr = refrozen_rain + refrozen_snow
           end if 
   
  -        ! 3) potential runoff
  +        !> 6. Runoff
           now%runoff = 0.0_dp
           now%runoff = now%melt + now%rf - refrozen_rain
   
  -        ! 4) accumulation: sum of all incoming solid water (just diagnostic, here)
  +        !> 7. Accumulation: sum of all incoming solid water (just diagnostic, here)
           if (.not. bnd%acc) then
               now%acc = now%sf - now%subl/rhow + now%refr
           end if
           
  -        ! 5) surface mass balance
  +        !> 8. Surface mass balance of snow
           now%smb_snow = now%sf - now%subl/rhow - now%melted_snow
   
           where (now%mask == 0)
               now%hsnow = 0.0_dp
           else where
  -            ! update snow height
  +            !> 9. 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)
  +        !> 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 + 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
   
  +        !> 11. Total surface mass balance
           if (.not. bnd%smb) then
               where (now%mask == 2)
                   now%smb = now%smb_snow + now%smb_ice - snow_to_ice/par%tstic
  
positive.code:         end if
      end subroutine
  
      elemental subroutine longwave_upward(ts,sigma,lwout)
          double precision, intent(in) :: ts, sigma
          double precision, intent(out) :: lwout
          lwout = sigma*ts*ts*ts*ts
      end subroutine
  
      elemental subroutine diurnal_cycle(amp,tmean,above,below)
          ! Calculate analytical expression for above-/below-freezing
          ! temperatures for a given mean temperature.
          ! Diurnal cycle amplitude can be either fixed (recommended)
          ! or decrease with surface temperature/pressure (name list)
  
          double precision, intent(in) :: tmean
          double precision, intent(in) :: amp
          double precision, intent(out) :: above, below
          double precision :: tmp1, tmp2
  
          tmp1 = 0.0_dp
  
negative.code:  
  
  
      !> Allocation of all fields of the surface_physics::state_class
      subroutine surface_alloc(now,npts)
  
          implicit none 
  