In Large-Eddy Simulation (LES), eddies that are bigger than the grid spacing are explicitly resolved and the sub-grid scale eddies are modelled. The region near the solid surface, where boundary effects are significant, play a critical role in the evolution of turbulence and therefore a high-resolution is desirable in the surface-layer. The number of grid points required to resolve the energy containing eddies is found in the literature to be in the range of Re1.76. Reynolds numbers of order 107 are commonly encountered in the atmospheric boundary layer flows and this would correspond to a total number of grid points of order 1012. In spite of the radical increase in the power of computers, LES of atmospheric flows with high-resolution in the surface-layer are not feasible yet. One approach to reduce the required number of grid points is by employing a grid nesting technique; i.e in the region of interest, a fine grid (FG) is embedded within the coarse grid (CG) parent domain. Numerical weather prediction models widely use horizontal grid nesting. However, vertical grid nesting is not offered in many LES models. We implement an interactive two-way vertical grid nesting in Parallelized Large-Eddy Simulation Model (PALM). It is a Finite Difference solver for the incompressible Boussinesq equations. It takes advantage of two-dimensional domain decomposition to achieve good-load balancing between the processing elements. We adopt MPI-1 standards for the communication between the CG and the FG so that the grid nesting feature will be supported by wide range of parallel compilers. The available processor cores are split between CG and FG depending on the estimated computational load for each grid. The use of lateral periodic boundary conditions are kept possible by maintaining the CG and FG to have the same horizontal extent. The bottom boundary of both grids represents the solid surface. For the FG, Monin-Obukhov similarity is assumed between the surface and the first grid point in the direction normal to the surface. We allow the vertical extent of the FG to be varied as needed. To set the top boundary condition for the FG based on the CG values, a sponge layer is introduced at the interface between the two grids. The CG prognostic variables in the overlapping region are updated/anterpolated based on FG values. Finally, We validate the LES with vertical nesting against a high-resolution LES under identical boundary and flow conditions. Our method allows higher resolution in the surface-layer without the need to have a quasi-uniform grid spacing in the entire domain, and consequently reduces the computational cost. Possible applications include virtual measurements in the surface-layer, e.g Eddy Covariance tower with height of order 10 m.