anchor.diff:          
           ! physical parameters used in the different parameterizations
           double precision :: tmin, tmax, albr, albl, alb_smax, alb_smin, hcrit, rcrit, &
  -                            ceff, tstic, amp, clh, csh, &
  +                            ceff, tstic, amp, clh, csh, shf_enh, lhf_enh, &
                               pdd_sigma, pdd_b, &
                               itm_cc, itm_t, &
                               tau_a, tau_f, w_crit, mcrit, &
                               mmfac, Pmaxfrac, &
  -                            Ts_wt(2), afac,tmid
  +                            Ts_wt(2), afac,tmid 
           double precision :: tsticsub
           integer :: ksub, n_ksub
   
  
positive.code:         sf     = dom%now%sf
          rf     = dom%now%rf
          
          ! assign values from previous timestep 
          melt    = dom%now%melt 
          massbal = dom%now%massbal 
  
          hsnow_max    = 10.d0 
          melt_snow    = 0.d0 
          melt_ice     = 0.d0 
          massbal_snow = 0.d0 
  
negative.code:         if (.not. bnd%hsnow)   dom%now%hsnow   = hsnow
          dom%now%hice    = new_ice + dom%now%hice
          dom%now%alb_snow = alb_snow
          dom%now%massbal_snow = massbal_snow
          dom%now%massbal_ice  = massbal_ice
  
      end subroutine surface_mass_balance
  
      elemental subroutine sensible_heat_flux(ts,ta,wind,rhoa,csh,enh,cap,shf)
          double precision, intent(in) :: ts, ta, wind, rhoa
          double precision, intent(in) :: csh, cap, enh
          double precision, intent(out) :: shf
          double precision :: coeff
  
          ! enhancement over land
          if ( wind*(ts-ta) > 0.0_dp ) then
              coeff = enh*csh
          else
              coeff = csh
          end if
          
          shf = coeff*cap*rhoa*wind*(ts-ta)
      end subroutine
  
      elemental subroutine latent_heat_flux(ts, wind, shum, sp, rhoatm, lmask, clh, enh, eps, cls, clv, lhf, subl, evap)
          double precision, intent(in)  :: ts, shum, sp, rhoatm, wind
          double precision, intent(in)  :: clh, eps, cls, clv, enh
          integer, intent(in)  :: lmask
          double precision, intent(out) :: lhf, subl, evap
          double precision :: esat_sur, shum_sat, coeff
  
          subl = 0.0_dp
          evap = 0.0_dp
          lhf  = 0.0_dp
          if (lmask == 1) then
              coeff = enh*clh
          else
              coeff = clh
          end if
          if (ts < t0) then
              esat_sur = ei_sat(ts)
              ! specific humidity at surface (assumed to be saturated) is
              shum_sat = esat_sur*eps/(esat_sur*(eps-1.0_dp)+sp)
              ! sublimation/deposition depends on air specific humidity
              subl = coeff*wind*rhoatm*(shum_sat - shum)
              lhf = subl*cls
          else
              esat_sur = ew_sat(ts)
              ! evaporation/condensation
              ! specific humidity at surface (assumed to be saturated) is
              shum_sat = esat_sur*eps/(esat_sur*(eps-1.0_dp)+sp)
              evap = coeff*wind*rhoatm*(shum_sat - shum)
              lhf = evap*clv
          end if
      end subroutine
  