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

Last change on this file since 11577 was 11577, checked in by nicolasmartin, 13 months ago

New LaTeX commands \nam and \np to mention namelist content
(Partial commit to serve as a backup before other large edits)
In order to benefit of the syntax highlighting and to have a simpler syntax for
citing namelist block (\nam) and parameter (\np) with an optional variable assignment (\forcode{…}),
at this time the only viable solution I found is to require a double marker for
what it looks like the same item:

  1. Marker with the real name: 'tra_adv' block or 'ln_flx' parameter
  2. Marker with underscore character escaping: 'tra\_adv' block or 'ln\_flx' parameter

Despite many searches and attempts, I did not find a workaround to edit on-the-fly one or
the other marker.
In fact, the problem is on one side that the LaTeX index interprets '_' as a switch for lowering like
in math mode while on the other hand the backslash is considered for Pygments as a typo in Fortran
(red box).

For instance, \nam and \np have as of now the aforementioned 2 mandatory arguments in
the previous order (between braces) + an optional argument for \np when the parameter is defined
(between brackets at the first position):

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