\documentclass[gmd, manuscript]{copernicus} \usepackage{acro} \DeclareAcronym{GCM}{ short = GCM, long = General Circulation Model } \DeclareAcronym{AGCM}{ short = AGCM, long = Atmospheric Global Climate Model } \DeclareAcronym{AOGCM}{ short = AOGCM, long = Atmospheric Ocean Global Climate Model } \DeclareAcronym{IPCC}{ short = IPCC, long = Intergovernmental Panel on Climate Change } \DeclareAcronym{ESM}{ short = ESM, long = Earth System Models } \DeclareAcronym{GFDL}{ short = GFDL, long = Geophysical Fluid Dynamic Laboratory } \DeclareAcronym{PBL}{ short = PBL, long = Planetary Boundary Layer } \DeclareAcronym{KE}{ short = KE, long = kinetic energy } \DeclareAcronym{ISA}{ short = ISA, long = International Standard Atmosphere } \DeclareAcronym{3D}{ short = 3D, long = three-dimensional } \DeclareAcronym{2D}{ short = 2D, long = two-dimensional } \DeclareAcronym{DCMIP}{ short = DCMIP, long = Dynamical Core Intercomparison Project } \DeclareAcronym{HS}{ short = HS, long = Held \& Suarez } \def\myhalf{1/2} \begin{document} \title{WAVETRISK adaptive dynamical core coupled to dry physics} \Author[]{G}{Ching-Johnson} \Author[2]{NK-R}{Kevlahan} \Author[3]{T}{Dubos} \Author[4]{F}{Hourdin} \affil[1]{Department of Meteorology, University of Reading, Reading RG6 7BE, United Kingdom} \affil[2]{Department of Mathematics and Statistics, McMaster University, Hamilton ON L8S 4K1, Canada} \affil[3]{Laboratoire de M\'et\'erologie Dynamique, \'Ecole Polytechnique, 91128 Palaiseau, France} \affil[4]{Laboratoire de M\'et\'eorologie Dynamique, Institut Pierre Simon Laplace, Paris, France} \correspondence{Gabrielle Ching-Johnson (g.n.ching-johnson@pgr.reading.ac.uk)} \runningtitle{Adaptive dry physics climate model} \runningauthor{G Ching-Johnson, N Kevlahan, T Dubos, F Hourdin} \received{} \pubdiscuss{} \revised{} \accepted{} \published{} \firstpage{1} \maketitle \begin{abstract} General circulation models (GCMs) of the atmosphere couple a dynamical core for fluid dynamics that are primarily two-dimensional with a physics model for processes like radiation and convection that are primarily one-dimensional. In this paper we fully describe the dry physics model of~\cite{Frederic_these}, which includes a soil column, radiation, vertical diffusion and convective adjustment. The climatology of the dry physics coupled with the {\sc wavetrisk-atmosphere} dynamical core is described. We show that the dynamically adaptive version of {\sc wavetrisk-atmosphere} can be stably and accurately coupled to the dry physics without scale-aware parameterizations or modifications to the dynamics-based adaptivity criteria. The turbulence generated by {\sc wavetrisk}--dry physics GCM is characterized by longitude-latitude projections of vorticity, three-dimensional vorticity isosurfaces and vertical profiles of energy spectrum power laws, which vary between $-2.8$ in the boundary layer to $-5.2$ in the stratosphere for the rotational mode. \end{abstract} \introduction \label{intro} The coupling of atmospheric and oceanic \acp{GCM} constitutes the backbone of global climate models (for which the acronym GCM is also used). Pioneered in the 1950's \citep{Phillips1956}, \acp{GCM} are specific numerical tools based on separating and coupling two fundamental components: \begin{enumerate} \item A ``dynamical core" that computes the ``large scale" fluid motions using approximate versions of the compressible or incompressible Navier--Stokes equations and transport of a thermodynamic prognostic variable, such as potential temperature or buoyancy. \item A set of ``physics" models, called ``parameterizations'', that account for physics other than fluid dynamics (e.g., radiation, phase changes, subsurface thermal conduction and hydrology). The physics includes subgrid scale (SGS) models to approximate the effect of unresolved fluid motions on the resolved scales (e.g., turbulent diffusion). \end{enumerate} The parameterization of subgrid scale motions relies on Reynolds decomposition, which separates scales into a ``large scale", or average, resolved component one one side, and an unresolved turbulent component, or perturbation on the other. In the Reynolds decomposition framework, the large scale is defined formally as the ensemble average of random realizations of the full turbulent flow. In this decomposition, the statistics of the SGS variables are assumed to be horizontally homogeneous and ergodicity is invoked to assume that ensemble averages are equivalent to both horizontal space and time averages (either a running mean, defined in each point of the continuous world, or a grid box average). Similarly, the plane-parallel approximation is used for radiation. The radiation model assumes that all quantities are locally uniform horizontally, above a plane surface. For both SGS motions and radiation, the assumption of horizontal homogeneity results in physics models which are one-dimensional in the vertical. \acp{GCM}, based on the separation of three-dimensional (fluid) dynamics and one-dimensional physics parameterizations, have proved to be very useful for understanding the climate machine, and in particular the link between convective and cloud processes and climate. From the beginning they have also been used successfully for prediction. An early success, in the 1960s, was weather prediction over a time horizon of a few days. In the late 1970s climate models based on coupled atmospheric and oceanic \acp{GCM} successfully predicted the subsequent global warming. They played a great role in the raising concern about the risks of global warming. Today climate models are used systematically to model the climate at seasonal, decadal and centennial time scales. They are also used to evaluate the sensitivity of the climate to variations in input data, most importantly the level of greenhouse gas emissions. Although \acp{GCM} have reached a high degree of maturity and sophistication, many challenges remain. Lastly, a major bottleneck of climate models is the enormous computational power and associated computational time required to produce useful results. Any increase in complexity, through the inclusion of more processes and higher resolution, requires more computational power. With only limited computational power, modellers are required to make numerous trade-offs between CPU time, physical and numerical accuracy, and the number of physical processes represented. Choosing horizontal grid resolution $\Delta x$ is the fundamental trade-off when designing and running a climate model. More resolution provides higher accuracy, but greatly increases the computational cost. Since the time step $\Delta t\propto 1/\Delta x$, doubling the resolution from $\Delta x$ to $\Delta x/2$ means a simulation year takes 8~times longer to compute on a given CPU (assuming memory is not a limitation), and in any case 8~times more operations must be done~ This is the main motivation for dynamically adaptive grids, in which the resolution is refined locally depending on the time-dependent structure of the flow. Note that in practice only the horizontal resolution is adapted, even if Lagrangian vertical coordinates with adaptive remapping and dormant layers could in principle provide vertical resolution adaptivity. Often, high resolution is not required everywhere and so uniformly increasing resolution wastes computational resources: some regions are over-resolved at some times and others remain under-resolved. Dynamically adaptive climate simulations have been proposed as a way of using computational resources optimally and ensuring uniform accuracy~\citep[e.g.,][]{dynam_adapt_Jablonowski,Kevlahan/Dubos:2019}. A dynamically adaptive GCM locally refines and coarsens the grid as a function of time and location based on dynamical, physical or numerical accuracy criteria~\citep{dynam_adapt_Jablonowski}. Beyond capturing the evolution of the small scale dynamics, adaptivity aids in capturing the effects of the interactions between small-scale and large scale processes~\citep{dynam_adapt_Jablonowski}. In theory, adaptive grids could improve the realism and accuracy of model simulations, but there are known complications and challenges that come with adaptivity. Most notably, as mentioned above, SGS physics parameterizations are often designed for a particular horizontal resolution and may not perform properly when the grid is refined and coarsened during the simulation~\citep{Collins13_grid_adapt}. This is expected to be the case for moist deep convection, when the grid resolution reaches a few tens of kilometres, or for dry or cloudy boundary layer convection for resolutions of 1 or 2~km \citep{Honnert2016}. The question of whether physics processes should be ``scale-aware" in this ``grey zones" of convection is a challenge that many researchers are actively exploring~\citep[e.g.,][]{ scale_aware2,scale_aware1}. Scale-aware physics parameterizations adjust themselves based on the current local resolution of the grid. At its most basic, a scale-aware physics parameterization adjusts its internal model appropriately in response to changes in local grid resolution. A simple example is resolution dependent grid-scale viscosity, where the non-dimensional viscosity $\nu\Delta t/\Delta x^2$ is kept constant when the grid resolution $\Delta x$ changes. In the most extreme case where the resolution is refined enough to allow for the process to be fully resolved, the SGS parameterization is deactivated. Specifically, in the case of dynamic adaptivity, a refinement criteria needs to be met before local refinement (or coarsening) can occur. A central challenge of adaptivity is whether these criteria needs to be based on both the physics and the dynamics, or where adapting on the dynamics alone is sufficient. The answer likely depends on the type of physics being modeled. In current models the criteria are dynamics-based. For example, based on gradients or vorticity \citep[e.g.,][]{FlowCriteria}, or based on the numerical approximation~\citep[e.g.,][]{ErrorCriteria, Adaptive_Mesh_criteria,Kevlahan/Dubos:2019}. In order to make progress addressing these questions, researchers need well-understood physics models that can work in both the adaptive and non-adaptive models. One strategy to tackle a particular bottleneck and test new approaches such as adaptive grids, is indeed to simplify other components or aspects of the model in order to avoid the complexity of running and analyzing results of a full \acp{GCM}. For example, models may be simplified by reducing the number or complexity of the components included. \acp{AGCM} simulates the general circulation of the atmosphere, while an \ac{AOGCM} couples the ocean and atmosphere sub-models to simulate the general circulation. The atmosphere and ocean \acp{GCM} can also be more or less complex, and include more or fewer physical processes. When addressing bottlenecks that concern the atmospheric dynamics, the ``physics" can be simplified further. The most basic atmosphere \acp{GCM} are hydrostatic and neglect topography, moist physics and ocean--atmosphere coupling. The land surface model is reduced to a basic radiation boundary condition and turbulent boundary layer model. Such highly simplified models are used to better understand physics--dynamics coupling and the effect of resolution. A classical way to develop new dynamical cores, or models, is to replace the atmospheric physics by a relaxation toward a prescribed latitude-altitude temperature field and include Rayleigh friction drag in the lowest layers to replace turbulent energy dissipation in the planetary boundary layer~\citep{HeldSuarez}. In this paper, we use a simple and versatile dry physics model developed by \citet{Frederic_these} to explore various atmospheric general circulation regimes for various planetary conditions. This model includes sub-models for radiation, convection, turbulence, the atmospheric boundary layer and the soil. This simple physics (with a few parameters modified) was used to investigate some aspects of the terrestrial atmospheric circulation. These included Hadley circulation and its sensitivity to planetary rotation rate and radiative properties, as well as the atmospheric circulation of Mars, Venus and Titan. It was also used to investigate atmospheric ``superrotation'' \citep{Hourdin1995}, a phenomenon dominates the circulation of Venus and Titan. We chose this particular dry physics model to explore physics--dynamics coupling in a fairly realistic, yet simple, dry physics model with all assumptions and approximations fully described. The lack of moisture keeps the physics simple, but is still complex enough for us to investigate basic physics-dynamics coupling questions. By using dry physics, without phase change, the scale at which parameterizations of SGS motions start to be questionable is reduced to a few kilometres. A strength of this physics model is that it can be easily coupled with different dynamical cores and adapted to model different planets, including gas giants. The final reason for focusing on this particular model is to fully document, implement and test a package that can be coupled to both adaptive and non-adaptive dynamical cores. This SGS physics parameterization package allows investigation of dynamic adaptivity, scale-awareness and the scale-dependence of basic physics processes, and could help reduce the need of climate model tuning with increases in resolution. With these goals in mind, this paper presents the simple dry physics model in detail and evaluates its climatology in the non-adaptive and adaptive cases by coupling it to an existing dynamical core. This initial characterization of the simple physics model provides guidance on physics-dynamics coupling for adaptive GCMs, and introduces the simple physics package as a way of comparing different dynamical cores by coupling them to the same physics package. Specifically, this paper addresses the following key questions: % \begin{enumerate} \item What is the climatology of the {\sc wavetrisk}--simple physics climate model? How does it differ from the \cite{HeldSuarez} climate model? \item Is the dynamically adaptive \ac{GCM} numerically stable and and able to capture the emergence of small scale structure? Is adaptivity efficient for seasonal climate simulations? \item Is scale-aware dry physics and/or physics dependent adaptivity required for stable and accurate dynamically adaptive climate simulations? \item Does the adaptive simulation capture the full range of active turbulent length scale in the atmosphere, as well as the correct power law scaling? \item Can the {\sc wavetrisk}--simple physics climate model be restarted successfully from a checkpoint with significantly increased resolution? \end{enumerate} The paper is organized as follows. Section~\ref{SimplePhysics} describes in detail the simplified dry physics model, explaining its components and discussing its assumptions, limiting features and potential improvements. Section~\ref{Configuration} describes the coupling of the dry physics model with the adaptive dynamical core {\sc wavetrisk}~\citep{wavetrisk,Kevlahan/Lemarie:2022} in both the adaptive and non-adaptive cases. Section~\ref{Results} presents and interprets the results of the physics-dynamics coupling in the non-adaptive and adaptive cases, exploring the climatology and the effects of seasons, and Section~\ref{Conclusions} summarizes our conclusions and proposes future research directions. \section{Simple Dry Physics Model}\label{SimplePhysics} This section introduces the simple dry physics model. We begin by summarizing the model's assumptions and then describe in detail each component of the model (radiation, turbulent diffusion, convection, planetary boundary layer, soil column), with a particular focus on boundary conditions and time integration. Lastly we describe the limitations of the package and its potential improvements. As with most physics models, the simple dry physics is implemented as a time integration split step and works on a set of vertical columns (of varying cross-section in the adaptive case). The fact that the dynamics works on individual horizontal layers and the physics works on individual vertical columns is a fundamental approximation of GCMs, and is well-suited to horizontal grid adaptivity. \subsection{Physics model assumptions} To properly understand a physics model it is essential to explicitly state its assumptions and approximations. This not only allows for full transparency, but also avoids inconsistent assumptions when the physics and dynamics are coupled~\citep{Lauritzen/etal:2022}. The basic assumptions of this physics package are: dry thermodynamics, compressible flow (i.e.\ atmosphere, not ocean), hydrostatic pressure, and a floating (Lagrangian) vertical hybrid coordinate. {\bf No moisture and compressibility approximations.} A fundamental approximation of this package is the assumption that moisture (e.g., water in its three phases) is not taken into account. In particular, condensible substances and their effects on the atmospheric dynamics are absent. Such models are called dry physics models. Secondly, the model assumes compressible flow, which means it is suitable for the atmosphere, not the ocean. A fluid is compressible if changes in pressure causes significant changes in density. When modeling the atmosphere, the density of dry air changes significantly with changes in atmospheric pressure. Therefore, these two basic assumptions of no moisture and compressibility need to be ensured when coupling with a dynamical core. {\bf Hydrostatic approximation.} Another fundamental approximation in this simple physics model is the hydrostatic vertical pressure. The coupled dynamical model must therefore also assume hydrostatic balance. The hydrostatic approximation constrains the vertical coordinates in the model and simplifies the vertical equation of motion. An important implication of the hydrostatic approximation is that dry convection must be parameterized and cannot be resolved (even if the local grid resolution permits it). {\bf Floating vertical coordinate.} There are many possible spatial coordinate systems that can be used when modeling a planet. For example, the Cartesian coordinate system, with coordinates $(x, y , z)$, where $z$ is the vertical height above mean sea level. In the case of simplified \acp{AGCM}, where the dynamics and the physics are coupled, the dynamics works on two-dimensional layers and the physics works on one-dimensional vertical columns. While $z$ might be one's first choice for a vertical coordinate, other possible vertical coordinates include pressure, mass and sigma (a normalized pressure coordinate). Vertical fluxes of heat causes air parcels to expand and compress, and in turn change the density $\rho$ (and associated buoyancy), which is represented in the mass continuity equation by % \begin{equation} \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial z}(w\rho) = 0, \end{equation} % \noindent where $w$ is the vertical velocity. When altitude $z$ is the vertical coordinate there are vertical fluxes across the horizontal layers and the dynamics must include vertical velocity as a prognostic variable (in addition to horizontal velocities in each layer). The vertical velocity must be diagnosed from the density, which can be a tedious process. An easier option, implemented in this physics package, is to use a Lagrangian vertical coordinate, sometimes known as a floating vertical coordinate, so there is no mass flux across horizontal layers. In the case of a Lagrangian vertical coordinate $\zeta$, the $z$ coordinate has the form $z = z(\zeta,t)$ and there is no explicit vertical velocity (and no mass flux across vertical layer interfaces) since the horizontal layers move vertically with the local vertical velocity. Therefore, the incremental mass $dM$ is % \begin{equation} dM = \rho dz = m\, d\zeta = \rho\frac{\partial z}{\partial \zeta}\partial \zeta \end{equation} % and $dM/dt = 0$, where $m = \partial M/\partial \zeta = \rho\partial z/\partial \zeta$. % Note that since the thickness of the horizontal layers changes over time, the vertical coordinates must be remapped periodically onto the initial vertical coordinates, or optimized target vertical coordinates. \subsection{Model parameterizations} % The physics package includes parameterizations for radiation, vertical turbulent diffusion of velocity and potential temperature, a planetary boundary layer surface scheme, dry convective adjustment, soil column (for realistic surface flux boundary conditions), and diurnal radiation forcing. All models are one-dimensional (vertical) in space, and therefore are applied independently to each one-dimensional vertical column on the sphere as an implicit split step after the dynamics time step (i.e., the physics package takes as input a set of vectors). Note that none of the physics sub-models explicitly require the horizontal cross-section of the column (i.e, its horizontal resolution), although horizontal grid resolution could implicitly modify the physics models and physics--dynamics coupling. This is one of the questions we are interested in. The time scheme for the physics step uses Strang splitting: \begin{enumerate} \item Explicit forward Euler for soil, surface flux, radiation, and dry convective adjustment. \item Implicit backwards Euler for vertical turbulent diffusion of velocity and potential temperature. \end{enumerate} In the following, we give full details of each sub-model, including their default parameters. More details on the physics time scheme is given in Section~\ref{sec:physics_time}. \subsubsection{Radiation} % Radiation is a process that has a major impact on the atmosphere. In general, since radiation is not a dynamical process and due to its complexity, it's parameterization involves a hierarchy of simplifying assumptions (e.g., local thermodynamic equilibrium, only vertical radiative transfer, simplified or neglected scattering, discretized spectral bands, \ldots). The radiation model included in the physics package is based on the two-stream approximation. The two-stream approximation equations make the plane-parallel assumption, in which radiation is a function of height and the atmosphere is horizontally homogeneous and stratified in flat, horizontal layers (at least on the horizontal scale of the column). The resulting downward and upward radiation fluxes at a specific frequency $\nu$ and height $z$ are % \begin{align} F_{+}(\nu, z) &= F_+(\nu, 0) \tau(\nu,0, z) + \int_{0}^{z'} \pi B(\nu,T(z')) K(\nu,z')\tau(\nu,z', z) \,dz' \nonumber\\ &= F_+(\nu, 0) \tau(\nu,0, z) + \int_{0}^{z'} \pi B(\nu,T(z')) \frac{\partial\tau(\nu,z', z)}{\partial z'} \,dz' \\ F_{-}(\nu, z) &= F_-(\nu, \infty) \tau(\nu, z, \infty) + \int_{z'}^{\infty} \pi B(\nu,T(z')) K(\nu,z')\tau(\nu,z, z') \,dz' \nonumber\\ &= F_-(\nu, \infty) \tau(\nu,z, \infty) - \int_{z'}^{\infty} \pi B(\nu,T(z')) \frac{\partial\tau(\nu,z, z')}{\partial z'} \,dz', \end{align}\label{two-stream} % where the $+$ subscript indicates an upward flux and the $-$ subscript indicates a downward flux. $B(\nu,T(z')$ is Planck's function: % \begin{equation} B(\nu,T(z')) = \frac{2h\nu^3}{c^2(e^{h\nu /kT(z')}-1)}, \end{equation} % where $h$ is Planck's constant. Planck's function describes the radiation emitted by a black body in thermal equilibrium at a given temperature $T(z') $. $\tau_{\nu}(z, z')$ is the transmission function, which represents the proportion of the incident radiation, at frequency $\nu$, which is emitted from layer $z$ to layer $z'$. Therefore, it is the proportion of the flux that is not absorbed and in this model it is crudely represented by % \begin{equation} \tau(z, z') = e^{-|\psi(z) - \psi(z')|^\alpha}, \end{equation} % where $\psi$ is a monotonic function and $\alpha$ represents the absorption weight. This is known as a band averaged transmission function, as it approximates the transmission function over a given frequency band. The particular transmission function differs depending on the radiation band, thus the monotonic function in the physics package depends on the frequency band. In general, transmission functions are multiplicative. However, due to band averaging of the transmission function this property is lost. $K(\nu,z')$ represents the absorption, at layer $z'$, of radiation at frequency $\nu$, and is dependent on the infrared absorbers (i.e.\ gas) in the atmosphere, % \begin{equation} K(\nu,z') = \sum_{i}{} \rho_i \kappa_i(p, T, \nu), \end{equation} % where $\kappa(p, T, \nu)$ is the absorption coefficient of gas $i$ and $\rho$ is the density of the gas $i$. The simple physics model assumes that the absorption coefficients depend only on height. The simple physics model assumes constant absorption coefficients for the longwave (visible) and shortwave (visible) bands. The longwave (infrared) absorption coefficient is $\kappa_L = 0.08$\,\unit{m^{-1}} and the shortwave (visible) absorption coefficient is $\kappa_S =0.99$\,\unit{m^{-1}} . In general, the transmission function $\tau_{\nu}$ is a function of $K(\nu,z')$. Due to the band averaging of the transmission function, $K(\nu,z')$ is not explicitly defined. However the band averaged transmission function does describe strong or weak absorption of a frequency band through $\alpha$. A strong absorption produces a steep decay of transmission and strong absorption of radiation over shorter paths~\citep{pierrehumbert_2010}. This means that for longer paths the absorption in the first part of the path is stronger, while a absorption is weaker further along in the path. Therefore, strong absorption is usually used for frequencies that are absorbed first in a path. To solve for the net upward and downward irradiance, the fluxes need to be integrated over all frequencies. The physics packages uses the dual band two-stream approximations by crudely integrating the fluxes over just two spectral bands, longwave (L) and shortwave (S). Therefore, the net flux of radiation is % \begin{equation} F_{rad} = \int_{0}^{\infty} F_-(\nu,z) - F_+(\nu,z) \,d\nu = F_L(z) + F_S(z), \end{equation} % where $F_L$ and $F_S$ are respectively the net longwave and shortwave radiation fluxes. \subsubsection{Shortwave radiation} % The shortwave net flux is the net result of the upward and downward fluxes evaluated for shortwave radiation, % \begin{equation} F_S(z) = F_{S+}(z) - F_{S-}(z). \end{equation} % The simple physics case neglects atmospheric emission of shortwave radiation, leaving the single initial terms as the upward and downward flux, % \begin{align} F_{S-}(z) &= \tau_S(z,\infty) \mu S \label{downShort}\\ F_{S+}(z) &= \tau_S(0,z)F_{S+}(0) \label{upShort} \end{align} % where $S$ is the downward flux at the top of the atmosphere, which in the case of shortwave radiation is the solar flux and $\mu$ is the cosine of the zenith angle. The solar flux represents the amount of flux from the solar star impinging on the outer atmosphere. The upward flux at height $z$ (Equation \ref{upShort}) represents the proportion of the upward flux at the ground transmitted through the atmosphere, while the downward flux represents the proportion of the incident solar flux transmitted through the layer at $z$. The boundary upward flux of radiation, $F_{S+}(0)$, is the portion of downward shortwave radiation that is reflected by the ground, % \begin{equation} F_{S+}(0) = -(1-\alpha)F_{S-}(0), \end{equation} % where $\alpha$ is the albedo, which depends on the reflective properties of the ground ($\alpha\approx 1$ for fresh snow, $\alpha\ll 1$ for forests). The albedo $\alpha=0.112$ for our simulations, independent of position, although the physics package allows location-dependent albedo (e.g., different for land and ocean regions). Therefore there is no shortwave radiation emitted by the ground itself, however the ground does emit longwave (heat) radiation and is taken into consideration in the longwave fluxes described below. The transmission factor $\tau_S$ for the shortwave spectral band depends on height according to % \begin{align} \tau_S(z,\infty) &= e^{-c_S \frac{\psi(z)}{\mu}}, \nonumber\\ \tau_S(0,z) &= e^{-c_S\frac{\psi(0)-\psi(z)}{\mu_0}}, \nonumber\\ \psi(z) &= \frac{p(z)}{p_{\mathrm rad}} \label{eq:trans_factor}, \end{align} % where $p_{\mathrm rad}=1000$\,\unit{hPa} is a reference pressure height and $c_S = 0.005$ is the attenuation of shortwave radiation at pressure height $p_{\mathrm rad}$. $\mu$ is the cosine of the zenith angle and $\mu_0$ is an average of $\mu$, $\mu_0=3/5$. \subsubsection{Longwave radiation} % The longwave radiative fluxes in the simple physics package are represented using the gray gas approximation of the two-stream equations. As in the case of the shortwave net flux, the longwave net flux is determined by % \begin{equation} F_L(z) = F_{L+}(z) - F_{L-}(z), \end{equation} % where $F_{L+}(z)$ and $F_{L-}(z)$ represent the upward and downward longwave fluxes respectively. A further assumption in the gray gas model is the neglect of longwave solar emission meaning $F_{L-}(\infty) = 0$. Also, since the longwave radiative fluxes are the result of the two-stream approximations (\ref{two-stream}) integrated over infrared frequencies, the integration over infrared frequencies of the second term is well approximated by the Stefan--Boltzmann law~\citep{pierrehumbert_2010}. More specifically, % \begin{equation} \int \pi B \frac{\partial}{\partial z'}\tau(\nu,z,z') \,d\nu \approx \sigma T(z)^4 \frac{\partial}{\partial z'}\tau_{L}(z,z'), \end{equation} % where $\sigma=5.67\times 10^{-8}$~W m$^{-2}$~K$^{-4}$ is the Stefan--Boltzmann constant. Therefore, the downward and upward flux of longwave radiation can be approximated as % \begin{align} F_{L-}(z) &= \int F_{-} \,d\nu = -\int \sigma T(z')^4\frac{\partial \tau_L(z,z')}{\partial z'}dz'\\ F_{L+}(z) &= F_{L+}(z=0)\tau_L(0,z) + \int \sigma T(z')^4\frac{\partial \tau_L(z',z)}{\partial z'}dz'. \end{align} % The transmission factor $\tau_L$ for the shortwave spectral band depends on height according to (\ref{eq:trans_factor}), but with strong absorption coefficient $C_L= 2.53$. The upward surface boundary radiation flux is comprised of the heat flux emitted from the ground itself and the fraction of the downward atmospheric longwave radiation that is reflected by the ground (i.e., not absorbed by the ground). Therefore, % \begin{equation}\label{boundaryuprad} F_{L+}(0) = \epsilon \sigma T_s^4 + (1-\epsilon)F_{L-}(0), \end{equation} % where $\epsilon = 1 $ and $T_s$ are the soil emissivity and the surface temperature. The emissivity and the albedo, mentioned in the shortwave section are both characteristics of the ground, but they describe the effects for different wavelengths. The albedo $\alpha$ is the fraction of shortwave radiation reflected by a surface, while the emissivity is the effectiveness of a surface, the ground in this case, to emit longwave (heat) radiation. Both parameters are an input to the initialization of the physics package. Rearranging (\ref{boundaryuprad}) gives the net longwave flux at the boundary, % \begin{equation} F_{L}(0) = F_{L+}(0) - F_{L-}(0) = \epsilon (\sigma T_s^4 - F_{L-}(0)). \end{equation} The transmission factors for the longwave spectral band are represented as % \begin{align} \tau_L(z,z') &= e^{-c_{LW} \sqrt{|\psi(z')-\psi(z)|}}\\ \psi(z) &=\left(\frac{p(z)}{p_{\mathrm rad}}\right)^2, \end{align} % where $c_{LW} = 2.53$ is the attenuation of the longwave radiation at pressure height $p_{\mathrm rad}=1000$\,\unit{hPa}. (A weak absorption option is also available for longwave radiation with $\psi(z) = p(z)/p_{\mathrm rad}^2$ .) Upon discretization, the longwave equations are computed as, % \begin{align} F_{L-}(l) &= \sum_{k>l}{} (\tau_L(l,k_-) - \tau_L(l,k_+))\sigma T(k)^4\\ F_{L+}(l) &= F_{L+}(z=0)\tau_L(0,z) + \sum_{k Ri_c$ a minimum eddy viscosity is imposed~\citep{Frederic_these}, % \begin{equation} \label{emin_coef} K_z^u = l(z)\sqrt{e_{\mathrm min}}. \end{equation} % where $e_{\mathrm min}$ is interpreted as a lower bound on the TKE. The minimum TKE is set to the small value $e_{\mathrm min}=10^{-16}$\,\unit{m^2\,s^{-2}} to avoid spurious vertical diffusion in inactive regions. The maximum diffusivity $K^u_{\mathrm max}$ is set to 1~\,\unit{m^2\,s^{-1}}. The physics model for vertical diffusion is therefore simply given by Fick's second law, % \begin{equation} \frac{\partial \textbf{u}}{\partial t} = \frac{\partial}{\partial\zeta} \left( K_z^u \frac{\partial \textbf{u}}{\partial\zeta}\right), \label{turb_momentum} \end{equation} % where the eddy viscosity $K_z^u$ is computed using the mixing length model (\ref{eq:eddy_viscosity}) and $\zeta$ is a floating vertical coordinate. As mentioned above, the diffusion equation dissipates resolved kinetic energy $K$ provided by the dynamics step at the rate $\epsilon$ and diffuses this kinetic energy vertically. The SGS physics model for heat includes a radiative flux term and a source term $\epsilon$ due to heating (due to the dissipation of TKE and, indirectly, from the dissipation of $K$) in addition to the turbulent diffusion term, % \begin{equation} h_\theta\frac{\partial \theta}{\partial t} = \frac{\partial}{\partial\zeta} \left( h_\theta K_z^\theta \frac{\partial \theta}{\partial\zeta}\right) - \frac{\partial F_{rad}}{\partial \zeta} + \epsilon, \label{turb_heat} \end{equation} % where $\theta$ is the potential temperature, $h$ is the enthalpy, $h_\theta = \partial h/\partial \theta$. Recall that we assume that $K_z^\theta = K_z^u$. The equation for resolved kinetic energy $K = 1/2 \|{\textbf u}\|^2$ is derived from the prognostic equation for velocity~(\ref{turb_momentum}) by taking its inner product with $\textbf{u}$ and rearranging, % \begin{equation} \label{TKE} \frac{\partial K}{\partial t} + \frac{1}{\rho} \frac{\partial \textbf{F}_K}{\partial \zeta} = - \epsilon \end{equation} % where $\textbf{F}_K = \textbf{u} \cdot (\frac{1}{\rho} F_u \hat{{\textbf z}})$ is the vertical flux of $K$. This provides an expression for the TKE dissipated to heat, $\epsilon$, which is included in the turbulent diffusion equation for potential temperature (\ref{TKE}), % \begin{equation} \epsilon = - (\frac{1}{\rho} F_u \hat{{\textbf z}}) \cdot \frac{\partial \textbf{u}}{\partial z} = K_z^u \left|\left|\frac{\partial \textbf{u}}{\partial z}\right|\right|^2. \end{equation} % Since our turbulence model is only a first order closure approximation, there is no need for a model equation for $\epsilon$: it is diagnosed directly from the prognostic variables. It is assumed that the part of the resolved kinetic energy that is dissipated at rate $\epsilon$ is instantaneously converted to heat (rather than first being converted to TKE). The heat equation (\ref{turb_heat}) must therefore include a corresponding source term $\epsilon$. In contrast, in the second-order $k$-$\epsilon$ model, the eddy viscosity $K_z^u \propto k^2/\epsilon$, where $k$ is the (unresolved) small scale TKE. $K_z^u$ therefore depends on two unresolved quantities that must be modeled. They are computed by solving two coupled time-dependent model equations for their tendencies. With a block cell vertical discretization, as explained in section Section~\ref{discrete}, the tendency equation for momentum, discretized at cell layer center $k$ is % \begin{align} D_t\textbf{u}_k &= -\frac{1}{\rho_k} \frac{\textbf{F}_{k+1/2} - \textbf{F}_{k-1/2}}{\Delta z_k}, \\ \textbf{F}_{k+1/2} &= \left(K^u \frac{\partial u}{\partial z}\right)_{k+1/2} = K_{k+1/2}^u \frac{u_{k+1}-u_k}{\frac{1}{2}(\Delta z_{k+1} + \Delta z_k)}\\ \textbf{F}_{k-1/2} &= \left(K^u \frac{\partial u}{\partial z}\right)_{k-1/2} = K_{k-1/2}^u \frac{u_k-u_{k-1}}{\frac{1}{2}(\Delta z_k + \Delta z_{k-1})}, \end{align} % where $D_t\textbf{u}_k$ is the material derivative of momentum and $\textbf{F}_{k\pm 1/2}$ are the vertical turbulent momentum fluxes at the upper and lower interfaces $k\pm1/2$. The tendency equation for heat follows the same scheme as the momentum equation, % \begin{align} h_{\theta} D_t\theta_k &= -\frac{1}{\rho_k} \frac{\textbf{F}_{k+1/2} - \textbf{F}_{k-1/2}}{\Delta z_k} \\ F_{\theta ,k+1/2} &= \left(K^u \frac{\partial \theta}{\partial z}\right)_{k+1/2} = K^u_{k+1/2} \frac{\theta_{k+1}-\theta_k}{\frac{1}{2}(\Delta z_{k+1} + \Delta z_k)}\\ F_{\theta ,k-1/2} &= \left(K^u \frac{\partial \theta}{\partial z}\right)_{k-1/2} = K^u_{k-1/2} \frac{\theta_k-\theta_{k-1}}{\frac{1}{2}(\Delta z_k + \Delta z_{k-1})}, \end{align} % where the turbulent dissipation rate $\epsilon$ is neglected and $\textbf{F}$ is the total heat flux, incorporating both the turbulent and radiative heat flux. These equations are completed by Neumann (flux) boundary conditions at the bottom and top of each vertical column. The flux at the top of the atmosphere is assumed to be zero for both velocity and potential temperature, $F_{u,N} = F_{\theta,N} = 0$, while the fluxes at the ground, $F_{u,0}, F_{\theta,0}$ are provided by the surface flux parameterization scheme discussed in the following section. The radiative heat flux is computed by the radiative transfer scheme and is an input to the turbulence model. During a time step the inputs to the model are the surface temperature and surface heat flux at the surface boundary, while the outputs of the model are the tendencies for velocity and potential temperature due to the parameterized SGS processes. \subsubsection{Surface flux scheme} % The surface heat and momentum drag flux are computed using the inputs of $z, \theta$ and $u$, at the lowest (surface) layer center of the column ($k=1$). The parameterization for the fluxes are the bulk formulae of~\citet{Louis_1979}, which are based on the Monin--Obukhov similarity theory for boundary layer turbulence that describes vertical profiles of velocity and potential temperature as a function of height $z$ and the Obukhov length $L$, $z/L$. The flux model uses the following parameters, derived by \citet{Louis_1979}: % \begin{minipage}{0.35\linewidth} \begin{align*} a &= \frac{\kappa}{\log{(z/z_0)}}\\ b &= 4.7\\ C^u &= 7.4\\ C^\theta &= 5.3\\ Ri_B &= \frac{gz(\theta - \theta_0)}{\theta (u^2+e_{min})}\\ R &= 0.74, \end{align*} \end{minipage} \begin{minipage}{0.5\linewidth} \begin{align*} c^* &= 2bC^*a^2\left(\frac{z}{z_0}\right)^{1/2}\\ f^* &= 1 - \frac{2b Ri_B}{1 + c^*\sqrt{-Ri_B}} \tag{\textrm{if $Ri_B < 0$}}\\ f^* &= \frac{1}{(1+bRi_B)^2} \tag{\textrm{if $Ri_B > 0$}}\\ C_D^u &= a^2uf^u\\ C_D^\theta &= \frac{a^2u}{R}f^\theta \end{align*} \end{minipage}\\ % \noindent where $z_0=0.1$~m is roughness length, $\theta_0$ is the potential temperature at $z_0$, $Ri_B$ is the bulk Richardson's number, $\kappa=0.4$ is the Von~Karman constant. The exponent * replaces $u$ or $\theta$ and $C_D^\theta$ are drag coefficients for velocity and potential temperature. $a$ represents the drag coefficient in the neutral surface condition case, while the coefficient calculated for $Ri_B<0$ and $Ri_B>0$ models the drag coefficients in the unstable and stable cases respectively. $R$ is an empirical parameter that relates $C_D^\theta$ to $C_D^\theta$. The surface fluxes are therefore computed as surface drags using the drag coefficients computed above, \begin{align} \textbf{F}_u &= -\rho C_D^u \textbf{u},\\ F_\theta &= -\rho C_D^\theta (\theta - \theta_0). \end{align} % \subsubsection{Soil}\label{soilexplain} The inclusion of a multi-layer soil model is rare in simple physics models. However, our model includes a soil column which stores and releases heat over a wide range of time scales in order to provide more accurate and realistic conditions at the surface boundary. In this model the only prognostic variable of the soil is its depth-dependent temperature $T_{s}(z,t)$, with $z<0$. The heat flux in the soil is % \begin{equation} \label{heatflux} F_s = - k \frac{\partial T_s}{\partial t} \end{equation} % where $k$ is the thermal conductivity in W [K m]$^{-1}$. From Fick's second law the temperature tendency in the soil is given by the diffusion equation~\citep{Frederic_these}, % \begin{equation} \frac{\partial T_s}{\partial t} = -\frac{1}{I} \frac{\partial F_s}{\partial z} = \frac{k}{I} \frac{\partial^2 T_s}{\partial z^2}, \end{equation} % where $k$ is the thermal conductivity and the thermal inertia $I = 3000$\,\unit{J\,m^2\,K^{-1}\,s^{-1/2}}. The soil is discretized by $N_S=10$ layers, where the heat fluxes are computed at layer interfaces $l=-N_S,\ldots,0$, where $l=0$ is the surface. The interfaces are located at generalized depths % \begin{align} \zeta'_l = \frac{\alpha^{l} - 1}{\alpha - 1}\sqrt{T}, \end{align} % where the time scale $T = 2000$\,\unit{s}$/\pi$ and exponent $\alpha=2$ are model parameters. Depth is therefore measured as in terms of time, where deep layers correspond to increasingly slow heat diffusion time scales. Including more vertical layers add more slow time scales to heat diffusion in the soil column. Table~\ref{tab:layer_depth} shows the time scales corresponding to each soil layer for the $N_S=10$ soil layer case considered in this paper. % \begin{table}[h!] \centering \begin{tabular}{|c|r|r|} \hline \textbf{Layer} & \textbf{Period} & \textbf{Depth [m]} \\ \hline 1 & 5.6 hours & 0.11 \\ 2 & 2.1 days & 0.33 \\ 3 & 11 days & 0.76 \\ 4 & 52 days & 1.6 \\ 5 & 220 days & 3.4 \\ 6 & 2.5 years & 6.9 \\ 7 & 10 years & 14 \\ 8 & 41 years & 28 \\ 9 & 170 years & 56 \\ 10 & 660 years & 110 \\ \hline \end{tabular} \medskip \caption{The time periods of each soil layer and their equivalent depth. 10 soil layers is sufficient for climate simulations of up to about 1000 years. The calculated depth in the table uses the heat capacity $C_p=2.2 \times 10^6$\,\unit{J\,m^{-2}\,K^{-1}} and thermal conductivity $\lambda=4.09$\,\unit{J\,s^{-1}\,m^{-2}\,K^{-1}}.} \label{tab:layer_depth} \end{table} The first interface below the surface is ($\zeta'_1 = \sqrt{T}$) and it includes the damping depth of the diurnal wave. The damping depth is the depth of a wave of period $P$ whose amplitude is reduced by $e$~\citep{IntroMicrometeorology}. More generally, the depth of penetration of a wave of period $P$, in days, is $\sqrt{P}$~\citep{Frederic_these}. More layers included in the soil model allows for the release of energy at more (longer) time scales, for example a diurnal, monthly and annual time scales. Fluxes at interface $l$ are computed by: % \begin{align} F_{s,l} = I \frac{T_{s,l+1/2} - T_{s,l-1/2}}{\frac{1}{2}(\Delta \zeta'_{l+1/2} + \Delta \zeta'_{l-1/2})}, \end{align} % where $\Delta \zeta'_{l+1/2} = \zeta'_{l+1} - \zeta'_{l}$. Therefore, the soil temperature tendency of soil layer $l+1/2$ is defined by % \begin{align} \frac{\partial T_s}{\partial t} = -\frac{1}{I} \frac{F_{s,l} - F_{s,l-1}}{\zeta'_l - \zeta'_{l-1}}. \end{align} % If there are $L_s$ layers, the Neumann boundary conditions for the soil column are % \begin{align} F_{s,0} & = F_{r} + F_H, \tag{\textrm{Flux at ground boundary, $l=0$}}\\ F_{s,L_s} &= 0, \tag{\textrm{Flux of bottom soil interface} $l=L_s$} \end{align} % where $F_{r}$ is the surface radiative flux and $F_H$ is the surface heat flux. \subsubsection{Dry Convective Adjustment} % Dry convection must be parameterized because the dynamics makes the hydrostatic approximation. The vertical turbulence model only simulates turbulent diffusion, not convection. It could therefore produce an unstable vertical temperature stratification, characterized by $\partial \theta/\partial z < 0$~\citep{Frederic_these}. The dry convective adjustment scheme of the physics model mitigates any instability by instantaneously relaxing the temperature to an adiabatic profile. A layer is unstable if its potential temperature is greater than that of the layers above. The convection parameterization scheme checks each column to see if a layer is unstable. If a particular layer of the one-dimensional vertical column is unstable, local mixing is applied. This local mixing is modeled by a temperature relaxation to the mass weighted temperature average of the unstable layers. This scheme also takes into account the transport of momentum that occurs with convection by mixing the momentum in the unstable layers. For each unstable layer in a column only a proportion $\alpha$ of the cell is mixed, where % \begin{equation} \alpha = \frac{\int |\overline{\theta} - \theta| \,dp}{\int \theta \,dp}, \end{equation} % and $\alpha < 0$, which is always verified \citep{Frederic_these}. \subsection{Physics package usage, input and output}\label{Package_Inputs} The vertical physics model is intended to be coupled with a dynamical core. This coupling is usually done via a split step method, where the ``dynamics'' and ``physics'' components take alternate times steps. The prognostic variables of the physics package are temperature, zonal velocity and meridional velocity. When a physics step is taken, the package outputs the tendencies of the prognostic variables. Note that each of the the physics sub-models (radiation, turbulent diffusion, dry convective adjustment, diurnal cycle) can be turned on or off separately. One of the major objectives of the package is its versatility in modeling different planets. In order to accommodate this goal, the package initialization allows users to read in input parameters to tailor the simulation to the desired planet. Beyond these parameters, calls to initialization routines set up the logistics of the grid. % \subsection {Numerical integration and boundary coupling} % \subsubsection{Time integration\label{sec:physics_time}} % To solve the system of differential equations at the next time step a numerical integration method is employed. The explicit Euler method is a first-order method commonly used to numerically integrate a system of differential equations. Given an autonomous system of ODEs, % \begin{equation} \frac{d\mathbf{S}}{dt}=\mathbf{F}(\mathbf{S}(t)), \end{equation} % where $\mathbf{F}(\mathbf{S}(t))$ is called the ``tendency''. The stencil of the Explicit Euler method is given by the linear Taylor series approximation of $\mathbf{S}(t*) = \mathbf{S}(t+\tau) + O(\tau^2)$ % \begin{equation} \mathbf{S}(t*) = \mathbf{S}_t + \tau \mathbf{F}(\mathbf{S}_t) \end{equation} % where: % \begin{itemize} \item $t$ is the current time \item $\tau$ is the time step (assumed fixed) \item $t*$ = t + $\tau$ is the time at the next discrete time \item $\mathbf{S}_t$ is the numerical approximation of $\mathbf{S}(t)$, $\mathbf{S}_t\approx\mathbf{S}(t)$. \end{itemize} A variety of time integration methods have been developed from the Explicit Euler method by making different approximations for the tendency of $\mathbf{S}$, $\mathbf{F}(\mathbf{S}$), (e.g.\ Runge-Kutta schemes). These methods differ in their order of accuracy and numerical stability properties. In the Explicit Euler method the tendency is simply approximated as $\mathbf{F}(\mathbf{S}_t)$, where the values of $\mathbf{S}$ are at time t. The implicit Euler method improves numerical stability by using the tendency for \textbf{S} at the the next time step ($\mathbf{F}(\mathbf{S}_{t*})$). This requires more computation since a linear system must be solved at each time step, in general. However, Implicit Euler has a much wider stable domain compared to the Explicit Euler method, allowing for a larger time step and increasing its desirability. The radiation model uses Explicit Euler, while the soil and turbulence models use Implicit Euler. %The variation of Euler's method that is used in this physics package is semi-implicit Euler's method where the tendency is $\mathbf{F}(\mathbf{S}_{t}, \mathbf{S}_{t*})$. This is more complex as the tendency is a function of both $\mathbf{S}_t$ and $\mathbf{S}_{t*}$. \subsubsection{Boundary coupling} % Since the sub-components of the physics are air and soil, it is essential to specify how they are coupled in the time integration scheme. The air and soil columns are coupled via their boundaries, which is the horizontal interface between the soil and air (i.e., the ground), and the boundary conditions need to be taken into consideration and coupled when evaluating the tendencies of the two components. The simplest boundary coupling is explicit, and the boundary conditions are evaluated at time $t$ from the previous state. Therefore, the time integration equation includes the boundary conditions $\mathbf{q_t}$, % \begin{equation} \mathbf{S}_{t*} = \mathbf{S}_t + \tau \mathbf{F}(\mathbf{S}_t,\mathbf{q}_t). \end{equation} % Since the air and soil are strongly coupled, coupling the boundary conditions implicitly is desirable for improved stability. (Note that both the soil and air physics strongly affect temperature.) Therefore, the boundary conditions are evaluated at time $t*$. The tendency is approximated as % \begin{equation} \mathbf{F}(\mathbf{S}(t),\mathbf{q}(t)) \approx \mathbf{F}(\mathbf{S}_t,\mathbf{S}_{t*},\mathbf{q}_{t},\mathbf{q}_{t*}) \end{equation} % \begin{equation} \mathbf{S}_{t*} = \mathbf{S}_t + \tau \mathbf{F}(\mathbf{S}_t,\mathbf{q}_t) + \frac{\partial \mathbf{S}_{t*}}{\partial \mathbf{q}_{t*}} (\mathbf{q}_{t*} - \mathbf{q}_t), \end{equation} % where $\partial \mathbf{S}_{t*}/\partial \mathbf{q}_{t*}$ represents the dependence of $\mathbf{S}_{t*}$ on the boundary conditions at the new time $t^*$. \subsection{Spatial discretization}\label{discrete} % The physics package models the effects of the small scale unresolved physics on each 1D vertical column on the sphere, taking the current state of the prognostic dynamical variables as inputs. The vertical discretization is a block cell discretization, with all prognostic variables located at the center of each layer cell, while the fluxes are located at the interfaces between each layer. Each column has $N$ layers, which is an input parameter to the package from the dynamics. There are a total of $N + 1$ interfaces, where the surface (ground) interface has index 0 and top interface (top of atmosphere) has index $N$. The soil column is considered separately, with negative indices. With the chosen floating vertical coordinates, the height of each cell layer ($\Delta z_k$) changes over time. Technical details are provided in~\ref{sec:software}. \subsection{Limitations, Decisions and Possible Extensions} \label{limits_decisions} % The simple dry physics package makes approximations to ensure it is relatively simple, stable and easy to use and understand. However, these choices necessarily mean the physics is idealized and not sufficiently realistic for many applications. First, a dry climate is assumed by the package which neglects condensible substances (e.g., water or methane) on the atmosphere. The lack of moisture is a major limitation, especially for planets with condensible substances that have a major effect on the atmosphere, like Earth. However, assuming a dry climate removes the added complexity of the effects of condensible substances and makes the physics much simpler. Due to this limitation, we cannot expect accurate results when modeling Earth. For example, the temperature of the surface temperature should be much higher than that of Earth. This is due to the fact that the lapse rate of the simulation is constrained to be close to the dry adiabatic lapse rate, to keep a physically consistent dry model, while the lapse rate for Earth is closer to the moist lapse rate. Secondly, the surface flux parameterization in the physics package was created empirically using Earth observation data, which includes moisture. Therefore the scheme may not be suitable for modeling dry planets and might require tuning when modeling planets whose condensible substances are not water. Finally, the current radiation parameterization uses only two frequency bands: shortwave (visible) and longwave (infrared). In general, radiation models are expensive, and therefore having only two frequency bands keeps the cpu time per step minimal. However, with only two radiation bands, when modeling Earth, the temperature will not be accurate in the stratosphere as UV light has a major effect on the temperature. Nevertheless, unlike for moisture, it is relatively simple to include additional radiation bands. Adding a third absorption band for UV light would improve the realism of the vertical temperature profile in the stratosphere, producing the expected trend of an increase in temperature with height. However, this addition would increase the computational cost of the radiation step. \section{Configuration and Numerical Implementation\label{Configuration}} The simple dry physics package needs to be coupled with a dynamical core to simulate the atmospheric general circulation of a planet. In this section we describe the {\sc wavetrisk} dynamical core~\cite{wavetrisk} used to evaluate the characteristics of the physics package, the configuration of the physics-dynamics coupling and its numerical implementation. We also present the base simulation of the~\cite{HeldSuarez} dry physics model, which is to be used for future comparison of the simple dry physics package. \subsection{{\sc wavetrisk} dynamical core} Since we want to better understand the need for scale-aware physics, and how dynamic adaptivity combines with a physics package, we use the {\sc wavetrisk} dynamical core. {\sc wavetrisk} is a global three-dimensional dynamically adaptive hydrostatic dynamical core that can be run in both adaptive and non-adaptive modes (as well as in atmosphere or ocean configurations). This allows us to directly compare equivalent adaptive and non-adaptive runs. {\sc wavetrisk} utilizes a hexagonal-triangular C-grid with the TRiSK horizontal discretization scheme, which is well known for its mimetic properties and stability. Dynamic horizontal grid adaptivity is implemented using a multiscale wavelet approach with user-specified minimum (coarsest) and maximum (finest) resolutions. Nodes and edges are added or removed at each scale as necessary to achieve a user-specified relative error tolerance for each prognostic variable (i.e., pseudo-density $\mu = \rho_0 \Delta z$, mass-weighted potential temperature $\Phi = \mu\theta$, and edge velocities $u,v,w$). A node/edge is retained if it is retained in any vertical layer. The computational grid is adapted after each full time step as follows: \begin{enumerate} \item Compute wavelet coefficients at each horizontal position/vertical layer/horizontal resolution scale for all prognostic atmosphere variables. \item If a wavelet coefficient is above the dimensional tolerance threshold for any prognostic variable at any vertical layer, keep the associated grid point. If not, delete it. \item Add to the active grid the nearest neighbour grid points at the same scale and at the finer scale. This allows for a local scale-dependent CFL criterion $u\Delta t/\Delta x\leq 1$ and accounts for quadratic nonlinearities (which can generate scales twice as small in a single time step). \item Interpolate all prognostic variables (including soil temperature) onto the new adapted grid. \end{enumerate} This process ensures that the grid adaptation algorithm conserves mass and various mimetic properties. The resulting adaptive grid is a set of columns of varying cross-section. {\sc wavetrisk} may be run either adaptively with a range of horizontal resolutions, or non-adaptively with a fixed horizontal grid resolution. Full details on the {\sc wavetrisk} model are given in~\cite{Kevlahan/Dubos:2019}. \subsection{Model configurations} % Time discretization uses an explicit low storage fourth-order Runge--Kutta step for {\sc wavetrisk} and an implicit Euler step for the physics package. Both the dynamical core and the physics package have the same time step, with {\sc wavetrisk} determining the time step and supplying it to the physics. The adaptive simulation uses the time step determined by the finest resolution. We use horizontal bi-Laplacian hyperdiffusion for vorticity $\nabla\times\mathbf{u}$, velocity divergence $\nabla\cdot\mathbf{u}$ and scalars to reduce grid scale noise. The hyperviscosity $\nu_4$ is defined non-dimensionally as a proportion of the maximum stable value for the given time step $\Delta t$ and minimum average grid size $\Delta x_{\mathrm{min}}$ ($\Delta t$ is based on $\Delta x_{\mathrm{min}}$). The maximum stable viscosity is the analytical value for the hexagonal/pentagonal C-grid on the plane, modified with empirical correction factors for the irregular multiscale hexagonal C-grid on the sphere. (The empirical corrections are conservative by 15\% to 33\%, and are based on Gershgorin circle theorem estimates evaluated on the actual multiscale grids.) Note that the TRisK discretization~\citep{Ringler/etal:2010} is stable without horizontal diffusion on a fixed resolution grid, but there is some grid scale noise. {\sc wavetrisk} uses Lagrangian (floating) vertical coordinates with periodic conservative remapping of the prognostic variables to the initial hybrid pressure grid. Note that vertical remapping introduces some vertical diffusion. The data structure combines the atmosphere layers and soil layers into a single vertical column, but only the atmosphere layers require remapping. All simulations have 30 hybrid pressure layers, with the top of the atmosphere at 225~Pa. The $A,B$ hybrid pressure coefficients are given in~\cite{Ullrich2012_DCMIP_TestCases} Table~XVIII. Layer~6 is at 850~hPa and layer 12 is at 533~hPa. There is no orography. There are 10 soil layers for all simulations. The physics uses Earth astronomical values and boundary drag parameterizations, including seasons (see Table~\ref{tab:astronomical}). The initial conditions for prognostic variables are based on an isothermal temperature profile of $250$\,\unit{K} and zero velocity. While it would be more realistic to have an initial temperature that depends on latitude, a horizontally uniform isothermal profile is a simple initial condition that has been used in Dynamical Core Intercomparison Projects (DCMIP) (e.g. \citet{DCMIP_ic}) and the climate dynamics develops within a few days. The reference resolution used to characterize the climatology of the {\sc wavetrisk}--simple physics coupling is a relatively coarse non-adaptive $2^\circ$ ($\sim$~240\,\unit{km}) simulation. A resolution of $2^\circ$ is considered sufficient for coarse climate simulations. The adaptive simulation uses a minimum resolution of $2^\circ$ ($\sim$~240\,\unit{km}) and a maximum resolution of $ 0.5^\circ$ (60\,\unit{km}). This means that the maximum possible refinement is 4 times and the maximum grid compression is $4^2 = 16$ (i.e., about 16 times faster than a non-adaptive simulation with the maximum resolution). The simulations are initialized at the coarsest $2^\circ$ resolution and refined and coarsened after each time step with a relative error tolerance of 2.5\% for based on the $\ell_2$ norms for each vertical layers. All simulations are run for 5 (Earth) years. This period includes two years for the means of the prognostic variables in all layers to converge to within about 2\% for runs with seasons. The value of 2\% was chosen (rather than a smaller value) since the climate variables have significant intrinsic variance due to turbulence. The basic parameters and labels for all runs discussed below are given in Table~\ref{tab:runs}. % \begin{table} \begin{tabular}{lcccc}\hline Run label & Resolution [$^\circ$] & Resolution [km] & Physics Model & Seasons \\ \hline HS~2.0 & $2^\circ$ & 240 & Held and Suarez & no \\ UNIF~2.0~NS & $2^\circ$ & 240 & Simple physics & no \\ UNIF~2.0 & $2^\circ$ & 240 & Simple physics & yes \\ UNIF~1.0 & $1^\circ$ & 120 & Simple physics & yes \\ UNIF~0.5 & $0.5^\circ$ & 60 & Simple physics & yes \\ ADAPT~0.5 & $2.0^\circ$,\, $1^\circ$,\, $0.5^\circ$ & 240, 120, 60 & Simple physics & yes \\ ADAPT~0.25 & $2.0^\circ$,\, $1^\circ$,\, $0.5^\circ$, $0.25^\circ$ & 240, 120, 60, 30 & Simple physics & yes\\ \hline\medskip \end{tabular} \medskip \caption{\label{tab:runs} Parameters for the simulation runs discussed in the following sections. All runs use 30 vertical pressure layers. All simple physics runs also include 10~soil layers. All simulations, except ADAPT~0.25, are run for 5~years, with data saved every 5 days. ADAPT~0.25 is a special case that is restarted from UNIF~2.0 with 8 times higher resolution and then run for an additional 100~days. All adaptive simulation use a relative tolerance of $\varepsilon = 2.5$\%, except ADAPT~0.25 which uses 3\%. The relative tolerance $\varepsilon$ determines the local horizontal grid resolution at each time step.} \end{table} The climatology of the physics--dynamics coupled simulations is characterized by 5-year zonal first- and second-order statistics: mean temperature, mean zonal velocity, eddy kinetic energy, eddy heat flux temperature variance and eddy momentum flux. In addition, we include the vertical profile of temperature, averaged over the sphere and instantaneous horizontal slices of temperature and vorticity. The results for the $0.5^\circ$ adaptive simulation include a visualization of instantaneous three-dimensional vorticity isosurfaces, zonally averaged vertical velocity and energy spectrum results. \section{Results\label{Results}} This section presents and discusses results from the uniform and adaptive simulation runs summarized in Table~\ref{tab:runs}. We begin in \S\ref{HS-base} by comparing the climatology from the Held--Suarez physics model (HS~2.0) with the climatology of simple physics without seasons (UNIF~2.0~NS). This provides a baseline to appreciate the new features provided by the more complex simple physics model. \S~\ref{unif-climatol-with-seasons} presents low resolution climatology for the simple physics with seasons (UNIF~2.0). These results form the basis for comparing with the high resolution uniform (UNIF~0.5) and adaptive (ADAPT) results in \S~\ref{unif-climatol-with-seasons}. This section shows how the climatology and instantaneous prognostic fields change as the horizontal resolution is increased by a factor of 4; it also validates and characterizes the ability of the adaptive simulation to capture the climatology and energy spectra of the non-adaptive high resolution simulation. Finally, \S~\ref{adapt_performance} quantifies the computational performance of the adaptive runs ADAPT~0.5and ADAPT~0.25 compared to equivalent uniform runs. Key metrics include the grid compression ratio, stability and the ability to restart from a uniform low resolution run at much higher adaptive resolution. \subsection{Held--Suarez (HS~2.0) compared to low resolution simple physics without seasons (UNIF~2.0~NS) \label{HS-base}} % The Held--Suarez physics model \citep{HeldSuarez} is the simplest dry physics model that can be coupled to dynamical cores and a standard benchmark used for comparison of other simplified physics models. Therefore our initial comparison for the simple dry physics package will be to this model which includes only a relaxation of temperature to a radiative equilibrium value and a layer dependent linear damping of velocities in the surface boundary layer. HS~2.0 is a non-adaptive {\sc wavetrisk} simulation, using the initial conditions of \citealt{Jablonowski2006}. Figure~\ref{HS_climatology} presents the zonal statistics from HS~2.0. It is important to note that, although it is a standard benchmark, and therefore an interesting comparison, the Held and Suarez physics model is fundamentally different from the simple physics model, and so we should not expect qualitatively or quantitatively similar statistics. In the Held and Suarez physics model temperature is relaxed exponentially to a specified latitude and height-dependent profile, based on Earth climatology, and Rayleigh friction damps velocity near the surface (i.e., for pressures greater than 700\,\unit{hPa}) with a 1~day relaxation time. There are no other physical processes represented and the physics does not directly modify the dynamics, except by through the temperature relaxation and boundary friction. In addition, the simple physics model includes diurnal radiation forcing and there is no diurnal scale in the \ac{HS} temperature relaxation (the relaxation time scale is about 40~days in the stratosphere). It is helpful to begin by comparing the simple physics climatology without seasons to the most basic physics model of~\cite{HeldSuarez}. To do this, we ran a 5 year {\sc wavetrisk} non-adaptive simulation coupled to simple physics without seasons (i.e., obliquity of $0^\circ$). We refer to this run as UNIF~2.0~NS. (Note that \cite{Kevlahan/Dubos:2019} found that the Held and Suarez results are relatively insensitive to resolution.) The Held--Suarez and simple physics results shown in Figure~\ref{fig:J5_ZonalStats_Noseasons} are broadly similar, although their structures differ in position, shape and intensity. The biggest differences are in the eddy \ac{KE}, which has a different zonal jet shape, and in the zonal velocity which has its highest velocity near the top of the atmosphere, rather than at a similar pressure as the highest eddy \ac{KE}. Furthermore, while the eddy heat flux and temperature variance differs between the simple physics and the \ac{HS} simulations, the overall structure of the two statistics are similar, with high values near the mid-latitudes and at a normalized pressure of $\sim0.85$. \begin{figure*}[t] \includegraphics[width=\linewidth]{Figures/HSJ5Z30_climatology.png}\\ (a) HS~2.0 climatology. \\ \bigskip % \includegraphics[width=\linewidth]{Figures/SimpleJ5Z30_climatology_noseasons.png}\\ (b) UNIF~2.0~NS climatology (no seasons). % \caption{Comparison of the climatology of the zonal statistics for the low resolution HS~2.0 and UNIF~2.0 NS simulations. All statistics were averaged every 5 days over 5 years. Although the statistics are broadly similar, their features differ in position, shape and intensity.} \label{fig:J5_ZonalStats_Noseasons} \end{figure*} Figure~\ref{TempProfileClimatologyNonAdapt} presents the temperature profile for the non-adaptive $2^\circ$ low resolution simple physics simulation without seasons. compared to the \ac{ISA} model profile and the equilibrium state of the \ac{HS} model~\citep{HeldSuarez}, explained in the previous section. The \ac{ISA} is a standard temperature profile for Earth. The simple physics profile, plotted in red, is spatially averaged over the sphere at each layer and temporally averaged over five years. The temperature profile decreases monotonically with height, without a well-defined tropopause. In comparison to the standard atmosphere and \ac{HS}, the simple physics temperature profile has a significantly warmer surface and cooler upper stratosphere. A hotter surface is expected to the lack of moisture in the atmosphere, while the continued decrease in upper atmosphere temperature is due to the lack of a UV radiation absorption band in the stratosphere, as mentioned in section \ref{limits_decisions}. % \begin{figure*}[t] \centering \includegraphics[width=0.6\textwidth]{Figures/TempProfile_NoSeasons.pdf} \caption[Climatology Vertical Profile of Temperature Averaged over five years]{Temperature profile averaged over five years from the UNIF~2.0~NS run (simple physics with no seasons) and the HS~2.0 run (Held and Suarez physics). The figure also displays the International Standard Atmosphere temperature profile (ISA).} \label{TempProfileClimatologyNonAdapt} \end{figure*} \subsection{Non-adaptive low resolution $2^\circ$ simulation climatology with seasons (UNIF~2.0) \label{unif-climatol-with-seasons}} % Figure~\ref{fig:J5_ZonalStats_Seasons} shows the 5~year seasonal zonal first- and second-order climatology statistics for UNIF~2.0. Unsurprisingly, the statistics are markedly different with seasons compared to the equivalent results without seasons shown in Figure~\ref{fig:J5_ZonalStats_Noseasons}. As pressure decreases with height, the temperature decreases, but there is also minimal zonal temperature variance near the top of the atmosphere. The large temperature variance at the poles is simply due to the effect of seasons. Furthermore, eddy \ac{KE} has a zonal jet at the mid latitudes with a peak variance at $P/P_S\sim~0.25$. Overall, simple physics with seasons has qualitatively and quantitatively distinct climatology from Held--Suarez and Simple physics without seasons. Only the mean zonal temperature and high altitude peaks of eddy kinetic energy at $\pm 40$~N/S and $P/P_S = 0.25$ qualitatively similar. An interesting feature of the simple physics with seasons climatology is the peaks in eddy \ac{KE}, eddy momentum and heat flux at $15^\circ$~N/S near the surface. This feature does not appear in the simple physics simulation without seasons shown in Figure~\ref{fig:J5_ZonalStats_Noseasons}. This suggests that seasonal radiation forcing greatly changes the climatology, both (unsurprisingly) at the poles, but also in the surface boundary layer near the equator. The Held--Suarez model cannot capture this effect. % \begin{figure*}[t] \centering % \includegraphics[width=\textwidth]{Figures/SimpleJ5Z30_climatology_summer.png}\\ {(a) Summer climatology: UNIF~2.0.}\\ \bigskip % \includegraphics[width=\textwidth]{Figures/SimpleJ5Z30_climatology_fall.png}\\ {(b) Fall Climatology: UNIF~2.0.}\\ \bigskip % \caption{Summer and fall zonal climatology statistics of the low resolution UNIF~2.0 simulation. The effect of seasonality is very strong. Extreme values are larger, there is strong N/S asymmetry and more pronounced fine scale structure compared the equivalent non-seasonal simple physics climatology shown in Figure~\ref{fig:J5_ZonalStats_Noseasons}.} \label{fig:J5_ZonalStats_Seasons} \end{figure*} In the following section we present results for the high resolution non-adaptive (UNIF~0.5) and adaptive (ADAPT~0.5) runs. We assess the ability of the dynamically adaptive wavelet technique to capture important features of UNIF~0.5, and compare the climatology of the high resolution simulations with the low resolution UNIF~2.0 results discussed here. \subsection{Comparison of adaptive (ADAPT~0.5) and uniform (UNIF~0.5) high resolution simulations \label{high_res}} % Dynamically adaptive GCMs have the potential to significantly improve accuracy and computational efficiency by adjusting the local horizontal grid resolution at each time step to ensure a given tolerance. In this section we investigate the key questions from the Introduction~(Section~\ref{intro}) regarding the dynamically adaptive \ac{GCM}. We focus on the adaptive physics-dynamics coupling's stability, ability to capture small-scale structure and the complete energy spectra. Our main way of validating ADAPT~0.5 is to compare its results directly with an equivalent uniform high resolution simulation, UNIF~0.5. We also compare with the low resolution UNIF~2.0 results presented in the previous section. The simple physics package coupled to {\sc wavetrisk-atmosphere} is run adaptively with a coarsest average grid resolution of 2$^\circ$ and a finest average grid resolution of 0.5$^\circ$ (i.e., local horizontal grid refinement of up to 16 times). There are 30 vertical atmosphere layers and 10 soil layers. The relative error tolerance is $2.5\%$ for all prognostic variables with the $\ell_2$ tolerance normalization determined separately for each atmosphere layer. The grid is adapted explicitly only on the prognostic dynamical variables, not on the physics (e.g., tendencies of the vertical diffusion or radiation sub-models). Simulations are run for 5~years and climatology statistics are computed over the full five years. The adaptive results ADAPT~0.5 are compared with an equivalent non-adaptive run UNIF~0.5. \subsubsection{High resolution climatology statistics \label{ref:highres_climatology}} % The high resolution first- and second-order climatology statistics of ADAPT~0.5 and UNIF~0.5 for summer and fall are compared in Figures~\ref{fig:ZonalStats_Summer} and \ref{fig:ZonalStats_Fall}. The adaptive and non-adaptive results are qualitatively very similar, although the extrema of the second-order statistics are slightly smaller in ADAPT~0.5. However, a comparison with the low resolution results shown in Figures~\ref{fig:J5_ZonalStats_Seasons} shows that the high resolution climatology statistics are not qualitatively very different from the low resolution statistics. We will therefore examine the {\em differences\/} between the ADAPT~0.5 and UNIF~0.5 climatology statistics. We also discuss instantaneousness temperature and vorticity fields as well as the rotational and divergent components of the energy spectra to get a complete statistical, dynamical and spectral understanding of of the climate model. These more sensitive diagnostics provide different measures of the differences between the adaptive and uniform simulations, as well as how the low and high resolution results differ. % \begin{figure*}[t] \centering % \includegraphics[width=\textwidth]{Figures/SimpleJ5J7Z30_climatology_summer.png}\\ {(a) Summer climatology: ADAPT~0.5.}\\ \bigskip % \includegraphics[width=\textwidth]{Figures/SimpleJ7Z30_climatology_summer.png}\\ {(b) Summer climatology: UNIF~0.5.} % \caption{High resolution summer climatology of the adaptive and non-adaptive simulations. Statistics were averaged data saved every 5~days computed over the last 3~years of the 5~year run.} \label{fig:ZonalStats_Summer} \end{figure*} % \begin{figure*}[t] \centering \includegraphics[width=\textwidth]{Figures/SimpleJ5J7Z30_climatology_fall.png} {(a) Fall climatology: ADAPT~0.5.}\\ \bigskip % \includegraphics[width=\textwidth]{Figures/SimpleJ7Z30_climatology_fall.png}\\ {(b) Fall climatology: UNIF~0.5.} % \caption{High resolution fall climatology of the adaptive and non-adaptive simulations. Statistics were averaged data saved every 5~days computed over the last 3~years of the 5~year run.} \label{fig:ZonalStats_Fall} \end{figure*} Figure~\ref{fig:ZonalStats_diff_SummerFall} shows the differences in climatology statistics between ADAPT~0.5 and UNIF~0.5 computed over the final 3~years of the simulations. These differences are greatest at high latitudes, where the adaptive grid is coarsest, especially around the solstices, (see Figure~\ref{fig:Adaptive_grid}), and at high altitudes where the kinetic energy density $\frac{1}{2}\rho u^2$ is very low. The differences are small between -$60^\circ$S and $60^\circ$N. The exception is eddy heat flux, which has relatively large differences between -$60^\circ$S and $30^\circ$S in the summer. The relative $l_2$ norm differences between ADAPT~0.5 and UNIF~0.5 are summarized in Table~\ref{tab:adapt_errors}. The differences are in the range of 4--17\% (except for temperature, which is much smaller, 0.5\%). These quantitative errors are larger than the 2.5\% relative tolerance on the prognostic variables, which is perhaps unsurprising given the turbulent and unsteady nature of the climate dynamics. % \begin{figure*}[t] \centering \includegraphics[width=\textwidth]{Figures/Simple_climatology_diff_summer.png} {(a) Summer climatology: difference UNIF~0.5 - ADAPT~0.5.}\\ \bigskip % \includegraphics[width=\textwidth]{Figures/Simple_climatology_diff_fall.png} {(b) Fall climatology: difference UNIF~0.5 - ADAPT~0.5.} % \caption{Difference between high resolution simulations UNIF~0.5 - ADAPT~0.5 for the summer and fall zonal climatology statistics. Note difference in scales compared to the climatology shown in Figures~\ref{fig:ZonalStats_Summer} and \ref{fig:ZonalStats_Fall}. The largest differences are at high altitudes and latitudes.} \label{fig:ZonalStats_diff_SummerFall} \end{figure*} % \begin{table}[h] \begin{tabular}{llll} \hline Field & Relative $l_2$ difference (summer) & Relative $l_2$ difference (fall) & Relative $l_2$ difference (full year) \\ \hline Temperature & 0.49\% & 0.12\% & 0.54\% \\ Zonal velocity & 15\% & 12\% & 13\% \\ Eddy heat flux & 8.7\% & 13\% & 18\% \\ Eddy kinetic energy & 13\% & 16\% & 10\% \\ Eddy momentum flux & 17\% & 16\% & 9.8\% \\ Temperature variance & 11\% & 4.1\% & 3.7\% \\ \hline \end{tabular} \bigskip\\ \caption{\label{tab:adapt_errors} Relative $l_2$ zonal climatology differences between the UNIF~0.5 and ADAPT~0.5 runs. The corresponding figures for summer and fall are shown in Figure~\ref{fig:ZonalStats_diff_SummerFall}. } \end{table} In addition to the $l_2$ errors between UNIF~0.5 and ADAPT~0.5, it is helpful to directly compare the meridional profile of vertically integrated second order statistics. Figure~\ref{fig:climat_vert_int} shows vertically integrated eddy momentum flux for summer, fall and the entire year for UNIF~2.0, UNIF~0.5 and ADAPT~0.5. These results confirm that that ADAPT~0.5 is much closer to UNIF~0.5 than the low resolution UNIF~2.0 is. % \begin{figure} \centering \begin{tabular}{ccc} \includegraphics[width=0.3\linewidth]{Figures/Summer_EddyMomentumFlux.pdf}& \includegraphics[width=0.3\linewidth]{Figures/Fall_EddyMomentumFlux.pdf}& \includegraphics[width=0.3\linewidth]{Figures/All_EddyMomentumFlux.pdf}\\ {(a) Summer.} & {(b) Fall.} & {(c) Full year.} \end{tabular} \caption{Vertically integrated eddy momentum flux for the UNIF~2.0, UNIF~0.5 and ADAPT~0.5 runs for summer, fall and full year. } \label{fig:climat_vert_int} \end{figure} \subsubsection{Instantaneous temperature and vorticity dynamics\label{instant_field}} % Climatology statistics do not directly provide any information about the dynamics, or ``weather", of the coupled dynamics--physics climate model. It is revealing to explore the instantaneous temperature and vorticity fields that determine the climate. In particular, we are interested in how the (turbulent) vorticity field changes with height and resolution and how well the adaptive method reproduces its small scale intermittent structure. Figure~\ref{fig:vort_temp_layers} shows the instantaneous vorticity and temperature on 24 June of year 4 for layers 1, 9 and 27 (pressure levels 993~hPa, 713~hPa, 0.011~hPa respectively) for UNIF~0.5. The turbulence during the summer solstice is qualitatively different in the two hemispheres, with a narrow band of band of small scale turbulence centred around $30^\circ$~N and large scale less intense turbulence extending over the entire southern hemisphere. This figure shows how the flow structure changes with height, with intense filamentary vorticity in the lower part of the atmosphere and smoother jet-like structures in the upper atmosphere. The temperature fields illustrate how the temperature is advected and mixed by the vorticity (like a passive scalar), with strong meridional temperature gradients. % \begin{figure} \centering \begin{tabular}{cc} \includegraphics[width=0.45\linewidth]{Figures/Vorticity_J7Z27_T311.png} & \includegraphics[width=0.45\linewidth]{Figures/Temperature_J7Z01_T311.png} \\ \multicolumn{2}{c}{\raisebox{20pt}{(a) Layer 27 (0.011~hPa).}}\\ \includegraphics[width=0.45\linewidth]{Figures/Vorticity_J7Z09_T311.png} & \includegraphics[width=0.45\linewidth]{Figures/Temperature_J7Z09_T311.png} \\ \multicolumn{2}{c}{\raisebox{20pt}{(b) Layer 9 (713 hPa).}}\\ \includegraphics[width=0.45\linewidth]{Figures/Vorticity_J7Z01_T311.png} & \includegraphics[width=0.45\linewidth]{Figures/Temperature_J7Z27_T311.png}\\ \multicolumn{2}{c}{\raisebox{20pt}{(c) Layer 1 (993~hPa).}} \end{tabular} \caption{Vorticity and temperature near the summer solstice for UNIF~0.5 on 24 June of year 4 for three representative layers.} \label{fig:vort_temp_layers} \end{figure} To evaluate how well the adaptive simulation captures the weather dynamics, Figure~\ref{fig:Adapt_nonAdapt_Vorticity} at layer 7 (860~hPa) on 24 June of year~4, when the the weather is very different in the northern and southern hemispheres. ADAPT~0.5 is qualitatively similar to UNIF~0.5, capturing the qualitative features of the fine scale vorticity field. Note that we have set the relative tolerance to 2.5\%, and the tolerance could be decreased (or increased) to balance computational cost and fidelity of the weather as desired. It is also important to note that because the flow is turbulent we cannot expect the details of the flow to match (e.g., shape and position of particular vortex filaments, number of eddies etc.). The adaptive grid shown in Figure~\ref{fig:Adaptive_grid} (bottom) shows that this difference in turbulence activity between the northern and southern hemispheres is reflected in the adaptive dynamical grid which has essentially no refinement north of about $50^\circ$. Comparing the instantaneous vorticity and temperature fields with the climatology statistics shown in Figures~\ref{fig:ZonalStats_Summer} gives a better understanding of how the instantaneous dynamics determines the climate. % \begin{figure*}[t] \centering \includegraphics[width=0.45\textwidth]{Figures/Vorticity_J5Z07_T311.png}\\ {(a) UNIF~2.0 vorticity.}\\ \bigskip \begin{tabular}{cc} \includegraphics[width=0.45\textwidth]{Figures/Vorticity_J7Z07_T311.png}& \includegraphics[width=0.45\textwidth]{Figures/Vorticity_J5J7Z07_T311.png}\\ {(b) UNIF~0.5 vorticity.}& {(c) ADAPT~0.5 vorticity.} \end{tabular} \caption{Comparison of vorticity on 24 June of year 4 for layer 7 ($P = 860$~hPa). The adaptive simulation ADAPT~0.5 captures the qualitative small scale structure of the non-adaptive UNIF~0.5 vorticity field. Note that since the dynamics are chaotic and non-stationary the locations of vorticity structures are not the same.} \label{fig:Adapt_nonAdapt_Vorticity} \end{figure*} % Finally, Figure~\ref{fig:Vorticity3D_J5J7Z30} shows isosurfaces of vorticity for the adaptive simulation on March~22 of year~5. The vorticity is primarily columnar, although there is some evidence of the intense atmospheric boundary layer structure visible in the eddy kinetic energy, eddy heat flux and eddy momentum flux shown in Figures~\ref{fig:ZonalStats_Summer},\ref{fig:ZonalStats_Fall}. The primarily vertical vortices suggests that adapting only horizontally is efficient for these climate dynamics. % \begin{figure*}[t] \centering \includegraphics[width=0.8\textwidth]{Figures/Vorticity3D_J5J7Z30.png} \caption{Isosurfaces of vorticity on March 22 of year 5 ($\omega = \pm 10^{-4}$\,\unit{s^{-1}}) showing vertical structure of vorticity.} \label{fig:Vorticity3D_J5J7Z30} \end{figure*} The preceding statistical and instantaneousness results confirm that the dynamically adaptive \ac{GCM} is stable and able to capture the emergence of small scale turbulence and the transitions between different seasonal dynamics. Furthermore, the similarity of the high resolution adaptive climatology to the low resolution non-adaptive climatology suggests that adapting only on the prognostic variables is sufficient, at least for dry physics. The implicit effect of the physics on the dynamics is sufficient to resolve time- and scale-dependent dry physics phenomena. We now show that the adaptive wavelet method reproduces the energy spectrum of the uniform high resolution simulation at all active length scales. \subsubsection{Energy spectra} % A key feature of any \ac{GCM} is the energy budget, and in order to properly represent turbulent atmosphere dynamics we need to ensure that the energy spectrum of the adaptive simulation is similar to that of the finest non-adaptive simulation is important over all length scales. Because the wavelet multiscale approach adapts in both position and length scale, one of its strengths is its ability to capture the full turbulent energy spectrum with a much smaller number of computational elements than an equivalent spectral method~\citep{goldstein2005cvs_scales}. We now verify that the adaptive simulation accurately captures the complete turbulent energy spectrum of the non-adaptive simulation. Since the flow is compressible, the rotational and divergent parts of the flow are both relevant. Recall that turbulence is characterized by a power law energy spectrum, where the power law $p, E(k)\propto k^p$ varies depending on the type of turbulence. Figure \ref{fig:Spectra_scaling} shows the power laws for each vertical layer for the $0.5^\circ$ non-adaptive simulation. The slope $p$ is found from linear fits on a $\log-\log$ of energy spectrum for length scales between 600\,\unit{km} and 2000\,\unit{km}. The energy spectra were averaged over five years. The power laws for the $0.5^\circ$ resolution simulations are between $-1.8$ and $-3$ for the divergent mode and between $-2.8$ and $-5.2$. For comparison, incompressible homogeneous isotropic three-dimensional turbulence has a power law of $-5/3$, while two-dimensional compressible turbulence has various spectra between $-2$ and $-4$. Overall, the divergent component has a shallower slope, meaning that it is more turbulent. Furthermore, the flow is the most turbulent in the atmospheric boundary layer for $P(z) > 600$\,\unit{hPa}. Figure~\ref{fig:Spectra} compares the rotational and divergent components of the energy spectra of the non-adaptive (left) and adaptive (right) simulation. Results are shown for layers 1 (993\,\unit{hPa}), 9 (713\,\unit{hPa}) and 27 (0.0111\,\unit{hPa}). Figure~\ref{fig:vort_temp_layers} shows an example of the instantaneous vorticity and temperature structure for these layers. The energy spectra for the adaptive simulation are very similar to the non-adaptive spectra, but slightly smoother. Most importantly ADAPT~0.5 captures the full range of active length scales and the correct turbulent slopes but with about 3 times fewer computational modes (see \S\ref{adapt_performance} below). \begin{figure*}[t] \centering \includegraphics[width=0.5\textwidth]{Figures/EnergyPowerLaws_J7Z30.pdf} \caption{Energy spectra data averaged over year 5 for the non-adaptive $0.5^\circ$ resolution case: spectral slope $p$ for each layer obtained from linear fits for scales between 600\,\unit{km} and 2000\,\unit{km}. Note that the divergent spectra are much shallower than the rotational spectra and the flow is most turbulent (shallowest slope) at about 700\,\unit{hPa}.} \label{fig:Spectra_scaling} \end{figure*} \begin{figure*}[t] \centering \begin{tabular}{ccc} \multicolumn{3}{c}{Rotational (divergence-free) component of energy spectrum.}\\ \includegraphics[width=0.3\textwidth]{Figures/Rot_Layer_1.pdf}& \includegraphics[width=0.3\textwidth]{Figures/Rot_Layer_9.pdf}& \includegraphics[width=0.3\textwidth]{Figures/Rot_Layer_27.pdf}\\ (a) Layer 1 (993~hPa). & (b) Layer 7 (713 hPa). & (c) Layer 27 (0.011 hPa). \bigskip\\ \multicolumn{3}{c}{Divergent (curl-free) component of energy spectrum.}\\ \includegraphics[width=0.3\textwidth]{Figures/Div_Layer_1.pdf}& \includegraphics[width=0.3\textwidth]{Figures/Div_Layer_9.pdf}& \includegraphics[width=0.3\textwidth]{Figures/Div_Layer_27.pdf}\\ (d) Layer 1 (993~hPa). & (e) Layer 7 (713 hPa). & (f) Layer 27 (0.011 hPa). \end{tabular} \caption{Energy spectra averaged over year 5 for the three representative layers shown in Figure~\ref{fig:vort_temp_layers}. We compare the low resolution uniform $2^\circ$ simulation (UNIF~2.0) , the high resolution uniform $0.5^\circ$ simulation (UNIF~0.5) and the adaptive simulation (ADAPT~0.5). The adaptive simulation results have been interpolated to the finest resolution ($0.5^\circ$) for the spectral analysis. The top row shows the rotational (divergence-free) component and the bottom row shows the divergent (curl-free) component of the energy spectrum. The adaptive results are consistent with the uniform high resolution results, apart from differences at the smallest scales. In particular, the adaptive simulation matches the spectral slope of the high resolution uniform simulation and resolves the full range of active scales.} \label{fig:Spectra} \end{figure*} \subsection{Computational performance of dynamic grid adaptivity \label{adapt_performance}} % We begin by summarizing the quantitative and qualitative features of the dynamic grid adaptivity. The adaptive grid has an average compression ratio of $2.9$, meaning that ADAPT~0.5 uses on average 2.9~times fewer active cells than UNIF~0.5. Over a simulation year the compression ratio varies depending on the season, with the highest compression at the solstices and the lowest compression at the equinoxes, due to the changing weather patterns of the season. The maximum compression ratios of about 3.7 occur at the solstices and the minimum compression ratios of about 2.4 occur at the equinoxes. Figure~\ref{fig:Adaptive_grid} shows the adaptive grids with the lowest and highest compression ratios during the spring equinox on March~22 of simulation year 5 and during the summer solstice on June 25 of simulation year 4 respectively. At the summer solstice the northern (hotter) hemisphere is less active than the southern (cooler) hemisphere. Overall, the atmosphere is much less turbulent at the solstices than at the equinoxes. The annual pattern of grid compression is clear in Figure~\ref{fig:compression}. These results confirm the ability of the adaptive grid to refine and coarsen, both locally in space and following the seasons, as required to capture the dynamically important features of the flow. \begin{figure*}[t] \centering \begin{tabular}{cc} \includegraphics[width=0.45\textwidth]{Figures/AdaptiveGrid_J5J7_T365.png}& \includegraphics[width=0.45\textwidth]{Figures/AdaptiveGrid_J5J7_T311.png}\\ (a) Spring equinox. & (b) Summer ssolstice. \end{tabular} \caption{Adaptive dual triangular grids from ADAPT~0.5 simulation with tolerance $2.5\%$. (a) On 5 March 22 of year 5, during an equinox. (b) On June 25 of year 4, during a summer solstice. At the summer solstice the northern (hotter) hemisphere is less active than the southern (cooler) hemisphere, leading to a higher compression ratio.} \label{fig:Adaptive_grid} \end{figure*} \begin{figure*}[t] \centering \includegraphics[width=0.6\textwidth]{Figures/compression.pdf} \caption{Compression ratio for the ADAPT~0.5 run with relative tolerance 2.5\% showing seasonal variation over 5 years. The solstices and equinoxes are labelled. The highest compression ratios are at the solstices and the lowest compression ratios are at the equinoxes.} \label{fig:compression} \end{figure*} Finally, we investigate whether a low resolution simulation can be stably restarted from a saved checkpoint at much higher dynamically adaptive resolution. This allows investigation of climate dynamics at higher resolutions for shorter periods of time and permits quickly spinning up the model at low resolutions. We restarted from a non-adaptive UNIF~2.0 checkpoint at 5~years (March 22), with maximum resolution increased by 8 times to $0.25^\circ$. This corresponds to a 64 times increase in horizontal resolution. We ran the simulation for an additional one year. This run is referred to as ADAPT~0.25. Figure~\ref{fig:J5J8_compression}~(top) shows the grid compression ratio as a function of time and Figure~\ref{fig:J5J8_compression}~(middle) shows the adapted horizontal grid 95 days after restarting. The relative tolerance is the same as for ADAPT~0.5, $\varepsilon = 2.5\%$. The restarted high resolution simulation ran stably with an average grid compression ratio of 5.3. After start-up, the grid compression varied between 3.3 and 9.5, depending on the season. Adding one more resolution level achieves an average compression ratios almost twice as large as ADAPT~0.5, with maximum compression ratios 2.6 times higher. The corresponding vorticity field shown in Figure~\ref{fig:J5J8_compression}~(bottom) illustrates how the adaptive grid tracks the intense vorticity filaments. This result shows the potential of adaptivity to both achieve relatively high grid compression (and corresponding speed-up) and the ability to investigate higher horizontal resolutions for shorter periods. % \begin{figure*}[t] \centering \includegraphics[width=0.6\textwidth]{Figures/J5J8Z30_compression.pdf}\\ (a) Compression ratio for ADAPT~0.25.\bigskip\\ \begin{tabular}{cc} \includegraphics[width=0.45\textwidth]{Figures/AdaptiveGrid_J5J8Z07.png}& \includegraphics[width=0.5\textwidth]{Figures/Vorticity_J5J8Z07.png}\\ (b) Adaptive grid for year 6 March 22. & (c) Layer 7 (860~hPa) vorticity field for year 6 March 22. \end{tabular} \caption{Results for ADAPT~0.25 simulation restarted after 5 years from UNIF~2.0 with relative tolerance $\varepsilon = 2.5\%$. The average compression ratio is 5.3.} \label{fig:J5J8_compression} \end{figure*} \conclusions\label{Conclusions} This paper presents the results of coupling a simple dry physics sub-grid scale model to the {\sc wavetrisk-atmosphere} dynamical core. We have four main goals: \begin{enumerate} \item Fully document the approximations and numerical implementation of the simple dry physics model introduced by \cite{Frederic_these}. \item Couple the simple physics package to the {\sc wavetrisk-atmosphere} multiscale adaptive dynamical core. To do this, the simple physics code has been rewritten in FORTRAN~2008 and simplified. \item Describe the 5-year climatology of the {\sc wavetrisk}--simple physics climate model, run non-adaptively at coarse resolution $2^\circ$ with 30 vertical atmosphere layers and 10 soil layers. The results without and with seasons are compared to an equivalent~\cite{HeldSuarez} climate model. \item Explore an adaptive {\sc wavetrisk}--simple physics climate model coupling and assess the need for scale aware parameterizations or grid adaptation criteria. The coarsest resolution is $2^\circ$ and we allow two layers of refinement: $1^\circ$ and $0.5^\circ$, for a maximum grid compression ratio of 16 times. The adaptive results are compared with an equivalent non-adaptive $0.5^\circ$ resolution simulation. \end{enumerate} % The basic approximations of the physical model are: dry atmosphere (i.e., no moisture), hydrostatic balance, compressibility and a floating vertical coordinate. The physics model is coupled to the dynamics as an implicit split step in time. It includes sub-grid scale parameterizations for radiation, vertical turbulent diffusion of velocity and potential temperature, a planetary boundary layer, surface fluxes, dry convection, diurnal cycle, seasons and, significantly, a soil column. The lack of moisture reduces model realism, but maintains simplicity and simplifies the interpretation of results. Model parameters, such as gravity and astronomical parameters, are set to Earth values. Compared to the Held \& Suarez and the International Standard Atmosphere, the temperature profile for the simple physics--{\sc wavetrisk} climate model decreases is significantly warmer near the surface of the atmosphere and decreases monotonically with height. These differences are likely due to two limitations of the simple dry physics package, specifically the lack of moisture and the simplicity of the radiation parameterization (e.g., no ultraviolet absorption band in the stratosphere). We presented zonal statistics of temperature, zonal velocity, eddy kinetic energy, eddy heat flux, temperature variance and eddy momentum flux. Without seasons, these statistics are qualitatively similar to those of the Held \& Suarez model, although the location of zonal jets is shifted and the eddy kinetic energy is about half as strong. This suggests the additional slow time scales introduced by the soil column do not significantly affect the climatology. Unsurprisingly, adding seasons dramatically modifies all statistical quantities. The adaptive simple physics--{\sc wavetrisk} climate model with horizontal resolutions $2^\circ, 1^\circ, 0.5^\circ$ is stable and accurate. With relative tolerance 2.5\% produces grid compression ratios of between about 1.8 times (at the equinoxes) and 3.0 times (at the solstices). This demonstrates the efficiency advantages of horizontal grid adaptivity, as it is able to track the active and inactive parts of the globe throughout the year. The zonal statistics are very similar to those of the $2\circ$ non-adaptive simulation and also agree well with the non-adaptive $0.5^\circ$ results. In addition to comparing the 5-year statistics, we also compared longitude-latitude projections of vorticity for the adaptive and non-adaptivity simulations on 25 June of year 4 when their is a strong north-south asymmetry in the turbulence. The weather is characterized by a narrow band of small scale turbulence centered around $30^\circ$\,N and a large scale less intense turbulence extending over the entire southern hemisphere. The adaptive simulation successfully captures the large and small scale vortical structures in each hemisphere. Finally, we presented isosurfaces of vorticity for March 22 of year 5. The vorticity is primarily columnar, but also includes intense atmospheric boundary layer structures, which are also visible in the statistics. To complement the physical space comparisons, we compared energy spectra for the rotational and divergent components for adaptive and non-adaptive simulations. These results confirm that the adaptive simulation accurately captures the full range of active turbulent scales, although there is some smoothing due to interpolation of the adaptive grid onto the full grid (needed to compute the spherical harmonics transform). Vertical profiles of the energy spectrum power law show strong dependence on depth, with the most turbulent region in the boundary layer, $P(z)>600$\,\unit{hPa}. Overall, we conclude that the simple dry physics model can be stably coupled to an adaptive dynamical core. Scale-aware parameterizations are not necessary, at least for dry physics. The climatology presented here could provide the basis for new dynamical core test cases since the simple dry physics can be easily coupled to different dynamical cores. The simple physics includes relatively sophisticated radiation, turbulence, convection and diurnal/seasonal forcing and thus fills a gap between the Held \& Suarez model and realistic moist physics models. Future work could explore an intercomparison of a set of dynamical cores coupled to the same simple dry physics package. Although scale-aware physics parameterization is not necessary for this dry physics model, it is quite likely to be necessary for moist physics. In future work, we will investigate coupling the adaptive {\sc wavetrisk} model to a basic moist physics, such as~\cite{kessler1969water}. This model has already been coupled to the Held \& Suarez model~\cite{ullrich2016moistheld}. We hope that the results presented here will stimulate interest in~\cite{Frederic_these}'s interesting dry physics model and has demonstrated the potential of adaptive dynamical cores to improve the efficiency and accuracy of climate modeling. Our results may also help climate modelers simplify the tuning required when grid resolution is increased by focusing their attention on sub-models not included in the simple dry physics. \codedataavailability{The {\sc wavetrisk} code, including the simple dry physics and Held--Suarez physics, are available at \url{https://github.com/kevlahan/wavetrisk_hydrostatic}. The test case {\tt ~/test/climate} is used to run the simulations. The data presented in this paper are available at this globus link: \href{https://globus.alliancecan.ca/file-manager?origin_id=2cd3329f-064f-4489-b634-eeeafa59f1a7\&origin_path=\%2F}{Globus} .} \appendix \section{Software Engineering Design Decisions and Challenges\label{sec:software}} %% Appendix A To couple the simple dry physics package to the dynamical core {\sc wavetrisk} certain design decisions were made in order to facilitate the physics-dynamics coupling and overcome challenges encountered due to dynamic adaptivity. This appendix section will summarize the general steps needed to couple the simple dry physics package to a dynamical core and the decisions made to couple the package to {\sc wavetrisk}. In order to utilize the physics package as a test case and couple it with a dynamical core an interface is used. This interface is the connector and facilitator of the entire model and will vary among climate models due to different dynamical core architectures and their methods of incorporating test cases. In terms of coupling the physics to {\sc wavetrisk}, we made the design decision to create an interface formatted as one of {\sc wavetrisk}'s test cases. The interface calls all physics initialization routines through a wrapper function which includes setting up the plugins of the physics and reading all constants and constraints set by the dynamics. As the facilitator, the interface couples the dynamics and physics by calling their associated time step routines, usually in a split step. When calling the physics time step a few aspects need to be taken into account, the precision of the both the dynamics and physics, the format of the main prognostic variables stored and the data structures of the prognostic variables, which is highly dependent on the grid. The physics package is written in single precision, while {\sc wavetrisk} is written in double precision. To avoid converting to and from both data types, before and after the physics step, we decided to compile the physics package with the Fortran flag, {\tt -freal-4-real-8}. This flag converts all real numbers (single precision) in the physics package to real (8) values (double precision), allowing for compatibility with {\sc wavetrisk}. The key input variables required for a physics package time step includes temperature, zonal velocity and meridional velocity at each vertical layer, along with the pressure and geopotential of the vertical columns. Calculating the pressure and geopotential requires the simple integration from the surface upwards for each column and the utilization of the hydrostatic approximation. In comparison, the temperature and velocities require a more in-depth conversion. {\sc wavetrisk} saves the mass weighted potential temperature as its prognostic variable. Therefore, when calling the physics, the interface converts the value, at each vertical layer of a column, to potential temperature and then to the required temperature input, in Kelvin, of the physics. Upon return from the physics call, the interface converts the temperature tendencies back to mass weighted potential temperature using the reverse process. Regarding the zonal and meridional velocities, wavetrisk does not store the velocities at the center of the node, instead it stores three velocity edge values. To retrieve the zonal and meridional velocities at the node, interpolation using the three edge values and neighbouring columns' edge values is required. The inverse procedure is also applied when the physics returns the zonal and meridional tendencies for each column. In general, the data structures storing the prognostic variables of both the dynamics and physics are different. {\sc wavetrisk} utilizes a hybrid data structure which is optimal for the grid chosen and for adaptivity and load balancing associated with the use of parallelism (see~\citep{wavetrisk_datastruct} for more on the {\sc wavetrisk} grid architecture and data structure). Conversely, the physics package utilizes a regular data structure to store the values for each column. The regular data structure is a 2D array, with each row representing a different vertical column and each column of the array represents a vertical layer. Upon initialization of the physics package, the number of vertical columns to be sent at each time step is fixed. At each physics step, the package takes in the fixed number of vertical columns, the time step, the simulation day, the fraction of the day and the columns' prognostic variable values, pressure and geopotential at each vertical layer. A major overhead of the simulation, even in the case of parallelism, is the retrieval and conversion of all vertical columns from the hybrid data structure to the required regular structure for the input of the prognostic variables, as well as the conversion of the output tendencies from the regular to hybrid data structure. The use of a 2D data structure constrains the package as it expects the user to call the physics for all vertical columns at once. However, another major obstacle this structure produces is in the case of parallelism and adaptivity. {\sc wavetrisk} utilizes {\tt mpi} and domain decomposition of the columns of the sphere. When adaptivity is used, a major concern is ensuring correct load balancing on each CPU. Therefore, {\sc wavetrisk} allows for the re-balancing of the domains on each CPU, to ensure approximately equal computational loads. The complication arising with the ability of rebalancing is that the physics package saves the soil and surface temperatures internally for each column. To reduce the overhead and overcome the rebalancing obstacle, we have modified the physics step to call the physics for vertical columns individually and the step is completed when the last column's tendency returns from the call. However, to allow for single column transmission, the physics package and {\sc wavetrisk} required some fine tuning and additions. The major addition to the physics package is the creation of the single column module. It allows the model to call the physics package for each column individually by providing the soil and surface temperatures to the dynamics. This module contains a wrapper subroutine that the interface calls, including all required inputs for the original physics call plus the added soil and surface temperatures. Once the wrapper subroutine is called, it updates the soil and surface temperatures for the column and calls the original physics routine. The subroutine outputs the original tendencies and the new updated soil and surface temperatures for the interface to save. The module also contains a subroutine that allows for the change in latitude and longitude required when sending a different column at each physics call. A major change to the physics package is the input of the number of desired soil layers when initializing the parameters of the grid in the package. Originally this was fixed to 10, but a variable number of layers is desired to test the importance of the soil model. The only major change to {\sc wavetrisk} was in the grid initialization. To send single columns to the physics with the possibility of rebalancing, {\sc wavetrisk} must save the soil temperatures and surface temperature of each column. The vertical layers were therefore extended, below the surface, employing negative and zero indexes for each vertical column to represent the soil layers ($k<0$) and surface ($k=0$) temperature respectively. The interface facilitates this by setting the number of soil layers, before all grids are initialized. The only drawback of extending the grid is that all prognostic variables' data structures needs to be extended, even though only the mass-weighted potential temperature negative indexes are used to save the soil and surface values. Therefore, there is additional unnecessary memory overhead. In future, it is recommended that only the required variable's data structure be extended, but this will require a more extensive update, not only to the architecture of the grid, but also to the architecture of the calls for the dynamics step. Another {\sc wavetrisk} design decision was the format and storage of the extra temperature values. It was decided that the soil and surface temperatures would be stored in the mass-weighted potential temperature's data structure in the negative and 0 indexes respectively. More importantly, the temperatures would not be converted to mass-weighted potential temperature, but would be retained as temperature (units of Kelvin). This decision was made to avoid the extra overhead that comes with the conversion, as the values are not used by the dynamics. Finally, while the main decision to send individual grid columns to the physics required extra temperature values to be stored by the dynamics, no other physics variables need to be stored. In many models, as in the {\sc wavetrisk-ocean} model configuration, a new dynamics variable tke is required to keep track of the turbulent kinetic and the tke evolves over time with the physics' k-epsilon model. However, since the simple physics model's vertical diffusion is a first order closure model, no further variables need to be stored and the interpolation of any internal physics variables onto the adaptive grid is not required. \subsection{The workflow of the interface and physics call} % To illustrate the sequence of steps taken by the interface Fig.~\ref{InterfaceProgram} displays the workflow of the interface written to couple {\sc wavetrisk} to the physics package. Furthermore, to illustrate the physics-dynamics coupling challenges mentioned above, Fig.~\ref{PhysicsStepWorkflow} displays a flowchart of the algorithm used by the interface and required for the a single physics (split) step. The steps in black are the extra steps required to change prognostic variables to the required form, for example mass-weighted potential temperature to temperature or vice versa. The yellow step represents the call to the wrapper subroutine in the single column module added to physics package to allow the architecture of sending individual columns. \begin{figure}[ht] \centering \includegraphics[scale=.27]{Figures/Interface_Program.png} \caption[Interface Program Work Flow]{A flow chart of the interface (physics-dynamics) coupling the physics and dynamics. The workflow shows the the steps of the simulation until the end.} \label{InterfaceProgram} \end{figure} \begin{figure}[ht] \centering \includegraphics[scale=.18]{Figures/Physics_SplitStep_Workflow.png} \caption[Workflow of the {\sc wavetrisk} Physics Split Step Call]{A flow chart of the physics step. Steps in black contain overhead concerning the conversion of the prognostic variables from one type to another (i.e.\ mass-weighted potential temperature to temperature ($K$)). The yellow process represents the physics call of the added single column module of the physics.} \label{PhysicsStepWorkflow} \end{figure} \noappendix \authorcontribution{F Hourdin developed the original simple physics model in his PhD thesis. GCJ, NKRK prepared the paper with contributions from TD. NKRK developed the model code with contributions from GCJ and TD for the dry physics module. GCJ and NKRK performed the model simulations using {\sc wavetrisk-atmosphere}.} \competinginterests{The contact author has declared that none of the authors has any competing interests.} \begin{acknowledgements} NKRK and GCJ were supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) under grant no. RGPIN-2024-05282. \end{acknowledgements} \bibliographystyle{copernicus} \bibliography{bibliography} \end{document}