docs: def _initialize_model_diagnostics(self):
          """Extra diagnostics for two-layer model"""
  
          self.add_diagnostic('entspec',
              description='barotropic enstrophy spectrum',
              function= (lambda self:
                        np.abs(self.del1*self.qh[0] + self.del2*self.qh[1])**2.)
          )
  
          self.add_diagnostic('APEflux',
              description='spectral flux of available potential energy',
              function= (lambda self:
                self.rd**-2 * self.del1*self.del2 *
                np.real((self.ph[0]-self.ph[1])*np.conj(self.Jptpc)) )
          )
  
          self.add_diagnostic('KEflux',
              description='spectral flux of kinetic energy',
              function= (lambda self:
                np.real(self.del1*self.ph[0]*np.conj(self.Jpxi[0])) +
                np.real(self.del2*self.ph[1]*np.conj(self.Jpxi[1])) )
          )
          
          self.add_diagnostic('APEgenspec',
              description='spectrum of APE generation',
              function= (lambda self: self.U[:,np.newaxis] * self.rd**-2 * self.del1 * self.del2 *
                         np.real(1j*self.k*(self.del1*self.ph[0] + self.del2*self.ph[1]) *
                                    np.conj(self.ph[0] - self.ph[1])) )
          )
  
          self.add_diagnostic('APEgen',
              description='total APE generation',
              function= (lambda self: self.U * self.rd**-2 * self.del1 * self.del2 *
                         np.real((1j*self.k*
                              (self.del1*self.ph[0] + self.del2*self.ph[1]) *
                              np.conj(self.ph[0] - self.ph[1])).sum()
                                +(1j*self.k[:,1:-2]*
                              (self.del1*self.ph[0,:,1:-2] + self.del2*self.ph[1,:,1:-2]) *
                              np.conj(self.ph[0,:,1:-2] - self.ph[1,:,1:-2])).sum()) /
                              (self.M**2) )
          )