In order to improve the resolution of traveltime seismic tomography surveys, P and Swave traveltimes can be used complementary under a joint inversion scheme. Usually both P and S arrivals can be obtained with minimum additional effort; therefore joint inversion of these datasets offers an inexpensive and efficient way to improve their interpretation. A joint inversion algorithm of P and S traveltimes can maximize useful information from existing datasets and enhance the effectiveness of seismic refraction tomography surveys. In this work we present a joint inversion algorithm that inverts the two different datasets subject to a velocity ratio linking constraint. This constraint could be constant for all model parameters or could vary spatially, depending on the available information. The compressional to shear waves velocity ratio value can be roughly estimated from existing a-priori information, based on well logging or lab measurements or even can be estimated from geological information. A critical issue regarding joint inversion schemes is the significance of the constraint equations in the overall inversion procedure. If this constraint is too lose then the cross-correlation between the parameters of different type is negligible and the two models vary independently, degrading the joint inversion scheme in two separate inversions. On the contrary if too heavy weighting is applied, the solution is strongly biased towards the a-priori assumption, neglecting the actual data information. The participation of linking equation is controlled by a Lagrangian multiplier, which is usually defined empirically (i.e. extracted from a trial and error procedure) and is adopted as uniform for the whole model area. We introduce a spatially variable Lagrangian multiplier vector instead of a scalar one that scales differently the ratio constraint for each pair of parameters. The weighting is based on the parameter resolution matrices of the two models and spread function analysis. For highly resolvable parameters, a small value of the Lagrangian multiplier is assigned and the parameters are allowed to vary almost independently. For poorly resolved parameters, the Lagrangian multiplier that is assigned is large and due to the lack of information the specific parameters are forced to follow the ratio constraint. This method allows areas of the models with high information density values to vary based on this information rather than the ratio constraint equations, preserving the information content and the contrast that separate inversions would provide. At the same time the areas of the models which are lacking of resolving power and introduce instability under independent inversions, are constrained under the joint inversion scheme. Additional regularization is used and spatially variable smoothness is also applied to the models parameters. Testing with synthetic and real data suggest that the presented joint inversion algorithm can lead to improved results, stabilize the inversion process and reduce the non-uniqueness of the problem.