Global sensitivity analysis (GSA) of functional-output models is usually performed by combining statistical techniques, such as basis expansions, metamodeling and sampling based estimation of sensitivity indices. By neglecting truncation error from basis expansion, two main sources of errors propagate to the final sensitivity indices: the metamodeling related error and the sampling-based, or pick-freeze (PF), estimation error. This work provides an efficient algorithm to estimate these errors in the frame of Gaussian processes (GP), based on the approach of Le Gratiet et al. [16]. The proposed algorithm takes advantage of the fact that the number of basis coefficients of expanded model outputs is significantly smaller than output dimensions. Basis coefficients are fitted by GP models and multiple conditional GP trajectories are sampled. Then, vector-valued PF estimation is used to speed-up the estimation of Sobol indices and generalized sensitivity indices (GSI). We illustrate the methodology on an analytical test case and on an application in non-Newtonian hydraulics, modelling an idealized dam-break flow. Numerical tests show an improvement of 15 times in the computational time when compared to the application of Le Gratiet et al. [16] algorithm separately over each output dimension.