Is energy extensivity necessary in thermodynamics? Given a partition of a system into two smaller systems, the energy $U$ is devided into $U_1$ and $U_2$, with 
$$U=\mathcal{P}(U_1,U_2):=U_1+U_2,$$ 
so that $U_2$ is given by $U-U_1$. Here the operation $\mathcal{P}$ denotes the partition rule. 
What I've written above is basically the extensivity of energy and it somewhat enforces the idea of packages of energy, or particles. The idea of extensivity is closely related to additivity. A small energy flow from system one to system two is then given via 
$$U_1\rightarrow \mathbb{fl}^-(U_1):= U_1-\epsilon,
\\ U_2\rightarrow \mathbb{fl}^+(U_2):=U_2+\epsilon.$$
What we don't want to drop is energy conservation, so that 
$$U=\mathcal{P}(U_1,U_2)=\mathcal{P}(\mathbb{fl}^-(U_1),\mathbb{fl}^+(U_2)),$$
where $\mathbb{fl}^\pm$ define what a flow of energy is supposed to mean, and due to energy conservation these functions in part determined by $\mathcal{P}$.

Are there any good arguments why in thermodynamics or in classical statistical physics, the implicit laws for this kind of divistion should be formulated with a plus sign (other than the fact that there are good theories in which this idea works fulfill)?

After all, there are some models for which the notions of extensivity or additivity are not so clear for all quantities (e.g. some of the entropy-definitions of the last decades).

The relations
$$\mathcal{P}(U_2,U_1)=\mathcal{P}(U_1,U_2),$$
$$\mathcal{P}(\mathbb{fl}^-(U_1),\mathbb{fl}^+(U_2))=\mathcal{P}(\mathbb{fl}^+(U_1),\mathbb{fl}^-(U_2)),$$ 
also seem very natural - although this might not be totally necessary, as the values of $U_1$ and $U_2$ are elements of the reals and therefore ordered/distinguishable. 
Also, I think the order of putting subsystems together can't matter
$$\mathcal{P}(U_1,\mathcal{P}(U_2,U_3))=\mathcal{P}(\mathcal{P}(U_1,U_2),U_3).$$
There should be more of these kind of restrictions.

So for example, a first ansatz would be that one could consider $$U=\mathcal{P}(U_1,U_2):=U_1U_2,$$
so that 
$$\mathbb{fl}^-(U_1)\equiv\tfrac{U_2}{\mathbb{fl}^+(U_2)}U_1.$$
One realization of that (which is constructed such that the change via $\mathbb{fl}^-$ is the same as in the additive case) would be 
$$U_1\rightarrow \mathbb{fl}^-(U_1):= \left(1-\tfrac{\epsilon}{U_1}\right)U_1 =U_1-\epsilon,
\\ U_2\rightarrow \mathbb{fl}^+(U_2):= \left(1-\tfrac{\epsilon}{U_1}\right)^{-1}U_2 =U_2+\left(\tfrac{U_2}{U_1}\right)\epsilon+ O\left((\tfrac{\epsilon}{U_1})^2\right).$$
 A: From the statistical mechanics perspective, the fact that energy is conserved  forces it to be one of the extensive variables in any thermodynamic description.
A: I'm not sure to what extent this answers your question, but here are some of my thoughts on this. If energy were not extensive then it would not cause problems for the formalism of thermodynamics, but we would no longer have the fact that all temperatures become equal at equilibrium. Let's examine your example where $U = U_1 U_2$. We'll say the entropy is still extensive, so that $S = S_1(U_1) + S_2(U_2)$. Let's consider putting these two bodies in contact and allowing them to go to equilibrium. This happens when $$\frac{\partial S_1}{\partial U_1}dU_1 = \frac{\partial S_2}{\partial U_2}dU_2,$$
or
$$\frac{dU_1}{T_1} = \frac{dU_2}{T_2},$$
according to the definition of temperature. 
In regular thermodynamics we have that $dU_1 = -dU_2$, which allows us to conclude that $1/T_1 = 1/T_2$ when the bodies are in equilibrium. However, in this example we instead have that $dU_1/U_1 = dU_2/U_2$. (You can see this by noting that $\log(U_1) + \log(U_2)$ is the conserved quantity.) This means that, at equilibrium we have $U_1/T_2 = U_2/T_1$ instead of $T_1 = T_2$, so the temperatures won't be equal in general.
This sort of thing happens all the time in chemistry: if you're dealing with a reaction like $X\leftrightarrow Y$ then the two chemical potentials become equal at equilibrium, whereas for something like $X\leftrightarrow 2Y$ they do not, because $X+Y$ is not conserved.
On a more physical level, I should note that the extensivity of energy is an approximation, and it fails in particular for very small systems. This is because there can be energy in the interaction between the two systems. In particular, for two strongly interacting particles, it really doesn't make sense to try and divide the Hamiltonian $H$ into $H_1+H_2$. It's only as the systems become large that the internal energy terms dominate over the interaction terms and energy becomes approximately extensive. As long as the two systems are interacting there has to be an interaction term, i.e. $H=H_1 + H_2 + H_\text{interaction}$, but $H_\text{interaction}$ can be very small compared to the other terms. Extensivity can also fail to happen if gravitational attraction between the systems is important, which I think is ultimately why gravitationally-bound systems have negative heat capacities.
The situation for entropy is similar, in that it only becomes extensive when interaction terms can be neglected, and this tends to happen for macroscopically-sized systems near equilibrium. It's only when we start to reach this large-system limit that statistical mechanics starts to turn into thermodynamics. 
Finally, I want to note that thermodynamics could be defined to work just fine, even if energy didn't exist at all. The fundamental equation is usually written
$$
dU = TdS - pdV + \sum_i\mu_i dN_i + \text{terms for things like charge and momentum}
$$
but what you really get from statistical mechanics is
$$
dS = \lambda_U dU + \lambda_V dV + \sum_i \lambda_{N_i}dN_i + \dots,
$$
where $\lambda_U = 1/T$, $\lambda_V = p/T$, $\lambda_{N_i} = -\mu_i/T$, etc. I  regard this as much more fundamental, because it puts the "special" quantity $S$ on the left hand side and all the conserved quantities on the right. When thermodynamics is formulated this way, energy loses what appeared to be its special role, and it behaves exactly like all of the other conserved quantities. If energy didn't exist at all (but the microscopic dynamics still obeyed Liouville's theorem) then we would simply lose the first term on the right-hand side of this equation, but everything else would stay the same, as long as we worked in terms of $\lambda_V$ and $\lambda_{N_i}$ instead of $p$ and $\mu_i$.
A: What happened to U if U1=U2 ? You loose energy conservation...
Thermodynamics is not a fundamental theory. It can be derived from kinetic theory which is itself derived from mechanics. Energy in mechanics is an additive function, and unless you want to modify the rules for integrals (Chasles) it remains additive in kinetic theory and thermodynamics.
A: I did some research about your idea of non extensive statistical theory and have found a positive answer to your question, though the theory is quite involved. To make things simple: Gibbs theorem presupposes ergodicity. We have extensive evidence that systems far from equilibrium exhibit both self organisation at any scale and chaotic regions: a behaviour incomptible with ergodicity. To take into account these observations Tsallis has introduced a non extensive entropy which asymptoticaly matches the Gibbs Boltzmann entropy. The corresponding non equilibrium theory involves fractal dynamics and fractal operators.
You can find an excellent presentation of the theory in this paper:
http://arxiv.org/ftp/arxiv/papers/1203/1203.4003.pdf
