Development of BFMCOUPLER ( v 1 . 0 ) , the coupling scheme that links the MITgcm and BFM models for ocean biogeochemistry simulations

In this paper, we present a coupling scheme between the Massachusetts Institute of Technology general circulation model (MITgcm) and the Biogeochemical Flux Model (BFM). The MITgcm and BFM are widely used models for geophysical fluid dynamics and for ocean biogeochemistry, respectively, and they benefit from the support of active developers and user communities. The MITgcm is a state-of-the-art general circulation model for simulating the ocean and the atmosphere. This model is fully three dimensional (including the non-hydrostatic term of momentum equations) and 15 includes a finite-volume discretization and a number of additional features enabling simulations from global (O(10)m) to local scales (O(100)m). The BFM is a complex biogeochemical model that simulates the cycling of a number of constituents and nutrients within marine ecosystems. The coupler presented in this paper links the two models through an efficient scheme that manages communication and memory sharing between the models. We also test specific model options to balance the numerical accuracy against the computational performance. The coupling scheme allows us to solve several 20 processes that are not considered by each of the models alone, including light attenuation parameterizations along the water column, phytoplankton and detritus sinking, external inputs, and surface and bottom fluxes. Moreover, this new coupled hydrodynamic-biogeochemical model has been configured and tested against an idealized problem (a cyclonic gyre in a midlatitude closed basin) and a realistic case study (central part of the Mediterranean Sea in 2006-2012). The numerical results are consistent with the expected theoretical and observed behaviour of both the idealized system and the Mediterranean 25 domain, thus demonstrating the applicability of the new coupled model to a wide range of ocean biogeochemistry problems.

Abstract.In this paper, we present a coupling scheme between the Massachusetts Institute of Technology general circulation model (MITgcm) and the Biogeochemical Flux Model (BFM).The MITgcm and BFM are widely used models for geophysical fluid dynamics and for ocean biogeochemistry, respectively, and they benefit from the support of active developers and user communities.The MITgcm is a state-of-the-art general circulation model for simulating the ocean and the atmosphere.This model is fully 3-D (including the non-hydrostatic term of momentum equations) and is characterized by a finite-volume discretization and a number of additional features enabling simulations from global (O(10 7 ) m) to local scales (O(100) m).The BFM is a biogeochemical model based on plankton functional type formulations, and it simulates the cycling of a number of constituents and nutrients within marine ecosystems.The online coupling presented in this paper is based on an open-source code, and it is characterized by a modular structure.Modularity preserves the potentials of the two models, allowing for a sustainable programming effort to handle future evolutions in the two codes.We also tested specific model options and integration schemes to balance the numerical accuracy against the computational performance.The coupling scheme allows us to solve several processes that are not considered by each of the models alone, including light attenuation parameterizations along the water column, phytoplankton and detritus sinking, external inputs, and surface and bottom fluxes.
Moreover, this new coupled hydrodynamic-biogeochemical model has been configured and tested against an idealized problem (a cyclonic gyre in a mid-latitude closed basin) and a realistic case study (central part of the Mediterranean Sea in [2006][2007][2008][2009][2010][2011][2012].The numerical results consistently reproduce the interplay of hydrodynamics and biogeochemistry in both the idealized case and Mediterranean Sea experiments.The former reproduces correctly the alternation of surface bloom and deep chlorophyll maximum dynamics driven by the seasonal cycle of winter vertical mixing and summer stratification; the latter simulates the main basin-wide and mesoscale spatial features of the physical and biochemical variables in the Mediterranean, thus demonstrating the applicability of the new coupled model to a wide range of ocean biogeochemistry problems.

Introduction
Coupling different models that have been specifically developed to study only limited aspects of the Earth's systems is becoming increasingly common due to the need to simulate different environmental components -and their interactions -simultaneously (Heavens et al., 2013).As regards numerical oceanography, coupled hydrodynamic-biogeochemical models are widely used to investigate and predict the physical, biogeochemical, and ecological properties of marine ecosystems across a wide range of scales and provide useful tools that support environmental management and policies.
The numerical implementation of a coupling framework between 3-D hydrodynamic models and biogeochemical models is not a trivial task (Bruggeman and Bolding, 2014) because every model focuses on processes that occur on different temporal and spatial scales and uses different numerical parameterizations and schemes.Additionally, these models might be coded in different languages or follow different coding "philosophies" with respect to memory allocation, computational schemes, and code workflow.Furthermore, hydrodynamic and biogeochemical models are often developed by different and highly specialized scientific groups, whereas coupling requires interdisciplinary expertise.
In recent decades, the increasing availability of significant computational resources has allowed substantial improvements in hydrodynamic and biogeochemical models in terms of both temporal and spatial resolution of the simulations, which required new specific programming and coding expertise (i.e.code optimization and parallel programming).In addition, biogeochemical model complexity has increased through the inclusion of new variables and processes (Robson, 2014), and model development has become a cooperative and multidisciplinary task rather than an individual effort.A large number of generic, open-source models are utilized by the scientific community, and they can be customized to match the users' specific applications.A nonexhaustive list of the main state-of-the-art, hydrodynamic community models includes the MITgcm (Adcroft et al., 2016), GOTM (Burchard et al., 2006), ROMS (Haidvogel et al., 2000), and NEMO (Madec, 2014), whereas examples of community biogeochemical models include the BFM (Vichi et al., 2015), ERSEM (Butenschön et al., 2016), PISCES (Aumont et al., 2015), and ERGOM (Neumann, 2000).
Hydrodynamic and biogeochemical models can be coupled by merging their codes into a single larger new code, in which the original parts are intertwined.In this case, biological models are inserted into the workflow of the existing hydrodynamic model code (Burchard et al., 2006;Follows et al., 2006) because, in general, hydrodynamic models have already been developed to solve the partial differential equation of tracers and provide the coding infrastructure to handle the spatial-temporal properties of the simulations (i.e.bathymetry, boundaries, computational domain discretization).Alternatively, a modular approach can be adopted: each component preserves its own peculiarities, the coupling is performed only on localized portions of the code, and there are clear application programming interfaces (APIs).The separation of the two coupled components facilitates the maintenance of each code within its development community, avoids possible large efforts in solving the language differences between models, and eliminates the need to keep models up to date with respect to the parent model.As an example, Bruggeman and Bolding (2014) proposed a set of programming interfaces (FABM) that allows communication between different hydrodynamic and biogeochemical models.
In this paper, we present a coupling scheme between the MITgcm hydrodynamic model and the BFM biogeochemical model for ocean biogeochemical simulations.The two models are widely used, as described in the next sections, and have already been coupled with several other models.For example, the MITgcm has already been coupled to low- (Parekh et al., 2005;Follows et al., 2006) or intermediatecomplexity (Hauck et al., 2013;Cossarini et al., 2015a) biogeochemical models for a few specific applications and to a specific high-complexity model (Dutkiewicz et al., 2009) to explore the theoretical aspects of intraspecific competition in plankton communities.On the other side, the BFM has already been coupled to POM (Polimene et al., 2006), NEMO (Vichi and Masina, 2009;Epicoco et al., 2016), and the offline OGSTM, an upgraded version of OPA (Lazzari et al., 2012).A direct coupling between MITgcm and BFM has not been implemented yet.Thus, we developed a dedicated online modular coupler linking them.The new coupler is open source, and allows us to exploit the high potentiality of the two models, to preserve the sustainability of the programming effort, and to handle the future evolution of the two codes.Further, the online coupling of hydrodynamic and biogeochemical models allows us to drive the biogeochemistry at the same frequency of the hydrodynamic processes, avoiding the use of large files where hydrodynamic variables are saved at high frequency.It also ensures the use of consistent differential operators (advection and diffusion) for hydrodynamic and biogeochemical variables, and would eventually provide a framework to describe possible feedbacks from biogeochemistry to hydrodynamics.
We demonstrate that the new online coupled model provides reliable results when simulating different marine ecosystems by correctly reproducing the interplay between physical, chemical, and biological processes and components.The coupled model also runs with good computational performance and preserves the numerical accuracy of the solution.We consider that the MITgcm-BFM model represents a promising tool for investigating marine biogeochemistry at different spatial and temporal scales.
This paper is organized as follows.After a brief presentation of the two models (Sect.2), we focus on the technical aspects of the coupling algorithm.In the subsequent section (Sect.3), we describe the testing of the new coupled hydrodynamic-biogeochemical model against the idealized case of a cyclonic circulation in a closed basin and against a real case study in the central Mediterranean Sea.The paper closes with a discussion of the key issues of the coupling and future perspectives.A manual of the new code package is detailed in the Appendix.

Formulation of the hydrodynamic-biogeochemical coupling
A coupled hydrodynamic-biogeochemical model is composed of three main elements: a hydrodynamic sub-model, which solves the governing equations for oceanic flows; a tracer transport sub-model, which solves for the transport (advection and diffusion) of biogeochemical variables (commonly called tracers); and a biogeochemical sub-model, which describes the relationships (i.e.biogeochemical reactions) among the biogeochemical variables.Following the common practice in which the biological feedback on transport is negligible, one can assume that changes in biogeochemical properties do not affect the water velocity, density, or other physical properties; therefore, modifying the standard equations that underpin hydrodynamic models is unnecessary.We adopted such an assumption for this numerical coupling framework; however, this coupler was developed, in principle, to also handle biological feedbacks on hydrodynamics.The coupled model solves the set of partial differential equations specified below: Momentum conservation equations, Eqs. ( 1)-( 2), continuity and density equations, Eqs. ( 3)-( 4), and activetracer equations (for potential temperature θ and salinity S), Eqs. ( 5)-( 6), are formulated according to the semicompressible Boussinesq approximation.In the equations, v H = (u, v) is the horizontal component of velocity, w is the vertical velocity, f is the Coriolis parameter, ρ c is a constant reference density, and p is the pressure term.The right-hand side (RHS) terms in Eqs.(1), (2), (5), and (6) correspond to the forcing and dissipation terms, including the diffusion, which acts on the momentum (F H and F V in Eqs.1-2) and on the temperature and salinity (Q θ and Q S in Eqs.5-6).Similarly, Eq. ( 7), which stands for a system of partial differential equations of tracers (C), encompasses the forcing and dissipation terms for biogeochemical tracers, Q C , and the biogeochemical reactions that occur in the sea, R bio .
Equation ( 8) is an equation of state that calculates the modulation of irradiance PAR (photosynthetic active radiation) with depth starting from short-wave surface radiation fields (Q sw ).The total derivative accounts for the partial derivative in time and the advection term, which is related to the flow field, Eq. ( 9).
By adopting a more explicit formulation and commonly used assumptions based on scale analysis (see Crise et al., 1999), Eq. ( 7) can be rewritten as follows: The first three terms on the RHS of Eq. ( 10) represent the advection (first term) and diffusion (second -horizontal -and third -vertical -terms) of biogeochemical tracers, where K H and K V are the horizontal and vertical diffusivities, respectively, which are considered separately because they have different spatial scales.The remaining terms describe the sinking processes that affect biological particles (fourth term) and biogeochemical reactions (fifth term).
Within a coupled model, Eqs. ( 1)-( 6) are solved by the hydrodynamic sub-model, whereas Eq. ( 10) is solved partly by the transport sub-model, which is usually embedded in the hydrodynamic code, and partly by the biogeochemical submodel.The other components, such as Eq. ( 8), the biogeochemical tracers' forcing terms (Q C , e.g.surface and bottom boundary conditions), and the sinking terms, can be handled by either the hydrodynamic or biogeochemical model, according to the specific processes and the features of the codes.
A coupler is defined as the interface that transfers the hydrodynamic information from Eqs. ( 1)-( 6) to Eq. ( 10) and controls the communication between the different terms of Eq. ( 10).In this study, the sub-models coupled are the MITgcm (managing both hydrodynamics and transport) and BFM (for the biogeochemistry) models, which are described in Sects.2.2 and 2.3.The algorithm used to construct the fully coupled system is detailed in Sect.2.4.

Nomenclature and units
Throughout the text, we used the following convention.In equations and text, C refers to the concentration (mass per unit volume) of biogeochemical model state variables, which are referred to as pTracer (passive tracer) in the MITgcm nomenclature.As regards BFM, the chemical components in the subscript are in blackboard style (C: carbon; N: nitrogen; P: phosphorus; S: silica).The pieces of code and the names of the routines and files are in typewriter font.Appendix B reports a list of all symbols and variables used throughout the text.

MITgcm
The MITgcm (Massachusetts Institute of Technology general circulation model; Marshall et al., 1997) is a 3-D, finitevolume, general circulation model used by a broad community of researchers.It can be customized to create different simulation set-ups by modifying its packages and parameters accordingly (Adcroft et al., 2016) and it has already been successfully applied to a wide range of case studies for the world's ocean at various spatial and temporal scales.The code and documentation of the MITgcm are under continuous development.The modular Fortran77 code is open source (copyright ©2016 MITgcm Developers and Contributors), and it can be downloaded from the MITgcm website (http://mitgcm.org/)as a TAR file or using a CVS pserver.The most recent online documentation can be found at http://mitgcm.org/public/r2_manual/latest/.The MITgcm, which is designed to run on high-performance computing (HPC) platforms, can solve fully non-hydrostatic and hydrostatic equations and can handle different freesurface formulations.Subgrid-scale turbulence in both the horizontal and vertical directions can be parameterized by using different types of closure schemes with either constant or variable coefficients (e.g.Gent-McWilliams, Redi, Leith, Smagorinsky, KPP, and GGL90).The MITgcm code is composed of several packages, and, depending on the selected experiment, the compiled packages can be enabled or disabled during the runtime by specifying the flag use PACKAGENAME=.TRUE./.FALSE. in the data.pkginput namelist.The MITgcm's implementation used in this paper was based on the Release 1 -Checkpoint 65 k (April 2015) version of the code.Among the different available customization options, we adopted the fully implicit barotropic time stepping for the free surface, which is unconditionally stable.The vertical diffusion and viscosity terms in the horizontal momentum equations were treated implicitly in time and were solved by using the Euler backward method.The terms that were evaluated explicitly in time were discretized by using the third-order Adams-Bashforth method for the momentum equations and the Euler forward-in-time method for the transport equations.
A native transport sub-model for passive tracers (the Passive TRACERS -PTRACERS -package; according to the MITgcm's jargon, a passive tracer is a generic tracer that has no influence on the hydrodynamics -e.g. by changing the density and/or viscosity) is included in the MITgcm code.This sub-model solves the first three terms on the RHS of Eq. (10) (transport of a generic passive tracer).This transport is calculated by adopting a direct space-time discretization method for the advection-diffusion part of the tracer equations and a non-linear, third-order advection scheme with a Sweby flux limiter (Sweby, 1984) to avoid spurious oscillations in the model output fields.When employing the direct space-time method and the flux-limited schemes, the Euler forward time stepping is adopted rather than Adams-Bashforth.
Because of the different length scales, horizontal and vertical turbulent processes are treated separately and are solved by adopting a selected subset of several available parameterizations: in this study, we chose a mixed Leith-Smagorinsky scheme for the horizontal processes (second term on the RHS of Eq. 10) and the K-profile parameterization (KPP, Large et al., 1994) for the vertical processes (third term on the RHS of Eq. 10).
This code was compiled onto a Linux cluster that was equipped with Intel Xeon Ivy Bridge processors by using both the native GNU compiler (gfortran with openmpi libraries) and the Intel compiler (ifort: Intel Composer XE 2013 SP1) and by adopting the optimization levels -O3 and -O2, respectively.Overall, the model performance increased by approximately 10 % when using the Intel compiler.The results in this paper were obtained using the Intel compiler with the optimization set to -O2.Further details on the custom model installation are available in the MITgcm's online documentation.

BFM
The Biogeochemical Flux Model (BFM) is an open-source, modular Fortran90 numerical model that was designed to describe the dynamics of the major biogeochemical processes that occur in marine ecosystems (Vichi et al., 2015).The standard configuration of the BFM solves the cycles of carbon, phosphorus, nitrogen, silica, and oxygen in the water-dissolved phase and in the plankton, detritus, and benthic compartments.Plankton dynamics are parameterized by considering a number of plankton functional groups, each representing a class of taxa.The BFM's plankton functional groups are subdivided into producers (phytoplankton), consumers (zooplankton), and decomposers (bacteria).These broad functional classifications are further partitioned into functional subgroups to create a planktonic food web (e.g.diatoms, picophytoplankton, microzooplankton).The structure of the plankton functional types is modular and can be adapted to specific needs.In fact, the BFM's code is organized into several modules devoted to several plankton function types: Phytodynamics (for the phytoplankton functional types), Mesozoodynamics and Microzoodynamics (for the zooplankton functional types), and PelBacDynamics (for bacteria).The two modules OxygenReaerationDynamics and PelChemDynamics for the oxygen and carbonate system dynamics, respectively, complete the pelagic system (subroutine PelagicSystemDynamics).The EcologyDynamics interface routine manages the memory allocation of the biogeochemical state variables and derivatives, and the external information that is required to calculate the biological equations: temperature, salinity, presence of ice, wind, position of the cell with respect to the surface or bottom, and atmospheric CO 2 partial pressure.The code and a full description of the model equations and parameterizations are freely available at http://bfm-community.eu.
For this application, we adopted version v2 (Lazzari et al., 2012(Lazzari et al., , 2016;;Teruzzi et al., 2013;Melaku Canu et al., 2015;Cossarini et al., 2015b), which can be downloaded upon request from the BFM consortium website (http:// BFM-consortium.eu)under the GNU GPL license.The current BFM version uses a 0-D data structure for the biogeochemical state variables.The present BFM includes four components (C, N, P, and S); four phytoplankton groups; four zooplankton groups; one group each of bacteria, detritus, labile, and semilabile organic matter; and additional variables, such as dissolved oxygen and alkalinity (Fig. 1).In addition, chlorophyll is solved as a prognostic variable according to the formulation of Geider et al. (1997), and the carbonate system is solved by using the OCMIP formulation (Melaku Canu et al., 2015;Cossarini et al., 2015b).

The coupler
In this coupling scheme, we adopted a modular approach by considering the high complexity of the two models that were employed.The size of the codes according to the SLOC-Count tool (Wheeler, 2015) is approximately 400 000 code lines for the MITgcm and approximately 20 000 for the BFM.The coupler is a package that handles the interface (APIs) between the host code (MITgcm) and the BFM to solve Eqs. ( 7)-( 8) and to efficiently manage the matrices that contain the variables and tendencies shared by the two models and the flow of information among the different sub-model components.
The MITgcm-BFM coupling (Fig. 2) was achieved by upgrading a few routines of the MITgcm GCHEM package (Geo-CHEMistry, details in Appendix A), which handles the evolution of tracers, and by developing an additional package, BFMCOUPLER, which was specifically designed as the interface with the BFM model.The BFM is called by the MITgcm as an external library; therefore, the BFM was compiled separately using the same compiler used for the MITgcm (additional details on the compilation options and instructions are provided in Appendix A).
The BFMCOUPLER package (dashed box in Fig. 2) manages the initialization and memory usage of the BFM.This package also calls the BFM core routines and solves several processes that are not included in either model.The interfaces among the different components of the coupled model were designed so that the tracer transport sub-model (MITgcm PTRACERS package) uses the u, v, and w components As an interface, the BFMCOUPLER manages the transfer of information that is required by the BFM from both the hydrodynamic and transport sub-models of the MITgcm, and provides the integration solver (a MITgcm package) with the biogeochemical surface and bottom forcing and the sink/source terms originating from the BFM (gTracer bio ).The values of the tracers are derived from the transport sub-model.Moreover, the hydrodynamic sub-model supplies the temperature, salinity, and photosynthetic active radiation (PAR) values as well as additional forcing parameters (presence of ice, wind speed, and air partial pressure of CO 2 ) and information, such as the position of the specific grid cell within the water column (surface, intermediate, or bottom), which activates specific processes (e.g.surface air-sea gas transfer or bottom sediment fluxes).Then, the BFM calculates the biological partial derivative of Eq. ( 10) (fourth term), and the BFMCOUPLER returns this term to the time integration package, which integrates the transport and biogeochemical derivative terms to solve Eq. ( 10).Certain information used by the BFM, such as the PAR and wind values, can be calculated directly from the internal variables of the hydrodynamic sub-model (such as short-wave radiation (Q sw ) or other atmospheric fields) or managed from external sources.

Integration scheme, operator splitting, and LONGSTEP options
We considered several coupling strategies according to the MITgcm's code structure (Fig. 3).Within each time step of the model integration, which is coded in the FORWARD_STEP routine, the MIT- gcm solves the hydrodynamic equations (Eqs.1-6) through several routines: DO_ATMOSPHERIC_PHYS, DO_OCEAN_PHYS, DYNAMICS, TEMP_INTEGRATE, and SALT_INTEGRATE; further adjustments for temperature and salinity (e.g.filters) are applied in TRACERS_CORRECTION_STEP (Adcroft et al., 2016).Different options can be used to solve the evolution of tracers (Eq.10), which can be controlled by the gchem_separate_forcing pre-compilation option.When this option is false (#ifndef in Fig. 3), a direct integration scheme is adopted; therefore, the transport (gTracer trsp ) and biogeochemical (gTracer bio ) tendencies of each tracer are calculated by using the same (current) values of the physical and biogeochemical variables: where v, θ , S, and F C are the hydrodynamic variables and the additional biogeochemical forcing and t is the time discretization, which is the same adopted by the hydrodynamic sub-model.The biogeochemical tendency, which is solved by calling the BFM through the BFMCOUPLER_CALC_TENDENCY routine, is temporarily stored in the gchemTendency matrix, which is then summed to the overall tracer tendency, gTracer, in the PTRACER_APPLY_FORCING routine, along with the tendency term from the evaporation minus precipitation minus runoff effect (surfPtracers).The transport terms of the tracers (which update the gTracer matrix) are subsequently calculated within the PTRACER_INTEGRATE routine by several subroutines (GAD_ADVECTION, GAD_CALC_RHS, IMPLDIFF, and others) according to the options and numerical schemes selected in the specific MITgcm simulation set-up.The TIMESTEP_TRACER routine calculates the integration of Eq. ( 10) by providing a new state for the tracers.However, when the MITgcm set-up is prescribed with an implicit vertical diffusion scheme, an update of the state of tracers is solved within the IMPLDIFF routine according to the specific parameterization of the vertical diffusion (e.g.KPP, GGL90).Finally, if open boundary conditions are prescribed in the MITgcm set-up, the OBCS_APPLY_PTRACER routine applies the updated values of the tracers at the boundaries.The calculation of the derivative of the transport processes (gTracer trsp ) involves several MITgcm packages (GENERIC_ADVDIFF, PTRACERS, GCHEM, OBCS, KPP, and EXF) and options (choice of the advection scheme, viscosity and diffusivity coefficients, parameterization of surface dilution/concentration of tracers from evaporation, precipitation, and runoff), which are exhaustively described in the MITgcm documentation (Adcroft et al., 2016).
For the second coupling option, an operator splitting scheme is selected when gchem_separate_forcing is true (#ifdef).In this case, the biogeochemical tendency (gTracer bio ) is calculated after the state of the tracers has been updated by the transport equation terms.An intermediate value of the tracers, C n+1 , is passed to BFMCOUPLER_CALC_TENDENCY along with the values of the updated hydrodynamic variables (Eq.12).
This option allows for the development of an integration scheme with different time steps for the hydrodynamic and transport parts on one side and for the biological processes on the other.
A third option is an operator splitting algorithm, which involves the MITgcm LONGSTEP package (Adcroft et al., 2016) and adopts different time steps for the hydrodynamic and transport-biogeochemical components, thus increasing the computational performance of a coupled simulation.In particular, the tracer time step is set as a multiple (LS n ) of the main (hydrodynamic) time step, t trc = LS n • t, whereas the terms v, θ , S, and F C in Eq. ( 11) are replaced by suitable averages of the physical variables.The calculation of the averages is controlled by the LS_when_to_sample parameter, which defines the position within the code workflow in which the hydrodynamic variables are sampled (longstep_average, Fig. 3) and biogeochemical tracer tendencies are calculated (LONGSTEP_THERMODYNAMICS, Fig. 3).To activate this option, the LONGSTEP package code must be modified properly: the LONGSTEP_THERMODYNAMICS routine must be modified by adding a call to the modified GCHEM_CALC_TENDENCY routine.
This third method is preferred over the previous one as a possible method of decoupling the numerical biogeochemistry solution from the hydrodynamic solution.We tested the model to verify the trade-off between the increase in computational performance and the loss of accuracy in the model results as a function of the extension of the time step for the tracer equations (LS n ; see Sect.3.1.3).

BFMCOUPLER processes
The core of the present coupling scheme is the new BFMCOUPLER_CALC_TENDENCY routine, which is called by GCHEM_CALC_TENDENCY or by GCHEM_ FORCING_STEP.The approach adopted in this coupling scheme is to loop in space and to call the BFM as a subroutine to calculate the derivative terms of each biogeochemical tracer for each computational grid point (gTracer bio in Fig. 2).The derivatives of the chemical and biological processes are calculated by the BFM model via an Euler forward scheme through the BFM0D_ECOLOGY_DYNAMICS routine (a BFM routine) and stored in the 4-D MITgcm gchemTendency matrix.Additionally, the contributions of other processes, which are not explicitly coded in the BFM, are solved within BFMCOUPLER_CALC_TENDENCY, namely, the light penetration formulation, the sinking of phytoplankton and detritus, and the exchange processes in the surface and bottom layers of the water column.
In particular, the BFMCOUPLER package calculates the vertical profile of PAR along the water column, starting from the surface PAR, which is read from an external file or by using the short-wave radiation field (Q sw ), which is converted into PAR by a standard bulk formula if the native MITgcm atmospheric forcing package EXF is active (Britton and Dodd, 1976): www.geosci-model-dev.net/10/1423/2017/Geosci.Model Dev., 10, 1423-1445, 2017 where PAR s is the PAR at the sea surface, conv is a conversion factor of 4.6 µEin m −2 s −1 (W m −2 ) −1 , and pfrac is the fraction of the radiation in the visible band, which equals 0.4.The calculation of PAR along the water column, Eq. ( 14), is performed in the cell centre according to the Lambert-Beer formulation for the light exponential decay with depth and the shading of detritus and phytoplankton: where Chl j is the chlorophyll concentration of the j th phytoplankton functional type (PFT), R C is the carbon concentration of the detritus or optically active organic matter, and Kp j and K R are the corresponding extinction factors.K ext represents a background value set constant and is equal to 0.035 m −1 (considering pure water), or it can be read from an external file.In the latter case, the external file contains maps of the background extinction factor, which can be built a priori to incorporate the contributions of different unparameterized processes (e.g. the pattern distribution of yellow substances).
BFMCOUPLER solves the sinking processes and is activated for the phytoplankton groups and detritus (R C,N,P variables in Fig. 1): where w s is the sinking velocity (m s −1 ), which is provided both as a constant value or as a diagnostic result produced by the BFM model based on the nutrient stress conditions of the phytoplankton cells (Lazzari et al., 2012).The equation is solved numerically based on an Euler forward scheme.
A second module of BFMCOUPLER was designed to easily handle the boundary conditions at the surface and bottom.At the surface, air deposition can constitute an important source of nutrients in oligotrophic systems.Furthermore, when the runoff and nutrient discharge from rivers cannot be incorporated into the MITgcm OBCS package (as in Cossarini et al., 2015a), incorporating these factors as external surface forcings may be necessary (i.e. as localized runoff).Therefore, such contributions are prescribed as additional terms in gchemTendency in the surface layer by reading timevarying 2-D maps from external files: The coupled MITgcm-BFM model includes a simple parameterization of the fluxes at the water-sediment interface, which includes the burial of detritus (e.g. a net export flux from the ecosystem) and an incoming flux of nutrients into the deepest cell of the water column.Burial is parameterized as the first-order kinetics of the carbon (C), nitrogen (N), and phosphorus (P) content in the detritus (R (3) C,N,P variables), which is exported out from the bottom grid cell: In the same grid cell, the nutrient (for C equal to N (1) , N (2) and N (3) in Fig. 1) bottom fluxes are set either as a constant rate over the entire domain or as time-varying 2-D maps that can be read from an external file or provided by the benthic module of the BFM (which is foreseen in an ongoing development of the model): The BFMCOUPLER involves the use of several external surface forcing fields, such as the surface photosynthetic active radiation, background light extinction factor, sediment fluxes, and partial pressure of atmospheric carbon dioxide, which are used by the BFM to calculate the air-sea CO 2 exchanges.These fields are managed by the BFMCOUPLER package through the BFMCOUPLER_FIELDS_LOAD routine, which is a specifically modified replica of the EXTERNAL_FIELDS_LOAD routine of the MITgcm (Adcroft et al., 2016).This reading of external fields is controlled by two parameters: the timespan (forcingCycle) and the frequency (forcingPeriod) of external forcing, which are specified in the BFMCOUPLER namelist (additional details in the Appendix).

Compilation and set-up
The MITgcm and BFM must be compiled with the same compiler.We tested the code by using both the GNU and Intel compilers on several HPC platforms.Here, we report the results obtained by running the model (compiled with Intel) on a Linux cluster.The BFM is compiled as an independent library by using the following option of the BFM makefile: mkmf -p $BFM_LIB, and by configuring the config_BFM.shcompiling bash script with the appropriate compilation options (modules, optimization, and compiler), which are also used for the MITgcm compilation.Then, the build_option file for generating the MITgcm makefile must be modified by adding the path to the BFM compiled library and include files.Additional details are given in the manual of the BFMCOUPLER package (Appendix A).

Results
We tested the new coupled hydrodynamic-biogeochemical model against two case studies: an idealized experiment (a cyclonic gyre in a mid-latitude closed domain) and a realistic configuration (central Mediterranean Sea).In the first case study, which was released along with the code and the manual (https://github.com/gcossarini/BFMCOUPLER),we aimed to test the coherence of the model with the expected dynamics based on theoretical considerations and to test the model's performance under different coupling configurations.The second application was not meant to produce a thoroughly validated description of the dynamics of the area, but has been designed to show that the coupled model (once run in a realistic set-up) can be used to investigate a wide range of processes from coastal areas to open ocean.A thorough quantitative validation of the Mediterranean model output and an exploration of the results for analyses of the biogeochemical dynamics in the area are beyond the scope of this paper.

Idealized case study
This experiment was based on a simplified case study that consisted of an idealized domain (2 • × 2 • × 280 m closed box) that was forced by steady winds and a seasonal cycle of surface heat (downward long-wave and short-wave radiation) and mass (precipitation) fluxes.The horizontal shear in the surface wind field maintained a permanent cyclonic gyre, whereas the surface heat fluxes acted on the thermohaline properties of the water column, inducing a yearly cycle (summer stratification -winter mixing).This simulation was run for several years to reach steady-state conditions (perpetual year simulation).

Numerical configuration
This domain was discretized by adopting a uniform grid spacing (1/32 • ) in the horizontal direction, creating 64 grid cells in both directions.All the peripheral grid points of the bathymetry were land points (closed box), whereas the bottom of the domain was a bowl-shaped pit.In the vertical direction, the model was composed of 30 layers with non-uniform thickness (from 1.5 to 21 m).The time step equalled 300 s.External forcing fields were introduced via the MITgcm native EXF package.The meteorological forcing consisted of nine surface fields, namely, the 2 m air temperature (atemp), 2 m specific humidity (aqh), 10 m zonal and meridional wind (uwind, vwind), precipitation (precip), long-and short-wave incident radiation (lwdown, swdown), air pressure (apressure), and surface runoff (runoff).The wind stress and total heat flux were calculated via standard bulk formulae.The experiment was designed with no open boundaries to verify the mass conservation of chemical elements and simulate the effect of free surface dynamics on the distribution of tracers in the surface layer, which can be important for certain processes, such as the effects of concentration and dilution on the carbonate system variables.We chose the pre-compilation option, which allows for the presence of mass sources/sinks of fluid in the domain (3-D generalization of the oceanic real-freshwater flux option: #define ALLOW_ADDFLUID).With this option enabled, the net contribution of precipitation, evaporation, and runoff can be considered in the total mass budget.In particular, we activated the "exact conservation" of fluid in the free-surface formulation (#define EXACT_CONSERV) so that the temporal evolution of the free surface height exactly equalled the divergence of the volume transport.We allowed the use of the non-linear free-surface option so that the surface level thickness (hFactor) could vary with time (#define NONLIN_FRSURF).The tests were run by adopting the following runtime options (in the "data" namelist) for the free-surface formulation and the volume conservation constraints: &PARM01 implicitFreeSurface=.TRUE., exactConserv=.TRUE., useRealFreshwaterFlux=.TRUE., selectAddFluid=1, linFSConserveTr=.FALSE., nonlinFreeSurf=4, &END When configuring the options for the passive tracers package (PTRACERS), we set the concentrations of the tracers in the surface mass fluxes (evaporation minus precipitation minus runoff) to always equal zero (PTRACERS_EvPrRn(tracer_number) = 0.0).We used the same advection scheme (third order and direct space-time with a Sweby flux limiter) for active and passive tracers (tracerAdvScheme = 33).Biogeochemical variables were initialized with suitable vertical profiles for winter conditions all over the domain.The BFMCOUPLER package was configured without external forcing both at the surface and at the bottom for nutrients, so that a closed system is simulated and mass conservation is checked.PAR was converted from short-wave radiation, and the light extinction factor was calculated considering a background value (K ext = 0.035 m −1 ) and the shading effect of phytoplankton groups.All details of this experiment along with namelists and input files are given in Appendix A.

Results of the simulation
The model simulated a realistic cyclonic circulation with associated mesoscale variability from vertical thermohaline stratification and flow instability.Relatively well-mixed thermohaline conditions in the winter induced a more unstable cyclonic gyre with small-scale mesoscale eddies (Fig. 4a), whereas a more stable and energetic cyclonic circulation occurred from stratified thermohaline conditions in the summer (Fig. 4b).
Figure 5 shows the evolution of several physical properties and biological components within the central part of the gyre.The coupled model simulated the evolution of the therwww.geosci-model-dev.net/10/1423/2017/Geosci.Model Dev., 10, 1423-1445, 2017  mocline and nutricline and the effect of winter vertical mixing on the temperature and nutrient profiles (Fig. 5a and b). Figure 5 also shows the formation of surface phytoplankton blooms during early winter (Fig. 5c), the formation of the deep chlorophyll maximum (DCM) during summer (as a trade-off between the light penetration and the depth of the nutricline), and the effect of the erosion of the stratification during autumn on the biogeochemical properties of the basin (deepening of mixed layer depth -MLD -Fig.5a).Net pri-mary production (NPP, contour plot in Fig. 5d) showed the highest values in the proximity of the DCM during spring, although high primary productivity was also simulated in the upper part of the water column, where the high level of irradiance stimulated carbon fixation, especially for small-sized phytoplankton groups (not shown), even in the presence of low phytoplankton biomass.
The region close to the DCM was the most active biological area, i.e. the concentrations of all of the living variables (small and mesozooplankton groups and bacteria; Fig. 5e  and f) were the highest and the fluxes fuelled the so-called classic food chain (Legendre and Rassoulzadegan, 1995).Nevertheless, significant bacterial biomass was also simulated in the upper part of the water column, where bacteria consumed the labile organic matter, which was sideproduced by phytoplankton in the well-lit upper levels.Small zooplankton (sum of micro-and hetero-trophic nanoflagellate groups) took advantage of the bacterial biomass, triggering the so-called microbial food web (Legendre and Rassoulzadegan, 1995), which dominated the upper part of the water column during summer.Oxygen (Fig. 5d) was higher in the upper part of the water column during winter because of the high level of NPP and the effect of re-aeration processes with the atmosphere.Bacterial production and the predominance of respiration over phytoplankton photosynthesis caused the autumn minimum.

Application of the LONGSTEP option
The computational cost of a 1-year simulation was approximately 5 h when adopting an MPI configuration that featured 16 Ivy Bridge cores.The code profiling (Fig. 6) indicated that most of the CPU time (i.e. up to 85 %) was devoted to solving the differential equation for the high number of tracers (51).Solving the transport part (Tracers trsp in Fig. 6) of the tracer equation (Eq.10) accounted for 50 % of the overall computational cost, whereas solving the biological part (Tracers bio in Fig. 6) accounted for 35 %.The cost of solving tracer transport increased linearly with the number of tracers (e.g.Tracers trsp is almost 25 times larger than the time used to solve for temperature and salinity; Fig. 6), whereas the cost of the BFM calculations was primarily dependent on the solution of the carbonate system, although the complexity of the relationships among the biogeochemical variables (results not shown) was also a factor.The use of the LONGSTEP MITgcm package caused an almost exponential reduction in the computational cost for the integration of the tracer equation (Fig. 6).With a LS n set to 8 (a time step for tracers, t trc , equals 2400 s), the runtime was reduced by more than 80 % with respect to the reference run.Assuming this optimization, the fraction that was devoted to the solution of the hydrodynamic and MPI routines accounted for 45 %, whereas the remaining part (55 %) was devoted to solving the transport-biogeochemical part.Within the tracer equation, 60 % of the quota was allotted for transport and 40 % of the quota was allotted for the BFM and the other biogeochemical processes.
The use of a coarser time resolution for the solution of the tracer equations implied errors with respect to the reference solution (Fig. 6).The errors were calculated as the root mean square of the difference of the integrated 0-200 m chlorophyll between the reference run (LS n = 1, i.e. no LONGSTEP) and the run with increased LS n .The magnitude of the mean annual error increased almost linearly with the coarsening of t and equalled 0.0025 mg m −3 at LS n = 8.Within a simulation, the largest differences between the reference run and the coarser time discretization run were registered during periods with the highest chlorophyll tendency, such as during autumn vertical mixing events along the entire water column and during the deep chlorophyll maximum formation in the spring (not shown).The errors became relevant (> 0.01 mg m −3 ) when larger values for LS n (e.g.LS n ≥ 10) were adopted.

Mass budget
The reference run was also used to verify the mass conservation of the coupled hydrodynamic-biogeochemical model by considering that the model configuration (i.e.non-linear free surface) was set to properly simulate the effects of freesurface dynamics on the concentrations of the biogeochemical variables at the surface.Figure 7 shows the time series of the sea surface height (SSH) averaged over the entire basin.The results indicated the prevalence of rain over evaporation for the first part of the year and vice versa from May to October.For example, the evolution of alkalinity, which is a key variable for resolving the ocean carbonate system (Follows et al., 2006), was correctly anti-correlated with the derivative of www.geosci-model-dev.net/10/1423/2017/Geosci.Model Dev., 10, 1423-1445, 2017 SSH in the surface layer because the effects of concentration and dilution at the surface are dependent on the water mass balance.This model feature was provided along with the mass conservation capability for tracers (Fig. 7).The errors in mass conservation over time were small (O(10 −9 )) and they were caused by the computation of the time average of the model output.The coupled MITgcm-BFM model, which was configured with the non-linear free-surface option, allowed us to efficiently simulate the dilution-concentration dynamics while preserving the ability to calculate the budget of the chemical elements with a high level of accuracy.This feature is indeed important considering the dynamics of variables like alkalinity, whose spatial patterns at the surface were dominated by the regional-spatial-scale distribution of the water mass budget (Cossarini et al., 2015b).

Adriatic-Ionian system case study
The coupled model was also used to simulate a realistic domain: the central Mediterranean Sea.This area, which encompasses the Adriatic and Ionian seas (Fig. 8), was chosen because it is characterized by a wide range of interconnected ecosystems that span coastal areas, which are influenced by river discharges, and offshore regions, which are characterized by open-sea dynamics.Indeed, the northern part of the Adriatic is a continental shelf area influenced by terrestrial input (Solidoro et al., 2009;Cossarini et al., 2015a).This area is a site of dense water formation (Gačić et al., 2001;Querin et al., 2013) and represents one of the most productive areas of the Mediterranean Sea (Mangoni et al., 2008).
The southern Adriatic Sea is characterized by an almost permanent geostrophic gyre modulated by deep winter mixing episodes (Gačić et al., 2002;Bensi et al., 2014), and it is connected to the Ionian Sea via the Otranto Strait.The Ionian Sea is the deepest sub-basin of the Mediterranean, and it is characterized by basin-scale circulation patterns and smaller mesoscale eddies.This sea is influenced by oligotrophic and salty waters originating from the Levantine basin and by the relatively fresh Atlantic water masses that flow from the west.The hydrodynamics of the area have been simulated by the Adriatic-Ionian implementation of the MITgcm (ADriatic IOnian System model (ADIOS), Querin et al., 2016), which we used in this study.The aim of this experiment is to show the ability of the new coupled model to properly simulate the effects of hydrodynamics on biogeochemistry within a wide range of oceanographic and ecological processes that range from a few kilometres to hundreds of kilometres and from oligotrophic to high-level trophic conditions.

Domain and model set-up
The model domain was delimited by the Sicily channel (lon 12.2 • E) on the western side and by the Cretan Passage (lon 22.7 • E) on the eastern side.The Strait of Messina and the Gulf of Corinth were excluded in this study.The horizontal resolution was 1/32 • (approximately 3 km), whereas the vertical grid consisted of 72 z levels; therefore, the ADIOS model could be easily nested in the 1/16 • Coperni-  The model set-up only considered the main rivers that flow into the Adriatic Sea, whereas the minor contributions that flow into the Ionian Sea were neglected.River contributions were introduced as local boundary conditions, imposing observed daily freshwater flow rates for the major rivers (e.g.Po) and climatological annual flow rates for the others, with spring and autumn maxima and winter and summer minima (Querin et al., 2013;Janeković et al., 2014).The tracer concentrations at the river mouths were constant in space and time (Table 1), and the mass fluxes were calculated by multiplying the concentrations by the flow rate of each river.
The boundary conditions along the Sicily Channel and along the Cretan Passage were derived from the CMEMS MED-MFC system (Tonani et al., 2008;Lazzari et al., 2010) for both the hydrodynamic and biogeochemical variables (OBC and OBC C in Fig. 2).The output of the 1999-2012 reanalysis (Salon et al., 2015) was downloaded from the Copernicus web portal (http://marine.copernicus.eu).The present model configuration adopted a finer horizontal resolution (from 1/16 • to 1/32 • ) with respect to the CMEMS MED-MFC system, whereas the vertical spacing was the same; hence, interpolating/extrapolating the hydrodynamic and biogeochemical fields in the vertical direction was unnecessary.Furthermore, both the CMEMS MED-MFC system and ADIOS adopted the BFM biogeochemical model; therefore, changes or conversions to the biogeochemical variables were not required.The initial conditions for the hydrodynamic and biogeochemical variables were also derived from the CMEMS MED-MFC system by linearly interpolating the original fields from 1/16 • to 1/32 • .Additional details on the ADIOS model set-up are provided by Querin et al. (2016).
Surface meteorological forcing was derived from the Regional Climate Model (RegCM) developed at the International Centre for Theoretical Physics (ICTP) in Trieste.We used the 12 km horizontal resolution version with 3 h output frequency (as in Querin et al., 2016).The heat fluxes (Q θ in Fig. 2) at the air-sea interface were calculated using standard bulk formulae (via the MITgcm native EXF package); the air temperature, specific humidity, precipitation, incoming radiation, and wind speed values were interpolated from the meteorological model; the sea surface temperature was provided by the oceanographic model.The 3 h temporal resolution can highlight the daily variability in the physical and biogeochemical properties of the uppermost layers of the water column (daily cycling of the PAR, temporal variability in the temperature and wind).
The specific settings for the BFMCOUPLER package were specified as follows.The background water light extinction factor was set considering a longitudinal negative gradient according to Lazzari et al. (2012) and the coefficient for the self-shading effect was set to 10 −3 and 8 × 10 −3 m 2 mg −1 chlorophyll for diatoms and the other three phytoplankton groups, respectively.The nutrient surface forcing (air deposition) was set to 0.00096 and 0.057 mmol m −2 d −1 for phosphorus and nitrate, respectively (Lazzari et al., 2012 and reference therein), whereas we assumed that the atmospheric carbon dioxide (pCO atm 2 ) linearly increased from 380 to 395 in the period 2006-2012 according to the trend that was reported in Artuso et al. (2009).No bottom forcing was prescribed for the biogeochemistry.

Results of the simulation
The simulation covered the period from January 2006 to December 2012 at a time step of 200 s.In the following analysis, we disregarded the first 2 years of the simulation, which we considered a spin-up period for the biogeochemical variables from the CMEMS's coarser resolution fields.The MPI domain decomposition consisted of 16×14 subdomains run on 224 Intel Xeon Ivy Bridge cores of a Linux cluster, and the computational cost of the simulation was 65.8 h per year.The runtime was significantly reduced by adopting the LONGSTEP option (Table 2).The wall clock time progressively decreased by increasing the LONGSTEP factor (LS n ) from 1 to 9.Then, time steps that were higher than 30 min www.geosci-model-dev.net/10/1423/2017/Geosci.Model Dev., 10, 1423-1445, 2017 substantially decreased the accuracy without further reducing the computational cost (Table 2).We present the results for the ADIOS case study to demonstrate the ability of the new MITgcm-BFM coupled model to investigate closely interconnected hydrodynamic and biogeochemical processes for both coastal and open-sea ecosystems.
In the western coastal areas of the Adriatic Sea, the maps in Fig. 9 correctly display the patterns of low salinity, southward currents, high nitrate and chlorophyll concentrations, and strong primary production, which are all typical fingerprints of the Western Adriatic Current (WAC) system in the Adriatic Sea.The effect of the input from the northern rivers and the basin-scale cyclonic circulation generates a frontal system along the Italian coast.As is commonly observed in satellite chlorophyll maps (Barale et al., 2008), the width of the WAC frontal system decreases southwards, whereas weaker recirculation patterns are also visible in the central Adriatic Sea (Fig. 9).Other river-influenced coastal areas are simulated along the south-eastern areas of the Adriatic Sea, where the input from the Neretva and other south-eastern rivers triggers small-scale chlorophyll a signals along those areas, as reported by Marini et al. (2010).The northward flow of salty and oligotrophic water, which enters through the Otranto Strait, confines the river's fertilization to a narrow coastal strip.
The coastal to open-sea gradients of nutrients were accurately simulated by the coupled model.As an example, Fig. 9 shows that the nitrate patterns display a longitudinal gradient along the Adriatic and northern Ionian seas, and these results are consistent with the current climatologies (Cossarini et al., 2012;Solidoro et al., 2009;Zavatarelli et al., 1998).In the open-sea area of the Ionian Sea, the surface circulation is dominated by large mesoscale structures and a basin-scale anticyclone in the middle, and the downwelling area is characterized by minimal nitrate and chlorophyll concentrations (Fig. 9).This pattern is consistent with the climatology of Manca et al. (2004), even if the nitrate concentrations are slightly higher in the eastern Ionian Sea, which is related to overestimated eastern boundary values.
If we focus on the open-sea sub-surface dynamics, we can analyse how vertical processes affect the biogeochemistry.The vertical profiles of chlorophyll and phosphate for the two sites in Fig. 8 are depicted in Fig. 10.One site is located in the centre of the southern Adriatic gyre, which is characterized by strong winter vertical mixing, whereas the second is located in the centre of the large anticyclonic gyre in the Ionian Sea.A comparison between the two sites shows the ability of the coupled model to simulate the different regimes in the two areas.The southern Adriatic Sea presents a much higher mixed layer depth in winter, a shallower nutricline than the Ionian Sea, more intense inter-annual variability in the cyclic alternation of winter vertical mixing phases, and the onset of summer stratification.
The intense vertical mixing in the southern Adriatic area during winter drives the upwelling of nutrient-rich water, which contributes to a shallow nutricline (up to the depth of the DCM) during summer.However, winter ventilation in the Ionian Sea's open areas rarely reaches a depth of 250 m; consequently, nutrient-rich water remains confined to the deepest layers (below 200 m).The two areas are characterized by different biological regimes because of the different depths of the nutricline and the superimposed longitudinal gradient of the background light extinction factor (according to Lazzari et al., 2012).
Another interesting coupled hydrodynamicbiogeochemical feature is displayed along the southern coast of Sicily, where the entrance of modified Atlantic water (MAW, low-saline water mass in Fig. 9a) and the simulated coastal upwelling from westerly winds induce vertical transport of nutrients, consistent with the findings of Patti et al. (2010) and Rinaldi et al. (2014).Intense vertical dynamics trigger the high concentrations of nutrients and chlorophyll and the strong primary production simulated in the upper layer of the northern Sicily channel (Fig. 9b), and these results are consistent with the typical patterns observed in satellite chlorophyll maps (Volpe et al., 2012).
The computation and diagnostics of the transport components for the tracers (e.g.zonal and meridional advection and diffusion, vertical advection and implicit and explicit diffusion) are already implemented in the native PTRACERS and DIAGNOSTIC packages of the MITgcm.This feature, which is complemented by the ability to calculate the surface and lateral fluxes at the boundaries through the BFMCOUPLER package, allows us to calculate the budget of the simulated chemical elements in marine ecosystems.As an example, we evaluated the meridional transport across the Otranto Strait for the carbon components along with other fluxes  , Z (1,2,3,4) C , and B C ; and detritus, R C ) is mainly confined to the surface layer for both the inflow and outflow.A barely visible flux of organic carbon toward the Ionian Sea is depicted along the western slope below a depth of 200 m (mainly because of the sinking of detritus).The northward and southward fluxes of DIC along the surface (Fig. 11) are characterized by the same organic carbon pattern and nearly balanced.Additionally, an outflow (blue) area at a depth of 300-900 m along the left flank of the strait indicates DIC transport associated with the Adriatic Dense Water Outflow Current (DWOC, Gačić et al., 2001).This carbon flux represents the export term that closes the budget of the Adriatic Sea and replenishes the layer of the Ionian Sea below the depth of the Levantine Intermediate Water, which suggests a possible mechanism for the long-term carbon sequestration in the Mediterranean Sea.

Discussion and conclusions
In this paper, we presented a coupling between two widely used models, the MITgcm and BFM, and we showed the potential of the new coupled model.These two models were developed by two different scientific communities that are actively and constantly involved in improving the codes.When one model is directly embedded in another, code developments might represent an issue because of the constant and tedious work of keeping one code updated with respect to the other.Therefore, the coupling in this paper was designed to preserve the independence of the two models as much as possible.The number of modifications that were required for the two original codes was limited, and changes could be easily managed should each single model be upgraded.In our solution, the MITgcm remained the host code, the BFM was compiled and linked as an independent library, and the new BFMCOUPLER package handled all the coupling procedures and concentrated all the coding effort.The upgrades to the MITgcm enumerated less than 10 new code lines in a few routines (in the GCHEM and LONGSTEP packages) and the list of available diagnostics (in the DIAGNOSTIC package).On the BFM side, several "include" files contained a list of newly added variables.The order of the variables in the BFM's include files and in the MITgcm's data.ptracerfile must be consistent (see Appendix A).This feature is important because the BFM (Vichi et al., 2015) can be customized in terms of the number of state variables and processes, thus increasing the flexibility of the new coupled model for a wider range of applications.
Despite the growth of computational resources, the efficiency of coupled codes can still be an issue because of the large size of the computational grids (Blom and Verwer, 2000).Domain decomposition and parallelization tools are available in several coupling environments (e.g.FABM, Bruggeman and Bolding, 2014;MESSy, Jöckel et al., 2008).Likewise, our coupling scheme has been thought to fully exploit the parallelization efficiency of the MITgcm (Marshall et al., 1997), and no additional coding effort (in terms of parallelization) is required by the users.
Other biogeochemical models of various complexity have already been embedded in the MITgcm (Dutkiewicz et al., 2009;Hauck et al., 2013;Cossarini et al., 2015a).Nevertheless, the BFM in this new coupled model has a biological complexity and a number of features (Lazzari et al., 2016) that increase the attractiveness of the model for many marine applications.
The MITgcm-BFM coupling scheme was primarily designed by considering the direct integration scheme because this framework has the highest level of numerical accuracy.The use of the LONGSTEP option reformulated the coupling as an operator splitting algorithm that allows for different time steps for hydrodynamics and coupled transport-biogeochemistry at the cost of accuracy.When using the LONGSTEP option, the results (Fig. 6 and Table 2) show that the loss of accuracy remained negligible only for a limited increase in the tracer time step.Furthermore, the coupling framework could handle a separate solution of hydrodynamics and transport processes from the biogeochemical processes through the use of the gchem_separate_forcing option (Fig. 3).However, this approach would require a wider modification of the GCHEM package to introduce independent integration steps for the transport and biogeochemical parts of the tracers.Then, a more detailed analysis of the sensitivity (e.g.similar to what was proposed in Butenschön et al., 2012) of the biogeochemical model's results to the different coupling schemes and time steps should be performed for each specific application.
A direct integration scheme might be more appropriate for investigating the feedback of the biogeochemistry on the hydrodynamics of the system.An example is the calculation for the sinking of certain phytoplankton groups, which is a physical 1-D process solved within BFMCOUPLER and related to the sinking velocity calculated by the BFM.Furthermore, the shading effect on light penetration caused by phytoplankton and other suspended matter currently only affects the PAR vertical profile (Eq.14).However, this factor could be introduced as an extra term in the routine that calculates seawater thermodynamics (in the SWFRAC and EXTERNAL_FORCING routines).A new parameterization of the penetration of solar radiation could be used to estimate the biological effects on the seawater temperature, which might be an interesting issue in highly productive areas, such as the northern Adriatic Sea and the coastal strip along the Italian coast reached by the Western Adriatic Current (WAC).A realistic simulation of light absorption with depth could reduce the model errors when estimating temperature, which is affected by many other sources of uncertainty originating from the surface forcing data, the heat flux bulk formulation, the vertical resolution, and the parameterization of vertical turbulent processes.The design of our coupler, which is characterized by the sharing of biogeochemical variables and their tendencies in the host model's memory structure, allows for the future implementation of the feedback effects of biology on hydrodynamics.Furthermore, the new coupling scheme was designed to foster development towards a full Earth system modelling approach, in which a wide range of processes among the Earth's spheres can be simulated online and the interactions and feedback effects can be directly considered.For example, the BFM has already been coupled with other ecosystem components (e.g.online coupling with high-trophic-level model Ecopath with Ecosym, Akoglu et al., 2015).Moreover, the parameterization of Eqs. ( 16) and ( 17) can be easily substituted by a call to a benthic model function, which solves the processes that occur in a single-layer sediment model and calculates the exchanges between the pelagic environment and the sediment.
Similarly, the MITgcm has already been coupled with atmospheric models.For example, the MITgcm has been coupled online with the RegCM atmospheric model in the Mediterranean Sea region (Giorgi et al., 2006) using the OASIS coupling framework (Artale et al., 2010).Therefore, our coupling scheme can act as a link between atmospherehydrosphere models and biosphere models.This coupler could be successfully used to study ocean-atmosphere interactions, such as the effects of climate scenarios on hightrophic-level ecosystem components or the feedback of ocean carbon pumps on the climate.
Finally, the results of the two test cases show that the new coupled model provides a realistic representation of a wide range of marine processes from costal to open-sea ecosystems, where the interplay of hydrodynamics and biogeochemistry is crucial.The effects of river plumes, coastal upwelling, and different vertical mixing regimes on phytoplankton dynamics were reasonably reproduced by the model and found to be consistent with both theoretical knowledge (Mann and Lazier, 2006) and published experimental findings for the Mediterranean Sea.
Code availability.The code can be downloaded from the link: https://github.com/gcossarini/BFMCOUPLER/tree/release-1.0 (SHA-1 hash 02ce96dc), which corresponds to version v1.0 described in the present paper.The numerical experiment described in Sect.3.1 consists of an idealized domain forced by steady wind and a seasonal cycle of surface heat and mass fluxes.This case study simulates a permanent cyclonic gyre with a yearly cycle of thermohaline and biogeochemical properties.The input files along with the MITgcm and BFM namelists of the experiment are available at https://github.com/gcossarini/BFMCOUPLER/tree/release-1.0/Input.The modified MITgcm files, located in the /input/modified_MITgcm_files directory, should be linked through the -mods option of the MITgcm builder (see Sect. 3.4 of the manual) in order to override the original MITgcm source files with the modified ones, when the code is built.
New diagnostic quantities are listed in the namelist in the data.diagnosticsparameter file, which specifies the frequency and type of output, the number of levels, and the names of all the separate output files.
The coupled MITgcm-BFM model can use a large number of tracers; therefore, increasing the ndiagMax parameter in diagnostics_size.hmay be necessary.The initialization of BFMCOUPLER diagnostics is provided by adding a call statement to BFMCOUPLER_INIT_FIXED.F in the GCHEM_INIT_FIXED.F routine.

A2.5 LONGSTEP
The LONGSTEP MITgcm package allows the tracer time step to be longer than the time step used by the hydrodynamic model.When this package is activated along with the BFMCOUPLER package, a new specifically developed version of the LONGSTEP_THERMODYNAMICS.F routine has to be used.The new version of this routine includes a call to BFMCOUPLER_CALC_TENDENCY.The BFMCOUPLER routines use the hydrodynamic variables stored in the LONGSTEP variables, which are either the averages or temporal sub-samplings of the variables of the master hydrodynamic model depending on the when_to_sample parameter set in the data.longstepnamelist file.

A2.6 Compilation and compile time flags
The BFM is a Fortran95 code and must be compiled separately as an external library in advance ($BFM_LIB/lib/libbfm.a).According to the BFM's manual, a compiled library version is obtained by customizing the BFM makefile (mkmf -p $BFM_LIB).The config_BFM.sh compiling bash script must contain build options (modules, optimization options, and compiler) that are consistent with those of the MITgcm compilation.
When the MITgcm is compiled, the build_options file must be modified and the following lines must be added: The subroutines of the new BFMCOUPLER package must be included in the /MITgcm/pkg/BFMCOUPLER folder, which can be added to the original source tree of the code.BFMCOUPLER must be specified in the packages_conf compile configuration file.
Several specific compile time flags are set in BFMcoupler_OPTIONS.h:USE_QSW: use Q sw from the MITgcm to calculate the photosynthetic active radiation (PAR).READ_PAR: read the PAR from a file set in data.bfmcoupler.USE_SHADE: include the role of phytoplankton and detritus in the calculation of the vertical profile of the PAR.READ_xESP: read the background light extinction factor from a file set in data.bfmcoupler.USE_SINK: use the calculation for the sinking of phytoplankton and detritus.USE_BURIAL: calculate the contribution of burial for detritus tendency at the bottom.USE_BOT_FLUX: use input sediment fluxes for nutrients at the bottom.BFMCOUPLER_DEBUG: activate a control on the tendencies calculated by the BFM.

A3 Do's and dont's
This package must be run with both PTRACERS and GCHEM enabled.This package is configured for a number of biogeochemical variables specified by the BFM model.Therefore, data.ptracersmust be configured accordingly (order of tracers equals what is specified in the ModuleMem.F90 file from the BFM code).This package must also be run with diagnostics enabled.

Figure 2 .
Figure 2. Description of the MITgcm-BFM coupling and interfaces among the different components of the coupled model.

Figure 3 .
Figure 3. Workflow of the MITgcm FORWARD_STEP routine.The boxes indicate the routines (darker colours denote dependency).The matrix of the tracers' state variables (pTracer), the overall tendency of the tracer (gTracer), and the tendency for the biogeochemistry only (gchemTendency) are also specified.The blue boxes indicate modifications to either the MITgcm code or the BFMCOUPLER routines, whereas the green boxes indicate BFM routines.The pre-compilation options (#) and omitted parts (. . . ) are also shown.

Figure 4 .
Figure 4. Idealized case study (circulation in a 2 • × 2 • × 280 m closed domain).Horizontal component of velocity: current speed (colour) and direction (vectors) at 12 m depth.Five-day average in winter (a) and summer (b), and yearly average (c).

Figure 5 .
Figure 5. Hovmöller diagrams of the (a) temperature and evolution of the mixed-layer depth (MLD), (b) phosphate and PAR, (c) chlorophyll (sum of the chlorophyll content in the four phytoplankton functional groups) and phytoplankton expressed in carbon biomass, (d) oxygen and net primary production (NPP), (e) small zooplankton (small zoopl) and mesozooplankton (mesozoopl), and (f) bacteria.

Figure 6 .
Figure 6.Wall clock time of the main MITgcm routines clustered in selected groups (left axes) as a function of the number of hydrodynamic time steps between tracer time steps (LS n ): initialization and I/O (sum of the routines: MODEL_I/O, DO_STATEVARS_DIAGS, LOAD_FIELDS_DRIVER, MONITOR, DO_THE_MODEL_IO, and DO_WRITE_PICKUP), hydrodynamics (sum of DYNAMICS, SOLVE_FOR_PRESSURE, INTEGR_CONTINUITY, and other routines); temperature and salinity (sum of the routines: TEMP_INTEGRATE and SALT_INTEGRATE); MPI tasks (BLOCKING_EXCHANGES routine); tracers bio (GCHEM_CALC_TENDENCY) and tracers trsp (PTRACER_INTEGRATE).The root mean square difference of the integrated chlorophyll (right axis) is shown as a function of LS n .

Figure 7 .
Figure 7. Evolution of SSH (blue line) and alkalinity (red line) at the surface layer together with the relative variation of total alkalinity mass (M) with respect to the initial condition (M 0 ) over the whole domain (black line).The total alkalinity mass was obtained by multiplying the daily average model output by the domain volume, which included the time-varying SSH at the surface layer.

Figure 8 .
Figure 8. Bathymetry (depth in metres) of the Adriatic-Ionian model.The plot also indicates the location of the major rivers (arrows), the Otranto Strait, and the position of the two sites (circles) that were selected to display the Hovmöller diagrams in Fig. 10.
Forecasting system (CMEMS MED-MFC; Lazzari et al., 2010), which shares the same bathymetry along the open boundary of ADIOS.

Figure 9 .
Figure 9. (a) Map of the surface currents (arrows) and salinity.(b) Map of the surface chlorophyll and contours (solid black lines) of nitrate concentration in the upper layer (0-20 m), and contours (dotted blue lines) of the annually averaged and vertically integrated (0-200 m) net primary production (NPP).

Figure 10 .Figure 11 .
Figure 10.Hovmöller diagrams of chlorophyll (colour) and phosphate (contour, mmol m −3 ) and plots of the mixed layer depth (dashed lines, m) for the southern Adriatic Pit (a) and the Ionian offshore area (b).

Table 1 .
Concentrations of tracers in the rivers.

Table 2 .
Computational cost as a function of the LONGSTEP factor (LS n ) and the mean error of the integrated chlorophyll.The error was the annual average of the rms of the differences between the LONGSTEP simulations and the reference (Ref) simulation.The error was normalized using the reference simulation.