The manuscripts describes a C-grid implementation of CICE. It is a well written paper, although I would question the advantages of the C-grid as compared to the variational implementation on the B-grid from purely numerical side -- the implementation of rheology is of course possible on a quadrilateral C-grid, as demonstrated previously and in the present manuscript, but is not optimal in  the sense that diagonal and off-diagonal components of the strain rate are naturally placed at T and U locations, while the field of Delta  needs all them at the same location. Additionally, the variational (in essence bilinear finite-element) B-grid implementation is more accurate. Transport in narrow channels (not coupling, which is trivial) is a valid argument, but what is the accuracy of the divergence of internal stresses in such channels in a general case?

Leaving this discussion aside, I certainly recommend the manuscript as it be helpful for those users who prefer the C-grid implementation. The central place in the manuscript is devoted to the description of the adjustments of the incremental remapping methods used in CICE (EFA method), which is of obvious interest. It is, however, mentioned that this method was already available in CICE, which immediately leads to the question what is new? I would recommend to state more clearly what is the new part, and what is well known. 

I do not think the analysis presented in appendix B is relevant as it starts from a problematic assumption: for u=u_0+u'=u' as taken further in the analysis, (B1) is not reduced to sigma=-P. Furthermore, the manuscript deals with the EVP method, and uses the replacement pressure, which is not reflected in the analysis, but is crucial for the behavior of perturbations. The answer is simpler. Generally one is not allowed to interpolate transport velocities, as the discrete divergence is only consistent in a certain sense. The Fourier symbol of C-grid divergence is d=i*[u*sin(ka/2)+v*sin(la/2)](2/a), where a is the cell side. The u and v velocities are not arbitrary, but given by the dynamics. Averaging to U-points and computing the divergence after will introduce (1+cos(la))/2 with u and (1+cos(ka))/2 with v, i.e. will create a spurious additional divergence d_sp=i*[u*sin(ka/2)*(cos(la)-1)+v*sin(la/2)*(cos(ka)-1)](1/a), which is simply a source of spurious modulation of scalar fields in (B7). The divergence d in the experiment shown in Fig. 3(b) has a large scale structure, with the implication that the u and v contributions in d nearly cancel for k,l approaching to grid cutoff wavenumbers. But this would imply a substantial d_sp at grid scales. This is a common issue, to avoid it one reconstructs scalars at faces or in control volumes, but not the advecting velocities. The EFA method restores the consistent divergence and removes spurious sources. 

Minor things:

line 22  'The C-grid ...' -- but the B-grid may have some other advantages. I would prefer to see a more critical analysis here.

line 31 - 40. I would say here that the manuscript presents an EFA version of the incremental remapping which is crucial for C-grid velocities. The history should be in the text, otherwise the story is told twice.

line 50 ... (3) is hardly the main contribution -- it is an element of history, the main contribution is the statement that the EFA method has to be used; As mentioned, I do not think the analysis (4) is rigorous enough to keep it in its present form in this manuscript.   

line 99 Why is E_0 less than 1? In the limit dt->0 this would require to increase the number  of substeps.

       

line 205 ' mathematical mode' --> numerical mode 

Figure 5: In a way this illustration shows that the C-remap in not really an accurate scheme (although better than the first order upwind). How is it related to a slope or flux limited advection (e.g. with superbee)? Which accuracy is needed?

line 325 Could the reason be that 'effective' diffusion associated with the advection is sufficiently high on coarse meshes? 

line 360  I think that the point is not in the null space created by averaging but in spurious divergence introduced by this averaging. 

line 450 Formulas below are written under an assumption that xi_1 and xi_2 parameterize a flat manifold. In reality, a covariant derivative of velocities will contain a vertical contribution on the  sphere (try to compute a vector Laplacian in the same way as A13-A16 to see that some terms are lost)

line 611 For the velocity field used further D_d is an oscillating function, which makes this statement contradictory.

753 Why 'in CICE'? This is a general expression.

The paper describes the implementation of a C-grid discretization in CICE. It explains how the C-grid discretization can be coupled with an incremental mapping scheme originally developed for a B-grid. The authors follow the approach of Kimmritz 2018 and Bouillon 2013 in implementing the C-grid for the momentum equation. The main part of the paper describes how the existing incremental remapping scheme in CICE can be modified to work with the C-grid approximation of the momentum equation. To my understanding, the coupling (modification) of the remapping is the only scientifically new development/method that is not documented in another paper. However, the coupling presented is not discussed in relation to existing approaches in the field/literature. It seems inconsistent to me to implement a C-grid for the momentum equation, but adapting a B-grid transport (instead of using C-grid schemes e.g., the second order DG, which is more accurate). Furthermore, the proposed EFA scheme has problems on non-uniform grids.

Figure 5. shows the transport in a one cell wide channel and compares the C-upwind to the C-remap. Both approaches are diffusing. There are also DG type transport schemes which work in one cell wide channel.  These methods are able to preserve the shape e.g. The neXtSIM-DG dynamical core: A Framework for Higher-order Finite Element Sea Ice Modeling.I think this should be discussed and acknowledged in the paper.  

The authors state several times that the most difficult part is the approximation of the rheology term. But how well are the stresses and the strain rates approximated compared to the existing B-grid approach in CICE? To address the issue, I would solve a 2D equation

Div( nabla v+nabla v^ T)=f,

calculate the analytical solution for this and check under mesh refinement how the error convergences for both approaches, C-grid and B-gird. When it comes to the fully coupled system how large the coupling error?  Can you represent LKFs with the same quality with the C-grid?

Minor things:

Introduction: I miss a general overview of how other models with C-grid deal with the advection and how they couple to the advection. For example, the MITgcm. Is your proposed coupling more or less accurate?

Introduction: Could you give more details on the EFA method.  Who developed the EFA method? Is it used in other research fields?

L.144 What is the difference to Bouillon and Kimmritz? If there is no important difference why do repeat it in the appendix?

L 213 “ Fortunately”   It sounds like you are happy to avoid programming. I would remove the word.

166 How much accuracy do you lose by the averaging?

257 How does the EFA method ensure that the divergence is consistent with the divergence calculated by the dynamical solver?

343 Here are to many dots

381-3.85 I think it is trivial that this is fulfilled since your using C-grid for the momentum equation that has already been used by other models.