anchor.diff:    ! using n-layer physics it is constrained to balance the layers
     ! above it.
     if (RedGrav) then
  -    do j=1, ny
  -      do i=1, nx
  +    do j = 1, ny
  +      do i = 1, nx
           dhdt_GM(i,j,layers) = &
               ah(layers)*(h(i+1,j,layers)*wetmask(i+1,j)   &
                 + (1d0 - wetmask(i+1,j))*h(i,j,layers)     & ! boundary
  
positive.code:   dhdt = 0d0
  
    do k = 1, layers
      do j=1, ny
        do i=1, nx
          dhdt(i,j,k) = &
              dhdt_GM(i,j,k) & ! horizontal thickness diffusion
              - ((h(i,j,k)+h(i+1,j,k))*u(i+1,j,k) &
  
negative.code:     ! 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
  
      unew = u+dt*(23d0*dudt - 16d0*dudtold + 5d0*dudtveryold)/12d0
      vnew = v+dt*(23d0*dvdt - 16d0*dvdtold + 5d0*dvdtveryold)/12d0
      hnew = h+dt*(23d0*dhdt - 16d0*dhdtold + 5d0*dhdtveryold)/12d0
  
      ! Apply the boundary conditions
      ! - Enforce no normal flow boundary condition
  