Implementation of the Community Earth System Model (cesm) Version 1.2.1 as a New Base Model into Version 2.50 of the Messy Framework

The Community Earth System Model (CESM1), maintained by the United States National Centre for Atmospheric Research (NCAR) is connected with the Modular Earth Submodel System (MESSy). For the MESSy user community, this offers many new possibilities. The option to use the Community Atmosphere Model (CAM) atmospheric dynamical cores, especially the state-of-the-art spectral element (SE) core, as an alternative to the ECHAM5 spectral transform dynamical core will provide scientific and computational advances for atmospheric chemistry and climate modelling with MESSy. The well-established finite volume core from CESM1(CAM) is also made available. This offers the possibility to compare three different atmospheric dynamical cores within MESSy. Additionally, the CESM1 land, river, sea ice, glaciers and ocean component models can be used in CESM1/MESSy simulations, allowing the use of MESSy as a comprehensive Earth system model (ESM). For CESM1/MESSy setups , the MESSy process and diagnostic submodels for atmospheric physics and chemistry are used together with one of the CESM1(CAM) dynamical cores; the generic (infrastructure) submodels support the atmospheric model component. The other CESM1 component models, as well as the coupling between them, use the original CESM1 infrastructure code and libraries; moreover, in future developments these can also be replaced by the MESSy framework. Here, we describe the structure and capabilities of CESM1/MESSy, document the code changes in CESM1 and MESSy, and introduce several simulations as example applications of the system. The Supplements provide further comparisons with the ECHAM5/MESSy atmospheric chemistry (EMAC) model and document the technical aspects of the connection in detail.


Introduction
Increasing scientific and societal interest in understanding and forecasting the state of the atmosphere, oceans, land and ice has led to the development of Earth system models (ESMs).The Community Earth System Model (CESM1; Hurrell et al., 2013) is a fully coupled global climate model, which has integrated individual Earth system component models, using a coupler and a generic IO library, but otherwise modifying the component models as little as possible.CESM1 has shown to be a very useful tool for many types of studies; see, for example, the special issue on CCSM (Community Climate System Model) and CESM in the Journal of Climate. 1 The Modular Earth Submodel System (MESSy) uses a different approach.The code is organized in four layers: a base model of any level of complexity is complemented by a base model interface layer.A further interface layer to the submodels makes it possible to keep process submodels as distinct as possible in the submodel core layer.For the A. J. G. Baumgaertner et al.: CESM1/MESSy ECHAM5/MESSy atmospheric chemistry (EMAC) model, the base model ECHAM5 provides only the dynamical core, including advection; all physics parametrizations have been recoded or replaced by submodels, and infrastructure code has been recoded or replaced by generic infrastructure submodels.For a list of available submodels, see Table 1 in Jöckel et al. (2010) or the MESSy website. 2ere, we have implemented CESM1 (version 1.2.1) as an additional base model for MESSy (implemented into MESSy version 2.50), similar to the implementation of ECHAM5.Note, however, that CESM1 provides a much larger number of process descriptions of all components of the Earth than ECHAM5.This means that much larger portions of the CESM1 code are still used in a CESM1/MESSy simulation.Here, we present test simulations using MESSy atmospheric physics and chemistry submodels for the atmosphere, with execution and data handling done by MESSy generic interface submodels, using one of the CESM1(CAM5) atmospheric dynamical cores, and CESM1 component models for ocean, land, ice and rivers.
The code integration can be seen from a MESSy or CESM user point of view.For MESSy users, CESM1/MESSy offers additional state-of-the art atmospheric dynamical cores, as well as the ability to couple with other component models.
As the development was aimed at MESSy users, the code structure, set-up design, configuration and script environment are analogous to ECHAM5/MESSy.For CESM users, CESM1/MESSy offers the opportunity to use an independent physics and chemistry suite, replacing the CAM physics and chemistry.

The Modular Earth Submodel System
The MESSy (Jöckel et al., 2005(Jöckel et al., , 2010)), maintained by the MESSy consortium, defines a strategy for building comprehensive ESMs from process-based modules, the so-called submodels.Technically, MESSy comprises standard interfaces to couple the different components, a simple coding standard and a set of submodels coded accordingly.The code is organized into four different layers: -The base model layer (BML) can be a model of arbitrary complexity starting from a global climate model (GCM) (as CESM1 or ECHAM5), to regional climate models (RCMs; such as COSMO) to models spanning the basic entity of the process (i.e. a box model for atmospheric chemistry or a column model for a convection model).
-The base model interface layer (BMIL) comprises the base-model-specific implementation of the MESSy infrastructure.
-The submodel interface layer (SMIL) represents the connector of a specific process to the infrastructure (BMIL).
-The submodel core layer (SMCL) comprises the basemodel-independent implementation of a specific process in the Earth system, or of a diagnostic tool of the model system.It uses data provided via its SMIL and returns data back via its SMIL to other submodels and/or the base model.
Coupled to the base model ECHAM5, MESSy has proven as a useful framework for atmospheric chemistry and physics studies.An up-to-date list of publications using the model is available at http://messy-interface.org.The layer structure described above makes comparisons of physics parametrizations a straightforward task; see, for example, Tost et al. (2006b).
For the second MESSy development cycle, which is comprehensively documented by Jöckel et al. (2010), complete independence of ECHAM5 was achieved by several new generic submodels.This has been exploited, for example, by the COSMO/MESSy development (Kerkweg and Jöckel, 2012a, b), for CMAT/MESSy (Baumgaertner et al., 2013a), and is also used here to connect to the CESM1 Earth system model.The CESM1 code was implemented into MESSy version 2.50, yielding an intermediate version 2.50+.The modifications will be made available in upcoming versions.

The Community Earth System Model
The Earth system model CESM1 (version 1.2.1) is a fully coupled global climate model.The physics-based models that serve for the different Earth system components are the Community Atmosphere Model (CAM), the Community Land Model (CLM), the sea ice model Community Ice CodE (CICE), the ocean model Parallel Ocean Program (POP), the land-ice model Community Ice Sheet Model (Glimmer-CISM), and the River Transport Model (RTM).As an alternative to the physics-based models, climatological data models are provided for each component.The models are coupled through the CESM1 coupler (CPL7), which uses the Model Coupling Toolkit (MCT).For a specific simulation, the user can choose a so-called component set, which describes the used model, model version as well as specific settings for each component.
The atmosphere component, CAM5, provides a set of physics parametrizations, and several dynamical cores, which also include advection.While CAM5 provides four different cores, we describe only the cores implemented in CESM1/MESSy, the CAM5 default finite volume (FV) core and the new spectral element (SE) core.The FV dynamics were initially developed by NASA's Data Assimilation Office (DAO).The discretization is local and entirely in physical space.In the horizontal, it uses a flux-form semi-Lagrangian scheme (Lin andRood, 1996, 1997), whereas the vertical discretization is quasi-Lagrangian.For more details, see the CAM5 description,3 Sect.3.1.
The SE dynamical core originates from the High-Order Method Modeling Environment (HOMME; Dennis et al., 2005).More specifically, SE uses a continuous Galerkin spectral finite element method (Taylor et al., 2009;Fournier et al., 2004;Thomas and Loft, 2005;Wang et al., 2007;Taylor and Fournier, 2010).It is currently implemented for a cubed-sphere grid, although the core can in principle be employed for fully unstructured quadrilateral meshes.The main advantages compared to traditional approaches are its scalability up to 10 5 compute cores, which is useful for current and future computing architectures, and local energy conservation on top of mass and potential vorticity conservation.Also, no polar filters are required since the grid is quasi-uniform.A detailed description and further references are given in the CAM5 description (Sect.3.2).A recent publication by Bacmeister et al. (2014) discusses some improvements, but also some problems at very high-resolution (0.23 • latitude × 0.31 • longitude) simulations.
CESM1 time stepping (so-called run alarms) can be chosen through the driver namelist, but most component sets use 30 min for all components except for the ice sheet model.For CAM, the 30 min time step applies to the physics parametrization, whereas the dynamical cores can have shorter time steps, depending on the horizontal resolution.This is achieved through substepping within the coupling to the core.The coupling is performed in a time-split manner for both FV and SE.For details see Sect. 2 in the CAM5 description.

Technical implementation of CESM1/MESSy
The development of CESM1/MESSy was driven by two goals: first, to provide the state-of-the art SE dynamical core to the MESSy user community, and second to provide further components (land, ice, etc.) to MESSy simulations, making it a comprehensive Earth system model.The strategy chosen to achieve both goals was to implement the entire CESM1 code as a base model into MESSy, analagous to the implementation of the base model ECHAM5.A diagram of the CESM1/MESSy structure is shown in Fig. 1.It indicates the MESSy layer structure as described above, the basics of the call structure between CESM1 and MESSy submodels, and basics of the data exchange.
The entire CESM1 repository is taken over as part of MESSy, which makes updates to newer versions of CESM1 straight forward.All changes to the CESM1 Fortran code are encapsulated using preprocessor commands: The CESM1 model components including the coupler can still be used in the CESM1/MESSy configuration; only the CAM5 process parametrizations are disabled and replaced by the MESSy atmospheric physics and chemistry.
The MESSy main control interface is called from the CCSM driver module ccsm_comp_mod, the CAM module atm_comp_mct and for the row loop in physpkg.The module atm_comp_mct is the outermost module in CAM, and also takes care of the coupling to the other component models.Most calls could also be moved to the ccsm_comp module, which controls the CESM1 time stepping and call the different component models, but since MESSy currently only replaces the CAM5 atmospheric physics and chemistry, atm_comp_mct is the most straightforward place in the code.For an overview of the call structure, see Fig. 1 in the Supplement "Implementation Documentation".
For MESSy, the submodel core layer remains unchanged, but the generic BMIL, as well as the SMIL, is modified.For submodels with a generic SMIL the modifications are encapsulated using preprocessor statements (#ifdef CESM1).For most SMIL modules no changes or very minor adjustments were necessary.For the remaining submodels4 that are more base-model-specific, new SMIL modules were created based on the respective ECHAM5 SMIL.
The following subsections provide an overview of these changes in MESSy and CESM1.

Time integration
CESM1/MESSy employs an explicit Euler time integration for the atmosphere with long time steps for the physics and chemistry, and higher-order types of integration (e.g.Runge-Kutta for SE) in the dynamical cores.The dynamical cores use sub-cycling for shorter integration times.Note that this is different to ECHAM5/MESSy, which uses leapfrog integration and a time filter.Sub-time stepping in MESSy is used for chemistry submodels such as MECCA and SCAV, whereas longer time steps (n • t) are used for radiation; i.e. the radiation submodel is called less frequently.
For CESM1/MESSy, the CAM time-integration scheme was adopted.Note however that while CAM performs a time integration after every individual physics process, allowing to use the state x for each process, MESSy performs a time integration at the end of every time step, but explicitly integrates required variables in every submodel, x + dx/dt • t.When using the SE core, the CESM1/MESSy integration is applied to temperature, winds, specific humidity, cloud water (liquid and ice), and trace gas mixing ratios.The coupling between the physics and dynamics is a time-split cou- pling, where physical and dynamical core time-integration components are calculated sequentially.This is equivalent to the coupling of the FV and SE cores with the CAM physics, which is described in more detail in Sect. 2 of the CAM5 description.

Data representation, input/output
MESSy uses representations (see Jöckel et al., 2010, for an explanation of the terminology) that describe the geometric structure of data objects based on dimensions.For CESM1/MESSy, representations analogous to the ECHAM5/MESSy grid point (or Eulerian) representations are used for all atmosphere data for both the FV and SE cores.All data are stored in CHANNEL objects, which contain the data fields, the object's representation, and metadata.The CHANNEL infrastructure module (Jöckel et al., 2010) also controls the model output and writing of restart files.
A namelist file gives the user full control over the output data.
For data import from files, MESSy provides the infrastructure submodel IMPORT.IMPORT is namelist controlled, and provides the data regridded to the required representation as channel objects, which every submodel can access through coupling with the respective channel objects.For CESM1/MESSy, this infrastructure is used for all data import.The TRACER submodel (Jöckel et al., 2008), which provides the handling of atmospheric trace gas variables, directly uses the NCREGRID (Jöckel, 2006) or GRID_TRAFO submodels for initialization of the tracers.Note that currently for the SE core, which employs an unstructured grid, all im-ported data, including those for tracer initialization, have to be provided on the grid used for the simulation.
In CESM1, explicit-shape arrays are used, such that the horizontal and vertical resolution as well as the number of tracers have to be selected before compilation.MESSy, in contrast, applies a dynamical memory management at run time.However, the replacement of CESM1 explicit-shape arrays by pointers in the dynamical cores has so far only been implemented for the tracers.The horizontal and vertical resolution have to be specified when MESSy is configured; for example, CESM1HRES=1.9x2.5 CESM1VRES=26 have to be added to the call of configure.
For the grid point representation, each process (MPI task) has its own set of rows and columns.The only difference is that for ECHAM the number of columns in the last row is in general different to the other rows, whereas in CAM the number of columns can be different for all rows.For the base model interfaces and submodel interfaces, this requires a distinction as detailed in the documentation Supplement.

Coupling to other component models through MCT
CESM1 uses the open-source MCT (Larson et al., 2005;Jacob et al., 2005), maintained by the Argonne National Laboratory.For CESM1/MESSy, this coupling is left in place, although in the future a coupling through the MESSy Multi-Model Driver (MMD, Kerkweg and Jöckel, 2012b) is anticipated.The MESSy channel objects for the atmospheric component are coupled to the data of the other component model analogously to CAM coupling.For a list of variables and the technical documentation, see the Supplement.

Parallelization
CESM1 is structured to have all component models handle their parallelization separately, giving each component model its own set of processors, which can be controlled via the namelist drv_in.The CAM physics and dynamical cores also have separate parallelization, depending on the employed grid.Due to the similarity of the MESSy and CAM physics data representation, the parallelization routines of the CAM physics are employed also for MESSy submodels.Technically, this means that the MPI infrastructure submodel uses the spmd_utils and phys_grid modules from CAM for the low-level gather/scatter routines.Specifically, the parallel data types, gather (gather_chunk_to_field) and scatter (scatter_field_to_chunk) subroutines available from spmd_utils, which directly uses the MPI library, are employed.In comparison, for ECHAM5/MESSy simulations the MPI submodel uses ECHAM5's mo_mpi low-level routines.

Namelists and scripts
Similar to CESM1, CESM1/MESSy also offers a large variety of set-up possibilities.In CESM1, there are a number of evaluated set-ups, so-called component sets (see Sect. 2.2).MESSy also offers several set-ups that the user can choose for a simulation, and that can be easily modified depending on the scientific requirements.
A variety of scripts support the CESM1 model set-up, which generate for instance the makefiles and namelists.MESSy uses autoconf/configure/make utilities, and a single script for runcontrol (xmessy_mmd).Run-time options are set in well-documented namelist files directly.The model comes with several namelist set-ups for different model configurations.
Instead of the automatic namelist generation in CESM1, the MESSy namelist set-ups contain some variables that are replaced by the runscript, for example, for resolutiondependent filenames, or start/stop dates.

Trace constituents and mixing ratios
In general, atmospheric air masses can be treated to include (wet) or exclude (dry) water vapour.Both in CAM and MESSy, specific humidity is treated as wet mass mixing ratio, i.e. water mass with respect to total air mass [kg kg −1 = (kg H 2 O)/(kg total air)].Also, in both CAM and MESSy cloud liquid and ice are treated as mass mixing ratios with respect to dry air [kg kg −1 = (kg H 2 O)/(kg dry air)].In MESSy, other trace constituents are treated as dry volume mixing ratio, i.e. [mol (mol of dry air) −1 ].The dynamical cores FV and SE both expect wet mass mixing ratios for advection.Therefore, advected trace constituents are converted before and after the advection through the dynamical core.

Vertical diffusion
The current suite of MESSy physical parametrization submodels does not include a submodel for vertical diffusion.For ECHAM5/MESSy, vertical diffusion is treated by the ECHAM5 base model.For CESM1/MESSy, the vertical diffusion code of CAM5 was restructured as a MESSy submodel (VERTDIFF).However, both models use a similar approach.In both models, the free atmosphere diffusion coefficients are estimated using the gradient Richardson number.For the boundary layer, they both use a Monin-Obukov similarity approach.The vertical diffusion equation is solved using an implicit method.For details of the implementation, see the VERTDIFF documentation in the Supplement.
3. maCMAC-FV: CESM1/MESSy with finite volume core at 1.9 The trace gas emissions and prescribed mixing ratios of long-lived trace gases (TNUDGE; see Kerkweg et al., 2006) are all from the year 2000.All simulations were performed for one model year, without spin-up using initializations from existing simulations.Note that the maEMAC simulation contains a more complete set of trace gas emissions than the CESM1/MESSy simulations.The respective namelist setups are provided in the Supplement.Baumgaertner (2015) contains a comparison of these set-ups for all major output variables.The following subsections present several evaluation examples.

Using the global electric circuit for model evaluation
The global electric circuit (GEC) is a system of currents spanning the globe.The currents are generated by thunderstorms and electrified clouds, whereas the spatial and temporal distribution of conductivity determines the potential and current distribution in the fair-weather atmosphere.For a recent review on the GEC, see Williams and Mareev (2014).
The physical state of the atmosphere determines the current generation as well as conductivity.Therefore, for a model to simulate the state and variability of the GEC correctly depends on its ability to reproduce temperature, humidity, air density, cloud cover, trace gas transport and a correct representation of convection.Modelling studies on the GEC with CESM1 are presented by Lucas et al. (2015) and Baumgaertner et al. (2013b).
We use the GEC current generation as well as conductivity as a way to collectively evaluate the operation and coupling amongst the various submodels involved in CESM1/MESSy simulations.Since the derived variables combine several ba-  sic aspects such as temperature, pressure and tracer transport, the GEC offers a way to evaluate several variables at the same time.Of course, this does not substitute a full evaluation, but rather presents an example application.
Both current generation parametrization and the conductivity have been implemented as a diagnostic MESSy submodel named GEC.
We parametrize current generation analogously to Kalb et al. (2016), who found that convection updraft mass flux averaged between 200 and 800 hPa is correlated with measured electrified cloud and thunderstorm occurrence.The MESSy submodel CONVECT offers eight different convection schemes, all providing updraft mass flux.Here, we show results from several additional CESM1/MESSy and EMAC sensitivity simulations that use the Tiedtke scheme (Tiedtke, 1989) with Nordeng closure (Nordeng, 1994), and the Bechtold scheme (Bechtold et al., 2001), respectively.The most critical aspect of GEC source current is the diurnal cycle, referred to as the Carnegie curve from electric field measurements in fair-weather regions.Figure 2 shows the total current composite mean, averaged over 45 • S to 45 • N as a function of universal time, using hourly stored data for one simulation year, as well as the Carnegie E-field measurements, provided by Harrison (2013).In general, the simulations reproduce a diurnal cycle similar to the Carnegie data.However, the current peaks too early in the day for all simulations, which is a common problem with convection parametrizations (see e.g.Lucas et al., 2015).Only the simulation using the Bechtold convection scheme (blue) has its maximum at 18:00 UT, close to the peak in the Carnegie data.
Conductivity is calculated similar to the approach described by Baumgaertner et al. (2013b), B13 hereafter, who used CESM1(Whole Atmosphere Community Climate Model -WACCM) to study spatial and temporal conductivity variability.Conductivity is proportional to ion pair concentrations, n, and positive/negative ion mobilities, µ +/− , and is defined as where e is the elementary charge, and positive and negative ion concentrations are assumed to be equal.Ion concentration is given by where the ion production rate is q, the ion-ion recombination rate α and the effective loss of ions by aerosol particles with rate i, r β(r i )S(i, r).
Here, we use the same parametrizations for galactic cosmic ray (GCR) ion production, mobility, and ion-ion recombination as described by B13.Lower atmosphere ionization sources include 222 Rn (Radon), obtained from the DRADON submodel, and further radioactive decay sources, also parametrized in the same way as presented by B13.While the aerosol attachment rate could be calculated using MESSy aerosol submodels, for consistency with B13 we use the same input data sets from CESM1(WACCM) simulations with CARMA (Community Aerosol and Radiation Model for Atmospheres).Note that clouds are not introduced as additional resistors in the present study.Column resistance is defined as the vertical integral of the reciprocal of conductivity (see e.g.B13 and references therein): where dz is the model layer thickness, which depends on height and geographic location.Figure 3 presents January column resistance from the maCMAC-FV (left) and maEMAC (right) simulations.Higher resistance at low latitudes, specifically at low geomagnetic latitudes, is due to the smaller GCR ionization.Mountains lead to a decrease in column resistance because there is less atmosphere between the mountain and the upper boundary.Terrestrial emissions of Radon decrease column resistance over land compared to ocean.Radon has a halflife of approximately 4 days, therefore advection of Radon from land to ocean can lead to elevated ionization rates near the coasts, so the transition is usually smooth.

Trace constituents and atmospheric chemistry
As a further example, we compare surface-tropospheric hydroxyl (OH), an important atmospheric cleaning agent, as well as stratospheric ozone concentrations.Note that the chosen variables and types of comparisons have no scientific justification for a full model evaluation, but are only example applications.
Zonal mean surface OH number concentrations are shown in Fig. 4 for the CMAC-FV (left), CMAC-SE (middle) and maEMAC (right) simulations for 1 year.As the CESM1/MESSy simulations are free running, different synoptic meteorologies lead to some differences on timescales of weeks, but overall the expected annual variations are present in all three simulations.This confirms the functionality of the emission, boundary condition and chemistry integration scheme.Tropospheric OH concentrations are important for the tropospheric methane lifetime (τ CH 4 ).With τ CH 4 =7.61 years, CMAC-FV is more reactive than maEMAC (τ CH 4 = 8.24 years), whereas CMAC-SE is less reactive (τ CH 4 = 10.46 years).This finding highlights the large influence of the dynamical core.
Figure 5 depicts the zonal mean ozone (top panel) and the column ozone (bottom panel) between 60 and 90 • S. Again, agreement is found between the maCMAC-FV (left) and maEMAC (right) simulations, showing the principal functionality especially of the dynamics, transport, and chemistry systems.However, the expected polar spring (September/October) ozone loss around 50 hPa is only shown by maEMAC.There is also more column ozone evident in maCMAC-FV than in maEMAC.Note that for low and midlatitudes the ozone column is very similar with no discernible bias (not shown).

Conclusions
CESM1 is connected to the Modular Earth Submodel System (MESSy) as a new base model.This allows MESSy users the option to utilize either the state-of-the art spectral element dynamical core or the finite volume core of CESM1.Additionally, this makes several other component models available to MESSy users.As example applications, an initial evaluation with respect to the global electric circuit, which offers a unique opportunity for evaluating a range of atmospheric parameters under a single scientific aspect, was performed.Good agreement between the CESM1/MESSy simulations and ECHAM5/MESSy is found.Similarly, an exemplary comparison of surface OH and Antarctic ozone shows the principal functionality of the atmospheric chemistry in the model.A broader evaluation will be published elsewhere.The developments and experiences will be useful also for further MESSy extensions, for example with the new ICON (Icosahedral non-hydrostatic) GCM (Zängl et al., 2015).
Further technical work on CESM1/MESSy is likely to include the following: -The coupling between the CESM1 component models with MCT can be replaced by the MESSy infrastructure.
-The CESM1 component models can be adapted to use the MESSy CHANNEL infrastructure submodel for memory management and data output.
-The CAM5 physical parametrizations can be implemented as MESSy submodels such that they can be used as alternative submodels for the current parametrization suite.
-The new MESSy infrastructure submodel GRID (Kerkweg and Jöckel, 2015) for regridding can be adapted for handling the SE data.

Figure 5 .
Figure 5. Top panel: zonal mean ozone (µmol mol −1 ) averaged between 60 and 90 • S for the year 2000 from the maCMAC-FV (left) and maEMAC (right) simulations.Bottom panel: column ozone for the same region.

Table 1 .
Jöckel et al. (2010)diagnostic submodels used in the simulations presented here.For a full list of available submodels, see Table1inJöckel et al. (2010)or the MESSy website (http://www.messy-interface.org/).