The numerical aspects of physical parametrization are discussed mainly in the context of the ECMWF Integrated Forecasting System. Two time integration techniques are discussed. With parallel splitting the tendencies of all the parametrized processes are computed independently of each other. With sequential splitting, tendencies of the explicit processes are computed first and are used as input to the subsequent implicit fast process. It is argued that sequential splitting is better than parallel splitting for problems with multiple time scales, because a balance between processes is obtained during the time integration. It is shown that sequential splitting applied to boundary layer diffusion in the ECMWF model leads to much smaller time truncation errors than does parallel splitting. The so called Semi-Lagrangian Averaging of Physical Parametrizations (SLAVEPP), as implemented in the ECMWF model, is explained. The scheme reduces time truncation errors compared to standard first order methods, although a few implementation questions remain. In the scheme fast and slow processes are handled differently and it remains a research topic to find the optimal way of handling convection and clouds. Process specific numerical issues are discussed in the context of the ECMWF parametrization package. Examples are the non-linear stability problems in the vertical diffusion scheme, the stability related mass flux limit in the convection scheme and the fast processes in the cloud microphysics. Vertical resolution in the land surface scheme is inspired by the requirement to represent diurnal to annual time scales. Finally, a coupling strategy between atmospheric models and land surface schemes is discussed. It allows for fully implicit coupling also for tiled land surface schemes.