The parallelization of a fully implicit and stable finite element algorithm with relative low memory requirements for the accurate simulation of time-dependent, free-surface flows of multimode viscoelastic liquids is presented. It is an extension of our multi-stage sequential solution procedure which is based on the mixed finite element method for the velocity and pressure fields, an elliptic grid generator for the deformation of the mesh, and the discontinuous Galerkin method for the viscoelastic stresses [Dimakopoulos and Tsamopoulos [12,14]]. Each one of the above subproblems is solved with the Newton-Rapshon technique according to its particular characteristics, while their coupling is achieved through Picard cycles. The physical domain is graphically partitioned into overlapping subdomains. In the process, two different kinds of parallel solvers are used for the solution of the distributed set of flow and mesh equations: a multifrontal, massively parallel direct one (MUMPS) and a hierarchical iterative parallel one (HIPS), while viscoelastic stress components are independently calculated within each finite element. The parallel algorithm retains all the advantages of its sequential predecessor, related with the robustness and the numerical stability for a wide range of levels of viscoelasticity. Moreover, irrespective of the deformation of the physical domain, the mesh partitioning remains invariant throughout the simulation. The solution of the constitutive equations, which constitutes the largest portion of the system of the governing, non-linear equations, is performed in a way that does not need any data exchange among the cluster's nodes. Finally, indicative results from the simulation of an extensionally thinning polymeric solution, demonstrating the efficiency of the algorithm are presented. (C) 2010 Elsevier B.V. All rights reserved.