source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_OBS.tex @ 11693

Last change on this file since 11693 was 11693, checked in by nicolasmartin, 12 months ago

Macros renaming

File size: 57.5 KB
Line 
1\documentclass[../main/NEMO_manual]{subfiles}
2
3\begin{document}
4
5\chapter{Observation and Model Comparison (OBS)}
6\label{chap:OBS}
7
8%\subsubsection*{Changes record}
9%\begin{tabular}{l||l|m{0.65\linewidth}}
10%    Release   & Author        & Modifications \\
11%    {\em 4.0} & {\em D. J. Lea} & {\em \NEMO\ 4.0 updates}  \\
12%    {\em 3.6} & {\em M. Martin, A. Ryan} & {\em Add averaging operator, standalone obs oper} \\
13%    {\em 3.4} & {\em D. J. Lea, M. Martin, ...} & {\em Initial version}  \\
14%    {\em --\texttt{"}--} & {\em ... K. Mogensen, A. Vidard, A. Weaver} & {\em ---\texttt{"}---}  \\
15%\end{tabular}
16
17\thispagestyle{plain}
18
19\chaptertoc
20
21\paragraph{Changes record} ~\\
22
23{\footnotesize
24  \begin{tabularx}{\textwidth}{l||X|X}
25    Release & Author(s) & Modifications \\
26    \hline
27    {\em   4.0} & {\em ...} & {\em ...} \\
28    {\em   3.6} & {\em ...} & {\em ...} \\
29    {\em   3.4} & {\em ...} & {\em ...} \\
30    {\em <=3.4} & {\em ...} & {\em ...}
31  \end{tabularx}
32}
33
34\clearpage
35
36The observation and model comparison code, the observation operator (OBS), reads in observation files
37(profile temperature and salinity, sea surface temperature, sea level anomaly, sea ice concentration, and velocity) and calculates an interpolated model equivalent value at the observation location and nearest model time step.
38The resulting data are saved in a ``feedback'' file (or files).
39The code was originally developed for use with the NEMOVAR data assimilation code,
40but can be used for validation or verification of the model or with any other data assimilation system.
41
42The OBS code is called from \mdl{nemogcm} for model initialisation and to calculate the model equivalent values for observations on the 0th time step.
43The code is then called again after each time step from \mdl{step}.
44The code is only activated if the \nam{obs}{obs} namelist logical \np{ln_diaobs}{ln\_diaobs} is set to true.
45
46For all data types a 2D horizontal interpolator or averager is needed to
47interpolate/average the model fields to the observation location.
48For {\em in situ} profiles, a 1D vertical interpolator is needed in addition to
49provide model fields at the observation depths.
50This now works in a generalised vertical coordinate system.
51
52Some profile observation types (\eg\ tropical moored buoys) are made available as daily averaged quantities.
53The observation operator code can be set-up to calculate the equivalent daily average model temperature fields using
54the \np{nn_profdavtypes}{nn\_profdavtypes} namelist array.
55Some SST observations are equivalent to a night-time average value and
56the observation operator code can calculate equivalent night-time average model SST fields by
57setting the namelist value \np{ln_sstnight}{ln\_sstnight} to true.
58Otherwise (by default) the model value from the nearest time step to the observation time is used.
59
60The code is controlled by the namelist \nam{obs}{obs}.
61See the following sections for more details on setting up the namelist.
62
63In \autoref{sec:OBS_example} a test example of the observation operator code is introduced, including
64where to obtain data and how to setup the namelist.
65In \autoref{sec:OBS_details} some more technical details of the different observation types used are introduced, and we
66also show a more complete namelist.
67In \autoref{sec:OBS_theory} some of the theoretical aspects of the observation operator are described including
68interpolation methods and running on multiple processors.
69In \autoref{sec:OBS_sao} the standalone observation operator code is described.
70In \autoref{sec:OBS_obsutils} we describe some utilities to help work with the files produced by the OBS code.
71
72%% =================================================================================================
73\section{Running the observation operator code example}
74\label{sec:OBS_example}
75
76In this section an example of running the observation operator code is described using
77profile observation data which can be freely downloaded.
78It shows how to adapt an existing run and build of \NEMO\ to run the observation operator. Note also the observation operator and the assimilation increments code are run in the ORCA2\_ICE\_OBS SETTE test.
79
80\begin{enumerate}
81\item Compile \NEMO.
82
83\item Download some EN4 data from \href{http://www.metoffice.gov.uk/hadobs}{www.metoffice.gov.uk/hadobs}.
84  Choose observations which are valid for the period of your test run because
85  the observation operator compares the model and observations for a matching date and time.
86
87\item Compile the OBSTOOLS code in the \path{tools} directory using:
88\begin{cmds}
89./maketools -n OBSTOOLS -m [ARCH]
90\end{cmds}
91
92replacing \texttt{[ARCH]} with the build architecture file for your machine. Note the tools are checked out from a separate location of the repository (under \path{/utils/tools}).
93
94\item Convert the EN4 data into feedback format:
95\begin{cmds}
96enact2fb.exe profiles_01.nc EN.4.1.1.f.profiles.g10.YYYYMM.nc
97\end{cmds}
98
99\item Include the following in the \NEMO\ namelist to run the observation operator on this data:
100\end{enumerate}
101
102Options are defined through the \nam{obs}{obs} namelist variables.
103The options \np{ln_t3d}{ln\_t3d} and \np{ln_s3d}{ln\_s3d} switch on the temperature and salinity profile observation operator code.
104The filename or array of filenames are specified using the \np{cn_profbfiles}{cn\_profbfiles} variable.
105The model grid points for a particular observation latitude and longitude are found using
106the grid searching part of the code.
107This can be expensive, particularly for large numbers of observations,
108setting \np{ln_grid_search_lookup}{ln\_grid\_search\_lookup} allows the use of a lookup table which
109is saved into an \np{cn_gridsearch}{cn\_gridsearch} file (or files).
110This will need to be generated the first time if it does not exist in the run directory.
111However, once produced it will significantly speed up future grid searches.
112Setting \np{ln_grid_global}{ln\_grid\_global} means that the code distributes the observations evenly between processors.
113Alternatively each processor will work with observations located within the model subdomain
114(see \autoref{subsec:OBS_parallel}).
115
116A number of utilities are now provided to plot the feedback files, convert and recombine the files.
117These are explained in more detail in \autoref{sec:OBS_obsutils}.
118Utilities to convert other input data formats into the feedback format are also described in
119\autoref{sec:OBS_obsutils}.
120
121%% =================================================================================================
122\section{Technical details (feedback type observation file headers)}
123\label{sec:OBS_details}
124
125Here we show a more complete example namelist \nam{obs}{obs} and also show the NetCDF headers of
126the observation files that may be used with the observation operator.
127
128\begin{listing}
129  \nlst{namobs}
130  \caption{\forcode{&namobs}}
131  \label{lst:namobs}
132\end{listing}
133
134The observation operator code uses the feedback observation file format for all data types.
135All the observation files must be in NetCDF format.
136Some example headers (produced using \mbox{\textit{ncdump~-h}}) for profile data, sea level anomaly and
137sea surface temperature are in the following subsections.
138
139%% =================================================================================================
140\subsection{Profile feedback file}
141
142\begin{clines}
143netcdf profiles_01 {
144dimensions:
145     N_OBS = 603 ;
146     N_LEVELS = 150 ;
147     N_VARS = 2 ;
148     N_QCF = 2 ;
149     N_ENTRIES = 1 ;
150     N_EXTRA = 1 ;
151     STRINGNAM = 8 ;
152     STRINGGRID = 1 ;
153     STRINGWMO = 8 ;
154     STRINGTYP = 4 ;
155     STRINGJULD = 14 ;
156variables:
157     char VARIABLES(N_VARS, STRINGNAM) ;
158          VARIABLES:long_name = "List of variables in feedback files" ;
159     char ENTRIES(N_ENTRIES, STRINGNAM) ;
160          ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
161     char EXTRA(N_EXTRA, STRINGNAM) ;
162          EXTRA:long_name = "List of extra variables" ;
163     char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
164          STATION_IDENTIFIER:long_name = "Station identifier" ;
165     char STATION_TYPE(N_OBS, STRINGTYP) ;
166          STATION_TYPE:long_name = "Code instrument type" ;
167     double LONGITUDE(N_OBS) ;
168          LONGITUDE:long_name = "Longitude" ;
169          LONGITUDE:units = "degrees_east" ;
170          LONGITUDE:_Fillvalue = 99999.f ;
171     double LATITUDE(N_OBS) ;
172          LATITUDE:long_name = "Latitude" ;
173          LATITUDE:units = "degrees_north" ;
174          LATITUDE:_Fillvalue = 99999.f ;
175     double DEPTH(N_OBS, N_LEVELS) ;
176          DEPTH:long_name = "Depth" ;
177          DEPTH:units = "metre" ;
178          DEPTH:_Fillvalue = 99999.f ;
179     int DEPTH_QC(N_OBS, N_LEVELS) ;
180          DEPTH_QC:long_name = "Quality on depth" ;
181          DEPTH_QC:Conventions = "q where q =[0,9]" ;
182          DEPTH_QC:_Fillvalue = 0 ;
183     int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
184          DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
185          DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
186     double JULD(N_OBS) ;
187          JULD:long_name = "Julian day" ;
188          JULD:units = "days since JULD_REFERENCE" ;
189          JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
190          JULD:_Fillvalue = 99999.f ;
191     char JULD_REFERENCE(STRINGJULD) ;
192          JULD_REFERENCE:long_name = "Date of reference for julian days" ;
193          JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
194     int OBSERVATION_QC(N_OBS) ;
195          OBSERVATION_QC:long_name = "Quality on observation" ;
196          OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
197          OBSERVATION_QC:_Fillvalue = 0 ;
198     int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
199          OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
200          OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
201          OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
202     int POSITION_QC(N_OBS) ;
203          POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
204          POSITION_QC:Conventions = "q where q =[0,9]" ;
205          POSITION_QC:_Fillvalue = 0 ;
206     int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
207          POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
208          POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
209          POSITION_QC_FLAGS:_Fillvalue = 0 ;
210     int JULD_QC(N_OBS) ;
211          JULD_QC:long_name = "Quality on date and time" ;
212          JULD_QC:Conventions = "q where q =[0,9]" ;
213          JULD_QC:_Fillvalue = 0 ;
214     int JULD_QC_FLAGS(N_OBS, N_QCF) ;
215          JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
216          JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
217          JULD_QC_FLAGS:_Fillvalue = 0 ;
218     int ORIGINAL_FILE_INDEX(N_OBS) ;
219          ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
220          ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
221     float POTM_OBS(N_OBS, N_LEVELS) ;
222          POTM_OBS:long_name = "Potential temperature" ;
223          POTM_OBS:units = "Degrees Celsius" ;
224          POTM_OBS:_Fillvalue = 99999.f ;
225     float POTM_Hx(N_OBS, N_LEVELS) ;
226          POTM_Hx:long_name = "Model interpolated potential temperature" ;
227          POTM_Hx:units = "Degrees Celsius" ;
228          POTM_Hx:_Fillvalue = 99999.f ;
229     int POTM_QC(N_OBS) ;
230          POTM_QC:long_name = "Quality on potential temperature" ;
231          POTM_QC:Conventions = "q where q =[0,9]" ;
232          POTM_QC:_Fillvalue = 0 ;
233     int POTM_QC_FLAGS(N_OBS, N_QCF) ;
234          POTM_QC_FLAGS:long_name = "Quality flags on potential temperature" ;
235          POTM_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
236          POTM_QC_FLAGS:_Fillvalue = 0 ;
237     int POTM_LEVEL_QC(N_OBS, N_LEVELS) ;
238          POTM_LEVEL_QC:long_name = "Quality for each level on potential temperature" ;
239          POTM_LEVEL_QC:Conventions = "q where q =[0,9]" ;
240          POTM_LEVEL_QC:_Fillvalue = 0 ;
241     int POTM_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
242          POTM_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on potential temperature" ;
243          POTM_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
244          POTM_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
245     int POTM_IOBSI(N_OBS) ;
246          POTM_IOBSI:long_name = "ORCA grid search I coordinate" ;
247     int POTM_IOBSJ(N_OBS) ;
248          POTM_IOBSJ:long_name = "ORCA grid search J coordinate" ;
249     int POTM_IOBSK(N_OBS, N_LEVELS) ;
250          POTM_IOBSK:long_name = "ORCA grid search K coordinate" ;
251     char POTM_GRID(STRINGGRID) ;
252          POTM_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
253     float PSAL_OBS(N_OBS, N_LEVELS) ;
254          PSAL_OBS:long_name = "Practical salinity" ;
255          PSAL_OBS:units = "PSU" ;
256          PSAL_OBS:_Fillvalue = 99999.f ;
257     float PSAL_Hx(N_OBS, N_LEVELS) ;
258          PSAL_Hx:long_name = "Model interpolated practical salinity" ;
259          PSAL_Hx:units = "PSU" ;
260          PSAL_Hx:_Fillvalue = 99999.f ;
261     int PSAL_QC(N_OBS) ;
262          PSAL_QC:long_name = "Quality on practical salinity" ;
263          PSAL_QC:Conventions = "q where q =[0,9]" ;
264          PSAL_QC:_Fillvalue = 0 ;
265     int PSAL_QC_FLAGS(N_OBS, N_QCF) ;
266          PSAL_QC_FLAGS:long_name = "Quality flags on practical salinity" ;
267          PSAL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
268          PSAL_QC_FLAGS:_Fillvalue = 0 ;
269     int PSAL_LEVEL_QC(N_OBS, N_LEVELS) ;
270          PSAL_LEVEL_QC:long_name = "Quality for each level on practical salinity" ;
271          PSAL_LEVEL_QC:Conventions = "q where q =[0,9]" ;
272          PSAL_LEVEL_QC:_Fillvalue = 0 ;
273     int PSAL_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
274          PSAL_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on practical salinity" ;
275          PSAL_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
276          PSAL_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
277     int PSAL_IOBSI(N_OBS) ;
278          PSAL_IOBSI:long_name = "ORCA grid search I coordinate" ;
279     int PSAL_IOBSJ(N_OBS) ;
280          PSAL_IOBSJ:long_name = "ORCA grid search J coordinate" ;
281     int PSAL_IOBSK(N_OBS, N_LEVELS) ;
282          PSAL_IOBSK:long_name = "ORCA grid search K coordinate" ;
283     char PSAL_GRID(STRINGGRID) ;
284          PSAL_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
285     float TEMP(N_OBS, N_LEVELS) ;
286          TEMP:long_name = "Insitu temperature" ;
287          TEMP:units = "Degrees Celsius" ;
288          TEMP:_Fillvalue = 99999.f ;
289
290// global attributes:
291          :title = "NEMO observation operator output" ;
292          :Convention = "NEMO unified observation operator output" ;
293}
294\end{clines}
295
296%% =================================================================================================
297\subsection{Sea level anomaly feedback file}
298
299\begin{clines}
300netcdf sla_01 {
301dimensions:
302     N_OBS = 41301 ;
303     N_LEVELS = 1 ;
304     N_VARS = 1 ;
305     N_QCF = 2 ;
306     N_ENTRIES = 1 ;
307     N_EXTRA = 1 ;
308     STRINGNAM = 8 ;
309     STRINGGRID = 1 ;
310     STRINGWMO = 8 ;
311     STRINGTYP = 4 ;
312     STRINGJULD = 14 ;
313variables:
314     char VARIABLES(N_VARS, STRINGNAM) ;
315          VARIABLES:long_name = "List of variables in feedback files" ;
316     char ENTRIES(N_ENTRIES, STRINGNAM) ;
317          ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
318     char EXTRA(N_EXTRA, STRINGNAM) ;
319          EXTRA:long_name = "List of extra variables" ;
320     char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
321          STATION_IDENTIFIER:long_name = "Station identifier" ;
322     char STATION_TYPE(N_OBS, STRINGTYP) ;
323          STATION_TYPE:long_name = "Code instrument type" ;
324     double LONGITUDE(N_OBS) ;
325          LONGITUDE:long_name = "Longitude" ;
326          LONGITUDE:units = "degrees_east" ;
327          LONGITUDE:_Fillvalue = 99999.f ;
328     double LATITUDE(N_OBS) ;
329          LATITUDE:long_name = "Latitude" ;
330          LATITUDE:units = "degrees_north" ;
331          LATITUDE:_Fillvalue = 99999.f ;
332     double DEPTH(N_OBS, N_LEVELS) ;
333          DEPTH:long_name = "Depth" ;
334          DEPTH:units = "metre" ;
335          DEPTH:_Fillvalue = 99999.f ;
336     int DEPTH_QC(N_OBS, N_LEVELS) ;
337          DEPTH_QC:long_name = "Quality on depth" ;
338          DEPTH_QC:Conventions = "q where q =[0,9]" ;
339          DEPTH_QC:_Fillvalue = 0 ;
340     int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
341          DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
342          DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
343     double JULD(N_OBS) ;
344          JULD:long_name = "Julian day" ;
345          JULD:units = "days since JULD_REFERENCE" ;
346          JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
347          JULD:_Fillvalue = 99999.f ;
348     char JULD_REFERENCE(STRINGJULD) ;
349          JULD_REFERENCE:long_name = "Date of reference for julian days" ;
350          JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
351     int OBSERVATION_QC(N_OBS) ;
352          OBSERVATION_QC:long_name = "Quality on observation" ;
353          OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
354          OBSERVATION_QC:_Fillvalue = 0 ;
355     int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
356          OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
357          OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
358          OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
359     int POSITION_QC(N_OBS) ;
360          POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
361          POSITION_QC:Conventions = "q where q =[0,9]" ;
362          POSITION_QC:_Fillvalue = 0 ;
363     int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
364          POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
365          POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
366          POSITION_QC_FLAGS:_Fillvalue = 0 ;
367     int JULD_QC(N_OBS) ;
368          JULD_QC:long_name = "Quality on date and time" ;
369          JULD_QC:Conventions = "q where q =[0,9]" ;
370          JULD_QC:_Fillvalue = 0 ;
371     int JULD_QC_FLAGS(N_OBS, N_QCF) ;
372          JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
373          JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
374          JULD_QC_FLAGS:_Fillvalue = 0 ;
375     int ORIGINAL_FILE_INDEX(N_OBS) ;
376          ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
377          ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
378     float SLA_OBS(N_OBS, N_LEVELS) ;
379          SLA_OBS:long_name = "Sea level anomaly" ;
380          SLA_OBS:units = "metre" ;
381          SLA_OBS:_Fillvalue = 99999.f ;
382     float SLA_Hx(N_OBS, N_LEVELS) ;
383          SLA_Hx:long_name = "Model interpolated sea level anomaly" ;
384          SLA_Hx:units = "metre" ;
385          SLA_Hx:_Fillvalue = 99999.f ;
386     int SLA_QC(N_OBS) ;
387          SLA_QC:long_name = "Quality on sea level anomaly" ;
388          SLA_QC:Conventions = "q where q =[0,9]" ;
389          SLA_QC:_Fillvalue = 0 ;
390     int SLA_QC_FLAGS(N_OBS, N_QCF) ;
391          SLA_QC_FLAGS:long_name = "Quality flags on sea level anomaly" ;
392          SLA_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
393          SLA_QC_FLAGS:_Fillvalue = 0 ;
394     int SLA_LEVEL_QC(N_OBS, N_LEVELS) ;
395          SLA_LEVEL_QC:long_name = "Quality for each level on sea level anomaly" ;
396          SLA_LEVEL_QC:Conventions = "q where q =[0,9]" ;
397          SLA_LEVEL_QC:_Fillvalue = 0 ;
398     int SLA_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
399          SLA_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on sea level anomaly" ;
400          SLA_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
401          SLA_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
402     int SLA_IOBSI(N_OBS) ;
403          SLA_IOBSI:long_name = "ORCA grid search I coordinate" ;
404     int SLA_IOBSJ(N_OBS) ;
405          SLA_IOBSJ:long_name = "ORCA grid search J coordinate" ;
406     int SLA_IOBSK(N_OBS, N_LEVELS) ;
407          SLA_IOBSK:long_name = "ORCA grid search K coordinate" ;
408     char SLA_GRID(STRINGGRID) ;
409          SLA_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
410     float MDT(N_OBS, N_LEVELS) ;
411          MDT:long_name = "Mean Dynamic Topography" ;
412          MDT:units = "metre" ;
413          MDT:_Fillvalue = 99999.f ;
414
415// global attributes:
416          :title = "NEMO observation operator output" ;
417          :Convention = "NEMO unified observation operator output" ;
418}
419\end{clines}
420
421To use Sea Level Anomaly (SLA) data the mean dynamic topography (MDT) must be provided in a separate file defined on
422the model grid called \ifile{slaReferenceLevel}.
423The MDT is required in order to produce the model equivalent sea level anomaly from the model sea surface height.
424Below is an example header for this file (on the ORCA025 grid).
425
426\begin{clines}
427dimensions:
428        x = 1442 ;
429        y = 1021 ;
430variables:
431        float nav_lon(y, x) ;
432                nav_lon:units = "degrees_east" ;
433        float nav_lat(y, x) ;
434                nav_lat:units = "degrees_north" ;
435        float sossheig(y, x) ;
436                sossheig:_FillValue = -1.e+30f ;
437                sossheig:coordinates = "nav_lon nav_lat" ;
438                sossheig:long_name = "Mean Dynamic Topography" ;
439                sossheig:units = "metres" ;
440                sossheig:grid = "orca025T" ;
441\end{clines}
442
443%% =================================================================================================
444\subsection{Sea surface temperature feedback file}
445
446\begin{clines}
447netcdf sst_01 {
448dimensions:
449     N_OBS = 33099 ;
450     N_LEVELS = 1 ;
451     N_VARS = 1 ;
452     N_QCF = 2 ;
453     N_ENTRIES = 1 ;
454     STRINGNAM = 8 ;
455     STRINGGRID = 1 ;
456     STRINGWMO = 8 ;
457     STRINGTYP = 4 ;
458     STRINGJULD = 14 ;
459variables:
460     char VARIABLES(N_VARS, STRINGNAM) ;
461          VARIABLES:long_name = "List of variables in feedback files" ;
462     char ENTRIES(N_ENTRIES, STRINGNAM) ;
463          ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
464     char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
465          STATION_IDENTIFIER:long_name = "Station identifier" ;
466     char STATION_TYPE(N_OBS, STRINGTYP) ;
467          STATION_TYPE:long_name = "Code instrument type" ;
468     double LONGITUDE(N_OBS) ;
469          LONGITUDE:long_name = "Longitude" ;
470          LONGITUDE:units = "degrees_east" ;
471          LONGITUDE:_Fillvalue = 99999.f ;
472     double LATITUDE(N_OBS) ;
473          LATITUDE:long_name = "Latitude" ;
474          LATITUDE:units = "degrees_north" ;
475          LATITUDE:_Fillvalue = 99999.f ;
476     double DEPTH(N_OBS, N_LEVELS) ;
477          DEPTH:long_name = "Depth" ;
478          DEPTH:units = "metre" ;
479          DEPTH:_Fillvalue = 99999.f ;
480     int DEPTH_QC(N_OBS, N_LEVELS) ;
481          DEPTH_QC:long_name = "Quality on depth" ;
482          DEPTH_QC:Conventions = "q where q =[0,9]" ;
483          DEPTH_QC:_Fillvalue = 0 ;
484     int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
485          DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
486          DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
487     double JULD(N_OBS) ;
488          JULD:long_name = "Julian day" ;
489          JULD:units = "days since JULD_REFERENCE" ;
490          JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
491          JULD:_Fillvalue = 99999.f ;
492     char JULD_REFERENCE(STRINGJULD) ;
493          JULD_REFERENCE:long_name = "Date of reference for julian days" ;
494          JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
495     int OBSERVATION_QC(N_OBS) ;
496          OBSERVATION_QC:long_name = "Quality on observation" ;
497          OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
498          OBSERVATION_QC:_Fillvalue = 0 ;
499     int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
500          OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
501          OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
502          OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
503     int POSITION_QC(N_OBS) ;
504          POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
505          POSITION_QC:Conventions = "q where q =[0,9]" ;
506          POSITION_QC:_Fillvalue = 0 ;
507     int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
508          POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
509          POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
510          POSITION_QC_FLAGS:_Fillvalue = 0 ;
511     int JULD_QC(N_OBS) ;
512          JULD_QC:long_name = "Quality on date and time" ;
513          JULD_QC:Conventions = "q where q =[0,9]" ;
514          JULD_QC:_Fillvalue = 0 ;
515     int JULD_QC_FLAGS(N_OBS, N_QCF) ;
516          JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
517          JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
518          JULD_QC_FLAGS:_Fillvalue = 0 ;
519     int ORIGINAL_FILE_INDEX(N_OBS) ;
520          ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
521          ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
522     float SST_OBS(N_OBS, N_LEVELS) ;
523          SST_OBS:long_name = "Sea surface temperature" ;
524          SST_OBS:units = "Degree centigrade" ;
525          SST_OBS:_Fillvalue = 99999.f ;
526     float SST_Hx(N_OBS, N_LEVELS) ;
527          SST_Hx:long_name = "Model interpolated sea surface temperature" ;
528          SST_Hx:units = "Degree centigrade" ;
529          SST_Hx:_Fillvalue = 99999.f ;
530     int SST_QC(N_OBS) ;
531          SST_QC:long_name = "Quality on sea surface temperature" ;
532          SST_QC:Conventions = "q where q =[0,9]" ;
533          SST_QC:_Fillvalue = 0 ;
534     int SST_QC_FLAGS(N_OBS, N_QCF) ;
535          SST_QC_FLAGS:long_name = "Quality flags on sea surface temperature" ;
536          SST_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
537          SST_QC_FLAGS:_Fillvalue = 0 ;
538     int SST_LEVEL_QC(N_OBS, N_LEVELS) ;
539          SST_LEVEL_QC:long_name = "Quality for each level on sea surface temperature" ;
540          SST_LEVEL_QC:Conventions = "q where q =[0,9]" ;
541          SST_LEVEL_QC:_Fillvalue = 0 ;
542     int SST_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
543          SST_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on sea surface temperature" ;
544          SST_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
545          SST_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
546     int SST_IOBSI(N_OBS) ;
547          SST_IOBSI:long_name = "ORCA grid search I coordinate" ;
548     int SST_IOBSJ(N_OBS) ;
549          SST_IOBSJ:long_name = "ORCA grid search J coordinate" ;
550     int SST_IOBSK(N_OBS, N_LEVELS) ;
551          SST_IOBSK:long_name = "ORCA grid search K coordinate" ;
552     char SST_GRID(STRINGGRID) ;
553          SST_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
554
555// global attributes:
556          :title = "NEMO observation operator output" ;
557          :Convention = "NEMO unified observation operator output" ;
558}
559\end{clines}
560
561%% =================================================================================================
562\section{Theoretical details}
563\label{sec:OBS_theory}
564
565%% =================================================================================================
566\subsection{Horizontal interpolation and averaging methods}
567
568For most observation types, the horizontal extent of the observation is small compared to the model grid size and so
569the model equivalent of the observation is calculated by interpolating from
570the four surrounding grid points to the observation location.
571Some satellite observations (\eg\ microwave satellite SST data, or satellite SSS data) have a footprint which
572is similar in size or larger than the model grid size (particularly when the grid size is small).
573In those cases the model counterpart should be calculated by averaging the model grid points over
574the same size as the footprint.
575\NEMO\ therefore has the capability to specify either an interpolation or an averaging
576(for surface observation types only).
577
578The main namelist option associated with the interpolation/averaging is \np{nn_2dint}{nn\_2dint}.
579This default option can be set to values from 0 to 6.
580Values between 0 to 4 are associated with interpolation while values 5 or 6 are associated with averaging.
581\begin{itemize}
582\item \np[=0]{nn_2dint}{nn\_2dint}: Distance-weighted interpolation
583\item \np[=1]{nn_2dint}{nn\_2dint}: Distance-weighted interpolation (small angle)
584\item \np[=2]{nn_2dint}{nn\_2dint}: Bilinear interpolation (geographical grid)
585\item \np[=3]{nn_2dint}{nn\_2dint}: Bilinear remapping interpolation (general grid)
586\item \np[=4]{nn_2dint}{nn\_2dint}: Polynomial interpolation
587\item \np[=5]{nn_2dint}{nn\_2dint}: Radial footprint averaging with diameter specified in the namelist as
588  \texttt{rn\_[var]\_avglamscl} in degrees or metres (set using \texttt{ln\_[var]\_fp\_indegs})
589\item \np[=6]{nn_2dint}{nn\_2dint}: Rectangular footprint averaging with E/W and N/S size specified in
590  the namelist as \texttt{rn\_[var]\_avglamscl} and \texttt{rn\_[var]\_avgphiscl} in degrees or metres
591  (set using \texttt{ln\_[var]\_fp\_indegs})
592\end{itemize}
593Replace \texttt{[var]} in the last two options with the observation type (sla, sst, sss or sic) for
594which the averaging is to be performed (see namelist example above).
595The \np{nn_2dint}{nn\_2dint} default option can be overridden for surface observation types using
596namelist values \texttt{nn\_2dint\_[var]} where \texttt{[var]} is the observation type.
597
598Below is some more detail on the various options for interpolation and averaging available in \NEMO.
599
600%% =================================================================================================
601\subsubsection{Horizontal interpolation}
602
603Consider an observation point ${\mathrm P}$ with longitude and latitude (${\lambda_{}}_{\mathrm P}$, $\phi_{\mathrm P}$) and
604the four nearest neighbouring model grid points ${\mathrm A}$, ${\mathrm B}$, ${\mathrm C}$ and ${\mathrm D}$ with
605longitude and latitude ($\lambda_{\mathrm A}$, $\phi_{\mathrm A}$),($\lambda_{\mathrm B}$, $\phi_{\mathrm B}$) etc.
606All horizontal interpolation methods implemented in \NEMO\ estimate the value of a model variable $x$ at point $P$ as
607a weighted linear combination of the values of the model variables at the grid points ${\mathrm A}$, ${\mathrm B}$ etc.:
608
609\begin{align*}
610  {x_{}}_{\mathrm P} =
611\frac{1}{w} \left( {w_{}}_{\mathrm A} {x_{}}_{\mathrm A} +
612{w_{}}_{\mathrm B} {x_{}}_{\mathrm B} +
613{w_{}}_{\mathrm C} {x_{}}_{\mathrm C} +
614{w_{}}_{\mathrm D} {x_{}}_{\mathrm D} \right)
615\end{align*}
616
617where ${w_{}}_{\mathrm A}$, ${w_{}}_{\mathrm B}$ etc. are the respective weights for the model field at
618points ${\mathrm A}$, ${\mathrm B}$ etc., and $w = {w_{}}_{\mathrm A} + {w_{}}_{\mathrm B} + {w_{}}_{\mathrm C} + {w_{}}_{\mathrm D}$.
619
620Four different possibilities are available for computing the weights.
621
622\begin{enumerate}
623\item {\bfseries Great-Circle distance-weighted interpolation.}
624  The weights are computed as a function of the great-circle distance $s(P, \cdot)$ between $P$ and
625  the model grid points $A$, $B$ etc.
626  For example, the weight given to the field ${x_{}}_{\mathrm A}$ is specified as the product of the distances
627  from ${\mathrm P}$ to the other points:
628
629  \begin{alignat*}{2}
630    {w_{}}_{\mathrm A} = s({\mathrm P}, {\mathrm B}) \, s({\mathrm P}, {\mathrm C}) \, s({\mathrm P}, {\mathrm D})
631  \end{alignat*}
632
633  where
634
635  \begin{alignat*}{2}
636    s\left({\mathrm P}, {\mathrm M} \right) & = & \hspace{0.25em} \cos^{-1} \! \left\{
637               \sin {\phi_{}}_{\mathrm P} \sin {\phi_{}}_{\mathrm M}
638             + \cos {\phi_{}}_{\mathrm P} \cos {\phi_{}}_{\mathrm M}
639               \cos ({\lambda_{}}_{\mathrm M} - {\lambda_{}}_{\mathrm P})
640                   \right\}
641   \end{alignat*}
642
643   and $M$ corresponds to $B$, $C$ or $D$.
644   A more stable form of the great-circle distance formula for small distances ($x$ near 1)
645   involves the arcsine function (\eg\ see p.~101 of \citet{daley.barker_bk01}:
646
647   \begin{alignat*}{2}
648     s\left( {\mathrm P}, {\mathrm M} \right) = \sin^{-1} \! \left\{ \sqrt{ 1 - x^2 } \right\}
649   \end{alignat*}
650
651   where
652
653   \begin{alignat*}{2}
654     x = {a_{}}_{\mathrm M} {a_{}}_{\mathrm P} + {b_{}}_{\mathrm M} {b_{}}_{\mathrm P} + {c_{}}_{\mathrm M} {c_{}}_{\mathrm P}
655   \end{alignat*}
656
657   and
658
659   \begin{alignat*}{3}
660   & {a_{}}_{\mathrm M} & = && \quad \sin {\phi_{}}_{\mathrm M}, \\
661   & {a_{}}_{\mathrm P} & = && \quad \sin {\phi_{}}_{\mathrm P}, \\
662   & {b_{}}_{\mathrm M} & = && \quad \cos {\phi_{}}_{\mathrm M} \cos {\phi_{}}_{\mathrm M}, \\
663   & {b_{}}_{\mathrm P} & = && \quad \cos {\phi_{}}_{\mathrm P} \cos {\phi_{}}_{\mathrm P}, \\
664   & {c_{}}_{\mathrm M} & = && \quad \cos {\phi_{}}_{\mathrm M} \sin {\phi_{}}_{\mathrm M}, \\
665   & {c_{}}_{\mathrm P} & = && \quad \cos {\phi_{}}_{\mathrm P} \sin {\phi_{}}_{\mathrm P}.
666   \end{alignat*}
667
668\item {\bfseries Great-Circle distance-weighted interpolation with small angle approximation.}
669  Similar to the previous interpolation but with the distance $s$ computed as
670  \begin{alignat*}{2}
671    s\left( {\mathrm P}, {\mathrm M} \right)
672    & = & \sqrt{ \left( {\phi_{}}_{\mathrm M} - {\phi_{}}_{\mathrm P} \right)^{2}
673                                      + \left( {\lambda_{}}_{\mathrm M} - {\lambda_{}}_{\mathrm P} \right)^{2}
674                                      \cos^{2} {\phi_{}}_{\mathrm M} }
675  \end{alignat*}
676  where $M$ corresponds to $A$, $B$, $C$ or $D$.
677
678\item {\bfseries Bilinear interpolation for a regular spaced grid.}
679  The interpolation is split into two 1D interpolations in the longitude and latitude directions, respectively.
680
681\item {\bfseries Bilinear remapping interpolation for a general grid.}
682  An iterative scheme that involves first mapping a quadrilateral cell into
683  a cell with coordinates (0,0), (1,0), (0,1) and (1,1).
684  This method is based on the \href{https://github.com/SCRIP-Project/SCRIP}{SCRIP interpolation package}.
685
686\end{enumerate}
687
688%% =================================================================================================
689\subsubsection{Horizontal averaging}
690
691For each surface observation type:
692\begin{itemize}
693\item The standard grid-searching code is used to find the nearest model grid point to the observation location
694  (see next subsection).
695\item The maximum number of grid points required for that observation in each local grid domain is calculated. Some of these points may later turn out to have zero weight depending on the shape of the footprint.
696\item The longitudes and latitudes of the grid points surrounding the nearest model grid box are extracted using
697  existing MPI routines.
698\item The weights for each grid point associated with each observation are calculated,
699  either for radial or rectangular footprints.
700  For grid points completely within the footprint, the weight is one;
701  for grid points completely outside the footprint, the weight is zero.
702  For grid points which are partly within the footprint the ratio between the area of the footprint within
703  the grid box and the total area of the grid box is used as the weight.
704\item The weighted average of the model grid points associated with each observation is calculated,
705  and this is then given as the model counterpart of the observation.
706\end{itemize}
707
708Examples of the weights calculated for an observation with rectangular and radial footprints are shown in
709\autoref{fig:OBS_avgrec} and~\autoref{fig:OBS_avgrad}.
710
711\begin{figure}
712  \centering
713  \includegraphics[width=0.66\textwidth]{OBS_avg_rec}
714  \caption[Observational weights with a rectangular footprint]{
715    Weights associated with each model grid box (blue lines and numbers)
716    for an observation at -170.5\deg{E}, 56.0\deg{N} with a rectangular footprint of 1\deg\ x 1\deg.}
717  \label{fig:OBS_avgrec}
718\end{figure}
719
720\begin{figure}
721  \centering
722  \includegraphics[width=0.66\textwidth]{OBS_avg_rad}
723  \caption[Observational weights with a radial footprint]{
724    Weights associated with each model grid box (blue lines and numbers)
725    for an observation at -170.5\deg{E}, 56.0\deg{N} with a radial footprint with diameter 1\deg.}
726  \label{fig:OBS_avgrad}
727\end{figure}
728
729%% =================================================================================================
730\subsection{Grid search}
731
732For many grids used by the \NEMO\ model, such as the ORCA family, the horizontal grid coordinates $i$ and $j$ are not simple functions of latitude and longitude.
733Therefore, it is not always straightforward to determine the grid points surrounding any given observational position.
734Before the interpolation can be performed, a search algorithm is then required to determine the corner points of
735the quadrilateral cell in which the observation is located.
736This is the most difficult and time consuming part of the 2D interpolation procedure.
737A robust test for determining if an observation falls within a given quadrilateral cell is as follows.
738Let ${\mathrm P}({\lambda_{}}_{\mathrm P} ,{\phi_{}}_{\mathrm P} )$ denote the observation point,
739and let ${\mathrm A}({\lambda_{}}_{\mathrm A} ,{\phi_{}}_{\mathrm A} )$, ${\mathrm B}({\lambda_{}}_{\mathrm B} ,{\phi_{}}_{\mathrm B} )$,
740${\mathrm C}({\lambda_{}}_{\mathrm C} ,{\phi_{}}_{\mathrm C} )$ and ${\mathrm D}({\lambda_{}}_{\mathrm D} ,{\phi_{}}_{\mathrm D} )$
741denote the bottom left, bottom right, top left and top right corner points of the cell, respectively.
742To determine if P is inside the cell, we verify that the cross-products
743\begin{align*}
744  \begin{array}{lllll}
745    {{\mathbf r}_{}}_{\mathrm PA} \times {{\mathbf r}_{}}_{\mathrm PC}
746    & = & [({\lambda_{}}_{\mathrm A}\; -\; {\lambda_{}}_{\mathrm P} )
747          ({\phi_{}}_{\mathrm C}   \; -\; {\phi_{}}_{\mathrm P} )
748          - ({\lambda_{}}_{\mathrm C}\; -\; {\lambda_{}}_{\mathrm P} )
749          ({\phi_{}}_{\mathrm A}   \; -\; {\phi_{}}_{\mathrm P} )] \; \widehat{\mathbf k} \\
750    {{\mathbf r}_{}}_{\mathrm PB} \times {{\mathbf r}_{}}_{\mathrm PA}
751    & = & [({\lambda_{}}_{\mathrm B}\; -\; {\lambda_{}}_{\mathrm P} )
752          ({\phi_{}}_{\mathrm A}   \; -\; {\phi_{}}_{\mathrm P} )
753          - ({\lambda_{}}_{\mathrm A}\; -\; {\lambda_{}}_{\mathrm P} )
754          ({\phi_{}}_{\mathrm B}   \; -\; {\phi_{}}_{\mathrm P} )] \; \widehat{\mathbf k} \\
755    {{\mathbf r}_{}}_{\mathrm PC} \times {{\mathbf r}_{}}_{\mathrm PD}
756    & = & [({\lambda_{}}_{\mathrm C}\; -\; {\lambda_{}}_{\mathrm P} )
757          ({\phi_{}}_{\mathrm D}   \; -\; {\phi_{}}_{\mathrm P} )
758          - ({\lambda_{}}_{\mathrm D}\; -\; {\lambda_{}}_{\mathrm P} )
759          ({\phi_{}}_{\mathrm C}   \; -\; {\phi_{}}_{\mathrm P} )] \; \widehat{\mathbf k} \\
760    {{\mathbf r}_{}}_{\mathrm PD} \times {{\mathbf r}_{}}_{\mathrm PB}
761    & = & [({\lambda_{}}_{\mathrm D}\; -\; {\lambda_{}}_{\mathrm P} )
762          ({\phi_{}}_{\mathrm B}   \; -\; {\phi_{}}_{\mathrm P} )
763          - ({\lambda_{}}_{\mathrm B}\; -\; {\lambda_{}}_{\mathrm P} )
764          ({\phi_{}}_{\mathrm D}  \;  - \; {\phi_{}}_{\mathrm P} )] \; \widehat{\mathbf k} \\
765  \end{array}
766  % \label{eq:OBS_cross}
767\end{align*}
768point in the opposite direction to the unit normal $\widehat{\mathbf k}$
769(\ie\ that the coefficients of $\widehat{\mathbf k}$ are negative),
770where ${{\mathbf r}_{}}_{\mathrm PA}$, ${{\mathbf r}_{}}_{\mathrm PB}$, etc. correspond to
771the vectors between points P and A, P and B, etc..
772The method used is similar to the method used in the \href{https://github.com/SCRIP-Project/SCRIP}{SCRIP interpolation package}.
773
774In order to speed up the grid search, there is the possibility to construct a lookup table for a user specified resolution.
775This lookup table contains the lower and upper bounds on the $i$ and $j$ indices to
776be searched for on a regular grid.
777For each observation position, the closest point on the regular grid of this position is computed and
778the $i$ and $j$ ranges of this point searched to determine the precise four points surrounding the observation.
779
780%% =================================================================================================
781\subsection{Parallel aspects of horizontal interpolation}
782\label{subsec:OBS_parallel}
783
784For horizontal interpolation, there is the basic problem that
785the observations are unevenly distributed on the globe.
786In \NEMO\ the model grid is divided into subgrids (or domains) where
787each subgrid is executed on a single processing element with explicit message passing for
788exchange of information along the domain boundaries when running on a massively parallel processor (MPP) system.
789
790For observations there is no natural distribution since the observations are not equally distributed on the globe.
791Two options have been made available:
7921) geographical distribution;
793and 2) round-robin.
794
795%% =================================================================================================
796\subsubsection{Geographical distribution of observations among processors}
797
798\begin{figure}
799  \centering
800  \includegraphics[width=0.66\textwidth]{OBS_obsdist_local}
801  \caption[Observations with the geographical distribution]{
802    Example of the distribution of observations with
803    the geographical distribution of observational data}
804  \label{fig:OBS_local}
805\end{figure}
806
807This is the simplest option in which the observations are distributed according to
808the domain of the grid-point parallelization.
809\autoref{fig:OBS_local} shows an example of the distribution of the {\em in situ} data on processors with
810a different colour for each observation on a given processor for a 4 $\times$ 2 decomposition with ORCA2.
811The grid-point domain decomposition is clearly visible on the plot.
812
813The advantage of this approach is that all information needed for horizontal interpolation is available without
814any MPP communication.
815This is under the assumption that we are dealing with point observations and only using a $2 \times 2$ grid-point stencil for
816the interpolation (\eg\ bilinear interpolation).
817For higher order interpolation schemes this is no longer valid.
818A disadvantage with the above scheme is that the number of observations on each processor can be very different.
819If the cost of the actual interpolation is expensive relative to the communication of data needed for interpolation,
820this could lead to load imbalance.
821
822%% =================================================================================================
823\subsubsection{Round-robin distribution of observations among processors}
824
825\begin{figure}
826  \centering
827  \includegraphics[width=0.66\textwidth]{OBS_obsdist_global}
828  \caption[Observations with the round-robin distribution]{
829    Example of the distribution of observations with
830    the round-robin distribution of observational data.}
831  \label{fig:OBS_global}
832\end{figure}
833
834An alternative approach is to distribute the observations equally among processors and
835use message passing in order to retrieve the stencil for interpolation.
836The simplest distribution of the observations is to distribute them using a round-robin scheme.
837\autoref{fig:OBS_global} shows the distribution of the {\em in situ} data on processors for
838the round-robin distribution of observations with a different colour for each observation on a given processor for
839a 4 $\times$ 2 decomposition with ORCA2 for the same input data as in \autoref{fig:OBS_local}.
840The observations are now clearly randomly distributed on the globe.
841In order to be able to perform horizontal interpolation in this case,
842a subroutine has been developed that retrieves any grid points in the global space.
843
844%% =================================================================================================
845\subsection{Vertical interpolation operator}
846
847Vertical interpolation is achieved using either a cubic spline or linear interpolation.
848For the cubic spline, the top and bottom boundary conditions for the second derivative of
849the interpolating polynomial in the spline are set to zero.
850At the bottom boundary, this is done using the land-ocean mask.
851
852For profile observation types we do both vertical and horizontal interpolation. \NEMO\ has a generalised vertical coordinate system this means the vertical level depths can vary with location. Therefore, it is necessary first to perform vertical interpolation of the model value to the observation depths for each of the four surrounding grid points. After this the model values, at these points, at the observation depth, are horizontally interpolated to the observation location.
853
854%\usepackage{framed}
855
856%% =================================================================================================
857\section{Standalone observation operator}
858\label{sec:OBS_sao}
859
860%% =================================================================================================
861\subsection{Concept}
862
863The observation operator maps model variables to observation space. This is normally done while the model is running, i.e. online, it is possible to apply this mapping offline without running the model with the \textbf{standalone observation operator} (SAO). The process is divided into an initialisation phase, an interpolation phase and an output phase.
864During the interpolation phase the SAO populates the model arrays by
865reading saved model fields from disk. The interpolation and the output phases use the same OBS code described in the preceding sections.
866
867There are two ways of exploiting the standalone capacity.
868The first is to mimic the behaviour of the online system by supplying model fields at
869regular intervals between the start and the end of the run.
870This approach results in a single model counterpart per observation.
871This kind of usage produces feedback files the same file format as the online observation operator.
872The second is to take advantage of the ability to run offline by calculating
873multiple model counterparts for each observation.
874In this case it is possible to consider all forecasts verifying at the same time.
875By forecast, we mean any method which produces an estimate of physical reality which is not an observed value.
876
877% sao.exe
878
879%% =================================================================================================
880\subsection{Using the standalone observation operator}
881
882%% =================================================================================================
883\subsubsection{Building}
884
885In addition to \emph{OPA\_SRC} the SAO requires the inclusion of the \emph{SAO\_SRC} directory.
886\emph{SAO\_SRC} contains a replacement \mdl{nemo} and \mdl{nemogcm} which
887overwrites the resultant \textbf{nemo.exe}.
888Note this a similar approach to that taken by the standalone surface scheme \emph{SAS\_SRC} and the offline TOP model \emph{OFF\_SRC}.
889
890% Running
891%% =================================================================================================
892\subsubsection{Running}
893
894The simplest way to use the executable is to edit and append the \textbf{sao.nml} namelist to
895a full \NEMO\ namelist and then to run the executable as if it were nemo.exe.
896
897% Configuration section
898%% =================================================================================================
899\subsection{Configuring the standalone observation operator}
900The observation files and settings understood by \nam{obs}{obs} have been outlined in the online observation operator section.
901In addition is a further namelist \nam{sao}{sao} which used to set the input model fields for the SAO
902
903%% =================================================================================================
904\subsubsection{Single field}
905
906In the SAO the model arrays are populated at appropriate time steps via input files.
907At present, \textbf{tsn} and \textbf{sshn} are populated by the default read routines.
908These routines will be expanded upon in future versions to allow the specification of any model variable.
909As such, input files must be global versions of the model domain with
910\textbf{votemper}, \textbf{vosaline} and optionally \textbf{sshn} present.
911
912For each field read there must be an entry in the \nam{sao}{sao} namelist specifying
913the name of the file to read and the index along the \emph{time\_counter}.
914For example, to read the second time counter from a single file the namelist would be.
915
916\begin{forlines}
917!----------------------------------------------------------------------
918!       namsao Standalone obs_oper namelist
919!----------------------------------------------------------------------
920!   sao_files    specifies the files containing the model counterpart
921!   nn_sao_idx   specifies the time_counter index within the model file
922&namsao
923   sao_files = "foo.nc"
924   nn_sao_idx = 2
925/
926\end{forlines}
927
928%% =================================================================================================
929\subsubsection{Multiple fields per run}
930
931Model field iteration is controlled via \textbf{nn\_sao\_freq} which
932specifies the number of model steps at which the next field gets read.
933For example, if 12 hourly fields are to be interpolated in a setup where 288 steps equals 24 hours.
934
935\begin{forlines}
936!----------------------------------------------------------------------
937!       namsao Standalone obs_oper namelist
938!----------------------------------------------------------------------
939!   sao_files    specifies the files containing the model counterpart
940!   nn_sao_idx   specifies the time_counter index within the model file
941!   nn_sao_freq  specifies number of time steps between read operations
942&namsao
943   sao_files = "foo.nc" "foo.nc"
944   nn_sao_idx = 1 2
945   nn_sao_freq = 144
946/
947\end{forlines}
948
949The above namelist will result in feedback files whose first 12 hours contain the first field of foo.nc and
950the second 12 hours contain the second field.
951
952%\begin{framed}
953\textbf{Note} Missing files can be denoted as "nofile".
954%\end{framed}
955
956A collection of fields taken from a number of files at different indices can be combined at
957a particular frequency in time to generate a pseudo model evolution.
958If all that is needed is a single model counterpart at a regular interval then
959the standard SAO is all that is required.
960However, just to note, it is possible to extend this approach by comparing multiple forecasts, analyses, persisted analyses and
961climatologies with the same set of observations.
962This approach is referred to as \emph{Class 4} since it is the fourth metric defined by the GODAE intercomparison project. This requires multiple runs of the SAO and running an additional utility (not currently in the \NEMO\ repository) to combine the feedback files into one class 4 file.
963
964%% =================================================================================================
965\section{Observation utilities}
966\label{sec:OBS_obsutils}
967
968For convenience some tools for viewing and processing of observation and feedback files are provided in
969the \NEMO\ repository.
970These tools include OBSTOOLS which are a collection of \fortran\ programs which are helpful to deal with feedback files.
971They do such tasks as observation file conversion, printing of file contents,
972some basic statistical analysis of feedback files.
973The other main tool is an IDL program called dataplot which uses a graphical interface to
974visualise observations and feedback files.
975OBSTOOLS and dataplot are described in more detail below.
976
977%% =================================================================================================
978\subsection{Obstools}
979
980A series of \fortran\ utilities is provided with \NEMO\ called OBSTOOLS.
981This are helpful in handling observation files and the feedback file output from the observation operator. A brief description of some of the utilities follows
982
983%% =================================================================================================
984\subsubsection{corio2fb}
985
986The program corio2fb converts profile observation files from the Coriolis format to the standard feedback format.
987It is called in the following way:
988
989\begin{cmds}
990corio2fb.exe outputfile inputfile1 inputfile2 ...
991\end{cmds}
992
993%% =================================================================================================
994\subsubsection{enact2fb}
995
996The program enact2fb converts profile observation files from the ENACT format to the standard feedback format.
997It is called in the following way:
998
999\begin{cmds}
1000enact2fb.exe outputfile inputfile1 inputfile2 ...
1001\end{cmds}
1002
1003%% =================================================================================================
1004\subsubsection{fbcomb}
1005
1006The program fbcomb combines multiple feedback files produced by individual processors in
1007an MPI run of \NEMO\ into a single feedback file.
1008It is called in the following way:
1009
1010\begin{cmds}
1011fbcomb.exe outputfile inputfile1 inputfile2 ...
1012\end{cmds}
1013
1014%% =================================================================================================
1015\subsubsection{fbmatchup}
1016
1017The program fbmatchup will match observations from two feedback files.
1018It is called in the following way:
1019
1020\begin{cmds}
1021fbmatchup.exe outputfile inputfile1 varname1 inputfile2 varname2 ...
1022\end{cmds}
1023
1024%% =================================================================================================
1025\subsubsection{fbprint}
1026
1027The program fbprint will print the contents of a feedback file or files to standard output.
1028Selected information can be output using optional arguments.
1029It is called in the following way:
1030
1031\begin{cmds}
1032fbprint.exe [options] inputfile
1033
1034options:
1035     -b            shorter output
1036     -q            Select observations based on QC flags
1037     -Q            Select observations based on QC flags
1038     -B            Select observations based on QC flags
1039     -u            unsorted
1040     -s ID         select station ID
1041     -t TYPE       select observation type
1042     -v NUM1-NUM2  select variable range to print by number
1043                      (default all)
1044     -a NUM1-NUM2  select additional variable range to print by number
1045                      (default all)
1046     -e NUM1-NUM2  select extra variable range to print by number
1047                      (default all)
1048     -d            output date range
1049     -D            print depths
1050     -z            use zipped files
1051\end{cmds}
1052
1053%% =================================================================================================
1054\subsubsection{fbsel}
1055
1056The program fbsel will select or subsample observations.
1057It is called in the following way:
1058
1059\begin{cmds}
1060fbsel.exe <input filename> <output filename>
1061\end{cmds}
1062
1063%% =================================================================================================
1064\subsubsection{fbstat}
1065
1066The program fbstat will output summary statistics in different global areas into a number of files.
1067It is called in the following way:
1068
1069\begin{cmds}
1070fbstat.exe [-nmlev] <filenames>
1071\end{cmds}
1072
1073%% =================================================================================================
1074\subsubsection{fbthin}
1075
1076The program fbthin will thin the data to 1 degree resolution.
1077The code could easily be modified to thin to a different resolution.
1078It is called in the following way:
1079
1080\begin{cmds}
1081fbthin.exe inputfile outputfile
1082\end{cmds}
1083
1084%% =================================================================================================
1085\subsubsection{sla2fb}
1086
1087The program sla2fb will convert an AVISO SLA format file to feedback format.
1088It is called in the following way:
1089
1090\begin{cmds}
1091sla2fb.exe [-s type] outputfile inputfile1 inputfile2 ...
1092
1093Option:
1094     -s            Select altimeter data_source
1095\end{cmds}
1096
1097%% =================================================================================================
1098\subsubsection{vel2fb}
1099
1100The program vel2fb will convert TAO/PIRATA/RAMA currents files to feedback format.
1101It is called in the following way:
1102
1103\begin{cmds}
1104vel2fb.exe outputfile inputfile1 inputfile2 ...
1105\end{cmds}
1106
1107%% =================================================================================================
1108\subsection{Building the obstools}
1109
1110To build the obstools use in the tools directory use ./maketools -n OBSTOOLS -m [ARCH].
1111
1112%% =================================================================================================
1113\subsection{Dataplot}
1114
1115An IDL program called dataplot is included which uses a graphical interface to
1116visualise observations and feedback files. Note a similar package has recently developed in python (also called dataplot) which does some of the same things that the IDL dataplot does. Please contact the authors of the this chapter if you are interested in this.
1117
1118It is possible to zoom in, plot individual profiles and calculate some basic statistics.
1119To plot some data run IDL and then:
1120
1121\begin{minted}{idl}
1122IDL> dataplot, "filename"
1123\end{minted}
1124
1125To read multiple files into dataplot,
1126for example multiple feedback files from different processors or from different days,
1127the easiest method is to use the spawn command to generate a list of files which can then be passed to dataplot.
1128
1129\begin{minted}{idl}
1130IDL> spawn, 'ls profb*.nc', files
1131IDL> dataplot, files
1132\end{minted}
1133
1134\autoref{fig:OBS_dataplotmain} shows the main window which is launched when dataplot starts.
1135This is split into three parts.
1136At the top there is a menu bar which contains a variety of drop down menus.
1137Areas - zooms into prespecified regions;
1138plot - plots the data as a timeseries or a T-S diagram if appropriate;
1139Find - allows data to be searched;
1140Config - sets various configuration options.
1141
1142The middle part is a plot of the geographical location of the observations.
1143This will plot the observation value, the model background value or observation minus background value depending on
1144the option selected in the radio button at the bottom of the window.
1145The plotting colour range can be changed by clicking on the colour bar.
1146The title of the plot gives some basic information about the date range and depth range shown,
1147the extreme values, and the mean and RMS values.
1148It is possible to zoom in using a drag-box.
1149You may also zoom in or out using the mouse wheel.
1150
1151The bottom part of the window controls what is visible in the plot above.
1152There are two bars which select the level range plotted (for profile data).
1153The other bars below select the date range shown.
1154The bottom of the figure allows the option to plot the mean, root mean square, standard deviation or
1155mean square values.
1156As mentioned above you can choose to plot the observation value, the model background value or
1157observation minus background value.
1158The next group of radio buttons selects the map projection.
1159This can either be regular longitude latitude grid, or north or south polar stereographic.
1160The next group of radio buttons will plot bad observations, switch to salinity and
1161plot density for profile observations.
1162The rightmost group of buttons will print the plot window as a postscript, save it as png, or exit from dataplot.
1163
1164\begin{figure}
1165  \centering
1166  \includegraphics[width=0.66\textwidth]{OBS_dataplot_main}
1167  \caption{Main window of dataplot}
1168  \label{fig:OBS_dataplotmain}
1169\end{figure}
1170
1171If a profile point is clicked with the mouse button a plot of the observation and background values as
1172a function of depth (\autoref{fig:OBS_dataplotprofile}).
1173
1174\begin{figure}
1175  \centering
1176  \includegraphics[width=0.66\textwidth]{OBS_dataplot_prof}
1177  \caption[Profile plot from dataplot]{
1178    Profile plot from dataplot produced by right clicking on a point in the main window}
1179  \label{fig:OBS_dataplotprofile}
1180\end{figure}
1181
1182\subinc{\input{../../global/epilogue}}
1183
1184\end{document}
Note: See TracBrowser for help on using the repository browser.