anchor.diff:  subroutine read_input_fileH_2D(name, array, default, nx, ny)
     implicit none
   
  -  character(30) name
  -  integer nx, ny
  -  double precision array(0:nx+1, 0:ny+1), default
  +  character(30), intent(in) :: name
  +  double precision, intent(out) :: array(0:nx+1, 0:ny+1)
  +  double precision, intent(in) :: default
  +  integer, intent(in) :: nx, ny
  +
     double precision array_small(nx, ny)
   
     if (name.ne.'') then
  
positive.code: subroutine read_input_fileH(name, array, default, nx, ny, layers)
    implicit none
  
    character(30) name
    integer nx, ny, layers, k
    double precision array(0:nx+1, 0:ny+1, layers), default(layers)
    double precision array_small(nx, ny, layers)
  
    if (name.ne.'') then
      open(unit=10, form='unformatted', file=name)
  
negative.code: subroutine calc_boundary_masks(wetmask, hfacW, hfacE, hfacS, hfacN, nx, ny)
    implicit none
  
    double precision, intent(in)  :: wetmask(0:nx+1, 0:ny+1)
    double precision, intent(out) :: hfacW(0:nx+1, 0:ny+1)
    double precision, intent(out) :: hfacE(0:nx+1, 0:ny+1)
    double precision, intent(out) :: hfacN(0:nx+1, 0:ny+1)
    double precision, intent(out) :: hfacS(0:nx+1, 0:ny+1)
    integer, intent(in) :: nx !< number of grid points in x direction
    integer, intent(in) :: ny !< number of grid points in y direction
  
    double precision temp(0:nx+1, 0:ny+1)
    integer i, j
  
    hfacW = 1d0
  