Our generalized neural closure models (gnCMs) based on unified neural partial differential equations (PDEs) are applied to ocean, sea ice, and chaotic systems. We augment existing/low-fidelity dynamical models directly in their PDE forms with both Markovian and non-Markovian neural network (NN) closures. The melding of the existing models with NNs in the continuous spatiotemporal space followed by numerical discretization automatically allows for generalizability. The Markovian term is designed to enable extraction of its analytical form and thus provides interpretability. The non-Markovian terms allow accounting for inherently missing time delays needed to represent the real world. Our flexible gnCMs provide full autonomy for the design of the unknown closure terms such as using any linear-, shallow-, or deep-NN architectures, selecting the span of the input function libraries, and using either or both Markovian and non-Markovian closure terms, all in accord with prior knowledge. We apply the gnCMs to learning experiments with advecting nonlinear waves, shocks, ocean acidification, ocean submesoscales, and sea ice models. We highlight applications to chaotic systems, emphasizing the need for adaptive learning schemes. Our learned gnCMs discover missing chaotic physics, find leading numerical error terms, discriminate among candidate functional forms in an interpretable fashion, achieve generalization, and compensate for the lack of complexity in simpler models.