anchor.diff:          !> 2. surface_physics::latent_heat_flux \n
           !! bulk formulation of latent heat flux (W/m^2), only accounts for sublimation/deposition,
           !! not for evaporation/condensation (would require estimate of liquid water content)
  -        if (.not. bnd%lhf) then  
  +        if (.not. bnd%lhf) then
               now%subl = 0.0_dp
               now%evap = 0.0_dp
               now%lhf  = 0.0_dp
               call latent_heat_flux(now%tsurf, now%wind, now%qq, now%sp, now%rhoa, now%mask, &
                                     par%clh, eps, cls, clv, now%lhf, now%subl, now%evap)
           end if
  +        now%subl = now%subl/rhow ! sublimation was calculated as kg/(s m2) -> m/s
   
           !> 3. surface_physics::longwave_upward \n
           !! outgoing long-wave flux following Stefan-Boltzmann's law (W/m^2)
  
positive.code:         character(len=*)  :: step                        !< current step:
                                                           !! "init", "step", or "end"
          double precision, optional, intent(in) :: nt     !< number of total steps
          
          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
  
  
negative.code:             where (now%mask == 0)
                  now%alb = 0.06_dp
              end where
  
              if (trim(par%alb_scheme) .eq. "alex") then
                  now%alb_snow = par%alb_smin + (par%alb_smax - par%alb_smin)*(0.5_dp*tanh(par%afac*(now%t2m-par%tmid))+0.5_dp) 
                  now%alb      = now%alb_snow
              end if
  
          end if
  
          ! store residual energy
          where (now%mask == 0)
              now%qmr_res = 0.0_dp
          elsewhere
              now%qmr_res = now%qmr
          end where
  
      end subroutine mass_balance
  
  