the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A simplified atmospheric boundary layer model for an improved representation of air–sea interactions in eddying oceanic models: implementation and first evaluation in NEMO (4.0)
Florian Lemarié
Guillaume Samson
JeanLuc Redelsperger
Hervé Giordani
Théo Brivoal
Gurvan Madec
A simplified model of the atmospheric boundary layer (ABL) of intermediate complexity between a bulk parameterization and a threedimensional atmospheric model is developed and integrated to the Nucleus for European Modelling of the Ocean (NEMO) general circulation model. An objective in the derivation of such a simplified model, called ABL1d, is to reach an apt representation in oceanonly numerical simulations of some of the key processes associated with air–sea interactions at the characteristic scales of the oceanic mesoscale. In this paper we describe the formulation of the ABL1d model and the strategy to constrain this model with largescale atmospheric data available from reanalysis or realtime forecasts. A particular emphasis is on the appropriate choice and calibration of a turbulent closure scheme for the atmospheric boundary layer. This is a key ingredient to properly represent the air–sea interaction processes of interest. We also provide a detailed description of the NEMOABL1d coupling infrastructure and its computational efficiency. The resulting simplified model is then tested for several boundarylayer regimes relevant to either ocean–atmosphere or seaice–atmosphere coupling. The coupled system is also tested with a realistic 0.25^{∘} resolution global configuration. The numerical results are evaluated using standard metrics from the literature to quantify the wind–seasurfacetemperature (a.k.a. thermal feedback effect), wind–current (a.k.a. current feedback effect), and ABL–seaice couplings. With respect to these metrics, our results show very good agreement with observations and fully coupled ocean–atmosphere models for a computational overhead of about 9 % in terms of elapsed time compared to standard uncoupled simulations. This moderate overhead, largely due to I/O operations, leaves room for further improvement to relax the assumption of horizontal homogeneity behind ABL1d and thus to further improve the realism of the coupling while keeping the flexibility of oceanonly modeling.
Owing to advances in computational power, global oceanic models used for research or operational purposes are now configured with increasingly higher horizontal and vertical resolution, thus resolving the baroclinic deformation radius in the tropics (e.g., Deshayes et al., 2013; Metzger et al., 2014; von Schuckmann et al., 2018). Meanwhile finescale local models are routinely used to simulate submesoscales, which occur on scales on the order of 0.1–20 km horizontally, and their impact on larger scales (e.g., Marchesiello et al., 2011; McWilliams et al., 2019). By increasing the oceanic model resolution, smallscale features are explicitly resolved, but an apt representation of the associated processes also requires the relevant scales to be present in the surface forcings including the proper interaction with the lowlevel atmosphere.
1.1 Historical context
Historically, oceanic general circulation models (OGCMs) were forced by specified wind stress and thermal boundary conditions (from observations or reanalysis) independent from the oceanic state, thus often leading to important drifts in model sea surface properties. To minimize such drifts, a flux correction in the form of a restoration of seasurface temperature and salinity toward climatological values can be added (e.g., Haney, 1971; Barnier et al., 1995). To overcome the shortcomings of the forcing with specified flux, Takano et al. (1973) proposed to use a parameterization of the atmospheric surface layer (ASL) constrained by largescale meteorological data and by the sea state (essentially the seasurface temperature and sometimes the surface currents) to compute the turbulent components of air–sea fluxes. Currently, whatever the target applications, such a technique is widely used in the absence of a concurrently running atmospheric model. Such parameterization of the ASL (known as bulk parameterization, e.g., Beljaars, 1995; Large, 2006), which corresponds to a generalization of the classical neutral wall law to stratified conditions (Monin and Obukhov, 1954), is expected to be valid in the first few tens of meters in the atmosphere. In practice, unless a fully coupled ocean–atmosphere model is used, atmospheric quantities at 10 m, either from existing numerical simulations of the atmosphere or from observations, are prescribed as input to the bulk parameterization. Throughout the paper, this approach will be referred to as “ASL forcing strategy”. A problem with such methodology is that the fast component of the system (the atmosphere) is specified to force the slow component (the ocean), whereas the inertia is in the latter. Indeed, a change in wind stress or heat flux will affect 10 m winds and temperature more strongly than sea surface currents and temperature. In the “ASL forcing strategy”, the key marine atmospheric boundary layer (MABL) processes are not taken into account, and thus feedback loops between the MABL and the upper ocean are not represented.
1.2 Air–sea interactions at oceanic mesoscales
An increasing number of studies based either on observational studies and/or on air–sea coupled simulations have unambiguously shown the existence of air–sea interactions at oceanic mesoscales (e.g., Giordani et al., 1998; Bourras et al., 2004; Chelton and Xie, 2010; Frenger et al., 2013; Schneider and Qiu, 2015; Oerder et al., 2016). Those interactions affect the mass, heat, and momentum exchange between the atmosphere and the ocean. We focus in this work on the dynamical response of the surface wind stress to the seasurface properties (sea surface temperature (SST) and currents) which directly affects lowlevel winds, temperature, and humidity. Several mechanisms responsible for the surface windstress response to SST and oceanic currents can be invoked:
 i.
Downward momentum mixing. SSTinduced changes in the stratification produce significant changes of wind speed and turbulent fluxes throughout the MABL with an increase (decrease) in wind speed over warm water (cold water). As the wind blows over warm water, the MABL becomes more unstable, which leads to an increased vertical mixing, resulting in a downward mixing of momentum from the upper atmosphere to the surface strengthening surface winds on the warm side of an SST front (e.g., Wallace et al., 1989). This mechanism results in a proportional relationship between windstress intensity and SST mesoscale anomalies which has been identified in observations and coupled simulations (e.g., O'Neill et al., 2010; Oerder et al., 2016). Considering spatial derivatives of this proportional relationship leads to a correlation between windstress divergence (curl) and downwind (crosswind) SST–gradient (e.g., Chelton and Xie, 2010; Schneider and Qiu, 2015).
 ii.
Atmospheric pressure adjustment. This mechanism corresponds to an adjustment of the atmospheric pressure gradient to the underlying SST, which manifests itself as a linear relation between horizontal wind divergence and the Laplacian of SST (Lindzen and Nigam, 1987; Minobe et al., 2008; Lambaerts et al., 2013).
 iii.
Oceanic current feedback. The momentum exchange between the ocean and the atmosphere is also largely affected by a dynamical coupling through the dependence of surface wind stress on oceanic surface currents (e.g., Dewar and Flierl, 1987). This coupling results in a drag exerted by the air–sea interface on the ocean which leads to a systematic reduction of the wind power input to the oceanic circulation.
Even if these three mechanisms are mainly active at oceanic eddy scales, they can induce significant effects at larger scales in regions with large SST gradients and/or surface currents (Hogg et al., 2009; Bryan et al., 2010; Renault et al., 2016a). They jointly leave their imprint on the wind divergence, and identifying the relative importance of each mechanism on the momentum balance is difficult because it depends on the dynamical regime and on the spatial and temporal scales of interest (Schneider and Qiu, 2015; Ayet and Redelsperger, 2019).
In the ASL coupling strategy the pressure adjustment mechanism is absent, and only a small fraction of the downward momentum mixing mechanism is accounted for through the modification of the surface drag coefficient depending on the ASL stability (Businger and Shaw, 1984; Chelton and Xie, 2010). As far as the current feedback is concerned, Renault et al. (2016b) showed that the reduction of wind power input to the ocean is systematically overestimated in oceanic simulations based on an ASL forcing strategy compared to air–sea coupled simulations. A simulation that neglects the MABL adjustment to the current feedback cannot represent the partial reenergization of the ocean by the atmosphere and hence overestimates the drag effect by more than 30 % (e.g., Renault et al., 2016b, 2019a). The ASL forcing strategy used in most oceanic models will thus overestimate the current feedback effect and underestimate the downward momentum mixing.
1.3 The proposed approach and focus for this paper
The various aspects discussed so far suggest that a relevant coupling at the characteristic scales of the oceanic mesoscales requires nearly the same horizontal resolution in the ocean and the atmosphere (since the atmosphere must “see” oceanic eddies and fronts) as well as an atmospheric component more complete than a simple ASL parameterization to estimate air–sea fluxes. This assessment raises numerous questions on current practices to force oceanic models across all scales^{1} in the absence of an interactive atmospheric model. The computational cost associated with the systematic use of fully coupled ocean–atmosphere models of similar horizontal resolution is generally unaffordable and comes with practical issues like the proper definition of initial conditions via data assimilation techniques (e.g., Mulholland et al., 2015) and the proper choice of a parameterization set. Moreover, in the fully coupled case at basin or global interannual scales the temporal consistency with the observed variability is generally lost unless a nudging toward observations or reanalysis is done in the atmosphere above the MABL (e.g., Bielli et al., 2009). There is thus clearly room for improvement in the methodology to compute the surface boundary conditions for an ocean model. Alternatives to the ASL forcing strategy have already been suggested by Kleeman and Power (1995), Seager et al. (1995), and Deremble et al. (2013). They proposed a vertically integrated thermodynamically active and dynamically passive MABL model where the wind and the MABL height are specified as in the current practices. Such a model allows a better feedback between SSTs and lowlevel air temperature and humidity because the latter two are prognostic (Abel, 2018). However, by construction, such models do not reproduce the various aforementioned coupling mechanisms affecting the surface wind stress. Their focus is on the improvement of the largescale thermodynamics while ours is on the improvement of the eddyscale momentum exchanges. In the present study we propose an alternative methodology to improve the representation of the downward momentum mixing and of the current feedback effect in oceanonly simulations, leaving aside the pressure gradient adjustment from now on. Our aim is to account for the modulation of atmospheric turbulence by anomalies in seasurface properties in the air–sea flux computation, which is thought to be the main coupling mechanism at the characteristic scales of the oceanic mesoscales. As a step forward beyond the ASL forcing strategy we propose to complement the ASL parameterization with an ABL parameterization while keeping a singlecolumn frame. By construction our approach excludes horizontal advection whose effect can be important in the vicinity of strong SST fronts (e.g., Kilpatrick et al., 2014; Ayet and Redelsperger, 2019). However, we considered that finding a simple and efficient MABL parameterization is the top priority to start investigating the viability of our approach in terms of practical implementation and computational cost. Indeed there exists a large variety of parameterization schemes to represent the effects of subgridscale turbulent mixing in the ABL (see LeMone et al., 2019, and references therein). The schemes based on a diagnostic or prognostic turbulent kinetic energy (TKE) are very popular for operational and research purposes despite wellidentified shortcomings (e.g., Baklanov et al., 2011). For our purposes we do not need the full complexity of the schemes used in practice in atmospheric models because aspects like cloud processes and complex terrains are outside our scope. For this reason, the guideline in this paper is the development and the testing of a simplified version of the TKEbased scheme proposed by Cuxart et al. (2000) for overwater and overseaice conditions. Note that the singlecolumn approximation for our simplified model selected in this study is only a temporary choice to provide evidence on the viability of the whole approach. More advanced formulations allowing a more realistic momentum balance (i.e., including advection) to be recovered will be studied in future work.
1.4 Content
The objective of the present study is to introduce a simplified model of the MABL of intermediate complexity between a bulk parameterization and a full threedimensional atmospheric model and to describe its integration to the Nucleus for European Modelling of the Ocean (NEMO) general circulation model (Madec, 2012). This approach will be referred to as the “ABL coupling strategy”. A constraint in the conception of such a simplified model is to allow an apt representation of the downward momentum mixing mechanism and partial reenergization of the ocean by the atmosphere while keeping the computational efficiency and flexibility inherent to oceanonly modeling. The paper is organized as follows. In Sect. 2, we describe the continuous formulation of the simplified model called ABL1d, including the parameterization scheme used to represent vertical turbulent mixing in the MABL and the strategy to constrain this model with largescale atmospheric conditions. Section 3 provides the description of the dicretization and of the practical implementation of the ABL1d model in the NEMO framework. In Sects. 4 and 5 numerical results obtained for some atmosphereonly simplified test cases available in the literature and for a coupled NEMOABL1d simulation in a global configuration are shown. Finally, our conclusions and perspectives are summarized and discussed in Sect. 6.
In this section we first provide some basic elements on model reduction to motivate our approach and mention possible alternatives (Sect. 2.1). Then we detail the continuous formulation of the ABL1d model and discuss the assumptions made. In particular the governing equations and necessary boundary conditions are given in Sect. 2.2 and the turbulence closure scheme for the MABL in Sect. 2.3. Finally in Sect. 2.4 we discuss the methodology to relax the ABL1d prognostic quantities toward largescale data.
2.1 Motivations and proposed approach
Global oceanic models can be run at higher resolution than global atmospheric models because of their affordable computational cost. From an oceanic perspective, we generally simulate at high resolution (in space and time) ocean fields ${\mathit{\varphi}}_{\mathrm{HR}}^{\mathrm{oce}}(x,y,z,t)$ over a time interval $t\in [\mathrm{0},T]$ over which only largescale atmospheric data ${\mathit{\varphi}}_{\mathrm{LS}}^{\mathrm{atm}}(x,y,z,t)$ are known from the integration of a model ℳ_{atm} using lowerresolution surface oceanic data ${\mathit{\varphi}}_{\mathrm{LS}}^{\mathrm{oce}}(x,y,z=\mathrm{0},t)$ to compute its surface boundary conditions, namely
Instead of directly using ${\mathit{\varphi}}_{\mathrm{LS}}^{\mathrm{atm}}(x,y,z=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{m},t)$ to constrain the oceanic model as in the ASL forcing strategy, our objective is to estimate (without running the full atmospheric model again) the correction to the 10 m largescale atmospheric data associated both with the fine resolution in the oceanic surface fields and with the twoway air–sea coupling. Somehow we aim to find a methodology to get a cheap estimate ${\stackrel{\mathrm{\u0303}}{\mathit{\varphi}}}_{\mathrm{HR}}^{\mathrm{atm}}(x,y,z=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{m},t)$ of the solution that would have been obtained using a coupling of ℳ_{atm} and the oceanic model at high resolution. To do so we could imagine several approaches: (i) estimate $\frac{\partial {\mathcal{M}}_{\mathrm{atm}}}{\partial {\mathit{\varphi}}_{\mathrm{LS}}^{\mathrm{oce}}}$ (i.e., the derivatives of the atmospheric solution with respect to the oceanic parameters) via sensitivity analysis which would require to have the possibility to operate ℳ_{atm}; (ii) build a surrogate model via learning strategies which would require a huge amount of data and computing time; (iii) select the feedback loops of interest and define a simplified model to mimic the underlying physical mechanisms. Following the terminology of Razavi et al. (2012), the first two approaches enter the class of statistical or empirical datadriven models emulating the original model responses while the third one enters the class of lowfidelity physically based surrogates which are built on a simplified version of the original system of equations. In the present study we consider this latter approach, in the spirit of Giordani et al. (2005), who derived a simplified oceanic model by degenerating the primitive equation system and prescribing geostrophic currents into the momentum equation in substitution of the horizontal pressure gradient. In this model, a simple 1D oceanic mixed layer is threedimensionalized via advective terms to couple the vertical columns with each other. The idea here is to translate this idea to the MABL context. In the rest of this section we describe the continuous formulation of our simplified MABL model which will be referred to as ABL1d.
2.2 Formulation of a singlecolumn approach
The formulation of the ABL1d model is derived under the following assumptions: (i) horizontal homogeneity (i.e., ${\partial}_{x}\cdot ={\partial}_{y}\cdot =\mathrm{0}$); (ii) the atmosphere in the computational domain being transparent (i.e., ${\partial}_{z}\mathcal{I}=\mathrm{0}$ with ℐ the radiative flux) meaning that cloud physics is ignored and solar radiation and precipitations at the air–sea interface are specified as usual from observations (e.g., Large and Yeager, 2009); (iii) vertical advection being neglected. Such assumptions prevent the model from prognostically accounting for the SSTinduced adjustment of the atmospheric horizontal pressure gradient and for horizontal advective processes associated with a higher resolution boundary condition at the air–sea interface. The focus here is on the proper representation of the modulation of the MABL turbulent mixing by the air–sea feedback, which is thought to be the main coupling mechanism at the characteristic scales of the oceanic mesoscales impacting ${\mathit{\varphi}}^{\mathrm{atm}}(z=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{m},t)$ and hence air–sea fluxes. This mechanism is expected to explain most of the eddyscale wind–SST and wind–current interactions and is key to properly downscaling largescale atmospheric data produced by a coarseresolution GCM to the oceanic resolution. At a given location in space, the ABL1d model for the Reynoldsaveraged profiles of horizontal velocities u_{h}(z,t), potential temperature θ(z,t), and specific humidity q(z,t), given a suitable initial condition, is
for the height z between a lower boundary z_{sfc} and an upper boundary z_{top}, which will be considered horizontally constant because only the ocean and seaicecovered areas are of interest. In Eq. (1), $\mathit{k}=(\mathrm{0},\mathrm{0},\mathrm{1}{)}^{t}$ is a vertical unit vector, f is the Coriolis parameter, K_{m} and K_{s} are the eddy diffusivity for momentum and scalars respectively, the subscript LS is used to characterize largescale quantities known a priori, λ_{s}(z,t) is the inverse of a relaxation timescale, and R_{LS} denotes a largescale forcing for the momentum equation. R_{LS} can either represent a forcing by geostrophic winds u_{G} (i.e., ${\mathbf{R}}_{\mathrm{LS}}=f\mathit{k}\times {\mathbf{u}}_{\mathrm{G}}$) or equivalently by a horizontal pressure gradient (i.e., ${\mathbf{R}}_{\mathrm{LS}}={\left(\frac{\mathrm{1}}{{\mathit{\rho}}_{\mathrm{a}}}{\mathbf{\nabla}}_{\mathrm{h}}p\right)}_{\mathrm{LS}}$) combined with a standard Newtonian relaxation (i.e., ${\mathbf{R}}_{\mathrm{LS}}={\mathit{\lambda}}_{\mathrm{m}}({\mathbf{u}}_{\mathrm{LS}}{\mathbf{u}}_{\mathrm{h}})$). Because of the simplifications made to derive the ABL1d model the R_{LS} term and a nonzero λ_{s} are necessary to prevent the prognostic variables from drifting very far away from the largescale values used to “guide” the model. By itself a relaxation term does not directly represent any real physical process, but the rationale is that it accounts for the influence of largescale threedimensional circulation processes not explicitly represented in a simple 1D model. Note that this methodology is currently used to evaluate GCM parameterizations in 1D column models. Once the turbulent mixing and the Coriolis term have been computed to provide a provisional prediction ${\mathit{\varphi}}^{n+\mathrm{1},\star}$ at time n+1 for any ABL1d prognostic variable ϕ, the relaxation term provides a weighting between this prediction and the largescale quantities:
with Δt the increment of the temporal discretization. Above the boundary layer, the ABL1d formulation is unable to properly represent the physics; therefore the λ parameter should be large, while in the first tens of meters near the surface we expect the ABL1d model to accurately represent the interaction with the fineresolution oceanic state, and thus the relaxation toward ϕ_{LS} should be small. The exact form of the λ_{s} and λ_{m} coefficients is discussed later in Sect. 2.4. Note that because of the relaxation term, threedimensional atmospheric data for u_{LS}, θ_{LS}, q_{LS}, and possibly ${\left(\frac{\mathrm{1}}{{\mathit{\rho}}_{\mathrm{a}}}{\mathbf{\nabla}}_{\mathrm{h}}p\right)}_{\mathrm{LS}}$ sampled between z_{sfc} and z_{top} must be provided to the oceanic model instead of the twodimensional data (usually at 10 m) necessary for an ASL forcing strategy. Since the ABL1d model does not include any representation of radiative processes and microphysics, the radiative fluxes and precipitation at the air–sea interface are similar to the one provided for a standard uncoupled oceanic simulation. The model requires boundary conditions for the vertical mixing terms which are computed via a standard bulk formulation:
For the sake of consistency, it is preferable to use a bulk formulation as close as possible to the one used to compute the threedimensional largescale atmospheric data ${\mathit{\varphi}}_{\mathrm{LS}}^{\mathrm{atm}}$. Because in the present study the plan is to use a largescale forcing from ECMWF reanalysis products, we use the IFS (Integrated Forecasting System: https://www.ecmwf.int/en/forecasts/documentationandsupport/changesecmwfmodel/ifsdocumentation, last access: 20 January 2021) bulk formulation such as implemented in the AeroBulk (https://brodeau.github.io/aerobulk/, last access: 20 January 2021) package (Brodeau et al., 2017) to compute C_{D}, C_{θ}, and C_{q} in realistic simulations (see Sect. 5). Note that for an ASL forcing strategy u_{h}(z_{sfc}) and ϕ(z_{sfc}) in Eq. (3) would be equal to u_{LS}(z=10 m) and ϕ_{LS}(z=10 m) respectively, while in the ABL coupling strategy those variables are provided prognostically by an ABL1d model. As far as the boundary conditions at z=z_{top} are concerned, Dirichlet boundary conditions u_{h}(z_{top})=u_{LS}(z_{top}) and ϕ(z_{top})=ϕ_{LS}(z_{top}) are prescribed. Model (1) is a first step before evolving toward a more advanced surrogate model including horizontal advection and finescale pressure gradients in the future. A particular focus of the present study is on the appropriate choice of a closure scheme to diagnose the eddy diffusivities K_{m} and K_{s}. This is a key step to properly represent the downward mixing process.
2.3 Turbulence closure scheme
This subsection describes the turbulence scheme used to compute the eddy diffusivity for momentum and scalars. Those eddy diffusivities are responsible for a vertical mixing of atmospheric variables due to turbulent processes. The turbulence scheme we have implemented in our ABL1d model is very similar to the socalled CBR1d scheme of Cuxart et al. (2000), which is used operationally at Meteo France (Bazile et al., 2012). We chose to recode the parameterization from scratch for several reasons: computational efficiency, consistency with the NEMO coding rules, use of a geopotential vertical coordinate, and flexibility to add elements specific to the marine atmospheric boundary layer. CBR1d is a oneequation turbulence closure model based on a prognostic TKE and a diagnostic computation of appropriate length scales. The prognostic equation for the TKE $e=\frac{\mathrm{1}}{\mathrm{2}}\left(\u2329{u}^{\prime}{u}^{\prime}\u232a+\u2329{v}^{\prime}{v}^{\prime}\u232a+\u2329{w}^{\prime}{w}^{\prime}\u232a\right)$ (with 〈⋅〉 the Reynolds averaging operator) is
where horizontal terms and vertical advection are neglected, as usually done in mesoscale atmospheric models. Here θ_{v} is the virtual potential temperature, ρ_{a} is atmospheric density, and ε is a dissipation term. In order to express the evolution of e in terms of Reynoldsaveraged atmospheric variables we consider the standard closure assumptions for the first order turbulent fluxes (Cuxart et al., 2000) to obtain the classical TKE prognostic equation
where l_{ε} is a dissipative length scale, c_{ϵ} a constant, and N^{2} is the moist Brunt–Väisälä frequency computed as ${N}^{\mathrm{2}}=(g/{\mathit{\theta}}_{\mathrm{v}}^{\mathrm{ref}})\left({\partial}_{z}\u2329\mathit{\theta}\u232a+\mathrm{0.608}\phantom{\rule{0.25em}{0ex}}{\partial}_{z}\left(\u2329\mathit{\theta}\u232a\u2329q\u232a\right)\right)$ with ${\mathit{\theta}}_{\mathrm{v}}^{\mathrm{ref}}=\mathrm{288}\phantom{\rule{0.25em}{0ex}}\mathrm{K}$. The eddy diffusivities for momentum K_{m}, TKE K_{e}, and scalars K_{s} all depend on e and on a mixing length scale l_{m}:
with $({C}_{\mathrm{m}},{C}_{\mathrm{s}},{C}_{\mathrm{e}})$ a triplet of constants and ϕ_{z} a stability function proportional to the inverse of a turbulent Prandtl number, given by
The ϕ_{z} function is bounded not to exceed ${\mathit{\varphi}}_{z}^{max}=\mathrm{2.2}$ as done in the Arpege model of Meteo France (e.g., Bazile et al., 2012). Assuming that the minimum of ϕ_{z} is attained in the linearly stratified limit (i.e., for ${l}_{\mathrm{m}}={l}_{\mathit{\epsilon}}=\sqrt{\mathrm{2}e/{N}^{\mathrm{2}}}$), values of the maximum Prandtl number ${\mathrm{Pr}}_{\mathrm{t}}={C}_{\mathrm{m}}/\left({C}_{\mathrm{s}}{\mathit{\varphi}}_{z}\right)$ are given in Table 1. Constant values for C_{m}, C_{s}, C_{e}, c_{ε}, and C_{1} can be determined from different methods, leading to nearly similar values. The traditional way is to use the inertial–convective subrange theory of locally isotropic turbulence (Lilly, 1967; Deardorff, 1974). Another way relies on a theoretical turbulence model partly based on renormalization group methods (see Cheng et al., 2002). For the present study, the sets proposed by Cuxart et al. (2000) and Cheng et al. (2002) will be considered (Table 1). A major difference between the two sets concerns the value of C_{m}. This difference is explained by a reevaluation of the energy redistribution among velocity components by pressure fluctuations, whose magnitude is assumed to be proportional to the degree of energy anisotropy as initially introduced by Rotta (1951). Note that the constant set of Cheng et al. (2002) is now used by default in both research and operational Meteo France models.
Cuxart et al. (2000)Cheng et al. (2002)The Dirichlet boundary condition for TKE applied at the top z=z_{top} is $e(z={z}_{\mathrm{top}})={e}_{min}={\mathrm{10}}^{\mathrm{6}}$ m^{2} s^{−2} and at the bottom z=z_{sfc} we have
with u_{⋆} and w_{⋆} the friction and convective velocities given by the bulk formulation. The value for e_{min} has been chosen empirically as well as background values ${K}_{\mathrm{m}}^{min}={\mathrm{10}}^{\mathrm{4}}$ m^{2} s^{−1} and ${K}_{\mathrm{s}}^{min}={\mathrm{10}}^{\mathrm{5}}$ m^{2} s^{−1} for eddy diffusivities.
The minimum value for l_{m} is simply set as ${l}_{\mathrm{m}}^{min}=\frac{{K}_{\mathrm{m}}^{min}}{{C}_{\mathrm{m}}\sqrt{{e}_{min}}}$. There are multiple options to compute the mixing lengths l_{m} and l_{ε} (this point will be discussed later in Sect. 3.2.2), but all options have identical boundary conditions ${l}_{\mathrm{m}}(z={z}_{\mathrm{top}})={l}_{\mathrm{m}}^{min}$ and
The value of L_{sfc} results from the similarity theory in the neutrally stratified surface layer (Sect. 4.1 in Redelsperger et al., 2001, and Appendix A). In Eq. (8), κ is the von Karman constant and z_{0} a roughness length computed within the bulk algorithm. The way e_{sfc} and L_{sfc} are obtained is detailed in Appendix A.
Our current implementation of boundary layer subgrid processes is an eddydiffusivity approach which does not include any explicit representation of boundarylayer convective structures. This could be done via a massflux representation (e.g., Hourdin et al., 2002; Soares et al., 2004) or the introduction of a countergradient term (e.g., Troen and Mahrt, 1986). This point is left for future developments of the ABL1d model.
2.4 Processing of largescale forcing and Newtonian relaxation
As mentioned earlier, the ABL1d model (1) requires threedimensional $(x,y,z)$ largescale atmospheric variables ${\mathit{\varphi}}_{\mathrm{LS}}^{\mathrm{atm}}$, while existing uncoupled oceanic forcing strategies require only twodimensional (x,y) atmospheric variables. This is a difficulty for efficiency reasons since it substantially increases the number of I/Os but also for practical reasons because it requires the development of a dedicated tool to extract largescale atmospheric data and interpolate them on prescribed geopotential heights from their native vertical grid, which can be either pressure based or arbitrary Lagrangian Eulerian. Such tools have been developed specifically to work with ERAInterim, ERA5, and operational IFS datasets and are described in Appendix B.
Beyond the particular values of ${\mathit{\varphi}}_{\mathrm{LS}}^{\mathrm{atm}}$, the form of the relaxation timescale has a great impact on model solutions. The vertical profile for the λ_{m} and λ_{s} coefficients in Eq. (1) is chosen to nudge strongly above the MABL and moderately in the MABL with a smooth transition between its minimum and maximum value to avoid large vertical gradients in λ_{m} and λ_{s}, which would result in artificially large vertical gradients in atmospheric variables. In practice the λ_{m}(z) and λ_{s}(z) functions depend on the following parameters:

$({\mathit{\lambda}}_{\mathrm{m}}^{\mathrm{max}},{\mathit{\lambda}}_{\mathrm{m}}^{\mathrm{min}})$ and $({\mathit{\lambda}}_{\mathrm{s}}^{\mathrm{max}},{\mathit{\lambda}}_{\mathrm{s}}^{\mathrm{min}})$, which define the maximum and minimum of the nudging coefficient for momentum and scalars respectively. Following Eq. (2), a guideline to set reasonable values for those parameter values would be to make sure that $\mathrm{\Delta}t{\mathit{\lambda}}_{\mathrm{s}}^{max}\approx \mathrm{1}$ (i.e., the largescale value is imposed above the boundary layer) and choose ${\mathit{\lambda}}_{\mathrm{s}}^{min}$ based on the typical adjustment timescale of the ABL to surface perturbations. Broadly speaking the ABL can be defined as the region that responds to surface forcings with a timescale of about an hour (e.g., LeMone et al., 2019). In the realistic numerical experiments shown in Sect. 5, we used ${\mathit{\lambda}}_{\mathrm{s}}^{min}=\frac{\mathrm{1}}{\mathrm{90}\left[\mathrm{min}\right]}$, which, for an oceanic dynamical time step Δt=1080 s, would lead to $\mathrm{\Delta}t{\mathit{\lambda}}_{\mathrm{s}}^{min}=\mathrm{0.2}$. (i.e., the boundary layer values are the result of a weighting with a weight 0.8 for the ABL1d prediction and 0.2 for the largescale value).

$({\mathit{\beta}}_{min},{\mathit{\beta}}_{max})$, which defines the extent of the transition zone separating the maximum and minimum of the nudging coefficient
We considered the following general form for λ_{s}(z) and λ_{m}(z), with h_{bl} the boundary layer height whose value is diagnosed using an integral Richardson number criteria (Sect. 3.2 and 3.3 in Lemarié et al., 2012) with a critical value equal to C_{1}:
where four α_{m} coefficients are necessary to guarantee the continuity of λ_{s}(z) and its derivative ∂_{z}λ_{s} at $z={\mathit{\beta}}_{min}{h}_{\mathrm{bl}}$ and $z={\mathit{\beta}}_{max}{h}_{\mathrm{bl}}$. We easily find
The value of h_{bl} is bounded beforehand to guarantee that at least 3 grid points are such that $z\le {\mathit{\beta}}_{min}{h}_{\mathrm{bl}}$ and $z\ge {\mathit{\beta}}_{max}{h}_{\mathrm{bl}}$. A typical profile of the λ_{s}(z) is shown in Fig. 1a.
When the model is forced by the largescale pressure gradient (or the geostrophic winds), the parameter λ_{m}(z) should be theoretically zero at high and middle latitudes. However, for the equatorial region, a Newtonian relaxation toward the largescale winds should be maintained. To do so, the coefficient λ_{m}(z) is multiplied by a coefficient r_{eq}, which is a function of the Coriolis parameter f. The r_{eq} coefficient is equal to zero for large values of $\leftf\right$ and increases to 1 when approaching the Equator. The following form satisfies those constraints (see also Fig. 1b):
We have introduced so far the continuous formulation of the ABL1d model. In this section we describe the discretization methods used and how this model is included in the NEMO modeling framework. In particular, the discretization of the Coriolis term and of the TKE Eq. (6) and associated mixing lengths are described in Sect. 3.1 and 3.2 respectively. Details about the practical implementation in NEMO are given in Sect. 3.3 for the coupling aspects and 3.4 for the computational aspects. The ABL1d model (1) is discretized in time with an Euler backward scheme for the vertical diffusion terms, semiimplicitly for the Coriolis term, and explicitly for the relaxation term, which means that the model is stable as long as λ_{s}Δt≤1. The variables are defined on a nonstaggered grid in the horizontal (a.k.a Arakawa Agrid). Because we consider a computational domain exclusively over water or sea ice, topography is not considered and vertical levels are flat and fixed in time which, among other things, allows the largescale data ϕ_{LS} to be interpolated on the vertical grid offline. The position of the various quantities introduced so far on the computational grid is given in Fig. 2.
3.1 Coriolis term treatment
Since in our implementation the horizontal velocity components are collocated, the discretization of the Coriolis term is straightforward and is energetically neutral. In the event the ABL1d is integrated with a time step much larger than the oceanic time step, specific care must be given to the stability of the Coriolis term time stepping. A semiimplicit scheme with weighting parameter γ reads ${\mathbf{u}}_{\mathrm{h}}^{n+\mathrm{1},\star}=\left(f\mathrm{\Delta}t\right)\mathit{k}\times \left((\mathrm{1}\mathit{\gamma}){\mathbf{u}}_{\mathrm{h}}^{n}+\mathit{\gamma}{\mathbf{u}}_{\mathrm{h}}^{n+\mathrm{1},\star}\right)$, where the exponent ⋆ is used here to emphasize that ${\mathbf{u}}_{\mathrm{h}}^{n+\mathrm{1},\star}$ is a temporary value at time n+1 before vertical diffusion and Newtonian relaxation are applied. For a given grid cell with index (i,j), the semiimplicit scheme can be written in a more compact way as
The associated amplification factor modulus is $\left{\mathcal{A}}_{\mathrm{cor}}\right=\sqrt{\frac{\mathrm{1}+(\mathrm{1}\mathit{\gamma}{)}^{\mathrm{2}}(f\mathrm{\Delta}t{)}^{\mathrm{2}}}{\mathrm{1}+{\mathit{\gamma}}^{\mathrm{2}}(f\mathrm{\Delta}t{)}^{\mathrm{2}}}}$ meaning that unconditional stability is obtained as long as $\mathit{\gamma}\ge \mathrm{1}/\mathrm{2}$. For the numerical results obtained below in Sects. 4 and 5 we used γ=0.55, which is deliberately slightly dissipative.
3.2 Discretization of TKE equation
In Sect. 2.3 we have presented the continuous formulation of the TKEbased turbulence closure of the ABL1d model. In the following we describe how the positivity of TKE can be preserved and how the mixing lengths l_{m} and l_{ε} are computed. We provide a substantial discussion on the latter aspect because numerical results are very sensitive to the choices made.
3.2.1 TKE positivity preservation
The TKE equation is discretized using a backward Euler scheme in time with a linearization of the dissipation term $\frac{{c}_{\mathit{\epsilon}}}{{l}_{\mathit{\epsilon}}}{e}^{\mathrm{3}/\mathrm{2}}$, which is discretized as $\frac{{c}_{\mathit{\epsilon}}}{{l}_{\mathit{\epsilon}}}\sqrt{{e}^{n}}{e}^{n+\mathrm{1}}$. However, such discretization is not unconditionally positivitypreserving for TKE, which could give rise to unphysical solutions (e.g., Burchard, 2002b). Ignoring the diffusion term, the TKE prognostic Eq. (6) can be written as an ordinary differential equation (ODE) of the form
where the last term can be seen as a damping term. For ODEs like Eq. (11) it can be shown that for an initial condition e(0)≥0 and $S({\mathbf{u}}_{\mathrm{h}},{N}^{\mathrm{2}})\ge \mathrm{0}$, the solution e(t) keeps the same sign as e(0) whatever the sign of the damping coefficient D(e,t). Assuming that S(u_{h},N^{2}) and D(e,t) are positive, a backward Euler discretization of the damping term in Eq. (11) would lead to ${e}^{n+\mathrm{1}}=\frac{{e}^{n}+\mathrm{\Delta}tS({\mathbf{u}}_{\mathrm{h}},{N}^{\mathrm{2}})}{\mathrm{1}+\mathrm{\Delta}tD(e,t)}$, which preserves positivity since for e^{n}≥0 we obtain ${e}^{n+\mathrm{1}}\ge \mathrm{0}$. However, there is no guarantee that the forcing term S(u_{h},N^{2}) is positive, in particular when the shear is weak and the stratification is large. When S(u_{h},N^{2}) is negative a specific treatment (known as the “Patankar trick”; see Deleersnijder et al., 1997; Burchard, 2002b) is required. In the event of a negative S(u_{h},N^{2}), the idea is to move the buoyancy term from S to D after dividing it by e^{n}, such that $S({\mathbf{u}}_{\mathrm{h}},{N}^{\mathrm{2}})={K}_{\mathrm{m}}\Vert {\partial}_{z}{\mathbf{u}}_{\mathrm{h}}{\Vert}^{\mathrm{2}}$ is now strictly positive and $D(e,t)=\frac{{c}_{\mathit{\epsilon}}}{{l}_{\mathit{\epsilon}}}\sqrt{{e}^{n}}+{K}_{\mathrm{s}}\frac{{N}^{\mathrm{2}}}{{e}^{n}}$. Such a procedure is a sufficient condition to preserve the positivity of the TKE without ad hoc clipping of negative values. Moreover our discretization of the shear and buoyancy terms in the TKE equation is done in an energetically consistent way following Burchard (2002a).
3.2.2 Mixing length computation
Another challenging task when implementing a TKE scheme is the discretization of the mixing lengths. As mentioned earlier, four different discretizations of l_{m} (l_{ε}) have been coded. All discretizations consider the boundary conditions given in Eq. (8). The values of l_{m} and l_{ε} are traditionally computed from two intermediate length scales l_{up} and l_{dwn}, which respectively correspond to the maximum upward and downward displacement of a parcel of air with a given initial kinetic energy. Once l_{up} and l_{dwn} have been estimated by one of the methods described below, the dissipative and mixing length scales l_{m} and l_{ε} are computed as
where $a\approx \frac{\mathrm{3}}{\mathrm{2}}$ for CBR00 and $a\approx \frac{\mathrm{6}}{\mathrm{7}}$ for CCH02 (see Appendix A). The impact of the weighting between l_{up} and l_{dwn} to compute l_{m} can be significant for idealized experiments like the ones presented in Sect. 4.2 but for more realistic cases results are weakly sensitive and equivalent to the ones obtained with the simpler weighting ${l}_{\mathrm{m}}=\sqrt{{l}_{\mathrm{up}}{l}_{\mathrm{dwn}}}$.
In the following we provide the continuous form of the various ways to compute l_{up} and l_{dwn} implemented in the ABL1d model. The discretization aspects are detailed in Appendix C.

Bougeault and Lacarrère (1989) length scale. A classical approach in atmospheric models is the use of the Bougeault and Lacarrère (1989) mixing length (see also Bougeault and André, 1986) which defines l_{up} and l_{dwn} as
$$\begin{array}{}\text{(13)}& \begin{array}{rl}& \underset{z}{\overset{z+{l}_{\mathrm{up}}}{\int}}{N}^{\mathrm{2}}\left(s\right)(sz)\mathrm{d}s=e\left(z\right),\\ & \phantom{\rule{1em}{0ex}}\underset{z{l}_{\mathrm{dwn}}}{\overset{z}{\int}}{N}^{\mathrm{2}}\left(s\right)(zs)\mathrm{d}s=e\left(z\right).\end{array}\end{array}$$By construction such mixing lengths are bounded by the distance to the bottom and the top of the computational domain. It is worth noting that for constant values of N^{2}, Eq. (13) gives $\frac{{l}_{\mathrm{up}}^{\mathrm{2}}{N}^{\mathrm{2}}}{\mathrm{2}}=e\left(z\right)$ and $\frac{{l}_{\mathrm{dwn}}^{\mathrm{2}}{N}^{\mathrm{2}}}{\mathrm{2}}=e\left(z\right)$ respectively, which is equivalent to the Deardorff (1980) length scale. In the remainder we will note l_{BL89}, the mixing length obtained from Eq. (13).

Adaptation of NEMO's length scale. The standard NEMO algorithm (Sect. 10.1.3 in Madec, 2012) is simple and efficient compared to Eq. (13). This algorithm is based on the Deardorff (1980) length scale ${l}_{\mathrm{D}\mathrm{80}}=\sqrt{\mathrm{2}e\left(z\right)/{N}^{\mathrm{2}}}$. l_{up} and l_{dwn} are first initialized to ${l}_{\mathrm{up}}={l}_{\mathrm{dwn}}={l}_{\mathrm{D}\mathrm{80}}$. The resulting length scales are then limited not only by the distance to the surface and to the top but also by the distance to a strongly stratified portion of the air column. This limitation amounts to control of the vertical gradients of l_{up}(z) and l_{dwn}(z) such that they are not larger than the variations of altitude. The resulting mixing length will be simply referred to as l_{D80}. Note that the Taylor expansion of the integral in Eq. (13) is
$$}\underset{z}{\overset{z+{l}_{\mathrm{up}}}{\int}}{N}^{\mathrm{2}}\left(s\right)(sz)\mathrm{d}s\approx {\displaystyle \frac{{N}^{\mathrm{2}}\left(z\right){l}_{\mathrm{up}}^{\mathrm{2}}}{\mathrm{2}}}+{\displaystyle \frac{\frac{d{N}^{\mathrm{2}}}{\mathrm{d}z}{l}_{\mathrm{up}}^{\mathrm{3}}}{\mathrm{3}}}+\mathcal{O}\left({l}_{\mathrm{up}}^{\mathrm{4}}\right),$$which shows that the l_{D80} mixing length is an approximation of l_{BL89}, which is obtained by retaining only the leading order term in the Taylor expansion.

Rodier et al. (2017) length scale. Recently, Rodier et al. (2017) proposed a modification of the Bougeault and Lacarrère (1989) mixing length. This modification turns out to improve results for stably stratified boundary layers typical of areas covered by ice. They propose to add a shearrelated term to Eq. (13) such that the definition of l_{up} and l_{dwn} becomes
$$\begin{array}{}\text{(14)}& \begin{array}{rl}& \underset{z}{\overset{z+{l}_{\mathrm{up}}}{\int}}\left[{N}^{\mathrm{2}}\left(s\right)(sz)+{c}_{\mathrm{0}}\sqrt{e\left(s\right)}\Vert {\partial}_{\mathrm{s}}{\mathbf{u}}_{\mathrm{h}}\Vert \right]\mathrm{d}s=e\left(z\right),\\ & \phantom{\rule{1em}{0ex}}\underset{z{l}_{\mathrm{dwn}}}{\overset{z}{\int}}\left[{N}^{\mathrm{2}}\left(s\right)(zs)+{c}_{\mathrm{0}}\sqrt{e\left(s\right)}\Vert {\partial}_{\mathrm{s}}{\mathbf{u}}_{\mathrm{h}}\Vert \right]\mathrm{d}s=e\left(z\right),\end{array}\end{array}$$where c_{0} is a parameter whose value should be smaller than $\sqrt{{C}_{\mathrm{m}}/{c}_{\mathit{\epsilon}}}$. The value of c_{0} will be chosen based on numerical experiments presented in Sect. 4. In the following this mixing length will be referred to as l_{R17}.

A local buoyancy and shearbased length scale. For the sake of computational efficiency, we have derived a local version of the Rodier et al. (2017) length scale which is original to the present paper. Under the assumption that l_{up} (l_{dwn}) is small compared to the spatial variations of N^{2}, e, and $\Vert {\partial}_{z}{\mathbf{u}}_{\mathrm{h}}\Vert $, we end up with the following secondorder equation for l_{up}:
$$\frac{{N}^{\mathrm{2}}\left(z\right)}{\mathrm{2}}}{l}_{\mathrm{up}}^{\mathrm{2}}+{c}_{\mathrm{0}}\sqrt{e\left(z\right)}\Vert {\partial}_{z}{\mathbf{u}}_{\mathrm{h}}\Vert {l}_{\mathrm{up}}=e\left(z\right),$$whose unique positive solution is
$$}{l}_{\mathrm{D}\mathrm{80}}^{\star}\left(z\right)={\displaystyle \frac{\mathrm{2}\sqrt{e\left(z\right)}}{{c}_{\mathrm{0}}\Vert {\partial}_{z}{\mathbf{u}}_{\mathrm{h}}\Vert +\sqrt{{c}_{\mathrm{0}}^{\mathrm{2}}\Vert {\partial}_{z}{\mathbf{u}}_{\mathrm{h}}{\Vert}^{\mathrm{2}}+\mathrm{2}{N}^{\mathrm{2}}\left(z\right)}}}.$$We easily find that ${l}_{\mathrm{D}\mathrm{80}}^{\star}={l}_{\mathrm{D}\mathrm{80}}$ for $\Vert {\partial}_{z}{\mathbf{u}}_{\mathrm{h}}\Vert =\mathrm{0}$, and ${l}_{\mathrm{D}\mathrm{80}}^{\star}=\frac{\sqrt{e\left(z\right)}}{{c}_{\mathrm{0}}\Vert {\partial}_{z}{\mathbf{u}}_{\mathrm{h}}\Vert}$ for N^{2}=0, which is consistent with the shearbased length scale of Wilson and Venayagamoorthy (2015). Once ${l}_{\mathrm{D}\mathrm{80}}^{\star}$ has been computed we apply the same algorithmic approach as in the l_{D80} case.
The performance of those four length scales for various physical flows is discussed in Sect. 4.
3.3 Coupling with ocean and sea ice
For the practical implementation of the ABL coupling strategy within a global oceanic model, a proper coupling method is required for stability and consistency purposes (e.g., Beljaars et al., 2017; Renault et al., 2019a), and the ABL1d must have the ability to handle grid cells partially covered by sea ice. For the coupling strategy, a socalled implicit flux coupling which is unconditionally stable (Appendix B in Beljaars et al., 2017) and asymptotically consistent for Δt→0 (Renault et al., 2019a) is used. Because vertical diffusion in ABL1d is handled implicitly in time, the boundary conditions (Eqs. 3 and 4) should be provided at time n+1. The implicit flux coupling amounts to discretize the boundary conditions Eqs. (3) and (4) as
where ${\stackrel{\mathrm{\u0303}}{\mathbf{u}}}_{\mathrm{oce}}$ and ${\stackrel{\mathrm{\u0303}}{\mathit{\varphi}}}_{\mathrm{oce}}$ are either the instantaneous values at time n if NEMO and ABL1d have the same time step or an average over the successive oceanic substeps otherwise.
Particular care has also been given to the compatibility between the ABL1d model and SI^{3} (Sea Ice model Integrated Initiative) the seaice component of NEMO. SI^{3} is a multicategory model whose state variables relevant for our study are the ice surface temperature ${T}_{l}^{\mathrm{ice}}$ with associated fractional area a_{l} (for the lth category), and the ice velocity u^{ice} (same for all categories). Note that the values of the exchange coefficients over sea ice ${C}_{D}^{\mathrm{ice}}$, ${C}_{\mathit{\theta}}^{\mathrm{ice}}$, and ${C}_{q}^{\mathrm{ice}}$ are different from their oceanic counterparts but are the same over all seaice categories. At this point there are several strategies for the ABL1d∕SI^{3} coupling:

Run the ABL1d model over the whole ABL for each category l and then average atmospheric variables weighted by a_{l}.

Run a single ABL1d model with a categoryaveraged surface flux. In the current version of NEMO ${C}_{\mathit{\theta}}^{\mathrm{ice}}$ is a function of the averaged temperature T^{ice} which means that it is equivalent to compute a flux over each category before averaging them and to compute a single flux using the averaged surface temperature, indeed
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}\sum _{l}{a}_{l}\left[{C}_{\mathit{\theta}}^{\mathrm{ice}}\Vert {\mathbf{u}}_{\mathrm{h}}\left({z}_{\mathrm{sfc}}\right){\mathbf{u}}^{\mathrm{ice}}\Vert \left(\mathit{\theta}\right({z}_{\mathrm{sfc}}){T}_{l}^{\mathrm{ice}})\right]=\\ {\displaystyle}& {\displaystyle}{C}_{\mathit{\theta}}^{\mathrm{ice}}\Vert {\mathbf{u}}_{\mathrm{h}}\left({z}_{\mathrm{sfc}}\right){\mathbf{u}}^{\mathrm{ice}}\Vert \left(\mathit{\theta}\left({z}_{\mathrm{sfc}}\right)\sum _{l}{a}_{l}{T}_{l}^{\mathrm{ice}}\right).\end{array}$$
The second option has been preferred because it is much easier to implement and more computationally efficient. It amounts to consider an ice surface temperature averaged over all categories ${T}^{\mathrm{ice}}={\sum}_{l=\mathrm{1}}^{{n}_{\mathrm{cat}}}{a}_{l}{T}_{l}^{\mathrm{ice}}$ for the computation of ice–atmosphere turbulent fluxes (T^{ice} also enters in the computation of q_{ice}). Noting F_{oce} the fraction of open water (lead), the boundary condition (15) and (16) are modified in
Because the dynamics of sea ice is computed before the thermodynamics (see Fig. 1 in Rousset et al., 2015), the ABL1d–SI^{3} coupling follows these different steps:

compute surface fluxes over ice and ocean and integrate the ABL1d model for given values ${F}_{\mathrm{oce}}^{n}$ and ${a}_{l}^{n}$,

compute the dynamics of sea ice,

update ${F}_{\mathrm{oce}}^{n}$ and ${a}_{l}^{n}$ in ${F}_{\mathrm{oce}}^{\star}$ and ${a}_{l}^{\star}$ because of step 2,

distribute the fluxes over each ice category considering the updated values ${a}_{l}^{\star}$ (Sect. 3.6 in Rousset et al., 2015)

compute the thermodynamics of sea ice.
3.4 Computational aspects
As described in Maisonnave and Masson (2015), the NEMO source code is organized to separate the ocean routines on one side and the routines responsible for the surface boundary conditions computation (including sea ice and the coupling interfaces) on the other side. This makes a clear separation between the standard ocean model (OCE component) and the socalled surface module (SAS component). As schematically described in Fig. 3, the ABL1d model has been implemented within the SAS component, which allows the following useful features:

The ABL1d model can be run in standalone mode (coupled or not with sea ice) with prescribed oceanic surface fields.

The ABL1d model can be run in detached mode; i.e., the OCE and SAS components run on potentially separate processors and computational grids communicating via the OASIS3MCT coupling library (Craig et al., 2017).
An other capability implemented within the NEMO modeling framework is the possibility to interpolate forcing fields on the fly. This is particularly useful for the ABL coupling strategy since threedimensional atmospheric data must be interpolated on the ABL1d computational grid. As the current implementation of the onthefly interpolation only works in the horizontal, the vertical interpolation of largescale atmospheric data on the ABL1d vertical grid is done offline. Nevertheless it means that the size of input data compared to an ASL forcing strategy is N times larger with N the number of vertical levels in the ABL. A possibility to improve the efficiency for the reading of input data would be to take advantage of the parallel I/O capabilities provided by the XIOS library (XMLIOServer; Meurdesoif et al., 2016) which is currently used in NEMO only for writing output data. This technical development is left for future work. This is a key aspect because, as discussed later in Appendix D, the main source of computational overhead associated with the ABL coupling strategy is due to the time spent waiting for input files to be read.
4.1 Sensitivity experiments and objectives
To check the relevance of our ABL1d model for idealized atmospheric situations typical of the atmospheric boundary layer over water or sea ice, we performed a set of singlecolumn experiments. Each of those experiments are evaluated with benchmark large eddy simulations (LESs). Moreover, we use standardized test cases from the literature to allow our results to be crosscompared with other wellestablished ABL schemes. In the following we consider a neutrally stratified (Sect. 4.2) and a stably stratified (Sect. 4.3) case as well as a case with a transition from stable to unstable stratification representative of an atmospheric flow over an SST front (Sect. 4.4). All ABL1d simulations presented here have been performed directly within the SAS component of the NEMO modeling framework and can be reproduced using the code available at https://zenodo.org/record/3904518 (last access: 20 January 2021) (Lemarié and Samson, 2020), which also includes the scripts to generate the figures. An objective of the present section is to illustrate the type of sensitivity we can expect from the ABL1d model and discriminate between the various options available in the code. The experiments showed in Sect. 4.2 and 4.3 are meant to investigate the impact of (i) the set of constant coefficients (CBR00 vs. CCH02), (ii) the various formulations of l_{m} and l_{ε} among the algorithms described in Sect. 3.2.2, and (iii) the parameter value c_{0} in the l_{R17} and ${l}_{\mathrm{D}\mathrm{80}}^{\star}$ mixing length computation. Those experiments will allow several options to be discarded. The ability of the remaining options to represent the downward mixing mechanism, discussed in Sect. 1.2, is then evaluated in Sect. 4.4. The robustness of the results to the bulk formulation and to the nudging coefficient is also checked. For each experiment we explicitly provide the initial and boundary conditions as well as all the necessary parameter values (see Table 2) so that the experiments can be reproduced easily by other modeling groups.
4.2 Neutral turbulent Ekman layer
We first propose to investigate the simulation of a neutrally stratified atmosphere analogous to a classical turbulent Ekman layer. The selected case is based on the setup described in Andren et al. (1994). The initial conditions for this experiment are not defined analytically; they are given by Table A1 in Andren et al. (1994)^{2}. This test case is mainly used to check the adequacy of our surface boundary conditions with similarity theory and the proper calibration of the parameter c_{0} in the ${l}_{\mathrm{D}\mathrm{80}}^{\star}$ and l_{R17} formulations of the mixing lengths. In theory, the l_{D80} and l_{BL89} mixing lengths do not support the asymptotic limit N^{2}=0 but for the integrity of numerical results a minimum threshold ${N}_{\mathit{\epsilon}}^{\mathrm{2}}$ on the stratification is imposed in the code. In this case the procedure to compute those mixing lengths as described in Appendix C will provide identical results, namely ${l}_{\mathrm{up}}={z}_{\mathrm{top}}z$ and ${l}_{\mathrm{dwn}}=z{z}_{\mathrm{sfc}}$ (i.e., the distance from the top and from the bottom of the computational domain). We test here the ${l}_{\mathrm{D}\mathrm{80}}^{\star}$ and l_{R17} introduced in Sect. 3.2.2. The reference solution is taken from Cuxart et al. (2000) (panels a and b in their Fig. 16). Results are obtained using the ABL1d model with either the CBR00 (Fig. 4a–d) or the CCH02 (Fig. 4e–h) set of parameters. All experiments have been done with c_{0}=0.15 and c_{0}=0.2. All simulations are able to reproduce the overall behavior of the LES case. The main outcomes are as follows:

The best agreement is obtained when using the CCH02 constants along with ${l}_{\mathrm{D}\mathrm{80}}^{\star}$ mixing length and c_{0}=0.2.

The results obtained for l_{D80} and l_{BL89} are identical and close to the l_{R17} results with c_{0}=0.15 (not shown).

All simulations with the CCH02 set of parameters show reasonable results.
4.3 Stably stratified boundary layer (GABLS1)
Within the Global Energy and Water Exchanges (GEWEX) atmospheric boundary layer study (GABLS), idealized cases for stable surface boundary layers have been investigated (e.g., Cuxart et al., 2006). Such conditions are typical of areas covered with sea ice. Here we consider the GABLS1 case, whose technical description is available at http://turbulencia.uib.es/gabls/gabls1d_desc.pdf (last access: 20 January 2021). This experiment is particularly interesting as significant differences generally exist between solutions obtained from LES and singlecolumn simulations, for example when the Bougeault and Lacarrère (1989) length scale is used (e.g., Cuxart et al., 2006; Rodier et al., 2017). A largescale geostrophic wind is imposed as well as a cooling of the surface temperature θ_{s}(t) given by ${\mathit{\theta}}_{\mathrm{s}}\left(t\right)=\mathrm{263.5}\mathrm{0.25}(t/\mathrm{3600}\phantom{\rule{0.125em}{0ex}}\mathrm{s})$. The parameter values for this test are reported in Table 2 and the initial conditions are ${\mathbf{u}}_{\mathrm{h}}(z,t=\mathrm{0})={\mathbf{u}}_{\mathrm{G}}$, and
The solutions after 9 h of simulation are shown in Fig. 5 (left panels) for CBR00 parameter values and in Fig. 5 (right panels) for CCH02 parameter values. The reference solution is taken from Rodier et al. (2017) LESs. As expected, solutions based on a mixing length ignoring the contribution from the vertical shear exhibit a boundary layer that is too thick and a wind speed maximum located too high in altitude. Using a buoyancy and shearbased mixing length mitigates the issue and provides very good agreement with reference solutions when the CCH02 model constants are used. The best results are obtained for ${l}_{\mathrm{D}\mathrm{80}}^{\star}$ with c_{0}=0.2 and l_{R17} with c_{0}=0.15. Solutions obtained with the CBR00 model constants systematically predict larger turbulent kinetic energy and mixing lengths, resulting in large values of K_{s} in the first 100 m near the surface (not shown). The mismatch in terms of TKE is partially explained by the difference in boundary conditions since with CBR00 constants we have ${e}_{\mathrm{sfc}}=\mathrm{4.628}\phantom{\rule{0.25em}{0ex}}{u}_{\star}^{\mathrm{2}}$ while with CCH02 constants we get ${e}_{\mathrm{sfc}}=\mathrm{3.065}\phantom{\rule{0.25em}{0ex}}{u}_{\star}^{\mathrm{2}}$ from Eq. (7). Note that the proper calibration of the c_{0} constant jointly with the c_{ε} is the subject of several ongoing studies. Since our simulations reproduce the known sensitivity to those parameters, the ABL1d model could directly benefit from new findings on that topic.
The main outcomes are as follows:

The CCH02 set of parameters provides results of better quality than the CBR00 constants. For the sake of simplicity, we will retain only the CCH02 parameters for the numerical results shown in the remainder.

The buoyancy and verticalshearbased mixing lengths l_{R17} and ${l}_{\mathrm{D}\mathrm{80}}^{\star}$ are superior to the buoyancybased mixing lengths l_{D80} and l_{BL89} for stable boundary layers.
4.4 Winds across a midlatitude SST front
4.4.1 Setup and reference solutions
An idealized experiment particularly relevant for the coupling of the MABL with mesoscale oceanic eddies (and potentially submesoscale fronts) was initially suggested by Spall (2007) and then revised by Kilpatrick et al. (2014). More recently Ayet and Redelsperger (2019) derived an analytical model based on a similar setup. The geometry of the problem is twodimensional x–z with an SST front along the x axis:
where Δθ=3 K, L_{θ}=100 km, and $x\in [\mathrm{1800}\phantom{\rule{0.125em}{0ex}}\mathrm{km},\mathrm{1800}\phantom{\rule{0.125em}{0ex}}\mathrm{km}]$. As indicated in Table 2, a zonal geostrophic wind of 15 m s^{−1} is prescribed, balanced by a vertically homogeneous meridional pressure gradient. The wind thus flows over cold water before reaching a warm SST anomaly, which is 3 K warmer. We consider a dry case, in which the model is initialized ∀x with
where ${N}^{\mathrm{2}}={\mathrm{10}}^{\mathrm{4}}$ s^{−2} and θ^{ref}=288 K.
The velocities are systematically initialized with geostrophic winds. All simulations are run for 36 h when the flow reaches a quasiequilibrium state. For this configuration the reference solution is obtained from the mesoscale nonhydrostatic model (MesoNH) v5.3.0 (Lafore et al., 1998; Lac et al., 2018), where microphysics and radiation packages have not been activated. The horizontal resolution is Δx=1 km and the model is discretized with 91 vertical levels from the surface to 20 km height. The vertical resolution near the surface is Δz=10 m and around Δz=100 m at 2000 m height. The turbulence scheme is the 1.5order closure of Cuxart et al. (2000) in its onedimensional form with the l_{BL89} mixing length and CCH02 set of parameters. Sea surface fluxes are computed using the bulk parameterization COARE3.0, which is also available in NEMO from the AeroBulk package. As far as the ABL1d model is concerned, the top of the computational domain is z_{top}=2000 m and the vertical grid is stretched with a typical resolution of 20 m near the surface and 100 m near z=z_{top} with a first grid point located at z=10 m. In the horizontal, the resolution is Δx=6 km.
4.4.2 Numerical results
For this configuration, results will be mostly evaluated in terms of 10 m winds u_{10} and temperature θ_{10}. As an illustration of the type of result we get, we first compare the MesoNH solution and the ABL1d solution obtained with the l_{D80} and l_{BL89} mixing lengths in Fig. 6. It is worth noting that the MesoNH solution closely compares with the solution of Kilpatrick et al. (2014) (their Fig. 2) with a shallow boundary layer height (around 400 m) before the front and a thicker one (around 800 m) after the front where momentum mixing is enhanced. Over the front, as noted by Ayet and Redelsperger (2019) with a similar setup, the effect of advection is predominant for meridional winds, thus explaining the differences seen with the ABL1d simulations. Indeed with ABL1d, whatever the numerical options, the atmospheric column will locally adjust to the underlying oceanic conditions since horizontal advection is neglected. This explains the absence of horizontal lag when passing over the front in the ABL1d solution compared to the MesoNH solution. However, away from the SST front the solutions are very similar in terms of boundary layer height and vertical wind structure. In anticipation of a coupling with an oceanic model, the most important quantities to look at are the 10 m atmospheric variables rather than the full 3D vertical structure of the MABL. In Fig. 7, the 10 m wind components and temperature when the ABL1d model reaches a quasiequilibrium state are shown for different mixing length options, as well as the MesoNH results. First, the results obtained with the l_{R17} are very different from the expected behavior, and we will focus the discussion on other options. In terms of zonal 10 m wind the buoyancybased l_{BL89} and l_{D80} mixing lengths provide a good agreement with the MesoNH solution, which could be expected as the MesoNH solution has been generated using the l_{BL89} mixing length. As soon as the mixing length is a function of buoyancy and vertical shear (as is the case for ${l}_{\mathrm{D}\mathrm{80}}^{\star}$) the simulated winds are weaker because the boundary layer is thinner. This leads to improved results in the stably stratified case shown earlier, but in the present case, which is more representative of realistic configuration in the MABL, it leads to a mixing that is too weak. However, compared to the l_{R17} mixing length the ${l}_{\mathrm{D}\mathrm{80}}^{\star}$ still performs reasonably well but the winds on the warm side of the front are about 1 m s^{−1} weaker than the MesoNH winds for c_{0}=0.15 and become weaker and weaker as c_{0} increases.
The main outcomes are as follows:

In the frontal region the effect of horizontal advection is predominant and the ABL1d model cannot reproduce the horizontal lag seen in the reference solution when passing over the front.

The ABL1d model reproduces the downward momentum mixing mechanism correctly. The best results are obtained with the buoyancybased l_{BL89} and l_{D80} mixing lengths.

The l_{R17} mixing length will be discarded from the comparison.
Although relevant for the present study this 2D x–z setup is not fully representative of realistic conditions because the air column has time to adjust to the underlying oceanic state, which is kept frozen in time.
4.4.3 A singlecolumn version
An alternative to the x–z setup would be to formulate the test case as a Lagrangian advection of an air column over an SST front by prescribing a temporal evolution of sea surface temperature θ_{s}(t) as
In this case the air column does not necessarily have time to adjust to the underlying oceanic conditions. Initial conditions are the same as the ones of the 2D x–z case. For this test case we do not have a reference solution, but it is expected that the temporal evolution of the solution should be relatively similar to the spatial evolution in the MesoNH 2D x–z case studied in previous subsection. This can be seen from Fig. 8, where there is a clear similarity between the time vs. height sections obtained with the ABL1d simulations and the x vs. height sections shown for MesoNH in Fig. 6. The ABL1d solution shows a temporal lag analogous to the horizontal lag in the reference solution for the 2D x–z case. In addition to that, we also use this test case to investigate the sensitivity of the solutions to the bulk formulation and to the Newtonian relaxation which was absent in simulations discussed so far. We consider the COARE and IFS bulk formulations, which are relatively close to each other, to check the robustness of the results to small perturbations in surface fluxes. We also consider simulations with a relaxation of the temperature variable toward the initial condition with a fast relaxation timescale ${\mathit{\lambda}}_{\mathrm{s}}=\frac{\mathrm{1}}{\mathrm{6}\left[\mathrm{h}\right]}$ above the ABL and a slower one ${\mathit{\lambda}}_{\mathrm{s}}=\frac{\mathrm{1}}{\mathrm{48}\left[\mathrm{h}\right]}$ in the ABL. This is meant to check that the relaxation does not completely overwrite the physics of the coupling we aim to represent with the ABL1d model. Results from those sensitivity experiments are shown in Figs. 8 and 9. In particular in Fig. 9 the evolution of the 10 m winds across the SST front closely resembles the one shown in Fig. 7 (dashed red lines) for MesoNH. Moreover, the results in Fig. 9 are robust to a change of bulk formulation to compute the surface fluxes. Reassuringly, adding a relaxation toward largescale data which did not see the SST front does not deteriorate the realism of the solutions, as can be seen from Fig. 8c and f and 9 (gray lines).
The main outcomes are as follows:

The response of the ABL1d model to evolving oceanic conditions is not local in time (it shows a temporal lag).

The good representation of the downward momentum mixing process is not sensitive to the bulk formulation.

Adding a relaxation term toward largescale data does not deteriorate the realism of the solutions significantly.
Based on the results reported in this section, the best balance between efficiency and physical relevance is obtained when using the parameter values from CCH02 and the modified Deardorff (1980) mixinglength formulation l_{D80} or ${l}_{\mathrm{D}\mathrm{80}}^{\star}$. In particular we could imagine using the l_{D80} formulation over water and the ${l}_{\mathrm{D}\mathrm{80}}^{\star}$ formulation over sea ice.
Using atmosphereonly experiments, we have been focused so far on the good representation of the downward momentum mixing mechanism and of the stable boundary layers typical of areas covered with sea ice. In the following we check that those two aspects are still adequately represented in a realistic coupled NEMOABL1d simulation. This simulation will also be used to look at the wind–current interaction, which was left aside so far. We performed a 5year global simulation using the ORCA025 configuration. Details and illustrations are given hereafter.
5.1 Coupled NEMOABL1d configuration
We use here a global ORCA025 configuration at a 0.25^{∘} horizontal resolution (Barnier et al., 2006) with 75 vertical z levels forced by the ECWMF ERAInterim 6 h analysis (Dee et al., 2011). This configuration is identical to the one described in Couvelard et al. (2020) (see their Sect. 4.1.1). The ABL1dNEMO coupled simulation is carried out with the same numerical options as in a standard ASL forcing strategy. However, in the ABL coupling strategy, the twodimensional nearsurface air temperature, humidity, and winds used in the usual ASL forcing are replaced by threedimensional atmospheric variables sampled between the surface and 2000 m preprocessed following the different steps described in Appendix B. The largescale pressure gradient computed during the preprocessing is used as a geostrophic forcing for the ABL1d model dynamics. Threedimensional atmospheric variables are generated over the 2014–2018 period and vertically interpolated on 50 levels between 10 and 2000 m with a vertical resolution increasing with height. Grid resolution is about 20 m near the air–sea interface and reaches 70 m at the top of the ABL1d domain. The choice of a vertical extent of 2000 m and 50 vertical levels in the ABL1d model is somewhat arbitrary and the robustness of numerical results to these choices will be investigated in future studies. For the simulations presented here, the same horizontal grid and time step (Δt=1200 s) are chosen in the ABL1d and NEMO models. The options associated with the ABL coupling available through the NEMO standard name list are reported in Table 3, and a detailed profiling of the code is presented in Appendix D in order to assess the overhead associated with the ABL coupling strategy vs. ASL coupling strategy. This profiling shows that the overhead associated with the ABL1d (when using the l_{D80} mixing length) is on the order of 4 % and the one associated with the input part of the I/O operations is 5 %. Overall there is an increase of 9 % in elapsed time compared to the standard ASL forcing strategy.
5.2 Numerical results
In this section, we evaluate the ABL coupling strategy in a realistic context for a set of relevant metrics. The objective is not to conduct a thorough physical analysis of the numerical results but to illustrate the potential of the ABL coupling strategy and its proper implementation in NEMO. To evaluate our numerical results, we use standard metrics from the literature to quantify the wind–SST (a.k.a. thermal feedback effect), wind–currents (a.k.a. current feedback effect), and MABL–seaice couplings (e.g., Bryan et al., 2010; Renault et al., 2019b).
5.2.1 Thermal feedback effect
To quantify the surface wind response to SST, we show in Fig. 10a a global map of the temporal correlation between the highpassfiltered 10 m wind speed from the first vertical level in the ABL1d model and the SST. The same correlation is shown in Bryan et al. (2010) from satellite observations (their Fig. 1d) and from coupled numerical experiments between a 0.1^{∘} ocean and a 0.25^{∘} atmospheric model (their Fig. 1c). Consistent with observations and fully coupled models, the correlation obtained from the coupled NEMOABL1d simulation shows large positive correlations over regions like the Southern Ocean, Kuroshio, and Gulf Stream extensions as well as in the Gulf of Guinea. Correlations are, however, weaker than observations in the northern and equatorial Pacific between 90 and 180^{∘} W. As the thermal feedback strength is related to the ocean model resolution (Bryan et al., 2010), we can expect a better agreement with observations using a higher resolution configuration such as ORCA12 ($\mathrm{1}/\mathrm{12}{}^{\circ}$ resolution). This coupling sensitivity to the oceanic resolution will be addressed in a future study.
5.2.2 Current feedback effect
Other processes of interest are those related to the coupling between oceanic surface currents, wind stress, and wind. Such coupling is responsible for a dampening of the eddy mesoscale activity in the ocean. In Renault et al. (2019b), two coupling coefficients called s_{w} and s_{τ} are defined to quantify this effect. s_{τ} is a measure of the sink of energy from the eddies and fronts to the ABL and s_{w} quantifies the partial reenergization of the ocean by the wind response to the wind–current coupling. This reenergization is absent in the ASL forcing strategy, which results in an excessive dampening of the oceanic eddy mesoscale activity. In practice, s_{τ} (s_{w}) corresponds to the slope of the linear relationship between highpassfiltered surface current vorticity and surface windstress (wind) curl. Global maps of s_{τ} and s_{w} computed from our coupled NEMOABL1d global simulation are shown in Fig. 10. Large negative values of s_{τ} indicate an efficient dampening of the eddy mesoscale activity by the current feedback (i.e., a large sink of energy from the ocean to the atmosphere), and the large positive values of s_{w} indicate an efficient wind response and reenergization of the mesoscale currents. Our numerical experiment provides results very consistent with the results obtained from coupled simulations between NEMO and the WRF shown in Renault et al. (2019b) (their Fig. 1b for s_{τ} and 2c for s_{w}). As mentioned earlier, with an ASL forcing strategy we would systematically have s_{w}=0 and stronger s_{τ} values.
5.2.3 MABL and seaice coupling
The last illustration of our implementation presented in this section is the coupling of ABL1d with sea ice. As described in Sect. 4.3, sea ice generally induces a shallow stably stratified boundary layer due to the nearsurface air cooling. This increased vertical stability tends to reduce atmospheric turbulence, producing shallower ABL heights over sea ice. This relationship between seaice concentration and ABL height is clearly visible from Fig. 11 on both Arctic and Antarctic domains, where the ABL height follows a progressive decrease from about 800 to 200 m in the transition zone between the open ocean and fully icecovered regions. This coupling between the ABL and sea ice have important effects on nearsurface wind, temperature, and humidity, and consequently on seaice concentration evolution, which will need to be specifically assessed in future ABLbased studies.
6.1 Summary
A simplified atmospheric boundary layer (ABL) model has been developed and integrated to an oceanic model. This development is made with the objective to improve the representation of air–sea interactions in eddying oceanic models compared to the standard forcing strategy where the 10 m height atmospheric quantities are prescribed. For this preliminary study, the simplified ABL model takes the form of a singlecolumn model including a turbulence scheme coupled to each oceanic grid point. A crucial hypothesis is that the dominant process at the characteristic scale of the oceanic mesoscale is the socalled downward mixing process which stems from a modulation of atmospheric turbulence by sea surface temperature (SST) anomalies. Our approach can be seen as an extended bulk approach: instead of prescribing atmospheric quantities at 10 m to compute air–sea fluxes via an atmospheric surface layer (ASL) parameterization, atmospheric quantities in the first few hundred meters are used to constrain an ABL model which provides 10 m atmospheric values to the ASL parameterization. An important point is that our modeling strategy keeps the computational efficiency and flexibility inherent to oceanonly modeling. Indeed, the overhead generally observed in terms of computational cost compared to the usual ASL forcing strategy is roughly 10 %, and half of this overhead is due to I/O operations since the ABL model is constrained by 3D atmospheric data. In this paper the key components of such an approach have been described. This includes the largescale forcing strategy, the coupling with the ocean and sea ice and last but not least the ABL turbulence closure scheme based on a prognostic equation for the turbulence kinetic energy. The resulting simplified model, called ABL1d, has been tested for several boundarylayer regimes relevant to either ocean–atmosphere or seaice–atmosphere coupling. Results have systematically been evaluated against large eddy simulations (LESs). Furthermore we have investigated the behavior of the model to several parameters including the formulation of the mixing length and the turbulence model constants. First results from a global ABL1dNEMO configuration show an excellent behavior in terms of wind–SST twoway coupling. A first analysis of the impact of the coupling with ABL1d from a physical viewpoint is presented in Brivoal et al. (2020).
6.2 Future work
Now that an adequate computational framework and an efficient turbulent scheme that can be operated for a reasonable computational overhead have been developed, the next step is to investigate the relevance of the singlecolumn representation of the ABL selected for the present paper. Indeed, several studies have already shown that momentum vertical turbulent mixing, pressure gradient, Coriolis, and nonlinear advection are all important to the momentum balance in the marine atmospheric boundary layer at the vicinity of oceanic fronts (see for example Spall, 2007; Small et al., 2008; O'Neill et al., 2010). It is well known that the relative importance of those terms depends on the wind regime: for strong winds the vertical mixing is the dominant mechanism while for weak winds the pressure adjustment mechanism dominates. The current singlecolumn approximation is based on the assumption that the balance is dominated by vertical turbulent mixing, and the effect of other terms is roughly represented by the geostrophic guide and/or the nudging term. The test case presented in Sect. 4.4 clearly illustrates the limitations of a singlecolumn approach ignoring advective effects. However, before moving to more advanced formulations, our rationale was that the two main bottlenecks in terms of computational cost inherent to the ABL coupling strategy are the reading of 3D atmospheric data and the choice of ABL turbulent scheme. As a first step, we focused on those two aspects to assess whether or not our approach can be a viable option. Even if the justification of our model is not beyond reproach, it already brings an improvement compared to the ASL forcing strategy.
Several ways to improve the methodology presented here are currently under investigation. At a practical level, ways to lower the computational overhead due to I/O operations will be investigated using the parallel I/O capabilities provided by the XIOS library which is currently used in NEMO only for outputs. At a more fundamental level, the continuous formulation of the ABL1d model will be completed to improve the representation of the momentum balance by integrating the effect of horizontal advection and finescale pressure gradients. Increasing the complexity of the model should allow the impact of the nudging term on the ABL solutions to be lowered. In the event our approach turns out to be physically sound for a reasonable complexity it could be useful not only for offline oceanic simulations but also in coupled simulations to downscale the information from a lowresolution atmospheric component to a highresolution oceanic component. A standalone ABL model of intermediate complexity could also play a role in coupled data assimilation where the current practice is generally to assimilate data separately in the ocean and the atmosphere, ignoring the air–sea interactions, which results in inconsistencies at the air–sea interface in the initial conditions, causing initial shocks in the coupled forecasts (e.g., Mulholland et al., 2015). We wish to conclude this study by clarifying that the framework we have developed within NEMO is general enough to allow alternative approaches (e.g., via modeldriven empirical models) to be seamlessly tested and confronted with the ABL coupling strategy.
In this appendix, following the methodology of Redelsperger et al. (2001) reexpressed with our notations, we quickly recall how the surface boundary conditions for the turbulent kinetic energy (TKE) and for the mixing lengths l_{m} and l_{ε} are determined via a matching between the subgrid turbulence scheme and the surfacelayer theory. The simplest form of atmospheric surface layer (ASL) theory, namely for neutral stratification, is considered. Under the quasiequilibrium hypothesis the evolution Eq. (6) for the TKE reduces to the equilibrium between turbulence production and dissipation ${K}_{\mathrm{m}}\Vert {\partial}_{z}{\mathbf{u}}_{\mathrm{h}}{\Vert}^{\mathrm{2}}=\frac{{c}_{\mathit{\epsilon}}}{{l}_{\mathit{\epsilon}}}{e}^{\mathrm{3}/\mathrm{2}}$ which, combined with ${K}_{\mathrm{m}}={C}_{\mathrm{m}}{l}_{\mathrm{m}}\sqrt{e}$, leads to
The similarity theory for the ASL in the neutral case is such that $\Vert {\partial}_{z}{\mathbf{u}}_{\mathrm{h}}\Vert =\frac{{u}_{\star}}{\mathit{\kappa}(z+{z}_{\mathrm{0}})}$ with κ the Von Karman constant and z_{0} the roughness length which can be combined with Eq. (A1) to get that in the surface layer:
with ${u}_{\star}^{\mathrm{2}}=\Vert \u2329{\mathbf{u}}_{\mathrm{h}}^{\prime}{w}^{\prime}\u232a\Vert $ the friction velocity. Moreover, enforcing the consistency between the eddy diffusivity for momentum given by the ASL theory (${K}_{\mathrm{m}}=\mathit{\kappa}{u}_{\star}(z+{z}_{\mathrm{0}})$) and the one given by the TKE closure (${K}_{\mathrm{m}}={C}_{\mathrm{m}}{l}_{\mathrm{m}}\sqrt{e}$) leads to
We thus have two relations (A2) and (A3) for three unknowns (e, l_{m}, and l_{ε}). At this point our derivation will differ from Redelsperger et al. (2001) as we will assume that ${l}_{\mathrm{m}}={l}_{\mathit{\epsilon}}={L}_{\mathrm{sfc}}$ in the ASL.
Under this assumption, combining Eqs. (A2) and (A3) we easily obtain
for z≤δ_{asl}, where δ_{asl} corresponds to the extent of the ASL.
The expression of L_{sfc}(z) for z≤δ_{asl} is also used as a constraint to define the weighted average needed to determine l_{m} from l_{up} and l_{dwn}:
(equivalent to Eq. 12a). In the ASL we further assume that l_{dwn}≈0 and l_{up}≈δ_{asl}; for l_{m}(z=δ_{asl}) to be consistent with L_{sfc}(z=δ_{asl}) we should have
Considering that δ_{asl}≫z_{0}, Eq. (A4) is satisfied for
Considering the CBR00 model constants we obtain $a=\mathrm{1.4796}\approx \mathrm{3}/\mathrm{2}$ and $a=\mathrm{0.860834}\approx \mathrm{6}/\mathrm{7}$ for the CCH02 constants (see Table 1).
B1 Altitude of IFS vertical levels
The ABL1d model is discretized on fixed in time and space geopotential levels while the IFS model uses a pressurebased sigma coordinate. A first step is to recover the altitude associated with each sigma level. The pressure p_{k+½} defined at cell interfaces between two successive vertical layers is given by
where A_{k+½} (Pa) and B_{k+½} (dimensionless) are constants given by a smooth analytical function defining the vertical grid stretching. Typical values of the altitude of grid points in the vertical for a standard 60level grid (L60) and a surface pressure of 1013 hPa are given in Table B1. Once the values of p_{k+½} and p_{s} are known, the altitude of cell interfaces can be computed by integrating the hydrostatic equilibrium
vertically. In Eq. (B1), ϕ is the geopotential, T_{v} the virtual temperature, and R_{d} the specific gas constant for dry air. At a discrete level we get
which gives
Once the layer thicknesses $\mathrm{e}\mathrm{3}{\mathrm{t}}_{k}^{\mathrm{ifs}}$ are known, horizontal wind components, potential temperature, and specific humidity can be interpolated on the ABL1d vertical levels. Under the constraint that ${\int}_{{z}_{\mathrm{sfc}}}^{{z}_{\mathrm{top}}}{\mathit{\psi}}^{\mathrm{ifs}}\left(z\right)\phantom{\rule{0.25em}{0ex}}\mathrm{d}z={\int}_{{z}_{\mathrm{sfc}}}^{{z}_{\mathrm{top}}}\mathit{\psi}\left(z\right)\phantom{\rule{0.25em}{0ex}}\mathrm{d}z$ for any IFS quantity ψ^{ifs} to be interpolated. Wind components are interpolated using a fourthorder compact scheme while tracers are interpolated using a WENOlike PPM scheme (Alexander Shchepetkin, personal communication, 2001) which is monotonic.
B2 Filtering in the presence of boundaries
Because of the IFS numerical formulation and of the postprocessing of output data, the solutions sometimes contain highfrequency oscillations at the vicinity of the landsea interface. This problem is further compounded when the nearshore topography is steep. The atmospheric fields over water thus need to be smoothed horizontally to specifically remove the 2Δx noise. We use a standard twodimensional Shapiro filter which, in the absence of lateral boundaries, can be formulated as
where ${\mathit{\delta}}_{i+\mathrm{1}/\mathrm{2},j}^{\left(x\right)}={\mathit{\psi}}_{i+\mathrm{1},j}{\mathit{\psi}}_{i,j}$ and ${\mathit{\delta}}_{i,j+\mathrm{1}/\mathrm{2}}^{(y,\star )}={\mathit{\psi}}_{i,j+\mathrm{1}}^{\star}{\mathit{\psi}}_{i,j}^{\star}$. The amplification factor associated with this filter is
which guarantees that one iteration of the filter is sufficient to remove the gridscale noise since ${\mathcal{A}}_{\mathrm{shap}}(\mathit{\pi},\mathit{\pi})={\mathcal{A}}_{\mathrm{shap}}(\mathit{\pi},{\mathit{\theta}}_{y})={\mathcal{A}}_{\mathrm{shap}}({\mathit{\theta}}_{x},\mathit{\pi})=\mathrm{0}$ and that 𝒜_{shap}≤1 (i.e., no waves are amplified). In the presence of solid boundaries we would like to retain those properties as much as possible. A straightforward approach would be to impose a nogradient condition at the coast, i.e., ${\mathit{\delta}}_{i+\mathrm{1}/\mathrm{2},j}^{\left(x\right)}=\mathrm{0}$ as soon as ${\mathrm{tmask}}_{i+\mathrm{1},j}\times {\mathrm{tmask}}_{i,j}=\mathrm{0}$ (${\mathit{\delta}}_{i,j+\mathrm{1}/\mathrm{2}}^{(y,\star )}=\mathrm{0}$ as soon as ${\mathrm{tmask}}_{i,j+\mathrm{1}}\times {\mathrm{tmask}}_{i,j}=\mathrm{0}$), with tmask the indicator function equal to 1 over water and 0 over land. Let us also consider the following alternative boundary conditions
and similar in the y direction. We do not elaborate on this choice but it can be shown theoretically that boundary conditions (B2) provide a better control of gridscale noise near the coast. To illustrate this point, in Fig. B1 the surface pressure gradients are shown for different boundary conditions. In particular it can be seen near the coast that the nogradient boundary condition (panels b and e) leaves some artificial patterns in gradients, especially in the Peru–Chile current system, while the boundary condition (B2) efficiently mitigates this issue. Note that it is particularly essential to make sure that the surface pressure field is sufficiently smooth because gradients of this field are used to compute geostrophic winds which are important for the largescale forcing of the ABL1d model.
B3 Largescale pressure gradient computation
The last aspect of the preprocessing of atmospheric data we would like to discuss is the computation of the largescale pressure gradient (or equivalently of the geostrophic wind components) The objective is to estimate the following terms:
where (⋅)_{z} denotes a gradient along constant geopotential height. Using the hydrostatic balance we have $\frac{\mathrm{1}}{{\mathit{\rho}}_{\mathrm{a}}}=g({\partial}_{z}p{)}^{\mathrm{1}}$, which leads to
Assuming a generalized vertical coordinate $s=s(x,y)$ the computation of gradients along constant height is not straightforward since $({\partial}_{x}p{)}_{z}=({\partial}_{x}p{)}_{\mathrm{s}}\left({\partial}_{z}p\right)({\partial}_{x}z{)}_{\mathrm{s}}$ leading to
In the particular case of the IFS coordinate s we have
and (∂_{x}z)_{s} can be estimated after integrating the hydrostatic balance. Starting from the layer interface height ${z}_{i,j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{ifs}}$, surface pressure (p_{s})_{i,j}, and parameter values ${A}_{k},{B}_{k},{A}_{k+\mathrm{1}/\mathrm{2}},{B}_{k+\mathrm{1}/\mathrm{2}}$ the different steps are the following:

compute Δx_{i,j} and Δy_{i,j} from latitudes and longitudes;

compute horizontal gradients ∂_{x}p_{s} and ∂_{y}p_{s} for surface pressure
$$\begin{array}{l}{\displaystyle}{\mathrm{FX}}_{i+\mathrm{1}/\mathrm{2},j}={\displaystyle \frac{\mathrm{2}\left\{({p}_{\mathrm{s}}{)}_{i+\mathrm{1},j}({p}_{\mathrm{s}}{)}_{i,j}\right\}}{\mathrm{\Delta}{x}_{i,j}+\mathrm{\Delta}{x}_{i+\mathrm{1},j}}},\\ {\displaystyle}\phantom{\rule{1em}{0ex}}{\mathrm{FY}}_{i,j+\mathrm{1}/\mathrm{2}}={\displaystyle \frac{\mathrm{2}\left\{({p}_{\mathrm{s}}{)}_{i,j+\mathrm{1}}({p}_{\mathrm{s}}{)}_{i,j}\right\}}{\mathrm{\Delta}{y}_{i,j}+\mathrm{\Delta}{y}_{i,j+\mathrm{1}}}};\end{array}$$ 
compute horizontal gradients (∂_{x}z)_{s} and (∂_{y}z)_{s}
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}{\mathrm{dZx}}_{i+\mathrm{1}/\mathrm{2},j,k}=\\ {\displaystyle}& {\displaystyle}\phantom{\rule{1em}{0ex}}{\displaystyle \frac{{z}_{i+\mathrm{1},j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}{z}_{i,j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}+{z}_{i+\mathrm{1},j,k\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}{z}_{i,j,k\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}}{\mathrm{\Delta}{x}_{i,j}+\mathrm{\Delta}{x}_{i+\mathrm{1},j}}},\\ {\displaystyle}& {\displaystyle}{\mathrm{dZy}}_{i,j+\mathrm{1}/\mathrm{2},k}=\\ {\displaystyle}& {\displaystyle}\phantom{\rule{1em}{0ex}}{\displaystyle \frac{{z}_{i,j+\mathrm{1},k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}{z}_{i,j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}+{z}_{i,j+\mathrm{1},k\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}{z}_{i,j,k\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}}{\mathrm{\Delta}{y}_{i,j}+\mathrm{\Delta}{y}_{i,j+\mathrm{1}}}};\end{array}$$ 
compute $\left({\partial}_{z}p{)}^{\mathrm{1}}\right({\partial}_{x}p{)}_{\mathrm{s}}$ via (B4)
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}{\mathrm{wrkX}}_{i,j,k}=\\ {\displaystyle}& {\displaystyle}\phantom{\rule{1em}{0ex}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\phantom{\rule{0.25em}{0ex}}{\displaystyle \frac{{B}_{k}\left({z}_{i,j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}{z}_{i,j,k\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}\right)\left({\mathrm{FX}}_{i+\mathrm{1}/\mathrm{2},j}+{\mathrm{FX}}_{i\mathrm{1}/\mathrm{2},j}\right)}{\left({p}_{s}{)}_{i,j}\right({B}_{k+\mathrm{1}/\mathrm{2}}{B}_{k\mathrm{1}/\mathrm{2}})+({A}_{k+\mathrm{1}/\mathrm{2}}{A}_{k\mathrm{1}/\mathrm{2}})}},\\ {\displaystyle}& {\displaystyle}{\mathrm{wrkY}}_{i,j,k}=\\ {\displaystyle}& {\displaystyle}\phantom{\rule{1em}{0ex}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\phantom{\rule{0.25em}{0ex}}{\displaystyle \frac{{B}_{k}\left({z}_{i,j,k+\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}{z}_{i,j,k\mathrm{1}/\mathrm{2}}^{\mathrm{cep}}\right)\left({\mathrm{FY}}_{i+\mathrm{1}/\mathrm{2},j}+{\mathrm{FY}}_{i\mathrm{1}/\mathrm{2},j}\right)}{\left({p}_{s}{)}_{i,j}\right({B}_{k+\mathrm{1}/\mathrm{2}}{B}_{k\mathrm{1}/\mathrm{2}})+({A}_{k+\mathrm{1}/\mathrm{2}}{A}_{k\mathrm{1}/\mathrm{2}})}};\end{array}$$ 
finalize (we get a minus sign in ${R}_{\mathrm{LS}}^{u}$ because the grid in the y direction is flipped in the raw data)
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}({R}_{\mathrm{LS}}^{u}{)}_{i,j,k}=\\ {\displaystyle}& {\displaystyle}\phantom{\rule{1em}{0ex}}g\left({\mathrm{wrkY}}_{i,j,k}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathrm{dZy}}_{i,j+\mathrm{1}/\mathrm{2},k}+{\mathrm{dZy}}_{i,j\mathrm{1}/\mathrm{2},k})\right),\\ {\displaystyle}& {\displaystyle}({R}_{\mathrm{LS}}^{v}{)}_{i,j,k}=\\ {\displaystyle}& {\displaystyle}\phantom{\rule{1em}{0ex}}g\left({\mathrm{wrkX}}_{i,j,k}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathrm{dZx}}_{i+\mathrm{1}/\mathrm{2},j,k}+{\mathrm{dZx}}_{i\mathrm{1}/\mathrm{2},j,k})\right).\end{array}$$
In the following we describe the discrete algorithms used to provide the mixing lengths l_{up} and l_{dwn} given in Sect. 3.2.2. Four different ways to compute those quantities have been implemented in the ABL1d model.
C1 Bougeault and Lacarrère (1989) length scale
The Bougeault and Lacarrère (1989) mixing length defines l_{up} and l_{dwn} as
By construction such mixing lengths are bounded by the distance to the bottom and the top of the computational domain and revert to the Deardorff (1980) length scale (i.e., ${l}_{\mathrm{up}}={l}_{\mathrm{dwn}}=\sqrt{\mathrm{2}e\left(z\right)/{N}^{\mathrm{2}}}$) for N^{2}=cste. An objective is to also satisfy this last property at a discrete level. Considering a simple trapezoidal rule to approximate the integral in Eq. (C1) over each grid cell, the procedure for the computation of ${l}_{\mathrm{up}}\left({z}_{k+\mathrm{1}/\mathrm{2}}\right)$ is given in Algorithm 1. In the case ${N}^{\mathrm{2}}\left({z}_{p+\mathrm{1}/\mathrm{2}}\right)={N}^{\mathrm{2}}\left({z}_{p\mathrm{1}/\mathrm{2}}\right)={N}_{\mathrm{cst}}^{\mathrm{2}}\phantom{\rule{0.25em}{0ex}}(\forall p)$, Algorithm 1 gives the following sequence:
As soon as $\mathrm{FC}\left({z}_{p+\mathrm{1}/\mathrm{2}}\right)$ changes sign we stop the procedure because l_{up} such that $e\left({z}_{k+\mathrm{1}/\mathrm{2}}\right)+{N}_{\mathrm{cst}}^{\mathrm{2}}{l}_{\mathrm{up}}^{\mathrm{2}}=\mathrm{0}$, which corresponds to the Deardorff (1980) length scale, has been found. We note l_{BL89} the mixing length corresponding to the Bougeault and Lacarrère (1989) algorithm.
C2 Adaptation of NEMO's length scale
The standard NEMO algorithm (Sect. 10.1.3 in Madec, 2012) is much easier to discretize. As a first step the Deardorff (1980) length scale l_{D80} is computed at cell interfaces, such that
with ${N}_{\mathit{\epsilon}}^{\mathrm{2}}$ the minimum stratification allowed whose value is set to the smallest positive real computer value. The vertical gradients of l_{D80} are then limited such that they stay smaller than the variations of height. This amounts to compute l_{up} and l_{dwn} as
with e3t_{k} the thickness of vertical layer k (Fig. 2). The resulting mixing length is simply referred to as l_{D80}.
C3 Rodier et al. (2017) length scale
Recently, Rodier et al. (2017) proposed to add a shearrelated term to Eq. (C1) such that the definition of l_{up} and l_{dwn} becomes
where c_{0} is a parameter whose value should be smaller than $\sqrt{{C}_{\mathrm{m}}/{c}_{\mathit{\epsilon}}}$. At a discrete level, the FC function in Algorithm 1 is replaced by
This mixing length will be referred to as l_{R17}.
C4 A local buoyancy and shearbased length scale
For the sake of computational efficiency, we have derived a local version of the Rodier et al. (2017) length scale (denoted as ${l}_{\mathrm{D}\mathrm{80}}^{\star}$) which is original to the present paper:
Once ${l}_{\mathrm{D}\mathrm{80}}^{\star}$ has been computed at cell interfaces $z={z}_{k}+\mathrm{1}/\mathrm{2}$ we apply the limitations (C2) and (C3) as in the NEMO algorithm.
To finalize our description of the implementation of the simplified atmospheric boundary layer model in NEMO, we assess in this appendix the computational efficiency of our approach. We compare the performance of two simulations: one with a coupling with the ABL1d model (with 50 vertical levels) which requires reading 3D atmospheric data in input files, and one with a standard ASL forcing strategy which necessitates reading only 2D atmospheric data. For the coupling with ABL1d, we consider the l_{D80} mixing length which gave robustly good results across the different numerical tests investigated earlier in the paper. The simulations are performed with NEMO version 4.0 for the ORCA025 configuration previously described on 128 cores (Intel(R) Xeon(R) E5 processors 2.6 GHz) compiled with ifort (v13.0.1) using the “$\mathrm{i}\mathrm{4}\phantom{\rule{0.25em}{0ex}}\mathrm{r}\mathrm{8}\phantom{\rule{0.25em}{0ex}}\mathrm{O}\mathrm{3}\phantom{\rule{0.25em}{0ex}}\mathrm{fp}\mathrm{model}\phantom{\rule{0.25em}{0ex}}\mathrm{precise}\phantom{\rule{0.25em}{0ex}}\mathrm{fno}\mathrm{alias}$” options. The I/Os are handled via the Lustre file system. Each MPI subdomain has 80×130 points in the horizontal and 75 points in the vertical. The various reports given below have been obtained from a builtin NEMO code instrumentation dedicated to calculation measurement (e.g., Maisonnave and Masson, 2019). As mentioned earlier, the outputs are done using the parallel I/O capabilities provided by the XIOS library. Thanks to XIOS, we do not expect any significant difference between the two simulations regarding the cost of output operations. However, the use of XIOS to handle input operations is still under development, and because of the significant amount of data to read in the ABL coupling strategy it makes sense to assess the associated overhead. We ran the ASL forced and ABL coupled NEMO simulations for 20 d such that the cost of the initialization step is no longer visible in the averaged cost per time step. Moreover, for the two simulations, the atmospheric data necessary for the computation of the turbulent components of air–sea fluxes are provided every 6 h.
We first show in Fig. D1 the elapsed time for each time step over the first 48 h of the simulations with different ways to specify the surface fluxes. For most time steps, the overhead associated with the ABL1d when using the l_{D80} mixing length is very small (on the order of 4 %); however, every 18 time steps (i.e., every 6 h), there is a larger overhead due to the input part of the I/O operations. To further refine our assessment, we report in Table D1 the elapsed and CPU time spent on average over all the processors in the 11 most expansive sections of the code. As expected, the CPU time is not significantly affected by the ABL1d model (increase of 4 %), but the elapsed time is increased by about 9 % because of the time spent in waiting for I/O operations. The overhead associated with input operations could be mitigated by reducing the number of vertical levels in the ABL1d model (we used 50 levels here to get an upper bound on the computational overhead) and either by using XIOS to handle input operations or by running ABL1d in detached mode as explained in Sect. 3.4 such that the time spent reading input files is covered by actual computations. Nonetheless the small increase in CPU time leaves room for further improvements of the ABL model to relax the horizontal homogeneity assumption.
The changes to the NEMO code have been made on the standard NEMO code (release 4.0). The code can be downloaded from the NEMO website (http://www.nemoocean.eu/, last access: 23 June 2020). The NEMO code modified to include the ABL1d model is available in the Zenodo archive (https://doi.org/10.5281/zenodo.3904518, Lemarié and Samson, 2020). The name lists and data used to produce the figures are also available in the Zenodo archive.
FL wrote the paper with the help of all the coauthors. FL, GS, and GM designed and developed a preliminary version of the ABL1d model within the NEMO 3.6 stable version. This original code was then ported to NEMO release 4.0. JLR and HG provided inputs in the design of the TKE closure scheme and of the numerical experiments. FL carried out the idealized numerical experiments, GS the realistic experiments, and JLR the MesoNH simulations.
The authors declare that they have no conflict of interest.
We thank Anton Beljaars and one anonymous reviewer whose efforts helped to improve earlier versions of this paper. We also thank Pascal Marquet and Sébastien Masson for useful discussions.
Florian Lemarié and JeanLuc Redelsperger also acknowledge the support by MercatorOcean and the Copernicus Marine Environment Monitoring Service (CMEMS) through contract 22GLOHR – Lot 2 (Highresolution ocean, waves, atmosphere interaction).
This research has been supported by the H2020 European Institute of Innovation and Technology (IMMERSE (grant no. 821926)).
This paper was edited by Christina McCluskey and reviewed by Anton Beljaars and one anonymous referee.
Abel, R.: Aspects of airsea interaction in atmosphereocean models, PhD thesis, Kiel University, 2018. a
Andren, A., Brown, A. R., Mason, P. J., Graf, J., Schumann, U., Moeng, C.H., and Nieuwstadt, F. T. M.: Largeeddy simulation of a neutrally stratified boundary layer: A comparison of four computer codes, Q. J. Roy. Meteor. Soc., 120, 1457–1484, 1994. a, b, c, d
Ayet, A. and Redelsperger, J.L.: An analytical study of the atmospheric boundary layer flow and divergence over a SST front, Q. J. Roy. Meteor. Soc., 145, 2549–2567, https://doi.org/10.1002/qj.3578, 2019. a, b, c, d
Baklanov, A. A., Grisogono, B., Bornstein, R., Mahrt, L., Zilitinkevich, S. S., Taylor, P., Larsen, S. E., Rotach, M. W., and Fernando, H. J. S.: The Nature, Theory, and Modeling of Atmospheric Planetary Boundary Layers, B. Am. Meteorol. Soc., 92, 123–128, https://doi.org/10.1175/2010BAMS2797.1, 2011. a
Barnier, B., Siefridt, L., and Marchesiello, P.: Thermal forcing for a global ocean circulation model using a threeyear climatology of ECMWF analyses, J. Mar. Res., 6, 363–380, https://doi.org/10.1016/09247963(94)000349, 1995. a
Barnier, B., Madec, G., Penduff, T., Molines, J.M., Treguier, A.M., Le Sommer, J., Beckmann, A., Biastoch, A., Böning, C., Dengg, J., Derval, C., Durand, E., Gulev, S., Remy, E., Talandier, C., Theetten, S., Maltrud, M., McClean, J., and De Cuevas, B.: Impact of partial steps and momentum advection schemes in a global ocean circulation model at eddypermitting resolution, Ocean Dynam., 56, 543–567, https://doi.org/10.1007/s1023600600821, 2006. a
Bazile, E., Marquet, P., Bouteloup, Y., and Bouyssel, F.: The Turbulent Kinetic Energy (TKE) scheme in the NWP models at Meteo France, in: Workshop on Workshop on Diurnal cycles and the stable boundary layer, 7–10 November 2011, ECMWF, Shinfield Park, Reading, 127–135, available at: https://www.ecmwf.int/node/8006 (last access: 20 January 2021), 2012. a, b
Beljaars, A.: The parametrization of surface fluxes in largescale models under free convection, Q. J. Roy. Meteor. Soc., 121, 255–270, 1995. a
Beljaars, A., Dutra, E., Balsamo, G., and Lemarié, F.: On the numerical stability of surface–atmosphere coupling in weather and climate models, Geosci. Model Dev., 10, 977–989, https://doi.org/10.5194/gmd109772017, 2017. a, b
Bielli, S., Douville, H., and Pohl, B.: Understanding the West African monsoon variability and its remote effects: An illustration of the grid point nudging methodology, Clim. Dynam., 35, 159–174, https://doi.org/10.1007/s0038200906678, 2009. a
Bougeault, P. and André, J.C.: On the Stability of the THIRDOrder Turbulence Closure for the Modeling of the StratocumulusTopped Boundary Layer, J. Atmos. Sci., 43, 1574–1581, https://doi.org/10.1175/15200469(1986)043<1574:OTSOTT>2.0.CO;2, 1986. a
Bougeault, P. and Lacarrère, P.: Parameterization of orographyinduced turbulence in a mesobetascale model, Mon. Weather Rev., 117, 1872–1890, 1989. a, b, c, d, e, f
Bourras, D., Reverdin, G., Giordani, H., and Caniaux, G.: Response of the atmospheric boundary layer to a mesoscale oceanic eddy in the northeast Atlantic, J. Geophys. Res., 109, https://doi.org/10.1029/2004JD004799, 2004. a
Brivoal, T., Samson, G., Giordani, H., BourdalléBadie, R., Lemarié, F., and Madec, G.: Impact of the current feedback on kinetic energy over the NorthEast Atlantic from a coupled ocean/atmospheric boundary layer model, Ocean Sci. Discuss. [preprint], https://doi.org/10.5194/os202078, in review, 2020. a
Brodeau, L., Barnier, B., Gulev, S. K., and Woods, C.: Climatologically Significant Effects of Some Approximations in the Bulk Parameterizations of Turbulent Air–Sea Fluxes, J. Phys. Oceanogr., 47, 5–28, https://doi.org/10.1175/JPOD160169.1, 2017. a
Bryan, F. O., Tomas, R., Dennis, J. M., Chelton, D. B., Loeb, N. G., and McClean, J. L.: Frontal Scale Air–Sea Interaction in HighResolution Coupled Climate Models, J. Climate, 23, 6277–6291, https://doi.org/10.1175/2010JCLI3665.1, 2010. a, b, c, d
Burchard, H.: Energyconserving discretisation of turbulent shear and buoyancy production, Ocean Modell., 4, 347–361, https://doi.org/10.1016/S14635003(02)000094, 2002a. a
Burchard, H.: Applied Turbulence Modelling in Marine Waters, Lecture Notes in Earth Sciences, Springer Berlin Heidelberg, 2002b. a, b
Businger, J. and Shaw, W.: The response of the marine boundary layer to mesoscale variations in seasurface temperature, Dynam. Atmos. Oceans, 8, 267–281, 1984. a
Chelton, D. B. and Xie, S.P.: Coupled oceanatmosphere interaction at oceanic mesoscales, Oceanography, 23, 52–69, 2010. a, b, c
Cheng, Y., Canuto, V. M., and Howard, A. M.: An improved model for the turbulent PBL, J. Atmos. Sci., 59, 1550–1565, 2002. a, b, c, d, e
Couvelard, X., Lemarié, F., Samson, G., Redelsperger, J.L., Ardhuin, F., Benshila, R., and Madec, G.: Development of a twowaycoupled ocean–wave model: assessment on a global NEMO(v3.6)–WW3(v6.02) coupled configuration, Geosci. Model Dev., 13, 3067–3090, https://doi.org/10.5194/gmd1330672020, 2020. a
Craig, A., Valcke, S., and Coquart, L.: Development and performance of a new version of the OASIS coupler, OASIS3MCT_3.0, Geosci. Model Dev., 10, 3297–3308, https://doi.org/10.5194/gmd1032972017, 2017. a
Cuxart, J., Bougeault, P., and Redelsperger, J.L.: A turbulence scheme allowing for mesoscale and largeeddy simulations, Q. J. Roy. Meteor. Soc., 126, 1–30, 2000. a, b, c, d, e, f, g, h, i
Cuxart, J., Holtslag, A. A. M., Beare, R. J., Bazile, E., Beljaars, A., Cheng, A., Conangla, L., Ek, M., Freedman, F., Hamdi, R., Kerstein, A., Kitagawa, H., Lenderink, G., Lewellen, D., Mailhot, J., Mauritsen, T., Perov, V., Schayes, G., Steeneveld, G.J., Svensson, G., Taylor, P., Weng, W., Wunsch, S., and Xu, K.M.: SingleColumn Model Intercomparison for a Stably Stratified Atmospheric Boundary Layer, Bound.Lay. Meteorol., 118, 273–303, 2006. a, b, c
Deardorff, J. W.: Threedimensional numerical study of turbulence in an entraining mixed layer, Bound.Lay. Meteorol., 7, 199–226, 1974. a
Deardorff, J. W.: Stratocumuluscapped mixing layers derived from a three dimensional model, Bound.Lay. Meteorol., 18, 495–527, 1980. a, b, c, d, e, f
Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., MongeSanz, B. M., Morcrette, J.J., Park, B.K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.N., and Vitart, F.: The ERAInterim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597, https://doi.org/10.1002/qj.828, 2011. a
Deleersnijder, E., Beckers, J.M., Campin, J.M., El Mohajir, M., Fichefet, T., and Luyten, P.: Some mathematical problems associated with the development and use of marine models, in: The Mathematics of Models for Climatology and Environment, edited by: Díaz, J. I., Springer Berlin Heidelberg, Berlin, Heidelberg, 39–86, 1997. a
Deremble, B., Wienders, N., and Dewar, W. K.: CheapAML: A Simple, Atmospheric Boundary Layer Model for Use in OceanOnly Model Calculations, Mon. Weather Rev., 141, 809–821, 2013. a
Deshayes, J., Tréguier, A.M., Barnier, B., Lecointre, A., Sommer, J. L., Molines, J.M., Penduff, T., BourdalléBadie, R., Drillet, Y., Garric, G., Benshila, R., Madec, G., Biastoch, A., Böning, C. W., Scheinert, M., Coward, A. C., and Hirschi, J. J.M.: Oceanic hindcast simulations at high resolution suggest that the Atlantic MOC is bistable, J. Geophys. Res., 40, 3069–3073, https://doi.org/10.1002/grl.50534, wOS:000321951300034, 2013. a
Dewar, W. K. and Flierl, G. R.: Some Effects of the Wind on Rings, J. Phys. Oceanogr., 17, 1653–1667, 1987. a
Frenger, I., Gruber, N., Knutti, R., and Munnich, M.: Southern Ocean Eddies Affect Local Weather, Nat. Geosci., 6, 608–612, 2013. a
Giordani, H., Planton, S., Bénech, B., and Kwon, B.H.: Atmospheric boundary layer response to sea surface temperatures during the Semaphore experiment, J. Geophys. Res., 103, 25047–25060, 1998. a
Giordani, H., Caniaux, G., and Prieur, L.: A Simplified 3D Oceanic Model Assimilating Geostrophic Currents: Application to the POMME Experiment, J. Phys. Oceanogr., 35, 628–644, https://doi.org/10.1175/JPO2724.1, 2005. a
Haney, R. L.: Surface Thermal Boundary Condition for Ocean Circulation Models, J. Phys. Oceanogr., 1, 241–248, https://doi.org/10.1175/15200485(1971)001<0241:STBCFO>2.0.CO;2, 1971. a
Hogg, A., Dewar, W., Berloff, P., Kravtsov, S., and Hutchinson, D. K.: The effects of mesoscale oceanatmosphere coupling on the largescale ocean circulation, J. Climate, 22, 4066–4082, 2009. a
Hourdin, F., Couvreux, F., and Menut, L.: Parameterization of the Dry Convective Boundary Layer Based on a Mass Flux Representation of Thermals, J. Atmos. Sci., 59, 1105–1123, https://doi.org/10.1175/15200469(2002)059<1105:POTDCB>2.0.CO;2, 2002. a
Kilpatrick, T., Schneider, N., and Qiu, B.: Boundary Layer Convergence Induced by Strong Winds across a Midlatitude SST Front, J. Climate, 27, 1698–1718, 2014. a, b, c
Kleeman, R. and Power, S.: A Simple Atmospheric Model of Surface Heat Flux for Use in Ocean Modeling Studies, J. Phys. Oceanogr., 25, 92–105, https://doi.org/10.1175/15200485(1995)025<0092:ASAMOS>2.0.CO;2, 1995. a
Lac, C., Chaboureau, J.P., Masson, V., Pinty, J.P., Tulet, P., Escobar, J., Leriche, M., Barthe, C., Aouizerats, B., Augros, C., Aumond, P., Auguste, F., Bechtold, P., Berthet, S., Bielli, S., Bosseur, F., Caumont, O., Cohard, J.M., Colin, J., Couvreux, F., Cuxart, J., Delautier, G., Dauhut, T., Ducrocq, V., Filippi, J.B., Gazen, D., Geoffroy, O., Gheusi, F., Honnert, R., Lafore, J.P., Lebeaupin Brossier, C., Libois, Q., Lunet, T., Mari, C., Maric, T., Mascart, P., Mogé, M., Molinié, G., Nuissier, O., Pantillon, F., Peyrillé, P., Pergaud, J., Perraud, E., Pianezze, J., Redelsperger, J.L., Ricard, D., Richard, E., Riette, S., Rodier, Q., Schoetter, R., Seyfried, L., Stein, J., Suhre, K., Taufour, M., Thouron, O., Turner, S., Verrelle, A., Vié, B., Visentin, F., Vionnet, V., and Wautelet, P.: Overview of the MesoNH model version 5.4 and its applications, Geosci. Model Dev., 11, 1929–1969, https://doi.org/10.5194/gmd1119292018, 2018. a
Lafore, J. P., Stein, J., Asencio, N., Bougeault, P., Ducrocq, V., Duron, J., Fischer, C., Héreil, P., Mascart, P., Masson, V., Pinty, J. P., Redelsperger, J. L., Richard, E., and VilàGuerau de Arellano, J.: The MesoNH Atmospheric Simulation System. Part I: adiabatic formulation and control simulations, Ann. Geophys., 16, 90–109, 1998. a
Lambaerts, J., Lapeyre, G., Plougonven, R., and Klein, P.: Atmospheric response to sea surface temperature mesoscale structures, J. Geophys. Res., 118, 9611–9621, 2013. a
Large, W. G.: Surface Fluxes for Practitioners of Global Ocean Data Assimilation, in: Ocean Weather Forecasting. An Integrated View of Oceanography, edited by: Chassignet, E. P. and Verron, J., chap. 9, Springer, 229–270, 2006. a
Large, W. G. and Yeager, S. G.: The global climatology of an interannually varying airsea flux data set, Clim. Dynam., 33, 341–364, https://doi.org/10.1007/s0038200804413, 2009. a
Lemarié, F. and Samson, G.: A simplified atmospheric boundary layer model for an improved representation of airsea interactions in eddying oceanic models: implementation and first evaluation in NEMO (v4.0)), Zenodo, https://doi.org/10.5281/zenodo.3904518, 2020. a, b
Lemarié, F., Kurian, J., Shchepetkin, A. F., Molemaker, M. J., Colas, F., and McWilliams, J. C.: Are there inescapable issues prohibiting the use of terrainfollowing coordinates in climate models?, Ocean Modell., 42, 57–79, https://doi.org/10.1016/j.ocemod.2011.11.007, 2012. a
LeMone, M. A., Angevine, W. M., Bretherton, C. S., Chen, F., Dudhia, J., Fedorovich, E., Katsaros, K. B., Lenschow, D. H., Mahrt, L., Patton, E. G., Sun, J., Tjernström, M., and Weil, J.: 100 Years of Progress in Boundary Layer Meteorology, Meteorol. Monogr., 59, 9.1–9.85, https://doi.org/10.1175/AMSMONOGRAPHSD180013.1, 2019. a, b
Lilly, D.: The representation of smallscale turbulence in numerical simulation experiments, in: Proc. IBM Sci. Comput. Symp. on Environmental Sci., 14–16 November 1966, Thomas J. Watson Res. Center, Yorktown Heights, N. Y., IBM Form 320–1951, 195–210, 1967. a
Lindzen, R. S. and Nigam, S.: On the role of sea surface temperature gradients in forcing lowlevel winds and convergence in the tropics, J. Atmos. Sci., 44, 2418–2436, 1987. a
Madec, G.: NEMO ocean engine, in: Note du Pole de modélisation No. 27, Institut PierreSimon Laplace (IPSL), France, 2012. a, b, c
Maisonnave, E. and Masson, S.: Ocean/seaice macro task parallelism in NEMO, in: Technical report, TR/GMGC/15/54, available at: https://www.cerfacs.fr/~maisonna/Reports/opa_sas_tr.pdf (last access: 20 January 2021), 2015. a
Maisonnave, E. and Masson, S.: NEMO 4.0 performance: how to identify and reduce unnecessary communications, in: Technical report, TR/CMGC/19/19, available at: https://cerfacs.fr/wpcontent/uploads/2019/01/GLOBCTR_MaisonnaveNemo2019.pdf (last access: 20 January 2021), 2019. a
Marchesiello, P., Capet, X., Menkes, C., and Kennan, S.: Submesoscale dynamics in tropical instability waves, Ocean Modell., 39, 31–46, https://doi.org/10.1016/j.ocemod.2011.04.011, 2011. a
McWilliams, J. C., Gula, J., and Molemaker, M. J.: The Gulf Stream North Wall: Ageostrophic Circulation and Frontogenesis, J. Phys. Oceanogr., 49, 893–916, https://doi.org/10.1175/JPOD180203.1, 2019. a
Metzger, E. J., Smedstad, O. M., Thoppil, P. G., Hurlburt, H. E., Cummings, J. A., Wallcraft, A. J., Zamudio, L., Franklin, D. S., Posey, P. G., Phelps, M. W., Hogan, P. J., Bub, F. L., and DeHaan, C. J.: US Navy Operational Global Ocean and Arctic Ice Prediction Systems, Oceanogr., 27, 32–43, https://doi.org/10.5670/oceanog.2014.66, 2014. a
Meurdesoif, Y., Caubel, A., Lacroix, R., Dérouillat, J., and Nguyen, M.: XIOS Tutorial, available at: http://forge.ipsl.jussieu.fr/ioserver/rawattachment/wiki/WikiStart/XIOStutorial.pdf (last access: 20 January 2021), 2016. a
Minobe, S., KuwanoYoshida, A., Komori, N., Xie, S.P., and Small, R. J.: Influence of the Gulf Stream on the troposphere, Nature, 452, 206–209, 2008. a
Monin, A. S. and Obukhov, A. M.: Basic laws of turbulent mixing in the surface layer of the atmosphere, Trudy Akademii Nauk SSSR Geofizicheskogo Instituta, 24, 163–187, 1954. a
Mulholland, D. P., Laloyaux, P., Haines, K., and Balmaseda, M. A.: Origin and Impact of Initialization Shocks in Coupled AtmosphereOcean Forecasts, Mon. Weather Rev., 143, 4631–4644, https://doi.org/10.1175/MWRD150076.1, 2015. a, b
Oerder, V., Colas, F., Echevin, V., Masson, S., Hourdin, C., Jullien, S., Madec, G., and Lemarié, F.: Mesoscale SST–wind stress coupling in the Peru–Chile current system: Which mechanisms drive its seasonal variability?, Clim. Dynam., 47, 2309–2330, 2016. a, b
O'Neill, L. W., Esbensen, S. K., Thum, N., Samelson, R. M., and Chelton, D. B.: Dynamical Analysis of the Boundary Layer and Surface Wind Responses to Mesoscale SST Perturbations, J. Climate, 23, 559–581, https://doi.org/10.1175/2009JCLI2662.1, 2010. a, b
Razavi, S., Tolson, B. A., and Burn, D. H.: Review of surrogate modeling in water resources, Water Resour. Res., 48, W07401, https://doi.org/10.1029/2011WR011527, 2012. a
Redelsperger, J. L., Mahé, F., and Carlotti, P.: A Simple And General Subgrid Model Suitable Both For Surface Layer And FreeStream Turbulence, Bound.Lay. Meteorol., 101, 375–408, 2001. a, b, c
Renault, L., Molemaker, M. J., Gula, J., Masson, S., and McWilliams, J. C.: Control and Stabilization of the Gulf Stream by Oceanic Current Interaction with the Atmosphere, J. Phys. Oceanogr., 46, 3439–3453, https://doi.org/10.1175/JPOD160115.1, 2016a. a
Renault, L., Molemaker, M. J., McWilliams, J. C., Shchepetkin, A. F., Lemarié, F., Chelton, D., Illig, S., and Hall, A.: Modulation of Wind Work by Oceanic Current Interaction with the Atmosphere, J. Phys. Oceanogr., 46, 1685–1704, 2016b. a, b
Renault, L., Lemarié, F., and Arsouze, T.: On the implementation and consequences of the oceanic currents feedback in ocean–atmosphere coupled models, Ocean Modell., 141, 101 423, https://doi.org/10.1016/j.ocemod.2019.101423, 2019a. a, b, c
Renault, L., Masson, S., Oerder, V., Jullien, S., and Colas, F.: Disentangling the Mesoscale OceanAtmosphere Interactions, J. Geophys. Res., 124, 2164–2178, https://doi.org/10.1029/2018JC014628, 2019b. a, b, c
Rodier, Q., Masson, V., Couvreux, F., and Paci, A.: Evaluation of a Buoyancy and Shear Based Mixing Length for a Turbulence Scheme, Front. Earth Sci., 5, 65, https://doi.org/10.3389/feart.2017.00065, 2017. a, b, c, d, e, f, g, h
Rotta, J.: Statistische theorie nichthomogener turbulenz, Z. Physik, 129, 547–572, 1951. a
Rousset, C., Vancoppenolle, M., Madec, G., Fichefet, T., Flavoni, S., Barthélemy, A., Benshila, R., Chanut, J., Levy, C., Masson, S., and Vivier, F.: The LouvainLaNeuve sea ice model LIM3.6: global and regional capabilities, Geosci. Model Dev., 8, 2991–3005, https://doi.org/10.5194/gmd829912015, 2015. a, b
Schneider, N. and Qiu, B.: The Atmospheric Response to Weak Sea Surface Temperature Fronts, J. Atmos. Sci., 72, 3356–3377, https://doi.org/10.1175/JASD140212.1, 2015. a, b, c
Seager, R., Blumenthal, M. B., and Kushnir, Y.: An advective atmospheric mixed layer model for ocean modeling purposes: global simulation of surface heat fluxes, J. Climate, 8, 1952–1964, 1995. a
Small, R. J., deSzoeke, S. P., Xie, S. P., O'Neill, L., Seo, H., Song, Q., Cornillon, P., Spall, M., and Minobe, S.: Airsea interaction over ocean fronts and eddies, Dynam. Atmos. Oceans, 45, 274–319, https://doi.org/10.1016/j.dynatmoce.2008.01.001, 2008. a
Soares, P. M. M., Miranda, P. M. A., Siebesma, A. P., and Teixeira, J.: An eddydiffusivity/massflux parametrization for dry and shallow cumulus convection, Q. J. Roy. Meteor. Soc., 130, 3365–3383, https://doi.org/10.1256/qj.03.223, 2004. a
Spall, M.: Midlatitude Wind Stress–Sea Surface Temperature Coupling in the Vicinity of Oceanic Fronts, J. Climate, 20, 3785–3801, https://doi.org/10.1175/JCLI4234.1, 2007. a, b
Takano, K., Mintz, Y., and Han, J.Y.: Numerical simulation of the world ocean circulation, Second Conf. on Numerical Weather Prediction, Monterey, CA, Amer. Meteor. Soc., 121–129, 1973. a
Troen, I. B. and Mahrt, L.: A simple model of the atmospheric boundary layer; sensitivity to surface evaporation, Bound.Lay. Meteorol., 37, 129–148, https://doi.org/10.1007/BF00122760, 1986. a
von Schuckmann, K., Traon, P.Y. L., Smith, N., Pascual, A., Brasseur, P., Fennel, K., Djavidnia, S., Aaboe, S., Fanjul, E. A., Autret, E., Axell, L., Aznar, R., Benincasa, M., Bentamy, A., Boberg, F., BourdalléBadie, R., Nardelli, B. B., Brando, V. E., Bricaud, C., Breivik, L.A., Brewin, R. J., Capet, A., Ceschin, A., Ciliberti, S., Cossarini, G., de Alfonso, M., de Pascual Collar, A., de Kloe, J., Deshayes, J., Desportes, C., Drévillon, M., Drillet, Y., Droghei, R., Dubois, C., Embury, O., Etienne, H., Fratianni, C., Lafuente, J. G., Sotillo, M. G., Garric, G., Gasparin, F., Gerin, R., Good, S., Gourrion, J., Grégoire, M., Greiner, E., Guinehut, S., Gutknecht, E., Hernandez, F., Hernandez, O., Høyer, J., Jackson, L., Jandt, S., Josey, S., Juza, M., Kennedy, J., Kokkini, Z., Korres, G., Kõuts, M., Lagemaa, P., Lavergne, T., le Cann, B., Legeais, J.F., LemieuxDudon, B., Levier, B., Lien, V., Maljutenko, I., Manzano, F., Marcos, M., Marinova, V., Masina, S., Mauri, E., Mayer, M., Melet, A., Mélin, F., Meyssignac, B., Monier, M., Müller, M., Mulet, S., Naranjo, C., Notarstefano, G., Paulmier, A., Gomez, B. P., Gonzalez, I. P., Peneva, E., Perruche, C., Peterson, K. A., Pinardi, N., Pisano, A., Pardo, S., Poulain, P.M., Raj, R. P., Raudsepp, U., Ravdas, M., Reid, R., Rio, M.H., Salon, S., Samuelsen, A., Sammartino, M., Sammartino, S., Sandø, A. B., Santoleri, R., Sathyendranath, S., She, J., Simoncelli, S., Solidoro, C., Stoffelen, A., Storto, A., Szerkely, T., Tamm, S., Tietsche, S., Tinker, J., Tintore, J., Trindade, A., van Zanten, D., Vandenbulcke, L., Verhoef, A., Verbrugge, N., Viktorsson, L., von Schuckmann, K., Wakelin, S. L., Zacharioudaki, A., and Zuo, H.: Copernicus Marine Service Ocean State Report, J. Oper. Oceanogr., 11, S1–S142, 2018. a
Wallace, J. M., Mitchell, T. P., and Deser, C.: The Influence of SeaSurface Temperature on Surface Wind in the Eastern Equatorial Pacific: Seasonal and Interannual Variability, J. Climate, 2, 1492–1499, https://doi.org/10.1175/15200442(1989)002<1492:TIOSST>2.0.CO;2, 1989. a
Wilson, J. M. and Venayagamoorthy, S. K.: A ShearBased Parameterization of Turbulent Mixing in the Stable Atmospheric Boundary Layer, J. Atmos. Sci., 72, 1713–1726, https://doi.org/10.1175/JASD140241.1, 2015. a
This remark is supported by the conclusions of the CLIVAR Working Group on Model Development following the Kiel meeting in April 2014: http://www.clivar.org/sites/default/files/documents/exchanges65_0.pdf (last access: 20 January 2021)
However, we did not find significant differences in numerical solutions when using the following initial conditions: ${\mathbf{u}}_{\mathrm{h}}(z,t=\mathrm{0})={\mathbf{u}}_{\mathrm{G}},\phantom{\rule{1em}{0ex}}e(z,t=\mathrm{0})={e}_{\mathrm{min}}.$
 Abstract
 Introduction
 Model equations
 Numerical discretization and implementation within NEMO
 Atmosphereonly numerical experiments
 Coupled numerical experiments
 Conclusions
 Appendix A: Surface boundary conditions for TKE and mixing lengths
 Appendix B: Preprocessing of atmospheric data from IFS
 Appendix C: Discrete algorithms to compute l_{up} and l_{dwn}
 Appendix D: Code performance and profiling
 Code and data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Model equations
 Numerical discretization and implementation within NEMO
 Atmosphereonly numerical experiments
 Coupled numerical experiments
 Conclusions
 Appendix A: Surface boundary conditions for TKE and mixing lengths
 Appendix B: Preprocessing of atmospheric data from IFS
 Appendix C: Discrete algorithms to compute l_{up} and l_{dwn}
 Appendix D: Code performance and profiling
 Code and data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References