assemble1D Subroutine

    subroutine assemble1D(points, cells, diff, vel, GlobalA)
        !* Assemble the global stiffness matrix based on element connectivity
        integer,  intent(in),   dimension(:,:)  :: cells    !! Element connectivity
        real(wp), intent(in)                    :: diff     !! Diffusivity coefficient [m/s^2]
        real(wp), intent(in)                    :: vel      !! Velocity [m/s]
        real(wp), intent(in),   dimension(:)    :: points   !! Array of nodal coordinates
        real(wp), intent(out),  dimension(:,:)  :: GlobalA  !! Global Stiffness matrix

        integer               :: ii, num_cells, num_pts
        real(wp), parameter   :: eps = epsilon(1e0)
        real(wp), allocatable :: xy(:), Ie(:,:)

        GlobalA   = 0._wp
        num_cells = size(cells, 1)

        ! Add elemental stiffness matrices to Global Stiffness Matrix
        !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(ii, xy, Ie) REDUCTION(+:GlobalA)
        !$OMP DO
        do ii = 1, num_cells

            num_pts   = size(cells(ii,:))

            ! Reallocate elemental stiffness matrix
            allocate(xy(num_pts), Ie(num_pts, num_pts))

            xy = points(cells(ii,:))

            Ie = assembleElementalMatrix1D(num_pts, 1, 1, xy)
            ! call r8mat_print(num_pts, num_pts, Ie, 'Elemental Stiffness Matrix:')

            GlobalA(cells(ii,:), cells(ii,:)) = &
                GlobalA(cells(ii,:), cells(ii,:)) - diff*Ie

            if ( abs(vel) .gt. eps ) then
                Ie = assembleElementalMatrix1D(num_pts, 0, 1, xy)
                ! call r8mat_print(num_pts, num_pts, Ie, 'Elemental Stiffness Matrix:')

                GlobalA(cells(ii,:), cells(ii,:)) = &
                    GlobalA(cells(ii,:), cells(ii,:)) - vel*Ie

            endif

            ! Deallocate elemental stiffness matrix after every loop
            deallocate(Ie, xy)
            ! stop
        enddo
        !$OMP END DO NOWAIT
        !$OMP END PARALLEL

        ! call r8mat_print(num_nodes, num_nodes, GlobalA, 'Global Stiffness Matrix:')
        ! stop
        return
    end subroutine assemble1D