-
-**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 FV3. 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.
-
diff --git a/docs/Chapter2.md b/docs/Chapter2.md
deleted file mode 100644
index 54310c2b7..000000000
--- a/docs/Chapter2.md
+++ /dev/null
@@ -1,52 +0,0 @@
-Cubed-sphere grid {#cube}
-=========================================
-
-##Chapter 2
-
-GFDL will provide the additional documentation by the end of May 2021.
-
-(This information is a reproduction of Appendix A from PL07)
-
-In Section 2 of PL07 the flux-form multidimensional transport scheme is discretized in general non-orthogonal curvilinear coordinates. The covariant and contra-variant wind vector components are presented in Eqs. (4) and (5) (of Section 2 of PL07) based on the local unit vectors \f$ (\vec{e_1},\vec{e_2}) \f$ of the coordinate system. Given the angle \f$(\alpha)\f$ between the two unit vectors
-
-\f[
- \cos \alpha = \vec{e_1} \cdot \vec{e_2} \\ \tag {2.1}
- \f]
-
-the covariant and contravariant components are related by the following relationships:
-
-\f[
- u = \tilde{u} + \tilde{v} \cos \alpha \\ \tag {2.2}
- \f]
-
-\f[
- v = \tilde{v} + \tilde{u} \cos \alpha \\ \tag {2.3}
- \f]
-
-or (solving for the contravaraint components)
-
-\f[
- \tilde{u} = \frac{1}{\sin^2 \alpha} [u - v \cos \alpha] \\ \tag {2.4}
- \f]
-
-\f[
- \tilde{v} = \frac{1}{\sin^2 \alpha} [v - u \cos \alpha] \\ \tag {2.5}
- \f]
-
-The winds on the cubed-sphere can be oriented to/from local coordinate orientation to a spherical latitude–longitude component form using the local unit vectors of the curvilinear coordinate system \f$(\vec{e_1},\vec{e_2})\f$ and the unit vector from the center of the sphere to the surface at the point of the vector location \f$(\vec{e_\lambda},\vec{e_\theta})\f$. Eqs.(2.6) and (2.7) represent the transformation from the spherical orientation \f$(u_{\lambda \theta},v_{\lambda \theta})\f$ to the local cubed-sphere form \f$(u,v)\f$ and the reverse transformation is presented in Eqs. (2.8) and (2.9).
-
-\f[
- u = (\vec{e_1} \cdot \vec{e_\lambda}) u_{\lambda \theta} + (\vec{e_1} \cdot \vec{e_\theta}) v_{\lambda \theta} \ \tag {2.6}
- \f]
-
-\f[
- v = (\vec{e_2} \cdot \vec{e_\lambda}) u_{\lambda \theta} + (\vec{e_2} \cdot \vec{e_\theta}) v_{\lambda \theta} \ \tag {2.7}
- \f]
-
-\f[
- u_{\lambda \theta} = \frac{ (\vec{e_2} \cdot \vec{e_\theta}) u - (\vec{e_1} \cdot \vec{e_\theta}) v }{ (\vec{e_1} \cdot \vec{e_\lambda}) (\vec{e_2} \cdot \vec{e_\theta}) - (\vec{e_2} \cdot \vec{e_\lambda}) (\vec{e_1} \cdot \vec{e_\theta}) } \\ \tag {2.8}
- \f]
-
-\f[
- v_{\lambda \theta} = \frac{ (\vec{e_2} \cdot \vec{e_\lambda}) u - (\vec{e_1} \cdot \vec{e_\lambda}) v }{ (\vec{e_1} \cdot \vec{e_\lambda}) (\vec{e_2} \cdot \vec{e_\theta}) - (\vec{e_2} \cdot \vec{e_\lambda}) (\vec{e_1} \cdot \vec{e_\theta}) } \\ \tag {2.9}
- \f]
diff --git a/docs/Chapter3.md b/docs/Chapter3.md
deleted file mode 100644
index 6f172c815..000000000
--- a/docs/Chapter3.md
+++ /dev/null
@@ -1,6 +0,0 @@
-Finite-volume formulation and flux evaluation {#flux}
-=========================================
-
-##Chapter 3
-
-GFDL will provide the additional documentation by the end of May 2021.
diff --git a/docs/Chapter4.md b/docs/Chapter4.md
deleted file mode 100644
index 3ce95f17d..000000000
--- a/docs/Chapter4.md
+++ /dev/null
@@ -1,75 +0,0 @@
-The vertical Lagrangian Solver {#lagrangian}
-=========================================
-
-##Chapter 4
-
-GFDL will provide the additional documentation by the end of May 2021.
-
-(This information is a reproduction of Sections 3 and 4 from HCZC20)
-
-###4.1 Lagrangian vertical coordinates
-
-A *Lagrangian* vertical coordinate is used in FV3. This coordinate uses the depth of each layer (in terms of mass or as geometric height) as a prognostic variable, allowing the layer interfaces to deform freely as the flow evolves. Further, the flow is constrained within the Lagrangian layers, with no flow across the layer interfaces (even for non-adiabatic flows). Instead, the flow deforms the layers themselves by advecting the layer thickness and by straining the layers by the vertical gradient of explicit vertical motion. This form is automatically consistent with the LR96 scheme, avoids the need for explicit calculation and dimensional splitting of verticaladvection, greatly reduces implicit vertical diffusion, and has no vertical Courant number restriction.
-
-FV3 uses a hybrid-pressure coordinate based on the hydrostatic surface pressure \f$p_s^*\f$:
-
-\f[
- p_k^* = a_k + b_kp_s^* \\ \tag {4.1}
- \f]
-
-where *k* is the vertical index of the layer interface, counting from the top down, and \f$a_k\f$, \f$b_k\f$ are pre-defined coefficients. Typically, the top interface is at a constant pressure \f$p_T\f$, so \f$a_0 = p_T\f$ and \f$b_0 = 0\f$. The spacing of the levels depends on the particular application, and is chosen depending on how high of a model top is desired, where additional vertical resolution is applied (typically in the boundary layer, but sometimes also near the tropical tropopause), and where to transition from hybrid \f$b_k > 0\f$ to pure pressure \f$b_k = 0\f$ coordinates.
-
-###4.2 Prognostic variables and governing equations
-
-The mass per unit area \f$\delta m\f$ can be expressed in terms of the difference in hydrostatic pressure \f$\delta p^*\f$ between the top and bottom of the layers; and, using the hydrostatic equation, can also be written in terms of the layer depth \f$\delta z^1\f$:
-
-\f[
- \delta m = \frac{\delta p^*}{g} = -p \delta z \\ \tag {4.2}
- \f]
-
-The continuous Lagrangian equations of motion, in a layer of finite depth \f$\delta z\f$ and mass \f$\delta p^*\f$, are then given as
-
-\f[
- D_L \delta p^* + \nabla \cdot (V \delta p^*) = 0 \\
- D_L \delta p^* \Theta_v + \nabla \cdot (V \delta p^* \Theta_v) = 0 \\
- D_L \delta p^*w + \nabla \cdot (V \delta p^*w) = -g \delta z \frac{\partial p'}{\partial z} \\
- D_L u = \Omega u - \frac{\partial}{\partial x} K - \frac{1}{p} \frac{\partial p}{\partial x} \biggr\rvert_{z} \\
- D_L v = - \Omega u - \frac{\partial}{\partial y} K - \frac{1}{p} \frac{\partial p}{\partial y} \biggr\rvert_{z} \\ \tag {4.3}
- \f]
-
-Note that these equations are exact: no discretization has been made yet, and the only change from the original differential form of Euler’s equations is to integrate over an arbitrary depth \f$\delta p^*\f$. The operator \f$D_L\f$ is the “vertically-Lagrangian” derivative, formally equal to \f$\frac{\partial \psi}{\partial t} + \frac{\partial}{\partial z}(w \psi)\f$ for an arbitrary scalar \f$\psi\f$. The flow is entirely along the Lagrangian surfaces, including the vertical motion (which deforms the surfaces as appropriate, an effect included in the semi-implicit solver).
-
-####Prognostic variables in FV3
-
-Variable | Description
-:--------------: | :----------
-\f$\delta p^*\f$ | Vertical difference in hydrostatic pressure, proportional to mass
-\f$u\f$ | D-grid face-mean horizontal x-direction wind
-\f$v\f$ | D-grid face-mean horizontal y-direction wind
-\f$\Theta_v\f$ | Cell-mean virtual potential temperature
-\f$w\f$ | Cell-mean vertical velocity
-\f$\delta z\f$ | Geometric layer height
-
-The vertical component of absolute vorticity is given as \f$\Omega\f$ and \f$p\f$ is the full nonhydrostatic pressure. The kinetic energy is given as \f$K = \frac{1}{2}(\tilde{u}u + \tilde{v}v)\f$: since FV3 does not assume that the horizontal coordinate system is orthogonal, we use the covariant (\f$u\f$ and \f$v\f$) components of the wind vector as prognostic variables and the contravariant (\f$\tilde{u}\f$ and \f$\tilde{v}\f$) components for advection, avoiding the need to explicitly include metric terms. See PL07 and HL13 for more information about covariant and contravariant components.
-
-The nonhydrostatic pressure gradient term in the \f$w\f$ equation is computed by the semi-implicit solver described section 5, which also computes the prognostic equation for \f$\delta z\f$. There is no projection of the vertical pressure gradient force into the horizontal; similarly, there is no projection of the horizontal winds \f$u\f$, \f$v\f$into the vertical, despite the slopes of the Lagrangian surfaces.
-
-Finally, the ideal gas law:
-
-\f[
- p = p^* + p' = \rho R_dT_v = - \frac{\partial p^*}{g \delta z} R_d T_v \\ \tag {4.4}
- \f]
-
-where
-
-\f[
- T_v = T(1 + \in q_v)(1 - q_{cond}) \\ \tag {4.5}
- \f]
-
-is the “condensate modified” virtual temperature, or density temperature, is used to close the system of equations. Here, \f$d_{cond}\f$ is the (moist) mixingratio of all of the liquid and solid-phase microphysical species, if present. When the gas law is used, the mass \f$p^*\f$ in this computation must be the mass of only the dry air and water vapor, and not including the mass of the condensate (non-gas) species. A rigorous derivation of the virtual and density temperatures is given in K. Emanuel, Atmospheric Convection (1994, Oxford), Sec. 4.3.
-
-These equations are also applicable to hydrostatic flow, in which \f$w\f$ is not prognosed and \f$p = p^*\f$ is entirely hydrostatic.
-
-
-______________________________
-1 In this document, to avoid confusion we write \f$\delta z\f$ as if it is a positive-definite quantity. In the solver itself, \f$\delta z\f$ is defined to be negative-definite, incorporating the negative sign into the definition of \f$\delta z\f$; this definition has the additional advantage of being consistent with how \f$\delta z\f$ is defined, being measured as the difference in hydrostatic pressure between the bottom and top of a layer.
diff --git a/docs/Chapter5.md b/docs/Chapter5.md
deleted file mode 100644
index 70ff389be..000000000
--- a/docs/Chapter5.md
+++ /dev/null
@@ -1,42 +0,0 @@
-The horizontal discretization and Lagrangian Dynamics {#horizontal}
-=========================================
-
-##Chapter 5
-
-GFDL will provide the additional documentation by the end of May 2021.
-
-(This information is a reproduction of Section 5 from HCZC20)
-
-###5.1 Dynamics along Lagrangian surfaces
-
-The layer-integrated equations (4.3) are discretized along the Lagrangian surfaces and integrated on the “acoustic” or “dynamical” time step \f$\delta t\f$ using forward-backward time-stepping as in LR97. The vertical velocity \f$w\f$ is a three-dimensional cell-mean value and partially advanced using the advection scheme. The geometric layer depth \f$\delta z\f$ is simply the difference of the heights of the successive layer interfaces, which with \f$\delta p^*\f$ defines the layer-mean density and the location of the Lagrangian surfaces. The air mass is the total air mass, including water vapor and condensate species.
-
-FV3 places the wind components using the D-grid (following Arakawa’s terminology), which defines the winds as face-tangential quantities. The D-grid permits us to compute the cell-mean absolute vorticity \f$\Omega\f$ exactly using Stokes’ theorem and a cell-mean value of the local Coriolis parameter \f$f\f$, without performing any averages or interpolations. The wind components themselves are face-mean values “along the cell edges” (not cell-mean values).
-
-Following the notation from L04, PL07, and HL13, we can write the discretized forms of (4.3), excluding the vertical components, as:
-
-\f[
- \delta p^{*(n + 1)} = \delta p^{*n} + F[\tilde{u^*}, \delta p_y] + G [\tilde{v^*}, \delta p_x] \\ \tag {5.1}
- \f]
-
-\f[
- \Theta^{n + 1} = \frac{1}{\delta p^{*(n+1)}} [\Theta^n \delta p^{*n} + F[X^*, \Theta_y] + G[Y^*, \Theta_x]] \\ \tag {5.2}
- \f]
-
-\f[
- w^* = \frac{1}{\delta p^{*(n+!)}} [w^n \delta p^{*n} + F[X^*, w_y] + G[Y^*, w_x]] \\ \tag {5.3}
- \f]
-
-\f[
- u^{n + 1} = u^n + \triangle \tau [Y(\tilde{v^*}, \Omega_x) - \delta_x (K^* -v \nabla^2 D) + \hat{P_x}] \\ \tag {5.4}
- \f]
-
-\f[
- v^{n + 1} = v^n + \triangle \tau [X(\tilde{u^*}, \Omega_y) - \delta_y (K^* -v \nabla^2 D) + \hat{P_y}] \\ \tag {5.5}
- \f]
-
-The quantities \f$\hat{P_x}\f$ \f$\hat{P_y}\f$ are the horizontal pressure-gradient force terms computed as in L97, the primary difference being that the forces due to hydrostatic and nonhydrostatic pressures arecomputed separately, and that the hydrostatic pressure gradient computation uses the log of the pressure to improve accuracy. The vertical nonhydrostatic pressure-gradient force is evaluated by the semi-implicit solver described in Chapter 5; only the forward advection of \f$w\f$ is performed during the Lagrangian dynamics, producing a partially-updated \f$w^*\f$.
-
-For stability, the pressure gradient force is evaluated backwards-in-time: the flux terms in the momentum, mass, and entropy equations are evaluated forward by the advection scheme, and the resulting updated fields are used to compute the pressure gradient force. This forward-backward time-stepping is stable without needing to use predictor-corrector or Runge-Kutta methods.
-
-In nonhydrostatic simulations it is recommended that the time off-centering for the horizontal pressure-gradient force be consistent with that used in the semi-implicit solver, which includes the vertical nonhydrostatic pressure-gradient force computation,to ensure consistency between the two. If the semi-implicit solver is run fully-implicit \f$(\alpha=1)\f$ then the pressure-gradient force should be evaluated fully backward \f$(\beta = 0)\f$; otherwise use \f$\beta = 1 - \alpha\f$
diff --git a/docs/Chapter6.md b/docs/Chapter6.md
deleted file mode 100644
index 2a6ada505..000000000
--- a/docs/Chapter6.md
+++ /dev/null
@@ -1,46 +0,0 @@
-The Nonhydrostatic Solver {#solver}
-=========================================
-
-##Chapter 6
-
-GFDL will provide the additional documentation by the end of May 2021.
-
-(This information is a reproduction of Section 6 from HCZC20)
-
-An equation for \f$z\f$ can be derived from the definition of \f$w\f$:
-
-\f[
- w = \frac{Dz}{Dt} = D_Lz + \vec{U} \cdot \nabla z \\ \tag {6.1}
- \f]
-
-The time-tendency of geopotential height is then the sum of the advective height flux along the Lagrangian interfaces and the vertical distortion of the surfaces by the gradient of \f$z\f$. Discretizing:
-
-\f[
- z^{n + 1} = z^n + F [\tilde{u^*}, \delta z_y] + G[\tilde{v^*}, \delta z_x] + w^{n + 1} \triangle t \\ \tag {6.2}
- \f]
-
-Since \f$w\f$ is solved for on the interfaces, we can then simply take the vertical difference to get \f$\delta z\f$.
-
-Recalling that the Lagrangian dynamics in (7) only performs the forward advection of the vertical velocity, yielding \f$w^*\f$, we then need to evaluate the vertical pressure-gradient force:
-
-\f[
- w^{n + 1} = w^* -g \delta z^{n + 1} \delta_z p'^{n + 1} \\ \tag {6.3}
- \f]
-
-The pressure perturbation \f$p'\f$ can be evaluated from the ideal gas law,
-
-\f[
- p' = p - p^* = \frac{\delta p}{g \delta z} R_d T_v - p^* \\ \tag {6.4}
- \f]
-
-requiring simultaneous solution of \f$w\f$, \f$p'\f$, and \f$\delta z\f$ using a tridiagonal solver.
-
-There is an option to off-center the semi-implicit solver toreduce implicit diffusion. The parameter \f$\alpha\f$ can be varied between 0.5 and 1 to control the amount of off-centering, with \f$\alpha = 1\f$ being fully-implicit. As discussed in Chapter 4 this off-centering parameter should be set to \f$\alpha = \beta - 1\f$, consistent with that used for the horizontal pressure-gradient force.
-
-The boundary conditions used are \f$p' = 0\f$ at the model top, and \f$\vec{U} \widehat{\cdot} n_s = 0\f$ at the lower boundary of \f$z = z_s\f$. This is the “free-slip” boundary condition, that the lower boundary is a streamline. The surface vertical velocity \f$w_s\f$ can be computed from (11) by advecting the surface height \f$z_s\f$:
-
-\f[
- w_s = \frac{z^*_s - z_s}{\triangle t} \\ \tag {6.5}
- \f]
-
-where \f$z_s^*\f$ is the advected value and \f$z_s\f$ the height of the topography.
diff --git a/docs/Chapter7.md b/docs/Chapter7.md
deleted file mode 100644
index 7261675be..000000000
--- a/docs/Chapter7.md
+++ /dev/null
@@ -1,112 +0,0 @@
-Stabilization and filtering options {#stabilization}
-=========================================
-
-##Chapter 7
-
-(This information is a reproduction of Chapter 2 from LPH17)
-
-###7.1 Divergence damping
-Horizontal divergence (along a Lagrangian surface) is computed as a cell integrated quantity on the dual grid:
-
-\f[
- D = \frac{1}{\Delta A_c} \left[ \delta_x (u_c \Delta y_c sin\alpha) + \delta_y (v_c \Delta x_c sin \alpha) \right] \\ \tag {7.1}
- \f]
-
-The Laplacian of D can also be computed as a cell-integrated quantity on the dual grid:
-
-\f[
- \nabla^2 D = \frac{1}{\Delta A_c} \left[ \delta_x \left(\frac{\delta_x D}{\Delta x} \Delta y_c sin\alpha \right) + \delta_y \left(\frac{\delta_y D}{\Delta y} \Delta x_c sin \alpha \right) \right] \\ \tag {7.2}
- \f]
-
-This operator can be applied on ∇2D instead of D to yield ∇4D. The damping is then applied when the forward time step is taken for the horizontal dynamics along vertically-Lagrangian surfaces:
-
-\f[
- u^{n+1} = u^n + . . . + \nu_D \frac{\delta_x \nabla^{2N} D}{\Delta x} \\ \tag {7.3}
- \f]
-
-\f[
- v^{n+1} = v^n + . . . + \nu_D \frac{\delta_y \nabla^{2N} D}{\Delta y} \\ \tag {7.4}
- \f]
-
-where N (equal to the namelist parameter `nord`) is 1 for fourth-order and 2 for sixth-order damping. The nondimensional damping coefficient is given
-
-\f[
- \nu_D = (d_4\Delta A_{min} )^{N+1} \\ \tag {7.5}
- \f]
-
-in which d4 is the parameter `d4_bg` in the namelist, and ΔAmin is the *global* minimum grid-cell area. It is recommended that this parameter be set to a value between 0.1 and 0.16, with instability likely for higher or lower values. Note that divergence damping is necessary as there is no implicit damping on the divergence in FV3. An optional second-order ∇2 damping, in addition the higher-order divergence damping, can be applied as well; in this case the added damping is of the form ν2D(δxD/ Δx), where ν2D = d2ΔAmin. Typically, the coefficient for d2 should be much smaller—by at least an order of magnitude—than the higher-order coefficient, if it is used at all, since the second-order damping is only weakly scale-selective and will significantly diffuse even resolved-scale features.
-
-The divergence damping can also be modified to add an approximate Smagorinsky-type damping, in which second-order divergence damping can be added to the flow dependent on the amount of stretching and dilation in thepflow. In this case, the d2 in the expression for ν2D is replaced by dSΔt(D2 + ζ2)1/2, where dS is the Smagorinsky coefficient (typically set to 0.2 if used) and ζ is the relative vorticity interpolated to cell corners so as to be co-located with the divergence. This form of the damping coefficient is more physical than the artificial damping discussed in the rest of this chapter, and will typically be very weak except in regions of very strong flow deformation.
-
-Divergence and flux damping (described in the next section) are both applied entirely within Lagrangian surfaces; there is no explicit diffusion applied across the surfaces. However, in regions of vertical motion the Lagrangian surfaces are distorted in the vertical as they follow the flow upward or downward. The amount of the distortion depends on the along-surface gradient of the vertical velocity; so where the distortion is largest is where there is the strongest horizontal shearing of the vertical velocity, which is also where ∇2n of the scalar fields should be the largest.
-
-
-###7.2 Hyperdiffusion (flux, or “vorticity”) damping
-
-Traditionally in FV3 computational noise is controlled through two means: the explicit divergence damping, and the implicit, nonlinear diffusion from the monotonicity constraint used in computing the fluxes. However, the implicit diffusion may be too heavily damping of marginally-resolved flow features, and for high-resolution applications it may be beneficial to use a nonmonotonic scheme and then add a user-controllable hyperdiffusion instead. This added hyperdiffusion need not be as strong as the divergence damping, since the non-monotonic advection schemes are still weakly diffusive (while no implicit diffusion is applied to the divergence), and often the hyperdiffusion coefficient df (`vtdm4`) should be much less than the divergence damping coefficient d4.
-
-In FV3 the hyperdiffusion is primarily meant to diffuse the kinetic energy of the rotational component of the flow, similarly to how the divergence damping dissipates the kinetic energy in the divergence component of the flow. (For this reason, the hyperdiffusion is sometimes called “vorticity” damping.) The diffusion is applied to the vorticity *flux*, allowing application of diffusion to only the rotational component fo the flow without needing to explicitly compute the rotational component.To maintain consistent advection of the other prognostic variables—w, θv, and δp*—the fluxes for these quantities are diffused as well, so that the potential vorticity and updraft helicity are still consistently advected as if they were scalars. (There is no way to add diffusion to the tracer mass fluxes, since higher-order diffusion cannot ensure monotonicity preservation without performing an additional flux limiting, adding even more implicit diffusion.)
-
-The hyperdiffusion is equal to the order of the divergence damping, unless eighth-order divergence damping is used, for which the hyperdiffusion remains sixth-order. The diffusion operator itself is second-order, but higher order diffusion is easily computed by repeatedly applying the diffusion operator, as is done for the divergence damping.
-
-Vertical vorticity is a cell-integrated quantity on the model grid. The vorticity fluxes vΩ/Δx and -uΩ/Δy are used to update the vector-invariant momentum equations. We can apply damping on the vorticity as well; to maintain consistent advection, the same damping is applied to the mass, heat, and vertical momentum fields, all of which are co-located with the vorticity. This additional damping is beneficial when using a non-monotonic advection scheme, which lacks the implicit diffusion of monotonic advection.
-
-Since the diffusion added to the vorticity fluxe is known explicitly, the loss of kinetic energy due to this explicit diffusion can be computed. The lost energy optionally can be added back as heat, after applying a horizontal smoother to the energy density field (so as not to restore the grid-scale noise which the damping was intended to remove). This can greatly improve the dynamical activity on marginally-resolved scales.
-
-###7.3 Energy-, momentum-, and mass-conserving 2Δz filter
-
-Local Richardson-number dynamic instabilities can create instabilities, especially in higher-top models, if the vertical turbulence scheme in the physical parameterizations is either disabled or insufficiently-strong to relieve these instabilities. These instabilities can grow large enough to crash the model. To relieve these instabilities, FV3 has the option to use a local (2Δz), vertical mixing to release the instability. This is similar to the Richardson number based subgrid-scale diffusion formulations of Lilly (1962, Tellus) and of Smagorinsky (1963), although their isotropic formulations have been simplified so as to only act on vertical gradients and perform diffusion in the vertical. This filter is completely local (2Δz), diagnosing and acting only on adjacent grid cells, and is typically applied only in the stratosphere and above to avoid interference with physical dynamic instabilities in the free troposphere or boundary layer which are more accurately-simulated by a planetary boundary layer scheme or the resolved dynamics. This filter is applied at the same frequency that the physical parameterizations are updated.
-
-We compute the local Richardson number on cell interfaces. Recall that k = 1 is the top layer of the domain and the index increases downward:
-
-\f[
- Ri_{k-\frac{1}{2}} = \frac{g\delta z\delta_z\theta_v}{(\theta^{k}_{v}+\theta^{k-1}_{v})((\delta_z u)^2+(\delta_z v)^2)} \\ \tag {7.6}
- \f]
-
-
-If Ri < 1, then mixing M is performed, scaled by Ri so that complete mixing occurs if R ≤ 0:
-
-\f[
- M = max(1, (1-Ri)^2) \frac{\delta p^{*k}\delta p^{*(k-1)}}{\delta p^{*k} +\delta p^{*(k-1)}} \\ \tag {7.7}
- \f]
-
-The mixing is applied to each variable, including the winds interpolated to the A-grid (necessary for application of the physical paramterization tendencies) on a timescale τ (namelist parameter `fv_sg_adj`) which should be larger than the physics timestep to avoid suppressing resolved convective motions. The mixing is applied to the momentum (δp* ua, δp* va), total energy, air mass, and all tracer masses, so that all of these quantities are conserved:
-
-\f[
- \frac{\partial \phi^k}{\partial t} = - \frac{M}{\delta p^{*k} } \left(\phi^k - \phi^{k-1} \right) \frac{1}{\tau} \\ \tag {7.8}
- \f]
-
-
-\f[
- \frac{\partial \phi^{k-1}}{\partial t} = + \frac{M}{\delta p^{*k} } \left(\phi^k - \phi^{k-1} \right) \frac{1}{\tau} \\ \tag {7.9}
- \f]
-
-where ϕ is a generic scalar. Note that since total energy and momentum are both conserved, lost kinetic energy automatically becomes heat.
-
-This mixing is most useful for removing instabilities caused by vertically propagating waves near the top of the domain. The namelist variable `n_sponge` controls the number of levels at the top of the domain to which the filter is applied.
-
-###7.4 Model-top sponge layer and energy-conserving Rayleigh damping
-
-Two forms of damping are applied at the top of the domain to absorb vertically propagating waves nearing the upper boundary. The first is a diffusive sponge layer, which applies second-order damping to the divergence and to the vertical-momentum flux, and optionally also to the vorticity and mass fluxes if the hyperdiffusion is enabled. (This differs from earlier versions of FV3, which instead of adding explicit damping applied first-order upwind advection in the sponge layer, the strength of which is flow-dependent and not user-controllable.) The damping is computed in the same way as described earlier, although typically a very strong value of d2 is used to ensure the vertically-propagating waves are sufficiently damped. The additional ∇2 sponge-layer damping is applied to the top two layers of the model, with a weaker damping also applied in the third layer if dk2 > 0.05. Since the model top is at a constant pressure, not constant height, it acts as a flexible lid, and therefore does not reflect as strongly as a rigid lid would.
-
-The second form of damping is a Rayleigh damping, which damps all three components of the winds to zero with a timescale which depends on the pressure. Given a minimum timescale τ0 and a cutoff pressure pc the damping timescale is:
-
-\f[
- \tau(p^*) = \tau_0 sin \left( \frac{\pi}{2} \frac{log(p_c / p^*)}{log(p_c / p_T} \right)^2 \\ \tag {7.10}
- \f]
-
-The strength of the applied damping is then determined by the magnitude of the cell-mean 3D vector wind U3D, including the vertical velocity, divided by a scaling velocity U0. The damping is only applied if the horizontal wind speed exceeds a latitude-dependent threshold (currently 25cosθ) or if the vertical velocity is larger than a given threshold. The damping is then applied, at the beginning of each large (physics) time step and before the Lagrangian dynamics is first called, by:
-
-\f[
- u \longleftarrow u(1 + \tau U_{3D} / U_0)^{-1} \\ \tag {7.11}
- \f]
-
-The dissipated kinetic energy can then be restored as heat:
-
-\f[
- T \longleftarrow T + \frac{1}{2} U_{3D} * (1 - (1 + \tau U_{3D} / U_0)^{-2}) / C_v \\ \tag {7.12}
- \f]
-
-
-
-
diff --git a/docs/Chapter8.md b/docs/Chapter8.md
deleted file mode 100644
index cac9512ae..000000000
--- a/docs/Chapter8.md
+++ /dev/null
@@ -1,65 +0,0 @@
-Physics-dynamics coupling {#physics}
-=========================================
-
-##Chapter 8
-
-(This information is a reproduction of Chapter 3 from LPH17)
-
-###8.1 Staggered wind interpolation
-
-The coupling to physical parameterizations and surface models is straight forward; the primary complication is interpolating between cell-centered, orthogonal winds used by most physics packages and the FV3 staggered non-orthogonal D-grid. The unstaggered orthogonal wind is defined by horizontal components in longitude and latitude; the staggered non-orthogonal D-grid wind is defined by the winds tangential to the grid-cell interfaces by the horizontal covariant components in the cubed-sphere curvilinear components. A two-stage horizontal re-mapping of the wind is used when coupuling the physics packages to the dynamics. The D-grid winds are first transformed to locally orthogonal and unstaggered wind components at cell centers, as input to the physics. After the physics returns its tendencies, the wind tendencies (du/dt, dv/dt) are then remapped (using high-order algorithm) back to the native grid for updating the prognostic winds. This procedure satisfies the “no data no harm” principle — that the interpolation/remapping procedure creates no tendency on its own if there are no external sources from physics or data assimilation.
-
-###8.2 Condensate loading and mass conservation
-
-The mass δm, and thereby also δp*, in FV3 is the total mass of both the dry air and of the water categories, including the vapor and condensate phases; the precise number N of water species is dependent upon the microphysics scheme used, or may be zero. This incorporates the effect of condensate loading into the air mass without a special parameterization for loading. The dry weight (per unit area) can be given as:
-
-\f[
- g \delta m_d = \delta p^* \left( 1 - \sum_{m=1}^N q_m \right) = \left( \delta p^* - \sum_{m=1}^N Q_m \right) \\ \tag {8.1}
- \f]
-
-where Qm = δp* qm is the tracer mass. Dry mass should be conserved by the physical parameterizations; here we will assume this to be the case, so δmd should be a constant in each grid cell during the physical tendency updates. The condition for dry mass conservation is then given by
-
-\f[
- \delta p^{*(n+1)} = \delta p^{*n} + \delta \tau \sum_{m=1}^N \frac{dQ_m}{dt} = \delta p^{*n} \Delta M \tag {8.2}
-\f]
-
-
-where ΔM = 1 + δτ ΣNm=1 dqi/dt. Physics packages usually return the rate of change in tracer mass dQm /dt, and so is independent of whether the solver uses total air mass or just dry air mass (as is done in many weather models). The tracer update is then done by:
-
-\f[
- Q^{n+1}_m = Q^n_m + \delta \tau \frac{dQ_m}{dt} \tag {8.3}
-\f]
-
-
-or, using (8.2)
-
-\f[
- q^{n+1}_m = \left( Q^n_m + \delta \tau \frac{dq_m}{dt} \delta p^{*n} \right) / \left( \delta p^{*(n+1)} \right) \\ \tag {8.4}
-\f]
-
-\f[
- = \left( Q^n_m + \delta \tau \frac{dq_m}{dt}\delta p^{*n} \right) / \left( \delta p^{*n} \Delta M \right) \\ \tag {8.5}
-\f]
-
-The full mass-conserving update algorithm is then:
-
-\f[
- q^{*}_m = q^n_m + \delta \tau \frac{dq_m}{dt} \tag {8.6}
-\f]
-
-\f[
- \Delta M = 1 + \delta \tau \sum_{m=1}^N \frac{dq_m}{dt} \tag {8.7}
-\f]
-
-\f[
- \delta p^{n+1} = \delta p^{n} \Delta M \tag {8.8}
-\f]
-
-\f[
- \delta q^{n+1}_m = q^{*}_m / \Delta M \tag {8.9}
-\f]
-
-Note that often the mass of non-water species, such as ozone or aerosols, are considered so small that they are not included in δm; however, since their mixing ratio is still the quotient of the tracer mass and the total air mass, if the effects of water species are included in the total air mass their mixing ratios must still be adjusted by (8.9).
-
-
-
diff --git a/docs/Chapter9.md b/docs/Chapter9.md
deleted file mode 100644
index 7551a60df..000000000
--- a/docs/Chapter9.md
+++ /dev/null
@@ -1,67 +0,0 @@
-Grid refinement techniques {#grid}
-=========================================
-
-##Chapter 9
-
-(This information is a reproduction of Chapter 4 from LPH17)
-
-There is a need for increasingly high-resolution numerical models for weather and climate simulation, but also an increasing need for coupling the newly resolved scales to the large and global-scale circulations, for which limited area models are only of limited use. However, uniformly-high resolution global models are not always practical on present-day computers. The solution to this problem is to locally refine a global grid, allowing for enhanced resolution over the area of interest while also representing the global grid. FV3 has two variable-resolution methods: a simple Schmidt transformation for grid stretching, and two-way regional-to-global nesting. These methods can be combined for maximum flexibility.
-
-FV3 can also be configured as a doubly-periodic solver, in which the cubed-sphere is replaced by a Cartesian-coordinate doubly-periodic horizontal grid; otherwise the solver is unchanged. This can be useful for idealized simulations at a variety of resolutions, including very high horizontal resolutions useful for studying explicit convection.
-
-###9.1 Grid stretching
-Here we follow the development of HLT16. A relatively simple variable resolution grid can be created by taking the original cubed-sphere grid and applying the transformation of F. Schmidt (Beitr. Atmos. Phys., 1977) to “pull” grid intersections towards a “target” point, corresponding to the center point of the high-resolution region. This is done in two steps: the grid is distorted towards the south pole to get the desired degree of refinement, and then the south pole is rotated to the target point using a solid-body rotation. Distorting to the south pole means that the longitudes of the points are not changed, only the latitudes, greatly simplifying the transformation.
-
-The transformation of the latitude θ to ϑ is given by:
-
-\f[
- sin\vartheta = \frac{D + sin\theta}{1 + Dsin\theta} \\ \tag {9.1}
- \f]
-
-where the distortion is a function of the stretching factor c, which can be any positive number:
-
-\f[
- D = \frac{1 - c^2}{1 + c^2} \\ \tag {9.2}
- \f]
-
-Using c = 1 causes no stretching. Note that other forms for the transformation could also be used without making any other changes to the solver.
-
-Although the grid has been deformed, the solver still uses the assumption that the grid cells are bounded by great-circle arcs, which are not strictly identical to a Schmidt transformation of the cubed-sphere arcs of the unstretched grid.
-
-###9.2 Grid nesting
-
-Using grid nesting can greatly increase the flexibility of grid refinement in the model, at the cost of greater complexity in the solver. The major strength of grid nesting is its ability to use independent configurations on each grid, including different time steps and physical parameterizations, most appropriate for that particular grid. The ability to use a longer time step on the coarse grid than on the nested grid can greatly improve the efficiency of a nested-grid model; and choosing parameterizations independently allows values appropriate for each resolution without needing compromise or “scale-aware” parameterizations.
-
-Here we follow the development of HL13, with additional updates necessary for the nonhydrostatic solver. Implementing two-way grid nesting involves two processes: interpolating the global grid variables to create boundary conditions for the nested-grid, and then updating the coarse-grid solution with nested-grid data in the region they overlap. The goal is to do so in as efficient of a manner consistent with the finite-volume methodology.
-
-A major feature of FV3’s nesting is to use concurrent nesting, in which the nested and coarse grids run simultaneously, akin to how coupled models run their atmosphere and ocean components at the same time on different sets of processors. This can greatly reduce the amount of load imbalance between the different processors.
-
-The entire nesting cycle is as follows, starting at the beginning of call to the solver:
-
-- For each `p_split` step:
- - Call solver
- - Fetch boundary condition data from coarse grid
- - In Lagrangian dynamics, update boundary conditions at each ¢t by extrapolating from two earlier coarse-grid states.
- - Perform tracer transport and vertical remapping
- - Perform two-way update
-- Call physics
-
-Note that we do not do a compile cycle every coarse-grid time step, unlike many regional nested-grid models. The cycling can be carried out multiple times per physics time step, if more frequent updates of the boundary conditions and of the two-way communication are considered necessary. There is also an option to perform the last two-way update after the physics, instead of before, which changes how the physical parameterizations interact with the nested-grid solution passed to the coarse grid. Performing the update before calling the physics has been found to yield better results in real-data forecasts.
-
-Currently, nested grids in FV3 are constrained to be a proper refinement of a subset of coarse-grid cells; that is, each coarse-grid cell in the nested grid region is subdivided into N x N nested-grid cells. This greatly simplifies the nested-grid boundary condition interpolation and the two-way updating. Nested grids are also static and constrained to lie within one coarse-grid face. However, the algorithm does not require an aligned, static grid in one cube face, and any of these conditions may be relaxed in the future.
-
-The nested-grid boundary conditions are implemented in a simple way. Coarse-grid data is interpolated from the coarse grid into the halo of the nested grid, thereby providing the nested-grid boundary conditions. Linear interpolation, although it is simple and and is not conserving, does have the advantage of not introducing new extrema in the interpolated field. The boundary conditions for staggered variables are interpolated directly from the staggered coarse grids. Boundary conditions are needed for each prognostic variable, including the tracers; also, boundary conditions are needed for the C-grid winds, available at each half-time step, and for the divergence when using fourth or higher-order divergence damping.
-
-Finally, boundary conditions for the layer-interface nonhydrostatic pressure anomalies are also needed to evaluate the pressure-gradient force. Instead of interpolating these interface values from the coarse grid, they are instead diagnosed and interpolated from the other boundary condition variables using the same methods as the semi-implicit solver.
-
-Most nested-grid models perform time-interpolation between two coarse grid states on each time step, but since the grids are integrated concurrently in FV3, interpolation is not possible. Instead, we can extrapolate between two earlier coarse-grid states. If interpolated coarse-grid variables are available at times t and t - Δ τ, where Δ τ = N Δt, then the extrapolation for a given variable φ at time t + nδt (n=1,...,N) is given by:
-
-\f[
- \phi^{t+n\delta t} = \left( 1 + \frac{n}{N} \right) \phi^t - \frac{n}{N} \phi^{t-\Delta \tau} \\ \tag {9.3}
- \f]
-
-The extrapolation is constrained for positive-definite scalars so that the value of the boundary condition at t + Δτ is non-negative, which is done by the substitution φt - Δτ → min(φt - Δτ , 2φt ).
-
-Two-way updates from the nested to the coarse grid are performed consistent with the finite-volume numerics. Scalars are updated to the coarse grid by performing an average of nested-grid cells, consistent with the values being cell-averages. The staggered horizontal winds are updated by averaging the winds on the faces of nested-grid cells abutting the coarse-grid cell being updated, so that the update preserves the average of the vorticity on the nested-grid cells. In FV3 only the three wind components and the temperature is updated to the coarse grid; the air and tracer masses are not updated, trivially conserving their masses on the nested grid, and reducing the amount of noise created through overspecification of the solution on the coarse grid. Since the air mass determines the vertical coordinate, which will differ between the two grids, the averaged nested-grid data is remapped onto the coarse-grid’s vertical coordinate.
-
-
diff --git a/docs/DoxygenLayout.xml b/docs/DoxygenLayout.xml
deleted file mode 100644
index 171a20ef8..000000000
--- a/docs/DoxygenLayout.xml
+++ /dev/null
@@ -1,194 +0,0 @@
-