For quantitative predictions the tokamak plasma has to be modelled as a consistent entity. Distinguishing features of fusion plasma theory are the simultaneous importance of a large number of effects, essential multidimensionality in geometrical and velocity space, and a high degree of nonlinearity in the interaction between these effects. Only computational methods, mainly based on linearization, iteration and on procedures for solving huge systems of linear equations, are widely applicable, and provide quantitative `point results' of the `numerical experiment' type. The practical limitations of computer capacity and cost of computing time impose severe limits on details which can be taken into account in computational plasma models. For plasmas of tokamaks such as JET, the models are set up as initial boundaryvalue problems on several time scales, and composed of a cluster of interdependent computer codes. The basic magnetic field-plasma configuration is determined from a one-fluid (magnetohydrodynamic) theory in two dimensions. The macroscopic stability of these configurations is checked. For stable plasmas the secular evolution of a sequence of equilibria is computed by `transport codes'. In these the balance equations for the conservation quantities (mass, momentum, energy) are solved for fluxes between, and sources and sinks on, the magnetic surfaces (as determined from equilibrium). These are multifluid equations in one spatial dimension for all charged plasma species including impurities. Important source terms such as electromagnetic wave heating, injected particles, and fusion $\alpha $-particles must be calculated kinetically with at least one additional velocity-space coordinate. At least at the plasma boundary, neutral atoms cannot be neglected. Their distribution is calculated by Monte-Carlo methods in three spatial dimensions (and velocity space). The bulk plasma is usually surrounded by magnetic surfaces or field lines that cross solid walls. Here also, for the transport of charged particles, propagation in two dimensions must be calculated. The consistent combination of these major elements is considered.