anchor.diff:              qcold = 0.0_dp
           end where
   
  -        ! 1) ablation: melt (m/s); potential melt resulting from available melt energy
  +        !> 4. Ablation: melt (m/s); potential melt resulting from available melt energy
           if (.not. bnd%melt) then
               ! potential melt
               now%melt = qmelt/(rhow*clm)
  
positive.code:         if (tsurf > t0) then
              tm = 1.0_dp
          end if
          ! In contrast to the formulation in their paper, I summed up alpha_nir
          ! and alpha_nir immediately (fewer parameters: alb_smax and alb_smin).
          alb = alb_smax - (alb_smax - alb_smin)*tm**3.0_dp
          ! snow cover fraction and gridpoint-averaged albedo
          !sfrac = tanh(hsnow/hcrit)
          !alb = albr*(1.-sfrac) + alb*sfrac
          !alb = albr + hsnow/(hsnow + hcrit)*(alb - albr)
          ! write albedo to domain object 'now'
      end subroutine albedo_slater
  
      elemental subroutine albedo_denby(alb,melt,alb_smax, alb_smin, mcrit)
          double precision, intent(out) :: alb
          double precision, intent(in) :: melt
          double precision, intent(in) :: alb_smax, alb_smin, mcrit
          alb = alb_smin + (alb_smax - alb_smin)*dexp(-melt/mcrit)
      end subroutine albedo_denby
  
      ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      ! Subroutine : a l b e d o _ r e m b o
      ! Author     : Alex Robinson, modified by M. Krapp
      ! Purpose    : Determine the current surface
      ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      elemental subroutine albedo_rembo(as, hsnow, albr,alb_smin,alb_smax, melt_in)
        
        implicit none
        
        double precision, intent(in out)       :: as
        double precision, intent(in)           :: hsnow
        double precision, intent(in)           :: albr, alb_smax, alb_smin
        double precision, intent(in), optional :: melt_in
        double precision :: as_snow, as_ground, depth, as_snow_planet
        double precision :: hsnow_critical, melt
        
        ! Determine the critical snow height based on number of PDDs,
        ! which correspond to vegetation growth
        ! 0-100 PDDs    = desert, 10mm
        ! 100-1000 PDDs = tundra (grass), linearly increasing between 10mm => 100mm
        hsnow_critical = 0.1_dp ! 10 mm
        
        ! Determine the scaled snowdepth
        depth = hsnow / hsnow_critical
        
        ! Determine the amount of melt to affect abledo calculation.
        ! Default is high melt everywhere, so that jump in albedo is avoided
        ! at onset of melt season.
        ! After calculating actual melt, this is included as an argument and
        ! albedo is recalculated.
        melt = 1.e-3_dp + 1.e-3_dp
        if ( present(melt_in) ) melt = melt_in
        
        ! Figure out what the surface albedo would be with no snow, 
        ! based on the type of ground underneath ( ice or land )
        as_ground = albr
        
        ! Determine current snow albedo: if melt gt eg, 1mm/day, then change albedo to melting snow!
        as_snow = alb_smax
        if (melt .gt. 1.e-3_dp) as_snow = alb_smin
        
        ! Determine snow albedo for use with REMBO (to get planetary albedo)
        ! Where PDDs > 1d3, forest snow
        as_snow_planet = as_snow
  
        ! First, calculate surface albedo to use for REMBO, 
        ! which is not necessarily identical to snow albedo
        ! minimum albedo is that of bare ground
        as = min(as_ground + depth*(as_snow_planet-as_ground), as_snow_planet)
        
        ! Get current surface albedo to be used for ITM melt scheme
        ! It will either be the maximum albedo (that of dry snow: as_snow)
        ! or the wet snow albedo plus a fraction depending on height of snow
        ! minimum albedo now should be that of wet snow minimum
        as = min(albr + depth*(as_snow-albr), as_snow)
        
        return
        
      end subroutine albedo_rembo
  
      elemental function ew_sat(t) result(fsat)
  
          double precision, intent(in) :: t
          double precision :: fsat
          fsat = 611.2_dp*dexp(17.62_dp*(t-t0)/(243.12_dp+t-t0))
      end function ew_sat
  
      elemental function ei_sat(t) result(fsat)
  
          double precision, intent(in) :: t
          double precision :: fsat
          fsat = 611.2_dp*dexp(22.46_dp*(t-t0)/(272.62_dp+t-t0))
      end function ei_sat
  
  
      !! Data management subroutines 
  
      subroutine surface_physics_average(ave,now,step,nt)
          implicit none 
  
          type(surface_state_class), intent(INOUT) :: ave
          type(surface_state_class), intent(IN)    :: now 
          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)
  
negative.code:                 now%qmr = 0.0_dp
              end where
          end if
          !> 7. Update 2m air temperature
          now%t2m   = now%t2m + (now%shf+now%lhf)*par%tsticsub/par%ceff  
  
      end subroutine
  
      !>  Main routine to calculate
      !!  1. surface mass balance (m/s)
      !!  2. surface (snow/ice) albedo
      !!  3. snow height (m)
      !! based on surface energy and mass balance.
      !! Forced by atmospheric fields of
      !!  * air temperature
      !!  * surface wind speed
      !!  * air humidity
      !!  * snow fall
      !!  * rain fall
      !!  * surface pressure
      !!  * air density
      !!  * longwave radiation
      !!  * shortwave radiation
      subroutine mass_balance(now,par,bnd,day,year)
  
          type(surface_state_class), intent(in out) :: now !< object with daily variables
          type(boundary_opt_class),  intent(in)     :: bnd !< object with boundary switches
          type(surface_param_class), intent(in)     :: par !< object with parameter values
  
          integer, intent(in) :: day !< current day
          integer, intent(in) :: year !< current year
          
          ! auxillary variables
          double precision, dimension(par%nx) :: qmelt, qcold, &
              below, above, f_rz, f_alb, refrozen_rain, refrozen_snow, snow_to_ice 
  
          !> 1. Calculate above-/below-freezing temperatures for a given mean temperature
          call diurnal_cycle(par%amp,now%tsurf-t0,above,below)
  
          where (now%mask >= 1)
              !> 2. Calculate Melt energy where temperature exceeds freezing (difference to heat at freezing)
              qmelt = dmax1(0.0_dp,(above)*par%ceff/par%tstic+now%qmr)
              !> 3. Calcukate "cold" content
              ! watch the sign
              qcold = dmax1(0.0_dp,abs(below)*par%ceff/par%tstic)
          else where
  