anchor.diff:  ! -----------------------------------------------------------------------------
   !> Check to see if there are any NaNs in the data field and stop the
   !! calculation if any are found.
  -subroutine break_if_NaN(data,nx,ny,layers,n)
  +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
  -      do i=1,nx
  -        if (data(i,j,k).ne.data(i,j,k))  then
  +  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
  +      do i=1, nx
  +        if (data(i,j,k) .ne. data(i,j,k)) then
             ! write a file saying so
             open(unit=10, file='NaN detected.txt', action="write", &
                 status="replace", form="formatted")
  -          write(10,1000) n
  +          write(10, 1000) n
   1000      format( "NaN detected at time step ", 1i10.10)
             close(unit=10)
             ! print it on the screen
  
positive.code:   if (.not. RedGrav) then
      ! Initialise arrays for pressure solver
      ! a = derivatives of the depth field
      do j=1,ny
        do i=1,nx
          a(1,i,j)=g_vec(1)*0.5*(depth(i+1,j)+depth(i,j))/dx**2
          a(2,i,j)=g_vec(1)*0.5*(depth(i,j+1)+depth(i,j))/dy**2
          a(3,i,j)=g_vec(1)*0.5*(depth(i,j)+depth(i-1,j))/dx**2
          a(4,i,j)=g_vec(1)*0.5*(depth(i,j)+depth(i,j-1))/dy**2
        end do
      end do
      do j=1,ny
        a(1,nx,j)=0.0
        a(3,1,j)=0.0
      end do
      do i=1,nx
        a(2,i,ny)=0.0
        a(4,i,1)=0.0
      end do
      do j=1,ny
        do i=1,nx
          a(5,i,j)=-a(1,i,j)-a(2,i,j)-a(3,i,j)-a(4,i,j)
        end do
      end do
  
negative.code: 
      ! Calculate Bernoulli potential
      if (RedGrav) then
        call evaluate_b_RedGrav(b, h, u, v, nx, ny, layers, g_vec)
      else
        call evaluate_b_iso(b, h, u, v, nx, ny, layers, g_vec, depth)
      end if
  
      ! Calculate relative vorticity
      call evaluate_zeta(zeta, u, v, nx, ny, layers, dx, dy)
  
      ! Calculate dhdt, dudt, dvdt at current time step
      call evaluate_dhdt(dhdt, h, u, v, ah, dx, dy, nx, ny, layers, &
          spongeHTimeScale, spongeH, wetmask, RedGrav)
  
      call evaluate_dudt(dudt, h, u, v, b, zeta, wind_x, fu, au, ar, slip, &
          dx, dy, hfacN, hfacS, nx, ny, layers, rho0, &
          spongeUTimeScale, spongeU, RedGrav, botDrag)
  
      call evaluate_dvdt(dvdt, h, u, v, b, zeta, wind_y, fv, au, ar, slip, &
          dx, dy, hfacW, hfacE, nx, ny, layers, rho0, &
          spongeVTimeScale, spongeV, RedGrav, botDrag)
  
      ! Use dh/dt, du/dt and dv/dt to step h, u and v forward in time with
      ! the Adams-Bashforth third order linear multistep method
  