anchor.diff:      ! Set up the simulation boundaries.
       bounds: do i = 1, 2
         j = 2*i
  -      ixi(j-1 : j) = floor( xi%x(j-1 : j) / FldQty % celsiz - FldQty % ost(i) ) + 1
  +      ixi(j-1 : j) = floor( xi % x(j-1 : j) / FldQty % celsiz - FldQty % ost(i) ) + 1
       end do bounds
   
       ! Loop through y.
  
positive.code:     ly: do j = max(ixi(3), 1), min( ixi(4), FldQty % whd(2) )
        ! Loop through x.
        ! Start from the maximum value between the desired x0 and 1 (for the cell index) up to the minimum value between the desired x1 and the simulation width.
        lx: do i = max(ixi(1), 1), min( ixi(2), FldQty % whd(2) )
          ! Write value to old value.
          if ( abs(FldQty % ReadVal(i, j)) < abs(value) ) call FldQty % WriteOldVal(i,j,value)
        end do lx
  
negative.code:     class(FluidSlv), intent(inout) :: FldSlv
      integer, allocatable, save     :: ndim ! Number of dimensions.
  
      ! Allocate ndim and set its value to the number of dimensions.
      if (.not. allocated(ndim)) allocate(ndim, source = size(FldSlv % u))
  
      ! Choose 3D or 2D.
      C2D3D: if (ndim == 2) then
        call FldSlv % BuildRhs2D
      else C2D3D
        call FldSlv % BuildRhs3D
      end if C2D3D
    end subroutine BuildRhs
  
    subroutine BuildRhs2D(FldSlv)
      implicit none
      class(FluidSlv), intent(inout) :: FldSlv
      real(dp)                       :: scale     ! Scaling factor.
      integer                        :: i, j, idx ! Counters and index.
  
      ! Calculate the scale.
      scale =  1./FldSlv % celsiz
  