assemble2D Module Subroutine

    module subroutine assemble2D(points, cells, diff, vel, GlobalA)
        integer,  intent(in),   dimension(:,:)  :: cells
        real(wp), intent(in)                    :: diff
        real(wp), intent(in),   dimension(2)    :: vel
        real(wp), intent(in),   dimension(:,:)  :: points
        real(wp), intent(out),  dimension(:,:)  :: GlobalA

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

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

        allocate(xy(num_pts, 2), Ie(num_pts, num_pts))

        ! NOTE: This is assuming that all elements have the same number of points
        !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(ii, xy, Ie) REDUCTION(+:GlobalA)
        !$OMP DO
        outer: do ii = 1, num_cells

            xy = points(cells(ii,:), :)

            ! *** Elemental matrix for diffusion in X and Y ***

            Ie  = assembleElementalMatrix(num_pts, 1, 1, xy) + &
                & assembleElementalMatrix(num_pts, 2, 2, xy)

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

            ! Add elemental matrix for velcity in X (1) and Y (2) to GlobalA
            adv: do jj = 1, 2
                if ( abs(vel(jj)) .gt. eps ) then

                    Ie = assembleElementalMatrix(num_pts, 0, jj, xy)

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

                endif
            enddo adv
        enddo outer
        !$OMP END DO NOWAIT
        !$OMP END PARALLEL

        deallocate(Ie, xy)
        return
    end subroutine assemble2D