Interactive comment on “ A Lagrangian Advection scheme with Shape Matrix ( LASM ) for solving advection problems ” by L .

Dear authors, In my role as Executive editor of GMD, I would like to bring to your attention our Editorial: http://www.geoscientific-model-development.net/gmd_journal_white_paper.pdf http://www.geosci-model-dev.net/6/1233/2013/gmd-6-1233-2013.html This highlights some requirements of papers published in GMD. In particular, please note in the abstract of the Editorial: “– The paper must be accompanied by the code, or means of accessing the code, for C1375


Introduction
Advection is one of the key problems in geophysical modeling research.Many tracers (or constituents) need to be advected by the flow simulated by the atmospheric or oceanic dynamical core and physical parameterization.In an atmospheric global circulation model (AGCM), the most important tracers are the water vapor and other phases of the water (e.g., cloud water and cloud ice), which participate in many parameterization processes (e.g., convection, microphysics, and radiation).The quality of the computed water substance distributions are vital to the successful simulation, although many other aspects also affect the results.
The numerical schemes can be divided into three categories: (1) Eulerian schemes, (2) semi-Lagrangian schemes, and (3) Lagrangian schemes.The Eulerian schemes (Yu, 1994;Lin and Rood, 1996;Putman and Lin, 2007) are constructed on the fixed meshes and compute the flux of mass through the Eulerian cell edges.The semi-Lagrangian schemes (Staniforth and Côté, 1991;Lauritzen et al., 2010) are also bound to the fixed meshes, but they track the mass in a cell and restart from grids each time step.Lauritzen et al. (2011) gave a thorough review on the atmospheric transport schemes emphasized on the semi-Lagrangian view on finite-volume discretizations.In comparison, the Lagrangian schemes use mobile parcels as the discretization units and track them all the time.The most significant advantage of the Lagrangian schemes over the first two is that the numerical diffusion is highly suppressed because the simulated tracer parcels adapt to the flow and do not restart from the fixed grids.This advantage is praised in Stenke et al. (2008), which used a Lagrangian scheme to advect water vapor and cloud water, and gained positive results compared with the semi-Lagrangian scheme in ECHAM4.
Besides the application in the AGCMs, the Lagrangian schemes are also suitable for the use in the atmospheric chemistry transport models (CTMs); see Brunner (2013) for an overview.With the increasing of the chemical tracer number (approaching hundreds), the multi-tracer efficiency becomes an important issue.Some Eulerian or semi-Lagrangian schemes need to advect each tracer species independently, so this low efficiency may become a significant Published by Copernicus Publications on behalf of the European Geosciences Union.
L. Dong et al.: Lagrangian advection scheme with shape matrix computational bottleneck.In Lagrangian schemes, however, the trajectory information and some weights are shared by all the species, so this efficiency is very high.In addition, the trajectory information of individual parcel provided in default is useful for addressing sources and transport pathways.For the nonlinear chemical reaction, the unphysical numerical diffusion is a problem which necessitates more efforts to control in the Eulerian and semi-Lagrangian schemes when the resolution is coarse, and it will disturb the reaction.On the other hand, when the resolution is refined, the time step size must be reduced substantially due to the Courant-Friedrichs-Lewy constraint for some Eulerian schemes.Contrarily, the Lagrangian schemes have no numerical diffusion, and the interparcel mixing needs to be implemented explicitly, but this mixing can be designed in a physical way, which is good for simulating chemical reaction realistically.
Several Lagrangian schemes have also been proposed in the last decade: (1) ATTILA (atmospheric tracer transport in a Lagrangian model, Reithmeier and Sausen, 2002); (2) CLaMS (Chemical Lagrangian Model of the Stratosphere; McKenna et al., 2002); (3) TTS-C/I (Trajectorytracking scheme -"C" and "I" mean the tracking object is the centroid and interface, respectively; Dong andWang, 2012, 2013); and (4) HEL (Hybrid Eulerian-Lagrangian; Kaas et al., 2013).The discretization unit in the Lagrangian schemes is the tracer parcel, which carries tracer information and is advected by the flow passively all the way without restarting from the fixed model grids each time step as in the semi-Lagrangian schemes.To be applied in an AGCM, the carried tracer information must be remapped onto the the model grids, and some other information like tendencies from the physical parameterization needs to be remapped back onto the parcels.In some Lagrangian schemes unrealistic "spotty" tracer distribution will occur when the flow deforms greatly as noted in TTS-C (see Fig. 13 in Dong and Wang, 2012) and HEL (see Sect. 1.3 in Kaas et al., 2013).McKenna et al. (2002) also noted that small (unmixed) structures will arise (see Fig. 8 there).This is caused by the fact that the parcel spatial extent or shape is not simulated explicitly, so the remapping is isotropic, e.g., the inverse distance weighting interpolation in TTS-C.When the parcel is deformed severely, the tracer mass will not be distributed to the right grids or resolved well, and the so-called spectral blocking or aliasing will occur (Kaas et al., 2013).The main motivation of this work is to address the problem of tracer mass remapping between the irregularly distributed parcels and the fixed model grids (e.g., regular latitude-longitude grids).
ATTILA, CLaMS, and HEL all introduced some kinds of physical mixing among parcels to solve the above unrealistic tracer distribution problem with different methodologies.In order to keep the mass of air associated with a parcel as a compact volume around the parcel centroid, ATTILA redefines the parcel boundaries by bringing the mass mixing ratio c of a species in a parcel closer to an average back-ground mixing ratio, or called "interaction by exchange with mean" (Henne et al., 2013).CLaMS uses a Lyapunov exponent to measure for the deformation in the flow, and inserts or merges parcels based on this exponent.HEL incorporates the deformation rate of the flow and mixes tracers in a directionally biased way.In 2-D cases, two passive auxiliary Lagrangian parcels are associated with the main parcel to identify the asymptotic dilatation axis caused by the shear.Only model grids close to the axis are assigned the tracer mass carried by the main parcel, and this mixing is mainly along the asymptotic dilatation axis.
TTS-I tried another path by computing the parcel interfaces explicitly, which are polygons in 2-D cases, so the shape of each parcel is simulated.The polygons are advected and deformed by the flow, and new vertices are inserted to capture the curvature changes of the real interfaces.This solved the remapping problem to a great extent, and no mixing concept was added at this pure advection stage, though physical mixing is still needed in the real application.The drawback of this approach is that the complex geometric calculation, especially on the sphere, hindered its application.Furthermore, in 3-D cases the parcel polyhedron is almost impossible to deal with.
In this work a new design of TTS is proposed with the consideration of pros and cons in TTS-C and TTS-I, and other schemes.To extend the new scheme to 3-D easily, the centroid discretization is used as in TTS-C and other schemes.The shape of a parcel is still computed explicitly as in TTS-I to avoid the aliasing error as much as possible, but it is represented by a linear transformation matrix as in Yserentant (1999) and Gauger et al. (2000) with some key modifications that bypass the costly geometric calculations in TTS-I.This transformation transforms the physical space into a remapping space after translation, rotation, and size change of the parcel.In this manner the flow deformation can play an important role in the remapping process.Several computation difficulties on the sphere have been addressed.Henne et al. (2013) stated that this type of parcel shape treatments allow for a more realistic description of the interparticle mixing and also enable the realization of a fully Lagrangian dynamical core, but these interesting and important topics are not covered in this work.A flow adaptive interparcel mixing algorithm is also constructed in the new scheme to ensure the approximation of the parcel shape valid in the nonlinearly deformational flow.In other words, when the flow deformation is comparatively linear, the mixing can be turned off.The new scheme is renamed LASM, which stands for Lagrangian advection scheme with shape matrix, to better represent its features.
The paper is organized as following: Sect. 2 presents the new LASM in detail, i.e., the updating of the deformation matrix, and the interparcel mixing.Section 3 gives the results from several test cases to show the performance of LASM.The conclusion is drawn and discussion presented in Sect. 4. The code availability is shown in Sect. 5.

LASM details
The advection equation in the Eulerian framework can be written as where V is the velocity of the flow and ρ m is the tracer density for the tracer species m.The mixing ratio of tracer (m > 0) relative to the background tracer or dry air (m = 0) is and will be used in the following tests.In the Lagrangian framework, this partial differential equation turns into a differential equation system as where x is the coordinate of a parcel centroid.So the continuous tracer field is discretized into tracer parcels with a finite number accordingly.The parcel carries the tracer densities and masses for each species, and is transported by the given flow.Equation (3) is the trajectory equation for the parcel centroid.Its computation on the sphere is elaborated in Dong and Wang (2012), which integrates Eq. (3) by using the fourth-order Runge-Kutta method numerically and projects the centroids onto the polar stereographic plane when they are near the poles.Equation (4) describes the evolution of tracer density, which is affected by the flow divergence.The divergence on the model grids is evaluated by the secondorder central finite difference in this study.The interpolation of velocity and divergence onto the parcel centroids is bilinear.
Each parcel is constructed in the following two steps.Firstly the shape of the parcel is expressed by a linear trans-formation (see Fig. 1a) as where x 0i is the centroid coordinate of parcel i, x is some spatial coordinate, y is the body coordinate in the remapping space with x 0i as the origin, and H i is the transformation matrix or deformation matrix.In 2-D cases the shape of the parcel is an ellipse, and in 3-D it is an ellipsoid.The parcels can overlap each other.The inverse transformation is The matrix H i is the key of this transformation, which controls the shape and size of the parcel.This idea is proposed by Harry Yserentant through a series of papers (Yserentant, 1997a(Yserentant, , b, 1999(Yserentant, , 2003;;Gauger et al., 2000) in the context of fully Lagrangian dynamics.Compared with the polygons in TTS-I which can not overlap each other, the linear transformation is more feasible in a real application.Secondly the mass distribution within the remapping space of a parcel is prescribed by a shape function or kernel function: which is a composite of a B-spline function (see Fig. 1b) as in Gauger et al. (2000): where d is the spatial dimension number and y k is the kth component of the body coordinate y.This form of a shape function is adopted in this work, since there is no driving need to replace it.
L. Dong et al.: Lagrangian advection scheme with shape matrix Figure 2. Four skeleton points (black points) associated with a parcel.Their spatial coordinates will change with the flow, but the body coordinates in remapping space are fixed.The red point is the major axis vertex of the parcel.

Update of deformation matrix
The core of the new LASM is the deformation matrix H i of parcel i, which is updated every time step according to the advection flow as in the following equation (Yserentant, 2003): where ∇V is the velocity gradient tension (a d × d matrix).But on the sphere, the solution of Eq. ( 9) meets some difficulties when parcels traverse the poles in the spherical coordinate system.In some spherical meshes like cubed-sphere mesh, the coordinate system is on the local (curvilinear) coordinate, so the difficulties may be bypassed.In spite of this, LASM needs to support the lat-lon mesh, so other solutions are sought.To avoid the pole problems, the transformation is conducted on the local stereographic projection plane with the parcel centroid as the origin.The spatial coordinate is labeled as x on that plane, and Eqs. ( 5) and ( 6) become The coordinate transformation between spherical coordinate and local stereographic projected coordinate is routine and omitted here.
In 2-D cases, four auxiliary skeleton points are associated with each parcel.Their body coordinates are set as (−1, 0), (0, −1), (1, 0), and (0, 1).Initially the parcels are the shape of a circle (see the left figure in Fig. 2), and the distance between centroid and skeleton point is chosen as 1.5 times the maximum local model grid interval so that all the grids are covered by some parcels.It is noteworthy that a larger distance will cause a smoother initial tracer density distribution when remapping between tracers and grids as in Sect.2.2.
The skeleton points will be advected along with the parcel centroid, and the relative locations of them will be changed if the flow is deformational (see Fig. 2).So we can say that the skeleton points sense the flow deformation, and they determine H i .From Eq. ( 10), the elements of H i can be calculated as where the superscripts ( * ) are the labels for the skeleton points, and h (1) , h (1) , h (2) The determinant of H i is where 22 , and similar for the other three matrices.It is obvious that H i is composed of four matrices: , and The matrix H i can be decomposed by using the SVD (singular value decomposition) technique, which provides useful information about the transformation for the subsequent computation: where S i is a positive diagonal d × d scaling matrix with the diagonal elements in descending order, and U i and V i are for rotation.The production of the diagonal elements of S i is the determinant of H i , which is also a representation of the parcel volume V i .On the other hand, after solving the density of the tracer from Eq. ( 4), we can also calculate the tracer volume by where m i and ρ i are the mass and density of one tracer species (the first species is used, since it should be no difference for other species).Generally there will be discrepancy between the determinant of H i and V i , so the determinant of H i is reset to V i by scaling S i each time step: where s k i is kth diagonal element of S i .Then H i is reconstructed by replacing S i in Eq. ( 13) so that the rotation part (i.e., U i and V i ) is kept.After updating H i , the shape of the parcel, which affects the tracer mass remapping, on the new time step under the deformation of the flow is obtained.The remapping will be anisotropic so that the flow deformation is respected, and this is the key: that LASM can reduce the aliasing error without using substantial mixing as in other Lagrangian schemes.Section 2.3 will discuss the interparcel mixing in more detail.Compared with TTS-I, the deformation matrix is more practical than the polygon with varying edge number.Additionally LASM is also 3-D ready, and no quasi-Lagrangian vertical coordinate is needed when applied to a real AGCM or oceanic general circulation model (OGCM).

Remapping between tracers and grids
Since other dynamical and physical processes of the model take place on the Eulerian grids, the tracer densities carried by the parcels need to be remapped onto the grids.
After the trajectory calculation and deformation matrix updating, the parcels will be connected with the grids that have non-zero shape function value.This involves the neighbor searching algorithm.The searching procedure is as follows: -the grids covered by the circle (or sphere in 3-D) with the parcel semi-major axis length as their radius are searched for; -the grid coordinates are transformed into the parcel's body coordinate system by Eq. ( 11); -the corresponding shape function values are calculated by Eq. ( 7); -the grids with non-zero shape function value connect with the parcel.
Figure 3 shows a result of the neighbor searching.There is a flaw with this procedure.Imagine a case when a grid (i.e., cell center) is not covered by a parcel, but some part of the cell is; then the grid does not get tracer mass from that parcel.This grid is called a partially covered grid.The effect of this flaw becomes visible when the parcel is stretched into a filament.Therefore an adaptive interparcel mixing algorithm is proposed in the next subsection to reduce this effect.
In the future, another improved searching procedure may be designed which can process the partially covered grids, and the shape function value is modified to a reasonable non-zero value.
The remapping density for the mth species on the grid I is where the subscripts i and I are for parcel and grid, respectively, and S I is the set of parcels that are connected with this grid.Since this is a linear remapping, the density on the grid will be in the range of the minimum and maximum density on the parcels.The remapping in the reverse direction is similar to Eq. ( 16) as where S i is the set of grids that are connected with parcel i.The density (not the mass) is chosen as the remapping quantity, because the remapping weight is based on "distance", not overlapping area between parcels and grid cells as in TTS-I.So when the areas of the model grid cells are highly nonuniform (e.g., lat-lon mesh), large density error on the grids will occur if the mass is remapped, especially near the poles on the lat-lon mesh.The grid density will be constrained directly through Eq. ( 16) when the density is remapped, but the drawback is that the total tracer mass on the grids can not be guaranteed as a constant without sources and sinks.This dilemma is also encountered by TTS-C and HEL.A mass correction must be done if the total mass on the grids must be conserved.Diamantakis and Flemming (2014) described several mass fixer algorithms for the semi-Lagrangian schemes, which can be used in LASM.In this study, the simplest fixer is used as in TTS-C and HEL.This may make people uncomfortable, but, as stated in Kaas et al. (2013), it is different from the traditional global methods in Eulerian and semi-Lagrangian schemes because the correct values are preserved in the Lagrangian parcels; in other words, the total mass on the parcel is exactly conserved, and therefore the total mass on the grids is well constrained.The profound impact on the model simulation should be analyzed in the real application.On the other hand, when more uniform mesh is used (e.g., cubed-sphere, icosahedral mesh), we can remap the mass and the total mass on the grids is conserved automatically.There may be occasions during which some grids are not connected with any parcel, which are called void grids as in ATTILA.When this happens, the void grids will be filled by the density value interpolated from its neighbor grids by using the inverse-distance weight.In the real application, due to the continuity condition, the void grids are conjectured to be very rare.
Finally because the remapping weight ψ i (I ) is the same for all the tracer species, the multi-tracer efficiency will be as good as HEL, which meets the needs of chemistry transport models with hundreds of tracer species.

Interparcel mixing
It would be ideal for all the parcels to move around and keep good shape, but in the deformational flows the large shear will elongate some parcels along the asymptotic dilation axis and the strong vortex will stir the parcels heavily.At some critical time points, on the one hand the determinant of H i may turn out to be negative if the parcel size is not limited, which is unphysical and means negative volume; on the other hand the linear deformation approximation of the parcel shape may be violated.These are caused by the finite par-cel resolution.To avoid problems, Klingler et al. (2007) proposed a restart procedure in the finite mass method, which replaces the old, deformed parcels (particles in their case) by more regularly distributed parcels (i.e., redistribution).But the global parcel redistribution will eventually cause unacceptable numerical diffusion, and some good features of the Lagrangian method will be lost (e.g., good preservation of tracer correlation).Contrarily a different path is taken in LASM, where a parcel i is mixed with its surrounding parcels (the connected parcels of the cell where parcel i is contained) when the ratio γ i of the major axis of the parcel shape to the minor one (i.e., s 1i /s 2i in 2-D) exceeds any given threshold γ m (normally 100 in this study).
The threshold γ m is adjusted by the flow.We define a disorder degree D i of parcel i by the angles between the major axes of it and its surrounding parcels.Firstly the angles are calculated on the local stereographic plane of parcel i. Secondly the ratio between the maximum and mean angles1 is defined as D i .When D i is small, the parcels are well organized so that they are free to sense the deformation of the flow, such as the case in Fig. 4a.When D i exceeds any given threshold D * (e.g., 1.05), the parcels are disordered and the flow is more nonlinearly deformational, such as in Fig. 4b, so the interparcel mixing needs to be more easily triggered to improve the linear deformation approximation by reducing γ m to a strict value γ * m (e.g., 5) for the involved parcels.It is noteworthy that the mixing in LASM is adapted to the flow deformation, not the tracer density gradient as in the Eulerian or semi-Lagrangian schemes, so the correlation between tracers will not be largely disturbed because the mixing degrees of all the tracers are the same.It should also be noted that the values and formulas of the above parameters are by no means the optimal ones.The current values are gained under the objective of making LASM perform well in both ideal and barotropic test cases.For example, when γ * m is larger, unacceptable results (more aliasing) will be observed, and when it is smaller, the excessive diffusion will appear.A better formula of D i could be designed to ensure the mixing is sufficient so that the aliasing is eliminated entirely.
When parcel i needs to be mixed, its major axis will be shrunk by a factor α = 0.05 to reduce γ i ; so the mixing is a gradual process.Firstly s 1i is reduced to (1 − α)s 1i .Secondly the new H i are reconstructed by using Eq. ( 13).The skeleton points are also reset according to the new H i .The constant value of α may be changed to a function of D i and time step size suggested by the reviewer Eigil Kaas, which will make LASM more adaptive.This will be addressed in the future works.Some part (controlled by α) of the tracer mass (omitting the species index m) carried by parcel i will be distributed to its surrounding parcels.The mixing weight for a surrounding parcel j is calculated on the local stereographic plane of parcel i as where d (1) j is the distance parallel to the major axis of parcel i, d (2) j is vertical to the major axis, parameters β 1 and β 2 are two scale factors, and S i is the set of the surrounding parcels.Weight (Eq.18) is similar to Eq. ( 29) in Kaas et al. (2013), so only parcels close to parcel i will be involved; see the color of the surrounding parcels in Fig. 4c (The redder the color is, the larger the weight is).Moreover we can control the lateral mixing by adjusting β 2 .The larger β 2 is, the less the lateral mixing is.In addition, the lateral mixing can also be controlled by D i , so when D i is high, β 2 is changed to β * 2 (e.g., 10) to allow more lateral mixing.Then the new tracer mass of parcel j is m * j = m j + w j αm i .The volume of parcel j will be increased accordingly (H j is updated) as so the tracer density will be which is clearly a weighted average and reveals the mixing nature of this operation.No new extreme of tracer density will be introduced.The total mass on the parcels is also conserved in this process.
Compared to the mixing in other Lagrangian schemes, LASM can avoid the mixing as much as possible when the flow deformation is quite linear because the linear deformation approximation is sufficient for capturing the parcel shape changes and then helping to remap the tracer mass onto the right grids.This is the reason why γ m can be set as large as 100 so that the evolution of the parcel shape will not be interrupted when the flow deformation is in order.Only when this approximation is violated will the interparcel mixing be triggered to compensate for such degeneration when the aliasing occurs.Additionally, extra mixing can be added by casting the corresponding formula based on the physical thoughts.

Implementation
LASM is implemented in C++ language and incorporates the object-oriented programming model.The Fortran interfaces for application in the existing models will be added in the future.The workflow is shown in Fig. 5 and the overall programming units are depicted in Fig. 6, which are divided into two libraries (LASM and GEOMTK) with different functionalities.The key implementation details can be summarized as follows.
-Parcel and tracer codes: the Parcel class contains the centroid coordinate, H, U, S, V, and the index to the mesh cell.It is the core of LASM.The Tracer class derives from Parcel and adds the members related to tracer species, such as the density and mass.The skeleton points (TracerSkeleton) are also contained in Tracer.Tracer objects are managed by TracerManager, which contains a list of them.The tracer quantities remapped onto the mesh and the connectivity between parcels and grids are organized in MeshAdaptor.
-Mesh and field codes: the mesh is specified by the model (AGCM or OGCM) and is not the core of LASM, so the mesh details are separated from LASM.Different mesh could be used by implementing the corresponding mesh, index, and regrid classes.The field class are bound with the mesh class to provide a convenient container.Currently, only lat-lon mesh is implemented, but other meshes will be supported in the future.
-Neighbor search codes: the efficient neighbor search codes should depend on the mesh type, but currently a general search algorithm is adopted to reduce the development burden.This algorithm comes from the MLPACK library 2 and is based on the tree data structure.The fixed grids are organized into a tree structure at the initialization stage.When searching the connected grids of a parcel, this tree is walked through and the Euclidean space distance is used as the judge.In the future, if the efficiency needs to be improved further, the exclusive search algorithms for each mesh could be implemented.
-Interparcel mixing codes: the interparcel mixing is conducted after the remapping from parcels to grids, because this can avoid modification of connectivity be-2 See MLPACK website: http://www.mlpack.org.
tween parcels and grids due to the shape change of parcels.
-Test case codes: several test cases for testing LASM are coded in the same interface AdvectionTestCase.
Currently there are three concrete test case classes: The current implementation is not fully optimized; only the serial run is tested.The parallel version by using OpenMP and MPI techniques will be provided in other works.

Test case results
In the following subsections, three test cases are performed to validate LASM.It should be noted that the trajectory calculation in LASM is the same as it will be in a real application, whereas some other schemes used semi-analytical trajectory (i.e., use the true velocity at half time step), so the error caused by the time discretization will be a little larger  in LASM relatively.In order to reduce the unnecessary computation in the lat-lon mesh, the initial distribution of parcels is chosen as a reduced lat-lon mesh with only four parcels along the zonal circle closest to the poles, and the latitude of the transition from the normal region to the reduced region is 45  the interesting regions in the ideal test cases3 .The impact of interparcel mixing with different parameter configurations is tested in Sect.3.2.4.A test case with a barotropic model as the driver is also included to demonstrate the performance of LASM in the more complex flow.

Solid body rotation
The traditional solid body rotation on the sphere test case is used to illustrate the basic performance of LASM, since no apparent shape changes will occur.The sphere radius R is chosen as 1, and the cosine bell radius R c is R/3.The time step size of all runs is 1800 s. Figure 7 shows the traditional normalized errors (l 1 , l 2 , and l ∞ ) of three runs with rotation parameter α as π/2, π/4, and 0, respectively.The pole-crossing (α = π/2) run has rel- Figure 13.The effects of interparcel mixing with the slotted cylinders initial condition.There are three parameter configurations as shown on the top of the each row.The tested parameter that controls the amount of tracer mass distributed to the surrounding parcels is α.
atively larger errors, whereas another run with rotation along the Equator (α = 0) has significant error reduction.This reveals that the pole-crossing trajectory calculation affects the overall accuracy.On the other hand, HEL constructed on the cubed-sphere mesh performs with high accuracy, so the main inaccuracy sources are from lat-lon mesh (e.g., trajectory calculation on the polar stereographic projection plane).In the future, a more elaborate pole-crossing trajectory calculation may be designed, or computed on other meshes without pole problems.

Deformation flow
The deformation flow test cases proposed in Nair and Lauritzen (2010) and Lauritzen et al. (2012) are utilized here with the emphasis on convergence rate, filamentation preservation, and correlation preservation diagnostics.These tests are also used in a scheme intercomparison study (Lauritzen et al., 2014).The analytical flows selected are the same as the previous studies (Kaas et al., 2013), including one non-divergent flow and one divergent flow.The trajectories implied by the two flows are quite challenging, since the initial compact parcels will be stretched and bent into long filaments.The deformation reaches maximum at half simulation duration and is reversed after that.Three initial conditions consisting of Gaussian hills, cosine hills, and slotted cylinders defined in Lauritzen et al. ( 2012) are used for different purposes.
In all the following test runs, the sphere radius R is 1, and the normalized flow period T is 5.The Courant number (CN) is fixed for different spatial resolutions as about 1 with time step size 5/600 for 1.5 • spatial resolution.The selection of time step size is only subject to the accuracy, not the stability, compared with other Eulerian schemes.The simulated spatial distribution of tracer mixing ratios in the two flows with different initial conditions is shown in Fig. 10, with the parameters listed in Table 1, which shows that LASM can preserve the discontinuity in tracers and perform well in the deformational flows.The results with the slotted cylinders initial condition of other schemes can be referred to in Figs.7-10 in Lauritzen et al. (2014).

Convergence rate and minimal resolution
The numerical convergence rate of LASM is evaluated with the Gaussian hills' initial condition in the non-divergent flow.Figure 8 shows the results for error norms l 2 and l ∞ , where the rate is calculated by the linear regression from the four error samples.Currently, the convergence rate of LASM is first order with CN = 1.The analytical velocity field is also used, but the convergence rate is still first order (not shown), because the spatial remapping between tracers and grids is only first order.In contrast, HEL uses a second-order numerical approximation to the gradient of the tracer density on the Eulerian mesh (see Eq. 15 in Kaas et al., 2013).In the future, ρ mi and ρ mI on the right-hand side of Eqs. ( 16) and ( 17) can be modifed as in HEL to improve the convergence rate of LASM.
The minimal resolution λ m is an absolute error measure that supplements the convergence rate.The initial condition is changed to the cosine hills, which is quasi but not infinitely smooth, to challenge the schemes.λ m is the evaluated resolution at which the l 2 error is 0.033; therefore a larger value is better.For unfiltered CSLAM (conservative semi-Lagrangian multi-tracer transport scheme) (Lauritzen et al., 2010)  CN = 1 (Lauritzen et al., 2012), λ m is about 1.40 • for LASM with CN = 1.

Filamentation preservation
This diagnostics (l f ) is used to indicate how well the advection scheme resolves the thin filaments that develop at t = T /2, which is defined as where A(τ, t) is the total area whose tracer mixing ratio value exceeds τ at time t.Nineteen values across the range [0.1, 1] of τ are selected for the cosine hill initial condition.In nondivergent flow, the the area spanned by mixing ratio values larger than any given threshold value should be preserved.Therefore l f of a good scheme should be close to 100 %.When evaluating l f on the parcels, the values would be almost 100 % as stated in Kaas et al. (2013), which is also true for LASM.But because we are concerned with the results on the model grids, l f is evaluated on the model grids even though LASM is fully Lagrangian.
Figure 9 shows the result of two resolutions (1.5 • and 0.75 • ), which indicates that LASM can preserve the filament or large gradient very well.All the values of l f are around 100 % with a little oscillation at 1.5 • resolution, and even better when the resolution is doubled.The results of other schemes can be found in Fig. 5 of Lauritzen et al. (2014), where most schemes exhibit diffusive trait; i.e., l f is larger than 100 % at smaller τ and vice versa.

Correlation preservation
The correlation preservation diagnostics proposed in Lauritzen and Thuburn (2012) and Lauritzen et al. (2012) evaluate the potential performance of schemes in the chemistry transport models, where the pre-existing correlations among tracer species are chemically significant and influence the chemical reaction rates and equilibria.In this test two tracers with initial nonlinear functional relation are advected in the non-divergence flow, where χ corresponds to the cosine hills' initial condition and ξ corresponds to the nonlinearly related hills (see Eqs. 21 and 22 in Lauritzen and Thuburn, 2012, for setup).
The scatterplots of the two tracers on the mesh at t = T /2 are depicted in Fig. 11, including three mixing diagnostics: real mixing l r , unmixing l u and overshooting l o .The order of magnitude of unmixing is 10 −8 .When evaluated on the parcels, the unmixing and overshooting are both 0 as HEL.The real mixing is also small ( 10 −4 ), so the scatter points are almost along the prescribed functional curve.The results of other schemes can be found in Figs.11-14 of Lauritzen

Impact of interparcel mixing
The above deformation test cases have not included apparent interparcel mixing, because the mixing only occasionally occurs near the poles, where the tracer density is constant, so the impact of the mixing are demonstrated in this section by adjusting the key parameters.The selected test case is the non-divergent flow with slotted cylinders initial condition.Firstly γ m is decreased to 10 (γ * m is 5).The results with different β 2 (β * 2 is 10) are in Fig. 12.The effects of mixing is obvious, especially for β 2 = 10, which causes intensive lateral mixing.By increasing β 2 , the lateral mixing is largely inhibited, so the simulation result at time t = T with β 2 = 1000 is acceptable for this test case.The results of HEL can be found in Fig. 7 of Kaas et al. (2013).The correlation preservation results at time t = T /2 are also expected (see left columns of Fig. 12).Secondly we test the impact of parameter α, which affects the amount of distributed tracer mass at each mixing event.As shown in Fig. 13, increasing α will cause some noise, so it is better to use a small α (0.05) to make the mixing occur gradually.
In the deformation test cases, LASM can already produce good results without introducing apparent interparcel mixing.There are two reasons: (1) the linear deformation approximation of the parcel shape is effective; (2) the two flows are ideal, so the local deformation can be captured by the approximation (the disorder degree D of parcels is low).But when the flow is complicated and the parcels are heavily elongated, this approximation will be violated (the disorder degree D of parcels is high), especially near the strong vor- tex.In the next, more realistic barotropic test, the adaptivity of the interparcel mixing will be shown.

Barotropic model test
To test the new LASM in a more realistic flow, this section shows the results from a barotropic dynamical model on the sphere.This dynamical model is solved by using a finite difference method with an implicit Euler midpoint time integrator that can conserve the total energy and total mass on the regular lat-lon mesh.The details about the numerical techniques in the model can be referred to in Wang and Ji (1994).
The barotropic equations are where H = √ φ, U = uH , V = vH , v * = v cos ϕ, and f * = 2 sin ϕ + u a tan ϕ. u, v, and φ are the zonal wind speed, meridional wind speed, and geopotential depth, respectively; φ s is the surface geopotential.The gravity acceleration g = 9.8 m s −2 , the Earth radius a is 6.371 × 10 6 m, and the rotation velocity of the Earth is 7.292 × 10 −5 s −1 .The model spatial resolution is 1.5 • × 1.5 • , and the time step size is 20 s to ensure stability.
The objective of this experiment is to test the effectiveness of the adaptive interparcel mixing.Three tracers are advected: -an initial uniform background tracer used to calculate mixing ratio of the step tracer; -the geopotential depth as in the barotropic model; -a step tracer (Fig. 14b).
Since the current LASM is not fully dynamical, the standard shallow water tests are not used, instead two subcases with different initial conditions and topography are conducted: -subcase 1: isolated cosine topography (Fig. 14a), constant initial geopotential depth, and wind field derived from geostrophic relation; -subcase 2: flat topography, initial wind field, and geopotential depth (Fig. 14c) from ERA reanalysis.Although the initial flow in subcase 1 is simple, a strong vortex will be generated above the topography.The interparcel mixing must be effective to control the parcel shape near the vortex to avoid intensive filamentation, which will degenerate the degree of approximation.The results of subcase 1 are shown in Fig. 15.Two resolutions of parcels are used (21 984 and 87 948).From the geopotential depth results, it can be concluded that the mixing is working, despite some noise around the vortex in LASM compared with the barotropic model.The noise level is reduced when the number of parcels is increased.Based on this observation, the adaptive refinement of parcels may be incorporated in the future to further improve the performance of LASM.From the step tracer results, it is obvious that LASM does a good job simulating the discontinuous tracer field, where the thin and rolled filament is preserved very well.The effects of the interparcel mixing adaptivity be further illustrated in Fig. 16.When D * is reasonable enough, the sensitivity on it is not very large, but when it is infinity -i.e., the adaptivity is turned off -the tracer density is in disorder and many spotty structures come out due to the degeneration of linear deformation approximation when parcels turn to be very long.The parcel distribution after 24 h is depicted in Fig. 15f.The parcels are stirred heavily by the vortex, and the mixing occurs frequently there to ensure the good approximation.
In subcase 2, a 200 hPa geopotential depth and wind field from ERA reanalysis are chosen as the initial conditions.After 24 h, the geopotential depth simulated by the barotropic model and LASM are compared in Fig. 17a-c.The discrepancy (i.e., wind-mass inconsistency) between the results of the barotropic model and LASM can not be avoided, but the difference in percentage is mostly within ±2 %.The step tracer in Fig. 17d also indicates the good performance on the discontinuous tracer of LASM, where the discontinuous tracer is deformed heavily, although no reference solution is given.The parcels after 24 h (Fig. 17e) are fairly chaotic.The time duration of this subcase is also extended to 10 days, and the step tracer is already deformed beyond recogition as shown in Fig. 15f (the animation can be found in the Supplement or viewed at http://dongli.github.io/dongli/assets/pictures/lasm.barotropic.240x120.step.cells.gif).
The sensitivity about γ m is also tested by reducing it to 10, and the results (not shown) of the two subcases are almost the same as the above, because the interparcel mixing is active to limit the parcel size in the nonlinearly deformational flow, so no long parcels survive.

Conclusions
A new Lagrangian advection scheme, LASM, is proposed in this work, which takes into account of the pros and cons of TTS-C/I and other Lagrangian schemes.The core of LASM is the linear transformation matrix that is used to describe the parcel shape.The remapping weight is then anisotropic so that the flow deformation is respected.The test results show that LASM can reduce the aliasing error that will occur in the Lagrangian schemes without introducing substantial mixing when the flow deformation is quite linear.The interparcel mixing is only triggered in the regions with very strong nonlinear deformation where the linear deformation approximation will be violated and the parcels will be shrunk to improve the approximation.The effects of such mixing is verified and acceptable as compared with HEL.Parameters γ m , α, β 1 , and β 2 can be used to control the degree of mixing.Furthermore the disorder degree D of the parcels is used to adjust γ m and β 2 .This adjustment is evaluated to be effective in the barotropic tests.The formulation of the parameters may be further polished to better control the mixing and judge the nonlinear deformation.Currently, the convergence rate of LASM is first order, and this can be improved by adding the second-order terms when remapping tracer mass.
LASM is ready to be extended into three-dimension with consideration of the vertical motion and boundaries.Other spherical meshes except for lat-lon mesh will be supported in the near future, and the pole-crossing performance will be better when the mesh does not have the pole problems.The grid searching procedure may be further polished to improve the remapping accuracy.The remaining challenges come from the parallel computing.In the Eulerian models (AGCM, OGCM, CTM), the spatial domain is decomposed into several parts, and each part is assigned to one CPU per node.For some Eulerian or semi-Lagrangian schemes, only boundary cells need to be communicated among the CPUs.But because Lagrangian parcels move around, the parcels in CPU1 initially may travel to the other CPU, and this may cause some trouble when parallelizing the codes.Efficient LASM for application in the real models will be implemented in other works.
Figure 1.(a)Linear transformation between physical space and remapping space.In this 2-D case, the black ellipse in physical space is transformed into a red circle in remapping space so that two points marked by 1 and 2 have the same distance to x 0 in remapping space, but their Euclidean distances to x 0 are different.(b) The deformed shape function is bell-like.

Figure 3 .
Figure3.The neighbor grids of the parcel schema.The black cross is the parcel centroid, and the black circle represents the parcel shape that is determined by the deformation matrix.The red filled circles are the neighbor grids connected with the parcel.

L.
Figure 4. (a)The low-disorder-degree case, where the black ellipses are the surrounding parcels of the red one.(b) The high-disorder-degree case.(c) The blue parcel is mixed with its surrounding parcels, whose colors indicate the mixing weights (the redder, the larger).The blue parcel is shrunk to the green one after mixing.

Figure 6 .Figure 7 .
Figure 6.The overall programming units of LASM and their relationship.Four external libraries (Boost, Armadillo, NetCDF, and MLPACK) are used.

Figure 8 .FilamentFigure 9 .
Figure8.Test of convergence for the deformation test case with Gaussian initial condition using error norms l 2 (left) and l ∞ (right).The Courant number is fixed.The upper reference line is for second order, and the lower one is for third order.

Figure 10 .
Figure 10.Results from deformation flow test cases.(a) is the initial slotted cylinders for the non-divergent flow test, and (b) is the initial cosine hills for the divergent flow test.(c) and (d) are the simulated distributions at time t = T /2 with maximum deformation.(e) and (f) are the final distributions at time t = T .The model grid resolution is 1.5 • × 1.5 • , and time step size is T /600.

Figure 11 .
Figure 11.Diagnostics for correlation preservation for the model grid resolution 1.5 and 0.75 • .

Figure 12 .
Figure 12.The effects of interparcel mixing with the slotted cylinders' initial condition.There are three parameter configurations as shown on the top of the each row.By increasing β 2 , the lateral mixing is decreased.The results for γ m = 100 are in Figs. 10 and 11.
Figure 14.(a) The surface geopotential in subcase 1.(b) The initial step tracer mixing ratio.(c) The initial geopotential depth in subcase 2 (wind field is omitted).

Figure 15 .
Figure 15.The simulation results after 24 h in subcase 1.

Figure 16 .
Figure 16.The effectiveness of the interparcel mixing adaptivity.Different disorder degree D * are tested.(a) is the same as Fig. 15b.

Figure 17 .
Figure 17.The simulation results after 24 h in subcase 2.