anchor.diff:      public :: surface_physics_par_load, surface_alloc, surface_dealloc
       public :: surface_boundary_define, surface_physics_average
       public :: print_param, print_boundary_opt
  +
  +    public :: sigm, cap, eps, cls, clv  
  +    public :: longwave_upward, latent_heat_flux, sensible_heat_flux 
   
   contains
   
  
positive.code:                     bnd%melt    = .TRUE. 
                  case("refr")
                      bnd%refr    = .TRUE. 
                  case("massbal")
                      bnd%massbal = .TRUE. 
                  case("acc")
                      bnd%acc     = .TRUE. 
                  case("lhf")
  
negative.code:     ! Define all parameters needed for the surface module
      type surface_param_class
          character (len=256) :: name, boundary(30), alb_scheme
          integer             :: nx        !< number of grid points
          integer             :: n_ksub    !< number of sub-daily time steps
          double precision    :: ceff      !< effective volumetric heat capacity of snow/ice [J/m3]
          double precision    :: albr      !< background albedo (bare ice) [no unit]
          double precision    :: albl      !< background albedo (bare land) [no unit]
          double precision    :: alb_smax  !< maximum snow albedo (fresh snow) [no unit]
          double precision    :: alb_smin  !< minimum snow albedo (old, wet snow) [no unit]
          double precision    :: hcrit     !< critical snow height for which grid cell is 50%snow covered [m]
          double precision    :: rcrit     !< critical snow height for which refreezing fraction is 50% [m]
          double precision    :: amp       !< Amplitude of diurnal cycle [K]
          double precision    :: csh       !< sensible heat exchange coefficient [no unit]
          double precision    :: clh       !< latent heat exchange coefficient [no unit]
          double precision    :: tmin      !< minimum temperature for which albedo decline becomes effective ("slater") [K]
          double precision    :: tmax      !< maximum temperature for which albedo decline becomes effective ("slater") [K]
          double precision    :: tstic     !< time step [s]
          double precision    :: shf_enh   !< sensible heat flux enhancement factor for positive u*(tsurf-t2m) [no unit]
          double precision    :: lhf_enh   !< latent heat flux enhancement factor over land (mask==1) [no unit]
          double precision    :: tau_a     !< dry albedo decline for "isba" albedo scheme [1/day]
          double precision    :: tau_f     !< wet albedo decline for "isba" albedo scheme [1/day]
          double precision    :: w_crit    !< critical liquid water content for "isba" albedo scheme [kg/m2]
          double precision    :: mcrit     !< critical melt rate for "isba" and "denby" albedo scheme [m/s]
          double precision    :: afac      !< param [no unit]
          double precision    :: tmid      !< param for "alex" albedo parametrization [K]
      end type
  
      type surface_state_class
          ! Model variables
          double precision, allocatable, dimension(:) :: t2m      !< 2m air temperature [K]
          double precision, allocatable, dimension(:) :: tsurf    !< surface temperature [K]
          double precision, allocatable, dimension(:) :: hsnow    !< snow pack height (water equivalent) [m]
          double precision, allocatable, dimension(:) :: hice     !< ice thickness (water equivalent) [m]
          double precision, allocatable, dimension(:) :: alb      !< grid-averaged albedo [no unit]
          double precision, allocatable, dimension(:) :: alb_snow !< snow albedo [no unit]
          double precision, allocatable, dimension(:) :: melt     !< surface melt [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 [??]
          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   !< surface runoff [m/s]
          ! 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 
          logical :: t2m, tsurf, hsnow, alb, melt, refr, smb, acc, lhf, shf, subl
      end type
  
      type surface_physics_class
  