New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
chap_SBC.tex in NEMO/trunk/doc/latex/NEMO/subfiles – NEMO

source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex @ 10354

Last change on this file since 10354 was 10354, checked in by nicolasmartin, 5 years ago

Vast edition of LaTeX subfiles to improve the readability by cutting sentences in a more suitable way
Every sentence begins in a new line and if necessary is splitted around 110 characters lenght for side-by-side visualisation,
this setting may not be adequate for everyone but something has to be set.
The punctuation was the primer trigger for the cutting process, otherwise subordinators and coordinators, in order to mostly keep a meaning for each line

File size: 85.1 KB
Line 
1\documentclass[../tex_main/NEMO_manual]{subfiles}
2\begin{document}
3% ================================================================
4% Chapter —— Surface Boundary Condition (SBC, ISF, ICB)
5% ================================================================
6\chapter{Surface Boundary Condition (SBC, ISF, ICB) }
7\label{chap:SBC}
8\minitoc
9
10\newpage
11$\ $\newline    % force a new ligne
12%---------------------------------------namsbc--------------------------------------------------
13
14\nlst{namsbc}
15%--------------------------------------------------------------------------------------------------------------
16$\ $\newline    % force a new ligne
17
18The ocean needs six fields as surface boundary condition:
19\begin{itemize}
20\item
21  the two components of the surface ocean stress $\left( {\tau _u \;,\;\tau _v} \right)$
22\item
23  the incoming solar and non solar heat fluxes $\left( {Q_{ns} \;,\;Q_{sr} } \right)$
24\item
25  the surface freshwater budget $\left( {\textit{emp}} \right)$
26\item
27  the surface salt flux associated with freezing/melting of seawater $\left( {\textit{sfx}} \right)$
28\end{itemize}
29plus an optional field:
30\begin{itemize}
31   \item the atmospheric pressure at the ocean surface $\left( p_a \right)$
32\end{itemize}
33
34Five different ways to provide the first six fields to the ocean are available which are controlled by
35namelist \ngn{namsbc} variables:
36an analytical formulation (\np{ln\_ana}\forcode{ = .true.}),
37a flux formulation (\np{ln\_flx}\forcode{ = .true.}),
38a bulk formulae formulation (CORE (\np{ln\_blk\_core}\forcode{ = .true.}),
39CLIO (\np{ln\_blk\_clio}\forcode{ = .true.}) or
40MFS \footnote { Note that MFS bulk formulae compute fluxes only for the ocean component}
41(\np{ln\_blk\_mfs}\forcode{ = .true.}) bulk formulae) and
42a coupled or mixed forced/coupled formulation (exchanges with a atmospheric model via the OASIS coupler)
43(\np{ln\_cpl} or \np{ln\_mixcpl}\forcode{ = .true.}).
44When used ($i.e.$ \np{ln\_apr\_dyn}\forcode{ = .true.}),
45the atmospheric pressure forces both ocean and ice dynamics.
46
47The frequency at which the forcing fields have to be updated is given by the \np{nn\_fsbc} namelist parameter.
48When the fields are supplied from data files (flux and bulk formulations),
49the input fields need not be supplied on the model grid.
50Instead a file of coordinates and weights can be supplied which maps the data from the supplied grid to
51the model points (so called "Interpolation on the Fly", see \autoref{subsec:SBC_iof}).
52If the Interpolation on the Fly option is used, input data belonging to land points (in the native grid),
53can be masked to avoid spurious results in proximity of the coasts as
54large sea-land gradients characterize most of the atmospheric variables.
55
56In addition, the resulting fields can be further modified using several namelist options.
57These options control
58\begin{itemize}
59\item
60  the rotation of vector components supplied relative to an east-north coordinate system onto
61  the local grid directions in the model;
62\item
63  the addition of a surface restoring term to observed SST and/or SSS (\np{ln\_ssr}\forcode{ = .true.});
64\item
65  the modification of fluxes below ice-covered areas (using observed ice-cover or a sea-ice model)
66  (\np{nn\_ice}\forcode{ = 0..3});
67\item
68  the addition of river runoffs as surface freshwater fluxes or lateral inflow (\np{ln\_rnf}\forcode{ = .true.});
69\item
70  the addition of isf melting as lateral inflow (parameterisation) or
71  as fluxes applied at the land-ice ocean interface (\np{ln\_isf}) ;
72\item
73  the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift
74  (\np{nn\_fwb}\forcode{ = 0..2});
75\item
76  the transformation of the solar radiation (if provided as daily mean) into a diurnal cycle
77  (\np{ln\_dm2dc}\forcode{ = .true.});
78  and a neutral drag coefficient can be read from an external wave model (\np{ln\_cdgw}\forcode{ = .true.}).
79\end{itemize}
80The latter option is possible only in case core or mfs bulk formulas are selected.
81
82In this chapter, we first discuss where the surface boundary condition appears in the model equations.
83Then we present the five ways of providing the surface boundary condition,
84followed by the description of the atmospheric pressure and the river runoff.
85Next the scheme for interpolation on the fly is described.
86Finally, the different options that further modify the fluxes applied to the ocean are discussed.
87One of these is modification by icebergs (see \autoref{sec:ICB_icebergs}),
88which act as drifting sources of fresh water.
89Another example of modification is that due to the ice shelf melting/freezing (see \autoref{sec:SBC_isf}),
90which provides additional sources of fresh water.
91
92
93% ================================================================
94% Surface boundary condition for the ocean
95% ================================================================
96\section{Surface boundary condition for the ocean}
97\label{sec:SBC_general}
98
99The surface ocean stress is the stress exerted by the wind and the sea-ice on the ocean.
100It is applied in \mdl{dynzdf} module as a surface boundary condition of the computation of
101the momentum vertical mixing trend (see \autoref{eq:dynzdf_sbc} in \autoref{sec:DYN_zdf}).
102As such, it has to be provided as a 2D vector interpolated onto the horizontal velocity ocean mesh,
103$i.e.$ resolved onto the model (\textbf{i},\textbf{j}) direction at $u$- and $v$-points.
104
105The surface heat flux is decomposed into two parts, a non solar and a solar heat flux,
106$Q_{ns}$ and $Q_{sr}$, respectively.
107The former is the non penetrative part of the heat flux
108($i.e.$ the sum of sensible, latent and long wave heat fluxes plus
109the heat content of the mass exchange with the atmosphere and sea-ice).
110It is applied in \mdl{trasbc} module as a surface boundary condition trend of
111the first level temperature time evolution equation
112(see \autoref{eq:tra_sbc} and \autoref{eq:tra_sbc_lin} in \autoref{subsec:TRA_sbc}).
113The latter is the penetrative part of the heat flux.
114It is applied as a 3D trends of the temperature equation (\mdl{traqsr} module) when
115\np{ln\_traqsr}\forcode{ = .true.}.
116The way the light penetrates inside the water column is generally a sum of decreasing exponentials
117(see \autoref{subsec:TRA_qsr}).
118
119The surface freshwater budget is provided by the \textit{emp} field.
120It represents the mass flux exchanged with the atmosphere (evaporation minus precipitation) and
121possibly with the sea-ice and ice shelves (freezing minus melting of ice).
122It affects both the ocean in two different ways:
123$(i)$  it changes the volume of the ocean and therefore appears in the sea surface height equation as
124a volume flux, and
125$(ii)$ it changes the surface temperature and salinity through the heat and salt contents of
126the mass exchanged with the atmosphere, the sea-ice and the ice shelves.
127
128
129%\colorbox{yellow}{Miss: }
130%
131%A extensive description of all namsbc namelist (parameter that have to be
132%created!)
133%
134%Especially the \np{nn\_fsbc}, the \mdl{sbc\_oce} module (fluxes + mean sst sss ssu
135%ssv) i.e. information required by flux computation or sea-ice
136%
137%\mdl{sbc\_oce} containt the definition in memory of the 7 fields (6+runoff), add
138%a word on runoff: included in surface bc or add as lateral obc{\ldots}.
139%
140%Sbcmod manage the ``providing'' (fourniture) to the ocean the 7 fields
141%
142%Fluxes update only each nf{\_}sbc time step (namsbc) explain relation
143%between nf{\_}sbc and nf{\_}ice, do we define nf{\_}blk??? ? only one
144%nf{\_}sbc
145%
146%Explain here all the namlist namsbc variable{\ldots}.
147%
148% explain : use or not of surface currents
149%
150%\colorbox{yellow}{End Miss }
151
152The ocean model provides, at each time step, to the surface module (\mdl{sbcmod})
153the surface currents, temperature and salinity. 
154These variables are averaged over \np{nn\_fsbc} time-step (\autoref{tab:ssm}), and
155it is these averaged fields which are used to computes the surface fluxes at a frequency of \np{nn\_fsbc} time-step.
156
157
158%-------------------------------------------------TABLE---------------------------------------------------
159\begin{table}[tb]   \begin{center}   \begin{tabular}{|l|l|l|l|}
160\hline
161Variable description             & Model variable  & Units  & point \\  \hline
162i-component of the surface current  & ssu\_m & $m.s^{-1}$   & U \\   \hline
163j-component of the surface current  & ssv\_m & $m.s^{-1}$   & V \\   \hline
164Sea surface temperature          & sst\_m & \r{}$K$      & T \\   \hline
165Sea surface salinty              & sss\_m & $psu$        & T \\   \hline
166\end{tabular}
167\caption{  \protect\label{tab:ssm}
168  Ocean variables provided by the ocean to the surface module (SBC).
169  The variable are averaged over nn{\_}fsbc time step,
170  $i.e.$ the frequency of computation of surface fluxes.}
171\end{center}   \end{table}
172%--------------------------------------------------------------------------------------------------------------
173
174%\colorbox{yellow}{Penser a} mettre dans le restant l'info nn{\_}fsbc ET nn{\_}fsbc*rdt de sorte de reinitialiser la moyenne si on change la frequence ou le pdt
175
176
177% ================================================================
178%       Input Data
179% ================================================================
180\section{Input data generic interface}
181\label{sec:SBC_input}
182
183A generic interface has been introduced to manage the way input data
184(2D or 3D fields, like surface forcing or ocean T and S) are specify in \NEMO.
185This task is archieved by \mdl{fldread}.
186The module was design with four main objectives in mind:
187\begin{enumerate}
188\item
189  optionally provide a time interpolation of the input data at model time-step, whatever their input frequency is,
190  and according to the different calendars available in the model.
191\item
192  optionally provide an on-the-fly space interpolation from the native input data grid to the model grid.
193\item
194  make the run duration independent from the period cover by the input files.
195\item
196  provide a simple user interface and a rather simple developer interface by
197  limiting the number of prerequisite information.
198\end{enumerate} 
199
200As a results the user have only to fill in for each variable a structure in the namelist file to
201define the input data file and variable names, the frequency of the data (in hours or months),
202whether its is climatological data or not, the period covered by the input file (one year, month, week or day),
203and three additional parameters for on-the-fly interpolation.
204When adding a new input variable, the developer has to add the associated structure in the namelist,
205read this information by mirroring the namelist read in \rou{sbc\_blk\_init} for example,
206and simply call \rou{fld\_read} to obtain the desired input field at the model time-step and grid points.
207
208The only constraints are that the input file is a NetCDF file, the file name follows a nomenclature
209(see \autoref{subsec:SBC_fldread}), the period it cover is one year, month, week or day, and,
210if on-the-fly interpolation is used, a file of weights must be supplied (see \autoref{subsec:SBC_iof}).
211
212Note that when an input data is archived on a disc which is accessible directly from the workspace where
213the code is executed, then the use can set the \np{cn\_dir} to the pathway leading to the data.
214By default, the data are assumed to have been copied so that cn\_dir='./'.
215
216% -------------------------------------------------------------------------------------------------------------
217% Input Data specification (\mdl{fldread})
218% -------------------------------------------------------------------------------------------------------------
219\subsection{Input data specification (\protect\mdl{fldread})}
220\label{subsec:SBC_fldread}
221
222The structure associated with an input variable contains the following information:
223\begin{forlines}
224!  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
225!             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
226\end{forlines}
227where
228\begin{description} 
229\item[File name]:
230  the stem name of the NetCDF file to be open.
231  This stem will be completed automatically by the model, with the addition of a '.nc' at its end and
232  by date information and possibly a prefix (when using AGRIF).
233  Tab.\autoref{tab:fldread} provides the resulting file name in all possible cases according to
234  whether it is a climatological file or not, and to the open/close frequency (see below for definition).
235
236%--------------------------------------------------TABLE--------------------------------------------------
237\begin{table}[htbp]
238\begin{center}
239\begin{tabular}{|l|c|c|c|}
240\hline
241                         & daily or weekLLL          & monthly                   &   yearly          \\   \hline
242\np{clim}\forcode{ = .false.} & fn\_yYYYYmMMdDD.nc  &   fn\_yYYYYmMM.nc   &   fn\_yYYYY.nc  \\   \hline
243\np{clim}\forcode{ = .true.}        & not possible                  &  fn\_m??.nc             &   fn                \\   \hline
244\end{tabular}
245\end{center}
246\caption{ \protect\label{tab:fldread}
247  naming nomenclature for climatological or interannual input file, as a function of the Open/close frequency.
248  The stem name is assumed to be 'fn'.
249  For weekly files, the 'LLL' corresponds to the first three letters of the first day of the week
250  ($i.e.$ 'sun','sat','fri','thu','wed','tue','mon').
251  The 'YYYY', 'MM' and 'DD' should be replaced by the actual year/month/day, always coded with 4 or 2 digits.
252  Note that (1) in mpp, if the file is split over each subdomain, the suffix '.nc' is replaced by '\_PPPP.nc',
253  where 'PPPP' is the process number coded with 4 digits;
254  (2) when using AGRIF, the prefix '\_N' is added to files, where 'N'  is the child grid number.}
255\end{table}
256%--------------------------------------------------------------------------------------------------------------
257 
258
259\item[Record frequency]:
260  the frequency of the records contained in the input file.
261  Its unit is in hours if it is positive (for example 24 for daily forcing) or in months if negative
262  (for example -1 for monthly forcing or -12 for annual forcing).
263  Note that this frequency must really be an integer and not a real.
264  On some computers, seting it to '24.' can be interpreted as 240!
265
266\item[Variable name]:
267  the name of the variable to be read in the input NetCDF file.
268
269\item[Time interpolation]:
270  a logical to activate, or not, the time interpolation.
271  If set to 'false', the forcing will have a steplike shape remaining constant during each forcing period.
272  For example, when using a daily forcing without time interpolation, the forcing remaining constant from
273  00h00'00'' to 23h59'59".
274  If set to 'true', the forcing will have a broken line shape.
275  Records are assumed to be dated the middle of the forcing period.
276  For example, when using a daily forcing with time interpolation,
277  linear interpolation will be performed between mid-day of two consecutive days.
278
279\item[Climatological forcing]:
280  a logical to specify if a input file contains climatological forcing which can be cycle in time,
281  or an interannual forcing which will requires additional files if
282  the period covered by the simulation exceed the one of the file.
283  See the above the file naming strategy which impacts the expected name of the file to be opened.
284
285\item[Open/close frequency]:
286  the frequency at which forcing files must be opened/closed.
287  Four cases are coded:
288  'daily', 'weekLLL' (with 'LLL' the first 3 letters of the first day of the week), 'monthly' and 'yearly' which
289  means the forcing files will contain data for one day, one week, one month or one year.
290  Files are assumed to contain data from the beginning of the open/close period.
291  For example, the first record of a yearly file containing daily data is Jan 1st even if
292  the experiment is not starting at the beginning of the year.
293
294\item[Others]:
295  'weights filename', 'pairing rotation' and 'land/sea mask' are associated with
296  on-the-fly interpolation which is described in \autoref{subsec:SBC_iof}.
297
298\end{description}
299
300Additional remarks:\\
301(1) The time interpolation is a simple linear interpolation between two consecutive records of the input data.
302The only tricky point is therefore to specify the date at which we need to do the interpolation and
303the date of the records read in the input files.
304Following \citet{Leclair_Madec_OM09}, the date of a time step is set at the middle of the time step.
305For example, for an experiment starting at 0h00'00" with a one hour time-step,
306a time interpolation will be performed at the following time: 0h30'00", 1h30'00", 2h30'00", etc.
307However, for forcing data related to the surface module,
308values are not needed at every time-step but at every \np{nn\_fsbc} time-step.
309For example with \np{nn\_fsbc}\forcode{ = 3}, the surface module will be called at time-steps 1, 4, 7, etc.
310The date used for the time interpolation is thus redefined to be at the middle of \np{nn\_fsbc} time-step period.
311In the previous example, this leads to: 1h30'00", 4h30'00", 7h30'00", etc. \\ 
312(2) For code readablility and maintenance issues, we don't take into account the NetCDF input file calendar.
313The calendar associated with the forcing field is build according to the information provided by
314user in the record frequency, the open/close frequency and the type of temporal interpolation.
315For example, the first record of a yearly file containing daily data that will be interpolated in time is assumed to
316be start Jan 1st at 12h00'00" and end Dec 31st at 12h00'00". \\
317(3) If a time interpolation is requested, the code will pick up the needed data in the previous (next) file when
318interpolating data with the first (last) record of the open/close period.
319For example, if the input file specifications are ''yearly, containing daily data to be interpolated in time'',
320the values given by the code between 00h00'00" and 11h59'59" on Jan 1st will be interpolated values between
321Dec 31st 12h00'00" and Jan 1st 12h00'00".
322If the forcing is climatological, Dec and Jan will be keep-up from the same year.
323However, if the forcing is not climatological, at the end of
324the open/close period the code will automatically close the current file and open the next one.
325Note that, if the experiment is starting (ending) at the beginning (end) of
326an open/close period we do accept that the previous (next) file is not existing.
327In this case, the time interpolation will be performed between two identical values.
328For example, when starting an experiment on Jan 1st of year Y with yearly files and daily data to be interpolated,
329we do accept that the file related to year Y-1 is not existing.
330The value of Jan 1st will be used as the missing one for Dec 31st of year Y-1.
331If the file of year Y-1 exists, the code will read its last record.
332Therefore, this file can contain only one record corresponding to Dec 31st,
333a useful feature for user considering that it is too heavy to manipulate the complete file for year Y-1.
334
335
336% -------------------------------------------------------------------------------------------------------------
337% Interpolation on the Fly
338% -------------------------------------------------------------------------------------------------------------
339\subsection{Interpolation on-the-fly}
340\label{subsec:SBC_iof}
341
342Interpolation on the Fly allows the user to supply input files required for the surface forcing on
343grids other than the model grid.
344To do this he or she must supply, in addition to the source data file, a file of weights to be used to
345interpolate from the data grid to the model grid.
346The original development of this code used the SCRIP package
347(freely available \href{http://climate.lanl.gov/Software/SCRIP}{here} under a copyright agreement).
348In principle, any package can be used to generate the weights, but the variables in
349the input weights file must have the same names and meanings as assumed by the model.
350Two methods are currently available: bilinear and bicubic interpolation.
351Prior to the interpolation, providing a land/sea mask file, the user can decide to remove land points from
352the input file and substitute the corresponding values with the average of the 8 neighbouring points in
353the native external grid.
354Only "sea points" are considered for the averaging.
355The land/sea mask file must be provided in the structure associated with the input variable.
356The netcdf land/sea mask variable name must be 'LSM' it must have the same horizontal and vertical dimensions of
357the associated variable and should be equal to 1 over land and 0 elsewhere.
358The procedure can be recursively applied setting nn\_lsm > 1 in namsbc namelist.
359Note that nn\_lsm=0 forces the code to not apply the procedure even if a file for land/sea mask is supplied.
360
361\subsubsection{Bilinear interpolation}
362\label{subsec:SBC_iof_bilinear}
363
364The input weights file in this case has two sets of variables:
365src01, src02, src03, src04 and wgt01, wgt02, wgt03, wgt04.
366The "src" variables correspond to the point in the input grid to which the weight "wgt" is to be applied.
367Each src value is an integer corresponding to the index of a point in the input grid when
368written as a one dimensional array.
369For example, for an input grid of size 5x10, point (3,2) is referenced as point 8, since (2-1)*5+3=8.
370There are four of each variable because bilinear interpolation uses the four points defining
371the grid box containing the point to be interpolated.
372All of these arrays are on the model grid, so that values src01(i,j) and wgt01(i,j) are used to
373generate a value for point (i,j) in the model.
374
375Symbolically, the algorithm used is:
376\begin{equation}
377f_{m}(i,j) = f_{m}(i,j) + \sum_{k=1}^{4} {wgt(k)f(idx(src(k)))}
378\end{equation}
379where function idx() transforms a one dimensional index src(k) into a two dimensional index,
380and wgt(1) corresponds to variable "wgt01" for example.
381
382\subsubsection{Bicubic interpolation}
383\label{subsec:SBC_iof_bicubic}
384
385Again there are two sets of variables: "src" and "wgt".
386But in this case there are 16 of each.
387The symbolic algorithm used to calculate values on the model grid is now:
388
389\begin{equation*} \begin{split}
390f_{m}(i,j) =  f_{m}(i,j) +& \sum_{k=1}^{4} {wgt(k)f(idx(src(k)))}     
391              +   \sum_{k=5}^{8} {wgt(k)\left.\frac{\partial f}{\partial i}\right| _{idx(src(k))} }    \\
392              +& \sum_{k=9}^{12} {wgt(k)\left.\frac{\partial f}{\partial j}\right| _{idx(src(k))} }   
393              +   \sum_{k=13}^{16} {wgt(k)\left.\frac{\partial ^2 f}{\partial i \partial j}\right| _{idx(src(k))} }
394\end{split}
395\end{equation*}
396The gradients here are taken with respect to the horizontal indices and not distances since
397the spatial dependency has been absorbed into the weights.
398
399\subsubsection{Implementation}
400\label{subsec:SBC_iof_imp}
401
402To activate this option, a non-empty string should be supplied in
403the weights filename column of the relevant namelist;
404if this is left as an empty string no action is taken.
405In the model, weights files are read in and stored in a structured type (WGT) in the fldread module,
406as and when they are first required.
407This initialisation procedure determines whether the input data grid should be treated as cyclical or not by
408inspecting a global attribute stored in the weights input file.
409This attribute must be called "ew\_wrap" and be of integer type.
410If it is negative, the input non-model grid is assumed not to be cyclic.
411If zero or greater, then the value represents the number of columns that overlap.
412$E.g.$ if the input grid has columns at longitudes 0, 1, 2, .... , 359, then ew\_wrap should be set to 0;
413if longitudes are 0.5, 2.5, .... , 358.5, 360.5, 362.5, ew\_wrap should be 2.
414If the model does not find attribute ew\_wrap, then a value of -999 is assumed.
415In this case the \rou{fld\_read} routine defaults ew\_wrap to value 0 and
416therefore the grid is assumed to be cyclic with no overlapping columns.
417(In fact this only matters when bicubic interpolation is required.)
418Note that no testing is done to check the validity in the model,
419since there is no way of knowing the name used for the longitude variable,
420so it is up to the user to make sure his or her data is correctly represented.
421
422Next the routine reads in the weights.
423Bicubic interpolation is assumed if it finds a variable with name "src05", otherwise bilinear interpolation is used.
424The WGT structure includes dynamic arrays both for the storage of the weights (on the model grid),
425and when required, for reading in the variable to be interpolated (on the input data grid).
426The size of the input data array is determined by examining the values in the "src" arrays to
427find the minimum and maximum i and j values required.
428Since bicubic interpolation requires the calculation of gradients at each point on the grid,
429the corresponding arrays are dimensioned with a halo of width one grid point all the way around.
430When the array of points from the data file is adjacent to an edge of the data grid,
431the halo is either a copy of the row/column next to it (non-cyclical case),
432or is a copy of one from the first few columns on the opposite side of the grid (cyclical case).
433
434\subsubsection{Limitations}
435\label{subsec:SBC_iof_lim}
436
437\begin{enumerate} 
438\item
439  The case where input data grids are not logically rectangular has not been tested.
440\item
441  This code is not guaranteed to produce positive definite answers from positive definite inputs when
442  a bicubic interpolation method is used.
443\item
444  The cyclic condition is only applied on left and right columns, and not to top and bottom rows.
445\item
446  The gradients across the ends of a cyclical grid assume that the grid spacing between
447  the two columns involved are consistent with the weights used.
448\item
449  Neither interpolation scheme is conservative. (There is a conservative scheme available in SCRIP,
450  but this has not been implemented.)
451\end{enumerate}
452
453\subsubsection{Utilities}
454\label{subsec:SBC_iof_util}
455
456% to be completed
457A set of utilities to create a weights file for a rectilinear input grid is available
458(see the directory NEMOGCM/TOOLS/WEIGHTS).
459
460% -------------------------------------------------------------------------------------------------------------
461% Standalone Surface Boundary Condition Scheme
462% -------------------------------------------------------------------------------------------------------------
463\subsection{Standalone surface boundary condition scheme}
464\label{subsec:SAS_iof}
465
466%---------------------------------------namsbc_ana--------------------------------------------------
467
468\nlst{namsbc_sas}
469%--------------------------------------------------------------------------------------------------------------
470
471In some circumstances it may be useful to avoid calculating the 3D temperature,
472salinity and velocity fields and simply read them in from a previous run or receive them from OASIS. 
473For example:
474
475\begin{itemize}
476\item
477  Multiple runs of the model are required in code development to
478  see the effect of different algorithms in the bulk formulae.
479\item
480  The effect of different parameter sets in the ice model is to be examined.
481\item
482  Development of sea-ice algorithms or parameterizations.
483\item
484  Spinup of the iceberg floats
485\item
486  Ocean/sea-ice simulation with both media running in parallel (\np{ln\_mixcpl}\forcode{ = .true.})
487\end{itemize}
488
489The StandAlone Surface scheme provides this utility.
490Its options are defined through the \ngn{namsbc\_sas} namelist variables.
491A new copy of the model has to be compiled with a configuration based on ORCA2\_SAS\_LIM.
492However no namelist parameters need be changed from the settings of the previous run (except perhaps nn{\_}date0).
493In this configuration, a few routines in the standard model are overriden by new versions.
494Routines replaced are:
495
496\begin{itemize}
497\item
498  \mdl{nemogcm}:
499  This routine initialises the rest of the model and repeatedly calls the stp time stepping routine (\mdl{step}).
500  Since the ocean state is not calculated all associated initialisations have been removed.
501\item
502  \mdl{step}:
503  The main time stepping routine now only needs to call the sbc routine (and a few utility functions).
504\item
505  \mdl{sbcmod}:
506  This has been cut down and now only calculates surface forcing and the ice model required.
507  New surface modules that can function when only the surface level of the ocean state is defined can also be added
508  (e.g. icebergs).
509\item
510  \mdl{daymod}:
511  No ocean restarts are read or written (though the ice model restarts are retained),
512  so calls to restart functions have been removed.
513  This also means that the calendar cannot be controlled by time in a restart file,
514  so the user must make sure that nn{\_}date0 in the model namelist is correct for his or her purposes.
515\item
516  \mdl{stpctl}:
517  Since there is no free surface solver, references to it have been removed from \rou{stp\_ctl} module.
518\item
519  \mdl{diawri}:
520  All 3D data have been removed from the output.
521  The surface temperature, salinity and velocity components (which have been read in) are written along with
522  relevant forcing and ice data.
523\end{itemize}
524
525One new routine has been added:
526
527\begin{itemize}
528\item
529  \mdl{sbcsas}:
530  This module initialises the input files needed for reading temperature, salinity and
531  velocity arrays at the surface.
532  These filenames are supplied in namelist namsbc{\_}sas.
533  Unfortunately because of limitations with the \mdl{iom} module,
534  the full 3D fields from the mean files have to be read in and interpolated in time,
535  before using just the top level.
536  Since fldread is used to read in the data, Interpolation on the Fly may be used to change input data resolution.
537\end{itemize}
538
539
540% Missing the description of the 2 following variables:
541%   ln_3d_uve   = .true.    !  specify whether we are supplying a 3D u,v and e3 field
542%   ln_read_frq = .false.    !  specify whether we must read frq or not
543
544
545
546% ================================================================
547% Analytical formulation (sbcana module)
548% ================================================================
549\section{Analytical formulation (\protect\mdl{sbcana})}
550\label{sec:SBC_ana}
551
552%---------------------------------------namsbc_ana--------------------------------------------------
553%
554%\nlst{namsbc_ana}
555%--------------------------------------------------------------------------------------------------------------
556
557The analytical formulation of the surface boundary condition is the default scheme.
558In this case, all the six fluxes needed by the ocean are assumed to be uniform in space.
559They take constant values given in the namelist \ngn{namsbc{\_}ana} by
560the variables \np{rn\_utau0}, \np{rn\_vtau0}, \np{rn\_qns0}, \np{rn\_qsr0}, and \np{rn\_emp0}
561($\textit{emp}=\textit{emp}_S$).
562The runoff is set to zero.
563In addition, the wind is allowed to reach its nominal value within a given number of time steps (\np{nn\_tau000}).
564
565If a user wants to apply a different analytical forcing,
566the \mdl{sbcana} module can be modified to use another scheme.
567As an example, the \mdl{sbc\_ana\_gyre} routine provides the analytical forcing for the GYRE configuration
568(see GYRE configuration manual, in preparation).
569
570
571% ================================================================
572% Flux formulation
573% ================================================================
574\section{Flux formulation (\protect\mdl{sbcflx})}
575\label{sec:SBC_flx}
576%------------------------------------------namsbc_flx----------------------------------------------------
577
578\nlst{namsbc_flx} 
579%-------------------------------------------------------------------------------------------------------------
580
581In the flux formulation (\np{ln\_flx}\forcode{ = .true.}),
582the surface boundary condition fields are directly read from input files.
583The user has to define in the namelist \ngn{namsbc{\_}flx} the name of the file,
584the name of the variable read in the file, the time frequency at which it is given (in hours),
585and a logical setting whether a time interpolation to the model time step is required for this field.
586See \autoref{subsec:SBC_fldread} for a more detailed description of the parameters.
587
588Note that in general, a flux formulation is used in associated with a restoring term to observed SST and/or SSS.
589See \autoref{subsec:SBC_ssr} for its specification.
590
591
592% ================================================================
593% Bulk formulation
594% ================================================================
595\section[Bulk formulation {(\textit{sbcblk\{\_core,\_clio,\_mfs\}.F90})}]
596         {Bulk formulation {(\protect\mdl{sbcblk\_core}, \protect\mdl{sbcblk\_clio}, \protect\mdl{sbcblk\_mfs})}}
597\label{sec:SBC_blk}
598
599In the bulk formulation, the surface boundary condition fields are computed using bulk formulae and atmospheric fields and ocean (and ice) variables.
600
601The atmospheric fields used depend on the bulk formulae used.
602Three bulk formulations are available:
603the CORE, the CLIO and the MFS bulk formulea.
604The choice is made by setting to true one of the following namelist variable:
605\np{ln\_core} ; \np{ln\_clio} or  \np{ln\_mfs}.
606
607Note:
608in forced mode, when a sea-ice model is used, a bulk formulation (CLIO or CORE) have to be used.
609Therefore the two bulk (CLIO and CORE) formulea include the computation of the fluxes over
610both an ocean and an ice surface.
611
612% -------------------------------------------------------------------------------------------------------------
613%        CORE Bulk formulea
614% -------------------------------------------------------------------------------------------------------------
615\subsection{CORE formulea (\protect\mdl{sbcblk\_core}, \protect\np{ln\_core}\forcode{ = .true.})}
616\label{subsec:SBC_blk_core}
617%------------------------------------------namsbc_core----------------------------------------------------
618%
619%\nlst{namsbc_core}
620%-------------------------------------------------------------------------------------------------------------
621
622The CORE bulk formulae have been developed by \citet{Large_Yeager_Rep04}.
623They have been designed to handle the CORE forcing, a mixture of NCEP reanalysis and satellite data.
624They use an inertial dissipative method to compute the turbulent transfer coefficients
625(momentum, sensible heat and evaporation) from the 10 metre wind speed, air temperature and specific humidity.
626This \citet{Large_Yeager_Rep04} dataset is available through
627the \href{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}{GFDL web site}.
628
629Note that substituting ERA40 to NCEP reanalysis fields does not require changes in the bulk formulea themself.
630This is the so-called DRAKKAR Forcing Set (DFS) \citep{Brodeau_al_OM09}.
631
632Options are defined through the  \ngn{namsbc\_core} namelist variables.
633The required 8 input fields are:
634
635%--------------------------------------------------TABLE--------------------------------------------------
636\begin{table}[htbp]   \label{tab:CORE}
637\begin{center}
638\begin{tabular}{|l|c|c|c|}
639\hline
640Variable desciption              & Model variable  & Units   & point \\    \hline
641i-component of the 10m air velocity & utau      & $m.s^{-1}$         & T  \\  \hline
642j-component of the 10m air velocity & vtau      & $m.s^{-1}$         & T  \\  \hline
64310m air temperature              & tair      & \r{}$K$            & T   \\ \hline
644Specific humidity             & humi      & \%              & T \\      \hline
645Incoming long wave radiation     & qlw    & $W.m^{-2}$         & T \\      \hline
646Incoming short wave radiation    & qsr    & $W.m^{-2}$         & T \\      \hline
647Total precipitation (liquid + solid)   & precip & $Kg.m^{-2}.s^{-1}$ & T \\   \hline
648Solid precipitation              & snow      & $Kg.m^{-2}.s^{-1}$ & T \\   \hline
649\end{tabular}
650\end{center}
651\end{table}
652%--------------------------------------------------------------------------------------------------------------
653
654Note that the air velocity is provided at a tracer ocean point, not at a velocity ocean point ($u$- and $v$-points).
655It is simpler and faster (less fields to be read), but it is not the recommended method when
656the ocean grid size is the same or larger than the one of the input atmospheric fields.
657
658The \np{sn\_wndi}, \np{sn\_wndj}, \np{sn\_qsr}, \np{sn\_qlw}, \np{sn\_tair}, \np{sn\_humi}, \np{sn\_prec},
659\np{sn\_snow}, \np{sn\_tdif} parameters describe the fields and the way they have to be used
660(spatial and temporal interpolations).
661
662\np{cn\_dir} is the directory of location of bulk files
663\np{ln\_taudif} is the flag to specify if we use Hight Frequency (HF) tau information (.true.) or not (.false.)
664\np{rn\_zqt}: is the height of humidity and temperature measurements (m)
665\np{rn\_zu}: is the height of wind measurements (m)
666
667Three multiplicative factors are availables:
668\np{rn\_pfac} and \np{rn\_efac} allows to adjust (if necessary) the global freshwater budget by
669increasing/reducing the precipitations (total and snow) and or evaporation, respectively.
670The third one,\np{rn\_vfac}, control to which extend the ice/ocean velocities are taken into account in
671the calculation of surface wind stress.
672Its range should be between zero and one, and it is recommended to set it to 0.
673
674% -------------------------------------------------------------------------------------------------------------
675%        CLIO Bulk formulea
676% -------------------------------------------------------------------------------------------------------------
677\subsection{CLIO formulea (\protect\mdl{sbcblk\_clio}, \protect\np{ln\_clio}\forcode{ = .true.})}
678\label{subsec:SBC_blk_clio}
679%------------------------------------------namsbc_clio----------------------------------------------------
680%
681%\nlst{namsbc_clio}
682%-------------------------------------------------------------------------------------------------------------
683
684The CLIO bulk formulae were developed several years ago for the Louvain-la-neuve coupled ice-ocean model
685(CLIO, \cite{Goosse_al_JGR99}).
686They are simpler bulk formulae.
687They assume the stress to be known and compute the radiative fluxes from a climatological cloud cover.
688
689Options are defined through the  \ngn{namsbc\_clio} namelist variables.
690The required 7 input fields are:
691
692%--------------------------------------------------TABLE--------------------------------------------------
693\begin{table}[htbp]   \label{tab:CLIO}
694\begin{center}
695\begin{tabular}{|l|l|l|l|}
696\hline
697Variable desciption           & Model variable  & Units           & point \\  \hline
698i-component of the ocean stress     & utau         & $N.m^{-2}$         & U \\   \hline
699j-component of the ocean stress     & vtau         & $N.m^{-2}$         & V \\   \hline
700Wind speed module             & vatm         & $m.s^{-1}$         & T \\   \hline
70110m air temperature              & tair         & \r{}$K$            & T \\   \hline
702Specific humidity                & humi         & \%              & T \\   \hline
703Cloud cover                   &           & \%              & T \\   \hline
704Total precipitation (liquid + solid)   & precip    & $Kg.m^{-2}.s^{-1}$ & T \\   \hline
705Solid precipitation              & snow         & $Kg.m^{-2}.s^{-1}$ & T \\   \hline
706\end{tabular}
707\end{center}
708\end{table}
709%--------------------------------------------------------------------------------------------------------------
710
711As for the flux formulation, information about the input data required by the model is provided in
712the namsbc\_blk\_core or namsbc\_blk\_clio namelist (see \autoref{subsec:SBC_fldread}).
713
714% -------------------------------------------------------------------------------------------------------------
715%        MFS Bulk formulae
716% -------------------------------------------------------------------------------------------------------------
717\subsection{MFS formulea (\protect\mdl{sbcblk\_mfs}, \protect\np{ln\_mfs}\forcode{ = .true.})}
718\label{subsec:SBC_blk_mfs}
719%------------------------------------------namsbc_mfs----------------------------------------------------
720%
721%\nlst{namsbc_mfs}
722%----------------------------------------------------------------------------------------------------------
723
724The MFS (Mediterranean Forecasting System) bulk formulae have been developed by \citet{Castellari_al_JMS1998}.
725They have been designed to handle the ECMWF operational data and are currently in use in
726the MFS operational system \citep{Tonani_al_OS08}, \citep{Oddo_al_OS09}.
727The wind stress computation uses a drag coefficient computed according to \citet{Hellerman_Rosenstein_JPO83}.
728The surface boundary condition for temperature involves the balance between
729surface solar radiation, net long-wave radiation, the latent and sensible heat fluxes.
730Solar radiation is dependent on cloud cover and is computed by means of an astronomical formula \citep{Reed_JPO77}.
731Albedo monthly values are from \citet{Payne_JAS72} as means of the values at $40^{o}N$ and $30^{o}N$ for
732the Atlantic Ocean (hence the same latitudinal band of the Mediterranean Sea).
733The net long-wave radiation flux \citep{Bignami_al_JGR95} is a function of
734air temperature, sea-surface temperature, cloud cover and relative humidity.
735Sensible heat and latent heat fluxes are computed by classical bulk formulae parameterised according to
736\citet{Kondo1975}.
737Details on the bulk formulae used can be found in \citet{Maggiore_al_PCE98} and \citet{Castellari_al_JMS1998}.
738
739Options are defined through the \ngn{namsbc\_mfs} namelist variables.
740The required 7 input fields must be provided on the model Grid-T and are:
741\begin{itemize}
742\item          Zonal Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windi})
743\item          Meridional Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windj})
744\item          Total Claud Cover (\%)  (\np{sn\_clc})
745\item          2m Air Temperature ($K$) (\np{sn\_tair})
746\item          2m Dew Point Temperature ($K$)  (\np{sn\_rhm})
747\item          Total Precipitation ${Kg} m^{-2} s^{-1}$ (\np{sn\_prec})
748\item          Mean Sea Level Pressure (${Pa}$) (\np{sn\_msl})
749\end{itemize}
750% -------------------------------------------------------------------------------------------------------------
751% ================================================================
752% Coupled formulation
753% ================================================================
754\section{Coupled formulation (\protect\mdl{sbccpl})}
755\label{sec:SBC_cpl}
756%------------------------------------------namsbc_cpl----------------------------------------------------
757
758\nlst{namsbc_cpl} 
759%-------------------------------------------------------------------------------------------------------------
760
761In the coupled formulation of the surface boundary condition,
762the fluxes are provided by the OASIS coupler at a frequency which is defined in the OASIS coupler,
763while sea and ice surface temperature, ocean and ice albedo, and ocean currents are sent to
764the atmospheric component.
765
766A generalised coupled interface has been developed.
767It is currently interfaced with OASIS-3-MCT (\key{oasis3}).
768It has been successfully used to interface \NEMO to most of the European atmospheric GCM
769(ARPEGE, ECHAM, ECMWF, HadAM, HadGAM, LMDz), as well as to \href{http://wrf-model.org/}{WRF}
770(Weather Research and Forecasting Model).
771
772Note that in addition to the setting of \np{ln\_cpl} to true, the \key{coupled} have to be defined.
773The CPP key is mainly used in sea-ice to ensure that the atmospheric fluxes are actually received by
774the ice-ocean system (no calculation of ice sublimation in coupled mode).
775When PISCES biogeochemical model (\key{top} and \key{pisces}) is also used in the coupled system,
776the whole carbon cycle is computed by defining \key{cpl\_carbon\_cycle}.
777In this case, CO$_2$ fluxes will be exchanged between the atmosphere and the ice-ocean system
778(and need to be activated in \ngn{namsbc{\_}cpl} ).
779
780The namelist above allows control of various aspects of the coupling fields (particularly for vectors) and
781now allows for any coupling fields to have multiple sea ice categories (as required by LIM3 and CICE).
782When indicating a multi-category coupling field in namsbc{\_}cpl the number of categories will be determined by
783the number used in the sea ice model.
784In some limited cases it may be possible to specify single category coupling fields even when
785the sea ice model is running with multiple categories -
786in this case the user should examine the code to be sure the assumptions made are satisfactory.
787In cases where this is definitely not possible the model should abort with an error message.
788The new code has been tested using ECHAM with LIM2, and HadGAM3 with CICE but
789although it will compile with \key{lim3} additional minor code changes may be required to run using LIM3.
790
791
792% ================================================================
793%        Atmospheric pressure
794% ================================================================
795\section{Atmospheric pressure (\protect\mdl{sbcapr})}
796\label{sec:SBC_apr}
797%------------------------------------------namsbc_apr----------------------------------------------------
798
799\nlst{namsbc_apr} 
800%-------------------------------------------------------------------------------------------------------------
801
802The optional atmospheric pressure can be used to force ocean and ice dynamics
803(\np{ln\_apr\_dyn}\forcode{ = .true.}, \textit{\ngn{namsbc}} namelist).
804The input atmospheric forcing defined via \np{sn\_apr} structure (\textit{namsbc\_apr} namelist)
805can be interpolated in time to the model time step, and even in space when the interpolation on-the-fly is used.
806When used to force the dynamics, the atmospheric pressure is further transformed into
807an equivalent inverse barometer sea surface height, $\eta_{ib}$, using:
808\begin{equation} \label{eq:SBC_ssh_ib}
809   \eta_{ib} = -  \frac{1}{g\,\rho_o}  \left( P_{atm} - P_o \right)
810\end{equation}
811where $P_{atm}$ is the atmospheric pressure and $P_o$ a reference atmospheric pressure.
812A value of $101,000~N/m^2$ is used unless \np{ln\_ref\_apr} is set to true.
813In this case $P_o$ is set to the value of $P_{atm}$ averaged over the ocean domain,
814$i.e.$ the mean value of $\eta_{ib}$ is kept to zero at all time step.
815
816The gradient of $\eta_{ib}$ is added to the RHS of the ocean momentum equation (see \mdl{dynspg} for the ocean).
817For sea-ice, the sea surface height, $\eta_m$, which is provided to the sea ice model is set to $\eta - \eta_{ib}$
818(see \mdl{sbcssr} module).
819$\eta_{ib}$ can be set in the output.
820This can simplify altimetry data and model comparison as
821inverse barometer sea surface height is usually removed from these date prior to their distribution.
822
823When using time-splitting and BDY package for open boundaries conditions,
824the equivalent inverse barometer sea surface height $\eta_{ib}$ can be added to BDY ssh data:
825\np{ln\_apr\_obc}  might be set to true.
826
827% ================================================================
828%        Surface Tides Forcing
829% ================================================================
830\section{Surface tides (\protect\mdl{sbctide})}
831\label{sec:SBC_tide}
832
833%------------------------------------------nam_tide---------------------------------------
834
835\nlst{nam_tide}
836%-----------------------------------------------------------------------------------------
837
838The tidal forcing, generated by the gravity forces of the Earth-Moon and Earth-Sun sytems,
839is activated if \np{ln\_tide} and \np{ln\_tide\_pot} are both set to \np{.true.} in \ngn{nam\_tide}.
840This translates as an additional barotropic force in the momentum equations \ref{eq:PE_dyn} such that:
841\begin{equation}     \label{eq:PE_dyn_tides}
842\frac{\partial {\rm {\bf U}}_h }{\partial t}= ...
843+g\nabla (\Pi_{eq} + \Pi_{sal})
844\end{equation} 
845where $\Pi_{eq}$ stands for the equilibrium tidal forcing and $\Pi_{sal}$ a self-attraction and loading term (SAL).
846 
847The equilibrium tidal forcing is expressed as a sum over the chosen constituents $l$ in \ngn{nam\_tide}.
848The constituents are defined such that \np{clname(1) = 'M2', clname(2)='S2', etc...}.
849For the three types of tidal frequencies it reads: \\
850Long period tides :
851\begin{equation}
852\Pi_{eq}(l)=A_{l}(1+k-h)(\frac{1}{2}-\frac{3}{2}sin^{2}\phi)cos(\omega_{l}t+V_{l})
853\end{equation}
854diurnal tides :
855\begin{equation}
856\Pi_{eq}(l)=A_{l}(1+k-h)(sin 2\phi)cos(\omega_{l}t+\lambda+V_{l})
857\end{equation}
858Semi-diurnal tides:
859\begin{equation}
860\Pi_{eq}(l)=A_{l}(1+k-h)(cos^{2}\phi)cos(\omega_{l}t+2\lambda+V_{l})
861\end{equation}
862Here $A_{l}$ is the amplitude, $\omega_{l}$ is the frequency, $\phi$ the latitude, $\lambda$ the longitude,
863$V_{0l}$ a phase shift with respect to Greenwich meridian and $t$ the time.
864The Love number factor $(1+k-h)$ is here taken as a constant (0.7).
865
866The SAL term should in principle be computed online as it depends on the model tidal prediction itself
867(see \citet{Arbic2004} for a discussion about the practical implementation of this term).
868Nevertheless, the complex calculations involved would make this computationally too expensive.
869Here, practical solutions are whether to read complex estimates $\Pi_{sal}(l)$ from an external model
870(\np{ln\_read\_load=.true.}) or use a ``scalar approximation'' (\np{ln\_scal\_load=.true.}).
871In the latter case, it reads:\\
872\begin{equation}
873\Pi_{sal} = \beta \eta
874\end{equation}
875where $\beta$ (\np{rn\_scal\_load}, $\approx0.09$) is a spatially constant scalar,
876often chosen to minimize tidal prediction errors.
877Setting both \np{ln\_read\_load} and \np{ln\_scal\_load} to false removes the SAL contribution.
878
879% ================================================================
880%        River runoffs
881% ================================================================
882\section{River runoffs (\protect\mdl{sbcrnf})}
883\label{sec:SBC_rnf}
884%------------------------------------------namsbc_rnf----------------------------------------------------
885
886\nlst{namsbc_rnf} 
887%-------------------------------------------------------------------------------------------------------------
888
889%River runoff generally enters the ocean at a nonzero depth rather than through the surface.
890%Many models, however, have traditionally inserted river runoff to the top model cell.
891%This was the case in \NEMO prior to the version 3.3. The switch toward a input of runoff
892%throughout a nonzero depth has been motivated by the numerical and physical problems
893%that arise when the top grid cells are of the order of one meter. This situation is common in
894%coastal modelling and becomes more and more often open ocean and climate modelling
895%\footnote{At least a top cells thickness of 1~meter and a 3 hours forcing frequency are
896%required to properly represent the diurnal cycle \citep{Bernie_al_JC05}. see also \autoref{fig:SBC_dcy}.}.
897
898
899%To do this we need to treat evaporation/precipitation fluxes and river runoff differently in the
900%\mdl{tra\_sbc} module.  We decided to separate them throughout the code, so that the variable
901%\textit{emp} represented solely evaporation minus precipitation fluxes, and a new 2d variable
902%rnf was added which represents the volume flux of river runoff (in kg/m2s to remain consistent with
903%emp).  This meant many uses of emp and emps needed to be changed, a list of all modules which use
904%emp or emps and the changes made are below:
905
906
907%Rachel:
908River runoff generally enters the ocean at a nonzero depth rather than through the surface.
909Many models, however, have traditionally inserted river runoff to the top model cell.
910This was the case in \NEMO prior to the version 3.3,
911and was combined with an option to increase vertical mixing near the river mouth.
912
913However, with this method numerical and physical problems arise when the top grid cells are of the order of one meter.
914This situation is common in coastal modelling and is becoming more common in open ocean and climate modelling
915\footnote{
916  At least a top cells thickness of 1~meter and a 3 hours forcing frequency are required to
917  properly represent the diurnal cycle \citep{Bernie_al_JC05}.
918  see also \autoref{fig:SBC_dcy}.}.
919
920As such from V~3.3 onwards it is possible to add river runoff through a non-zero depth,
921and for the temperature and salinity of the river to effect the surrounding ocean.
922The user is able to specify, in a NetCDF input file, the temperature and salinity of the river,
923along with the depth (in metres) which the river should be added to.
924
925Namelist variables in \ngn{namsbc\_rnf}, \np{ln\_rnf\_depth}, \np{ln\_rnf\_sal} and
926\np{ln\_rnf\_temp} control whether the river attributes (depth, salinity and temperature) are read in and used.
927If these are set as false the river is added to the surface box only, assumed to be fresh (0~psu),
928and/or taken as surface temperature respectively.
929
930The runoff value and attributes are read in in sbcrnf. 
931For temperature -999 is taken as missing data and the river temperature is taken to
932be the surface temperatue at the river point.
933For the depth parameter a value of -1 means the river is added to the surface box only,
934and a value of -999 means the river is added through the entire water column.
935After being read in the temperature and salinity variables are multiplied by the amount of runoff
936(converted into m/s) to give the heat and salt content of the river runoff.
937After the user specified depth is read ini,
938the number of grid boxes this corresponds to is calculated and stored in the variable \np{nz\_rnf}.
939The variable \textit{h\_dep} is then calculated to be the depth (in metres) of
940the bottom of the lowest box the river water is being added to
941(i.e. the total depth that river water is being added to in the model).
942
943The mass/volume addition due to the river runoff is, at each relevant depth level, added to
944the horizontal divergence (\textit{hdivn}) in the subroutine \rou{sbc\_rnf\_div} (called from \mdl{divcur}).
945This increases the diffusion term in the vicinity of the river, thereby simulating a momentum flux.
946The sea surface height is calculated using the sum of the horizontal divergence terms,
947and so the river runoff indirectly forces an increase in sea surface height.
948
949The \textit{hdivn} terms are used in the tracer advection modules to force vertical velocities.
950This causes a mass of water, equal to the amount of runoff, to be moved into the box above.
951The heat and salt content of the river runoff is not included in this step,
952and so the tracer concentrations are diluted as water of ocean temperature and salinity is moved upward out of
953the box and replaced by the same volume of river water with no corresponding heat and salt addition.
954
955For the linear free surface case, at the surface box the tracer advection causes a flux of water
956(of equal volume to the runoff) through the sea surface out of the domain,
957which causes a salt and heat flux out of the model.
958As such the volume of water does not change, but the water is diluted.
959
960For the non-linear free surface case (\key{vvl}), no flux is allowed through the surface.
961Instead in the surface box (as well as water moving up from the boxes below) a volume of runoff water is added with
962no corresponding heat and salt addition and so as happens in the lower boxes there is a dilution effect.
963(The runoff addition to the top box along with the water being moved up through
964boxes below means the surface box has a large increase in volume, whilst all other boxes remain the same size)
965
966In trasbc the addition of heat and salt due to the river runoff is added.
967This is done in the same way for both vvl and non-vvl.
968The temperature and salinity are increased through the specified depth according to
969the heat and salt content of the river.
970
971In the non-linear free surface case (vvl),
972near the end of the time step the change in sea surface height is redistrubuted through the grid boxes,
973so that the original ratios of grid box heights are restored.
974In doing this water is moved into boxes below, throughout the water column,
975so the large volume addition to the surface box is spread between all the grid boxes.
976
977It is also possible for runnoff to be specified as a negative value for modelling flow through straits,
978i.e. modelling the Baltic flow in and out of the North Sea.
979When the flow is out of the domain there is no change in temperature and salinity,
980regardless of the namelist options used,
981as the ocean water leaving the domain removes heat and salt (at the same concentration) with it.
982
983
984%\colorbox{yellow}{Nevertheless, Pb of vertical resolution and 3D input : increase vertical mixing near river mouths to mimic a 3D river
985
986%All river runoff and emp fluxes are assumed to be fresh water (zero salinity) and at the same temperature as the sea surface.}
987
988%\colorbox{yellow}{river mouths{\ldots}}
989
990%IF( ln_rnf ) THEN                                     ! increase diffusivity at rivers mouths
991%        DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + rn_avt_rnf * rnfmsk(:,:)   ;   END DO
992%ENDIF
993
994%\gmcomment{  word doc of runoffs:
995%
996%In the current \NEMO setup river runoff is added to emp fluxes, these are then applied at just the sea surface as a volume change (in the variable volume case this is a literal volume change, and in the linear free surface case the free surface is moved) and a salt flux due to the concentration/dilution effect.  There is also an option to increase vertical mixing near river mouths; this gives the effect of having a 3d river.  All river runoff and emp fluxes are assumed to be fresh water (zero salinity) and at the same temperature as the sea surface.
997%Our aim was to code the option to specify the temperature and salinity of river runoff, (as well as the amount), along with the depth that the river water will affect.  This would make it possible to model low salinity outflow, such as the Baltic, and would allow the ocean temperature to be affected by river runoff. 
998
999%The depth option makes it possible to have the river water affecting just the surface layer, throughout depth, or some specified point in between.
1000
1001%To do this we need to treat evaporation/precipitation fluxes and river runoff differently in the tra_sbc module.  We decided to separate them throughout the code, so that the variable emp represented solely evaporation minus precipitation fluxes, and a new 2d variable rnf was added which represents the volume flux of river runoff (in kg/m2s to remain consistent with emp).  This meant many uses of emp and emps needed to be changed, a list of all modules which use emp or emps and the changes made are below:
1002
1003%}
1004% ================================================================
1005%        Ice shelf melting
1006% ================================================================
1007\section{Ice shelf melting (\protect\mdl{sbcisf})}
1008\label{sec:SBC_isf}
1009%------------------------------------------namsbc_isf----------------------------------------------------
1010
1011\nlst{namsbc_isf}
1012%--------------------------------------------------------------------------------------------------------
1013Namelist variable in \ngn{namsbc}, \np{nn\_isf}, controls the ice shelf representation used.
1014\begin{description}
1015\item[\np{nn\_isf}\forcode{ = 1}]
1016  The ice shelf cavity is represented (\np{ln\_isfcav}\forcode{ = .true.} needed).
1017  The fwf and heat flux are computed.
1018  Two different bulk formula are available:
1019   \begin{description}
1020   \item[\np{nn\_isfblk}\forcode{ = 1}]
1021     The bulk formula used to compute the melt is based the one described in \citet{Hunter2006}.
1022     This formulation is based on a balance between the upward ocean heat flux and
1023     the latent heat flux at the ice shelf base.
1024   \item[\np{nn\_isfblk}\forcode{ = 2}]
1025     The bulk formula used to compute the melt is based the one described in \citet{Jenkins1991}.
1026     This formulation is based on a 3 equations formulation
1027     (a heat flux budget, a salt flux budget and a linearised freezing point temperature equation).
1028   \end{description}
1029   For this 2 bulk formulations, there are 3 different ways to compute the exchange coeficient:
1030   \begin{description}
1031   \item[\np{nn\_gammablk}\forcode{ = 0}]
1032     The salt and heat exchange coefficients are constant and defined by \np{rn\_gammas0} and \np{rn\_gammat0}
1033   \item[\np{nn\_gammablk}\forcode{ = 1}]
1034     The salt and heat exchange coefficients are velocity dependent and defined as
1035     \np{rn\_gammas0}$ \times u_{*}$ and \np{rn\_gammat0}$ \times u_{*}$ where
1036     $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn\_hisf\_tbl} meters).
1037     See \citet{Jenkins2010} for all the details on this formulation.
1038   \item[\np{nn\_gammablk}\forcode{ = 2}]
1039     The salt and heat exchange coefficients are velocity and stability dependent and defined as
1040     $\gamma_{T,S} = \frac{u_{*}}{\Gamma_{Turb} + \Gamma^{T,S}_{Mole}}$ where
1041     $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn\_hisf\_tbl} meters),
1042     $\Gamma_{Turb}$ the contribution of the ocean stability and
1043     $\Gamma^{T,S}_{Mole}$ the contribution of the molecular diffusion.
1044     See \citet{Holland1999} for all the details on this formulation.
1045   \end{description}
1046 \item[\np{nn\_isf}\forcode{ = 2}]
1047   A parameterisation of isf is used. The ice shelf cavity is not represented.
1048   The fwf is distributed along the ice shelf edge between the depth of the average grounding line (GL)
1049   (\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front
1050   (\np{sn\_depmin\_isf}) as in (\np{nn\_isf}\forcode{ = 3}).
1051   Furthermore the fwf and heat flux are computed using the \citet{Beckmann2003} parameterisation of isf melting.
1052   The effective melting length (\np{sn\_Leff\_isf}) is read from a file.
1053 \item[\np{nn\_isf}\forcode{ = 3}]
1054   A simple parameterisation of isf is used. The ice shelf cavity is not represented.
1055   The fwf (\np{sn\_rnfisf}) is prescribed and distributed along the ice shelf edge between
1056   the depth of the average grounding line (GL) (\np{sn\_depmax\_isf}) and
1057   the base of the ice shelf along the calving front (\np{sn\_depmin\_isf}).
1058   The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$.
1059 \item[\np{nn\_isf}\forcode{ = 4}]
1060   The ice shelf cavity is opened (\np{ln\_isfcav}\forcode{ = .true.} needed).
1061   However, the fwf is not computed but specified from file \np{sn\_fwfisf}).
1062   The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$.\\
1063\end{description}
1064
1065$\bullet$ \np{nn\_isf}\forcode{ = 1} and \np{nn\_isf}\forcode{ = 2} compute a melt rate based on
1066the water mass properties, ocean velocities and depth.
1067This flux is thus highly dependent of the model resolution (horizontal and vertical),
1068realism of the water masses onto the shelf ...\\
1069
1070$\bullet$ \np{nn\_isf}\forcode{ = 3} and \np{nn\_isf}\forcode{ = 4} read the melt rate from a file.
1071You have total control of the fwf forcing.
1072This can be useful if the water masses on the shelf are not realistic or
1073the resolution (horizontal/vertical) are too coarse to have realistic melting or
1074for studies where you need to control your heat and fw input.\\ 
1075
1076A namelist parameters control over how many meters the heat and fw fluxes are spread.
1077\np{rn\_hisf\_tbl}] is the top boundary layer thickness as defined in \citet{Losch2008}.
1078This parameter is only used if \np{nn\_isf}\forcode{ = 1} or \np{nn\_isf}\forcode{ = 4}.
1079
1080If \np{rn\_hisf\_tbl}\forcode{ = 0}., the fluxes are put in the top level whatever is its tickness.
1081
1082If \np{rn\_hisf\_tbl} $>$ 0., the fluxes are spread over the first \np{rn\_hisf\_tbl} m
1083(ie over one or several cells).\\
1084
1085The ice shelf melt is implemented as a volume flux with in the same way as for the runoff.
1086The fw addition due to the ice shelf melting is, at each relevant depth level, added to
1087the horizontal divergence (\textit{hdivn}) in the subroutine \rou{sbc\_isf\_div}, called from \mdl{divcur}.
1088See the runoff section \autoref{sec:SBC_rnf} for all the details about the divergence correction.
1089
1090
1091\section{Ice sheet coupling}
1092\label{sec:SBC_iscpl}
1093%------------------------------------------namsbc_iscpl----------------------------------------------------
1094
1095\nlst{namsbc_iscpl}
1096%--------------------------------------------------------------------------------------------------------
1097Ice sheet/ocean coupling is done through file exchange at the restart step.
1098NEMO, at each restart step, read the bathymetry and ice shelf draft variable in a netcdf file.
1099If \np{ln\_iscpl}\forcode{ = .true.}, the isf draft is assume to be different at each restart step with
1100potentially some new wet/dry cells due to the ice sheet dynamics/thermodynamics.
1101The wetting and drying scheme applied on the restart is very simple and described below for the 6 different cases:
1102\begin{description}
1103\item[Thin a cell down:]
1104  T/S/ssh are unchanged and U/V in the top cell are corrected to keep the barotropic transport (bt) constant
1105  ($bt_b=bt_n$).
1106\item[Enlarge  a cell:]
1107  See case "Thin a cell down"
1108\item[Dry a cell:]
1109  mask, T/S, U/V and ssh are set to 0.
1110  Furthermore, U/V into the water column are modified to satisfy ($bt_b=bt_n$).
1111\item[Wet a cell:] 
1112  mask is set to 1, T/S is extrapolated from neighbours, $ssh_n = ssh_b$ and U/V set to 0.
1113  If no neighbours along i,j and k, T/S/U/V and mask are set to 0.
1114\item[Dry a column:]
1115   mask, T/S, U/V are set to 0 everywhere in the column and ssh set to 0.
1116\item[Wet a column:]
1117  set mask to 1, T/S is extrapolated from neighbours, ssh is extrapolated from neighbours and U/V set to 0.
1118  If no neighbour, T/S/U/V and mask set to 0.
1119\end{description}
1120The extrapolation is call \np{nn\_drown} times.
1121It means that if the grounding line retreat by more than \np{nn\_drown} cells between 2 coupling steps,
1122the code will be unable to fill all the new wet cells properly.
1123The default number is set up for the MISOMIP idealised experiments.\\
1124This coupling procedure is able to take into account grounding line and calving front migration.
1125However, it is a non-conservative processe.
1126This could lead to a trend in heat/salt content and volume.
1127In order to remove the trend and keep the conservation level as close to 0 as possible,
1128a simple conservation scheme is available with \np{ln\_hsb}\forcode{ = .true.}.
1129The heat/salt/vol. gain/loss is diagnose, as well as the location.
1130Based on what is done on sbcrnf to prescribed a source of heat/salt/vol.,
1131the heat/salt/vol. gain/loss is removed/added, over a period of \np{rn\_fiscpl} time step, into the system.
1132So after \np{rn\_fiscpl} time step, all the heat/salt/vol. gain/loss due to extrapolation process is canceled.\\
1133
1134As the before and now fields are not compatible (modification of the geometry),
1135the restart time step is prescribed to be an euler time step instead of a leap frog and $fields_b = fields_n$.
1136%
1137% ================================================================
1138%        Handling of icebergs
1139% ================================================================
1140\section{Handling of icebergs (ICB)}
1141\label{sec:ICB_icebergs}
1142%------------------------------------------namberg----------------------------------------------------
1143
1144\nlst{namberg}
1145%-------------------------------------------------------------------------------------------------------------
1146
1147Icebergs are modelled as lagrangian particles in NEMO \citep{Marsh_GMD2015}.
1148Their physical behaviour is controlled by equations as described in \citet{Martin_Adcroft_OM10} ).
1149(Note that the authors kindly provided a copy of their code to act as a basis for implementation in NEMO).
1150Icebergs are initially spawned into one of ten classes which have specific mass and thickness as
1151described in the \ngn{namberg} namelist: \np{rn\_initial\_mass} and \np{rn\_initial\_thickness}.
1152Each class has an associated scaling (\np{rn\_mass\_scaling}),
1153which is an integer representing how many icebergs of this class are being described as one lagrangian point
1154(this reduces the numerical problem of tracking every single iceberg).
1155They are enabled by setting \np{ln\_icebergs}\forcode{ = .true.}.
1156
1157Two initialisation schemes are possible.
1158\begin{description}
1159\item[\np{nn\_test\_icebergs}~$>$~0]
1160  In this scheme, the value of \np{nn\_test\_icebergs} represents the class of iceberg to generate
1161  (so between 1 and 10), and \np{nn\_test\_icebergs} provides a lon/lat box in the domain at each grid point of
1162  which an iceberg is generated at the beginning of the run.
1163  (Note that this happens each time the timestep equals \np{nn\_nit000}.)
1164  \np{nn\_test\_icebergs} is defined by four numbers in \np{nn\_test\_box} representing the corners of
1165  the geographical box: lonmin,lonmax,latmin,latmax
1166\item[\np{nn\_test\_icebergs}\forcode{ = -1}]
1167  In this scheme the model reads a calving file supplied in the \np{sn\_icb} parameter.
1168  This should be a file with a field on the configuration grid (typically ORCA)
1169  representing ice accumulation rate at each model point.
1170  These should be ocean points adjacent to land where icebergs are known to calve.
1171  Most points in this input grid are going to have value zero.
1172  When the model runs, ice is accumulated at each grid point which has a non-zero source term.
1173  At each time step, a test is performed to see if there is enough ice mass to
1174  calve an iceberg of each class in order (1 to 10).
1175  Note that this is the initial mass multiplied by the number each particle represents ($i.e.$ the scaling).
1176  If there is enough ice, a new iceberg is spawned and the total available ice reduced accordingly.
1177\end{description}
1178
1179Icebergs are influenced by wind, waves and currents, bottom melt and erosion.
1180The latter act to disintegrate the iceberg.
1181This is either all melted freshwater,
1182or (if \np{rn\_bits\_erosion\_fraction}~$>$~0) into melt and additionally small ice bits
1183which are assumed to propagate with their larger parent and thus delay fluxing into the ocean.
1184Melt water (and other variables on the configuration grid) are written into the main NEMO model output files.
1185
1186Extensive diagnostics can be produced.
1187Separate output files are maintained for human-readable iceberg information.
1188A separate file is produced for each processor (independent of \np{ln\_ctl}).
1189The amount of information is controlled by two integer parameters:
1190\begin{description}
1191\item[\np{nn\_verbose\_level}] takes a value between one and four and
1192  represents an increasing number of points in the code at which variables are written,
1193  and an increasing level of obscurity.
1194\item[\np{nn\_verbose\_write}] is the number of timesteps between writes
1195\end{description}
1196
1197Iceberg trajectories can also be written out and this is enabled by setting \np{nn\_sample\_rate}~$>$~0.
1198A non-zero value represents how many timesteps between writes of information into the output file.
1199These output files are in NETCDF format.
1200When \key{mpp\_mpi} is defined, each output file contains only those icebergs in the corresponding processor.
1201Trajectory points are written out in the order of their parent iceberg in the model's "linked list" of icebergs.
1202So care is needed to recreate data for individual icebergs,
1203since its trajectory data may be spread across multiple files.
1204
1205
1206% ================================================================
1207% Miscellanea options
1208% ================================================================
1209\section{Miscellaneous options}
1210\label{sec:SBC_misc}
1211
1212% -------------------------------------------------------------------------------------------------------------
1213%        Diurnal cycle
1214% -------------------------------------------------------------------------------------------------------------
1215\subsection{Diurnal cycle (\protect\mdl{sbcdcy})}
1216\label{subsec:SBC_dcy}
1217%------------------------------------------namsbc_rnf----------------------------------------------------
1218%
1219\nlst{namsbc} 
1220%-------------------------------------------------------------------------------------------------------------
1221
1222%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1223\begin{figure}[!t]    \begin{center}
1224\includegraphics[width=0.8\textwidth]{Fig_SBC_diurnal}
1225\caption{ \protect\label{fig:SBC_diurnal}
1226  Example of recontruction of the diurnal cycle variation of short wave flux from daily mean values.
1227  The reconstructed diurnal cycle (black line) is chosen as
1228  the mean value of the analytical cycle (blue line) over a time step,
1229  not as the mid time step value of the analytically cycle (red square).
1230  From \citet{Bernie_al_CD07}.}
1231\end{center}   \end{figure}
1232%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1233
1234\cite{Bernie_al_JC05} have shown that to capture 90$\%$ of the diurnal variability of SST requires a vertical resolution in upper ocean of 1~m or better and a temporal resolution of the surface fluxes of 3~h or less.
1235Unfortunately high frequency forcing fields are rare, not to say inexistent.
1236Nevertheless, it is possible to obtain a reasonable diurnal cycle of the SST knowning only short wave flux (SWF) at
1237high frequency \citep{Bernie_al_CD07}.
1238Furthermore, only the knowledge of daily mean value of SWF is needed,
1239as higher frequency variations can be reconstructed from them,
1240assuming that the diurnal cycle of SWF is a scaling of the top of the atmosphere diurnal cycle of incident SWF.
1241The \cite{Bernie_al_CD07} reconstruction algorithm is available in \NEMO by
1242setting \np{ln\_dm2dc}\forcode{ = .true.} (a \textit{\ngn{namsbc}} namelist variable) when
1243using CORE bulk formulea (\np{ln\_blk\_core}\forcode{ = .true.}) or
1244the flux formulation (\np{ln\_flx}\forcode{ = .true.}).
1245The reconstruction is performed in the \mdl{sbcdcy} module.
1246The detail of the algoritm used can be found in the appendix~A of \cite{Bernie_al_CD07}.
1247The algorithm preserve the daily mean incoming SWF as the reconstructed SWF at
1248a given time step is the mean value of the analytical cycle over this time step (\autoref{fig:SBC_diurnal}).
1249The use of diurnal cycle reconstruction requires the input SWF to be daily
1250($i.e.$ a frequency of 24 and a time interpolation set to true in \np{sn\_qsr} namelist parameter).
1251Furthermore, it is recommended to have a least 8 surface module time step per day,
1252that is  $\rdt \ nn\_fsbc < 10,800~s = 3~h$.
1253An example of recontructed SWF is given in \autoref{fig:SBC_dcy} for a 12 reconstructed diurnal cycle,
1254one every 2~hours (from 1am to 11pm).
1255
1256%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1257\begin{figure}[!t]  \begin{center}
1258\includegraphics[width=0.7\textwidth]{Fig_SBC_dcy}
1259\caption{ \protect\label{fig:SBC_dcy}
1260  Example of recontruction of the diurnal cycle variation of short wave flux from
1261  daily mean values on an ORCA2 grid with a time sampling of 2~hours (from 1am to 11pm).
1262  The display is on (i,j) plane. }
1263\end{center}   \end{figure}
1264%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1265
1266Note also that the setting a diurnal cycle in SWF is highly recommended when
1267the top layer thickness approach 1~m or less, otherwise large error in SST can appear due to
1268an inconsistency between the scale of the vertical resolution and the forcing acting on that scale.
1269
1270% -------------------------------------------------------------------------------------------------------------
1271%        Rotation of vector pairs onto the model grid directions
1272% -------------------------------------------------------------------------------------------------------------
1273\subsection{Rotation of vector pairs onto the model grid directions}
1274\label{subsec:SBC_rotation}
1275
1276When using a flux (\np{ln\_flx}\forcode{ = .true.}) or
1277bulk (\np{ln\_clio}\forcode{ = .true.} or \np{ln\_core}\forcode{ = .true.}) formulation,
1278pairs of vector components can be rotated from east-north directions onto the local grid directions.
1279This is particularly useful when interpolation on the fly is used since here any vectors are likely to
1280be defined relative to a rectilinear grid.
1281To activate this option a non-empty string is supplied in the rotation pair column of the relevant namelist.
1282The eastward component must start with "U" and the northward component with "V". 
1283The remaining characters in the strings are used to identify which pair of components go together.
1284So for example, strings "U1" and "V1" next to "utau" and "vtau" would pair the wind stress components together and
1285rotate them on to the model grid directions;
1286"U2" and "V2" could be used against a second pair of components, and so on.
1287The extra characters used in the strings are arbitrary.
1288The rot\_rep routine from the \mdl{geo2ocean} module is used to perform the rotation.
1289
1290% -------------------------------------------------------------------------------------------------------------
1291%        Surface restoring to observed SST and/or SSS
1292% -------------------------------------------------------------------------------------------------------------
1293\subsection{Surface restoring to observed SST and/or SSS (\protect\mdl{sbcssr})}
1294\label{subsec:SBC_ssr}
1295%------------------------------------------namsbc_ssr----------------------------------------------------
1296
1297\nlst{namsbc_ssr} 
1298%-------------------------------------------------------------------------------------------------------------
1299
1300IOptions are defined through the \ngn{namsbc\_ssr} namelist variables.
1301On forced mode using a flux formulation (\np{ln\_flx}\forcode{ = .true.}),
1302a feedback term \emph{must} be added to the surface heat flux $Q_{ns}^o$:
1303\begin{equation} \label{eq:sbc_dmp_q}
1304Q_{ns} = Q_{ns}^o + \frac{dQ}{dT} \left( \left. T \right|_{k=1} - SST_{Obs} \right)
1305\end{equation}
1306where SST is a sea surface temperature field (observed or climatological),
1307$T$ is the model surface layer temperature and
1308$\frac{dQ}{dT}$ is a negative feedback coefficient usually taken equal to $-40~W/m^2/K$.
1309For a $50~m$ mixed-layer depth, this value corresponds to a relaxation time scale of two months.
1310This term ensures that if $T$ perfectly matches the supplied SST, then $Q$ is equal to $Q_o$.
1311
1312In the fresh water budget, a feedback term can also be added.
1313Converted into an equivalent freshwater flux, it takes the following expression :
1314
1315\begin{equation} \label{eq:sbc_dmp_emp}
1316\textit{emp} = \textit{emp}_o + \gamma_s^{-1} e_{3t}  \frac{  \left(\left.S\right|_{k=1}-SSS_{Obs}\right)}
1317                                             {\left.S\right|_{k=1}}
1318\end{equation}
1319
1320where $\textit{emp}_{o }$ is a net surface fresh water flux
1321(observed, climatological or an atmospheric model product),
1322\textit{SSS}$_{Obs}$ is a sea surface salinity
1323(usually a time interpolation of the monthly mean Polar Hydrographic Climatology \citep{Steele2001}),
1324$\left.S\right|_{k=1}$ is the model surface layer salinity and
1325$\gamma_s$ is a negative feedback coefficient which is provided as a namelist parameter.
1326Unlike heat flux, there is no physical justification for the feedback term in \autoref{eq:sbc_dmp_emp} as
1327the atmosphere does not care about ocean surface salinity \citep{Madec1997}.
1328The SSS restoring term should be viewed as a flux correction on freshwater fluxes to
1329reduce the uncertainties we have on the observed freshwater budget.
1330
1331% -------------------------------------------------------------------------------------------------------------
1332%        Handling of ice-covered area
1333% -------------------------------------------------------------------------------------------------------------
1334\subsection{Handling of ice-covered area  (\textit{sbcice\_...})}
1335\label{subsec:SBC_ice-cover}
1336
1337The presence at the sea surface of an ice covered area modifies all the fluxes transmitted to the ocean.
1338There are several way to handle sea-ice in the system depending on
1339the value of the \np{nn\_ice} namelist parameter found in \ngn{namsbc} namelist.
1340\begin{description}
1341\item[nn{\_}ice = 0]
1342  there will never be sea-ice in the computational domain.
1343  This is a typical namelist value used for tropical ocean domain.
1344  The surface fluxes are simply specified for an ice-free ocean.
1345  No specific things is done for sea-ice.
1346\item[nn{\_}ice = 1]
1347  sea-ice can exist in the computational domain, but no sea-ice model is used.
1348  An observed ice covered area is read in a file.
1349  Below this area, the SST is restored to the freezing point and
1350  the heat fluxes are set to $-4~W/m^2$ ($-2~W/m^2$) in the northern (southern) hemisphere.
1351  The associated modification of the freshwater fluxes are done in such a way that
1352  the change in buoyancy fluxes remains zero.
1353  This prevents deep convection to occur when trying to reach the freezing point
1354  (and so ice covered area condition) while the SSS is too large.
1355  This manner of managing sea-ice area, just by using si IF case,
1356  is usually referred as the \textit{ice-if} model.
1357  It can be found in the \mdl{sbcice{\_}if} module.
1358\item[nn{\_}ice = 2 or more]
1359  A full sea ice model is used.
1360  This model computes the ice-ocean fluxes,
1361  that are combined with the air-sea fluxes using the ice fraction of each model cell to
1362  provide the surface ocean fluxes.
1363  Note that the activation of a sea-ice model is is done by defining a CPP key (\key{lim3} or \key{cice}).
1364  The activation automatically overwrites the read value of nn{\_}ice to its appropriate value
1365  ($i.e.$ $2$ for LIM-3 or $3$ for CICE).
1366\end{description}
1367
1368% {Description of Ice-ocean interface to be added here or in LIM 2 and 3 doc ?}
1369
1370\subsection{Interface to CICE (\protect\mdl{sbcice\_cice})}
1371\label{subsec:SBC_cice}
1372
1373It is now possible to couple a regional or global NEMO configuration (without AGRIF)
1374to the CICE sea-ice model by using \key{cice}.
1375The CICE code can be obtained from \href{http://oceans11.lanl.gov/trac/CICE/}{LANL} and
1376the additional 'hadgem3' drivers will be required, even with the latest code release.
1377Input grid files consistent with those used in NEMO will also be needed,
1378and CICE CPP keys \textbf{ORCA\_GRID}, \textbf{CICE\_IN\_NEMO} and \textbf{coupled} should be used
1379(seek advice from UKMO if necessary).
1380Currently the code is only designed to work when using the CORE forcing option for NEMO
1381(with \textit{calc\_strair}\forcode{ = .true.} and \textit{calc\_Tsfc}\forcode{ = .true.} in the CICE name-list),
1382or alternatively when NEMO is coupled to the HadGAM3 atmosphere model
1383(with \textit{calc\_strair}\forcode{ = .false.} and \textit{calc\_Tsfc}\forcode{ = false}).
1384The code is intended to be used with \np{nn\_fsbc} set to 1
1385(although coupling ocean and ice less frequently should work,
1386it is possible the calculation of some of the ocean-ice fluxes needs to be modified slightly -
1387the user should check that results are not significantly different to the standard case).
1388
1389There are two options for the technical coupling between NEMO and CICE.
1390The standard version allows complete flexibility for the domain decompositions in the individual models,
1391but this is at the expense of global gather and scatter operations in the coupling which
1392become very expensive on larger numbers of processors.
1393The alternative option (using \key{nemocice\_decomp} for both NEMO and CICE) ensures that
1394the domain decomposition is identical in both models (provided domain parameters are set appropriately,
1395and \textit{processor\_shape~=~square-ice} and \textit{distribution\_wght~=~block} in the CICE name-list) and
1396allows much more efficient direct coupling on individual processors.
1397This solution scales much better although it is at the expense of having more idle CICE processors in areas where
1398there is no sea ice.
1399
1400% -------------------------------------------------------------------------------------------------------------
1401%        Freshwater budget control
1402% -------------------------------------------------------------------------------------------------------------
1403\subsection{Freshwater budget control (\protect\mdl{sbcfwb})}
1404\label{subsec:SBC_fwb}
1405
1406For global ocean simulation it can be useful to introduce a control of the mean sea level in order to
1407prevent unrealistic drift of the sea surface height due to inaccuracy in the freshwater fluxes.
1408In \NEMO, two way of controlling the the freshwater budget.
1409\begin{description}
1410\item[\np{nn\_fwb}\forcode{ = 0}]
1411  no control at all.
1412  The mean sea level is free to drift, and will certainly do so.
1413\item[\np{nn\_fwb}\forcode{ = 1}]
1414  global mean \textit{emp} set to zero at each model time step.
1415%Note that with a sea-ice model, this technique only control the mean sea level with linear free surface (\key{vvl} not defined) and no mass flux between ocean and ice (as it is implemented in the current ice-ocean coupling).
1416\item[\np{nn\_fwb}\forcode{ = 2}]
1417  freshwater budget is adjusted from the previous year annual mean budget which
1418  is read in the \textit{EMPave\_old.dat} file.
1419  As the model uses the Boussinesq approximation, the annual mean fresh water budget is simply evaluated from
1420  the change in the mean sea level at January the first and saved in the \textit{EMPav.dat} file.
1421\end{description}
1422
1423% -------------------------------------------------------------------------------------------------------------
1424%        Neutral Drag Coefficient from external wave model
1425% -------------------------------------------------------------------------------------------------------------
1426\subsection[Neutral drag coeff. from external wave model (\protect\mdl{sbcwave})]
1427            {Neutral drag coefficient from external wave model (\protect\mdl{sbcwave})}
1428\label{subsec:SBC_wave}
1429%------------------------------------------namwave----------------------------------------------------
1430
1431\nlst{namsbc_wave}
1432%-------------------------------------------------------------------------------------------------------------
1433
1434In order to read a neutral drag coefficient, from an external data source ($i.e.$ a wave model),
1435the logical variable \np{ln\_cdgw} in \ngn{namsbc} namelist must be set to \forcode{.true.}.
1436The \mdl{sbcwave} module containing the routine \np{sbc\_wave} reads the namelist \ngn{namsbc\_wave}
1437(for external data names, locations, frequency, interpolation and all the miscellanous options allowed by
1438Input Data generic Interface see \autoref{sec:SBC_input}) and
1439a 2D field of neutral drag coefficient.
1440Then using the routine TURB\_CORE\_1Z or TURB\_CORE\_2Z, and starting from the neutral drag coefficent provided,
1441the drag coefficient is computed according to stable/unstable conditions of the air-sea interface following
1442\citet{Large_Yeager_Rep04}.
1443
1444
1445% Griffies doc:
1446% When running ocean-ice simulations, we are not explicitly representing land processes,
1447% such as rivers, catchment areas, snow accumulation, etc. However, to reduce model drift,
1448% it is important to balance the hydrological cycle in ocean-ice models.
1449% We thus need to prescribe some form of global normalization to the precipitation minus evaporation plus river runoff.
1450% The result of the normalization should be a global integrated zero net water input to the ocean-ice system over
1451% a chosen time scale.
1452%How often the normalization is done is a matter of choice. In mom4p1, we choose to do so at each model time step,
1453% so that there is always a zero net input of water to the ocean-ice system.
1454% Others choose to normalize over an annual cycle, in which case the net imbalance over an annual cycle is used
1455% to alter the subsequent year�s water budget in an attempt to damp the annual water imbalance.
1456% Note that the annual budget approach may be inappropriate with interannually varying precipitation forcing.
1457% When running ocean-ice coupled models, it is incorrect to include the water transport between the ocean
1458% and ice models when aiming to balance the hydrological cycle.
1459% The reason is that it is the sum of the water in the ocean plus ice that should be balanced when running ocean-ice models,
1460% not the water in any one sub-component. As an extreme example to illustrate the issue,
1461% consider an ocean-ice model with zero initial sea ice. As the ocean-ice model spins up,
1462% there should be a net accumulation of water in the growing sea ice, and thus a net loss of water from the ocean.
1463% The total water contained in the ocean plus ice system is constant, but there is an exchange of water between
1464% the subcomponents. This exchange should not be part of the normalization used to balance the hydrological cycle
1465% in ocean-ice models.
1466
1467
1468\end{document}
Note: See TracBrowser for help on using the repository browser.