anchor.diff:      end do lz
     end subroutine BuildRhs3D
   
  +  subroutine GSPSlv2D(FldSlv, h, MaxIt, DErr)
  +    ! Gauss-Seidel pressure solve.
  +    class(FluidSlv), intent(inout) :: FldSlv
  +    real(dp)     , intent(in)    :: h      ! Timestep.
  +    integer      , intent(in)    :: MaxIt  ! Maximum number of iterations.
  +    real(dp)     , intent(in)    :: DErr   ! Desired error threshold.
  +    real(dp)                     :: scale, Err ! Scaling factor, iteration error.
  +    integer                      :: iter, i, j, idx ! Iteration, counters, index
  +    real(dp)                     :: Diag, OffDiag ! Diagonal and offdiagonal elements of the matrix.
  +    real(dp)                     :: NewP ! New pressure value.
  +
  +    ! Set scale.
  +    scale = h / (FldSlv % srho * FldSlv % celsiz * FldSlv % celsiz)
  +
  +    ! ITERation Loop.
  +    iterl: do iter = 1, MaxIt
  +      Err = 0.
  +      ! Loop over y.
  +      ly: do j = 1, FldSlv % whd(2)
  +        ! Loop over x.
  +        lx: do i = 1, FldSlv % whd(1)
  +          idx     = i + FldSlv % whd(1)*(j-1)
  +          Diag    = 0.
  +          OffDiag = 0.
  +
  +          ! Building the matrix as a five-point stencil. Solid borders, i.e. no fluid is outside of the simulation bounds.
  +
  +          ! The pressure contribution for each cell depends on its 4 nearest neighbours.
  +          ! For example, cell a22 will receive contributions from cells a12, a23, a33 and a21.
  +          ! Boundary cells will only receive contributions from the orthogonally placed cells that are found within the simulation bounds. For example, cell a13 will only receive contributions from cells a12 and a23, while cell a21 from cells a11, a22 and a31. This leads to the series of if statements found herein.
  +
  +          ! y-axis
  +          ! ^
  +          ! |__> x-axis
  +          !
  +          ! -------------------
  +          ! | a13 | a23 | a33 |
  +          ! -------------------
  +          ! | a12 | a22 | a32 |     --------->   [a11,a21,a31,a12,a22,a32,a13,a23,a33]
  +          ! -------------------
  +          ! | a11 | a21 | a31 |
  +          ! -------------------
  +
  +          ! Exclude Lower X Bound (Exclude cells a11, a12, a13 in the diagram).
  +          elxb: if (i > 1) then
  +            Diag    = Diag + scale
  +            OffDiag = OffDiag - scale * FldSlv % p(idx - 1)
  +          end if elxb
  +
  +          ! Exclude Lower X Bound (Exclude cells a13, a23, a33 in the diagram).
  +          euxb: if (i < FldSlv % whd(1)) then
  +            Diag    = Diag + scale
  +            OffDiag = OffDiag - scale * FldSlv % p(idx + 1)
  +          end if euxb
  +
  +          ! Exclude Upper Y Bound (Exclude cells a11, a12, 9 in the diagram).
  +          elyb: if (j > 1) then
  +            Diag    = Diag + scale
  +            OffDiag = OffDiag - scale * FldSlv % p(idx - FldSlv % whd(1))
  +          end if elyb
  +
  +          ! Exclude Lower Y Bound (Exclude cells 1, 4, 7 in the diagram).
  +          euyb: if (j < FldSlv % whd(2)) then
  +            Diag    = Diag + scale
  +            OffDiag = OffDiag - scale * FldSlv % p(idx + FldSlv % whd(1))
  +          end if euyb
  +
  +          NewP = (FldSlv % prhs(idx) - OffDiag) / Diag
  +
  +          Err = max(Err, abs(FldSlv % p(idx) - NewP))
  +
  +          FldSlv % p(idx) = NewP
  +
  +        end do lx
  +      end do ly
  +
  +      if (Err < Derr) return
  +    end do iterl
  +
  +  end subroutine GSPSlv2D
  +
     subroutine InitFluid(FldSlv, names, whd, scale)
       implicit none
       character(len=*),  dimension(:), intent(in) :: names
  -    integer           , intent(inout) :: whd(:) ! Width (i), height (j), depth (k)
  +    integer           , intent(in) :: whd(:) ! Width (i), height (j), depth (k)
       real(dp), optional, intent(in) :: scale
       class(FluidSlv),   intent(out) :: FldSlv ! Fluid
       real(dp), dimension(size(whd)) :: ost    ! Offset
  
positive.code:     class(FluidSlv), intent(inout) :: FldSlv
      integer, allocatable, save     :: ndim ! Number of dimensions.
  
      if (.not. allocated(ndim)) then
        allocate(ndim, source = size(FldSlv % u))
      end if
  
      if (ndim == 2) then
        call FldSlv % BuildRhs2D
      else
        call FldSlv % BuildRhs3D
      end if
    end subroutine BuildRhs
  
    subroutine BuildRhs2D(FldSlv)
      implicit none
      class(FluidSlv), intent(inout) :: FldSlv
      real(dp)                       :: scale
      integer                        :: i, j, idx
  
      ! Calculate the scale.
      scale =  1./FldSlv % celsiz
  
negative.code:     integer                     :: ix(2)
      real(dp)                    :: VtxVal(4) ! VALues at cell VERtices. VtxVal(1) = x00, VtxVal(2) = x10, VtxVal(3) = x01, VtxVal(4) = x11
      integer, parameter          :: ivtc(8) = [0,1,0,1,0,0,1,1] + 1 ! Index array of the cell vertex displacements.
  
      ! x01   -----------   x11
      !       |         |
      !       |         |
  