Merge pull request NCAR#8 from climbfuji/update_gsd_develop_from_mast…

…er_20201001

Update gsd/develop from NOAA-EMC dev/emc 2020/10/01

docs/AdvectionOperators.md

@@ -0,0 +1,19 @@
+A  Brief Guide to Advection Operators in FV³ {#advection}
+========================
+
+By:
+
+**Lucas Harris, Shian-Jiann Lin, and Xi Chen, GFDL**
+
+*17 August 2018*
+
+FV³ is unique among atmospheric dynamical cores in that it is formulated similarly to CFD solvers which treat most processes as the fluxes of some conserved quantity. With the exception of the pressure gradient force, all of the terms in the inviscid Euler equations, formulated along Lagrangian surfaces, can be expressed as fluxes. The fluxes are evaluated through the Lin & Rood (1996) FV advection scheme, extended to the cubed sphere by Putman and Lin (2007), based on a two-dimensional combination of one-dimensional flux operators. The operators in the most recent version of FV³ all use the piecewise-parabolic method (Collella and Woodward 1984), although this method gives great flexibility in choosing the subgrid reconstructions used to evaluate the fluxes.FV³
+
+Some understanding of finite-volume methods is necessary to understand the operators. A true finite-volume method treats all variables as cell-mean values, instead of point values. (Fluxes and winds are treated as face-mean values.) To compute the fluxes, a ​reconstruction​ function describing an approximation to the (unknown) subgrid distribution of each variable is needed. In PPM, this reconstruction is a quadratic function, specified by the cell mean value and the values on the interfaces of the cells. The cell-interface values are derived by computing a high-order interpolation from the cell-mean values to point values on the interfaces. The reconstructions can then be altered by various means to enforce different shape preservation properties, such as monotonicity or positivity, by looking at the cell-mean and cell-interface values of the surrounding cells. This process, called ​limiting​ is a ​non-linear​ method, allowing us to create a high-order monotone solution without violating Godunov’s Theorem. (This should not be confused with flux-corrected transport, or FCT, a simpler means of enforcing monotonicity)
+
+Monotonicity means that the reconstructions are chosen to ensure that no new extrema appear in the variable, in the absence of flow convergence or divergence. This is a property of scalar advection in real flows, and prevents the appearance of unphysical oscillations, or over- and under-shoots (which can create spurious condensates or unphysically rapid chemical reactions). Monotonicity also prevents the appearance of negative values in positive-definite mass scalars, which can rapidly lead to numerical instability. However, the stronger the enforcement of monotonicity, the more diffusive the solution is, so other techniques, such as schemes which only enforce positivity, or simple reconstruction filters.
+
+Here we briefly describe three PPM operators, all formally the same fourth-order accuracy but with different reconstruction limiters: An unlimited (also called linear) “fifth-order” operator (hord = 5), an unlimited operator with a 2dx filter (hord = 6), and the monotone Lin 2004 operator (hord = 8). We also describe the optional positive-definite limiter which can be applied to methods 5 and 6, which is essential for scalar advection (hord_tr). ​The number indicated by hord is only the selection amongst a series of schemes, ​similar to how different numbers in the physics correspond to different schemes. ​They do not change the order of accuracy of the advection, only the diffusivity and shape-preserving characteristics.Hord = 5 is the “linear” or “unlimited” scheme originally devised by Collella and Woodward, modified with a weak 2dx filter: if the cell-interface values show a 2dx signal (that is, that the interpolated interface values on opposing sides of a cell switch between greater than and less than the cell-mean value) the flux is evaluated as a first-order upwind value. The 2dx filter is necessary to suppress 2dx noise; without this filter the solver is more prone to instability. Hord = 6 uses a much stronger 2dx filter: the hord = 5 method is extended by reverting to first-order upwind flux if the difference in cell-interface values exceeds the mean of the two interface values by a tunable threshold (1.5x by default).
+
+The optional positive-definite limiter is activated by adding a minus sign to hord (for values <= 6). This forces all cell-interface values to be no more than 0, ensuring a positive-definite flux. The positive-definite hord = 5 further applies the positive-definite constraint of Lin and Rood (1996).
+

docs/Chapter1.md

@@ -0,0 +1,35 @@
+FV³ Introduction {#introduction}
+============
+
+##Chapter 1
+
+###1.1 A brief history of FV³
+
+
+FV³ is a natural evolution of the hydrostatic Finite-Volume dynamical core (FV core) originally developed in the 90s on the latitude-longitude grid. The FV core started its life at NASA/Goddard Space Flight Center (GSFC) during early and mid 90s as an offline transport model with emphasis on the conservation, accuracy, consistency (tracer to tracer correlation), and efficiency of the transport process. The development and applications of monotonicity preserving Finite-Volume schemes at GSFC were motivated in part by the need to have a solution for the noisy and unphysical negative water vapor and chemical species (L94, LR96). It subsequently has been used by several high-profile Chemistry Transport Models (CTMs), including the NASA community GMI model (Rotman et al., 2001, J. Geophys. Res.), GOCART (Chin et al., 2000, J. Geophys. Res.), and the Harvard University-developed GEOS-CHEM community model (Bey et al, 2001, , J. Geophys. Res.). This transport module has also been used by some climate models, including the ECHAM5 AGCM (Roeckner et al, 2003, MPI-Report No. 349). Motivated by the success of monotonicity-preserving FV schemes in CTM applications, a consistently-formulated shallow-water model was developed. This solver was first presented at the 1994 PDE on the Sphere Workshop, and published in LR97. The Lin-Rood algorithm for shallow-water equations maintains mass conservation and a key Mimetic property of “no false vorticity generation”, and for the first time in computational geophysical fluid dynamics, uses high-order monotonic advection consistently for momentum and all other prognostic variables, instead of the inconsistent hybrid finite difference and finite-volume approach used by practically all other “finite volume” models today. The full 3D hydrostatic dynamical core, the FV core, was constructed based on the LR96 transport algorithm and the Lin-Rood shallow-water algorithm (LR97). The pressure gradient force is computed by the L97 finite-volume integration method, derived from Green’s integral theorem and based directly on first principles, and demonstrated errors an order of magnitude smaller than other pressure-gradient schemes at the time. The vertical discretization is the “vertically Lagrangian” scheme described by L04, the most powerful aspect of FV³, which permits great computational efficiency as well as greater accuracy given the same vertical resolution. The FV core was implemented in the NCAR CESM model in 2001 (Rasch et al, 2006, J. Climate) and in the GFDL CM2 model in 2004 (Delworth et al, 2006, J. Climate).
+
+
+The need to scale to larger number of processors in modern massively parallel environments and the scientific advantage of eliminating the polar filter led to the development of FV&sup3;, the cubed-sphere version of FV. The cubed-sphere FV&sup3; is in use in all GFDL and GSFC global models, and the cubed-sphere version of the LR96 advection scheme is used by the most recent version of GEOS-CHEM. Most recently, a computationally efficient non-hydrostatic solver has been implemented using a traditional semi-implicit approach for treating the vertically propagating sound waves. A second option for the non-hydrostatic solver, using a Riemann solver to nearly exactly solve for vertical sound-wave propagation, is also available. The Riemann solver is highly accurate and is very efficient if the Courant number for vertical sound wave propagation is small, and so may be very useful for extremely high (< 1 km) horizontal resolutions. In July 2016 this non-hydrostatic core has been selected for the Next-Generation Global Prediction System (NGGPS), paving the way for the unification of not only weather and climate models but also potentially regional and global models.
+
+
+###1.2 Outline of the solver
+
+FV&sup3; is the solver that integrates the compressible, adiabatic Euler equations in a weather or climate model. The solver is modular and designed to be called as a largely independent component of a numerical model, consistent with modern standards for model design; however for best results it is recommended that a model using FV&sup3; as its dynamical core should use the provided application programming interface (API) to invoke the solver as well as to use the provided utility routines consistent with the dynamics, particularly for the initialization, updating the model state by time tendencies from the physics, and for incorporating increments from the data assimilation system.
+
+The leftmost column of Figure 1.1 shows the external API calls used during a typical process-split model integration procedure. First, the solver is called, which advances the solver a full “physics” time step. This updated state is then passed to the physical parameterization package, which then computes the physics tendencies over the same time interval. Finally, the tendencies are then used to update the model state using a forward-in-time evaluation consistent with the dynamics, as described in Chapter 3.
+
+There are two levels of time-stepping inside FV&sup3;. The first is the “remapping” loop, the orange column in Figure 1.1. This loop has three steps:
+-# Perform the Lagrangian dynamics, the loop shown in the green column of Figure 1.1
+-# Perform the sub-cycled tracer advection along Lagrangian surfaces,using accumulated mass fluxes from the Lagrangian dynamics. Subcycling is done independently within each layer to maintain local (within each layer) stability.
+-# Remap the deformed Lagrangian surfaces on to the reference, or “Eulerian”, coordinate levels.
+
+
+
+<img src="../image/FV3flowchart.png" alt="drawing" width="800"/>  
+
+**Figure 1.1:**  *FV3 structure, including subroutines and time-stepping. Blue represents external API routines, called once per physics time step; orange routines are called once per remapping time step; green routines once per acoustic time step.*
+
+This loop is typically performed once per call to the solver, although it is possible to improve the model’s stability by executing the loop (and thereby the vertical remapping) multiple times per solver call.
+
+The Lagrangian dynamics is the second level of time-stepping in FV&sup3;. This is the integration of the dynamics along the Lagrangian surfaces, across which there is no mass transport. Since the time step of the Lagrangian dynamics is limited by horizontal sound-wave processes, this is called the “acoustic” time step loop. (Note that the typical assumption that the advective wind speed is much slower than the sound wave speed is often violated near the poles, since the speed of the polar night jets can exceed two-thirds of the speed of sound.) The Lagrangian dynamics has two parts: the C-grid winds are advanced a half-time step, using simplified (but similarly constructed) core routines, which are then used to provide the advective fluxes to advance the D-grid prognostic fields a full time step. The integration procedure is similar for both grids: the along-surface flux terms (mass, heat, vertical momentum, and vorticity, and the kinetic energy gradient terms) are evaluated forward-in-time, and the pressure-gradient force and elastic terms are then evaluated backwards-in-time, to achieve enhanced stability.
+