A three-dimensional simulation of the formation of metachronal waves in rows of pulmonary cilia is presented. The cilia move in a two-layer fluid model. The fluid layer adjacent to the cilia bases is purely viscous while the tips of the cilia move through a viscoelastic fluid. An overlapping fixed-moving grid formulation is employed to capture the effect of the cilia on the surrounding fluid. In contrast with immersed boundary methods, this technique allows a natural enforcement of boundary conditions without the need for smoothing of singular force distributions. The fluid domains are discretized using a finite volume method. The 9 + 2 internal microtubule structure of an individual cilium is modeled using large-deflection, curved, finite-element beams. The microtubule skeleton is cross-linked to itself and to the cilium membrane through spring elements which model nexin links. The cilium membrane itself is considered to be elastic and subject to fluid stresses computed from the moving grid formulation as well as internal forces transmitted from the microtubule skeleton. A cilium is set into motion by the action of dynein molecules exerting forces between adjacent microtubules. Realistic models of the forces exerted by dynein molecules are extracted from measurements of observed cilia shapes.