This papers proposes an objective method to tune ESMs in a global approach. The authors proposes a mathematical equivalent to the way modellers tune their models so as to rationalize and automatize the tuning process, based on a Green’s functions approach. This is a very important topic for ESM developers that come to a limitation in the model developments due to the long process of human-made iterative tuning. This study is perfectly timely. The authors provides an example of application and discuss some of the assumptions made. However, some assumptions are not fully explained or would deserve more comments. I therefore recommend the manuscript to be accepted after minor revisions.

Main comments

Section 3: Explaining the methodology using a mathematical formalism is a very good idea, but in the present form, I have not the impression that it is sufficiently clear. Maybe I am not a good mathematician but I think that this manuscript should be well understood by modellers and improving the explanation would leverage the impact of the study. Moreover, the equations provided in the manuscript are already presented in Menemenlis et al. (2005), therefore, they do not provide anew outlook. I would suggest to better explain these equations here. It is not clear to me how equation (4) is established? Similarly, it is not clear to me what is the difference between G(ƞ) and G ƞ? Maybe the notations could be simplified. Additionally, the physical meaning of Q and R should be detailed as well as the methodology to set them. For instance in section 4.3, it is said that “R is chosen diagonal”: apart from being mathematically convenient, what are the physical hypothesis behind such a choice? What makes this choice sensible? Isn’t it a strong hypothesis? What happens if reality is far from diagonal? These points should be discussed with a physical point of view in the manuscript to better illustrate the advantages and flaws of the method used.

Similarly, section 5.1 is a bit hard to read. The physical idea behind the choice of functions is not clear enough. What is the physical meaning of Q choices proposed? This section need to be detailed to better catch the physical hypothesis that are made depending on each choice of covariance matrix.

There is also a crucial point that should be at least discussed in the paper: in such complex ESM there is a long spin-up process of several decades that has to be taken in to account. If only taking the first 10 years of ESM simulations for tuning, what is the result of such a sort term tuning on the long term equilibrium? How do the authors proposes to handle this problem?

Line 230: the fact that the surface temperature is less well reproduced after the tuning process is very interesting. Even if it reflects a flaw in the present choice of observational targets, the discussion on this should be enriched. What type of other observational targets could be proposed to avoid this?

Section 5.2: by increasing the number of parameters, there is a risk to over-fit the model to the set of observation chosen as discussed in Dommenget and Rezny, (2018). This risk should be discussed in the last section. Additionally, by increasing the number of parameters the final fit is necessary improved as more degree of freedom have been used. Do the authors could use a methods similar to the LASSO technique (Tibshirani, 1996, used in linear regression) to penalize the results according to the number of degree of freedom? This would help selecting only the more useful parameters and avoid over-fitting.

Ln 304-305: ”The largest cost reduction was found between 5 to 10 optimized parameters.” And what about the ratio projected cost/optimized cost? Maybe this information could be added on figure 6. Is this ratio stable depending on the number of parameters?

Section 5.4: what would be the cost if the first 5 years were used? Would it be different from the last five years? If yes, this would probably says that the spin-up is important. But then, what about taking the last 10 years over a 20 years simulation. The initial choice of 10 years appears as a bit short against the spin-up timescale of coupled models. The spin-up problem would be discussed more thoroughly. Maybe it would be nice to show the model initial drift.

I may have misunderstood the methodology, but reading again the introduction section, I am not convinced that the proposed methodology addresses all the points raised in the introduction.

- ln 49-52 :“An experiment with both perturbed would not reveal which had the greater impact on the solution, and experiments with only one at a time perturbed would not behave as an experiment with both perturbed. The sequential methodology, therefore, for this example, may result in a sub-optimal combination of these parameters.” I do not understand how the Green’s function approach do something different to “one at a time” perturbed experiment. Maybe it is just a question of formulation here?

- ln 58 :”The second drawback is the large number of needed simulations, since a new set of experiments is required for each set of observational targets”. Using the standard by hand method, developers usually have a look at a panel of diagnostics and do not need to redo the experiments for a new diagnostic, I do not see any difference on this point. This point is repeated on line 74, I would not emphasize this much.

- ln 71-71: “The interdependence is accounted for in the constraints on the error covariance for each parameter.” This points of interdependence is not well developed in the manuscript. I would promote to better illustrate this point.

Miscellaneous

Section 2: please indicate the ocean model time-step and better separate the ocean and atmosphere models description in two distinct paragraphs.

Table 4 and 5: these tables are large and difficult to interpret as is, the exact values not behind of strong interest. Please consider replacing them by plots.

Section 5.2: Concerning LW and SW radiation, is the change from downward to net? Why changing this point? Could the authors discuss a bit more this choice?

References :

Dommenget, D. and Rezny, M.: A Caveat Note on Tuning in the Development of Coupled Climate Models, J. Adv. Model. Earth Syst., 10, 78–97, https://doi.org/10.1002/2017MS000947, 2018.

Tibshirani, R.: Regression Shrinkage and Selection Via the Lasso, J. R. Stat. Soc. Ser. B Methodol., 58, 267–288, https://doi.org/10.1111/j.2517-6161.1996.tb02080.x, 1996.

In this work the authors used Green's function approach to optimize parameters in an earth system model. This optimization was achieved by minimizing a weighted least-squares distance. The numerical results demonstrated that the Green's functions estiamtion approach can be a good fine-tuning step in the model development process to adjust uncertain parameters of the model considered in this work. Also, sensitvity analysis results were presented to further investigate the effectiveness of the Greens's function approach. In general, as an appliction work, the manuscript is well organized and technical details are clear. 

I have a few comments on the approach itself although this work is focused on application.

1.  Eq(5) is the combination of Eq (1) and Eq (3), and G seems like a composite of M_i and H. How do we understand it as an "convolution of observation operator H and M_i"?

2. Eq(4) and Eq (9) are similar to some formulas in Kalman filter. Indeed, Eqs (2) and (3) are also the starting point of Kalman filter. What is the main difference between the Green's function method and Kalman filter? Or are they equivalent to each other to some extent or in some cases?

3. The impact of choosing different Q was discussed using Q=\infty and Q=I. How about the influence of  R? Here R is a diagonal matrix, which indicates that the errors for different observables are independent (which is a quite standard setup in Kalman filter). How about the errors are not independent? Any comments or intuitions?