A two-dimensional algorithm to determine the steady state thermal structure of the lithosphere that integrates thermal, gravity, and local isostasy analyses is presented. Gravity analyses together with seismic data are used to constrain spatial variations in density and crustal structure, while absolute elevation is used to determine the lithospheric mantle thickness. The calculation is performed using a finite element technique that links the different physical equations. The program optionally calculates the temperature at any material boundary and, with given rheological parameters, the strength distribution and the total lithospheric strength in selected columns. We apply the algorithm to the Northeastern Spanish Geotransect which extends from the Pyrenees to the Balearic Promontory and along which a strong variation in crustal and lithospheric thickness is evident. We assess the use of two different inferred density models for the lithospheric mantle: The first assumes a linear decrease in density with increasing temperature using the asthenospheric density as a reference; the second model assumes a constant density for the whole lithospheric mantle. Although conceptually the two hypotheses differ substantially, the results obtained do not show significant differences. Lithospheric thicknesses of 120--130 km below the Pyrenees, 60--65 km in the Valencia Trough, and 65--75 km below the Balearic Promontory are deduced. In all cases the mean lithospheric mantle density has to be 40--60 kg m-3 higher than the asthenospheric density. The algorithm is shown to be a powerful tool in lithospheric thermal modeling especially in areas where surface heat flow is poorly constrained because of the temperature-density-elevation relationship. ¿ American Geophysical Union 1994