anchor.diff:  subroutine calc_boundary_masks(wetmask, hfacW, hfacE, hfacS, hfacN, nx, ny)
     implicit none
   
  -  integer nx !< number of grid points in x direction
  -  integer ny !< number of grid points in y direction
  -  double precision wetmask(0:nx+1, 0:ny+1)
  +  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)
  -  double precision hfacW(0:nx+1, 0:ny+1)
  -  double precision hfacE(0:nx+1, 0:ny+1)
  -  double precision hfacN(0:nx+1, 0:ny+1)
  -  double precision hfacS(0:nx+1, 0:ny+1)
     integer i, j
   
     hfacW = 1d0
  
positive.code: subroutine break_if_NaN(data, nx, ny, layers, n)
    implicit none
  
    ! To stop the program if it detects at NaN in the variable being checked
  
    integer nx, ny, layers, n, i, j, k
    double precision data(0:nx+1, 0:ny+1, layers)
  
    do k = 1, layers
      do j = 1, ny
  
negative.code: subroutine read_input_fileH(name, array, default, nx, ny, layers)
    implicit none
  
    character(30), intent(in) :: name
    double precision, intent(out) :: array(0:nx+1, 0:ny+1, layers)
    double precision, intent(in) :: default(layers)
    integer, intent(in) :: nx, ny, layers
  
    double precision array_small(nx, ny, layers)
    integer k
  
    if (name.ne.'') then
      open(unit=10, form='unformatted', file=name)
  