A mathematical model is constructed for land glaciers with the thickness much less than the horizontal dimensions and radii of curvature of large bottom irregularities by means of the method of a thin boundary layer in dimensionless orthogonal coordinates. The dynamics are described by a statically determinate system of equations, so the solution for stresses is found. For the general non-isothermal case the interrelated velocity and temperature distributions are calculated by means of the iteration of solutions for velocity and for temperature. Temperature distribution is determined by a parabolic equation with a small parameter at the senior derivative. Its solution is reduced to the solution of a system of recurrent non-uniform differential equations of the first order by means of a series expansion of the small parameter. A relatively thin conducting boundary layer adjoins the upper and lower surfaces of a glacier, playing the role of a temperature damper in the ablation area. For ice divides, the statically indeterminate problem is solved, so the result for stresses depends on the temperature distribution.