Development of Integrated Land Simulator

Accurate simulations of land processes are crucial for many purposes, such as climate simulation, weather, flood, and drought prediction, and climate change impact assessment studies. In this paper, we present a new land simulator called the Integrated Land Simulator (ILS). The ILS consists of multiple models that represent processes related to land (hereafter, referred to as “land models”). They are coupled by a general-purpose coupler, Jcup, and executed using the Multiple Program Multiple Data approach. Currently, ILS includes a physical land surface model, the Minimal Advanced Treatments of Surface Interaction and Runoff model, and a hydrodynamic model, the Catchment-based Macro-scale Floodplain model, and the inclusion of additional land models is planned. We conducted several test simulations to evaluate the computational speed and scalability and the basic physical performance of the ILS. The results will become a benchmark for further development.


Introduction
Humans live on land, which is one of the most influential elements of the Earth system in the phenomenon of climate change. Therefore, it is important to model land processes as precisely and as detailed as possible to estimate the impact on and from the land. Land surface and hydrological modeling communities have exerted considerable effort in this direction (e.g., van den Hurk et al. 2011;Archfield et al. 2015;Decharme et al. 2019). However, progress in this direction has not been sufficient in terms of preciseness (i.e., more precisely represent geological, physical, chemical, and biological processes) or detailedness (i.e., representing land processes in higher spatiotemporal resolution). Some breakthroughs in land model development have been expected to fill the gap between societal expectations and reality (Wood et al. 2011). For example, there is a need to establish adaptation strategies for local municipalities (e.g., coastal fisheries, forestry, and water resources), and numerical experimental results of climate models are often used as reference. In such cases, the more precise and detailed simulations are for particular land regions, the more useful and beneficial the adaptation strategies are for research as well as for society as a whole.
Furthermore, there are common long-lasting biases in climate models. One such bias is the hot and dry summer over the mid-latitudinal semi-arid regions (Mueller and Seneviratne 2014;Hirota et al. 2016). The soil and snow hydrology elements of climate models seem to exhibit some inadequacies, generating surface fluxes that may trigger this bias; this issue has not yet been fully resolved. Therefore, improvements in land simulation are essential.
Climate sensitivity is one of the seven grand challenges of the World Climate Research Programme, and the land surface mechanism also impacts clouds and atmospheric circulation (Bony and Stevens 2012). Similarly, the land surface mechanism is known as a key process in other grand challenges, such as carbon feedbacks in the climate system, weather and climate extremes, and water which helps supply the food baskets of the world.
Why has land surface and hydrological model development not advanced sufficiently to help solve these issues? It is a complicated problem, as partly reported by Blöschl et al. (2019) with 23 unsolved problems in hydrology. However, particularly for the improvement of climate modeling, one technical element that prevents the smooth evolution of the land model is being one of the physical parameterizations of an atmospheric model. This affects both preciseness and detailedness, as explained below.
The preciseness of the land model is realized by better and more appropriate expressions for all types of land phenomena. Although there are independent modeling studies in hydrology, biology, geology, and other relevant fields, these studies typically require long periods for implementation in climate models or are never incorporated. This is because, historically, the land surface component in climate models was developed as one of the physical parameterizations of the atmospheric general circulation model to produce lower boundary conditions for the atmosphere and because the typical coordination of a land surface model is significantly different from those of the more precise independent models. Furthermore, climate models tend to prioritize globally impacted processes. In contrast, most land processes tend to impact the local scale; thus, they tend not to be prioritized in climate modeling.
The detailedness of the land model is even more technical than the preciseness. Although land has inherent strong heterogeneity, it cannot have a high spatial resolution because the land model is technically attached as a boundary condition of the atmospheric column, which involves a rather lower resolution. This problem essentially remains in higher-resolution atmospheric models because the heterogeneity of land is always greater than that of the atmosphere. Mosaic or tiled representation of land-type classification is a simple solution to partially mimic the behavior at a higher resolution. It should be noted that the spatial coordinates of the geophysical model are very important for expressing the target behavior. For example, ocean and atmospheric models usually use different spatial coordinates. For land models, there is also a spatial coordinate system that is more appropriate than simply using coordinates more optimized for atmospheric models, as most stand-alone land models (hydrological models, in particular) adopt basinunit coordinates.
In  (Takata et al. 2003) as the MIROC land component. It was originally developed as one of the physical parameterizations of the atmospheric general circulation model previously known as the Center for Climate System Research Atmospheric General Circulation Model, currently known as MIROC, with a vertical one-dimensional representation of the canopy, snow, and soil layers. ILS enables the improvement of preciseness in two ways. First, the MATSIRO code was rewritten in a modern way to enhance readability and userfriendliness. This makes it easier to implement new parameterizations. Second, new models can be coupled within the ILS. The extent of benefit from these improvements depend on the type of application; the former may be more beneficial for snow or soil parameterizations, while the latter may be more beneficial for independently developed vegetation models. The detailedness can be achieved by using appropriate resolution and grid systems. As the first step of this framework, we adopt the CaMa-Flood (Catchment-based Macro-scale Floodplain) hydrodynamic model (Yamazaki et al. 2011) and coupled these two models (MATSIRO and CaMa-Flood) in the ILS. The second section presents the ILS framework and components, followed by the experimental settings of the three experiments. The third section provides the results of the experiments. Finally, the last section provides the summary and conclusions.

Methods
In this section, we describe the concept of ILS, its physical processes, the input/output (I/O) component, coupling system, regridding table calculator, boundary condition creator, and benchmarking. Then, we explain three experiments to (1) examine computational speed and scalability, (2) show the impact of coupling concerning computational reproducibility, and (3) provide an example of benchmarking.
2.1 Framework and components 2.1.1 Concept of ILS Figure 1 shows a conceptual diagram of ILS. Multiple terrestrial models are coupled using a general-purpose coupler Jcup (Arakawa et al. 2020). Here, the ILS consists of two land models: a physical land surface model, MATSIRO (Takata et al. 2003;Nitta et al. 2014), and a hydrodynamic model, CaMa-Flood (Yamazaki et al. 2011). An I/O component, used to input and output files, is developed as an independent executable. The ILS employs the Multiple Program Multiple Data (MPMD) approach and executes these models and the I/O component in parallel. Each model can run in its preferred spatial coordinate and resolution. The physical land surface model, MATSIRO, is managed within the ILS because it models essential land processes. However, there are independent groups that have developed other land models. In such a case, the coupling procedure is as follows: (1) port the latest version of the component model, employ the Jcup Application Programming Interface (API) to be coupled with ILS, and (3) prepare regridding tables for spatial interpolation. If there is an updated version of that model, steps (1) and (2) are repeated. This approach allows the latest versions of each land model to be included in the ILS with minimal code modifications. We can use the same procedure for a new different land model, if necessary. As Jcup has been used to couple an atmospheric general circulation model with an oceanic general circulation model (Miyakawa et al. 2017), it is relatively easy to couple the ILS to the atmospheric and oceanic general circulation models. Previously, offline land simulations were conducted by manually modifying the code. However, ILS enables site experiments, global offline experiments, and coupled climate experiments in a single framework without any modification of the code. Our future goal is to contribute to improving climate model simulations by incorporating more sophisticated land models with rapid development and validation cycles. The details of each model of the ILS are described in the following subsections in this section.
To speed up the development and validation cycles, we decided to employ a software management environment. We manage the codes through version-control software and use GitLab (https://about.gitlab.com/) to share information via Wiki, manage issues, and test the code automatically when a merge is requested. Continuous testing is essential for model development. Currently, a 2-day simulation of the global offline experiment and a site experiment are performed, and we are planning to add benchmarking tools as well. The benchmark system is explained in a separate subsection of this section.

Physical processes
Currently, ILS includes two physical models: MATSIRO and CaMa-Flood. As explained above, MATSIRO has been adopted as the MIROC land surface component. It consists of a single-layer canopy and multiple userdefined layers of snow and soil and employs a tilescheme. MIROC participated in several phases of the Coupled Model Intercomparison Project. As one of the multi-models, many studies have evaluated its simulated land-related variables, such as surface air temperature and/or precipitation over land (Mueller and Seneviratne 2014;Miao et al. 2014), land-atmosphere interactions (Ukkola et al. 2018), evapotranspiration (Lian et al. 2018), surface albedo (Wang et al. 2016), snow cover (Thackeray et al. 2015), and permafrost (Koven et al. ). It has also been used in studies assessing water resources and flood impact, which combine offline experiments with observation-based meteorological data as inputs (Yoshimura et al. 2018). First, we rewrote the MIROC5 version of MATSIRO (Watanabe et al. 2010; hereafter referred to as MATS IRO5) to improve user-friendliness such that new schemes are easier to implement, and thus, the preciseness can be improved. We decoupled MATSIRO from MIROC5, removed all unnecessary parts produced by the decoupling, and separated the rest into MATSIRO's library and driver. The library includes only calculations of land processes, whereas the driver includes calendar management, parameter settings, reading and writing restart files, and subroutines for coupling with other land models. Users can modify the driver to conduct different types of simulations; ILS can handle any combination of land cover tiles, soil type tiles, and elevation tiles by modifying the driver. Users will choose which tiles to use and how to use them depending on their purposes. The original MATSIRO is a vertical one-dimensional model and has an IJ loop, a two-dimensional loop in space, in the innermost part of the code. We replaced this IJ loop with an expandable universal loop, which can be applied for multiple purposes, not just IJ but also ensemble members or land cover tiles. This simplifies the code, improves readability, and is advantageous for vector processors.
To confirm the reproducibility of the rewritten code, we performed two experiments with the same driving data provided by Plumbing of Land Surface Model (PALS; Best et al. 2015;Abramowitz 2012) using the ILS and MATSIRO5. The experimental settings are described in the experimental section. We examined seven main output variables: soil temperature and moisture in the first layer, snow water equivalent, canopy water, sensible heat flux, latent heat flux, and runoff for 20 sites with a 30-min temporal resolution. The results show bitidentical reproducibility in the single-point floatingpoint data type in most cases, and the differences are small enough compared to their standard deviations ( Fig. 2). We implemented iterative calculations in the flux calculation to solve the well-known instability of the skin temperature and surface fluxes. The surface heat balance equation can be written as follows: where Δ g and Δ c are the energy divergences at the ground surface and canopy, respectively; H g and H c are the respective sensible heat fluxes; R g and R c are the respective net radiations; E g , E s , E sn , E i , E is , and E t are the bare ground evaporation flux, bare ground sublimation flux, snow sublimation flux, latent heat flux from canopy interception, sublimation flux from canopy snow interception, and transpiration flux, respectively; and F g and F sn are the heat conduction fluxes to the soil and snow surface, respectively. If there is no melting of the ground surface, the above equations are solved by linearizing the skin temperatures T g and T c of the bare ground and the canopy, respectively. If the change in T g or T c from the previous time step is large, the error of the solution becomes large. In the initial verification of this study, there were some sites where the average daily variation fluctuated. We considered that one of the causes of the calculation instability might be that the skin temperature assumed an unrealistic value owing to an error in the heat balance calculation. Therefore, we implemented the Newton-Raphson method in MATSIRO.
CaMa-Flood (Yamazaki et al. 2011(Yamazaki et al. , 2013 is a global hydrodynamic model that uses daily river runoff simulated by land surface models as an input. It is currently one of the most widely used hydrodynamic models, owing to its unique method of representing flood-plain inundation using subgrid topographical parameters. In CaMa-Flood, an adaptive timestep based on the Courant-Friedrichs-Lewy condition is used. Aside from the implementation of modules necessary for the coupling system, the only modifications to the original CaMa-Flood model are related to acquiring daily river runoff data and outputting the simulated results. To input the river runoff data, Jcup is used, and only unit conversion is conducted in the modified model. Outputting the simulated results is conducted by the I/O component, rather than the original modules.

I/O component
In developing the ILS, we decided to develop an I/O component for reading and writing files independently. The I/O component is included in the MPMD program as an independent executable. Models can avoid the relatively time-consuming file input/output, such that the overall execution time can be reduced. The I/O component includes a user interface for temporal interpolations that users can implement freely. The grid interpolation function of the Jcup application programming interface (API) permits the use of input data with different spatial resolutions from that of the experiment. For example, we can use meteorological data with different spatial resolutions for an experiment. It is also possible to combine simulated results for only land grids with missing values for ocean grids before outputting the results. The function allows single-process execution and two options with multiple-process (MP) executions. The first MP option uses the Message Passing Interface (MPI) I/O function. In the second option, the root processor collects data from other processors and outputs them. The I/O component supports the binary and Network Common Data Form (NetCDF) formats. It is designed to be extendable to other output formats such as gtool. For NetCDF, metadata are prepared in JavaScript Object Notation (JSON) format. File names and metadata follow the Assistance for Land Surface Modeling Activities (ALMA) and Climate and Forecast (CF) conventions.

Coupling system
As described above, the ILS is designed to perform simulations by coupling a plurality of components. This section describes the coupling system built into the ILS. The coupling software used in ILS is Jcup (Arakawa et al. 2020). Jcup can couple multiple components combining both parallel and serial implementations. In addition, it has high versatility because there is no restriction on applicable grid coordinates, and users can freely implement interpolation schemes. However, because of its high flexibility, the interface is complicated, making it difficult for beginners to use. Therefore, an API with a more compact interface customized to the ILS was developed. The most distinguishable feature of Jcup is that it requires a correspondence between the sending-side grid point index, receiving-side grid point index, and interpolation coefficients in the interpolation calculation as input information (hereafter, referred to as "regridding table"). Figure 3 shows a typical use of the Jcup API subroutines. Before the time-integration loop, there are three types of subroutines (and moj_put_data), such as the initialization of the coupler, setting of grid information, and setting of the regridding table. In the time-integration loop, moj_set_time is called near the beginning of the loop. Data exchange and interpolation calculations were performed inside this subroutine. Subroutine moj_get_data obtains the data from the sending component, and moj_put_data is a subroutine that transfers the data to be sent to the receiving component to the coupler. These two subroutines can be called anywhere in the time-integration loop.
As described above, Jcup and its API require a regridding table as input information; therefore, the user needs to calculate a regridding table in advance. In previous cases, the calculation tool was developed individually for the models to be coupled. For example, to couple the icosahedral atmospheric model NICAM (Satoh et al. 2014) and the global ocean model COCO (Hasumi 2015), a program that searches for grid points included in each polygon of the opposite model grid and calculates the area of the overlapped region has been developed . For the seismic model Seism3D (Furumura 2005) and the structure model FrontISTR++ (Okuda 2019) coupling, the regridding table could be calculated easily because the grid structure of the seismic model was a simple rectangular grid with a constant grid spacing (Matsumoto et al. 2015).
Unlike these cases, the ILS has the possibility of coupling various models that have their own grid coordinates, and the combination of coupling should be flexible; therefore, a regridding table calculation tool applicable to a wide range of grid coordinates is required. Hence, a general-purpose regridding table calculation tool called Spheroidal Coordinates Regridding Interpolation Table Generator (SPRING) has been developed. The next section details this software.

Regridding table calculator
SPRING (Takeshima 2020) has been developed to couple models with any type of grid coordinate system, such as the latitude-longitude grid, polygonal grid, or grid defined by an aggregation of tiny latitude-longitude cells (Fig. 4). The basin-shaped coordinate system of the CaMa-Flood is of the last type. This tool provides regridding tables, a dataset of overlapping grid indices, and their interpolation weights for these coordinates. In conservative data communication in the ILS, the data value X IJ in the receiving-side grid point IJ is calculated by the following equation: where ij is the sending-side grid point index that overlaps with the grid IJ, c ij is the interpolation weight, x ij is the data value, s 0 ij is the overlapping area between grid IJ and grid ij, and S IJ is the area of the grid IJ. In the default configuration, MATSIRO and CaMa-Flood have basin-shaped coordinates, whose grids consist of tiny latitude-longitude cells, and the forcing data used in this study are expressed on a latitude-longitude grid. Therefore, s 0 ij and S IJ are calculated by the following equation regarding Earth as a perfect sphere with radius a where S is the area of a rectangular shape, θ 1 and θ 2 are the longitudes of the western/eastern boundaries, respectively, and φ 1 and φ 2 are the latitudes of the northern/southern boundaries, respectively. The detailed evaluation of regridding tables will be presented in a separate paper and is beyond the scope of this study. Fig. 4 Examples of grids to which SPRING is applicable. a Latitude-longitude grid, b NICAM grid (Satoh et al. 2008), a type of convex polygon grid, c concave polygon grid, and d basin-shaped grid, a type of grid defined by the aggregation of tiny cells

Boundary condition creator
The ILS includes a tool to create the boundary conditions of MATSIRO. They were manually created for different experiments. However, the information was not sufficiently shared, and sometimes the data included human errors. To avoid these issues, we manage the version of the codes to create boundary conditions. Thus, we can use the common parts repeatedly. Figure 5 shows a conceptual diagram of the boundary creation tool. By combining the original datasets and regridding tables described in the previous section, we can create new boundary condition data for any horizontal resolution, which include land cover type, soil type, soil albedo, standard deviation and gradient of elevation, and leaf area index (LAI). First, the land cover categories in the high-resolution original data are reclassified into categories to be used in MATSIRO. The regridding table is then used to calculate the dominant categories. For soil types, the proportions of sand, silt, and clay are interpolated using a regridding table. The soil type is then determined using the triangle diagram of the soil properties. The soil albedo is simply interpolated. The standard deviation and slope of elevation are calculated from high-resolution topographic data. The LAI is derived by a simple interpolation without the tiling scheme. If a land cover tiling scheme is used, the LAI is first classified by land use data and then averaged. The boundary tool creator allows the implementation of a new method, as the method varies depending on the source data.

Benchmarking
ILS incorporates benchmarking as part of the model development cycle. It is important to validate the model against a variety of terrestrial observations during each model development cycle and to evaluate how new parameterizations or models for preciseness and changes in resolution for detailedness can improve the simulated results. In this paper, we present the results of an experiment using PALS, which includes 20 Fluxnet sites (Fig. 6).
The experiment is suitable for benchmarking because the sites cover the locations with various conditions, and its computational cost is low. The experimental settings are shown in the next section.

Specification of three experiments
First, we conducted global offline land simulations and evaluated the computational performance. The ILS was forced using three-hourly meteorological datasets prepared by the Global Soil Wetness Project Phase 3 (GSWP3; Kim et al. 2019). GSWP3 forcing data is based on downscaling of the Twentieth-Century Reanalysis (Compo et al. 2011), including rainfall, snowfall, temperature, wind speed, specific humidity, shortwave and longwave radiation, and surface pressure. Cloud cover from the same downscaling experiment was also  Figure 7 shows the computational flow of the global offline experiment. The I/O component reads the meteorological forcing and boundary data, interpolates them temporally, and sends them to MATSIRO. Then, MATSIRO calculates land processes and sends simulated runoff to the CaMa-Flood routine, where the river processes are calculated. The output data from MATS IRO and CaMa-Flood are sent to the I/O component to be written to files. The experiment was run on a computer cluster with Intel Xeon processors at the University of Tokyo. We changed the number of cores for a component from 1 to 128 while fixing the number of cores (16 cores) for the other two components.
Second, we examined the difference between uncoupled (stand-alone) and coupled CaMa-Flood simulations to check whether coupling, including the regridding tables created by SPRING, works as expected. Note that we did not consider feedback from CaMa-Flood to MATSIRO in the coupled simulation and take the uncoupled simulation as the truth. We conducted a stand-alone river discharge simulation using the original CaMa-Flood with runoff data from MATSIRO in ILS (uncoupled simulation). We compared it with the simulation mentioned above (coupled simulation). This experiment intends to confirm if the coupled simulation is as similar to the uncoupled simulation by CaMa-Flood, in which its original regridding scheme is implemented.
Finally, we conducted site experiments for 20 Fluxnet sites from PALS, as shown in Fig. 6, to (1) validate the rewritten code and (2) show the results of the first benchmark. The sites include tropical, arid, temperate, and cold regions, as assigned by the Koppen climate classification (Peel et al. 2007), and forest and non-forest regions. The dataset includes meteorological forcing data: downward shortwave and longwave radiation, temperature, precipitation, specific humidity, wind speed, and surface pressure. In addition to these data, cloud cover from ERA5 (C3S 2017) was used. PALS also includes surface flux data from the observations and estimated by one-, two-, and three-variable regression equations for validation and benchmarking. The time steps of the input data and simulation were 30 min without any missing data. The simulation periods range from 2 to 10 years, depending on the site. We used the MIROC5 version of the ILS for the verification shown above. Then, we implemented the iteration of the flux calculation to MATS IRO as described above and evaluated the PALS experiment. Further observation stations, validation of global experiments, and increasing the amount of observation data are clearly in our future work. In this section, we show the results from the three experiments explained in the previous section. Then, we examined the total time spent on each component of MATSIRO: initial setup, main MATSIRO calculation, communication by the coupler, receiving the input data from the coupler, and sending the output data to the coupler. In this analysis, we only used the experiments with fixed 16 cores for the CaMa-Flood and I/ O components, respectively, and increasing the number of cores for MATSIRO from 1 to 128 cores (Fig. 8b). The initial setup with 64 cores took longer than with 32 cores; the initial setup with multiple nodes seems to take more time. The coupler communication cost is relatively fixed. In all settings, the MATSIRO main calculation is the most expensive component; it accounts for approximately 95% of the MATSIRO computation with 16 cores. However, it shows close to perfect scalability such that we would expect ILS to become faster with more cores. Although the computational performance depends on machines, we believe that high-resolution computation with a high-performance computer is promising. The validation of the global offline experiment is beyond the scope of this study.

Impact of coupling with respect to computational reproducibility
The impact of coupling CaMa-Flood to ILS with respect to computational reproducibility is displayed in Fig. 9. Figure 9a shows the standard deviation of river discharge, flood-plain discharge, river storage, flood-plain storage, and flooded area fraction from stand-alone CaMa-Flood simulation normalized by the mean for each 10-degree band from 60°S to 90°N. The maximum value is approximately 11, found in flood-plain discharge at 30-40°N. The minimum value is approximately 1, found in the flooded area fraction at 60-50°S. Figure 9b is the same as Fig. 9a, but for the difference between coupled and uncoupled simulations, also normalized by the mean. The maximum difference is found in floodplain discharge at 50-60°N, which is less than 0.001% of the standard deviation. The differences were due to regridding tables, I/O, and simulation framework and are negligible compared to the standard deviation.    Figure 10 summarizes the performance of the halfhourly latent and sensible heat fluxes using a Taylor plot. Blue indicates ILS simulations, and each plot shows each site. For reference, the results of linear regressions of the downward shortwave radiation and surface air temperature against surface fluxes from the PALS project are shown in gray. The statistics were normalized by the standard deviation of the reference data. Figure 10a shows the results of the latent heat flux. The median correlation coefficients were 0.78 and 0.83, and the median root-mean-square differences (RMSDs) were 0.66 and 0.59 for the ILS and empirical methods, respectively. For these two metrics, the empirical method showed better reproducibility. However, the median standard deviations were 0.93 and 0.74 for the ILS and the empirical methods, respectively. ILS results are closer to unity, which is the normalized standard deviation of the reference data. The empirical approach underestimated the standard deviation in most places. Figure 10b shows the results of the sensible heat flux. In the ILS and empirical methods, the median correlation coefficients were 0.84 and 0.89, the median RMSDs were 0.60 and 0.50, and the median standard deviations were 1.14 and 0.90, respectively. The empirical method tended to be better for all three metrics.

First result of benchmarking
Comparing the results of latent and sensible heat fluxes, the reproducibility of sensible heat fluxes tends to be better than that of latent heat fluxes. The result is more scattered in the sensible heat flux than the latent heat flux. The detailed analysis is beyond the scope of this study and will be shown in a separate validation paper.
The empirical method had a problem with the reproducibility of extreme values. For example, when the 99.9th percentile values of the latent heat and sensible heat flux at each site were examined for 30 min, the ILS tended to produce better results (Fig. 11). The mean absolute errors of latent heat flux were 77.1 W/m 2 and 203.7 W/m 2 for the ILS and empirical equation, respectively, and the mean absolute errors of sensible heat flux were 76.8 W/m 2 and 149.9 W/m 2 for the ILS and empirical equation, respectively. The ILS also showed a better tendency for the 99th and 99.99th percentile values.

Summary and conclusions
The improvement of the land surface model is expected not only by the climate community but also by society; however, there are two characteristics that are necessary for such a feat-preciseness and detailedness. Preciseness can be achieved by including geological, physical, chemical, and biological processes and their interactions into the comprehensive model. Detailedness can be achieved by computing a land model simulation with high spatiotemporal resolution, considering the inherent heterogeneous conditions of the land. A barrier to realizing these two targets is that the land model was developed as one of the physical parameterizations of the atmospheric model. Therefore, in this study, we developed a new framework to enhance the speed of the improvement of our land surface model. We aim to develop a new framework to facilitate the realization of improved preciseness and detailedness in a land model simulation, a framework named ILS. Note that this paper only addresses a way to construct the framework. Further application of the framework, for example, adding new land models to the ILS for better preciseness or executing the ILS in a much higher spatiotemporal resolution using heterogeneous boundary conditions, will be conducted in the future. First, we removed the land surface model, MATSIRO, from the climate model, MIROC. The vertical onedimensional core of MATSIRO is used as a component in the ILS. Then, the code was verified. The new hydrodynamic model, CaMa-Flood, was implemented in ILS. One of the advantages of the ILS is to port a new standalone model almost as it is, using the new generalpurpose coupler, Jcup. Furthermore, we developed ancillary but important tools, such as a component to handle input and output files, modules to create necessary boundary conditions automatically, and modules to generate regridding tables to interpolate the geospatial information. Next, we conducted some numerical experiments to confirm two aspects of performance evaluation-computational efficiency and validity. Computational efficiency was evaluated by conducting global simulations using different combinations of CPUs and ILS components, that is, MATSIRO, CaMa-Flood, and I/ O, and good scalability was demonstrated for simulations at very high resolutions with a many-core architecture. The validation of the land model by comparing it with reality is still preliminary; however, for selected sites used in the PALS project, the ILS was satisfactory.
As mentioned above, this paper presents the first description of the new framework. Applications of the ILS have already been planned in emerging fields, such as human activities (irrigation/pumping/damming; Hanasaki et al. 2008), sediment dynamics (Hatono and Yoshimura 2020), groundwater dynamics (Miura and Yoshimura 2020), water temperature and water quality in rivers (Tokuda et al. 2019), vegetation dynamics, subgrid representability of hillslope processes, and isotopic representations (Yoshimura et al. 2006). We believe that the basic development conducted in this study is needed to further advance land/climate models and improve the understanding and prediction of the earth system.