We report on the analysis of the first complete far-infrared spectrum (43-197 μm) of the Seyfert 2 galaxy NGC 1068 as observed with the Long Wavelength Spectrometer (LWS) on board the Infrared Space Observatory (ISO). In addition to the seven expected ionic fine-structure emission lines, the OH rotational lines at 79, 119, and 163 μm were all detected in emission, which is unique among galaxies with full LWS spectra, where the 119 μm line, when detected, is always in absorption. The observed line intensities were modeled together with ISO Short Wavelength Spectrometer (SWS) and optical and ultraviolet line intensities from the literature, considering two independent emission components: the active galactic nucleus (AGN) component and the starburst component in the circumnuclear ring of ∼3 kpc in size. Using the UV to mid-IR emission line spectrum to constrain the nuclear ionizing continuum, we have confirmed previous results: a canonical power-law ionizing spectrum is a poorer fit than one with a deep absorption trough, while the presence of a "big blue bump" is ruled out. Based on the instantaneous starburst age of 5 Myr constrained by the Brγ equivalent width in the starburst ring, and starburst synthesis models of the mid- and far-infrared fine-structure line emission, a low-ionization parameter (U = 10-3.5) and low densities (n = 100 cm-3) are derived. Combining the AGN and starburst components, we succeeded in modeling the overall UV to far-IR atomic spectrum of NGC 1068, reproducing the line fluxes to within a factor of 2.0 on average with a standard deviation of 1.3, and the overall continuum as the sum of the contribution of the thermal dust emission in the ionized and neutral components. The OH 119 μm emission indicates that the line is collisionally excited and arises in a warm and dense region. The OH emission has been modeled using spherically symmetric, nonlocal, non-LTE radiative transfer models. The models indicate that the bulk of the emission arises from the nuclear region, although some ...