Calculations using the cluster variation method (CVM) in the tetrahedron approximation have been applied to BCC Ising systems with spin number up to s=15/2. Ferromagnetic first neighbour interactions (J(1)0) of variable strength α=−J(2)/J(1) have been analysed. For small α values the system is completely ferromagnetic, showing a second order transition from the ferromagnetic to the paramagnetic state. The Curie temperature TC depends on the value of α. For a critical value αC a multicritical point appears at TC in the α×T phase diagram. For α values higher than αC, two transitions occur: at low temperatures, a first order transition from the ferromagnetic ground state to the antiferromagnetic state, followed by a second order transition from the antiferromagnetic to the paramagnetic state. The specific heat capacities calculated for these model systems present a remarkable resemblance to experimental data for HCP Tb (s=6) and Dy (s=15/2). These results suggest that the complex magnetic structures observed in these metals should be interpreted as an effect of the high spin number statistics.