This work presents a method to efficiently determine the dominant Karhunen--Loève (KL) modes of a random process with known covariance function. The truncated KL expansion is one of the most common techniques for the approximation of random processes, primarily because it is an optimal representation, in the mean squared error sense, with respect to the number of random variables in the representation. However, finding the KL expansion involves solving integral problems, which tends to be computationally demanding. This work addresses this issue by means of a work-subdivision strategy based on a domain decomposition approach, enabling the efficient computation of a possibly large number of dominant KL modes. Specifically, the computational domain is partitioned into smaller nonoverlapping subdomains, over which independent local KL decompositions are performed to generate local bases which are subsequently used to discretize the global modes over the entire domain. The latter are determined by means of a Galerkin projection. The procedure leads to the resolution of a reduced Galerkin problem, whose size is not related to the dimension of the underlying discretization space but is actually determined by the desired accuracy and the number of subdomains. It can also be easily implemented in parallel. Extensive numerical tests are used to validate the methodology and assess its serial and parallel performance. The resulting expansion is exploited in Part B to accelerate the solution of the stochastic partial differential equations using a Monte Carlo approach.