Recently, coupled tensor decomposition has been widely used in data fusion of a hyperspectral image (HSI) and a multispectral image (MSI) for hyperspectral super-resolution (HSR). However, exsiting works often ignore the inherent non-negative (NN) property of the image data, or impose the NN constraint via hard-thresholding which may interfere with the optimization procedure and cause the method to be sub-optimal. As such, we propose a novel NN coupled canonical polyadic decomposition (NN-C-CPD) algorithm, which makes use of the parametric method and nonlinear least squares (NLS) framework to impose the NN constraint into the C-CPD computation. More exactly, the NN constraint is converted into the squared relationship between the NN entries of the factor matrices and a set of latent parameters. Based on the chain rule for deriving the derivatives, the key entities such as gradient and Jacobian with regards to the latent parameters can be derived, thus the NN constraint is naturally integrated without interfering with the optimization procedure. Experimental results are provided to demonstrate the performance of the proposed NN-C-CPD algorithm in HSR applications.