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)

. A simpliﬁed model of the atmospheric boundary layer (ABL) of intermediate complexity between a bulk parameterization and a three-dimensional atmospheric model is developed and integrated to the Nucleus for European Mod-elling of the Ocean (NEMO) general circulation model. An objective in the derivation of such a simpliﬁed model, called ABL1d, is to reach an apt representation in ocean-only 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 large-scale atmospheric data available from re-analysis or real-time 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 NEMO-ABL1d coupling infrastructure and its computational efﬁciency. The resulting simpliﬁed model is then tested for several boundary-layer regimes relevant to either ocean–atmosphere or sea-ice–atmosphere coupling. The coupled system is also tested with a realistic 0 . 25 ◦ resolution global conﬁguration. The numerical results are evaluated using standard metrics from the literature to quantify the wind–sea-surface-temperature (a.k.a. thermal feed-back effect), wind–current (a.k.a. current feedback effect), and ABL–sea-ice 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 ﬂexibility of ocean-only modeling.

Instead of directly using φ atm LS (x, y, z = 10 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 large-scale atmospheric data associated both with the fine resolution in the oceanic surface fields and with the two-way air-sea coupling. Somehow we aim at finding a methodology to get a cheap estimate φ atm HR (x, y, z = 10 m, t) of the solution that would have been obtained using a coupling of M atm and the oceanic model at high resolution. To do so we could imagine several approaches: (i) estimate 110 ∂M atm ∂φ oce LS (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 M 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 data-driven models emulating the original model responses while the third one enters the class 115 of low-fidelity 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 equations 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 three-dimensionalized via advective terms to couple the vertical columns with each other. The idea here is to translate this idea to the MABL context. Our approach 120 has some similarities with the diagnostic meteorological model CALMET (Scire et al., 2000) based on objective analysis.
Using a wind field given by observations and/or a large scale model as an initial guess as well as fine resolution land-use and topography data, CALMET generates a new wind field which accounts for small scale processes like slope flows or terrain blocking effects, among other things. In our case the model will be prognostic and will be specialized for overwater conditions.
In the rest of this section we describe the continuous formulation of our simplified MABL model which will be referred to as 125 ABL1d.

Formulation of a single column approach
The formulation of the ABL1d model is derived under the following assumptions: (i) horizontal homogeneity (i.e. ∂ x · = ∂ y · = 0), (ii) the atmosphere in the computational domain is transparent (i.e. ∂ z I = 0 with I 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. 130 Large and Yeager, 2009) (iii) vertical advection is neglected. Such assumptions prevent the model to prognostically account for the SST-induced 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 φ atm (z = 10 m, t) hence air-sea fluxes. This mechanism is 135 expected to explain most of the eddy-scale wind-SST and wind-currents interactions and is key to properly downscale largescale atmospheric data produced by a coarse resolution GCM to the oceanic resolution.
At a given location in space, the ABL1d model for the Reynolds-averaged profiles of horizontal velocities u h (z, t), potential temperature θ(z, t), and specific humidity q(z, t), provided a suitable initial condition, is (1) 140 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 sea-ice covered areas are of interest. In (1), k = (0, 0, 1) t is a vertical unit vector, f the Coriolis parameter, K m and K s are the eddy diffusivity respectively for momentum and scalars, the subscript LS is used to characterize large-scale quantities known a priori, λ s (z, t) is the inverse of a relaxation timescale, and R LS denotes a large-scale forcing for the momentum equation. R LS can either represent a forcing by geostrophic winds u G (i.e. R LS = f k×u G ) or equivalently by a 145 horizontal pressure gradient (i.e. R LS = 1 ρa ∇ h p LS ) combined with a standard Newtonian relaxation (i.e. R LS = λ m (u LS − u 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 to drift very far away from the large-scale values used to "guide" the model. By itself a relaxation term does not represent directly any real physical process, but the rationale is that it accounts for the influence of large-scale threedimensional circulation processes not explicitly represented in a simple 1d model. Note that this methodology is currently used 150 to evaluate GCM parameterizations in 1D column model. Once the turbulent mixing and the Coriolis term have been computed to provide a provisional prediction φ n+1, at time n + 1 for any ABL1d prognostic variable φ, the relaxation term provides a weighting between this prediction and the large-scale quantities with ∆t the increment of the temporal discretization. Above the boundary layer, the ABL1d formulation is unable to properly 155 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 fine resolution oceanic state and thus the relaxation toward φ LS should be small. The exact form of the λ s and λ m coefficients is discussed later in Sec. 2.4. Note that because of the relaxation term, three-dimensional atmospheric data for u LS , θ LS , q LS , and possibly 1 ρa ∇ h p LS sampled between z sfc and z top must be provided to the oceanic model instead of the two-dimensional data (usually at 10 m) necessary for an ASL forcing strategy.

160
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 consistency, it is preferable to use a bulk formulation as close as possible to the one used to compute the threedimensional large-scale atmospheric data φ atm LS . Because in the present study the plan is to use a large-scale forcing from ECMWF reanalysis products, we use the IFS 2 bulk formulation such as implemented in the AeroBulk 3 package (Brodeau et al., 2017) to compute C D , C θ , and C q in realistic simulations (see Sec. 5). Note that for an ASL forcing strategy u h (z sfc ) and φ(z sfc ) in (3) would be respectively equal to u LS (z = 10 m) and φ LS (z = 10 m) while in the ABL coupling strategy 170 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 gradient 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.

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 so-called CBR-1d 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 180 reasons: computational efficiency, consistency with the NEMO coding rules, use of a geopotential vertical coordinate, and flexibility to add elements specific of the marine atmospheric boundary layer.
CBR-1d is a one-equation turbulence closure model based on a prognostic turbulent kinetic energy (TKE) and a diagnostic computation of appropriate length scales. The prognostic equation for the TKE e = 1 2 ( u u + v v + w w ) (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 ε a dissipation term. In order to express the evolution of e in terms of Reynolds averaged atmospheric variables, the following closure assumptions for the first order turbulent fluxes are made: where l ε is a dissipative length scale, c a constant, and N 2 is the moist Brunt-Väisälä frequency computed as 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 (K m , K s , K e ) = (C m , C s φ z , C e )l m √ e with (C m , C s , C 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 φ max z = 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 m = l ε = 2e/N 2 ), values of the maximum Prandtl number Pr t = C m /(C s φ z ) are given in Tab. 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 205 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 (Tab. 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 210 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) and Cheng et al. (2002). Prt = Cm/(Csφz) is the turbulent Prandtl number.
Dirichlet boundary conditions for TKE are applied at the top z = z top and at the bottom z = z sfc e(z = z top ) = e min = 10 −6 m 2 s −2 with u and w the friction and convective velocities given by the bulk formulation. The value for e min has been chosen 215 empirically as well as background values K min m = 10 −4 m 2 s −1 and K min s = 10 −5 m 2 s −1 for eddy diffusivities. The value of e sfc arise from the similarity theory in the surface layer under the assumption of an equilibrium between TKE shear production, buoyancy production/destruction, and dissipation (Redelsperger et al., 2001). The minimum value for l m is simply set as There are multiple options to compute the mixing lengths l m and l ε (this point will be discussed later in Sec.

3.2) but all options have identical boundary conditions
Again the value of L sfc results from the similarity theory in the neutrally stratified surface layer (Sec. 4.1 in Redelsperger et al., 2001). In (8), κ is the von Karman constant and z 0 a roughness length computed within the bulk algorithm.
Our current implementation of boundary layer subgrid processes is an eddy-diffuvisity approach which does not include any explicit representation of boundary-layer convective structures. This could be done via a mass-flux representation (e.g. Hourdin 225 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.

Processing of large-scale forcing and Newtonian relaxation
As mentioned earlier, the ABL1d model (1) requires three-dimensional (x, y, z) large-scale atmospheric variables φ atm LS while existing uncoupled oceanic forcing strategies require only two-dimensional (x, y) atmospheric variables. This is a difficulty 230 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 large-scale 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 ERA-Interim, ERA5 and operational IFS datasets and are described in App. A.  We considered the following general form for λ s (z) (resp. λ m (z)), with h bl the boundary layer height whose value is diagnosed using an integral Richardson number criteria (Sec. 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 = β min h bl and 250 z = β max h bl . We easily find The value of h bl is bounded beforehand to guarantee that at least 3 grid points are such that z ≤ β min h bl and z ≥ β max h bl . A typical profile of the λ s (z) is shown in Fig. 1a.
When the model is forced by the large-scale pressure gradient (or the geostrophic winds), the parameter λ m (z) should be 255 theoretically zero at high and mid latitudes. However for the equatorial region, a Newtonian relaxation toward the large scale 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 equal to zero for large values of |f | and increases to one when approaching the equator. The following form satisfies those constraints (see also Fig. 1b)

Coriolis term treatment
Since in our implementation the horizontal velocity components are collocated, the discretization of the Coriolis term is 270 straightforward and is energetically neutral. In the event the ABL1d is integrated with a time-step much larger than the oceanic time-step, a specific care must be given to the stability of the Coriolis term time-stepping. For a given grid cell with index (i, j), a semi-implicit scheme with weighting parameter γ reads The exponent is used here to emphasize that u n+1, are temporary values at time n + 1 before vertical diffusion 275 and Newtonian relaxation are applied. (11) can be written in a more compact way as The associated amplification factor modulus is |A cor | = 1 + (1 − γ) 2 (f ∆t) 2 1 + γ 2 (f ∆t) 2 meaning that unconditional stability is obtained as long as γ ≥ 1/2 . For the numerical results obtained below in Sec. 4 and 5 we used γ = 0.55 which is deliberately 280 slightly dissipative.

TKE positivity-preservation and mixing lengths computation
In Sec. 2.3 we have presented the continuous formulation of the TKE-based turbulence closure of the ABL1d model. The TKE equation is discretized using a backward Euler scheme in time with a linearization of the dissipation term c ε l ε e 3/2 which is discretized as c ε l ε √ e n e n+1 . However, such discretization is not unconditionally positivity-preserving for TKE which could 285 give rise to unphysical solutions (e.g. Burchard, 2002b). Ignoring the diffusion term, the TKE prognostic equation (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 (12) it can be shown that for an initial condition e(0) ≥ 0 and S(u h , N 2 ) ≥ 0, the solution e(t) keeps the same sign as e(0) whatever the sign of the damping coefficient D(e, t). Assuming 290 that S(u h , N 2 ) and D(e, t) are positive, a backward Euler discretization of the damping term in (12) would lead to which preserves positivity since for e n ≥ 0 we obtain e n+1 ≥ 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 the idea is to move the buoyancy term from S to D after dividing it by e n , such that S(u h , N 2 ) = K m ∂ z u h 2 is now strictly positive and D(e, t) = c ε l ε √ e n + K s N 2 e n . Such 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).
An other challenging task when implementing a TKE scheme is the discretization of the mixing lengths. As mentioned 300 earlier, 4 different discretizations of l m (resp. l ε ) have been coded. All discretizations consider the boundary conditions given in (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 method described below, the dissipative and mixing length scale l m and l ε are computed as where a ≈ − 3 2 for CBR00 and a ≈ − √ 3 2 for CCH02. 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 Sec. 4.1 but for more realistic cases results are weakly sensitive and equivalent to the ones obtained with the simpler weighting l m = l up l dwn .

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 which 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.

315
It is worth noting that for N 2 = cste, (15) gives respectively l 2 up N 2 2 = e(z) and l 2 dwn N 2 2 = e(z) which is equivalent to the Deardorff (1980) length scale. An objective is to satisfy this property also at a discrete level. Considering a simple trapezoidal rule to approximate the integral in (15) over each grid cells, the procedure for the computation of l up (z k+1/2 ) is given in Algorithm 1. In the case N 2 (z p+1/2 ) = N 2 (z p−1/2 ) = N 2 cst (∀p), Algorithm 1 gives the following sequence Algorithm 1 Procedure to compute the Bougeault and Lacarrère (1989) length scale l up (z k+1/2 ).
As soon as FC(z p+1/2 ) changes sign we stop the procedure because l up such that −e(z k+1/2 )+N 2 cst l 2 up = 0, which corresponds to the Deardorff (1980) length scale, has been found. In the remainder we will note l BL89 the mixing length corresponding to 325 the Bougeault and Lacarrère (1989) algorithm.

Adaptation of NEMO's length scale
The standard NEMO algorithm (Sec. 10.1.3 in Madec, 2012) is a simpler and more efficient version of the Bougeault and Lacarrère (1989) which is much easier to discretize. As a first step the Deardorff (1980) , l min with N 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 will be simply referred to as l D80 . Note that the Taylor expansion of the integral in the Bougeault and Lacarrère (1989) mixing length definition is which shows that the l D80 mixing length is an approximation of l BL89 which is obtained by retaining only the leading order 340 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). This modification turns out to improve results for stably stratified boundary layers typical of areas covered by ice. They propose to add a shear related term to (15) such that the definition of l up and l dwn becomes where c 0 is a parameter whose value should be smaller than C m /c ε . The value of c 0 will be chosen based on numerical experiments presented in Sec. 4. At a discrete level, the FC function in Algorithm 1 is replaced by In the following this mixing length will be referred to as l R17 .

A local buoyancy and shear-based 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 (resp. l dwn ) is small compared to the spatial variations of N 2 , e, and ∂ z u h , we end up with the following second-order equation for l up whose unique positive solution is .
We easily find that l D80 = l D80 for ∂ z u h = 0, and l D80 = e(z) c 0 ∂ z u h for N 2 = 0 which is consistent with the shear based length scale of Wilson and Venayagamoorthy (2015). Once l D80 has been computed we apply the limitations (16) and (17) as 360 in the NEMO algorithm.
The performance of those four length scales for various physical flows is discussed in Sec. 4.

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. Lemarié et al., 2015;Beljaars et al., 2017;Renault et al., 2019a) and the 365 ABL1d must have the ability to handle grid cells partially covered by sea-ice. For the coupling strategy, a so-called implicit flux where u oce and φ 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.
A particular care has also been given to the compatibility between the ABL1d model and SI 3 (Sea Ice model Integrated Initiative) the sea-ice component of NEMO. SI 3 is a multi-category model whose state variables relevant for our study are  (19) and (20) 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 the different steps 1. Compute surface fluxes over ice and ocean and integrate the ABL1d model for given values F n oce and a n l .
2. Compute the dynamics of sea-ice 3. Update F n oce and a n l in F oce and a l because of step 2. 5. Compute the thermo-dynamics of sea-ice

Additional details about the implementation
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)  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 proces-400 sors and computational grids communicating via the OASIS3 − MCT coupling library (Craig et al., 2017).
An other capability implemented within the NEMO modelling framework is the possibility to interpolate forcing fields on-thefly. This is particularly useful for the ABL coupling strategy since three-dimensional atmospheric data must be interpolated on the ABL1d computational grid. As the current implementation of the on-the-fly interpolation only works in the horizontal, the vertical interpolation of large-scale atmospheric data on the ABL1d vertical grid is done offline. Nevertheless it means that the 405 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 IO capabilities provided by the XIOS library (Xml-IO-Server, 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 Sec. 5.3, the main source of computational overhead associated with the ABL coupling strategy is due to the time spent waiting for input 410 files to be read.

Atmosphere-only numerical experiments
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 single column experiments. Each of those experiments have been run with a companion large-eddy simulation (LES) model whose solution is considered as the reference. In the following we consider a neutrally

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  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 2 ε on the stratification is imposed in the code. In this case the procedure to compute those mixing lengths as described in Sec. 3.2.1 and 3.2.2 will provide identical results, namely l up = z top − z and l dwn = z − z sfc (i.e. the distance from the top and from the bottom of the computational domain). We test

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/ 440 gabls1d_desc.pdf. This experiment is particularly interesting as significant differences generally exist between solutions obtained from LES and single-column 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 large scale geostrophic wind is imposed as well as a cooling of the surface temperature θ s (t) given by θ s (t) = 263.5 − 0.25(t/3600 s). The parameter values for this test are reported in Tab. 2 and the initial conditions are u h (z, t = 0) = u G , and (not shown). The mismatch in terms of TKE is partially explained by the difference in boundary condition since with CBR00 constants we have e sfc = 4.628 u 2 while with CCH02 constants we get e sfc = 3.065 u 2 from Eq. (7). Note that the proper 455 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.
So far the CCH02 set of parameters turn out to provide 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 SST front is centered at x = 0 km.

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-485 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 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   the expected behaviour and we will focus the discussion on other options. In terms of zonal 10m wind the buoyancy based 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 function of buoyancy and vertical shear (as is the case for l D80 ) the simulated winds are weaker because the boundary layer is thinner. This leads to improved 500 results in the stably stratified case shown earlier but in the present case more representative of realistic configuration in the MABL it leads to a too weak mixing. However, compared to the l R17 mixing length the l D80 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. For the last numerical experiments we will study, the l R17 mixing length will be discarded from the comparison.

505
Although relevant for the present study this 2D x-z experiment 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.

A single-column version
An alternative to the x-z setup would be to formulate the testcase as a Lagrangian advection of an air column over an SST front by prescribing a temporal evolution of sea surface temperature θ s (t) as 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 testcase 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 formulation 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 520 the temperature variable toward the initial condition with a fast restoring time scale λ s = 1 6[h] above the PBL and a slower one λ s = 1 48 [h] in the PBL. This is meant to check that the relaxation does not completely overwrites the physics of the coupling we aim at representing with the ABL1d model. Results from those sensitivity experiments are shown in Fig. 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, 525 adding a relaxation toward large-scale data which did not see the SST front does not deteriorate the realism of the solutions, as can be seen from Fig. 8 (rightmost panels) and 9 (gray lines).
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) mixing length formulation l D80 or l D80 .

530
In addition to idealized atmosphere-only simulations, we now show an illustration of the results obtained when the ABL1d model is coupled with NEMO to represent explicitly interactions between the upper ocean, sea-ice and the MABL. To do so, we performed a 5-years global simulation using the ORCA025 configuration. Details and illustrations are given hereafter.

Coupled NEMO-ABL1d configuration
We use here a global ORCA025 configuration at a 0.25 • horizontal resolution (Barnier et al., 2006) with 75 vertical z-levels 535 forced by ECWMF ERA-Interim 6-hours analysis (Dee et al., 2011). This configuration is identical to the one described in  Table 3. Namelist parameters in the NEMO(v4.0) to set in the namelist section namsbc_abl before running a simulation coupled with ABL1d.

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 550 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/sea-ice couplings (e.g. Bryan et al., 2010;Renault et al., 2019b). To quantify the surface wind response to SST, we show in Fig. 10 global maps of the temporal correlation between the high-pass filtered 10 m wind speed from the first vertical level in the ABL1d model and the SST. Same correlation is shown in Bryan et al. (2010) from satellite observations (their Fig. 1d) and from coupled 555 numerical experiments between a 0.1 • ocean and a 0.25 • atmospheric models (their Fig. 1c). Consistent with observations and fully-coupled models, the correlation obtained from the coupled NEMO-ABL1d 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 • W 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 (1/12 • resolution). This coupling sensitivity to the oceanic resolution will be addressed in a future study.
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 565 to the ABL and s w quantifies the partial re-energization of the ocean by the wind response to the wind/current coupling. This re-energization is absent in the ASL forcing strategy which results in an excessive dampening of the oceanic eddy mesoscale activity. In practice, s τ (resp. s w ) corresponds to the slope of the linear relationship between high-pass filtered surface current vorticity and surface wind-stress (resp. wind) curl. Global maps of s τ and s w computed from our coupled NEMO-ABL1d global simulation are shown in Fig. 11. Large negative values of s τ indicate an efficient dampening of the eddy mesoscale

575
The last illustration of our implementation presented in this section is the coupling of ABL1d with sea-ice. As described in Sec. 4.2, sea-ice generally induces a shallow stably stratified boundary layer due to the near-surface air cooling. This increased s w s ⌧ N s m 3 Figure 11. Global maps of the coupling coefficient between the surface current vorticity and the wind curl (sw, top) and between the surface current vorticity and the wind-stress curl (sτ , bottom) estimated from a 0.25 • resolution coupled NEMO-ABL1d global simulation. The fields are first temporally averaged using a 29-day running mean and spatially high-pass filtered. ice-covered regions. This coupling between the PBL and sea-ice have important effects on near-surface wind, temperature and humidity, and consequently on sea-ice concentration evolution which will need to be specifically assessed in future ABL-based studies.

Code performance and profiling
At this point in the manuscript, we anticipate that the interested reader has already wondered several times what could be the 585 computational cost of the proposed methodology. To finalize our description of the implementation of this simplified atmospheric boundary layer model in NEMO, we assess 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 590 tests investigated earlier in the manuscript. 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 "−i4 − r8 − O3 − fp − model precise − fno − alias" options. The I/Os are handled via 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 built-in NEMO code instrumentation dedicated to calculation measurement (e.g. Maisonnave and Masson, 2019). As 595 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 days such that the cost of the initialization step is no longer visible in the averaged cost per time-step. Moreover, 600 for the two simulations, the atmospheric data necessary for the computation of the turbulent components of air-sea fluxes are provided every 6 hours.
We first show in Fig. 13

615
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 meter height atmospheric quantities are prescribed. For this preliminary study, the simplified ABL model takes the form of a single column 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 so-called downward 620 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 meters 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 meter atmospheric values to the ASL parameterization. An important point is that our modeling strategy keeps the computational efficiency and flexibility inherent to ocean only modeling.

625
In this paper the key components of such an approach have been described. This includes the large-scale 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 sea-ice/atmosphere coupling. Results have systematically been evaluated configuration show an excellent behaviour in term of wind-SST two-way coupling. A thorough analysis of the impact of the coupling with ABL1d from a physical viewpoint will be presented in subsequent work (e.g. Brivoal et al., 2020).
Several ways to improve the methodology presented here are currently under investigation. The continuous formulation of the ABL1d model will be completed by adding a mass-flux representation in our turbulent closure scheme and and by 635 integrating the effect of horizontal advection and fine-scale pressure gradient. As mentioned several times in the paper ways to lower the computational overhead due to I/O inputs will be investigated. Moreover, following the work of Couvelard et al. (2019) it is planned to explicitly account for the effect of surface waves in the ABL coupling strategy.
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 model-driven empirical models) to be seamlessly tested and confronted with the ABL 640 coupling strategy.  of the table). The timing is averaged on all processors. The right most column provides a quick description of the task handled by the corresponding section. On top of the timing in seconds the percentage of the total CPU and elapsed associated to each section is reported in parentheses. The computational overhead associated to the ABL coupling strategy can be estimated from the sbc section and the elapsed/CPU time.
Code and data availability. 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.nemo-ocean.eu/, last access: 23 June 2020). The NEMO code modified to include the

A1 Altitude of IFS vertical levels
The ABL1d model is discretized on fixed in time and space geopotential levels while the IFS model uses a pressure-based sigma coordinate. A first step is to recover the altitude associated with each sigma level. The pressure p k+ 1 2 defined at cell interfaces between two successive vertical layers is given by where A k+ 1 2 (Pa) and B k+ 1 2 (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 60 levels grid (L60) and a surface pressure of 1013 hPa are given in Tab. A1.
Once the values of p k+ 1 2 and p s are known, the altitude of cell interfaces can be computed by integrating the hydrostatic equilibrium vertically. In (A1), φ is the geopotential, T v the virtual temperature and R d the specific gas constant for dry air. At a discrete level we get Once the the layer thicknesses e3t ifs k are known, horizontal wind components, potential temperature and specific humidity can be interpolated on the ABL1d vertical levels. Under the constraint that ztop z sfc ψ ifs (z) dz = ztop z sfc ψ(z) dz for any IFS quantity ψ ifs to be interpolated. Wind components are interpolated using a fourth-order compact scheme while tracers are interpolated using a WENO-like PPM scheme (A. Shchepetkin, personal communication) which is monotonic.

A2 Filtering in the presence of boundaries
Because of the IFS numerical formulation and of the post-processing of output data, the solutions sometimes contain high frequency oscillations at the vicinity of the land-sea 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.  Table A1. Altitude z k and layer thickness e3t k of the IFS L60 vertical grid in the first 2000 meters with respect to the parameter values A k and B k of a surface pressure ps = 1013 hPa.

680
A shap (θ x , π) = 0 and that A 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 no-gradient condition at the coast, i.e.
δ (x) i+1/2,j = 0 as soon as tmask i+1,j × tmask i,j = 0 (resp. δ (y, ) i,j+1/2 = 0 as soon as tmask i,j+1 × tmask i,j = 0 ) with tmask the indicator function equal to 1 over water and 0 over land. Let us also consider the following alternative boundary conditions    δ (x) i+1/2,j = −δ are shown for different boundary conditions. In particular it can be seen near the coast that the no gradient boundary condition (panels b and e) leaves some artifical patterns in gradients especially in the Peru-Chile current system while boundary condition (A2) efficiently mitigate this issue. Note that it is particularly essential to make sure that the surface pressure field is sufficiently 690 smooth because gradients of this field are used to compute geostrophic winds which are important for the large-scale forcing of the ABL1d model.
where (·) z denotes a gradient along constant geopotential height. Using the hydrostatic balance we have 1 ρa = −g(∂ z p) −1 which leads to Assuming a generalized vertical coordinate s = s(x, y) the computation of gradients along constant height is not straightfor-700 ward since (∂ x p) z = (∂ x p) s − (∂ z p)(∂ x z) 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.