Abstract: Stream function and velocity potential can represent the vorticity and divergence of flow fields, respectively, and have been widely used in studies of global and regional atmospheric and oceanic circulations, pollutant diffusion, and data assimilation for a long time. In recent years, accuracy has decreased sharply when the commonly used algorithms are applied to complex flows driven by mesoscale and storm-scale weather systems, particularly in boundary layers over complex terrains and/or heterogeneous underlying surfaces. This paper presents a comprehensive review of the algorithms developed since the 1950s in five categories, emphasizing their strengths and weakness from the perspective of mathematical principles and physical meanings and summarizing their scopes of applications. The previously developed harmonic cosine series expansion spectral approach is revisited and corrected for its ill-positioned solvability conditions to improve its suitability and accuracy in solving complex flow fields. Numerical experiments based on both idealized and real flow fields are performed to illustrate/highlight and summarize the applicability and accuracy of the algorithms in each category for different types of datasets with different spatial resolutions. The objective of this paper is to provide a solid scientific basis for the correct and efficient applications of stream function, velocity potential, and their derived variables in the diagnostic analysis and numerical prediction of extreme weather and climate events.