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_OBS.tex in NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/doc/latex/NEMO/subfiles – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/doc/latex/NEMO/subfiles/chap_OBS.tex @ 14644

Last change on this file since 14644 was 14644, checked in by sparonuz, 3 years ago

Merge trunk -r14642:HEAD

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