anchor.diff:          double precision, allocatable, dimension(:) :: runoff      !< potential surface runoff [m/s]
           double precision, allocatable, dimension(:) :: qmr         !< heat flux from melting/refreezing [W/m2]
           double precision, allocatable, dimension(:) :: qmr_res     !< residual heat flux from melting/refreezing(at end of time step) [W/m2]
  +        double precision, allocatable, dimension(:) :: amp         !< diurnal cycle amplitude [K]
           ! Forcing variables
           double precision, allocatable, dimension(:) :: sf          !< snow fall [m/s]
           double precision, allocatable, dimension(:) :: rf          !< rain fall [m/s]
  
positive.code:         character(len=*)  :: step                        !< current step:
                                                           !! "init", "step", or "end"
          double precision, optional, intent(in) :: nt     !< number of total steps
          
          call field_average(ave%t2m,     now%t2m,  step,nt)
          call field_average(ave%tsurf,   now%tsurf,step,nt)
          call field_average(ave%hsnow,   now%hsnow,step,nt)
          call field_average(ave%alb,     now%alb,  step,nt)
          call field_average(ave%melt,    now%melt, step,nt)
          call field_average(ave%refr,    now%refr, step,nt)
          call field_average(ave%smb,     now%smb,  step,nt)
          call field_average(ave%acc,     now%acc,  step,nt)
          call field_average(ave%lhf,     now%lhf,  step,nt)
          call field_average(ave%shf,     now%shf,  step,nt)
          call field_average(ave%lwu,     now%lwu,  step,nt)
  
          call field_average(ave%sf,      now%sf,   step,nt)
          call field_average(ave%rf,      now%rf,   step,nt)
          call field_average(ave%sp,      now%sp,   step,nt)
          call field_average(ave%lwd,     now%lwd,  step,nt)
          call field_average(ave%swd,     now%swd,  step,nt)
          call field_average(ave%wind,    now%wind, step,nt)
          call field_average(ave%rhoa,    now%rhoa, step,nt)
          call field_average(ave%qq,      now%qq,   step,nt)
  
          return
  
  
negative.code:             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)
              if (trim(par%alb_scheme) .eq. "none") then
                  now%alb_snow = par%alb_smax
              end if
              end if
              where (now%mask == 2)
                  now%alb = par%albi + f_alb*(now%alb_snow - par%albi)
  