anchor.diff:              !> 9. Update snow height
               now%hsnow = dmax1(0.0_dp, now%hsnow + now%smb_snow*par%tstic)
           end where
  -        
  +
           !> 10. Relax snow height to maximum (eg, 5 m)
           snow_to_ice   = dmax1(0.d0,now%hsnow-hsmax)
           now%hsnow   = now%hsnow - snow_to_ice
  
positive.code:                 now%melt = now%melted_snow + now%melted_ice
              elsewhere
                  now%melt = now%melted_snow
                  !now%qmr = now%qmr + now%melted_ice*rhow*clm
                  now%melted_ice = 0.0_dp
              end where
          end if
          ! energy needed for melting is added to residual energy
          !now%qmr = now%qmr + now%melt*rhow*clm
  
          !> 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
              !else where
                  !f_rz = now%hsnow/(now%hsnow+par%rcrit)
                  !f_rz = (1.0_dp+tanh(now%hsnow-par%rcrit))/2.0_dp
              !f_rz = dmin1(1.0_dp,now%hsnow/par%rcrit)
              !end where
              f_rz = 1.0_dp! - exp(-now%hsnow/par%rcrit)
              ! potential refreezing
              now%refr = qcold/(rhow*clm)
              refrozen_rain = dmin1(now%refr,now%rf)
  
negative.code:         write(*,'(a,l1)') 'shf     ', bnd%shf
          write(*,'(a,l1)') 'subl    ', bnd%subl
          write(*,'(a,l1)') 'smb     ', bnd%smb
          write(*,'(a,l1)') 'amp     ', bnd%amp
      end subroutine print_boundary_opt
  
  
  
  end module surface_physics
  