Abstract. The freeze out (FO) of the expanding systems, created in relativistic heavy ion collisions, is discussed. We start with kinetic FO model, which realizes complete physical FO in a layer of given thickness, and then combine our gradual FO equations with Bjorken type system expansion into a unified model. We shall see that the basic FO features, pointed out in the earlier works, are not smeared out by the expansion.PACS. 51.10.+y -24.10. At highest energies available nowadays in relativistic heavy ion collisions at RHIC the total number of the produced particles exceeds 6000, therefore one can expect that the produced system behaves as a "matter" and generates collective effects. Indeed strong collective flow patterns have been measured at RHIC, which suggests that the hydrodynamical models are well justified during the intermediate stages of the reaction: from the time when local equilibrium is reached until the freeze out (FO), when the hydrodynamical description breaks down. During this FO stage, the matter becomes so dilute and cold that particles stop interacting and stream towards the detectors freely, their momentum distribution freezes out. The FO stage is essentially the last part of a collision process and the main source for observables.Nowadays, FO is usually simulated in two extreme ways: A) FO on a hypersurface with zero thickness, B) FO described by volume emission model or hadron cascade, which in principle requires an infinite time and space for a complete FO. At first glance it seems that one can avoid troubles with FO modeling using hydro+cascade two module model [1], since in hadron cascades gradual FO is realized automatically. However, in a such a scenario there is an uncertain point, actually uncertain hypersurface, where one switchs from hydrodynamical to kinetic modeling. First of all it is not clear how to determine such a hypersurface. This hypersurface in general may have both time-like and space-like parts. Mathematically this problem is very similar to hydro to FO phase transition on the infinitely narrow FO hypersurface, therefore for example all the problems discussed for FO on the hypersurface with space-like normal vectors will take place here. Another complication is that while for the post FO domain we have mixture of non-interacting ideal gases, now for the hadron cascade we should generate distributions for the interacting hadronic gas of all possible species, as a starting point for the further cascade evolution. The volume emission models are based on the kinetic equations [2,3] defining the evolution of the distribution functions, and therefore these also require to generate initial distribution functions for the interacting hadronic species on some hypersurface.In this work we present FO model which allows us to study FO in a layer of any thichness, L, from 0 to ∞, and which connects the pre FO hydrodynamical quantities, like energy density, e, baryon density, n, with post FO distribution function in a relatively simple way. Many building blocks of the model are ...