How to constraint the Gaussian process model output ranging between 0 and 1? It requires the Gaussian process model output ranging between 0 and 1, is there any technique to constrain the Gaussian process model's output?

The simplest thing to do in such a case is to model a transformed version of $y$ by a Gaussian process. This is not specific to Gaussian processes; transformation is a standard instrument in the statistical toolbox.
Although many appropriate transformations exist, a good default is the logit transform. In stead of modelling
$$y(\cdot) \sim GP(m(\cdot), c(\cdot, \cdot))$$
Let $w = \text{logit}\,(y) = \log(y) - \log(1 - y)$.  You should then model
$$w(\cdot) \sim GP(m(\cdot), c(\cdot, \cdot))$$
and use the inverse-logit $y = \text{logit}^{-1}\,(w)  = \frac{1}{1 + \exp(-w)}$ to obtain predictions for $y$.
