anchor.diff:          ! energy needed for melting is added to residual energy
           !now%qmr = now%qmr + now%melt*rhow*clm
   
  -        !> 5. Refreezing as fraction of melt (increases with snow height)
  +        !> 5. Refreezing (m/s) as fraction of melt (increases with snow height)
           if (.not. bnd%refr) then
               !where (now%mask == 2)
               !    f_rz = 1.0_dp
  
positive.code:         
          !> 8. Surface mass balance of snow
          if (.not. bnd%smb) then
              now%smb_snow = now%sf - now%subl/rhow - now%melted_snow + refrozen_snow
          end if
  
          where (now%mask == 0)
  
negative.code:         double precision, allocatable, dimension(:) :: melt        !< potential surface melt [m/s]
          double precision, allocatable, dimension(:) :: melted_snow !< actual melted snow [m/s]
          double precision, allocatable, dimension(:) :: melted_ice  !< actual melted ice [m/s]
          double precision, allocatable, dimension(:) :: refr        !< refreezing [m/s]
          double precision, allocatable, dimension(:) :: smb         !< surface mass balance [m/s]
          double precision, allocatable, dimension(:) :: acc         !< surface accumulation [m/s]
          double precision, allocatable, dimension(:) :: lhf         !< latent heat flux [W/m2]
          double precision, allocatable, dimension(:) :: shf         !< sensible heat flux [W/m2]
          double precision, allocatable, dimension(:) :: lwu         !< upwelling longwave radiation [W/m2]
          double precision, allocatable, dimension(:) :: subl        !< sublimation [m/s]
          double precision, allocatable, dimension(:) :: evap        !< evaporation [??]
          double precision, allocatable, dimension(:) :: smb_snow    !< surface mass balance of snow [m/s]
          double precision, allocatable, dimension(:) :: smb_ice     !< surface mass balance of ice [m/s]
          double precision, allocatable, dimension(:) :: runoff      !< potential surface runoff [m/s]
          double precision, allocatable, dimension(:) :: qmr         !< heat flux from melting/refreezing [W/m2]
          double precision, allocatable, dimension(:) :: qmr_res     !< residual heat flux from melting/refreezing(at end of time step) [W/m2]
          ! Forcing variables
          double precision, allocatable, dimension(:) :: sf          !< snow fall [m/s]
          double precision, allocatable, dimension(:) :: rf          !< rain fall [m/s]
          double precision, allocatable, dimension(:) :: sp          !< surface pressure [Pa]
          double precision, allocatable, dimension(:) :: lwd         !< downwelling longwave radiation [W/m2]
          double precision, allocatable, dimension(:) :: swd         !< downwelling shortwave radiation [W/m2]
          double precision, allocatable, dimension(:) :: wind        !< surface wind speed [m/s]
          double precision, allocatable, dimension(:) :: rhoa        !< air density [kg/m3]
          double precision, allocatable, dimension(:) :: qq          !< air specific humidity [kg/kg]
          integer,          allocatable, dimension(:) :: mask        !< ocean/land/ice mask [0/1/2]
      end type
  
      type boundary_opt_class !< object to hold flags for overriding with external fields
          logical :: t2m   !< flag for 2m air temperature
          logical :: tsurf !< flag for surface temperature
          logical :: hsnow !< flag for snow height
          logical :: alb   !< flag for albedo
          logical :: melt  !< flag for melting
          logical :: refr  !< flag for refreezing
          logical :: smb   !< flag for surface mass balance
          logical :: acc   !< flag for accumulation
          logical :: lhf   !< flag for latent heat flux
          logical :: shf   !< flag for senible heat flux
          logical :: subl  !< flag for sublimation/refreezing
      end type
  
      type surface_physics_class !< container which holds all used objects
  