anchor.diff:      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
  +    ! 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
  +    else C2D3D
         call FldSlv % BuildRhs3D
  -    end if
  +    end if C2D3D
     end subroutine BuildRhs
   
     subroutine BuildRhs2D(FldSlv)
       implicit none
       class(FluidSlv), intent(inout) :: FldSlv
  -    real(dp)                       :: scale
  -    integer                        :: i, j, idx
  +    real(dp)                       :: scale     ! Scaling factor.
  +    integer                        :: i, j, idx ! Counters and index.
   
       ! Calculate the scale.
       scale =  1./FldSlv % celsiz
  
positive.code: 
  contains
  
  
    pure function Lerp1D(FldQty, yi, x)
      ! Linear intERPolate in 1D.
      ! 1D interpolation between yi(1) < y < yi(4) for 0 < x < 1
  
negative.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(1) )
            ! Write value to old value.
            if ( abs(FldQty % ReadVal(i, j, k)) < abs(value) ) call FldQty % WriteOldVal(i, j, k, value)
          end do lx
  