The atmospheric and oceanic turbulent circulations involve a variety of nonlinearly interacting physical processes spanning a broad range of spatial and temporal scales. To make simulations of these turbulent flows computationally tractable, processes with scales smaller than the typical grid size of weather/climate models have to be parameterized. Recently, there has been substantial interest (and progress) in using deep learning techniques to develop data-driven subgrid-scale (SGS) parameterizations for the climate system. Another approach that is rapidly gaining popularity is to learn the entire spatio-temporal variability of the climate system from data, i.e., developing fully data-driven forecast models or emulators. For either of these approaches to be useful and reliable in practice, a number of major challenges have to be addressed. These include: 1) instabilities or unphysical drifts, 2) learning in the small-data regime, 3) interpretability, and 4) extrapolation to different parameters. Using several setups of 2D turbulence, two-layer quasi-geostrophic turbulence, Rayleigh-Benard convection, and ERA5 reanalysis, we introduce methods to address (1)-(4). The key aspect of some of these methods is combining the spectral analyses of deep neural networks and turbulence/nonlinear physics, as well as leveraging recent advances in theory and applications of deep learning. . In the end, we will discuss scaling up these methods to more complex systems and real-world applications, e.g., for SGS modeling of atmospheric gravity waves and conducting short- and long-term weather forecasting.