A solution code of the Stokes flow system in three dimensions under a self-gravitating field was developed with a free surface in order to simulate the long-time scale evolution of the solid earth system on the massively parallel-vector super computer Earth Simulator 2 (ES2). The largely deformable surface of a planetary body is represented by a sharp boundary of scalar color function, which is advected by the CIP-CLSR method. The self-gravitating field is calculated by the Poisson equation of the gravity potential. A low viscous material (sticky air) surrounding the planet mimics its free surface motion. In order to solve the ill conditioned Stokes problem of the finite difference discretization on a staggered grid, I employ the iterative Stokes flow solver, robust to large viscosity jumps, using a strong Schur complement preconditioner and mixed precision arithmetic utilizing the double-double method. In this paper, I present an overview of the numerical algorithms and the current parallel programming strategy for the efficient performance on ES2. In addition, preliminary results of the core formation simulation are demonstrated in three dimensions as an example of the target Stokes flow system of this simulation code.