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:         ! 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)*exp(-tau_f*tstic/tau) + alb_smin
          alb_new = sf*tstic/(w_crn/rhow)*(alb_smax-alb_smin)
  
          ! dry/wet-averaged albedo
  
negative.code: 
          !> 5. Refreezing (m/s) as fraction of melt (increases with snow height)
          if (.not. bnd%refr) then
              f_rz = par%rcrit!1.0_dp - dexp(-now%hsnow/(par%rcrit+epsil))
              ! potential refreezing
              now%refr = qcold/(rhow*clm)
              refrozen_rain = dmin1(now%refr,now%rf)
              ! potential refreezing 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*dmin1(now%hsnow/par%tstic,refrozen_rain)
              refrozen_snow = f_rz*dmin1(now%hsnow/par%tstic,refrozen_snow)
              now%refr = refrozen_rain + refrozen_snow
              ! energy released during refreezing that has not been used
              ! is subtracted from residual energy
  