title: HelmHoltz decomposition using Python and avoiding spherical harmonic calculations
body: <p>I am trying to decompose a wind field (u,v) into the irrotational and non-divergent fields (uchi, vchi, upsi, vpsi). This field is global, and in the form (time, level,lat,lon), and works when I use the HelmHoltz decomposition from the pyspharm python module and a field that does not have any missing values--which is the case for publicly available data (reanalysis data that has interpolated to the surface). Pyspharm uses spherical harmonic calculations to calculate this.</p>
  <p>However, the issue arises when I try to use data I derived from my own personal model where I have missing values-due to mountains or high regions, as I can't interpolate to the surface. Once there is a single NaN on a given level, the whole level returns a NaN.</p>
  <p>I initially thought of interpolating horizontally to smooth over these NaN, but the regions are often too large, so I thought perhaps there is a different way to do the Helmholtz decomposition that would return NaN only near the initial NaN (and not entire level), such as finite differences?</p>
  <p><strong>Does anyone know how or if I can perform a finite differences for the HelmHoltz, and if so where I can get the code for this?THX</strong></p>