docs: def _initialize_background(self):
          """Set up background state (zonal flow and PV gradients)."""
  
          self.H = self.Hi.sum()
          if np.asarray(self.U).ndim == 2:
              self.Ubg = self.U * np.ones((self.ny))
          else:
              self.Ubg = np.expand_dims(self.U,axis=1) * np.ones((self.ny))
  
  
          if not (self.nz==2):
              self.gpi = self.g*(self.rhoi[1:]-self.rhoi[:-1])/self.rhoi[:-1]
              self.f2gpi = (self.f2/self.gpi)[:,np.newaxis,np.newaxis]
  
              assert self.gpi.size == self.nz-1, "Invalid size of gpi"
  
              assert np.all(self.gpi>0.), "Buoyancy jump has negative sign!"
  
              assert self.Hi.size == self.nz, self.logger.error('size of Hi does not' +
                       'match number of vertical levels nz')
  
              assert self.rhoi.size == self.nz, self.logger.error('size of rhoi does not' +
                       'match number of vertical levels nz')
  
              assert self.Ubg.size == self.nz * self.ny, self.logger.error('size of Ubg does not' +
                       'match number of vertical levels nz')
  
              assert self.Vbg.size == self.nz, self.logger.error('size of Vbg does not' +
                       'match number of vertical levels nz')
  
          else:
              self.f2gpi = np.array(self.rd**-2 *
                  (self.Hi[0]*self.Hi[1])/self.H)[np.newaxis]
  
  
          self._initialize_stretching_matrix()
  
          # the meridional PV gradients in each layer
          self.Qy = self.beta - np.dot(self.S, self.Ubg) + np.gradient(np.gradient(self.Ubg, self.dy, axis=1), self.dy, axis=1)
          self.Qx = np.dot(self.S,self.Vbg)
  
  
          # complex versions, multiplied by k, speeds up computations to precompute
          self.ikQy = self.Qy[:,:,np.newaxis]*1j*self.k
          self.ilQx = self.Qx[:,np.newaxis,np.newaxis]*1j*self.l