anchor.diff:      double precision, parameter :: cap  = 1000.0_dp  !< specific heat capacitiy of air [J/(kg K)]
       double precision, parameter :: rhow = 1000.0_dp  !< density of water [kg/m3]
       double precision, parameter :: hsmax= 5.0_dp     !< maximum snow height [m]
  +    double precision, parameter :: epsil= epsilon(1.0_dp) !< smallest double precision number
   
   
       type surface_param_class !< Define all parameters needed for the surface module
  
positive.code: 
          !> 5. Refreezing (m/s) as fraction of melt (increases with snow height)
          if (.not. bnd%refr) then
              f_rz = 0.0_dp
              where (now%hsnow>0.0_dp) f_rz = par%rcrit!1.0_dp - exp(-now%hsnow/par%rcrit)
              ! potential refreezing
              now%refr = qcold/(rhow*clm)
              refrozen_rain = dmin1(now%refr,now%rf)
              ! potential refeezing snow
              refrozen_snow = dmax1(now%refr-refrozen_rain,0.0_dp)
              ! actual refreezing snow
              refrozen_snow = dmin1(refrozen_snow,now%melted_snow)
              ! actual refreezing
              refrozen_rain =  f_rz*refrozen_rain
              refrozen_snow =  f_rz*refrozen_snow
              now%refr = refrozen_rain + refrozen_snow
              ! energy released during refreezing that has not been used
              ! is subtracted from residual energy
  
negative.code:         ! where no melting occurs, albedo decreases linearly
          alb_dry = alb - tau_a*tstic/tau
          !where melting occurs, albedo decreases exponentially
          alb_wet = (alb - alb_smin)*dexp(-tau_f*tstic/tau) + alb_smin
          alb_new = sf*tstic/(w_crn/rhow)*(alb_smax-alb_smin)
  
          ! dry/wet-averaged albedo
  