anchor.diff:  ! hard-coded numbers
         surface_physics_nml = 'semic.nml'
         forc_file = 'forcing.nc'
  +      mask_file = "glac.nc"
         nx = 192
         ny = 96
         nt = 365
  
positive.code:       dom%surface%now%tsurf = 260.0
        dom%surface%now%alb_snow = 0.8
  
        dom%surface%now%mask = 2.0 ! all points treated as ice point (0/1/2 = ocean/land/ice)
  
        ! read model parameters from namelist file
        dom%surface_physics_nml = surface_physics_nml
  
negative.code: 
        end subroutine update
  
        subroutine read_mask(mnm_in, nx, ny, mask)
  
          integer, intent(in) :: nx
          integer, intent(in) :: ny
          character(len=7), intent(in) :: mnm_in !< GLAC file
  
          integer, allocatable, dimension(:,:) :: tmp
          integer, dimension(nx*ny), intent(out) :: mask
  
          allocate(tmp(ny,ny))
          call nc_read(mnm_in,"glac", tmp)
          mask = reshape(tmp, shape(mask))
          deallocate(tmp)
  
        end subroutine read_mask
  
        subroutine read_forcing(fnm_in,nx,ny,nt,wind,tt,sf,rf,qq,rhoa,sp,lwd,swd)
            
            integer, intent(in) :: nx !< x dimension
  