source: NEMO/branches/2019/dev_r11233_obsasm_docfixes/doc/latex/NEMO/subfiles/chap_OBS.tex @ 11311

Last change on this file since 11311 was 11311, checked in by djlea, 15 months ago

#2297 Fix a few typos

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