def _initialize_background(self):
        """Set up background state (zonal flow and PV gradients)."""

        # background vel.
        if len(np.shape(self.U)) == 0:
          self.U = (self.U * np.ones((self.ny)))
        print(np.shape(self.U))
        self.set_U(self.U)

        # the meridional PV gradients in each layer
        self.Qy = (self.beta + np.gradient(np.gradient(self.U, self.dy), self.dy))[np.newaxis,...]


        # complex versions, multiplied by k, speeds up computations to pre-compute
        self.ikQy = np.expand_dims(self.Qy, axis=2) *  1j * self.k

        self.ilQx = 0.