A fully- or semi-implicit nonhydrostatic model requires that a solution to a Helmholtz equation be obtained at every time step. In the terrain-following coordinates employed by most models the Helmholtz operator assumes a general non-separable form, even when the model state is arbitrarily close to one of stably stratified rest, and the direct application of efficient Fourier transform methods is precluded. We examine a solution strategy in which the Helmholtz forcing is first interpolated to a Cartesian terrain- intersecting grid, so that the interior solution operator very closely approximates a separable, symmetric and horizontally homogeneous operator of algebraically simple and constant coefficients, which a Fourier method then efficiently solves. The reconciliation of the original Neumann-type boundary conditions over the terrain surface is then accomplished, to a very close approximation at least, by the construction, via a surface transform of the diagnosed Neumann residual from a preliminary trial solution, of a corrective forcing. A second iteration of the Fourier solver will now produce a solution in the Cartesian grid essentially consistent with the terrain boundary condition, so that a further interpolation back to the model's native grid achieves the original objective. By permitting the vertical portion of the grid for the solver to differ from that of the model we are free to choose from a much more general class of model vertical coordinates than would be the case if the choice were to be dictated by the requirement of numerical efficiency of the Helmholtz solver constrained to use this same grid. Implications for efficient semi-Lagrangian modeling will be discussed.