Towards in-situ visualization integrated earth system models : RegESM 1 . 1 regional modelling system

The data volume being produced by regional and global multi-component earth system models are rapidly increasing due to the improved spatial and temporal resolution of the model components, sophistication of the used numerical models in terms of represented physical processes and their non-linear complex interactions. In particular, very short time steps have to be defined in multi-component and multi-scale non-hydrostatic modelling systems to represent the evolution of the fastmoving processes such as turbulence, extra-tropical cyclones, convective lines, jet streams, internal waves, vertical turbulent 5 mixing and surface gravity waves. Consequently, the used small time steps cause extra computation and disk I/O overhead in the used modelling system even if today’s most powerful high-performance computing and data storage systems are being considered. Analysis of the high volume of data from multiple earth system model components at different temporal and spatial resolution also poses a challenging problem to efficiently perform integrated data analysis of the massive amounts of data by relying on the conventional post-processing methods available today. This study basically aims to explore the 10 feasibility and added value of integrating existing in-situ visualization and data analysis methods with the model coupling framework (ESMF) to increase interoperability between multi-component simulation code and data processing pipelines by providing easy to use, efficient, generic and standardized modeling environment for earth system science applications. The new data analysis approach enables simultaneous analysis of the vast amount of data produced by multi-component regional earth system models (atmosphere, ocean etc.) during the run process. The methodology aims to create an integrated modeling 15 environment for analyzing fast-moving processes and their evolution in both time and space to support better understanding of the underplaying physical mechanisms. The state-of-art approach can also be used to solve common problems in earth system model development workflow such as designing new sub-grid scale parametrizations (convection, air-sea interaction etc.) that requires inspecting the integrated model behavior in a higher temporal and spatial scale during the run or supporting visual debugging of the multi-component modeling systems, which usually are not facilitated by existing model coupling libraries 20 and modeling systems.


Introduction
The multi-scale and inherently coupled Earth System Models (ESMs) make them challenging to study and understand.Rapid developments in Earth system science, as well as in high-performance computing and data storage systems, have enabled fully coupled regional or global ESMs to better represent relevant processes, complex climate feedbacks, and interactions among the coupled components.In this context, regional ESMs are employed when the spatial and temporal resolution of the global climate models are not sufficient to resolve local features such as complex topography, land-sea gradients and the influence of human activities in a smaller spatial scale.Along with the development of the modeling systems, specialized software libraries for the model coupling become more and more critical to reduce the complexity of the coupled model development and increase the interoperability, reusability, and efficiency of the existing modeling systems.Currently, the existing model coupling software libraries have two main categories: couplers and coupling frameworks.
Couplers are mainly specialized in performing specific operations more efficiently and quickly such as coordination of components and interpolation among model components.For example, OASIS3 (Valcke, 2013) uses multiple executable approaches for coupling model components but sequentially performing internal algorithms such as sparse matrix multiplication (SMM) operation for interpolation among model grids become a bottleneck along with increased spatial resolution of the model components.To overcome the problem, OASIS4 uses parallelism in its internal algorithms (Redler et al., 2010), but OASIS3 coupler interfaced with the Model Coupling Toolkit (MCT; Jacob et al., 2005;Larson et al., 2005) to develop OASIS3-MCT that provides a parallel implementation of interpolation and data exchange.Besides generic couplers like OASIS, domainspecific couplers such as Oceanographic Multi-purpose Software Environment (OMUSE; Pelupessy et al., 2017) that aims to provide a homogeneous environment for ocean modeling to make verification of simulation models with different codes and numerical methods and Community Surface Dynamics Modeling System (CSDMS; Overeem et al., 2013) to develop integrated software modules for modeling of Earth surface processes are introduced.
A coupling framework is as an environment for coupling model components through a standardized calling interface and aims to reduce the complexity of regular tasks such as performing spatial interpolation across different computational grids and transferring data among model components to increase the efficiency and interoperability of multi-component modeling systems.Besides, the synchronization of the execution of individual model components, a coupling framework can simplify the exchange of metadata related to model components and exchanged fields through the use of existing conventions such as CF (Climate and Forecast) convention.The Earth System Modeling Framework (ESMF) is one of the most famous examples of this approach (Theurich et al., 2016).The ESMF consists of a standardized superstructure for coupling components of Earth system applications through a robust infrastructure of high-performance utilities and data structures that ensure consistent component behavior (Hill et al., 2004).The ESMF framework is also extended to include the National Unified Operational Prediction Capability (NUOPC) layer.The NUOPC layer simplifies component synchronization and run sequence by providing additional programming interface between coupled model and ESMF framework through the use of NUOPC "cap".In this case, a NUOPC "cap" is a Fortran module that serves as the interface to a model when it is used in a NUOPC-based coupled system.
The term "cap" is used because it is a small software layer that sits on top of a model code, making calls into it and exposing model data structures in a standard way.In addition to generic modeling framework like ESMF, Modular System for Shelves and Coasts (MOSSCO; Lemmen et al., 2018) creates a state-of-art domain and process coupling system by taking advantage of both ESMF and Framework for Aquatic Biogeochemical Models (FABM; Bruggeman and Bolding, 2014) for marine coastal Earth system community.
The recent study of Alexander and Easterbrook (2015) to investigate the degree of modularity and design of the existing global climate models reveals that the majority of the models use central couplers to support data exchange, spatial inter-polation, and synchronization among model components.In this approach, direct interaction does not have to occur between individual model components or modules, since the specific coupler component manages the data transfer.This approach is also known as the hub-and-spoke method of building a multi-component coupled model.A key benefit of using a hub-and-spoke approach is that it creates a more flexible and efficient environment for designing sophisticated multi-component modeling system regarding represented physical processes and their interactions.The development of the more complex and high-resolution modeling systems leads to an increased demand for both computational and data storage resources.In general, the high volume of data produced by the numerical modeling systems may not allow storing all the critical and valuable information to use later, despite recent advances in storage systems.As a result, the simulation results are stored in a limited temporal resolution (i.e., monthly averages), which are processed after numerical simulations finished (post-processing).The poor representation of the results of numerical model simulations prevents to analyze the fast-moving processes such as extreme precipitation events, convection, turbulence and non-linear interactions among the model components in a high temporal and spatial scale with the conventional post-processing approach.
The analysis of leading high-performance computing systems reveals that the rate of disk input-output (I/O) performance is not growing at the same speed as the peak computational power of the systems (Ahern, 2012;Ahrens, 2015).The recent report of U.S. Department of Energy (DOE) also indicates that the expected rate of increase in I/O bandwidth (100 times) will be slower than the peak system performance (500 times) of the new generations of exascale computers (Ashby et al., 2010).Besides, the movement of large volumes of data across relatively slow network bandwidth servers fails to match the ultimate demands of data processing and to archive tasks of the present high-resolution multi-component ESMs.As a result, the conventional post-processing approach has become a bottleneck in monitoring and analysis of fast-moving processes that require very high spatial resolution, due to the present technological limitations in high-performance computing and storage systems (Ahrens et al., 2014).In the upcoming computing era, state-of-art new data analysis and visualization methods are needed to overcome the above limitations evocatively.
Besides the conventional data analysis approach, the so-called in-situ visualization and co-processing approaches allow researchers to analyze the output while running the numerical simulations simultaneously.The coupling of computation and data analysis helps to facilitate efficient and optimized data analysis and visualization pipelines and boosts the data analysis workflow.Recently, a number of in-situ visualization systems for analyzing numerical simulations of Earth system processes have been implemented.For instance, the ocean component of Model for Prediction Across Scales (MPAS) has been integrated with an image-based in-situ visualization tool to examine the critical elements of the simulations and reduce the data needed to preserve those elements by creating a flexible work environment for data analysis and visualization (Ahrens et al., 2014;O'Leary et al., 2016).Additionally, the same modeling system (MPAS-Ocean) has been used to study eddies in large-scale, high-resolution simulations.In this case, the in-situ visualization workflow is designed to perform eddy analysis at higher spatial and temporal resolutions than available with conventional post-processing facing storage size and I/O bandwidth constraints (Woodring et al., 2016).Moreover, a regional weather forecast model (Weather Research and Forecasting Model; WRF) has been integrated with in-situ visualization tool to track cyclones based on an adaptive algorithm (Malakar et al., 2012).Despite the lack of generic and standardized implementation for integrating model components with in-situ visualization tools, the previous studies have shown that in-situ visualization can produce analyses of simulation results, revealing many details in an efficient and optimized way.It is evident that more generic implementations could facilitate smooth integration of the existing standalone and coupled ESMs with available in-situ visualization tools (Ahrens et al., 2005;Ayachit, 2015;Childs et al., 2012) and improve interoperability between such tools and non-standardized numerical simulation codes.
The main aim of this paper is to explore the added value of integrating in-situ analysis and visualization methods with a model coupling framework (ESMF) to provide in-situ visualization for easy to use, generic, standardized and robust scientific applications of Earth system modeling.The implementation allows existing ESMs coupled with the ESMF library to take advantage of in-situ visualization capabilities without extensive code restructuring and development.Moreover, the integrated model coupling environment allows sophisticated analysis and visualization pipelines by combining information coming from multiple ESM components (i.e., atmosphere, ocean, wave, land-surface) in various spatial and temporal resolutions.Detailed studies of fundamental physical processes and interactions among model components are vital to the understanding of complex physical processes and could potentially open up new possibilities for the development of ESMs.

The Design of the Modeling System
The RegESM (Regional Earth System Model; 1.1) modeling system can use five different model components to support many different modeling applications that might require detailed representation of the interactions among different Earth system processes (Fig. 1a-b).The implementation of the modeling system follows the hub-and-spoke architecture.The driver that is responsible for the orchestration of the overall modeling system resides in the middle and acts as a translator among model components (atmosphere, ocean, wave, river routing, and co-processing).In this case, each model component introduces its NUOPC cap to plug into the modeling system.The modeling system is validated in different model domains such as Caspian Sea (Turuncoglu et al., 2013), Mediterranean Basin (Surenkok and Turuncoglu, 2015;Turuncoglu and Sannino, 2017), and Black Sea Basin.

Atmosphere Models (ATM)
The flexible design of RegESM modeling system allows choosing a different atmospheric model component (ATM) in the configuration of the coupled model for a various type of application.Currently, two different atmospheric model is compatible with RegESM modeling system: (1) RegCM4 (Giorgi et al., 2012), which is developed by the Abdus Salam International Centre for Theoretical Physics (ICTP) and (2) the Advanced Research Weather Research and Forecasting (WRF) Model (ARW; Skamarock et al., 2005), which is developed and sourced from National Center for Atmospheric Research (NCAR).In this study, RegCM 4.6 is selected as an atmospheric model component because the current implementation of WRF coupling interface is still experimental and does not support coupling with co-processing component yet, but the next version of the modeling system (RegESM 1.2) will be able to couple WRF atmospheric model with co-processing component.The NUOPC cap of atmospheric model components defines state variables (i.e., sea surface temperature, surface wind components), rotates the winds relative to Earth, apply unit conversions and perform vertical interpolation to interact with the newly introduced co-processing component.

RegCM
The dynamical core of the RegCM4 is based on the primitive equation, hydrostatic version of the National Centre for Atmospheric Research (NCAR) and Pennsylvania State University mesoscale model MM5 (Grell, 1995).The latest version of the model (RegCM 4.6) also supports non-hydrostatic dynamical core to support applications with high spatial resolutions (< 10 km).The model includes two different land surface models: (1) Biosphere-Atmosphere Transfer Scheme (BATS; Dickinson et. al., 1989) and (2) Community Land Model (CLM), version 4.5 (Tawfik and Steiner, 2011).The model also includes specific physical parameterizations to define air-sea interaction over the sea and lake (one-dimensional lake model; Hostetler et al., 1993).The Zeng Ocean Air-Sea Parameterization (Zeng et al., 1998) is extended to introduce the atmosphere model as a component of the coupled modeling system.In this way, the atmospheric model can exchange both two and threedimensional fields with other model components such as an ocean, wave and river routing components that are active in an area inside of the atmospheric model domain as well as in-situ visualization component.

WRF
The WRF model consists of fully compressible non-hydrostatic equations, and the prognostic variables include the threedimensional wind, perturbation quantities of pressure, potential temperature, geo-potential, surface pressure, turbulent kinetic energy and scalars (i.e., water vapor mixing ratio, cloud water).The model is suitable for a broad range of applications and has a variety of options to choose parameterization schemes for the planetary boundary layer (PBL), convection, explicit moisture, radiation, and soil processes to support analysis of different Earth system processes.The PBL scheme of the model has a significant impact on exchanging moisture, momentum, and energy between air and sea (and land) due to the used alternative surface layer options (i.e., drag coefficients) in the model configuration.A few modifications are done in WRF (version 3.8.1)model itself to couple it with RegESM modeling system.These modifications include rearranging of WRF time-related subroutines, which are inherited from the older version of ESMF Time Manager API (Application Programming Interface) that was available in 2009, to compile model with the newer version of ESMF library (version 7.1.0)together with the older version that requires mapping of time manager data types between old and new versions.

Ocean Models (OCN)
The current version of the coupled modeling system supports two different ocean model components (OCN): (1) Regional Ocean Modeling System (ROMS revision 809; Shchepetkin and McWilliams, 2005;Haidvogel et al., 2008), which is developed and distributed by Rutgers University and (2) MIT General Circulation Model (MITgcm version c63s; Marshall et al., 1997a, b).In this case, ROMS and MITgcm models are selected due to their large user communities and different vertical grid representations.Although the selection of ocean model components depends on user experience and application, often the choice of vertical grid system has a determining role in some specific applications.For example, the ROMS ocean model uses terrain following (namely s-coordinates) vertical grid system that allows better representation of the coastal processes but MITgcm uses z levels and generally used for the applications that involve open oceans and seas.Similar to the atmospheric model component, both ocean models are slightly modified to allow data exchange with the other model components.In the current version of the coupled modeling system, there is no interaction between wave and ocean model components, which could be crucial for some applications (i.e., surface ocean circulation and wave interaction) that need to consider the two-way interaction between waves and ocean currents.The exchange fields defined in the coupled modeling system between ocean and atmosphere strictly depend on the application and the studied problem.In some studies, the ocean model requires heat, freshwater and momentum fluxes to be provided by the atmospheric component, while in others, the ocean component retrieves surface atmospheric conditions (i.e., surface temperature, humidity, surface pressure, wind components, precipitation) to calculate fluxes internally, by using bulk formulas (Turuncoglu et al., 2013).In the current design of the coupled modeling system, the driver allows selecting the desired exchange fields from the predefined list of the available fields.The exchange field list is a simple database with all known fields that can be exported or imported by the component.In this way, the coupled modeling system can be adapted to different applications without any code customizations in both the driver and individual model components.

ROMS
The ROMS is a three-dimensional, free-surface, terrain-following numerical ocean model that solves the Reynolds-averaged Navier-Stokes equations using the hydrostatic and Boussinesq assumptions.The governing equations are in flux form, and the model uses Cartesian horizontal coordinates and sigma vertical coordinates with three different stretching functions.The model also supports second, third and fourth order horizontal and vertical advection schemes for momentum and tracers via its preprocessor flags.

MITgcm
The MIT general circulation model (MITgcm) is a generic and widely used ocean model that solves the Boussinesq form of Navier-Stokes equations for an incompressible fluid.It supports both hydrostatic and non-hydrostatic applications with a spatial finite-volume discretization on a curvilinear computational grid.The model has an implicit free surface in the surface and partial step topography formulation to define vertical depth layers.The MITgcm model supports different advection schemes for momentum and tracers such as centered second order, third-order upwind and second-order flux limiters to support a variety of applications.The model used in the coupled modeling system is slightly modified by ENEA to allow data exchange with other model components.The detailed information about the regional applications of the MITgcm ocean model was initially described in the study of Artale et al. (2010) using PROTHEUS modeling system specifically developed for the Mediterranean Sea.

Wave Model (WAV)
Surface waves play a crucial role in the dynamics of PBL in the atmosphere and the currents in the ocean.Therefore, the wave component is included in the coupled modeling system to have a better representation of atmospheric PBL and surface conditions (i.e., surface roughness, friction velocity, wind speed).In this case, the wave component is based on WAM Cycle-4 (4.5.3-MPI).The WAM is a third-generation model without any assumption on the spectral shape (Monbaliu et al., 2000).It considers all the main processes that control the evolution of a wave field in deep water, namely the generation by wind, the nonlinear wave-wave interactions, and also white-capping.The model is initially developed by Helmholtz-Zentrum Geesthacht (GKSS, now HZG) in Germany.The original version of the WAM model is slightly modified to retrieve surface atmospheric conditions (i.e., wind speed components or friction velocity and wind direction) from the RegCM4 atmospheric model and send back calculated surface roughness.In the current version of the modeling system, wave component cannot be coupled with the WRF model due to the missing modifications in the WRF side.In the RegCM4, the received surface roughness is used to calculate air-sea transfer coefficients and fluxes over sea using Zeng ocean air-sea parameterization (Zeng et al., 1998).In this design, it is also possible to define a threshold for maximum roughness length (the default value is 0.02 m) and friction velocity (the default value is 0.02 m) in the configuration file of RegCM4 to ensure the stability of the overall modeling system.The initial results to investigate the added value of atmosphere-wave coupling in the Mediterranean Sea can be found in Surenkok and Turuncoglu (2015).

River Routing Model (RTM)
To simulate the lateral freshwater fluxes (river discharges) at the land surface and to provide river discharge to ocean model component, the RegESM modeling system uses Hydrological Discharge (HD, version 1.0.2) model developed by Max Planck Institute (Hagemann and Dumenil, 1998;Hagemann and Lydia, 2001).The model is designed to run in a fixed global regular grid with 0.5 • horizontal resolution using daily time series of surface runoff and drainage as input fields.In this case, the model uses the pre-computed river channel network to simulate the horizontal transport of the runoff within model watersheds using different flow processes such as overland flow, baseflow and riverflow.The river routing model (RTM) plays an essential role in the freshwater budget of the ocean model by closing the water cycle between the atmosphere and ocean model components.
The original version of the model is slightly modified to support interaction with the coupled model components.To close water cycle between land and ocean, model retrieves surface and sub-surface runoff from the atmospheric component (RegCM or WRF) and provides estimated river discharge to the selected ocean model component (ROMS or MITgcm).In the current design of the driver, rivers can be represented in two different ways: (1) individual point sources that are vertically distributed to model layers, and (2) imposed as freshwater surface boundary condition like precipitation (P) or evaporation minus precipitation (E-P).In this case, the driver configuration file is used to select the river representation type (1 or 2) for each river individually.The first option is preferred if river plumes need to be defined correctly by distributing river discharge vertically among the ocean model vertical layers.The second option is used to distribute river discharge to the ocean surface when there is a need to apply river discharge to a large areal extent close to the river mouth.In this case, a special algorithm implemented in NUOPC cap of ocean model components (ROMS and MITgcm) is used to find affected ocean model grids based on the effective radius (in km) defined in the configuration file of the driver.

The Driver: RegESM
The RegESM (version 1.1) is completely redesigned and improved version of the previously used and validated coupled atmosphere-ocean model (RegCM-ROMS) to study the regional climate of Caspian Sea and its catchment area (Turuncoglu et al., 2013).To simplify the design and to create more generic, extensible and flexible modeling system that aims to support easy integration of multiple model components and applications, the RegESM uses a driver to implement the hub-and-spoke approach.In this case, all the model components are combined using ESMF (version 7.1.0)framework to structure coupled modeling system.The ESMF framework is selected because of its unique online re-gridding capability, which allows the driver to perform different interpolation types (i.e., bilinear, conservative) over the exchange fields (i.e., sea surface temperature, heat and momentum fluxes) and the NUOPC layer.The NUOPC layer is a software layer built on top of the ESMF.It refines the capabilities of ESMF by providing a more precise definition of what it means for a model to be a component and how components should interact and share data in a coupled system.The ESMF also provides the capability of transferring computational grids in memory among the model components, which has critical importance in the integration of the modeling system with a co-processing environment (see also Sect.3).The RegESM modeling system also uses ESMF and NUOPC layer to support various configuration of component interactions such as defining multiple coupling time steps among the model components.
An example configuration of the four-component (ATM, OCN, RTM, and WAV) coupled modeling system can be seen in Fig. 2. In this case, the RTM component runs in a daily time step (slow) and interacts with ATM and OCN components, but ATM and OCN components can interact each other more frequently (fast) such as every three hours.
The interaction (also called as run sequences) among the model components and driver are facilitated by the connector components provided by NUOPC layer.Connector components are mainly used to create a link between individual model components and driver.In this case, the number of active components and their interaction determines the number of connector component created in the modeling system.The interaction between model components can be in two way: (1) bi-directional such as atmosphere and ocean coupled modeling system or (2) unidirectional such as atmosphere and co-processing modeling system.In the uni-directional case, the co-processing component does not interact with the atmosphere model and only process retrieved information; thus there is one connector component.The RegESM modeling system can be configured with two different type of time-integration scheme (coupling type) between the atmosphere and ocean components: (1) explicit and (2) semi-implicit (or leap-frog) (Fig. 3).In explicit type coupling, two connector components (ATM-OCN and OCN-ATM direction) are executed at every coupling time step and model components start and stop at the same model time (Fig. 3a).However, the ocean model receives surface boundary conditions from the atmospheric model at one coupling time step ahead of the current ocean model time in semi-implicit type coupling (Fig. 3b).
The implicit and semi-implicit coupling aimed at lowering the overall computational cost of a simulation by increasing stability for longer coupling time steps.The main difference between the implicit and semi-implicit coupling type is that the models interact on different time scales in implicit coupling scenarios.
As described earlier, the execution of the model components is controlled by the driver.Both sequential and concurrent execution of the model components is allowed in the current version of the modeling system.If the model components and the driver are configured to run in sequence on the same set of PETs (Persistent Execution Threads), then the modeling system executes in a sequential mode.This mode is a much more efficient way to run the modeling system in case of limited computing resources.In the concurrent type of execution, the model components run in mutually exclusive sets of PETs, but the NUOPC connector component uses a union of available computational resources (or PETs) of interacted model components.
By this way, the modeling system can support a variety of computing systems ranging from local servers to large computing systems that could include high-speed performance networks, accelerators (i.e., Graphics Processing Unit or GPU) and parallel I/O capabilities.The main drawback of concurrent execution approach is to assign correct amount of computing resource to individual model components, which is not an easy task and might require an extensive performance benchmark of a specific configuration of the model components, to achieve best available computational performance.In this case, a load-balancing analysis of individual components and driver play a critical role in the performance of the overall modeling system.For example, LUCIA (Load-balancing Utility and Coupling Implementation Appraisal) tool can be used to collect all required information such as waiting time, the calculation time of each system components for a load-balancing analysis in the OASIS3-MCT based coupled system.
In general, the design and development of the coupled modeling systems involve a set of technical difficulties that arise due to the usage of the different computational grids in the model components.One of the most common examples is the mismatch between the land-sea masks of the model components (i.e., atmosphere and ocean models).In this case, the unaligned land-sea masks might produce artificial or unrealistic surface heat and momentum fluxes around the coastlines, narrow bays, straits and seas.The simplest solution to this issue is to modify the land-sea masks of the individual model components manually to align them.However, the main disadvantage of this solution is the required time and the difficulty to fix the land-sea masks of the different model components (especially when the horizontal grid resolution is high).Besides, the same procedure needs to be repeated whenever the model domain (i.e., shift or change in the model domain) or horizontal grid resolution is changed.As a result, this approach is considered as application specific and very time-consuming.
Unlike manual editing of the land-sea masks, customized interpolation techniques that also include extrapolation support helps to create more generic and automatized solutions.The RegESM modeling system uses extrapolation approach to overcome the mismatched land-sea mask problem for the interaction between atmosphere, ocean and wave components.To perform extrapolation, the driver uses a specialized algorithm to find the mapped and unmapped ocean grid points in the interpolation stage for every coupling direction (Fig. 4).According to the algorithm, the mapped grid points have same land-sea mask type in both model components (i.e., both are sea or land).On the other hand, the land-sea mask type does not match completely in the case of unmapped grid points.The algorithm first interpolates the field from source to destination grid using grid points just over the sea and nearest-neighbor type interpolation (from Field_A to Field_B).In this case, if the source field (Field_A) belongs to the ATM component, then the nearest source to destination method is used.In other cases, the interpolation performed using the nearest destination to source method.Similarly, the same operation is also performed by using bilinear type interpolation between the source and destination grids (from Field_A to Field_C).Then, the results of both interpolation (Field_B and Field_C) is compared to find mapped and unmapped grid points and create a new modified mask for the exchange fields (Fig. 4).
After finding mapped and unmapped grid points, the exchange field can be interpolated from source to destination grid using two-step interpolation approach.In the first step, exchange field is interpolated from source to destination grid using bilinear interpolation and the original land-sea mask.Then, result field is used to fill unmapped grid points using nearestneighbor type interpolation that is performed in the destination grid (from mapped grid points to unmapped grid points).One of the main drawbacks of this method is that the result field might include unrealistic values and sharp gradients in the areas of complex land-sea mask structure (i.e., channels, straits).The artifacts around the coastlines can be fixed by applying a light smoothing after interpolation or using more sophisticated extrapolation techniques such as the sea-over-land approach (Kara et al., 2007;Dominicis et. al., 2014), which are not included in the current version of the modeling system.Also, the usage of the mosaic grid along with second-order conservative interpolation method, which gives smoother results when the ratio between horizontal grid resolutions of the source and destination grids are high, can overcome unaligned land-sea mask problem.The next major release of ESMF library (8.0) will include the creep fill strategy (Kara et al., 2007) to fill unmapped grid points.

Integration of RegESM Modeling System with Co-processing Component
The newly designed modeling framework is a combination of the ParaView co-processing plugin -which is called Catalyst (Fabian et. al., 2011) -and ESMF library that is specially designed for coupling different ESMs to create more complex regional and global modeling systems.In conventional co-processing enabled simulation systems (single physical model component such as atmosphere along with co-processing support), the Catalyst is used to integrate ParaView visualization pipeline with the simulation code to support in-situ visualization through the use of application-specific custom adaptor code (Malakar et al., 2012;Ahrens et al., 2014;O'Leary et al., 2016;Woodring et al., 2016).A visualization pipeline integrates a data flow network in which computation is described as a collection of executable modules that are connected in a directed graph representing how data moves between modules (Moreland, 2013).There are three types of modules in a visualization pipeline: sources (file readers and synthetic data generators), filters (transforms data), and sinks (file writers and rendering module that provide images to a user interface).The adaptor code acts as a wrapper layer and transforms information coming from NUOPC cap to the co-processing component in a compatible format that is defined using ParaView, Catalyst, and VTK (Visualization Toolkit) APIs.Moreover, the adaptor code is responsible for defining the underlying computational grid and associating them with the multi-dimensional fields.After defining computational grids and fields, the ParaView processes the received data to perform co-processing to create desired products such as rendered visualizations, added value information (i.e., spatial and temporal averages, derived fields) as well as writing raw data to the disk storage (Fig. 5a).
The implemented novel approach aims to create a more generic and standardized co-processing environment designed explicitly for Earth system science (Fig. 5b).By this approach, the existing ESMs, which are coupled with ESMF library using NUOPC interface, might benefit to use integrated modeling framework to analyze the data flowing from multi-component and multi-scale modeling system without extensive code development and restructuring.In this design, the adaptor code interacts with the driver through the use of NUOPC cap and provides an abstraction layer for the co-processing component.As discussed previously, the ESMF framework uses a standardized interface (initialization, run and finalize routines) to plug new model components into existing modeling system such as RegESM in an efficient and optimized way.To that end, the new approach will benefit from the standardization of common tasks in the model components to integrate co-processing component with the existing modeling system.In this case, all the information (grids, fields, and metadata associated with them) required by ParaView, Catalyst is received from the driver, and direct interaction between other model components and the co-processing component is not allowed (Fig. 5b).The implementation logic of the adaptor code is very similar to the conventional approach (Fig. 5a).However, in this case, it uses the standardized interface of the ESMF framework and NUOPC layer to define the computational grid and associated two and three-dimensional fields of model components.The adaptor layer maps the field (i.e., ESMF_Field) and grid (i.e., ESMF_Grid) objects to their VTK equivalents through the use of VTK and co-processing APIs, which are provided by ParaView and co-processing plugin (Catalyst).Along with the usage of the new approach, the interoperability between simulation code and in-situ visualization system are enhanced and standardized.The new design also ensures easy to develop, extensible and flexible integrated modeling environment for Earth system science.
The development of the adaptor component plays an essential role in the overall design and performance of the integrated modeling environment.The adaptor code mainly includes a set of functions for the initialization (defining computational grids and associated input ports), run and finalize the co-processing environment.Similarly, the ESMF framework also uses the same approach to plug new model components into the modeling system as ESMF components.In ESMF framework, the simulation code is separated into three essential components (initialization, run and finalize) and calling interfaces are triggered by the driver to control the simulation codes (i.e., atmosphere and ocean models).In this case, the initialization phase includes definition and initialization of the exchange variables, reading input (initial and boundary conditions) and configuration files and defining the underlying computational grid (step 1 in Fig. 6).The run phase includes a time stepping loop to run the model component in a defined period and continues until simulation ends (step 4 in Fig. 6).The time interval to exchange data between model and co-processing component can be defined using coupling time step just like the interaction among other model components.According to the ESMF convention, the model and co-processing components are defined as a gridded component while the driver is a coupler component.In each coupling loop, the coupler component prepares exchange fields according to the interaction among components by applying re-gridding (except coupling with co-processing component), performing a unit conversion and common operations over the fields (i.e., rotation of wind field).
In the new version of the RegESM modeling system (1.1), the driver is extended to redistribute two and three-dimensional fields from physical model components to allow interaction with the co-processing component.In the initialization phase, the numerical grid of ESMF components is transformed into their VTK equivalents using adaptor code (step 3 in Fig. 6).In this case, ESMF_Grid object is used to create vtkStructuredGrid along with their modified parallel two-dimensional decomposition configuration, which is supported by ESMF/NUOPC grid transfer capability (Fig. 7).According to the design, each model component transfers their numerical grid representation to co-processing component at the beginning of the simulation (step 1 in Fig. 6) while assigning independent two-dimensional decomposition ratio to the retrieved grid definitions.The example configuration in Figure 7 demonstrates mapping of 2x3 decomposition ratio (in x and y-direction) of ATM component to 2x2 in COP component.Similarly, the ocean model transfers its numerical grid with 4x4 decomposition ratio to co-processing component with 2x2 (Fig. 7).In this case, ATM and OCN model components do not need to have the same geographical domain.
The only limitation is that the domain of ATM model component must cover the entire OCN model domain for an ATM-OCN coupled system to provide the surface boundary condition for OCN component.The main advantage of the generic implementation of the driver component is to assign different computational resources to the components.The computational resource with accelerator support (GPU) can be independently used by co-processing component to do rendering (i.e., iso-surface extraction, volume rendering, and texture mapping) and processing the high volume of data in an efficient and optimized way.
The initialization phase is also responsible for defining exchange fields that will be transferred among the model components and maps ESMF_Field representations as vtkMultiPieceDataSet objects in co-processing component (step 2-3 in Fig. 6).Due to the modified two-dimensional domain decomposition structure of the numerical grids of the simulation codes, the adaptor code also modifies the definition of ghost regions -a small subset of the global domain that is used to perform numerical operations around edges of the decomposition elements.In this case, the ghost regions (or halo regions in ESMF convention) are updated by using specialized calls, and after that, the simulation data are passed (as vtkMultiPieceDataSet) to the co-processing component.During the simulation, the co-processing component of the modeling system also synchronizes with the simulation code and retrieves updated data (step 5 in Fig. 6) to process and analyze the results (step 6 in Fig. 6).The interaction between driver and the adaptor continues until the simulation ends (step 4, 5 and 6 in Fig. 6) and the driver continues to redistribute the exchange fields using ESMF_FieldRedist calls.The NUOPC cap of model components also supports vertical interpolation of the three-dimensional exchange fields to height (from terrain-following coordinates of RegCM atmosphere model) or depth coordinate (from s-coordinates of ROMS ocean model) before passing information to the co-processing component (COP).
Then, finalizing routines of the model and co-processing components are called to stop the model simulations and the data analysis pipeline that destroy the defined data structure/s and free the memory (step 7-8 in Fig. 6).

Use Case and Performance Benchmark
To test the capability of the newly designed integrated modeling system that is described briefly in the previous section, the three components (atmosphere, ocean, and co-processing) configuration of RegESM 1.1 modeling system is implemented to analyze category 5 Hurricane Katrina.Hurricane Katrina was the costliest natural disaster and has been named one of the five deadliest hurricanes in the history of the United States, and the storm is currently ranked as the third most intense United States landfalling tropical cyclone.After established in the southern Florida coast as a weak category 1 storm near 22:30 UTC 25 August 2005, it strengthened to a category 5 storm by 12:00 UTC 28 August 2005 as the storm entered the central Gulf of Mexico (GoM).The model simulations are performed between 27-30 Aug. 2005, which is the most intense period of the cyclone, for three days to observe the evolution of the Hurricane Katrina and understand the importance of air-sea interaction regarding its development and predictability.The next section mainly includes detailed information of three components configuration of the modeling system as well as used computing environment, the preliminary benchmark results that are done in limited computing resource (without GPU support) and analysis of the evolution of Hurricane Katrina.

Working Environment
The model simulations and performance benchmarks are done on a cluster (SARIYER) provided by the National Center for High-Performance Computing (UHeM) in Istanbul, Turkey.The CentOS 7.2 operating system installed in compute nodes are configured with a two Intel Xeon CPU E5-2680 v4 (2.40GHz) processor (total 28 cores) and 128 GB RAM.In addition to the compute nodes, the cluster is connected to a high-performance parallel disk system (Lustre) with 349 TB storage capacity.
The performance network, which is based on Infiniband FDR (56 Gbps) is designed to give the highest performance for the communication among the servers and the disk system.Due to the lack of GPU accelerators in the entire system, the in-situ visualization integrated performance benchmarks are done with the support of software rendering provided by Mesa library.
Mesa is an open-source OpenGL implementation that supports a wide range of graphics hardware each with its back-end called a renderer.Mesa also provides several software-based renderers for use on systems without graphics hardware.In this case, ParaView is installed with Mesa support to render information without using hardware-based accelerators.

Domain and Model Configurations
The Regional Earth System Model (RegESM 1.1) is configured to couple atmosphere (ATM; RegCM) and ocean (OCN; ROMS) models with newly introduced in-situ visualization component (COP; ParaView Catalyst version 5.4.1) to analyze the evolution of Hurricane Katrina and to assess the overall performance of the modeling system.In this case, two atmospheric model domains were designed for RegCM simulations using one-way nesting approach, as shown in Fig. 8.The outer atmospheric model domain (low-resolution; LR) with a resolution of 27-km is centered at 77.5 • W, 25.0 • N and covers almost entire the United States, the western part of Atlantic Ocean and north-eastern part of Pacific Ocean for better representation of the large-scale atmospheric circulation systems.The outer domain is enlarged as much as possible to minimize the effect of the lateral boundaries of the atmospheric model in the simulation results of the inner model domain (high-resolution; HR).The horizontal grid spacing of inner domain is 3-km and covers the entire GoM and the western Atlantic Ocean to provide highresolution atmospheric forcing for coupled atmosphere-ocean model simulations and perform cloud-resolving simulations.
Unlike the outer domain, the model for the inner domain is configured to use the non-hydrostatic dynamical core (available in RegCM 4.6) to allow better representation local scale vertical acceleration and essential pressure features.
The lateral boundary condition for the outer domain is obtained from European Centre for Medium-Range Weather Forecasts (ECMWF) latest global atmospheric reanalysis (ERA-Interim project; Dee et. al., 2011), which is available at 6-h intervals at a resolution of 0.75 • x0.75 • in the horizontal and 37 pressure levels in the vertical.On the other hand, the lateral boundary condition of the HR domain is specified by the results of the LR domain.Concerning cumulus convection, Massachusetts Institute of Technology-Emanuel convective parameterization scheme (MIT-EMAN; Emanuel, 1991;Emanuel and Zivkovic-Rothman, 1999) is used in outer model simulations.Along with selected cumulus convection parameterization, sub-grid explicit moisture (SUBEX; Pal et al., 2000) scheme is used to represent large-scale precipitation for LR domain.
As it can be seen in Fig. 8, the ROMS ocean model is configured to cover entire the GoM to allow better tracking of the Hurricane Katrina.In this case, the used ocean model configuration is very similar to the configuration used by Physical Oceanography Numerical Group (PONG), Texas A&M University (TAMU), in which the original model configuration can be accessed from their THREDDS server.The ocean model has a spatial resolution of 1/36 • , which corresponds to a nonuniform resolution of around 3 km (655 x 489 grid points) with highest grid resolution in the northern part of the domain.
The model has 60 vertical sigma layer (θ s = 10.0, θ b = 2.0) to provide detailed representation of the main circulation pattern of the region and vertical tracer gradients.The bottom topography data of the GoM is constructed using the ETOPO1 dataset (Amante and Eakins, 2009), and minimum depth (h c ) is set to 400 m.The bathymetry data are also modified so that the ratio of depths of any two adjacent grids does not exceed 0.25 to enhance the stability of the model and ensure hydrostatic consistency creation that prevents pressure gradient error.The Mellor-Yamada level 2.5 turbulent closure (MY; Mellor and Yamada, 1982) is used for vertical mixing, while rotated tensors of the harmonic formulation are used for horizontal mixing.The model coupling time step between atmosphere and ocean model component is set to 1 hour but 6 minutes coupling time step is used to provide one-way interaction with co-processing component to study Hurricane Katrina in a very high temporal resolution.In the coupled model simulations, the ocean model provides SST data to the atmospheric model in the region where their numerical grids overlap.In the rest of the domain, the atmospheric model uses SST data provided by ERA-Interim dataset (prescribed SST).The results of the performance benchmark also include additional tests with smaller coupling time step such as 3 minutes for the interaction with the co-processing component.In this case, the model simulations for the analysis of Hurricane Katrina runs over three days, but only one day of simulation length is chosen in the performance benchmarks to reduce the compute time.

Performance Benchmark
A set of simulations are performed with different model configurations to assess the overall performance of the coupled modeling system by focusing overhead of the newly introduced co-processing component (Table 1).The performance benchmarks include analysis of the extra overhead provided by the co-processing component, coupling interval between physical models and co-processing component under different rendering load such as various visualization pipelines (Table 1).Two different atmospheric model configurations (a low-resolution, LR and high-resolution HR) are defined to scale up to a large number of processors.The LR model domain includes around 900.000 grid points in the atmospheric model while the HR domain contains 25 million grid points.In both cases, the ocean model configuration is the same, and it has around 19 million grid points.Besides the change of the dynamical core of atmospheric model component in HR case (non-hydrostatic), the rest of the model configurations are preserved.To isolate the overhead of the driver from the overhead of the co-processing component, first individual model components (ATM and OCN) are run in standalone mode and then, the best-scaled model configurations regarding two-dimensional decomposition configuration are used in the coupled model simulations (CPL and COP).Due to the current limitation in the integration of the co-processing component, the coupled model only supports sequential type execution (see Section 2.5 for more information) when the co-processing component is activated, but this limitation will be removed in the future version of the modeling system (RegESM 2.0).As mentioned in the previous section, the length of the simulations is kept relatively short (1 day) in the benchmark analysis to perform many simulations with different model configurations (i.e., coupling interval, visualization pipelines and domain decomposition parameters).
The benchmark results of standalone model components (ATM and OCN) can be seen in Fig. 9.In this case, two different atmospheric model configurations are considered to see the effect of the domain size and non-hydrostatic dynamical core in the benchmark results (LR and HR; Fig. 8).The results show that the model scales pretty well and it is clear that the HR case shows better scaling results than LR configuration of the atmospheric component (ATM) as expected.It is also shown that around 588 processors, which is the highest available compute resource, the communication among the processors dominate the benchmark results and even HR case does not gain further performance (Fig. 9a).Similar to the atmospheric model component, the ocean model (OCN) is also tested to find the best two-dimensional domain decomposition configuration (tiles in x and y-direction).
As it can be seen from the Fig. 9b, the selection of the tile configuration affects the overall performance of the ocean model.
In general, model scales better if tile in the x-direction is bigger than tile in the y-direction, but this is more evident in the small number of processors.The tile effect in the results is mainly due to the memory management of Fortran programming language (column-major order) as well as the total number of active grid points (not masked as land) placed in each tile.On the other hand, the tile options must be selected carefully while considering the dimension of the model domain in each direction.
In some tile configuration, it is not possible to run the model due to the used underlying numerical solver and the required minimum ghost points.To summarize, the ocean model scales well until 588 cores with the best tile configurations indicated in Fig. 9b.
The performance of the two-component modeling system (CPL) can be investigated using the benchmark results of the standalone atmosphere and ocean models.In this case, the best two-dimensional decomposition parameters of the standalone ocean model simulations are used in the coupled model simulations (Fig. 9b).The comparison of the standalone and coupled model simulations show that the driver component introduces additional 5-10% (average is 5% for LR and 6% for HR cases) overhead in the total execution time, which slight increases along with the used total number of processors, which is acceptable when increased number of MPI communication between the components are considered (Fig. 9 and 10a).The extra overhead is mainly due to the interpolation (sparse matrix multiply performed by ESMF) and extrapolation along the coastlines to match land-sea masks of the atmosphere and ocean models and fill the unmapped grid points to exchange data (Fig. 4).
To investigate the overhead introduced by the newly designed co-processing component, the three-component modeling system (COP) is tested with three different visualization pipelines (P1, P2, and P3; Table 1) using two different atmospheric model configurations (LR and HR) and coupling interval (3 and 6 minutes with co-processing).In this case, the measured total execution time during the COP benchmark results also includes vertical interpolation (performed in ESMF cap) to map data from sigma coordinates to height (or depth) coordinates for both physical model components (ATM and OCN).As shown in Fig. 10b-d, the co-processing components require 10-40% extra execution time for both LR and HR cases depending on used visualization pipeline when it is compared with CPL simulations.The results also reveal that the fastest visualization pipeline is P3 and the slowest one is P1 for the HR case (Fig. 10b and d).Table 1 also includes the execution time of the single visualization pipeline (measured by using MPI_Wtime call) isolated from the rest of the tasks.In this case, each rendering task gets 2-4 seconds for P1 and P2 cases and 7-15 seconds for the P3 case in LR atmospheric model configuration.For HR case, P1 and P2 take around 17-80 seconds, and the P3 case is rendered in around 8-10 seconds.These results show that the time spent in co-processing component (sending data to ParaView, Catalyst, and rendering to create the output) fluctuates too much and do not show predictable, stable behavior.It might be due to the particular configuration of the ParaView, which is configured to use software-based rendering to process data in CPUs and load in the used high-performance computing system (UHeM) even if the benchmark tests are repeated multiple times.
In addition to the testing modeling system with various data processing load, a benchmark with increased coupling time step is also performed (see P23M in Fig. 10c).In this case, the coupling time step between physical model components and the co-processing component is increased (from 6 minutes to 3 minutes) to produce output in doubled frame rate, but coupling interval between physical model components (ATM and OCN) are kept same (1 hour).The benchmark results show that increased coupling time step also rises overhead due to the co-processing from 45% to 60% for HR case and pipeline P2 when it is compared with the results of two-component simulations (CPL; Fig. 10a).It is also shown that the execution time of co-processing enabled coupled simulations increase but the difference between P2 and P23M cases are reduced from 66% to 37% when the number of processors increased from 140 to 588.
In addition to the analysis of timing profiles of modeling system under different rendering load, the amount of data exchanged and used in the in-situ visualization case can be compared with the amount of data that would be required for offline visualization at the same temporal frequency to reveal the added value of the newly introduced co-processing component.For this purpose, the amount of data exchanged with co-processing component is given in Table 1 for three different visualization pipelines (P1, P2, and P3).In co-processing mode, the data retrieved from model components (single time step) through the use of the driver read from memory and passed to the ParaView, Catalyst for rendering.Besides, processing data concurrently with the simulation on co-processing component, the offline visualization (post-processing) consists of the computations that are done after the model is run and requires to store numerical results in a disk environment.For example, 3-days long simulation with 6-minutes coupling interval produces around 160 GB data (720 time-step) just for a single variable from high-resolution atmosphere component (P1 visualization pipeline) in case of using offline visualization.In the case of using co-processing, the same analysis can be done by applying same visualization pipeline (P1), which requires to process only 224 MB data stored in the memory, in each coupling interval.Moreover, storing results of three-day long, high-resolution simulation of RegCM atmosphere model (in netCDF format) for offline visualization requires around 1.5 TB data in case of using 6-minutes interval in the default configuration (7 x 3d fields and 28 x 2d fields).It is evident that the usage of co-processing component reduces the amount of data stored in the disk and allows more efficient data analysis pipeline.
Besides the minor fluctuations in the benchmark results, the modeling system with co-processing component scales pretty well to the higher number of processors (or cores) without any significant performance pitfalls in the current configuration.
On the other hand, the usage of accelerator enabled ParaView configuration (i.e., using NVIDIA EGL library) and ParaView plugins with accelerator support such as NVIDIA IndeX volume rendering plugin and new VTK-m filters to process data on GPU will improve the benchmark result.The NVIDIA IndeX for ParaView Plugin enables large-scale and high-quality volume the low-resolution atmospheric model component, underlying topography information, and vorticity of ocean surface that is provided by ocean model component.

Co-processing Mode
In addition to live visualization mode that is described briefly in the previous section, ParaView Catalyst also allows to render and store data using predefined co-processing pipeline (in Python) for further analysis.Co-processing mode can be used for three purposes: (1) the simulation output can be directed to the co-processing component to render data in batch mode and write image files to the disk, (2) added value information (i.e., vorticity from wind components, eddy kinetic energy from ocean current) can be calculated and stored in a disk for further analysis and (3) storing simulation output in a higher temporal resolution to process it later (post-processing) or create a representative dataset that can be used to design visualization pipeline for co-processing or live visualization modes.In this case, the newly designed modeling system can apply multiple visualization and data processing pipelines to the simulation results at each coupling time step to make a different set of analysis over the results of same numerical simulation for more efficient data analysis.The modeling system also facilitates multiple input ports to process data flowing from multiple ESM components.In this design, input ports are defined automatically by co-processing component based on activated model components (ATM, OCN, etc.) and each model components have two ports to handle two and three-dimensional grids (and fields) separately such as atm_input2d, atm_input3d, ocn_input2d and ocn_input3d.
To test the capability of the co-processing component, the evolution of Hurricane Katrina is investigated by using two different configurations of the coupled model (COP_LR and COP_HR) that are also used to analyze the overall computational performance of the modeling system (see Section 4.3).In this case, both model configuration uses the same configuration of OCN model component, but the different horizontal resolution of the ATM model is considered (27 km for LR and 3 km for HR cases).
Figure 12 shows 3-hourly snapshots of the model simulated clouds that are generated by processing three-dimensional relative humidity field calculated by the low-resolution version of the coupled model (COP_LR) using NVIDIA IndeX volume rendering plugin as well as streamlines of Hurricane Katrina, which is calculated using three-dimensional wind field.The visualization pipeline also includes sea surface height and surface current from the ocean model component to make an integrated analysis of the model results.Figure 12a-b shows the streamlines that are produced by extracting the hurricane using ParaView Threshold filter.In this case, extracted region is used as a seed to calculate backward and forward streamlines.In Figure 12c-e, sea surface height, sea surface current and surface wind vectors (10-meters) are shown together to give insight about the interaction of ocean-related variables with the atmospheric wind.Lastly, the hurricane reaches to the land and start to disappear due to increased surface roughness and lack of energy source (Fig. 12f).While low-resolution of atmosphere model configuration is used, the information produced by the new modeling system enabled to investigate the evolution of the hurricane in a very high temporal resolution, which was impossible before.A day-long animation that is also used to create Figure 12 can be found as a supplemental video (Turuncoglu, 2018a).
In addition to the analysis of low-resolution model results to reveal the evolution of the hurricane in a very high temporal resolution, low and high-resolution model results are also compared to see the added value of the increased horizontal reso-lution of the atmospheric model component regarding representation of the hurricane and its structure.To that end, a set of visualization pipelines are designed to investigate the vertical updraft in the hurricane, simulated track, precipitation pattern, and ocean state.In this case, two time snapshots are considered: (1) 28 August 2005 0000 UTC, when it is the early stage of the hurricane in Category 5 and (2) 29 August 2005 0000 UTC when it is just before Katrina makes its third and final landfall near Louisiana-Mississippi border, where the surface wind is powerful, and surface currents had a strong onshore component (McTaggart-Cowan et al., 2007a, b).In the analysis of vertical structure, the hurricane is isolated based on the criteria of surface wind speed that exceeds 20 ms −1 and the seed (basically set of points defined as vtkPoints) for ParaView StreamTracerWithCustomSource filter are defined dynamically using ProgrammableFilter as a circular plane with a radius of 1.2 • and points distributed with 0.2 • interval in both direction (x and y) around the center of mass of the isolated region.Then, forward and backward streamlines of vorticity are computed separately to see inflow at low and mid levels and outflow at upper levels for both low (COP_LR; Fig. 13a, b, d and e) and high-resolution (COP_HR; Fig. 14a, b, d and e) cases.The analysis of simulations reveal that the vertical air movement shows higher spatial variability in high-resolution simulation (COP_HR) case even if the overall structure of the hurricane is similar in both cases.As it is expected, the strongest winds occur in a region formed as a ring around the eyewall of the hurricane, which is the lowest surface pressure occurs.In addition, the analysis of cloud liquid water content also shows that low and mid-levels of the hurricane have higher water content in a decreasing trend with height and spatial distribution of precipitation is better represented in high resolution case (Fig. 14a-b and d-e), which is consistent with the previous modelling study of Trenberth et al. (2007).It is also seen that the realistic principal and secondary precipitation bands around the eye of the hurricane are more apparent and well structured in the high-resolution simulation while low-resolution case does not show those small-scale features (Fig. 13a-b and d-e).In addition to the analysis of the atmospheric model component, the loop current, which is a warm ocean current that flows northward between Cuba and the Yucatan Peninsula and moves north into the Gulf of Mexico, loops east and south before exiting to the east through the Florida Straits and joining the Gulf Stream, is well defined by the ocean model component in both cases (Fig. 13c and f; Fig. 14c and f).The track of the hurricane is also compared with the HURDAT2 second-generation North Atlantic (NATL) hurricane database, which is the longest and most complete record of tropical cyclone (TC) activity in any of the world's oceans (Landsea and Franklin, 2013).In this case, the eye of the hurricane is extracted as a region that surface pressure anomaly is greater than 15 millibar (shown as a circular region near the best track).As it can be seen from figures, Katrina move over in the central Gulf, which is mainly associated with the loop current and persistent warm and cold eddies, and intensified as it passed over the region due to the high ocean heat content in both simulation (Fig. 13c and f and Fig. 14c and f).The comparison of the low and high-resolution simulations also indicate that the diameter of hurricane-force winds at peak intensity is bigger in high-resolution simulation case at 29 August 2005 0000 UTC (Fig. 13f and Fig. 14f).An animation that shows the comparison of low and high-resolution model results can be found as a supplemental video (Turuncoglu, 2018b).
While the main aim of this paper is to give design details of the new in-situ visualization integrated modeling system and show its capability, the performance of the coupled modeling system to represent one of the most destructive hurricanes is very satisfactory especially for high-resolution case (COP_HR).Nonetheless, the individual components (atmosphere and ocean) of the modeling system can be tuned to have better agreement with the available observations and previous studies.Specifically for the analysis of the hurricane, a better storm tracking algorithm needs to be implemented using ParaView Programmable Filter by porting existing legacy Fortran codes for more accurate storm tracking in both live and co-processing mode.
6 Discussion of The Concepts Associated with Interoperability, Portability, and Reproducibility In the current design of the RegESM modeling system, the NUOPC cap of the co-processing component is designed to work with regional modeling applications that have specific horizontal grid (or mesh) types such as rectilinear and curvilinear grids.
Therefore, the newly introduced co-processing interface (NUOPC cap and adaptor code) need to be generalized to be compatible with other regional and global modeling systems coupled with ESMF and NUOPC layer.Specifically, the following issues need to be addressed to achieve better interoperability with existing modeling systems and model components: (1) redesigning the NUOPC cap of the newly introduced co-processing component to support various global and regional mesh types such as Cubed-Sphere and unstructured Voronoi meshes, (2) extending the adaptor code to represent mesh and exchange fields provided by NUOPC cap using VTK and ParaView Catalyst APIs, (3) adding support to co-processing interface for models with online nesting capability, and (4) adding support to have common horizontal grid definition in co-processing component to make integrated analysis of data (i.e., calculating air-sea temperature difference and correlation) produced by various model components.Moreover, the co-processing interface can be tightly integrated with the NUOPC layer to provide a simplified API for designing new in-situ visualization integrated modeling systems in an efficient and standardized way.Besides, addressing issues of interoperability, supporting different modeling systems and their applications to create a more generic in-situ visualization solution, the RegESM modeling system is also tested with different model configurations such as coupling RegCM, MITgcm, and co-processing component to investigate air-sea interaction in the Black Sea basin.The initial results show that the co-processing component can also process data flowing from different model configuration supported by RegESM rather than the configuration used in this study.
When the diverse nature of high-performance computing systems, their hardware infrastructure (i.e., performance networks and storage systems) and software stacks (i.e., operating systems, compilers, libraries for inter-node communication and their different versions) are considered, realizing fully portable modeling system becoming increasingly crucial for the scientific community.In this case, the detailed examination of possible configurations of the modeling system and exiting computing environments can help to improve the flexibility and portability of the developed modeling system.Specifically for RegESM modeling system, the use case application and benchmark simulations reveal that the single executable approach (combines all model components into one program) used in the design of the modeling system can cause a portability problem when visualization and simulation run on concurrent resources.In the case of a homogeneous computing environment (all nodes with or without GPU support), the in-situ enabled modeling system runs without any particular problem because each MPI (Message Passing Interface) processes have access to the same software and hardware resources.In contrast to the homogeneous computing environment, the computing system does not have a uniform view of underlying hardware and software stack (i.e., mixed servers with and without GPU support) in a heterogeneous computing environment.As a result, the simulation with in-situ visualization support fails due to missing shared software libraries specific to underlying GPU in all used nodes or servers.In this case, two approaches can be used to overcome the problem: (1) installation of required libraries to the entire system even servers do not have GPU support, and (2) restructuring modeling system to support two executable, one for co-processing and one for physical model components.The second approach is considered a more generic and flexible solution and enhances the portability of the modeling system.It also allows implementing a loosely coupled in-situ visualization system and enables the use of specialized hardware (GPU and more memory) for rendering (Rivi et al., 2012).The main drawback of the loosely coupled in-situ visualization approach is that it requires transferring data over the network.As a result, the network performance can be a bottleneck for the modeling system, especially for high-resolution multi-component modeling applications.
When the complexity of regional and global ESMs are considered, developing fully reproducible, and portable modeling system (i.e., RegESM) is a challenging task and requires significant human interaction to keep track of detailed metadata and provenance information about the model, simulation and computing environment (in both software and hardware level).
The use of scientific workflows in earth system science has demonstrated advantages in terms of metadata, provenance, error handling, and reproducibility in an automatized and standardized way (Turuncoglu et al., 2011(Turuncoglu et al., , 2012;;Turuncoglu, 2012).
Additionally, the rapid development in the software container technology can help to design flexible and portable computing environment for earth system modeling.Hence, the Docker container has been implemented to examine the feasibility of using the container approach for newly developed in-situ visualization integrated modeling system (RegESM).A container is a standard unit of software that helps to create a package for software and all its dependencies.By this way, the application can be ported from one computing environment to another without worrying underlying hardware infrastructure and software stack.It also enhances the numerical reproducibility of simulation results by creating a standardized computing environment isolated from the dependencies.In this study, the Docker is selected as a container environment because it is widely adopted across the software industry and has a very active user community.Despite the flexibility and easy to use nature of the Docker containers, using specialized hardware such as NVIDIA GPUs, which require kernel modules and user-level libraries to operate, is not supported natively.In this case, Docker container cannot access the underlying GPU resource to perform hardware level rendering for in-situ visualization.To enable portable GPU-based containers, NVIDIA developed a special container that mounts the user mode components of the driver and the GPUs into the container at lunch.As a part of this study, the newly developed RegESM modeling system tested with both Docker (software rendering through the use of Mesa library) and NVIDIA Docker (hardware based rendering) to evaluate the potential of using container technology to improve portability and reproducibility of the RegESM modeling system.The initial results show that the developed modeling system can take advantage of the container approach to create portable and reproducible modeling system in both in-situ and co-processing modes without considerable performance loss ( 5-10%).The added value of using NVIDIA Docker is that it enables to utilize the underlying GPU resource to perform rendering (i.e., representation of clouds using direct volume rendering method).More information about the creating Docker container for in-situ visualization enabled modeling system can be found in dedicated GitHub repository (see code availability section).
In this study, the newly developed state-of-art in-situ visualization integrated modeling system (RegESM 1.1) is used to demonstrate the feasibility and added value of the integrated modeling environment to analyze the high volume of data coming from multi-component ESM in an integrated way, which was not possible before.In this case, ParaView, Catalyst is used as a coprocessing component to process and render data.The results of the selected use case (evolution of Hurricane Katrina) show that the co-processing component provides easy to use and generic modeling and data analysis environment, which is independent of the underlying physical model components used.Moreover, it promotes the usage of co-processing capability with the existing ESMs, which is coupled using ESMF framework and NUOPC layer, without significant code restructuring and development and help to increase the interoperability between ESMs and ParaView, co-processing plugin (Catalyst).In the current implementation, the prototype version of the adaptor code acts as an abstraction layer to simplify and standardize the regular tasks to integrate the simulation code with in-situ visualization and analysis environment.The driver is also responsible for redistributing the data to co-processing component while preserving its numerical grid along with the support of vertical interpolation.The coupling of co-processing component with the generic driver facilitates to define custom data processing pipelines (defined by Python scripts) easily and allows analysis of data originated from different components (i.e., atmosphere and ocean models) of the RegESM modeling system in a very high temporal resolution.In this way, RegESM modeling system can be used to study various physical processes (i.e., extreme precipitation events, air-sea interaction, convection, and turbulence) that could not be analyzed with the conventional post-processing approaches.
While the results of the in-situ visualization integrated modeling system are encouraging, the co-processing component will be extended to support different regional and global computational grid representations supported by ESMF library such as unstructured meshes for having a generic adaptor for various model applications.Additionally, we are currently exploring: (1) the way to optimize the grid transfer feature and mapping exchange fields to enhance the overall performance of the modeling environment in terms of memory usage and computational efficiency especially for very high-resolution applications (< 3 km), (2) possibility of automatic detection of accelerators (GPUs) through the use of driver component and assigning available GPU resources automatically to the co-processing component for rendering, (3) improving modeling system and co-processing component to allow nested applications (both atmosphere and ocean), (4) developing more application of the integrated modeling environment (possibly with other ocean and atmosphere components such as WRF and MITgcm) to analyze different physical processes such as air-sea interactions in upwelling regions under extreme atmospheric forcing conditions.processing tool (ParaView, Catalyst) used in the modeling system are distributed mainly by their home institutes and might apply different licensing types.The reader that might want to get more information about the individual model components and their license type could refer to the websites of them.The release version 1.1 is permanently archived on Zenodo and accessible under the digital object identifier doi:10.5281/zenodo.1307212.The demo configuration of the modeling system that is used in NVIDIA GPU Technology Conference (GTC) 2018 is also permanently archived on Zenodo and accessible under the digital object identifier doi:10.5281/zenodo.1474753.The repository also includes detailed information about the installation of the individual components of the modeling system, third-party libraries, and commands to create Docker container.
Competing interests.The author declare that he has no conflict of interest.
The lateral boundary conditions for ROMS ocean model are provided by Naval Oceanographic Office Global Navy Coastal Ocean Model (NCOM) during 27-30 August 2005.

Figure 1 .
Figure 1.Design of the RegESM coupled modeling system: (a) model components including co-processing component, (b) their interactions (orange arrows represent the redistribution and green arrows shows interpolation).

Figure 2 .
Figure2.The run sequence of model components in case of explicit type coupling.In this case, the fast coupling time step is used for the interaction between the atmosphere, ocean and wave components.The slow coupling time step is only used to interact with the river routing component.

Figure 3 .
Figure 3. Schematic representation of (a) explicit and (b) semi-implicit model coupling between two model components (atmosphere and ocean).The numbers indicate the execution orders, which is initialized in each coupling interval.

Figure 4 .Figure 6 .Figure 7 .Figure 8 .
Figure 4. Processing flow chart of the algorithm to find mapped and unmapped grid points for two-step interpolation.

Figure 9 .
Figure 9. Benchmark results of standalone (a) atmosphere (ATM; both LR and HR) and (b) ocean (OCN) models.Note that timing results of the atmosphere model is in log axes to show both LR and HR cases in the same figure.The black lines represent measured wall clock times in second and red lines (and shaded envelope) show speed-up.The best two-dimensional decomposition parameters, timing results, and speed-up are shown as a line for the ocean model case.

Figure 10 .Figure 11 .
Figure 10.Benchmark results of (a) CPL simulations (b) COP simulations with P1 visualization pipeline, (c) COP simulations with P2 visualization pipeline and (d) COP simulations with P3 visualization pipeline.CPL represents the two-component modeling system (ATM and OCN), and COP indicates three-component modeling system (ATM, OCN and co-processing).Note that the HR case requires at least 140 cores to run and the speed-up results are given based on 140 cores.38

Figure 12 .Figure 13 .
Figure 12.Rendering of multi-component (ATM-OCN-COP) fully coupled simulation using ParaView.The temporal interval for the processed data is defined as 6-minutes.

Figure 14 .
Figure 14.Same with Fig. 13 but for COP_HR simulation.The comparison of low and high-resolution model results is shown in the supplemental video.