Eugene Kazantsev

Optimal Boundary Conditions for ORCA-2 Model

Keywords: Variational Data Assimilation, Boundary conditions, ORCA-2 model

A 4D-Var data assimilation technique is applied to a ORCA-2 configuration of the NEMO in order to identify the optimal parametrization of the boundary conditions on the lateral boundaries as well as on the bottom and on the surface of the ocean. The influence of the boundary conditions on the solution is analyzed as in the assimilation window and beyond the window. It is shown that optimal conditions for vertical operators allows to get stronger and finer jet streams (Gulf Stream, Kuroshio) in the solution. Analyzing the reasons of the jets reinforcement, we see that the major impact of the data assimilation is made on the parametrization of the bottom boundary conditions for lateral velocities u and v.

Automatic generation of the tangent and adjoint codes is also discussed. Tapenade software is shown to be able to produce the adjoint code that can be used after a memory usage optimization.

Introduction

Variational data assimilation technique, first proposed in [START_REF] Dimet | A general formalism of variational analysis[END_REF], [START_REF] Dimet | Variational algorithm for analysis and assimilation of meteorological observations. theoretical aspects[END_REF], is based on the optimal control methods [START_REF] Lions | Contrôle optimal de systèmes gouvernés pas des équations aux dérivées parielles[END_REF] and perturbations theory [START_REF] Marchuk | Formulation of theory of perturbations for complicated models[END_REF]. This technique allows us to retrieve an optimal data for a given model from heterogeneous observation fields. Since the early 1990's, many mathematical and geophysical teams are involved in the development of the data assimilation strategy. One can cite many papers devoted to this problem, as in the domain of development of different methods for the data assimilation and also in the domain of its applications to the atmosphere and oceans.

In the beginning, data assimilation methods were intended to identify and reconstruct an optimal initial state of the model. However, the idea that other model's parameters should also be identified by data assimilation has also been studied and discussed in numerous papers. One can cite several examples of using data assimilation to identify the bottom topography of simple models ( [START_REF] Losch | Bottom topography as a control variable in an ocean model[END_REF], [START_REF] Heemink | Inverse 3d shallow water flow modelling of the continental shelf[END_REF], [START_REF] Kazantsev | Identification of optimal topography by variational data assimilation[END_REF]), to control open boundary conditions in coastal and regional models ( [START_REF] Shulman | Local data assimilation in specification of open boundary conditions[END_REF], [START_REF] Shulman | Optimized boundary conditions and data assimilation with application to the m2 tide in the yellow sea[END_REF], [START_REF] Taillandier | Controlling boundary conditions with a four-dimensional variational data-assimilation method in a non-stratified open coastal model[END_REF], [START_REF] Brummelhuis | Identification of shallow sea models[END_REF]), boundary conditions on rigid boundaries ( [START_REF] Dimet | Retrieval of balanced fields: an optimal control method[END_REF], [START_REF] Leredde | On initial, boundary conditions and viscosity coefficient control for burgers' equation[END_REF] [START_REF] Lellouche | Boundary control on burgers equation. a numerical approach[END_REF], [START_REF] Kazantsev | Identification of an optimal boundary approximation by variational data assimilation[END_REF], [START_REF] Kazantsev | Optimal boundary discretisation by variational data assimilation[END_REF], [START_REF] Kazantsev | Boundary conditions control for a shallow-water model[END_REF]) and to determine other parameters of a model ( [START_REF] Zou | An optimal nudging data assimilation scheme using parameter estimation[END_REF], [START_REF] Panchang | On the Determination of Hydraulic Model Parameters Using the Strong Constraint Formulation Modeling Marine Systems[END_REF], [START_REF] Chertok | Variational data assimilation for a nonlinear hydraulic model[END_REF]).

It was pointed out in [START_REF] Navon | Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography[END_REF] that the problem of adjoint parameters identification is frequently ill-posed. This fact remains valid when boundary conditions are considered as a control parameter. It was shown in [START_REF] Kazantsev | Identification of an optimal boundary approximation by variational data assimilation[END_REF] that the presence of non-null kernel of the Hessian results in a non-unique choice of optimal boundary conditions. However, all sets in the kernel from the kernel are equivalent: they provide the same (or almost the same) cost function's value and almost the same evolution of the model's solution after the end of assimilation. Consequently, we can not talk about physical meaning of the optimal conditions. Indeed, any value from the kernel can be obtained as the result of the minimization procedure while only one of them may have a physical meaning. In this paper we take this ambiguity into account and consider the adjoint boundary condition estimation as a simple compensation of the model errors rather than real identification of model parameters.

Particular attention is payed to the design of the tangent and adjoint codes. As it has been shown in [START_REF] Kazantsev | Optimal boundary discretisation by variational data assimilation[END_REF], the derivative of the model with respect to boundary conditions is two or three times longer (as well as in terms of the development, the number lines of the code and the necessary CPU time) than the derivative used to control the initial point of the model. The usual tangent model describes a linear evolution of a small perturbation to a model solution while the derivative with respect to other parameter includes also a block that introduce the perturbation into the model. This block may be at least as long and complex as the whole derivative with respect to initial conditions.

That's why we focus our attention on the automated differentiators. In this paper tangent and adjoint codes have been obtained by Tapenade software [START_REF] Hascoët | Tapenade 2.1 user's guide[END_REF]. Current version of this software reveals to be able to produce the derivative of such a complex code as ORCA-2.

The paper is organized as follows. In the second section, 2 0 global ocean configuration (ORCA-2) 1of the Nucleus for European Modelling of the Ocean (NEMO)2 model is described. The third section is devoted to the data assimilation technique and to the Automatic Differentiation Engine Tapenade3 software to generate the adjoint and to the optimization of the generated code. Results and discussion are presented in the fourth section.

2 ORCA-2 configuration of the NEMO and its discretization.

ORCA-2 configuration of the Nucleus for European Modelling of the Ocean [START_REF] Madec | Nemo ocean engine[END_REF] is used in this paper. This is 2 0 global ocean configuration based on the OPA 8.2 [START_REF] Madec | Opa 8.1 ocean general circulation model reference manual[END_REF] primitive equations model.

Its domain extends from 78 0 S to 90 0 N. The model grid counts 182 × 149 × 31 nodes. Vertical discretization is performed on z levels with the partial step approximation of the bottom cell [START_REF] Barnier | Impact of partial steps and momentum advection schemes in a global ocean circulation model at eddy-permitting resolution[END_REF]. Vertical mixing is achieved using the turbulent kinetic energy (TKE) scheme described in [START_REF] Blanke | Variability of the tropical atlantic ocean simulated by a general circulation model with two different mixed layer physics[END_REF].

Spatially discretized equations of the ORCA-2 model are written as follows:

∂u ∂t = S x S y v S y (ω + f ) -D x S x u 2 + S y v 2 2 -S z
+ D zz A z u u Vertical diffusion + gD x (η + T c ∂ t η)
Surface Pres. Grad.

(1)

∂v ∂t = S x S y u S x (ω + f ) -D y S x u 2 + S y v 2 2 -S z S x wD z v Advection + D x A xy v ξ + D y A xy v ω Lateral diffusion + + g z surf ace D y S z ρ(x, y, ζ)dζ Hydrostatic Pressure grad. + D zz A z v v Vertical diffusion + gD y (η + T c ∂ t η) Surface Pres. Grad. (2) ∂T ∂t = -D x (uS x T ) -D y (vS y T ) -D z (wS z T ) Advection + A xy T D x D x T + D y D y T Lateral diffusion + + D zz A z T T Vertical diffusion +Solar Radiation + Geothermal Heating + BBL (3)
∂s ∂t = -D x (uS x s) -D y (vS y s) -D z (wS z s) Advection + A xy T D x D x s + D y D y s Lateral diffusion + + D zz A z T s Vertical diffusion +BBL (4) ∂η ∂t = w(x, y, z = surface) (5)
ξ = D x u + D y v, ω = D y u -D x v, w = z bottom ξ(x, y, ζ)dζ, ρ = ρ(T, s) (6)
where operators D and S are derivatives and interpolations on the Arakawa C-grid. These operators will be discussed below in details.

The set of variables in this system is composed of u, v and w that represent zonal, meridional and vertical velocity components, T and s denote the potential temperature and salinity, ξ and ω are horizontal divergence and vorticity, η is the sea surface elevation and ρ is the density anomaly that is defined as a function of the temperature and salinity by the state equation. As one can see, u, v, T, s, η are prognostic variables while w, ξ, ω and ρ are diagnostic ones.

Among parameters in these equations, we can see f = 2Ω sin(latitude) that represents the Coriolis parameter, gravity acceleration g = 9.81 m s 2 , lateral diffusion coefficients

A xy T = 2000 m 2 s , A xy u = 40000 m 2 s latitude > 15 0 2000 m 2 s latitude < 15 0 , .
To calculate coefficients of the vertical diffusion A z u , A z v , A z T , A z s we use the turbulent closure scheme accompanied by the double diffusive mixing and enhanced vertical diffusion approximations.

The turbulent closure scheme is applied to solve the problem of statically unstable density profiles. The vertical eddy viscosity and diffusivity coefficients are computed from a TKE turbulent closure model based on a prognostic equation for ē, the turbulent kinetic energy [START_REF] Hascoët | Tapenade 2.1 user's guide[END_REF], and a closure assumption for the turbulence length scales [START_REF] Madec | Opa 8.1 ocean general circulation model reference manual[END_REF].

∂ē ∂t = A z v ∂u ∂z 2 + ∂v ∂z 2 -A z T N 2 + ∂ ∂z A z u ∂ē ∂z -c ǫ ē3/2 l ǫ (7)
A z u = A z v = max(A z 0 , C k l k √ ē) (8)
A z T = A z u /P rt ( 9
)
where N is the local Brunt-Vaisälä frequency which is calculated as a function of T and s. Parameters l ǫ and l κ are the dissipation and mixing length scales, P rt is the Prandtl number. The constants C k = √ 2/2 and c ǫ = 0.1 are designed to deal with vertical mixing at any depth [START_REF] Gaspar | A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station papa and long-term upper ocean study site[END_REF]. Following [START_REF] Blanke | Variability of the tropical atlantic ocean simulated by a general circulation model with two different mixed layer physics[END_REF], P rt is defined as a function of the local Richardson number, R i :

P rt =    1 if R i ≤ 0.2 5 R i if 0.2 ≤ R i ≤ 2 10 if 2 ≤ R i
In frames of the enhanced vertical diffusion parameterization, we assign very big values to the vertical eddy mixing coefficients in regions where the stratification is unstable (i.e. when the Brunt-Vaisälä frequency is negative) [START_REF] Lazar | The deep interior downwelling, the veronis effect, and mesoscale tracer transport parameterizations in an ogcm[END_REF]. This is done on both momentum u, v and tracers T, s:

A z = { if N 2 < 0 then A z = 100} (10)
Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa. They contribute to diapycnal mixing in extensive regions of the ocean. The parameterization of such phenomena was included in a global ocean model and it was shown that it leads to relatively minor changes in circulation but exerts significant regional influences on temperature and salinity in [START_REF] Merryfield | A global ocean model with double-diffusive mixing[END_REF].

A z T = A z T +    0.7 10 -4 Rρ(1+(Rρ/1.6) 6 ) if R ρ > 1 A ddm if 0 < R ρ < 1 0 otherwise A z s = A z s +    10 -4 1+(Rρ/1.6) 6 if R ρ > 1 (1.85 R ρ -0.85) A ddm if 0.5 ≤ R ρ < 1 0.15R ρ A ddm if 0 < R ρ < 0.5 (11)
A z u = max(A z u , A z T , A z s ) A z v = max(A z v , A z T , A z s )
where A ddm = 1.3636 × 10 -6 e 4.6e -0.54(1/Rρ -1) and R ρ is the buoyancy ratio

R ρ ∼ ∂ z T /∂ z S.
The term T c ∂ t η in the equations ( 1) and ( 2) is introduced to dump the external gravity waves. These waves are fast so their timescale is short with respect to other processes described by the primitive equations. Explicit resolution of these waves requires too short time step which is unnecessary to all other physics. Consequently, the filter of temporally unresolved external gravity waves, proposed in [START_REF] Roullet | salt conservation, free surface, and varying levels: a new formulation for ocean general circulation models[END_REF], is introduced into the model. The cutoff time T c is equal to one time step of the model.

The model is discretized on the grid, that is the generalization to three dimensions of the well-known "C" grid in Arakawa's classification [START_REF] Mesinger | Numerical methods used in Atmospheric models[END_REF]. The arrangement of variables is the same in all directions. It consists of cells centered on scalar points (T , s, η, ρ) with vector points (u, v, w) defined in the center of each face of the cells. The relative and planetary vorticity, ω and f , are defined in the center of each vertical edge.

To perform interpolations on this grid and to calculate the derivatives, we introduce operators S and D in the discretized equations ( 1)- [START_REF] Griewank | Achieving logarithmic growth of temporal and spatial complexity in reverse automatic differentiation[END_REF]. These operators plays the key role in this study because they depend on the boundary conditions and introduce them into the model. Despite these operators are usually written as overbars and δ brackets (see [START_REF] Madec | Nemo ocean engine[END_REF], for example), we note them as letters with index in order to emphasize that these are operators under control in this paper. Interpolations are calculated as a weighted mean of two function values in the adjacent nodes. Weights are defined to be proportional to the grid steps of corresponding cells in order to ensure the second order interpolation of a grid function. These weights are described in details in the [START_REF] Madec | Nemo ocean engine[END_REF], but we omit them in this paper to concentrate our attention on the principal features. Thus, writing operators in a simplified way, we assume both the argument and the result of the interpolation operator and the derivative are multiplied by an appropriate weight.

At each grid-node in the ocean they we write interpolations S x , S y and S z in a common way:

(S x u) i+1/2,j,k = u i,j,k + u i+1,j,k 2 , (S y u) i,j+1/2,k = u i,j,k + u i,j+1,k 2 , (S z u) i,j,k+1/2 = u i,j,k+1 + u i,j,k 2 (12)
Following [START_REF] Kazantsev | Identification of an optimal boundary approximation by variational data assimilation[END_REF], [START_REF] Kazantsev | Boundary conditions control for a shallow-water model[END_REF], these expressions are modified in the grid-nodes adjacent to the boundary, i.e. near the continents for interpolations in x and y directions and near the bottom and the surface for vertical interpolation. Let us suppose the index i = 0 corresponds to left rigid boundary. That means the index i = 1 is the first grid node in the ocean. In this case, we use the formula

(S x u) 1/2,j,k = α Sux l 0 + α Sux l 1 u 0,j,k + α Sux l 2 u 1,j,k (13)
to calculate the interpolated value of u at the point 1/2, j, k, where scalar variables T, s, ρ are defined. The value of (S x u) N -1/2,j,k near the right boundary is calculated by the similar formula, but with different coefficients Similarly, interpolated value of T, s, ρ at points i = 1, i = N -1 are calculated by

α Sux r 0 , α Sux r 1 , α Sux r 2 . u 0 u 2 u 3 u N u N -2 u N -3 T 3/2 T 5/2 T N -3/2 T N -5/2 T 1/2 T N-1/2 u 1 u N-1
(S x T ) 1,j,k = α ST x l 0 + α ST x l 1 T 1/2,j,k + α ST x l 2 T 3/2,j,k (S x T ) N -1,j,k = α ST x r 0 + α ST x r 1 T N -1/2,j,k + α ST x r 2 T N -3/2,j,k
Coefficients α play the role of the control variables in this paper. Operators are allowed to change their properties near boundaries in order to find the best fit with requirements of the model and data. To assign all control variables we shall perform data assimilation procedure and find their optimal values.

The first coefficient, α 0 , is added to the interpolation formula to simulate non-uniform boundary conditions. Coefficient α 1 controls the contribution of the physical boundary value that is prescribed to u 0,j,k . It should be noted here, that if impermeability condition is imposed, u 0,j,k = 0, the interpolation (13) is controlled by two α only.

Similar modification of the interpolation formulas are performed near the north and the south boundaries in the y direction. Exception is made for the periodical conditions. They are applied on the 78 0 E longitude and near the North Pole. In this case no modification is made and no control is applied.

In the vertical direction we control the discretization at two points: one point near the surface and another one near the bottom (see fig. 2). To calculate the interpolated values of the temperature, for example, we use the following formula:

(S z T ) i,j,1 = α SzT s 0 + α SzT s 1 T i,j,1/2 + α SzT s 2 T i,j,3/2 (S z T ) i,j,K-1 = α SzT b 0 + α SzT b 1 T i,j,K-1/2 + α SzT b 2 T i,j,K-3/2 (14)
We note here, that the value of the temperature (and other variables defined at the temperature levels) is not extrapolated to nodes w 0 and w K where physical boundary conditions must be prescribed.

• • • w 0 hz 1/2 w 1 hz 3/2 w 2 hz 5/2 w 3 w K hz K-1/2 w K-1 hz K-3/2 w K-2 w K-3 T 1/2 hz 1 T 3/2 hz 2 T 5/2 hz 3 T 7/2 T K-1/2 hz K-1 T K-3/2 hz K-2 T K-5/2
Figure 2: Structure of the vertical grid. The number of layers depends on the longitude and latitude:

K = K i,j
These conditions participate in the interpolation from w levels to T levels

(S z w) i,j,1/2 = α Szw s 0 + α Szw s 1 w i,j,0 + α Szw s 2 w i,j,1 (S z w) i,j,K-1/2 = α Szw b 0 + α Szw b 1 w i,j,K + α Szw b 2 w i,j,K-1 (15)
Discretized derivatives near the boundary are defined in a similar way as for the horizontal derivatives and for the vertical ones.

(D x u) 1/2,j,k = α Dux l 0 -α Dxu l 1 u 0,j,k + α Dxu l 2 u 1,j,k (D z w) i,j,K-1/2 = α Dzw b 0 + α Dzw b 1 w i,j,K -α Dzw b 2 w i,j,K-1 (16)
As well as above, impermeability conditions are prescribed at lateral boundary for u and v. That means coefficients α 1 , being multiplied by the vanishing value of the function on the boundary, are not under control in the experiment. On the other hand, boundary condition on the surface of the ocean, prescribed for vertical derivatives, are really used in the model and controlled by corresponding coefficient α Dz 1 . Taking into account the fact that at different points optimal boundary conditions may be different, we accept the coefficients α may vary from point to point. Thus, α Dz and α Sz are considered as functions of the longitude and the latitude. Coefficients used in horizontal operators are also allowed to vary from one boundary point to another.

Along with the derivatives and interpolations, boundary conditions influence also the second derivatives in the vertical diffusion and reconstruction of the vertical velocity w from the horizontal divergence [START_REF] Griewank | Achieving logarithmic growth of temporal and spatial complexity in reverse automatic differentiation[END_REF]. Discretized equivalent of the (6) writes

w i,j,K-1 = α w b 0 -α w b 1 hz i,j,K-1 ξ i,j,K-1 w i,j,k-1 = w i,j,k -hz i,j,k-1 ξ i,j,k-1 ∀k : 2 ≤ k ≤ K -1 ( 17
)
w i,j,1 = w i,j,2 + α w s 0 -α w s 1 hz i,j,1 ξ i,j,1
One can see the impermeability condition on the bottom may be violated by the term α w b 0 in the approximation of the vertical velocity near the bottom w i,j,K-1 that may allow non-zero vertical velocity. Coefficient α w 1 allows to control the thickness of the surface and the bottom layers. Control of the boundary conditions for the second derivatives D zz in (1) is performed in a slightly different way. We add α 0 to the prescribed boundary conditions on the surface and on the bottom

∂u ∂z surf ace = α DzzU s 0 + τ x hz 1 ρ 0 , ∂v ∂z surf ace = α DzzU s 0 + τ y hz 1 ρ 0 , ∂T ∂z surf ace = ∂S ∂z surf ace = α DzzT s 0 u| bottom = v| bottom = α DzzU b 0 T | bottom = S| bottom = α DzzT b 0 ( 18
)
and modify the approximation of the second derivatives at the nodes, adjacent to the surface and to the bottom multiplying the finite differencing coefficients by α 1 and α 2 . Thus,for example, approximation of

∂ 2 u ∂z 2 writes D zz u i,j,1/2 = (A z u ) 1 hz 1 hz 1/2 (α DzzU s 2 u 3/2 -α DzzU s 1 u 1/2 ) ( 19
)
D zz u i,j,k-1/2 = 1 hz k-1/2 (A z u ) k hz k (u k+1/2 -u k-1/2 ) - - (A z u ) k-1 hz k-1 (u k-1/2 -u k-3/2 ) ∀k : 2 ≤ k ≤ K -1 D zz u i,j,K-1/2 = 1 hz K-1/2 α DzzU b 2 (A z u ) K-1 hz K-1 u K-1/2 - -α DzzU b 1 (A z u ) K hz K + (A z u ) K-1 hz K-1 u K-3/2
The value of hz k corresponds to the difference of depths between adjacent layers where u, v, T, S are defined (see fig. 2). The grid step hz k+1/2 = hz k + hz k+1 2 is the distance between layers where the vertical velocity w k and w k+1 are defined. The vertical grid is not uniform, so the gridsteps are not equal one to another.

As well can see, impermeability is no longer imposed on the bottom and the temperature and salinity of the bottom are modified. The wind tension on the top is controlled for velocities and an additional control flux of T and S is added on the surface. Parameters α 1 and α 2 helps also to control the position of the boundary modifying the depth of the first and the last layers.

Thus, controlling boundary conditions in this paper, we control, in fact, the discretization of operators near the boundary determined by the set of coefficients α. In this paper, we look for optimal values of α for the following 24 operators that approximate either derivatives or interpolations in the horizontal plane: We can see the number of control variables is comparable in both cases. These numbers are also comparable with the dimension of the system state (1 707 245 variables) that we need to control identifying the initial conditions of the model.

Argument u v T, s ω defined at Operator D x D x u divg. D x v vort. D x (u 2 + v 2 ) adv. D x ω dissip. D x (uS x T, uS x s) D x ρ press.grad. D y D y u vort. D x v divg. D y (u 2 + v 2 )

Data assimilation.

One of the principal purposes of variational data assimilation consists in the variation of control parameters in order to bring the model's solution closer to the observational data. This implies the necessity to measure the distance between the trajectory of the model and the data. Introducing the cost function, we define this measure. Generally speaking, the cost function is represented by some norm of the difference between the model's solution and observations accompanied by the difference with the background that is used as regularization term.

Cost function

To define the cost function we introduce dimensionless state vector φ that is composed of five variables of the model φ = {w u u, w v v, w T T, w s s, w η η} t weighted by coefficients w. The choice of these weights that are used to normalize the model variable will be discussed later. The distance between the model solution and observations is defined as the Euclidean norm of the difference

ξ 2 = ξ 2 (φ(p, t)) = m ((Hφ) m -φ obs m ) 2 = ( 20
)
where H is the operator that interpolates the model solution to the observation point and the sum is performed over all available observations at time t.

In this expression, we emphasize the implicit dependence of ξ on time and on the set of the control parameters p that is composed of

• the set of initial conditions of the model φ 0 = {u | t=0 , v | t=0 , T | t=0 , s | t=0 , η | t=0 },
• the set of the coefficients α xy that controls the discretizations of horizontal operators (Table 1) in the vicinity of the continents,

• the set of the coefficients α z that controls the discretizations of vertical operators (Table 2) near the bottom and near the surface,

• Vertical diffusion coefficients A z u , A z v , A z T , A z s
The cost function is composed of two terms. One of them measures the distance from the observations, and another one is added to avoid irregular solution that may occur due to lack of observational information and the ill-posedness of the problem. Indeed, the dimension of the model state is bigger than the quantity of available observations. That means we can not identify the model state in regions where observational information is absent. To avoid this problem, we add the background term in the cost function. This term allows us to determine the solution everywhere requiring that it must be close to the specified background. The regularization term we use has a form

I bgr (p) = m (p m -p bgr m ) 2 (21)
where p m is the set of parameters under control and p bgr m are background values of these parameters. In this paper we use initial guesses as background for all parameters. That means, when we look for the optimal initial state of the model, we use original unperturbed initial conditions as initial guess for the minimization and as the background as well. Looking for optimal α, we use combinations α 0 = 0, α 1 = α 2 = 1 for derivatives and α 0 = 0, α 1 = α 2 = 1/2 for interpolations both for initial guess and for the background.

The difference model-observations [START_REF] Lellouche | Boundary control on burgers equation. a numerical approach[END_REF] contribute to another component of the cost function. Taking into account the results obtained in [START_REF] Kazantsev | Optimal boundary discretisation by variational data assimilation[END_REF], we define the cost function as

I obs (p) = T 0 tξ 2 (φ(p, t))dt (22)
that gives a bigger importance to the difference ξ 2 at the end of assimilation interval. So far, we perform the data assimilation in order to make a forecast, we need a "better" estimate of the model state at the end of the assimilation window because this state is used as the initial point for the forecast that starts just after the assimilation. For this purpose, we require the model to go closer to observations at the end of the assimilation window increasing the weight of the distance in the cost function.

To search for a minimum of the cost function, we shall use its gradient with respect to control parameters. The gradient of the background term (21) can be calculated easily

(∇I) bgr n = ∂I bgr ∂p n = 2(p n -p bgr n ) ( 23
)
The gradient of the second component of the cost function [START_REF] Lions | Contrôle optimal de systèmes gouvernés pas des équations aux dérivées parielles[END_REF] can be calculated as a Gateaux derivative of an implicit function:

(∇I obs ) n = ∂I obs ∂p n = T 0 t ∂ξ 2 ∂p n dt = T 0 t m ∂ξ 2 ∂φ m ∂φ m ∂p n dt = = 2 T 0 t m ((Hφ) m -φ obs m ) ∂φ m ∂p n dt (24)
because the derivative ∂ξ 2 ∂φ m can easily be calculated from (20): ∂ξ 2 ∂φ m = 2((Hφ) m -φ obs m ). The second term in [START_REF] Madec | Opa 8.1 ocean general circulation model reference manual[END_REF], ∂φ m ∂p n , represents the matrix of the tangent linear model that relates the perturbation of mth component of the model state vector φ m to the perturbation of the parameter p n . This relationship, of course, is assumed in the linear approach, that means it is only valid for infinitesimal perturbations.

In the classical case, when initial conditions are considered as the only control variable, the derivative ∂φ(t) ∂p = ∂φ(t) ∂φ 0 is the classical tangent model that describes the temporal evolution of a small error in the initial model state. The matrix is a square matrix that is widely studied in numerous sensitivity analyses. Its singular values at infinite time limit are related to well known Lyapunov exponents that determine the model behavior (chaotic or regular) and the dimension of its attractor.

In our case, the matrix ∂φ(t) ∂p is rectangular in general. It describes the evolution of an infinitesimal error in any parameter (including initial state). However, we can study its properties in the similar way as we do with the classical tangent linear model. Its structure and composition is described in [START_REF] Kazantsev | Boundary conditions control for a shallow-water model[END_REF] for the case of using coefficients α as control parameters and in [START_REF] Kazantsev | Identification of optimal topography by variational data assimilation[END_REF] for the case when the bottom topography is used to control the model solution.

The product m ((Hφ) m -φ obs m ) ∂φ m ∂p n in (24) represents an unusual vector-matrix product. To calculate this product directly we would have to evaluate all the elements of the matrix. This would require as many tangent model runs as the size of the state vector is. So, instead of the tangent model, we shall use the adjoint one that allows us to get the result by one run of the model. Backward in time adjoint model integration that starts from (φ -φ obs ) provides immediately the product ∂φ ∂p * (φ -φ obs ) which is exactly equal to (φ -φ obs ) ∂φ ∂p in [START_REF] Madec | Opa 8.1 ocean general circulation model reference manual[END_REF]. Using these notations, we write

∇I obs = 2 T 0 t ∂φ(t) ∂p * (φ(p, t) -φ obs (t))dt ( 25
)
where the expression in the integral is the result of the adjoint model run from t to 0 starting from the vector (φ(p, t) -φ obs (t)).

Adjoint model

Tangent and adjoint models have been automatically generated by the Tapenade software [START_REF] Hascoët | Tapenade 2.1 user's guide[END_REF] developed by the TROPICS team in INRIA. This software analyses the source code of a nonlinear model and produces codes of its derivative ∂φ ∂p and of the adjoint ∂φ ∂p *

.

The obvious advantage of the automatic tangent and adjoint code generation consists in avoiding of huge development work. This fact is much appreciated especially in the case of control of internal model parameters other than initial conditions. As it has been shown in [START_REF] Kazantsev | Boundary conditions control for a shallow-water model[END_REF], the derivative of the model with respect to boundary conditions is composed by two blocks: one of them is composed of operators acting in the space of the model variables and another one is composed of operators linking the model variables and the controlled parameter. The first block is responsible for the evolution of a small perturbation by the models dynamics and is similar for any data assimilation. The second one determines the way how this perturbation is introduced into the model and is specific to the particular parameter under control. This block is absent when the goal is to identify the initial conditions of the model because the uncertainty in initial conditions is introduced only once, at the beginning of the model integration. But, when the uncertainty is presented in an internal model parameter, like in this case, the perturbation is introduced at each time step of the model.

Consequently, controlling internal model parameters we must write an additional block for the tangent model which is at least as complex as the whole tangent model for the initial point (see equations 5 and 6 in [START_REF] Kazantsev | Boundary conditions control for a shallow-water model[END_REF] where these blocks are compared for a shallow-water model). Moreover, controlling internal parameters we may be brought to the necessity to enlarge the control set, to study the influence of some other parameter. Simple automatic generation of tangent and adjoint codes helps us to be free in the choice of parameters to control.

On the other hand, automatically generated code frequently suffers from excessive requirements to computer memory. To be able to integrate the adjoint model backward in time, we have to keep the forward trajectory. There have been proposed in [START_REF] Tber | Building the tangent and adjoint codes of the ocean general circulation model OPA with the automatic differentiation tool tapenade[END_REF] to use Griewank and Walthers binomial checkpointing algorithm [START_REF] Griewank | Achieving logarithmic growth of temporal and spatial complexity in reverse automatic differentiation[END_REF], that assumes to keep forward model solution at several time-steps only and to recalculate all other time-steps during the adjoint model run. The spacing of time-steps where we keep the forward model solution should be chosen in order to minimize the calculation time. It has been shown in [START_REF] Tber | Building the tangent and adjoint codes of the ocean general circulation model OPA with the automatic differentiation tool tapenade[END_REF] that the slowdown of the adjoint model can be limited by the factor 7 for the assimilation window of 1000 time-steps while keeping only 27 instantaneous model states.

In this paper, we start from the analysis of the adjoint code obtained from the Tapenade software and removing from keeping an unnecessary data. So far, the time-stepping is performed by the leap-frog scheme, the adjoint code saves two time levels at each step: the nth and the n-1th. However, the n-1th layer is unnecessary to save because it was saved on the previous step. We avoid also saving the values of all variable on the continents because we do not need them at all. This divides the required memory approximately by two. We do not keep in memory all forcings: they have been read from files during the forward model run, they will also be read when necessary during the adjoint run. Despite the model is written with double precision, it reveals to be sufficient to keep the trajectory with the single precision only. Moreover, we can keep odd time-steps only, retrieving all even time-steps by interpolation. This interpolation perturbs only a little the gradient dividing the required memory by two.

Summarizing all the efforts, we can divide the memory, required for keeping the forward trajectory, by 25 ensuring that 10 days assimilation window of the ORCA2 model fits into 640 MB space. Taking into account possible use of the binomial checkpointing algorithm, we can state that automatic generation of the adjoint code by Tapenade can be used even for finer resolution models and for longer assimilation windows. In these cases we will be limited by a long computational time rather than by an excessive required memory. Massively parallel version of the adjoint code will, consequently, be necessary, but it is beyond the scope of this paper.

Discontinuous diffusion coefficients

Another problem we face trying to apply variational data assimilation to the ORCA-2 configuration consists in the discontinuous coefficients of the vertical diffusion. As we have seen in expressions ( 9), [START_REF] Janisková | Simplified and regular physical parameterizations for incremental four-dimensional variational assimilation[END_REF] and [START_REF] Kazantsev | Identification of optimal topography by variational data assimilation[END_REF], coefficients A z depend on the buoyancy ratio R ρ and on the Brunt-Vaisälä frequency N 2 in a discontinuous way. Infinitesimal variation of the Brunt-Vaisälä frequency from -ε to +ε in some grid point results in a drastic change of all diffusion coefficients A z from values of order 10 -5 ...10 -4 up to 100. Similar effect results from a small variations of the buoyancy ratio in the vicinity of 0.5 and 1. This implies the whole system is discontinuous and, hence, not differentiable. We can neither calculate the tangent model, nor its adjoint and the variational data assimilation can not be performed. This problem has already been reported in relation with the atmospheric model in [START_REF] Janisková | Simplified and regular physical parameterizations for incremental four-dimensional variational assimilation[END_REF].

So, we have to modify the diffusion coefficients in order to make the model differentiable. Unfortunately, a simple smoothing of the diffusion coefficients can not solve the problem. So far, the jumps of the coefficients values may exceed 10 6 , any reasonable smoothing modifies a lot low values of the coefficient in the adjacent points. Instead of values of order 10 -5 ...10 -4 we get 10 -2 ..10 -1 that dumps all the vertical physical processes.

To be able to perform variational data assimilation, we have to suppress completely the enhanced viscosity [START_REF] Janisková | Simplified and regular physical parameterizations for incremental four-dimensional variational assimilation[END_REF] and the double diffusive mixing [START_REF] Kazantsev | Identification of optimal topography by variational data assimilation[END_REF]. The turbulent closure scheme can, however, be smoothed replacing 8) by

A z v = max(A z 0 , C k l k √ ē) in (
A z v = C k l k √ ē + A z 0 -C k l k √ ē 1 + exp(5(C k l k √ ē -A z 0 )) (26)
Of course, this modifies the model solution. In particular, the suppression of the enhanced viscosity that dumps unstable situations with imaginary Brunt-Vaisälä frequency may lead to uncontrolled instability in the future. The question of the smooth diffusion coefficients that we can use in the variational data assimilation must be considered in details, but this long study is beyond the scope of this paper, where we focus our attention on the influence of the control parameters on the solution. In all the matter below, we shall use the model with smooth vertical diffusion coefficients as the basic model to be compared with models subjected to modified initial conditions or modified parametrizations of the boundary conditions. Controlling coefficients A z , we shall also start from the smooth ones comparing them to the obtained in the data assimilation process. fig. 3B, where the evolution of the distance ξ(t) [START_REF] Lellouche | Boundary control on burgers equation. a numerical approach[END_REF] is shown during the assimilation and 10 days after its end. The dotted line (the model with the optimal parametrization of the lateral boundary conditions) is almost indistinguishable from the solid line that corresponds to the original model.

This fact shows that lateral boundary conditions play no role in the ORCA-2 model. It is useless to control them hoping to improve the solution or to bring it closer to observations. This fact is in contradiction with the result obtained in [START_REF] Kazantsev | Boundary conditions control for a shallow-water model[END_REF][START_REF] Kazantsev | Sensitivity of a shallow-water model to parameters[END_REF] where the influence of the lateral boundary conditions on the solution of the shallow-water model is very important. But this contradiction can be evidently explained by different resolutions and different lateral dissipations. Indeed, the shallow water model discussed in [START_REF] Kazantsev | Boundary conditions control for a shallow-water model[END_REF][START_REF] Kazantsev | Sensitivity of a shallow-water model to parameters[END_REF] is an eddy resolving model with a developed turbulence while the grid step of ORCA-2 is approximately equal to 2 degrees (i.e. about 200 km) that can not be considered as eddyresolving at all. Moreover, the lateral dissipation coefficient used for the shallow water model was either 50 or 200 m 2 s , while in the ORCA-2 this parameter is equal to A xy u = 40000 m 2 s in extra-tropical region. Such a strong dissipation admits neither turbulence nor lateral boundary layers, i.e. configurations in which the lateral boundary conditions play an important role and worth to be controlled.

So far, the lateral boundary conditions are not important in this case, we shall not consider them as a controlled parameter below.

When we control either initial state φ 0 or vertical diffusion coefficients A z of the model we get similar results. Control of each parameter allows us to reduce the cost function value by more than 20% (two dashed lines with short dashes in fig. 3A). In the fig. 3B these two parameters show close behavior too (also two dashed lines with short dashes). The distance from observations is smaller than for the model with original parameters. That means the control of the initial state of the model and its vertical diffusion coefficients is important if we want to bring the trajectory closer to observations. We should note, that the trajectory remains closer to observations even after the end of assimilation. That means we can hope to get a better forecast with the optimal model. A. But the most influent parameter in this set of experiments is the parametrization of boundary conditions for vertical operators. If we control the conditions on the surface and on the bottom of the ocean, we get a rapid convergence of the cost function (dashed line with long dashes in fig. 3A) and the shortest distance to observations (the same line in fig. 3B). This indicates that vertical boundary conditions are very important in this case from the point of view of as the data reanalysis in the assimilation window and the forecast beyond the window.

Convergence of the

Modification of physical variables under control of vertical boundary conditions.

In order to examine the modification of the physical variables produced by the use of optimal model parameters, we plot first the sea surface elevation in two regions: North Atlantic and North Pacific. These regions are characterized by the presence of strong jet-streams (Gulf Stream in the North Atlantic and Kuroshio in the North Pacific) difficult to reproduce by low-resolution models like ORCA-2. The strength and the length of the jets are usually reduced by these models while the width is over-estimated. The sea surface elevation is the variable that reflects the integral impact of the jet on the model solution. The strength of the jet can be evaluated by the orthogonal derivative of the SSH while the length is represented by the length of the region with the big gradient.

A.

B.

C.

Figure 4: Sea surface elevation in the North Atlantic on the January, 30. Optimal initial conditions (A), Optimal Vertical boundary conditions (B) and Optimal vertical diffusion coefficients (C).

Three simulation have been performed for 30 days beginning from the 1st of January 2006. The sea surface elevation in the North Atlantic obtained as the final states are presented in the fig. 4. The first one (A) corresponds to the model that starts from optimal initial state. The second one is obtained with the optimal parametrization on boundary conditions of all vertical operators. Optimal coefficients for vertical diffusion were used to obtain the third result.

One can see that (A) and (C) parts in the fig. 4 look very similar. Optimal diffusion coefficients move the Gulf Stream slightly to the South increasing by 10 cm the positive anomaly on the South and reducing the negative anomaly on the North by the same 10 cm. The length and the strength of the Gulf Stream is the same in both pictures.

The modification of the circulation produced by the optimal parametrization of the vertical boundary conditions is much more visible in the fig. 4B. The positive anomaly is reinforced and moved to the North-East. The beginning of Gulf Stream is moved slightly to the North, becomes narrower and its force increased as well as its length. This case clearly represents much better solution from the physical point of view.

Similar effects can be observed in the North Pacific in fig. 5: slightly reinforced positive anomaly on the South of Kuroshio and reduced negative one on the North in (A) and (C) parts showing similar influence of the controlled parameters on the solution. And also, the fig. 5B obtained with the optimal boundary conditions for vertical operators is strongly modified. Like Gulf Stream, Kuroshio become longer, stronger and narrower. Moreover, it becomes also more tortuous that seems like a simulation of A.

B.

C. moving eddies. This is also a good result obtained in frames of a coarse resolution model.

The immediate idea one can have looking at these pictures that it is an artifact. Indeed, the sea surface elevation η is directly related to the vertical velocity on the surface (see [START_REF] Gilbert | Some numerical experiments with variable storage quasi-newton algorithms[END_REF]). We control w on the surface adding control coefficients α w s 0 and α w s 1 (see [START_REF] Dimet | A general formalism of variational analysis[END_REF]). So far, much observational data is collected by satellites and available for the sea surface elevation, the variable η might undergo a strong forcing pulling it toward observations. It is possible, consequently, that data assimilation modifies just coefficients α w s on the surface for w velocity and this modification is immediately translated to η with no influence on all other variables. However, analyzing the pattern4 of α w s we can not see any significant modification that can influence the jet streams. First, the magnitude of α w s 0 is of order 10 -7 while the velocity w on the surface is about 10 -5 ...10 -4 . Second, no characteristic patterns of jet streams can be observed in the fields α w s 0 and α w s 1 . And third, the perturbation of the vertical velocity by α w b 1 on the bottom is bigger than on the surface. Values of the velocity w in the vertical y -z section by the 60 0 W meridian are shown in fig. 6. This section is indicated by the thick vertical line in fig. 4B. One can see strong values of w, reaching 1.5 mm s , near the bottom of the ocean. These values can't be observed in the original model because of impermeability boundary condition imposed for w and no-slip condition for horizontal velocity components: u| bottom = v| bottom = 0.

Consequently, we shall state that it is the modification of the bottom boundary conditions that influences the jet streams near the surface.

An interesting question we may ask, whether it is the vertical velocity itself that is modified near the bottom, or the modification is due to divergent modification of lateral velocities u and v. In other words, whether coefficients α create fountains at the bottom that spray out jets of water, or they create three-dimensional eddies in which the divergence of lateral velocities ξ (6) is balanced by ∂w ∂z . To answer this question, we plot the u velocity in the same y -z section in fig. 7B (velocity u is orthogonal to the picture plane) and we see that the velocity u also has anomalies near the bottom in the same places as the velocity w, and namely in regions 18 0 ...21 0 N and 24 0 ...27 0 N. That indicates the source of the w velocity may lie in the divergent part of the lateral velocity components. Moreover, in Figure 6: Vertical velocity w obtained with optimal parametrization of vertical boundary conditions in y -z plane for x = 60 0 W.

A.

B.

Figure 7: Velocity component u in y -z plane for x = 60 0 W: obtained with optimal initial conditions (A) and optimal parametrization of vertical boundary conditions (B).

fig. 7B, one can see an important anomaly of u near the bottom (39 0 ...41 0 N) that is almost non divergent (the derivative ∂u ∂x is compensated by ∂v ∂y ) producing only a little anomaly in w. Moreover, in fig. 7B one can easily see the influence of anomalies of u and v near the bottom on the jet stream near the surface. The positive anomaly of u just below the jet stream (39 0 ...41 0 N) maintains the flux over the whole water column ensuring more than 10 cm s velocity at the 500 m depth. Comparing the velocity pattern with the velocity obtained in the experiment with optimal initial conditions presented in fig. 7A, we see that u can hardly reach the value of 5 cm s at the 500 m depth. That means the Gulf Stream become not only stronger (one can see the surface velocity exceeds 35 cm s instead of 25 cm s in fig. 7A) and longer, but also deeper.

This fact brings us to the conclusion that the major impact of the data assimilation is made on the parametrization bottom boundary conditions for lateral velocities. Consequently, the modification of the sea surface elevation in the jet-stream regions is not a simple artefact, but the result of optimization of the deep dynamics.

Finally, we should note that the temperature and salinity fields has been modified only a little in the data assimilation. We can see almost no difference with the original fields near the bottom, and a little difference near the surface, especially in the equatorial region.

As we have seen in fig. 3, optimal coefficients A z allows us to bring the model solution closer to observations. However, they don't amplify the stream jets fig. 4, fig. 5 like optimal boundary conditions. Analyzing the modifications made by the data assimilation, we can see numerous small regions (and even points) where the viscosity become vanishing. This result seems not to be physically interesting.

Of course, more pictures should be presented in this section in order to show in details the modification of different variables of the model. However, preserving the volume of the paper and taking into account the impossibility to present multimedia content in the Journal, the reader can find numerous color pictures and movies at http://www-ljk.imag.fr/membres/Kazantsev/orca2/index.html

Conclusion and perspectives

This paper shows possible advantages of extending of the now classical variational data assimilation procedure to the control of the model parameters other than initial conditions. Controlling parametrization of boundary conditions of vertical operators in the ORCA-2 model we have brought the solution closer to observations not only during the assimilation stage, but also during the forecast. We have shown also that the jet streams in the North Atlantic and the North Pacific are modified under control of the vertical boundary conditions. This result seems to be interesting because these jets are extremely difficult to simulate in frames of a low-resolution model. Another important particularity consists in the fact that it is the bottom boundary conditions for lateral velocities u and v were subjected to the major modifications. That shows the importance of the appropriate approximation of the bottom topography and processed in the bottom boundary layer.

Performing the control of parameters, it is not clear in advance what parameter is the most important for a given model. It was lateral boundary conditions in the case of the shallow-water model in [START_REF] Kazantsev | Boundary conditions control for a shallow-water model[END_REF][START_REF] Kazantsev | Sensitivity of a shallow-water model to parameters[END_REF], but they play no role in the ORCA-2 model. Taking into account the effort made to construct the vertical diffusion coefficients A z for this model (including an additional evolution equation for TKE ē(x, y, z, t)), one can suppose they might be extremely important in this configuration. However, their influence is not as strong as the influence of vertical boundary conditions.

This fact helps us to understand that boundary conditions are only important to optimize in the case when the corresponding operator is important for the model dynamics. The solution of the shallow-water model discussed in [START_REF] Kazantsev | Boundary conditions control for a shallow-water model[END_REF][START_REF] Kazantsev | Sensitivity of a shallow-water model to parameters[END_REF] was exhibiting either developed turbulence or strong under-resolved boundary currents showing that horizontal operators and their boundary conditions should be considered carefully. In the case of ORCA-2, almost all the horizontal dynamics is dumped by the strong dissipation that leads to a little influence of corresponding boundary conditions.

But, in the case of ORCA-2, the vertical dynamics is important. The dissipation and diffusion coefficients are chosen as low as possible (sometimes as low as A z = 10 -4 m 2 s compared to A xy = 4 × 10 +4 m 2 s ) in order to preserve the vertical currents. This fine adjustment points out the importance of the vertical dynamics and requires, consequently, careful treatment of its boundary condition. The present study confirms this fact, showing the optimizations of the vertical boundary conditions modifies much the currents.

Summarizing, we can say that the control of model parameters can help us

• to compensate the model errors and to bring the model solution closer to observational data improving the forecast;

• to see what parameters are more or less important for a given model and for a given configuration;

• to see what geographical region requires a particular attention in formulation of the boundary conditions.

The use of the automatic differentiation tool reveals to be extremely useful in this study helping us to avoid the huge coding and debugging work. This fact is very appreciated in the situation when we come out the frames of controlling the only classical model parameter -its initial state. Now, we can be brought to consider various parameters of the model as worthy to be optimized, facing the necessity to get the derivative of the model and its adjoint with respect to the chosen parameter.

The major shortcoming of automatic differentiation tool that consists in an excessive requirement of computer memory (hundreds or thousands times of the model code) can be avoided by the memory usage optimization and by the binomial checkpointing algorithm [START_REF] Griewank | Achieving logarithmic growth of temporal and spatial complexity in reverse automatic differentiation[END_REF]. These techniques bring the memory demand into a reasonable limit and the principal barrier to overcome becomes the usual computing time. Consequently, automatic differentiation of massively parallel codes and parallelization of the adjoints is necessary.

However, the present study can not be considered as a finalized control of the model parameters. Numerous points, that are extremely important from physical and mathematical points of view, have not been discussed in this paper.

Standing on the mathematical point of view, we must address two principal questions about uniqueness and stability of the parameter identification. These questions arise because we are solving an inverse problem, that is non-linear and probably ill-posed. The uniqueness of the solution determines whether it is at all possible to obtain the value of unknown parameter from observational data. The notion of stability determines whether small errors in data would cause serious errors in the identified parameter.

Both these questions have negative answer in this paper. The set of control coefficients was intentionally chosen to be too large to ensure an unique solution. As it has been shown in [START_REF] Kazantsev | Identification of an optimal boundary approximation by variational data assimilation[END_REF], exceedingly large set of α leads to a non-null kernel of the Hessian and results in a non-unique choice of optimal boundary conditions.

Moreover, the nonlinearity of the model solution with respect to boundary conditions is far from being quadratic. Consequently, the cost function I(p) [START_REF] Merryfield | A global ocean model with double-diffusive mixing[END_REF], along with the global minimum, may possess numerous local ones. Applying the minimization procedure, we shall find one of them instead of the global one that will result in a supplementary non-uniqueness of the solution. Incremental data assimilation technique may help us to improve the minimization by finding a deeper minimum, but it can not guarantee the minimum found is global.

And finally, we can not find a minimum at all because we can perform only a limited number (few) iterations of the minimization process in a reasonable computing time. Even with ORCA-2, one of the simplest models of the Global Ocean, we could perform some 40 iterations. Hence, the obtained result is far from even a local minimum.

Consequently, this study can not pretend to be an identification of the boundary conditions, but rather a way to compensate some model errors due to the lack of resolution and approximative parametrizations.

The same conclusion we can make from the physical point of view also. First of all, coefficients α are not physical parameters. They have been chosen as controls because any imaginable boundary condition can be approximated by a combination of these three α. Even unphysical conditions, like a permission to cross the continent contours, are intentionally accepted because they may point out that the model geometry is not in the agreement with the model dynamics. Moreover, we accept different α for different operators even dealing with the same variable assuming the possibility of different boundary conditions imposed to the same variable and increasing the probability to get a non-null kernel of the Hessian when one α compensates one another.

Thus, an exceedingly wide set of α is controlled instead of some physical set of boundary conditions in order to allow the data assimilation to modify any α bringing the model solution closer to observations.

Consequently, the values of α them-self are neither important, nor physically meaningful in this study. That's why they are not plotted in this paper. But they show the result we can potentially get optimizing boundary conditions, the regions where this optimization is particularly important and which operators and which variables of the model should attract a special attention. Analyzing the modification of the physical variables, we can further reduce the quantity of the controlled α and proceed to the control of the physical boundary conditions including the position of the boundary that may should be moved from the current position.

  S x wD z u Advection + D x A xy u ξ + D y A xy u x S z ρ(x, y, ζ)dζ Hydrostatic Pressure grad.

Figure 1 :

 1 Figure 1: Structure of the horizontal grid.

Figure 3 :

 3 Figure 3: Convergence of the cost function (27) during the minimization procedure (A) and the evolution of the distance ξ(t) (20) during and after the assimilation (B).

Figure 5 :

 5 Figure 5: Sea surface elevation in the North Pacific on the January, 30. Optimal initial conditions (A), Optimal Vertical boundary conditions (B) and Optimal vertical diffusion coefficients (C).

Table 1 :

 1 Controlled operators in x -y direction Total set of controlled α counts 2 000 808 coefficients. In the vertical direction, we control the discretization near the boundary of 4 derivatives and 5 interpolations (see table 2) accompanied by two approximations of the second derivative D zz u, v and D zz T, s according to[START_REF] Dimet | Variational algorithm for analysis and assimilation of meteorological observations. theoretical aspects[END_REF] and by the reconstruction of the vertical velocity w[START_REF] Dimet | A general formalism of variational analysis[END_REF]. The set of controlled coefficients counts 1 197 792 elements in this case.

	Argument			Operators
	T, s	D z u adv.	D z v adv.	S z T adv.	S z s	S z ρ
	w	D

adv. D y ω dissip. D y (vS y T, vS y s) D y ρ press.grad. S x S x u 2 Kin.en. S x v Coriolis S x s adv. S x (S y u) Coriolis S x T adv. S x ω adv. S y S y u Coriolis S y v 2 Kin.en. S y s adv. S y (S x v) Coriolis S y T adv. S y ω adv. z wT adv. D z ws adv. S z (wD z u) adv. S z (wD z v)

Table 2 :

 2 Controlled operators in z direction.

http://www.mercator-ocean.fr/eng/science/composantes-systemes/modelisation/orca2

http://www.nemo-ocean.eu/

http://www-tapenade.inria.fr:8080/tapenade/index.jsp

This picture is not presented. For more color pictures and movies, please see http://wwwljk.imag.fr/membres/Kazantsev/orca2/index.html

Acknowledgement: The author would like to express his gratitude to Arthur Vidard for providing the ORCA-2 configuration of the NEMO accompanied by forcings and observational data. All the contour pictures have been prepared by the Grid Analysis and Display System (GrADS) developed in the Centre for Ocean-Land-Atmosphere Interactions, Department of Meteorology, University of Maryland.

Assimilated data and background fields.

The observational data set for this model consists of the set of the ECMWF data representing sea surface anomalies issued from Jason-1 and Envisat altimetric missions and the set of vertical profiles of the Temperature and salinity from ENACT/ENSEMBLES data banque [START_REF] Ingleby | Quality control of ocean temperature and salinity profiles -Historical and real-time data[END_REF].

The sea level anomaly data set counts 112 thousands measured values distributed over 20 days interval (between the 1st and the 20th of January 2006). The distribution of the data is not uniform in time: we have about 6700 points per day during first 10 days and 4500 points per day after the 10th of January.

The temperature and salinity profiles counts about 200 000 measures during the same 20 days: more than 12000 observations per day before the the 10th of January, and about 6000 per day after this date. The quantity of data is not uniformly distributed within each day also. Some time steps of the model (1/15 day) receive more than 5000 data to assimilate, some other steps got 1500 data points only.

The background fields we use in the assimilation [START_REF] Leredde | On initial, boundary conditions and viscosity coefficient control for burgers' equation[END_REF] depend on the parameter under control. In this paper, the initial guess for the minimization and the background were always the same. Thus,

• the set of initial conditions of the model representing the model state on the January, 1, 2006, 00h GMT is used as initial guess and background when we control the initial point of the model

1 is used as initial guess and background when we controls the discretizations of horizontal operators (Table 1) in the vicinity of the continents;

1 that controls the discretizations of the vertical operators (Table 2) near the bottom and near the surface. This set is accompanied by α w 0 = 0, α w 1 = 1 used to calculate the vertical velocity w in [START_REF] Dimet | A general formalism of variational analysis[END_REF] and by α Dzz = 1 used to approximate the second vertical derivative in [START_REF] Dimet | Variational algorithm for analysis and assimilation of meteorological observations. theoretical aspects[END_REF];

• the instantaneous values of the vertical diffusion coefficients A z u , A z v , A z T , A z s on the 10th of January 0H GMT are used as initial guess and background when we look for optimal values of these coefficients.

Results

In this section we perform four experiments of data assimilation for identification of optimal values of four sets of parameters described above. All experiments are performed following the same algorithm and the same software. In each experiment one parameter is controlled, all other parameters are kept equal to their background states.

The data assimilation is performed by minimization of the total cost function

composed of ( 21) and ( 22) using the minimization procedure developed by Jean Charles Gilbert and Claude Lemarechal, INRIA [START_REF] Gilbert | Some numerical experiments with variable storage quasi-newton algorithms[END_REF]. This procedure uses the gradient of the cost function ( 23), [START_REF] Madec | Opa 8.1 ocean general circulation model reference manual[END_REF] in the the limited memory quasi-Newton method. We start from assimilation window T 5 days and allow the minimizer to make 20 calls of the routine that calculates the cost function and its gradient. Optimal values of the controlled parameters obtained as the final state of the minimization is used as initial guess in the second minimization, that is performed with the window T = 10 days and the minimizer is allowed to make 40 calls. The set of parameters obtained as this final result we shall call the optimal set.

Of course, this optimal set is optimal within the assimilation window only. However, we shall analyze the model behavior beyond the window, and namely 10 and 20 days later, on the 20th and 30th of January.

Comparison of different control parameters

Convergence of the cost function during the minimization procedure is shown in the fig. 3A. Despite 20 cost functions call were allowed in all experiments, the minimizer has made less iterations, i.e successful updates of the controlled parameter.

The worst convergence is obtained when the parametrization of the lateral boundary conditions is controlled (solid line in fig. 3A. Minimizer performs several iterations with both T = 5 and T = 10 days windows, but no visible reduction of the cost function is obtained. The same phenomenon we see in


Michaël Gabay

Gerd Finke

Nadia Brauner

Identical Coupled Task Scheduling Problem: The Finite Case

Keywords: Scheduling, Coupled tasks, Complexity, High multiplicity

Scheduling tasks on a radar system can be modeled as a multioperation scheduling problem where tasks are composed of two operations separated by a fixed duration. Tasks have to be scheduled on a single machine and the aim is to minimize the makespan. The complexity status of this problem has been settled for various cases but the one where all tasks are identical remains open. This paper aims at showing that this case raises original and difficult problems and is highly unlikely to be polynomial.

Introduction

A radar system usually has a single antenna which is both used to send and receive radio waves. Moreover, after the emission, a wave has to travel some distance, be reflected, and then travel back to the radar before it can be received and processed. Other operations can be scheduled on the radar system during the wave travel time.

In 1980, Shapiro [START_REF] Shapiro | Scheduling coupled tasks[END_REF] introduced coupled tasks to model radar operations. A coupled task i is a two-operations task. The operation lengths are a i (emit wave) and b i (receive wave) and these operations are separated by a fixed duration L i (travel time). A coupled task is illustrated Figure 1. The objective is to find a schedule that minimizes the makespan on a single machine (the radar). Shapiro [START_REF] Shapiro | Scheduling coupled tasks[END_REF] proved the N P-completeness of this problem and proposed 3 heuristics.

ai bi Li

Fig. 1. A single coupled task

Orman and Potts [START_REF] Orman | On the complexity of coupled-task scheduling[END_REF] studied this problem and settled its complexity for various cases (Table 1). However, they claim that the identical coupled task scheduling problem where ∀i, a i = a, L i = L, b i = b remains open. This special case differs from the other cases since the input consists of 4 integers: a, b, L, and the number of tasks n (w.l.o.g. assume a > b and a + b ≤ L). Hence the input length is O(log(L) + log(n)) and an algorithm that is polynomial in n would then be exponential. Therefore, the schedule is not a valid certificate. 1. Complexity of coupled task scheduling problems Such a problem is called a high multiplicity scheduling problem. Those problems discussed in [START_REF] Hochbaum | Minimizing the number of tardy job units under release time constraints[END_REF][START_REF] Hochbaum | Strongly polynomial algorithms for the high multiplicity scheduling problem[END_REF][START_REF] Clifford | Parallel machine scheduling with high multiplicity[END_REF][START_REF] Brauner | A framework for the complexity of high-multiplicity scheduling problems[END_REF] for instance.

Ahr et al. [START_REF] Ahr | An exact algorithm for scheduling identical coupled tasks[END_REF] proposed an exact algorithm using graphs for the identical coupled task scheduling problem. Its complexity is O(nr 2L ) where r ≤ a-1 √ a holds. Baptiste [START_REF] Baptiste | A note on scheduling identical coupled tasks in logarithmic time[END_REF] improved this result by proving that, for fixed a, b and L, the problem can be solved in O(log(n)) time. However, the constant is exponential in L and therefore this algorithm remains exponential.

Lehoux-Lebacque et al. [START_REF] Lehoux-Lebacque | Identical coupled task scheduling: polynomial complexity of the cyclic case[END_REF] solved the cyclic case. In this case, the aim is to maximize the throughput rate.

We worked on the finite case (n tasks) and showed that it is significantly different from the cyclic case. Especially, none of the properties crucial for the cyclic case hold. More details about the results exposed in this paper and proofs can be found in [START_REF] Gabay | Identical coupled task scheduling problem: the finite case[END_REF].

The Cyclic Case

Lehoux-Lebacque et al. [START_REF] Lehoux-Lebacque | Identical coupled task scheduling: polynomial complexity of the cyclic case[END_REF] proposed a polynomial time algorithm (complexity O(log(L) 2 )) to solve the cyclic case. They use a pattern to describe cycles and return a concise certificate. In order to prove their results, they point out and use the following properties that hold in the cyclic case:

1. We always place as many operations as possible on the separation time L; 2. The order of operations in a cycle does not matter; 3. The pattern of operations within each coupled task is an invariant; These properties allow to create a polynomial certificate in the length of the instance. This certificate is composed of 3 integers and its length is O(log(L)). An example of cycle is shown Figure 2. We started to work on the finite case with the idea that an optimal solution in this case may look like a cyclic solution surrounded by an extension to the left (to start the schedule) and an extension to the right (to finish it). We defined a pure strategy using the optimal cyclic pattern as the solution consisting to repeat this pattern (omitting the first b's as their a's were not scheduled) until n operations a have been scheduled and finish the last tasks. We proved the following theorem.

Theorem 1 (Asymptotic guarantee). The pure strategy using the optimal cyclic pattern is asymptotically optimal.

In fact, we proved that this strategy is a 1 + 2 k approximation where k is the number of times the cycle is repeated. Moreover, computing the makespan of this schedule can be done in polynomial time.

This raised the following question: is this pure-strategy-solution optimal for an n large enough ? We have shown that this is not the case: for some instances, for each value of n, there exists a larger n ′ such that the optimal makespan is strictly smaller than the pure-strategy solution.

Issues Raised

We have implemented the problem as an integer program and solved instances using the commercial solver CPLEX. Running extensive tests allowed to find counterexamples for all properties of the cyclic case. Moreover, some optimal solutions are really combinatorial and the optimal solutions of different instances are usually quite different. Because of these differences, we have not been able to identify patterns (which would have been a way to create a certificate) for solutions in the finite case. Moreover, even with fixed characteristics (same a, b and L), the optimal solutions may change a lot for different values of n. In concrete terms, none of the properties 1, 2, 3 in the cyclic case, section 2, remains valid. We constructed counterexamples for those three properties in the finite case.

Eventually, in some special cases we have been able to prove that weaker versions of some properties remain true (for instance: for some classes of instances, we always place as many operations as possible during the first idle time).

A Polynomial Case

We have worked on some special cases of this problem. When L a + 1 ≥ n, the greedy schedule is obviously optimal. We also found the optimal solutions for very specific values of n when a = b + 1. The result can be obtained and checked in polynomial time. In those cases, the optimal solutions are non trivial and require to use different patterns, connected by some irregular transitions as represented Figure 3.

Yet, all the other cases remain unsolved. There is in a way a lack of recursive properties for the optimal solution patterns (given a, b, L and increasing n).

Conclusion

We have shown that the finite case is significantly different from the cyclic case and that cyclic results and properties cannot be extended to the finite case. We demonstrated that a small subcase of the problem is polynomial but its solutions are very specific and elaborate. This strengthened the idea that the identical coupled task scheduling problem may not be polynomial and that we may not be able to find a valid certificate in the general case. The problem remains open but our results encourage us to think that it may not be in N P. It remains to prove this conjecture, for instance by reducing an EXPSPACE-complete problem to this one.

  Strongly N P-Hard aj ; Lj; bj aj = Lj = bj aj = bj = p; Lj aj = a; Lj = L; bj Open aj = a; Lj = L; bj = b Polynomial aj = Lj = p; bj aj = bj = p; Lj = L Table

Fig. 2 .

 2 Fig. 2. a = 5, b = 4, L = 15, an optimal cycle

Fig. 3 .

 3 Fig. 3. a = 10, b = 9, L = 53, n = 13, an optimal solution.


The model (Dynamics) ORCA2 Configuration : 182 × 149 × 31 nodes in curvilinear (x, y) coordinates with z levels. The model (Thermodynamics)

∂u ∂t = v(ω + f ) - ∂(u 2 + v 2 )/
∂T ∂t = - ∂uT ∂x - ∂vT ∂y - ∂wT ∂z Advection +A h T ∂ 2 T ∂x 2 + ∂ 2 T ∂y 2 + Hor.diffusion + ∂ ∂z A z T ∂T ∂z Vert.diffusion +Solar Radiation + Geothermal Heating + BBL + Surface ∂S ∂t = - ∂uS ∂x - ∂vS ∂y - ∂wS ∂z + A h T ∂ 2 S ∂x 2 + ∂ 2 S ∂y 2 + ∂ ∂z A z S ∂S ∂z + BBL + Surface A z T = Turbulent closure: ∼ max(A z 0 , C k l k √ ē) (1) ∂ē ∂t = A z v ∂u ∂z 2 + ∂v ∂z 2 -A z T N 2 + ∂ ∂z A z u ∂ē ∂z -cǫ ē3/2 lǫ ρ = ρ(T, S), N 2 = N 2 (T, S), R = R(T, S)
Eugene Kazantsev Boundary conditions control for ORCA2 page 3 of 21

∂u ∂t = SxSyv Sy(ω + f ) -Dx Sxu 2 + Syv 2 2 -Sz SxwDzu + DxA h u ξ + DyA h u ω + + g z 0 DxSzρ(x, y, ζ)dζ + Dzz(A z u u) + gDx(η + Tcφ) ∂T ∂t = -Dx(uSxT ) -Dy(vSyT ) -Dz(wSzT ) + A h T DxDxT + DyDyT + + Dzz(A z T T ) + Solar Radiation + Geothermal Heating + BBL ξ = Dxu + Dyv, ω = Dyu -Dxv, w = z H ξ(x, y, ζ)dζ; w(x, y, H) = 0
Interpolations and Derivatives (Sw

) k+1/2 = w k+1 + w k 2 k = 1, . . . , K -1 (DT ) k = T k+1/2 -T k-1/2 h k = 1, . . . , K -1
Eugene Kazantsev Boundary conditions control for ORCA2 page 4 of 21

∂u ∂t = SxSyv Sy(ω + f ) -Dx Sxu 2 + Syv 2 2 -Sz SxwDzu + DxA h u ξ + DyA h u ω + + g z 0 DxSzρ(x, y, ζ)dζ + Dzz(A z u u) + gDx(η + Tcφ) ∂T ∂t = -Dx(uSxT ) -Dy(vSyT ) -Dz(wSzT ) + A h T DxDxT + DyDyT + + Dzz(A z T T ) + Solar Radiation + Geothermal Heating + BBL ξ = Dxu + Dyv, ω = Dyu -Dxv, w = z H ξ(x, y, ζ)dζ; w(x, y, H) = 0
Interpolations and Derivatives Modified Near the boundary

(Sw) k+1/2 = w k+1 + w k 2 , k = 1, . . . , K -2, (Sw) 1/2 = α S 0 + α S 1 w 0 + α S 2 w 1 (DT ) k = T k+1/2 -T k-1/2 h , i = 2, . . . , K -2, (DT ) 1 = α D 0 + α D 1 T 1/2 + α D 2 T 3/2 h • • • w 0 w 1 w 2 w 3 w K w K-1 w K-2 w K-3 b T 1/2 b T 3/2 b T 5/2 b T 7/2 b T K-1/2 b T K-3/2 b T K-5/2

Eugene Kazantsev

Boundary conditions control for ORCA2 page 4 of 21

∂u ∂t = SxSyv Sy(ω + f ) -Dx Sxu 2 + Syv 2 2 -Sz SxwDzu + DxA h u ξ + DyA h u ω + + g z 0 DxSzρ(x, y, ζ)dζ + Dzz(A z u u) + gDx(η + Tcφ) ∂T ∂t = -Dx(uSxT ) -Dy(vSyT ) -Dz(wSzT ) + A h T DxDxT + DyDyT + + Dzz(A z T T ) + Solar Radiation + Geothermal Heating + BBL ξ = Dxu + Dyv, ω = Dyu -Dxv, w = z H ξ(x, y, ζ)dζ; w(x, y, H) = α 0 (x, y)
Vertical velocity

w i,j,K-1 = α w b 0 -α w b 1 hz i,j,K-1/2 ξ i,j,K-1/2 w i,j,k-1 = w i,j,k -hz i,j,k-1/2 ξ i,j,k-1/2 ∀k : 2 ≤ k ≤ K -1 w i,j,0 = w i,j,1 + α w s 0 -α w s 1 hz i,j,1/2 ξ i,j,1/2 • • • w 0 α w 0 α w 1 hz 1/2 w 1 w 2 w 3 w K +α w 0 α w 1 hz K-1/2 w K-1 w K-2 w K-3 Eugene Kazantsev
Boundary conditions control for ORCA2 page 4 of 21

u n+1 -u n-1 2dt = γ(u n+1 -2u n + u n-1 ) Asselin filter -SxSyv n Sy(ω n + f ) - -Dx Sx(u n ) 2 + Sy(v n ) 2 2 -Sz Sxw n Dzu n + + DxA h u ξ n-1 + DyA h u ω n-1 Explicit Euler + DzzA z u u n+1 Impl.Euler + + g z 0 DxSzρ n (x, y, ζ)dζ + gDx(η n+1 + Tcφ n+1 ) Implicit T n+1 -T n-1 2dt = γ(u n+1 -2u n + u n-1 ) Asselin filter -Dx(u n SxT n ) -Dy(v n SyT n ) - -Dz(w n SzT n ) + A h T DxxT n-1 + DyyT n-1 Explicit Euler + + Dzz(A z T )T n+1

Adjoint

The models solution depend on initial conditions and a number of parameters:

∂T ∂t = -Dx(uSxT )-Dy(vSyT )-D (α) z (wS (α) z T )+A h T DxxT +DyyT +Dzz (α) (A z T )T
The model x(t) = M 0,t (x 0 , α)

We calculate the derivatives and their adjoints with respect to Modified formula The control may violate the model physics;

w i,j,K-1 = α w b 0 -α w b 2 hz i,j,K-1/2 ξ i,j,K-1/2 w i,j,k-1 = w i,j,k -hz i,j,k-1/2 ξ i,j,k-1/2 ∀k : 1 ≤ k ≤ K -2 w i,j,0 = w i,j,1 + α w s 0 -α w s 2 hz i,
The physical meaning of α is difficult to understand;

The set of α is not unique;

The problem of identifiability is not addressed yet;

The problem of stability is not even posed.

Consequently:

It is not a parameter estimation study, but a way to compensate model errors showing the most influent parameter (vertical BC for ORCA-2, lateral BC for Shallow-Water).

Eugene Kazantsev Boundary conditions control for ORCA2 page 20 of 21

  Total reduction of required memory is up to 25 times. 10 hours window =⇒ 10 days window.

		TAPENADE 3.6 (Tropics team, INRIA)
		with the Memory usage optimization:
	search for push/pop CALL PUSHREAL8ARRAY(sold, nx*ny*nz) CALL PUSHREAL8ARRAY(told, nx*ny*nz) SSH, North Atlantic, January,1-30 2006. SSH, North Pacific, January,1-30 2006.	replace by call push uvts(u,v,t,s,ssh)
	CALL PUSHREAL8ARRAY(vold, nx*ny*nz)
	CALL PUSHREAL8ARRAY(uold, nx*ny*nz)
	CALL PUSHREAL8ARRAY(ssh, nx*ny)
	CALL PUSHREAL8ARRAY(s, nx*ny*nz)
	CALL PUSHREAL8ARRAY(t, nx*ny*nz)
	CALL PUSHREAL8ARRAY(v, nx*ny*nz)
	CALL PUSHREAL8ARRAY(u, nx*ny*nz)
	Procedure push/pop uvts(u,v,t,s,ssh):
	does not push n -1 step and pops appropriate values (divides the required memory by 2)
	does not push u, v, t, s in lower level routines
	does not push values on continents (divides by 2) Optimal IC Optimal BCz
	pushes values in Real*4 format (divides by 2) eventually pushes only odd timesteps and interpolate when poping (divides Optimal IC Optimal BCz
	by 2)
	Eugene Kazantsev Eugene Kazantsev Eugene Kazantsev	Boundary conditions control for ORCA2 Boundary conditions control for ORCA2 Boundary conditions control for ORCA2	page 6 of 21 page 6 of 21 page 11 of 21

x 0 , α by

TAPENADE 3.6 (Tropics team, INRIA)

. that allows us to avoid a HUGE development/coding (a double of the classical one, at least) to obtain immediately the derivative with respect to any parameter we want.

Eugene Kazantsev

Boundary conditions control for ORCA2 page 12 of 21

Data ECMWF data issued from Jason-1 and Envisat altimetric missions and ENACT/ENSEMBLES data banque. January, 1, 2006. Difference between observations and background during the 1st of January. Eugene Kazantsev Boundary conditions control for ORCA2 page 7 of 21 Data January, 1-20 2006. Number of observation per time step during the 1 -20 Jan.2006 Probability density function for the difference (observation-background).

Eugene Kazantsev Boundary conditions control for ORCA2 page 8 of 21

Vertical diffusion

Eugene Kazantsev

Boundary conditions control for ORCA2 page 4 of 21

Cost function

The model:


Jean -Paul Douzals
email: jean-paul.douzals@irstea.fr

Antoine Porte

Pierre Fernandez

Simulating CoV from Nozzles Spray Distribution: a necessity to investigate spray distribution quality with drift reducing surfactants

Keywords:

This paper aims at showing possibilities of spraying quality control through nozzle flow distribution (CoV), especially when a drift reducing surfactant is used. The effect of nozzle height combined with hydraulic pressure was achieved considering a homothetic translation of flow rate. In all cases of FF 110 02 (standard or air injection), the evolution of CoV with nozzle height shows periodic pattern where the effect of surfactant concentration revealed critical CoV variations. As a result, the use of such surfactant shall necessarily be preceded by flow distribution evaluation in order to adjust the boom height adequately. The use of a nozzle simulator may greatly help to develop boom adjustment guidelines.

Introduction

Flow distribution under a spray boom is one of the major macroscopic parameter of spraying quality. It is dependant on several parameters such as spray angle and the recovery between nozzles mainly depending on the boom height [START_REF] Debouche | Simulated and measured coefficients of variation for the spray distribution under a static spray boom[END_REF]. Transversal spray distribution was widely studied in the past due to the effect of boom movements (mainly yaw and roll) [START_REF] Chapple | Assessment of single-nozzle patternation and extrapolation to moving booms[END_REF], [START_REF] Lardoux | Test method for Boom Suspension Influence on spray distribution, Part I: Experimental study of Pesticide Application under a moving boom[END_REF]. Original work was also done on the capability of a single nozzle spray pattern to determine worn nozzles [START_REF] Maertens | Optimization of sequence and orientation for used nozzles based on few, full boom distribution measurements[END_REF]. A comprehensive literature is existing on the effects of surfactants or coformulants, especially those aiming at reducing spray drift, in terms of their impact on droplet sizes [START_REF] Miller | Effects of formulation on spray nozzle performance for applications from ground-based boom sprayers[END_REF], [START_REF] Stainier | Droplet size spectra and drift effect on two phenmediphanm formulations and four adjuvant mixtures[END_REF]. However, very few researches were dedicated to the effect of chemicals on spray angle and spraying quality. The simulation of flow distribution with regards to boom height and working pressure is studied in order to recover acceptable CoV accordingly to [EN 12761-2].

Materials and methods.

Flow distribution of a single nozzle placed at 50 cm height is measured with a 5 cm channel test bench compatible with ISO 5682-1. The maximum capacity of patternator is 3 m (60 channels). Flow rate collected by each channel is driven to a tube mounted on individual weighting cell. When the volume collected in the tube reaches 400 ml, all tubes are moved apart form channels, weighted and emptied. Each measurement consists of 3 repetitions at the same position. In case of flow distribution measurement under a complete boom of 10 m (20 nozzles), the patternator is moved by steps of 50 cm with 3 repetitions by step. FIGURE 1: Single nozzle replicationboom patterngeneric pattern FF02 2.5 bar A generic pattern of 50 cm length is calculated form the overlapping from beside nozzles (Fig. 1). As all nozzles are identical, the distribution Cov can be calculated directly on the generic pattern (10 values) with the following formula: CoV = standard deviation/mean. The effect of nozzle height is then simulated assuming a homothetic translation of distances with nozzle height and an inverse translation of flow rates (Fig. 2). FIGURE 2: Homothetic translation of flow rates distribution with nozzle height FF02 2.5 bar A flat fan 110 02 nozzle was tested at 30 up to 90 cm height and pressures of 2.5 and 5 bar with water. A model mixture with a highly effective surfactant was tested for concentrations from 0.03% to 0.12%.

Results

Feasibility of homothety.

CoV of 5 different nozzles of the same type (FF 110 02) were measured depending on the nozzle heights (40 -50 -60 -70 cm) for 5 pressures from 1.5 -2 -2.5 -3 -3.5 bar. Simulated CoV were calculated from values at 50 cm and simulated for other heights with a homothetic translation with acceptable correlation R² of 0.87. The simulation showed quite similar results for high values of nozzle heights and hydraulic pressure as results were sometimes less uncertain for low pressure and low heights. It is to say that the spatial resolution of the test bench is 5 cm. In order to get more accurate results, the simulation process needed values every cm calculated by linear interpolation.

3.2.

Evolution of CoV according to nozzle height and hydraulic pressure. 3.2.1. Pure water (table 1). The same evolution of CoV vs. nozzle height was found for water and surfactant mixture with a semi-periodical curve but the amplitude of CoV was much greater in the case of the surfactant with values up to 45%. As a practical consequence, the prediction of the negative impact of surfactant on the spray distribution and the benefit of optimizing boom height and working pressures are possible. In all cases, lowest CoV were obtained for heights around 37 cm or 80 cm as worse values of CoV were visible at 50 cm height. Practical consequences of this study rely on the necessity of the boom height adjustment in order to maintain an acceptable spraying quality.

FF In the case of an air injection nozzle (FF IA 02 2.5 bar), results were similar to the reference nozzle FF 02 for pure water (Fig. 4). However the combination of surfactant with the air injection nozzle amplified the effects on reduced spray angles as well as CoV degradation (Table 2).

FIGURE 4: Simulated CoV for FF AI 02 @ 2.5 bar with pure water

Conclusion

This paper aims at showing the effect of surfactant C1 on spray angle and the consequences in terms of spray quality degradation. This illustrates the notion of sprayability of surfactants and chemicals in general. With the help of a distribution CoV simulator based on homothetic translation of flow rate and distances with nozzle height, it is possible to predict a wide range of realistic situations and possible curative adjustment of boom heights. Moreover, this work shows that a correct flow distribution is required before any drift evaluation in the field or in a wind tunnel.

TABLE 1 :

 1 Combined effect of pressure and nozzle height on CoV values FF 02.One can see on Table1the effect of nozzle height on the periodic variation of CoV for a given pressure. The general trend is a decrease of CoV while nozzle height increases, showing the influence of overlapping in the quality of spraying. Table1shows also the cooperative effect of pressure with nozzle height. Increasing pressure allows acceptable CoV for nozzle heights from 30 to 90 cm. Improper combinations of pressure and height may involve less satisfying conditions. surfactant was used at different concentrations (0.03% to 0.12%) to simulate the effect of CoV. Single distributions at 50 cm height revealed immediate shift in the CoV values as the spray angle (Table2) was strongly narrowed, proportionally to the concentration of surfactant. A first measurement at 50 cm height showed CoV varying from 20 to 40%. The evolution of simulated CoV vs. nozzle height is shown on fig.3.

				Pressure
	Nozzle Height (cm)1.5 bar 2 bar	2.5 bar 3 bar	3.5 bar
	30	13.03% 11.04%	9.56%	8.84% 8.46%
	32.5	10.56%	7.93%	5.33%	4.14% 4.39%
	35	11.99%	6.43%	3.43%	2.72% 2.76%
	37.5	11.56%	5.62%	3.45%	3.89% 3.35%
	40	11.67%	4.67%	2.88%	3.77% 3.39%
	42.5	12.11%	4.60%	2.91%	3.51% 3.51%
	45	11.30%	3.62%	3.04%	3.86% 4.74%
	47.5	9.89%	3.20%	4.19%	4.95% 5.66%
	50	7.81%	4.42%	5.98%	6.55% 6.68%
	52.5	4.01%	5.92%	6.96%	6.90% 6.46%
	55	2.93%	7.10%	7.31%	6.94% 5.85%
	57.5	3.71%	8.15%	7.34%	6.71% 5.15%
	60	5.48%	8.19%	6.60%	5.74% 3.87%
	62.5	7.64%	8.01%	5.68%	4.46% 2.51%
	65	8.92%	6.97%	4.28%	2.99% 1.60%
	67.5	9.47%	5.71%	2.97%	2.13% 1.61%
	70	9.89%	4.57%	2.00%	2.18% 2.19%
	72.5	9.33%	3.25%	1.50%	2.76% 2.46%
	75	8.32%	2.24%	1.94%	3.49% 2.68%
	77.5	7.20%	1.68%	2.31%	3.67% 2.40%
	80	4.93%	2.11%	3.20%	4.24% 2.65%
	82.5	3.93%	2.20%	2.76%	3.47% 1.72%
	85	2.47%	2.43%	2.51%	2.94% 1.47%
	87.5	1.36%	2.97%	2.47%	2.63% 1.17%
	90	1.82%	2.90%	2.03%	2.10% 1.08%

TABLE 2 :

 2 Evolution of spray angle for 99,98% of total flow rate for FF 02 standard and Air Injection Height 50 cm -Resolution in angle definition is +/-4° FIGURE 3: Distribution CoV for water and surfactant as a function of nozzle height: 50cm measuredother heights ::simulated.

  Ageng 2012. International Conference on Agricultural Engineering, 8-12 Juillet 2012, Valence, Espagne 3.2.3. Case of Air injection nozzle.

							110 02			FF 110 02 AI
						Spray angle (°)	Spray angle (°)
		Water				94,7				103.2
		C1 0.03%				80				n.d.
		C1 0.06%				70.6				79.1
		C1 0.08%				66.9				75.7
		C1 0.1%				63				72.1
		C1 0.12%				63				64.5
	55%
	50%
														Water
	45%													C1 -0,03%
														C1 -0,06%
	40%													C1 -0,08%
														C1 -0,1%
	35%													C1 -0,12%
	30%
	CoV %
	25%
	20%
	15%
	10%
	5%
	0%
	25	30	35	40	45	50	55	60	65	70	75	80	85	90	95
							Nozzle Height (cm)

CemOA : archive ouverte d'Irstea / Cemagref CIGR-

  Ageng 2012. International Conference on Agricultural Engineering, 8-12 Juillet 2012, Valence, Espagne

		12.00%
		10.00%
		8.00%
	CoV	6.00%
		4.00%
		2.00%
		0.00%
		0	10	20	30	40	50	60	70	80	90	100
						Nozzle height (cm)

CemOA : archive ouverte d'Irstea / Cemagref CIGR-


CADP: A Toolbox for the Construction and Analysis of Distributed Processes

Hubert Garavel, Frédéric Lang, Radu Mateescu, Gwen Salaün, and Wendelin Serwe Inria Grenoble -Rhône-Alpes / Laboratoire d'Informatique de Grenoble 655, avenue de l'Europe, Inovallée Montbonnot, 38334 St Ismier Cedex, France Description Asynchronous concurrency is becoming increasingly present in a large spectrum of systems, spanning from the level of systems-and networks-on-chip, over multi-processor architectures, up to the level of grid and cloud computing. Due to the intrinsic complexity of asynchronous concurrency, the correct design of such systems is notoriously difficult, requiring the support of formal methods and verification tools.

CADP (Construction and Analysis of Distributed Processes) [START_REF] Garavel | CADP 2010: A Toolbox for the Construction and Analysis of Distributed Processes[END_REF] is a toolbox for the design, functional verification, and performance evaluation of asynchronous concurrent systems. Currently, CADP consists of about fifty interconnected tools and software libraries1 . The toolbox is distributed free of charge to academia and public research institutes2 , and is already used by more than 440 research institutions and companies worldwide in many application domains. Given the increasing number of systems featuring asynchronous concurrency, CADP could be used still more widely in research, industry, and education (in particular for teaching the concepts of concurrency theory).

Tutorial Objectives

This tutorial presents the architecture and main functionalities of CADP, with a twofold objective. On the one hand, the tutorial illustrates the application of CADP to the modeling, functional verification, and performance evaluation. On the other hand, the tutorial presents various input languages accepted as input by CADP, together with software libraries that enable users to develop their own analysis tools.

The well-known, but fundamental problem of mutual exclusion [START_REF] Mateescu | Model Checking and Performance Evaluation with CADP Illustrated on Shared-Memory Mutual Exclusion Protocols[END_REF] will serve as support to illustrate the principal functionalities of CADP: formal modeling of protocols, compositional state space generation, graph visualization, interactive step-by-step simulation, formulation and verification of temporal logic properties, as well as performance evaluation by compositional insertion of latency constraints and transformation into interactive Markov chains [START_REF] Hermanns | Interactive Markov Chains and the Quest for Quantified Quality[END_REF][START_REF] Coste | Ten Years of Performance Evaluation for Concurrent Systems Using CADP[END_REF].

Covered Topics

• Theoretical foundations: labeled transition systems (LTS), bisimulations, Boolean equation systems, etc.

See http://cadp.inria.fr.

Since July 2011, a CADP license is granted to any user who can provide an academic e-mail address and Web page, without the need to sign a paper contract.

• The LNT language for the formal modeling of asynchronous concurrent systems.

• State space exploration tools for step-by-step simulation, state space generation, graph visualization, etc.

• The MCL language for the expression of temporal logic properties extended with data.

• Verification tools for model checking, equivalence checking, etc.

• Models and tools for performance evaluation: interactive Markov chains, numerical analysis, steady-state simulation, etc.

Presenter

Radu Mateescu holds a PhD and an HDR (Habilitation à Diriger les Recherches) from Grenoble INP (National Polytechnic Institute). He is currently a senior researcher at Inria and leader of the joint Inria/LIG Convecs team (http://convecs.inria.fr), which is a follow-up of the Vasy team.

In the past, he gave courses on formal methods at ENSIMAG (Grenoble), ESIREM (Dijon), ESIA (Annecy), and the VTSA summer school of the Max Planck Institute (Saarbrücken, Germany).

Duration: half a day (about 3 hours).

Keywords: asynchronous concurrency, Boolean equation system, bisimulation, compositional verification, distributed systems, exhaustive state-space exploration, formal modeling languages, formal specification, interleaving semantics, Markov chains, modal µ-calculus, on-the-fly verification, performance evaluation, process algebra, temporal logic.