anchor.diff:          character(len=*)  :: step
           double precision, optional :: nt 
           
  -        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%massbal, now%massbal,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%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)
  +        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
   
  
positive.code:         tmin          = par%tmin      ! minimum albedo-affecting temperature
          tmax          = par%tmax      ! maximum albedo-affecting temperature (originally 273.15K)
          tstic         = par%tstic
          Ts_wt         = par%Ts_wt     ! Weight between surface and free-atmos temp.
          method        = par%method 
          alb_scheme    = par%alb_scheme 
          ! PDD
          pdd_sigma     = par%pdd_sigma
          pdd_b         = par%pdd_b
          ! ITM
          itm_cc        = par%itm_cc
          itm_t         = par%itm_t
          tau_a         = par%tau_a  
          tau_f         = par%tau_f  
          w_crit        = par%w_crit
  
negative.code:     end function ei_sat
  
  
      !! Data management subroutines 
  
      subroutine surface_physics_average(ave,now,step,nt)
  