The sub-crustal stress has been traditionally computed using the Runcorn's formulae. This method allows computing the stress field only with a very limited spectral resolution. To overcome this problem, we apply a new method of computing the sub-crustal stress components based on utilizing the stress function with a subsequent numerical differentiation. This method increases the (degree-dependent) convergence domain of the asymptotically-convergent series and consequently allows evaluating the stress components to a higher spectral resolution compatible with currently available global crustal models. This method also facilitates the variable Moho geometry, instead of assuming only a constant Moho depth in the Runcorn's formulae. The crustal thickness and the sub-crustal stress are then determined directly from gravity and (seismic) crustal structure models. The numerical result reveals that the largest intensity of the sub-crustal stress occurs mainly along seismically active convergent tectonic plate boundaries, particularly along oceanic subduction zones and continent-to-continent collision zones.