The main documentation of the Add_to_Matrix_Equation_Monomial Procedure contains additional explanation of this code listing.

```  subroutine Add_to_Matrix_Equation_Monomial (Monomial, ELLM, RHS_MV)

! Use association information.

use Caesar_Numbers_Module, only: one, zero

! Input variables.

type(Monomial_type), intent(inout) :: Monomial  ! Monomial to be added.

! Input/Output variables.

! ELL_Matrix to be incremented.
type(ELL_Matrix_type), intent(inout) :: ELLM
! RHS mathematic vector.
type(Mathematic_Vector_type), intent(inout) :: RHS_MV

! Internal variables.

! Arrays to manipulate values of the matrix.
type(real,1) :: Matrix_Values
type(integer,1) :: Matrix_Rows
type(integer,1) :: Matrix_Columns
type(integer) :: i                          ! Loop parameter.
type(integer) :: shift                      ! Index shift.
type(real,1) :: Volume_Cells                ! Volume of the cells.

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! Verify requirements.

VERIFY(Valid_State(Monomial),5) ! Monomial is valid.
VERIFY(Valid_State(ELLM),5)     ! ELLM is valid.
VERIFY(Valid_State(RHS_MV),5)   ! RHS_MV is valid.

! Create some working vectors.

call Initialize (Matrix_Values, NCells_PE(Monomial%Mesh))
call Initialize (Matrix_Rows, NCells_PE(Monomial%Mesh))
call Initialize (Matrix_Columns, NCells_PE(Monomial%Mesh))
call Initialize (Volume_Cells, NCells_PE(Monomial%Mesh))

! Set Matrix_Rows, Matrix_Columns.
!
! To spread out the rows and columns from this equation into a multiple
! equation matrix, translate local row numbers via:
!   i --> NEquations*(i - 1) + Equation
! and global column numbers via:
!   i + shift --> NEquations*(i + shift - 1) + Variable

shift = First_PE(Monomial%Structure) - 1
do i = 1, Length_PE(Monomial%Structure)
Matrix_Rows(i) = Monomial%NEquations*(i - 1) + Monomial%Equation
Matrix_Columns(i) = Monomial%NEquations*(i + shift - 1) + &
Monomial%Variable
end do

! Increment ELLM values.

call Get_Volume_Cells (Volume_Cells, Monomial%Mesh)
if (Monomial%Exponent /= zero) then
if (Monomial%Derivative) then
! Time derivative version.
Matrix_Values(:) = Monomial%Coefficient * Monomial%Exponent * &
Monomial%Phi_k ** (Monomial%Exponent - one) * &
Volume_Cells / Monomial%Delta_t
else
! Standard version.
Matrix_Values(:) = Monomial%Coefficient * Monomial%Exponent * &
Monomial%Phi_k ** (Monomial%Exponent - one) * &
Volume_Cells
end if
else
! Treat exponent=0 (monomial=constant) case separately to avoid
! division by zero if Phi_k=0 and to save a little time.
Matrix_Values(:) = zero
end if
call Add_Values (ELLM, Matrix_Values, Matrix_Rows, Matrix_Columns, &
Global=.false.)

! Increment RHS_MV values. Note that they have been multiplied by -1
! to move them to the RHS of the equation.

if (Monomial%Derivative) then
! Time derivative version.
Matrix_Values(:) = Monomial%Coefficient *  &
((Monomial%Exponent - one) * &
Monomial%Phi_k ** Monomial%Exponent + &
Monomial%Phi_n ** Monomial%Exponent) * &
Volume_Cells / Monomial%Delta_t
else
! Standard version.
Matrix_Values(:) = Monomial%Coefficient * &
(Monomial%Exponent - one) * &
Monomial%Phi_k ** Monomial%Exponent * &
Volume_Cells
end if
call Add_Values (RHS_MV, Matrix_Values, Matrix_Rows, Global=.false.)

! Finalize working vectors.

call Finalize (Matrix_Values)
call Finalize (Matrix_Rows)
call Finalize (Matrix_Columns)
call Finalize (Volume_Cells)

! Verify guarantees.

VERIFY(Valid_State(ELLM),5)     ! ELLM is valid.
VERIFY(Valid_State(RHS_MV),5)   ! RHS_MV is valid.

return