We propose a numerical self-consistent method for 3D classical lattice mode
ls, which optimizes the variational state written as a two-dimensional prod
uct of tensors. The variational partition function is calculated using the
corner transfer matrix renormalization group (CTMRG), which is a variant of
the density matrix renormalization group (DMRG). The numerical efficiency
of the method is exemplified in its application to the 3D Ising model.