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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Chap_OBS.tex in branches/2017/dev_r8657_UKMO_OBSoper/DOC/TexFiles/Chapters – NEMO

source: branches/2017/dev_r8657_UKMO_OBSoper/DOC/TexFiles/Chapters/Chap_OBS.tex @ 8921

Last change on this file since 8921 was 8921, checked in by mattmartin, 6 years ago

Updated documentation after testing the compilation of the latex.
File size: 63.1 KB
Line 
1\documentclass[NEMO_book]{subfiles}
2\begin{document}
3% ================================================================
4% Chapter observation operator (OBS)
5% ================================================================
6\chapter{Observation and model comparison (OBS)}
7\label{OBS}
8
9Authors: D. Lea, M. Martin, K. Mogensen, A. Vidard, A. Weaver, A. Ryan, ...   % do we keep that ?
10
11\minitoc
12
13
14\newpage
15$\ $\newline    % force a new line
16
17The observation and model comparison code (OBS) reads in observation files (profile
18temperature and salinity, sea surface temperature, sea level anomaly, sea ice concentration,
19and velocity) and calculates  an interpolated model equivalent value at the observation
20location and nearest model timestep. The resulting data are saved in a ``feedback'' file (or
21files). The code was originally developed for use with the NEMOVAR data assimilation code, but
22can be used for validation or verification of the model or with any other data assimilation system.
23
24The OBS code is called from \mdl{nemogcm} for model initialisation and to calculate the model
25equivalent values for observations on the 0th timestep. The code is then called again after
26each timestep from \mdl{step}. The code is only activated if the namelist logical \np{ln\_diaobs}
27is set to true.
28
29For all data types a 2D horizontal interpolator or averager is needed to interpolate/average the model fields to
30the observation location. For {\em in situ} profiles, a 1D vertical interpolator is needed in
31addition to provide model fields at the observation depths. This now works in a generalised vertical
32coordinate system.
33
34Some profile observation types (e.g. tropical moored buoys) are made available as daily averaged quantities.
35The observation operator code can be set-up to calculate the equivalent daily average model temperature fields
36using the \np{nn\_profdavtypes} namelist array. Some SST observations are equivalent to a night-time
37average value and the observation operator code can calculate equivalent night-time average model SST fields by
38setting the namelist value \np{ln\_sstnight} to true. Otherwise the model value from the nearest timestep to the
39observation time is used.
40
41The code is controlled by the namelist \textit{namobs}. See the following sections for more
42details on setting up the namelist.
43
44Section~\ref{OBS_example} introduces a test example of the observation operator code including
45where to obtain data and how to setup the namelist. Section~\ref{OBS_details} introduces some
46more technical details of the different observation types used and also shows a more complete
47namelist. Section~\ref{OBS_theory} introduces some of the theoretical aspects of the observation
48operator including interpolation methods and running on multiple processors.
49Section~\ref{OBS_ooo} describes the offline observation operator code.
50Section~\ref{OBS_obsutils} introduces some utilities to help working with the files
51produced by the OBS code.
52
53% ================================================================
54% Example
55% ================================================================
56\section{Running the observation operator code example}
57\label{OBS_example}
58
59This section describes an example of running the observation operator code using
60profile data which can be freely downloaded. It shows how to adapt an
61existing run and build of NEMO to run the observation operator.
62
63\begin{enumerate}
64\item Compile NEMO.
65
66\item Download some EN4 data from
67\href{http://www.metoffice.gov.uk/hadobs}{http://www.metoffice.gov.uk/hadobs}. Choose observations which are
68valid for the period of your test run because the observation operator compares
69the model and observations for a matching date and time.
70
71\item Compile the OBSTOOLS code using:
72\begin{verbatim}
73./maketools -n OBSTOOLS -m [ARCH].
74\end{verbatim}
75
76\item Convert the EN4 data into feedback format:
77\begin{verbatim}
78enact2fb.exe profiles_01.nc EN.4.1.1.f.profiles.g10.YYYYMM.nc
79\end{verbatim}
80
81\item Include the following in the NEMO namelist to run the observation
82operator on this data:
83\end{enumerate}
84
85%------------------------------------------namobs_example-----------------------------------------------------
86\namdisplay{namobs_example}
87%-------------------------------------------------------------------------------------------------------------
88
89Options are defined through the  \ngn{namobs} namelist variables.
90The options \np{ln\_t3d} and \np{ln\_s3d} switch on the temperature and salinity
91profile observation operator code. The filename or array of filenames are
92specified using the \np{cn\_profbfiles} variable. The model grid points for a
93particular  observation latitude and longitude are found using the grid
94searching part of the code. This can be expensive, particularly for large
95numbers of observations, setting \np{ln\_grid\_search\_lookup} allows the use of
96a lookup table which is saved into an ``xypos`` file (or files). This will need
97to be generated the first time if it does not exist in the run directory.
98However, once produced it will significantly speed up future grid searches.
99Setting \np{ln\_grid\_global} means that the code distributes the observations
100evenly between processors. Alternatively each processor will work with
101observations located within the model subdomain (see section~\ref{OBS_parallel}).
102
103A number of utilities are now provided to plot the feedback files, convert and
104recombine the files. These are explained in more detail in section~\ref{OBS_obsutils}.
105Utilites to convert other input data formats into the feedback format are also
106described in section~\ref{OBS_obsutils}.
107
108\section{Technical details}
109\label{OBS_details}
110
111Here we show a more complete example namelist  \ngn{namobs} and also show the NetCDF headers
112of the observation
113files that may be used with the observation operator
114
115%------------------------------------------namobs--------------------------------------------------------
116\namdisplay{namobs}
117%-------------------------------------------------------------------------------------------------------------
118
119The observation operator code uses the "feedback" observation file format for
120all data types. All the
121observation files must be in NetCDF format. Some example headers (produced using
122\mbox{\textit{ncdump~-h}}) for profile
123data, sea level anomaly and sea surface temperature are in the following
124subsections.
125
126\subsection{Profile feedback type observation file header}
127
128\begin{alltt}
129\tiny
130\begin{verbatim}
131netcdf profiles_01 {
132dimensions:
133     N_OBS = 603 ;
134     N_LEVELS = 150 ;
135     N_VARS = 2 ;
136     N_QCF = 2 ;
137     N_ENTRIES = 1 ;
138     N_EXTRA = 1 ;
139     STRINGNAM = 8 ;
140     STRINGGRID = 1 ;
141     STRINGWMO = 8 ;
142     STRINGTYP = 4 ;
143     STRINGJULD = 14 ;
144variables:
145     char VARIABLES(N_VARS, STRINGNAM) ;
146          VARIABLES:long_name = "List of variables in feedback files" ;
147     char ENTRIES(N_ENTRIES, STRINGNAM) ;
148          ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
149     char EXTRA(N_EXTRA, STRINGNAM) ;
150          EXTRA:long_name = "List of extra variables" ;
151     char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
152          STATION_IDENTIFIER:long_name = "Station identifier" ;
153     char STATION_TYPE(N_OBS, STRINGTYP) ;
154          STATION_TYPE:long_name = "Code instrument type" ;
155     double LONGITUDE(N_OBS) ;
156          LONGITUDE:long_name = "Longitude" ;
157          LONGITUDE:units = "degrees_east" ;
158          LONGITUDE:_Fillvalue = 99999.f ;
159     double LATITUDE(N_OBS) ;
160          LATITUDE:long_name = "Latitude" ;
161          LATITUDE:units = "degrees_north" ;
162          LATITUDE:_Fillvalue = 99999.f ;
163     double DEPTH(N_OBS, N_LEVELS) ;
164          DEPTH:long_name = "Depth" ;
165          DEPTH:units = "metre" ;
166          DEPTH:_Fillvalue = 99999.f ;
167     int DEPTH_QC(N_OBS, N_LEVELS) ;
168          DEPTH_QC:long_name = "Quality on depth" ;
169          DEPTH_QC:Conventions = "q where q =[0,9]" ;
170          DEPTH_QC:_Fillvalue = 0 ;
171     int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
172          DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
173          DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
174     double JULD(N_OBS) ;
175          JULD:long_name = "Julian day" ;
176          JULD:units = "days since JULD_REFERENCE" ;
177          JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
178          JULD:_Fillvalue = 99999.f ;
179     char JULD_REFERENCE(STRINGJULD) ;
180          JULD_REFERENCE:long_name = "Date of reference for julian days" ;
181          JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
182     int OBSERVATION_QC(N_OBS) ;
183          OBSERVATION_QC:long_name = "Quality on observation" ;
184          OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
185          OBSERVATION_QC:_Fillvalue = 0 ;
186     int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
187          OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
188          OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
189          OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
190     int POSITION_QC(N_OBS) ;
191          POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
192          POSITION_QC:Conventions = "q where q =[0,9]" ;
193          POSITION_QC:_Fillvalue = 0 ;
194     int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
195          POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
196          POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
197          POSITION_QC_FLAGS:_Fillvalue = 0 ;
198     int JULD_QC(N_OBS) ;
199          JULD_QC:long_name = "Quality on date and time" ;
200          JULD_QC:Conventions = "q where q =[0,9]" ;
201          JULD_QC:_Fillvalue = 0 ;
202     int JULD_QC_FLAGS(N_OBS, N_QCF) ;
203          JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
204          JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
205          JULD_QC_FLAGS:_Fillvalue = 0 ;
206     int ORIGINAL_FILE_INDEX(N_OBS) ;
207          ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
208          ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
209     float POTM_OBS(N_OBS, N_LEVELS) ;
210          POTM_OBS:long_name = "Potential temperature" ;
211          POTM_OBS:units = "Degrees Celsius" ;
212          POTM_OBS:_Fillvalue = 99999.f ;
213     float POTM_Hx(N_OBS, N_LEVELS) ;
214          POTM_Hx:long_name = "Model interpolated potential temperature" ;
215          POTM_Hx:units = "Degrees Celsius" ;
216          POTM_Hx:_Fillvalue = 99999.f ;
217     int POTM_QC(N_OBS) ;
218          POTM_QC:long_name = "Quality on potential temperature" ;
219          POTM_QC:Conventions = "q where q =[0,9]" ;
220          POTM_QC:_Fillvalue = 0 ;
221     int POTM_QC_FLAGS(N_OBS, N_QCF) ;
222          POTM_QC_FLAGS:long_name = "Quality flags on potential temperature" ;
223          POTM_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
224          POTM_QC_FLAGS:_Fillvalue = 0 ;
225     int POTM_LEVEL_QC(N_OBS, N_LEVELS) ;
226          POTM_LEVEL_QC:long_name = "Quality for each level on potential temperature" ;
227          POTM_LEVEL_QC:Conventions = "q where q =[0,9]" ;
228          POTM_LEVEL_QC:_Fillvalue = 0 ;
229     int POTM_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
230          POTM_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on potential temperature" ;
231          POTM_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
232          POTM_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
233     int POTM_IOBSI(N_OBS) ;
234          POTM_IOBSI:long_name = "ORCA grid search I coordinate" ;
235     int POTM_IOBSJ(N_OBS) ;
236          POTM_IOBSJ:long_name = "ORCA grid search J coordinate" ;
237     int POTM_IOBSK(N_OBS, N_LEVELS) ;
238          POTM_IOBSK:long_name = "ORCA grid search K coordinate" ;
239     char POTM_GRID(STRINGGRID) ;
240          POTM_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
241     float PSAL_OBS(N_OBS, N_LEVELS) ;
242          PSAL_OBS:long_name = "Practical salinity" ;
243          PSAL_OBS:units = "PSU" ;
244          PSAL_OBS:_Fillvalue = 99999.f ;
245     float PSAL_Hx(N_OBS, N_LEVELS) ;
246          PSAL_Hx:long_name = "Model interpolated practical salinity" ;
247          PSAL_Hx:units = "PSU" ;
248          PSAL_Hx:_Fillvalue = 99999.f ;
249     int PSAL_QC(N_OBS) ;
250          PSAL_QC:long_name = "Quality on practical salinity" ;
251          PSAL_QC:Conventions = "q where q =[0,9]" ;
252          PSAL_QC:_Fillvalue = 0 ;
253     int PSAL_QC_FLAGS(N_OBS, N_QCF) ;
254          PSAL_QC_FLAGS:long_name = "Quality flags on practical salinity" ;
255          PSAL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
256          PSAL_QC_FLAGS:_Fillvalue = 0 ;
257     int PSAL_LEVEL_QC(N_OBS, N_LEVELS) ;
258          PSAL_LEVEL_QC:long_name = "Quality for each level on practical salinity" ;
259          PSAL_LEVEL_QC:Conventions = "q where q =[0,9]" ;
260          PSAL_LEVEL_QC:_Fillvalue = 0 ;
261     int PSAL_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
262          PSAL_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on practical salinity" ;
263          PSAL_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
264          PSAL_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
265     int PSAL_IOBSI(N_OBS) ;
266          PSAL_IOBSI:long_name = "ORCA grid search I coordinate" ;
267     int PSAL_IOBSJ(N_OBS) ;
268          PSAL_IOBSJ:long_name = "ORCA grid search J coordinate" ;
269     int PSAL_IOBSK(N_OBS, N_LEVELS) ;
270          PSAL_IOBSK:long_name = "ORCA grid search K coordinate" ;
271     char PSAL_GRID(STRINGGRID) ;
272          PSAL_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
273     float TEMP(N_OBS, N_LEVELS) ;
274          TEMP:long_name = "Insitu temperature" ;
275          TEMP:units = "Degrees Celsius" ;
276          TEMP:_Fillvalue = 99999.f ;
277
278// global attributes:
279          :title = "NEMO observation operator output" ;
280          :Convention = "NEMO unified observation operator output" ;
281}
282\end{verbatim}
283\end{alltt}
284
285\subsection{Sea level anomaly feedback type observation file header}
286
287\begin{alltt}
288\tiny
289\begin{verbatim}
290netcdf sla_01 {
291dimensions:
292     N_OBS = 41301 ;
293     N_LEVELS = 1 ;
294     N_VARS = 1 ;
295     N_QCF = 2 ;
296     N_ENTRIES = 1 ;
297     N_EXTRA = 1 ;
298     STRINGNAM = 8 ;
299     STRINGGRID = 1 ;
300     STRINGWMO = 8 ;
301     STRINGTYP = 4 ;
302     STRINGJULD = 14 ;
303variables:
304     char VARIABLES(N_VARS, STRINGNAM) ;
305          VARIABLES:long_name = "List of variables in feedback files" ;
306     char ENTRIES(N_ENTRIES, STRINGNAM) ;
307          ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
308     char EXTRA(N_EXTRA, STRINGNAM) ;
309          EXTRA:long_name = "List of extra variables" ;
310     char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
311          STATION_IDENTIFIER:long_name = "Station identifier" ;
312     char STATION_TYPE(N_OBS, STRINGTYP) ;
313          STATION_TYPE:long_name = "Code instrument type" ;
314     double LONGITUDE(N_OBS) ;
315          LONGITUDE:long_name = "Longitude" ;
316          LONGITUDE:units = "degrees_east" ;
317          LONGITUDE:_Fillvalue = 99999.f ;
318     double LATITUDE(N_OBS) ;
319          LATITUDE:long_name = "Latitude" ;
320          LATITUDE:units = "degrees_north" ;
321          LATITUDE:_Fillvalue = 99999.f ;
322     double DEPTH(N_OBS, N_LEVELS) ;
323          DEPTH:long_name = "Depth" ;
324          DEPTH:units = "metre" ;
325          DEPTH:_Fillvalue = 99999.f ;
326     int DEPTH_QC(N_OBS, N_LEVELS) ;
327          DEPTH_QC:long_name = "Quality on depth" ;
328          DEPTH_QC:Conventions = "q where q =[0,9]" ;
329          DEPTH_QC:_Fillvalue = 0 ;
330     int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
331          DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
332          DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
333     double JULD(N_OBS) ;
334          JULD:long_name = "Julian day" ;
335          JULD:units = "days since JULD_REFERENCE" ;
336          JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
337          JULD:_Fillvalue = 99999.f ;
338     char JULD_REFERENCE(STRINGJULD) ;
339          JULD_REFERENCE:long_name = "Date of reference for julian days" ;
340          JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
341     int OBSERVATION_QC(N_OBS) ;
342          OBSERVATION_QC:long_name = "Quality on observation" ;
343          OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
344          OBSERVATION_QC:_Fillvalue = 0 ;
345     int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
346          OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
347          OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
348          OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
349     int POSITION_QC(N_OBS) ;
350          POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
351          POSITION_QC:Conventions = "q where q =[0,9]" ;
352          POSITION_QC:_Fillvalue = 0 ;
353     int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
354          POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
355          POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
356          POSITION_QC_FLAGS:_Fillvalue = 0 ;
357     int JULD_QC(N_OBS) ;
358          JULD_QC:long_name = "Quality on date and time" ;
359          JULD_QC:Conventions = "q where q =[0,9]" ;
360          JULD_QC:_Fillvalue = 0 ;
361     int JULD_QC_FLAGS(N_OBS, N_QCF) ;
362          JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
363          JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
364          JULD_QC_FLAGS:_Fillvalue = 0 ;
365     int ORIGINAL_FILE_INDEX(N_OBS) ;
366          ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
367          ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
368     float SLA_OBS(N_OBS, N_LEVELS) ;
369          SLA_OBS:long_name = "Sea level anomaly" ;
370          SLA_OBS:units = "metre" ;
371          SLA_OBS:_Fillvalue = 99999.f ;
372     float SLA_Hx(N_OBS, N_LEVELS) ;
373          SLA_Hx:long_name = "Model interpolated sea level anomaly" ;
374          SLA_Hx:units = "metre" ;
375          SLA_Hx:_Fillvalue = 99999.f ;
376     int SLA_QC(N_OBS) ;
377          SLA_QC:long_name = "Quality on sea level anomaly" ;
378          SLA_QC:Conventions = "q where q =[0,9]" ;
379          SLA_QC:_Fillvalue = 0 ;
380     int SLA_QC_FLAGS(N_OBS, N_QCF) ;
381          SLA_QC_FLAGS:long_name = "Quality flags on sea level anomaly" ;
382          SLA_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
383          SLA_QC_FLAGS:_Fillvalue = 0 ;
384     int SLA_LEVEL_QC(N_OBS, N_LEVELS) ;
385          SLA_LEVEL_QC:long_name = "Quality for each level on sea level anomaly" ;
386          SLA_LEVEL_QC:Conventions = "q where q =[0,9]" ;
387          SLA_LEVEL_QC:_Fillvalue = 0 ;
388     int SLA_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
389          SLA_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on sea level anomaly" ;
390          SLA_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
391          SLA_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
392     int SLA_IOBSI(N_OBS) ;
393          SLA_IOBSI:long_name = "ORCA grid search I coordinate" ;
394     int SLA_IOBSJ(N_OBS) ;
395          SLA_IOBSJ:long_name = "ORCA grid search J coordinate" ;
396     int SLA_IOBSK(N_OBS, N_LEVELS) ;
397          SLA_IOBSK:long_name = "ORCA grid search K coordinate" ;
398     char SLA_GRID(STRINGGRID) ;
399          SLA_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
400     float MDT(N_OBS, N_LEVELS) ;
401          MDT:long_name = "Mean Dynamic Topography" ;
402          MDT:units = "metre" ;
403          MDT:_Fillvalue = 99999.f ;
404
405// global attributes:
406          :title = "NEMO observation operator output" ;
407          :Convention = "NEMO unified observation operator output" ;
408}
409\end{verbatim}
410\end{alltt}
411
412The mean dynamic
413topography (MDT) must be provided in a separate file defined on the model grid
414 called {\it slaReferenceLevel.nc}. The MDT is required in
415order to produce the model equivalent sea level anomaly from the model sea
416surface height. Below is an example header for this file (on the ORCA025 grid).
417
418\begin{alltt}
419\tiny
420\begin{verbatim}
421dimensions:
422        x = 1442 ;
423        y = 1021 ;
424variables:
425        float nav_lon(y, x) ;
426                nav_lon:units = "degrees_east" ;
427        float nav_lat(y, x) ;
428                nav_lat:units = "degrees_north" ;
429        float sossheig(y, x) ;
430                sossheig:_FillValue = -1.e+30f ;
431                sossheig:coordinates = "nav_lon nav_lat" ;
432                sossheig:long_name = "Mean Dynamic Topography" ;
433                sossheig:units = "metres" ;
434                sossheig:grid = "orca025T" ;
435\end{verbatim}
436\end{alltt}
437
438\subsection{Sea surface temperature feedback type observation file header}
439
440\begin{alltt}
441\tiny
442\begin{verbatim}
443netcdf sst_01 {
444dimensions:
445     N_OBS = 33099 ;
446     N_LEVELS = 1 ;
447     N_VARS = 1 ;
448     N_QCF = 2 ;
449     N_ENTRIES = 1 ;
450     STRINGNAM = 8 ;
451     STRINGGRID = 1 ;
452     STRINGWMO = 8 ;
453     STRINGTYP = 4 ;
454     STRINGJULD = 14 ;
455variables:
456     char VARIABLES(N_VARS, STRINGNAM) ;
457          VARIABLES:long_name = "List of variables in feedback files" ;
458     char ENTRIES(N_ENTRIES, STRINGNAM) ;
459          ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
460     char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
461          STATION_IDENTIFIER:long_name = "Station identifier" ;
462     char STATION_TYPE(N_OBS, STRINGTYP) ;
463          STATION_TYPE:long_name = "Code instrument type" ;
464     double LONGITUDE(N_OBS) ;
465          LONGITUDE:long_name = "Longitude" ;
466          LONGITUDE:units = "degrees_east" ;
467          LONGITUDE:_Fillvalue = 99999.f ;
468     double LATITUDE(N_OBS) ;
469          LATITUDE:long_name = "Latitude" ;
470          LATITUDE:units = "degrees_north" ;
471          LATITUDE:_Fillvalue = 99999.f ;
472     double DEPTH(N_OBS, N_LEVELS) ;
473          DEPTH:long_name = "Depth" ;
474          DEPTH:units = "metre" ;
475          DEPTH:_Fillvalue = 99999.f ;
476     int DEPTH_QC(N_OBS, N_LEVELS) ;
477          DEPTH_QC:long_name = "Quality on depth" ;
478          DEPTH_QC:Conventions = "q where q =[0,9]" ;
479          DEPTH_QC:_Fillvalue = 0 ;
480     int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
481          DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
482          DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
483     double JULD(N_OBS) ;
484          JULD:long_name = "Julian day" ;
485          JULD:units = "days since JULD_REFERENCE" ;
486          JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
487          JULD:_Fillvalue = 99999.f ;
488     char JULD_REFERENCE(STRINGJULD) ;
489          JULD_REFERENCE:long_name = "Date of reference for julian days" ;
490          JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
491     int OBSERVATION_QC(N_OBS) ;
492          OBSERVATION_QC:long_name = "Quality on observation" ;
493          OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
494          OBSERVATION_QC:_Fillvalue = 0 ;
495     int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
496          OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
497          OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
498          OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
499     int POSITION_QC(N_OBS) ;
500          POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
501          POSITION_QC:Conventions = "q where q =[0,9]" ;
502          POSITION_QC:_Fillvalue = 0 ;
503     int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
504          POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
505          POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
506          POSITION_QC_FLAGS:_Fillvalue = 0 ;
507     int JULD_QC(N_OBS) ;
508          JULD_QC:long_name = "Quality on date and time" ;
509          JULD_QC:Conventions = "q where q =[0,9]" ;
510          JULD_QC:_Fillvalue = 0 ;
511     int JULD_QC_FLAGS(N_OBS, N_QCF) ;
512          JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
513          JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
514          JULD_QC_FLAGS:_Fillvalue = 0 ;
515     int ORIGINAL_FILE_INDEX(N_OBS) ;
516          ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
517          ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
518     float SST_OBS(N_OBS, N_LEVELS) ;
519          SST_OBS:long_name = "Sea surface temperature" ;
520          SST_OBS:units = "Degree centigrade" ;
521          SST_OBS:_Fillvalue = 99999.f ;
522     float SST_Hx(N_OBS, N_LEVELS) ;
523          SST_Hx:long_name = "Model interpolated sea surface temperature" ;
524          SST_Hx:units = "Degree centigrade" ;
525          SST_Hx:_Fillvalue = 99999.f ;
526     int SST_QC(N_OBS) ;
527          SST_QC:long_name = "Quality on sea surface temperature" ;
528          SST_QC:Conventions = "q where q =[0,9]" ;
529          SST_QC:_Fillvalue = 0 ;
530     int SST_QC_FLAGS(N_OBS, N_QCF) ;
531          SST_QC_FLAGS:long_name = "Quality flags on sea surface temperature" ;
532          SST_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
533          SST_QC_FLAGS:_Fillvalue = 0 ;
534     int SST_LEVEL_QC(N_OBS, N_LEVELS) ;
535          SST_LEVEL_QC:long_name = "Quality for each level on sea surface temperature" ;
536          SST_LEVEL_QC:Conventions = "q where q =[0,9]" ;
537          SST_LEVEL_QC:_Fillvalue = 0 ;
538     int SST_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
539          SST_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on sea surface temperature" ;
540          SST_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
541          SST_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
542     int SST_IOBSI(N_OBS) ;
543          SST_IOBSI:long_name = "ORCA grid search I coordinate" ;
544     int SST_IOBSJ(N_OBS) ;
545          SST_IOBSJ:long_name = "ORCA grid search J coordinate" ;
546     int SST_IOBSK(N_OBS, N_LEVELS) ;
547          SST_IOBSK:long_name = "ORCA grid search K coordinate" ;
548     char SST_GRID(STRINGGRID) ;
549          SST_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
550
551// global attributes:
552          :title = "NEMO observation operator output" ;
553          :Convention = "NEMO unified observation operator output" ;
554}
555\end{verbatim}
556\end{alltt}
557
558\section{Theoretical details}
559\label{OBS_theory}
560
561\subsection{Horizontal interpolation and averaging methods}
562
563For most observation types, the horizontal extent of the observation is small compared to the model grid size
564and so the model equivalent of the observation is calculated by interpolating from the four surrounding grid
565points to the observation location. Some satellite observations (e.g. microwave satellite SST data, or satellite SSS data)
566have a footprint which is similar in size or larger than the model grid size (particularly when the grid size is small).
567In those cases the model counterpart should be calculated by averaging the model grid points over the same size as the footprint.
568NEMO therefore has the capability to specify either an interpolation or an averaging (for surface observation types only).
569
570The main namelist option associated with the interpolation/averaging is \np{nn\_2dint}. This default option can be set to values from 0 to 6.
571Values between 0 to 4 are associated with interpolation while values 5 or 6 are associated with averaging.
572\begin{itemize}
573\item \np{nn\_2dint}=0: Distance-weighted interpolation
574\item \np{nn\_2dint}=1: Distance-weighted interpolation (small angle)
575\item \np{nn\_2dint}=2: Bilinear interpolation (geographical grid)
576\item \np{nn\_2dint}=3: Bilinear remapping interpolation (general grid)
577\item \np{nn\_2dint}=4: Polynomial interpolation
578\item \np{nn\_2dint}=5: Radial footprint averaging with diameter specified in the namelist as \np{rn\_???\_avglamscl} in degrees or metres (set using \np{ln\_???\_fp\_indegs})
579\item \np{nn\_2dint}=6: Rectangular footprint averaging with E/W and N/S size specified in the namelist as \np{rn\_???\_avglamscl} and \np{rn\_???\_avgphiscl} in degrees or metres (set using \np{ln\_???\_fp\_indegs})
580\end{itemize}
581The ??? in the last two options indicate these options should be specified for each observation type for which the averaging is to be performed (see namelist example above).
582The \np{nn\_2dint} default option can be overridden for surface observation types using namelist values \np{nn\_2dint\_???} where ??? is one of sla,sst,sss,sic.
583
584Below is some more detail on the various options for interpolation and averaging available in NEMO.
585
586\subsubsection{Horizontal interpolation}
587Consider an observation point ${\rm P}$ with
588with longitude and latitude $({\lambda_{}}_{\rm P}, \phi_{\rm P})$ and the
589four nearest neighbouring model grid points ${\rm A}$, ${\rm B}$, ${\rm C}$ 
590and ${\rm D}$ with longitude and latitude ($\lambda_{\rm A}$, $\phi_{\rm A}$),
591($\lambda_{\rm B}$, $\phi_{\rm B}$) etc.
592All horizontal interpolation methods implemented in NEMO
593estimate the value of a model variable $x$ at point $P$ as
594a weighted linear combination of the values of the model
595variables at the grid points ${\rm A}$, ${\rm B}$ etc.:
596\begin{eqnarray}
597{x_{}}_{\rm P} & \hspace{-2mm} = \hspace{-2mm} & 
598\frac{1}{w} \left( {w_{}}_{\rm A} {x_{}}_{\rm A} +
599                   {w_{}}_{\rm B} {x_{}}_{\rm B} +
600                   {w_{}}_{\rm C} {x_{}}_{\rm C} +
601                   {w_{}}_{\rm D} {x_{}}_{\rm D} \right)
602\end{eqnarray}
603where ${w_{}}_{\rm A}$, ${w_{}}_{\rm B}$ etc. are the respective weights for the
604model field at points ${\rm A}$, ${\rm B}$ etc., and
605$w = {w_{}}_{\rm A} + {w_{}}_{\rm B} + {w_{}}_{\rm C} + {w_{}}_{\rm D}$.
606
607Four different possibilities are available for computing the weights.
608
609\begin{enumerate}
610
611\item[1.] {\bf Great-Circle distance-weighted interpolation.} The weights
612  are computed as a function of the great-circle distance $s(P, \cdot)$ 
613  between $P$ and the model grid points $A$, $B$ etc. For example,
614  the weight given to the field ${x_{}}_{\rm A}$ is specified as the
615  product of the distances from ${\rm P}$ to the other points:
616  \begin{eqnarray}
617  {w_{}}_{\rm A} = s({\rm P}, {\rm B}) \, s({\rm P}, {\rm C}) \, s({\rm P}, {\rm D})
618  \nonumber
619  \end{eqnarray}
620  where
621  \begin{eqnarray}
622   s\left ({\rm P}, {\rm M} \right )
623     & \hspace{-2mm} = \hspace{-2mm} & 
624      \cos^{-1} \! \left\{ 
625               \sin {\phi_{}}_{\rm P} \sin {\phi_{}}_{\rm M}
626             + \cos {\phi_{}}_{\rm P} \cos {\phi_{}}_{\rm M} 
627               \cos ({\lambda_{}}_{\rm M} - {\lambda_{}}_{\rm P})
628                   \right\}
629   \end{eqnarray}
630   and $M$ corresponds to $B$, $C$ or $D$.
631   A more stable form of the great-circle distance formula for
632   small distances ($x$ near 1) involves the arcsine function
633   ($e.g.$ see p.~101 of \citet{Daley_Barker_Bk01}:
634   \begin{eqnarray}
635   s\left( {\rm P}, {\rm M} \right)
636     & \hspace{-2mm} = \hspace{-2mm} & 
637      \sin^{-1} \! \left\{ \sqrt{ 1 - x^2 } \right\}
638   \nonumber
639   \end{eqnarray}
640   where
641   \begin{eqnarray}
642    x & \hspace{-2mm} = \hspace{-2mm} & 
643      {a_{}}_{\rm M} {a_{}}_{\rm P} + {b_{}}_{\rm M} {b_{}}_{\rm P} + {c_{}}_{\rm M} {c_{}}_{\rm P}
644   \nonumber
645   \end{eqnarray}
646   and
647   \begin{eqnarray}
648      {a_{}}_{\rm M} & \hspace{-2mm} = \hspace{-2mm} & \sin {\phi_{}}_{\rm M},
649      \nonumber \\
650      {a_{}}_{\rm P} & \hspace{-2mm} = \hspace{-2mm} & \sin {\phi_{}}_{\rm P},
651      \nonumber \\
652      {b_{}}_{\rm M} & \hspace{-2mm} = \hspace{-2mm} & \cos {\phi_{}}_{\rm M} \cos {\phi_{}}_{\rm M},
653      \nonumber \\
654      {b_{}}_{\rm P} & \hspace{-2mm} = \hspace{-2mm} & \cos {\phi_{}}_{\rm P} \cos {\phi_{}}_{\rm P},
655      \nonumber \\
656      {c_{}}_{\rm M} & \hspace{-2mm} = \hspace{-2mm} & \cos {\phi_{}}_{\rm M} \sin {\phi_{}}_{\rm M},
657      \nonumber \\
658      {c_{}}_{\rm P} & \hspace{-2mm} = \hspace{-2mm} & \cos {\phi_{}}_{\rm P} \sin {\phi_{}}_{\rm P}.
659      \nonumber
660   \nonumber
661  \end{eqnarray}
662
663\item[2.] {\bf Great-Circle distance-weighted interpolation with small angle
664  approximation.} Similar to the previous interpolation but with the
665  distance $s$ computed as
666  \begin{eqnarray}
667    s\left( {\rm P}, {\rm M} \right)
668     & \hspace{-2mm} = \hspace{-2mm} & 
669      \sqrt{ \left( {\phi_{}}_{\rm M} - {\phi_{}}_{\rm P} \right)^{2} 
670      + \left( {\lambda_{}}_{\rm M} - {\lambda_{}}_{\rm P} \right)^{2}
671        \cos^{2} {\phi_{}}_{\rm M} }
672  \end{eqnarray}
673  where $M$ corresponds to $A$, $B$, $C$ or $D$.
674
675\item[3.] {\bf Bilinear interpolation for a regular spaced grid.} The
676  interpolation is split into two 1D interpolations in the longitude
677  and latitude directions, respectively.
678
679\item[4.] {\bf Bilinear remapping interpolation for a general grid.} An
680  iterative scheme that involves first mapping a quadrilateral cell
681  into a cell with coordinates (0,0), (1,0), (0,1) and (1,1). This
682  method is based on the SCRIP interpolation package \citep{Jones_1998}.
683 
684\end{enumerate}
685
686\subsubsection{Horizontal averaging}
687
688For each surface observation type:
689\begin{itemize}
690\item The standard grid-searching code is used to find the nearest model grid point to the observation location (see next subsection).
691\item The maximum number of grid points is calculated in the local grid domain for which the averaging is likely need to cover.
692\item The lats/longs of the grid points surrounding the nearest model grid box are extracted using existing mpi routines.
693\item The weights for each grid point associated with each observation are calculated, either for radial or rectangular footprints. For grid points completely within the footprint, the weight is one; for grid points completely outside the footprint, the weight is zero. For grid points which are partly within the footprint the ratio between the area of the footprint within the grid box and the total area of the grid box is used as the weight.
694\item The weighted average of the model grid points associated with each observation is calculated, and this is then given as the model counterpart of the observation.
695\end{itemize}
696
697Examples of the weights calculated for an observation with rectangular and radial footprints are shown in Figs.~\ref{fig:obsavgrec} and~\ref{fig:obsavgrad}.
698
699%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
700\begin{figure}      \begin{center}
701\includegraphics[width=0.90\textwidth]{Fig_OBS_avg_rec}
702\caption{      \label{fig:obsavgrec}
703Weights associated with each model grid box (blue lines and numbers) for an observation at -170.5E, 56.0N with a rectangular footprint of 1\deg x 1\deg.}
704\end{center}      \end{figure}
705%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
706
707%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
708\begin{figure}      \begin{center}
709\includegraphics[width=0.90\textwidth]{Fig_OBS_avg_rad}
710\caption{      \label{fig:obsavgrad}
711Weights associated with each model grid box (blue lines and numbers) for an observation at -170.5E, 56.0N with a radial footprint with diameter 1\deg.} 
712\end{center}      \end{figure}
713%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
714
715
716\subsection{Grid search}
717
718For many grids used by the NEMO model, such as the ORCA family,
719the horizontal grid coordinates $i$ and $j$ are not simple functions
720of latitude and longitude. Therefore, it is not always straightforward
721to determine the grid points surrounding any given observational position.
722Before the interpolation can be performed, a search
723algorithm is then required to determine the corner points of
724the quadrilateral cell in which the observation is located.
725This is the most difficult and time consuming part of the
7262D interpolation procedure.
727A robust test for determining if an observation falls
728within a given quadrilateral cell is as follows. Let
729${\rm P}({\lambda_{}}_{\rm P} ,{\phi_{}}_{\rm P} )$ denote the observation point,
730and let ${\rm A}({\lambda_{}}_{\rm A} ,{\phi_{}}_{\rm A} )$,
731${\rm B}({\lambda_{}}_{\rm B} ,{\phi_{}}_{\rm B} )$,
732${\rm C}({\lambda_{}}_{\rm C} ,{\phi_{}}_{\rm C} )$ 
733and
734${\rm D}({\lambda_{}}_{\rm D} ,{\phi_{}}_{\rm D} )$ denote
735the bottom left, bottom right, top left and top right
736corner points of the cell, respectively.
737To determine if P is inside
738the cell, we verify that the cross-products
739\begin{eqnarray}
740\begin{array}{lllll}
741{{\bf r}_{}}_{\rm PA} \times {{\bf r}_{}}_{\rm PC}
742& = & [({\lambda_{}}_{\rm A}\; -\; {\lambda_{}}_{\rm P} )
743      ({\phi_{}}_{\rm C}   \; -\; {\phi_{}}_{\rm P} )
744    - ({\lambda_{}}_{\rm C}\; -\; {\lambda_{}}_{\rm P} )
745      ({\phi_{}}_{\rm A}   \; -\; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
746{{\bf r}_{}}_{\rm PB} \times {{\bf r}_{}}_{\rm PA}
747& = & [({\lambda_{}}_{\rm B}\; -\; {\lambda_{}}_{\rm P} )
748      ({\phi_{}}_{\rm A}   \; -\; {\phi_{}}_{\rm P} )
749    - ({\lambda_{}}_{\rm A}\; -\; {\lambda_{}}_{\rm P} )
750      ({\phi_{}}_{\rm B}   \; -\; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
751{{\bf r}_{}}_{\rm PC} \times {{\bf r}_{}}_{\rm PD}
752& = & [({\lambda_{}}_{\rm C}\; -\; {\lambda_{}}_{\rm P} )
753      ({\phi_{}}_{\rm D}   \; -\; {\phi_{}}_{\rm P} )
754    - ({\lambda_{}}_{\rm D}\; -\; {\lambda_{}}_{\rm P} )
755      ({\phi_{}}_{\rm C}   \; -\; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
756{{\bf r}_{}}_{\rm PD} \times {{\bf r}_{}}_{\rm PB}
757& = & [({\lambda_{}}_{\rm D}\; -\; {\lambda_{}}_{\rm P} )
758      ({\phi_{}}_{\rm B}   \; -\; {\phi_{}}_{\rm P} )
759    - ({\lambda_{}}_{\rm B}\; -\; {\lambda_{}}_{\rm P} )
760      ({\phi_{}}_{\rm D}  \;  - \; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
761\end{array}
762\label{eq:cross}
763\end{eqnarray}
764point in the opposite direction to the unit normal
765$\widehat{\bf k}$ (i.e., that the coefficients of
766$\widehat{\bf k}$ are negative),
767where ${{\bf r}_{}}_{\rm PA}$, ${{\bf r}_{}}_{\rm PB}$,
768etc. correspond to the vectors between points P and A,
769P and B, etc.. The method used is
770similar to the method used in
771the SCRIP interpolation package \citep{Jones_1998}.
772
773In order to speed up the grid search, there is the possibility to construct
774a lookup table for a user specified resolution. This lookup
775table contains the lower and upper bounds on the $i$ and $j$ indices
776to be searched for on a regular grid. For each observation position,
777the closest point on the regular grid of this position is computed and
778the $i$ and $j$ ranges of this point searched to determine the precise
779four points surrounding the observation.
780
781\subsection{Parallel aspects of horizontal interpolation}
782\label{OBS_parallel}
783
784For horizontal interpolation, there is the basic problem that the
785observations are unevenly distributed on the globe. In numerical
786models, it is common to divide the model grid into subgrids (or
787domains) where each subgrid is executed on a single processing element
788with explicit message passing for exchange of information along the
789domain boundaries when running on a massively parallel processor (MPP)
790system. This approach is used by \NEMO.
791
792For observations there is no natural distribution since the
793observations are not equally distributed on the globe.
794Two options have been made available: 1) geographical distribution;
795and 2) round-robin.
796
797\subsubsection{Geographical distribution of observations among processors}
798
799%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
800\begin{figure}      \begin{center}
801\includegraphics[width=10cm,height=12cm,angle=-90.]{Fig_ASM_obsdist_local}
802\caption{      \label{fig:obslocal}
803Example of the distribution of observations with the geographical distribution of observational data.} 
804\end{center}      \end{figure}
805%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
806
807This is the simplest option in which the observations are distributed according
808to the domain of the grid-point parallelization. Figure~\ref{fig:obslocal}
809shows an example of the distribution of the {\em in situ} data on processors
810with a different colour for each observation
811on a given processor for a 4 $\times$ 2 decomposition with ORCA2.
812The grid-point domain decomposition is clearly visible on the plot.
813
814The advantage of this approach is that all
815information needed for horizontal interpolation is available without
816any MPP communication. Of course, this is under the assumption that
817we are only using a $2 \times 2$ grid-point stencil for the interpolation
818(e.g., bilinear interpolation). For higher order interpolation schemes this
819is no longer valid. A disadvantage with the above scheme is that the number of
820observations on each processor can be very different. If the cost of
821the actual interpolation is expensive relative to the communication of
822data needed for interpolation, this could lead to load imbalance.
823
824\subsubsection{Round-robin distribution of observations among processors}
825
826%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
827\begin{figure}     \begin{center}
828\includegraphics[width=10cm,height=12cm,angle=-90.]{Fig_ASM_obsdist_global}
829\caption{      \label{fig:obsglobal}
830Example of the distribution of observations with the round-robin distribution of observational data.}
831\end{center}     \end{figure}
832%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
833
834An alternative approach is to distribute the observations equally
835among processors and use message passing in order to retrieve
836the stencil for interpolation. The simplest distribution of the observations
837is to distribute them using a round-robin scheme. Figure~\ref{fig:obsglobal}
838shows the distribution of the {\em in situ} data on processors for the
839round-robin distribution of observations with a different colour for
840each observation on a given processor for a 4 $\times$ 2 decomposition
841with ORCA2 for the same input data as in Fig.~\ref{fig:obslocal}.
842The observations are now clearly randomly distributed on the globe.
843In order to be able to perform horizontal interpolation in this case,
844a subroutine has been developed that retrieves any grid points in the
845global space.
846
847\subsection{Vertical interpolation operator}
848
849Vertical interpolation is achieved using either a cubic spline or
850linear interpolation. For the cubic spline, the top and
851bottom boundary conditions for the second derivative of the
852interpolating polynomial in the spline are set to zero.
853At the bottom boundary, this is done using the land-ocean mask.
854
855\newpage
856
857% ================================================================
858% Offline observation operator documentation
859% ================================================================
860
861%\usepackage{framed}
862
863\section{Offline observation operator}
864\label{OBS_ooo}
865
866\subsection{Concept}
867
868The obs oper maps model variables to observation space. It is possible to apply this mapping
869without running the model. The software which performs this functionality is known as the
870\textbf{offline obs oper}. The obs oper is divided into three stages. An initialisation phase,
871an interpolation phase and an output phase. The implementation of which is outlined in the
872previous sections. During the interpolation phase the offline obs oper populates the model
873arrays by reading saved model fields from disk.
874
875There are two ways of exploiting this offline capacity. The first is to mimic the behaviour of
876the online system by supplying model fields at regular intervals between the start and the end
877of the run. This approach results in a single model counterpart per observation. This kind of
878usage produces feedback files the same file format as the online obs oper.
879The second is to take advantage of the offline setting in which multiple model counterparts can
880be calculated per observation. In this case it is possible to consider all forecasts verifying
881at the same time. By forecast, I mean any method which produces an estimate of physical reality
882which is not an observed value. In the case of class 4 files this means forecasts, analyses, persisted
883analyses and climatological values verifying at the same time. Although the class 4 file format
884doesn't account for multiple ensemble members or multiple experiments per observation, it is possible
885to include these components in the same or multiple files.
886
887%--------------------------------------------------------------------------------------------------------
888% offline_oper.exe
889%--------------------------------------------------------------------------------------------------------
890
891\subsection{Using the offline observation operator}
892
893\subsubsection{Building}
894
895In addition to \emph{OPA\_SRC} the offline obs oper requires the inclusion
896of the \emph{OOO\_SRC} directory. \emph{OOO\_SRC} contains a replacement \textbf{nemo.f90} and
897\textbf{nemogcm.F90} which overwrites the resultant \textbf{nemo.exe}. This is the approach taken
898by \emph{SAS\_SRC} and \emph{OFF\_SRC}.
899
900%--------------------------------------------------------------------------------------------------------
901% Running
902%--------------------------------------------------------------------------------------------------------
903\subsubsection{Running}
904
905The simplest way to use the executable is to edit and append the \textbf{ooo.nml} namelist to
906a full NEMO namelist and then to run the executable as if it were nemo.exe.
907
908\subsubsection{Quick script}
909
910A useful Python utility to control the namelist options can be found in \textbf{OBSTOOLS/OOO}. The
911functions which locate model fields and observation files can be manually specified. The package
912can be installed by appropriate use of the included setup.py script.
913
914Documentation can be auto-generated by Sphinx by running \emph{make html} in the \textbf{doc} directory.
915
916%--------------------------------------------------------------------------------------------------------
917% Configuration section
918%--------------------------------------------------------------------------------------------------------
919\subsection{Configuring the offline observation operator}
920The observation files and settings understood by \textbf{namobs} have been outlined in the online
921obs oper section. In addition there are two further namelists wich control the operation of the offline
922obs oper. \textbf{namooo} which controls the input model fields and \textbf{namcl4} which controls the
923production of class 4 files.
924
925\subsubsection{Single field}
926
927In offline mode model arrays are populated at appropriate time steps via input files.
928At present, \textbf{tsn} and \textbf{sshn} are populated by the default read routines.
929These routines will be expanded upon in future versions to allow the specification of any
930model variable. As such, input files must be global versions of the model domain with
931\textbf{votemper}, \textbf{vosaline} and optionally \textbf{sshn} present.
932
933For each field read there must be an entry in the \textbf{namooo} namelist specifying the
934name of the file to read and the index along the \emph{time\_counter}. For example, to
935read the second time counter from a single file the namelist would be.
936
937\begin{alltt}
938\tiny
939\begin{verbatim} 
940!----------------------------------------------------------------------
941!       namooo Offline obs_oper namelist
942!----------------------------------------------------------------------
943!   ooo_files    specifies the files containing the model counterpart
944!   nn_ooo_idx   specifies the time_counter index within the model file
945&namooo
946   ooo_files = "foo.nc"
947   nn_ooo_idx = 2
948/
949\end{verbatim} 
950\end{alltt}
951
952\subsubsection{Multiple fields per run}
953
954Model field iteration is controlled via \textbf{nn\_ooo\_freq} which specifies
955the number of model steps at which the next field gets read. For example, if
95612 hourly fields are to be interpolated in a setup where 288 steps equals 24 hours.
957
958\begin{alltt}
959\tiny
960\begin{verbatim} 
961!----------------------------------------------------------------------
962!       namooo Offline obs_oper namelist
963!----------------------------------------------------------------------
964!   ooo_files    specifies the files containing the model counterpart
965!   nn_ooo_idx   specifies the time_counter index within the model file
966!   nn_ooo_freq  specifies number of time steps between read operations
967&namooo
968   ooo_files = "foo.nc" "foo.nc"
969   nn_ooo_idx = 1 2
970   nn_ooo_freq = 144
971/
972\end{verbatim} 
973\end{alltt}
974
975The above namelist will result in feedback files whose first 12 hours contain
976the first field of foo.nc and the second 12 hours contain the second field.
977
978%\begin{framed}
979\textbf{Note} Missing files can be denoted as "nofile".
980%\end{framed}
981
982It is easy to see how a collection of fields taken fron a number of files
983at different indices can be combined at a particular frequency in time to
984generate a pseudo model evolution. As long as all that is needed is a single
985model counterpart at a regular interval then namooo is all that needs to
986be edited. However, a far more interesting approach can be taken in which
987multiple forecasts, analyses, persisted analyses and climatologies are
988considered against the same set of observations. For this a slightly more
989complicated approach is needed. It is referred to as \emph{Class 4} since
990it is the fourth metric defined by the GODAE intercomparison project.
991
992%--------------------------------------------------------------------------------------------------------
993% Class 4 file section
994%--------------------------------------------------------------------------------------------------------
995\subsubsection{Multiple model counterparts per observation a.k.a Class 4}
996
997A generalisation of feedback files to allow multiple model components per observation. For a single
998observation, as well as previous forecasts verifying at the same time there are also analyses, persisted
999analyses and climatologies.
1000
1001
1002The above namelist performs two basic functions. It organises the fields
1003given in \textbf{namooo} into groups so that observations can be matched
1004up multiple times. It also controls the metadata and the output variable
1005of the class 4 file when a write routine is called.
1006
1007%\begin{framed}
1008\textbf{Note: ln\_cl4} must be set to \emph{.TRUE.} in \textbf{namobs} 
1009to use class 4 outputs.
1010%\end{framed}
1011
1012\subsubsection{Class 4 naming convention}
1013
1014The standard class 4 file naming convention is as follows.
1015
1016\noindent
1017\linebreak
1018\textbf{\$\{prefix\}\_\$\{yyyymmdd\}\_\$\{sys\}\_\$\{cfg\}\_\$\{vn\}\_\$\{kind\}\_\$\{nproc\}.nc}
1019
1020\noindent
1021\linebreak
1022Much of the namelist is devoted to specifying this convention. The
1023following namelist settings control the elements of the output
1024file names. Each should be specified as a single string of character data.
1025
1026\begin{description}
1027\item[cl4\_prefix]
1028Prefix for class 4 files e.g. class4
1029\item[cl4\_date]
1030YYYYMMDD validity date
1031\item[cl4\_sys]
1032The name of the class 4 model system e.g. FOAM
1033\item[cl4\_cfg]
1034The name of the class 4 model configuration e.g. orca025
1035\item[cl4\_vn]
1036The name of the class 4 model version e.g. 12.0
1037\end{description}
1038
1039\noindent
1040The kind is specified by the observation type internally to the obs oper. The processor
1041number is specified internally in NEMO.
1042
1043\subsubsection{Class 4 file global attributes}
1044
1045Global attributes necessary to fulfill the class 4 file definition. These
1046are also useful pieces of information when collaborating with external
1047partners.
1048
1049\begin{description}
1050\item[cl4\_contact]
1051Contact email for class 4 files.
1052\item[cl4\_inst]
1053The name of the producers institution.
1054\item[cl4\_cfg]
1055The name of the class 4 model configuration e.g. orca025
1056\item[cl4\_vn]
1057The name of the class 4 model version e.g. 12.0
1058\end{description}
1059
1060\noindent
1061The obs\_type,
1062creation date and validity time are specified internally to the obs oper.
1063
1064\subsubsection{Class 4 model counterpart configuration}
1065
1066As seen previously it is possible to perform a single sweep of the
1067obs oper and specify a collection of model fields equally spaced
1068along that sweep. In the class 4 case the single sweep is replaced
1069with multiple sweeps and a certain ammount of book keeping is
1070needed to ensure each model counterpart makes its way to the
1071correct piece of memory in the output files.
1072
1073\noindent
1074\linebreak
1075In terms of book keeping, the offline obs oper needs to know how many
1076full sweeps need to be performed. This is specified via the
1077\textbf{cl4\_match\_len} variable and is the total number of model
1078counterparts per observation. For example, a 3 forecasts plus 3 persistence
1079fields plus an analysis field would be 7 counterparts per observation.
1080
1081\begin{alltt}
1082\tiny
1083\begin{verbatim}
1084   cl4_match_len = 7
1085\end{verbatim}
1086\end{alltt}
1087
1088Then to correctly allocate a class 4 file the forecast axis must be defined. This
1089is controlled via \textbf{cl4\_fcst\_len}, which in out above example would be 3.
1090
1091\begin{alltt}
1092\tiny
1093\begin{verbatim}
1094   cl4_fcst_len = 3
1095\end{verbatim}
1096\end{alltt}
1097
1098Then for each model field it is necessary to designate what class 4 variable and
1099index along the forecast dimension the model counterpart should be stored in the
1100output file. As well as a value for that lead time in hours, this will be useful
1101when interpreting the data afterwards.
1102
1103\begin{alltt}
1104\tiny
1105\begin{verbatim}
1106   cl4_vars = "forecast" "forecast" "forecast" "persistence" "persistence"
1107              "persistence" "best_estimate"
1108   cl4_fcst_idx = 1 2 3 1 2 3 1
1109   cl4_leadtime = 12 36 60
1110\end{verbatim}
1111\end{alltt}
1112
1113In terms of files and indices of fields inside each file the class 4 approach
1114makes use of the \textbf{namooo} namelist. If our fields are in separate files
1115with a single field per file our example inputs will be specified.
1116
1117\begin{alltt}
1118\tiny
1119\begin{verbatim}
1120   ooo_files = "F.1.nc" "F.2.nc" "F.3.nc" "P.1.nc" "P.2.nc" "P.3.nc" "A.1.nc"
1121   nn_ooo_idx = 1 1 1 1 1 1 1
1122\end{verbatim}
1123\end{alltt}
1124
1125When we combine all of the naming conventions, global attributes and i/o instructions
1126the class 4 namelist becomes.
1127
1128\begin{alltt}
1129\tiny
1130\begin{verbatim}
1131!----------------------------------------------------------------------
1132!       namooo Offline obs_oper namelist
1133!----------------------------------------------------------------------
1134!   ooo_files    specifies the files containing the model counterpart
1135!   nn_ooo_idx   specifies the time_counter index within the model file
1136!   nn_ooo_freq  specifies number of time steps between read operations
1137&namooo
1138   ooo_files = "F.1.nc" "F.2.nc" "F.3.nc" "P.1.nc" "P.2.nc" "P.3.nc" "A.1.nc"
1139   nn_ooo_idx = 1 1 1 1 1 1 1
1140/
1141!----------------------------------------------------------------------
1142!       namcl4 Offline obs_oper class 4 namelist
1143!----------------------------------------------------------------------
1144!
1145!  Naming convention
1146!  -----------------
1147!  cl4_prefix    specifies the output file prefix
1148!  cl4_date      specifies the output file validity date
1149!  cl4_sys       specifies the model counterpart system
1150!  cl4_cfg       specifies the model counterpart configuration
1151!  cl4_vn        specifies the model counterpart version
1152!  cl4_inst      specifies the model counterpart institute
1153!  cl4_contact   specifies the file producers contact details
1154!
1155!  I/O specification
1156!  -----------------
1157!  cl4_vars      specifies the names of the output file netcdf variable
1158!  cl4_fcst_idx  specifies output file forecast index
1159!  cl4_fcst_len  specifies forecast axis length
1160!  cl4_match_len specifies number of unique matches per observation
1161!  cl4_leadtime  specifies the forecast axis lead time
1162!
1163&namcl4
1164   cl4_match_len = 7
1165   cl4_fcst_len = 3
1166   cl4_fcst_idx = 1 2 3 1 2 3 1
1167   cl4_vars = "forecast" "forecast" "forecast" "persistence" "persistence"
1168              "persistence" "best_estimate"
1169   cl4_leadtime = 12 36 60
1170   cl4_prefix = "class4"
1171   cl4_date = "20130101"
1172   cl4_vn = "12.0"
1173   cl4_sys = "FOAM"
1174   cl4_cfg = "AMM7"
1175   cl4_contact = "example@example.com"
1176   cl4_inst = "UK Met Office"
1177/
1178\end{verbatim}
1179\end{alltt}
1180
1181\subsubsection{Climatology interpolation}
1182
1183The climatological counterpart is generated at the start of the run by restarting
1184the model from climatology through appropriate use of \textbf{namtsd}. To override
1185the offline observation operator read routine and to take advantage of the restart
1186settings, specify the first entry in \textbf{cl4\_vars} as "climatology". This will then
1187pipe the restart from climatology into the output class 4 file. As in every other
1188class 4 matchup the input file, input index and output index must be specified.
1189These can be replaced with dummy data since they are not used but they must be
1190present to cycle through the matchups correctly.
1191
1192\subsection{Advanced usage}
1193
1194In certain cases it may be desirable to include both multiple model fields per
1195observation window with multiple match ups per observation. This can be achieved
1196by specifying \textbf{nn\_ooo\_freq} as well as the class 4 settings. Care must
1197be taken in generating the ooo\_files list such that the files are arranged into
1198consecutive blocks of single match ups. For example, 2 forecast fields
1199of 12 hourly data would result in 4 separate read operations but only 2 write
1200operations, 1 per forecast.
1201
1202\begin{alltt}
1203\tiny
1204\begin{verbatim}
1205   ooo_files = "F1.nc" "F1.nc" "F2.nc" "F2.nc"
1206...
1207   cl4_fcst_idx = 1 2
1208\end{verbatim}
1209\end{alltt}
1210
1211The above notation reveals the internal split between match up iterators and file
1212iterators. This technique has not been used before so experimentation is needed
1213before results can be trusted.
1214
1215
1216
1217
1218\newpage
1219
1220\section{Observation Utilities}
1221\label{OBS_obsutils}
1222
1223Some tools for viewing and processing of observation and feedback files are provided in the
1224NEMO repository for convenience. These include OBSTOOLS which are a collection of Fortran
1225programs which are helpful to deal with feedback files. They do such tasks as observation file
1226conversion, printing of file contents, some basic statistical analysis of feedback files. The
1227other tool is an IDL program called dataplot which uses a graphical interface to visualise
1228observations and feedback files. OBSTOOLS and dataplot are described in more detail below. 
1229
1230\subsection{Obstools}
1231
1232A series of Fortran utilities is provided with NEMO called OBSTOOLS. This are helpful in
1233handling observation files and the feedback file output from the NEMO observation operator.
1234The utilities are as follows
1235
1236\subsubsection{c4comb}
1237
1238The program c4comb combines multiple class 4 files produced by individual processors in an
1239MPI run of NEMO offline obs\_oper into a single class 4 file. The program is called in the following way:
1240
1241\begin{alltt}
1242\footnotesize
1243\begin{verbatim}
1244c4comb.exe outputfile inputfile1 inputfile2 ...
1245\end{verbatim}
1246\end{alltt}
1247
1248\subsubsection{corio2fb}
1249
1250The program corio2fb converts profile observation files from the Coriolis format to the
1251standard feedback format. The program is called in the following way:
1252
1253\begin{alltt}
1254\footnotesize
1255\begin{verbatim}
1256corio2fb.exe outputfile inputfile1 inputfile2 ...
1257\end{verbatim}
1258\end{alltt}
1259
1260\subsubsection{enact2fb}
1261
1262The program enact2fb converts profile observation files from the ENACT format to the standard
1263feedback format. The program is called in the following way:
1264
1265\begin{alltt}
1266\footnotesize
1267\begin{verbatim}
1268enact2fb.exe outputfile inputfile1 inputfile2 ...
1269\end{verbatim}
1270\end{alltt}
1271
1272\subsubsection{fbcomb}
1273
1274The program fbcomb combines multiple feedback files produced by individual processors in an
1275MPI run of NEMO into a single feedback file. The program is called in the following way:
1276
1277\begin{alltt}
1278\footnotesize
1279\begin{verbatim}
1280fbcomb.exe outputfile inputfile1 inputfile2 ...
1281\end{verbatim}
1282\end{alltt}
1283
1284\subsubsection{fbmatchup}
1285
1286The program fbmatchup will match observations from two feedback files. The program is called
1287in the following way:
1288
1289\begin{alltt}
1290\footnotesize
1291\begin{verbatim}
1292fbmatchup.exe outputfile inputfile1 varname1 inputfile2 varname2 ...
1293\end{verbatim}
1294\end{alltt}
1295
1296
1297\subsubsection{fbprint}
1298
1299The program fbprint will print the contents of a feedback file or files to standard output.
1300Selected information can be output using optional arguments. The program is called in the
1301following way:
1302
1303\begin{alltt}
1304\footnotesize
1305\begin{verbatim}
1306fbprint.exe [options] inputfile
1307
1308options:
1309     -b            shorter output
1310     -q            Select observations based on QC flags
1311     -Q            Select observations based on QC flags
1312     -B            Select observations based on QC flags
1313     -u            unsorted
1314     -s ID         select station ID 
1315     -t TYPE       select observation type
1316     -v NUM1-NUM2  select variable range to print by number
1317                      (default all)
1318     -a NUM1-NUM2  select additional variable range to print by number
1319                      (default all)
1320     -e NUM1-NUM2  select extra variable range to print by number
1321                      (default all)
1322     -d            output date range
1323     -D            print depths
1324     -z            use zipped files
1325\end{verbatim}
1326\end{alltt}
1327
1328\subsubsection{fbsel}
1329
1330The program fbsel will select or subsample observations. The program is called in the
1331following way:
1332
1333\begin{alltt}
1334\footnotesize
1335\begin{verbatim}
1336fbsel.exe <input filename> <output filename>
1337\end{verbatim}
1338\end{alltt}
1339
1340\subsubsection{fbstat}
1341
1342The program fbstat will output summary statistics in different global areas into a number of
1343files. The program is called in the following way:
1344
1345\begin{alltt}
1346\footnotesize
1347\begin{verbatim}
1348fbstat.exe [-nmlev] <filenames>
1349\end{verbatim}
1350\end{alltt}
1351
1352\subsubsection{fbthin}
1353
1354The program fbthin will thin the data to 1 degree resolution. The code could easily be
1355modified to thin to a different resolution. The program is called in the following way:
1356
1357\begin{alltt}
1358\footnotesize
1359\begin{verbatim}
1360fbthin.exe inputfile outputfile
1361\end{verbatim}
1362\end{alltt}
1363
1364\subsubsection{sla2fb}
1365
1366The program sla2fb will convert an AVISO SLA format file to feedback format. The program is
1367called in the following way:
1368
1369\begin{alltt}
1370\footnotesize
1371\begin{verbatim}
1372sla2fb.exe [-s type] outputfile inputfile1 inputfile2 ...
1373
1374Option:
1375     -s            Select altimeter data_source
1376\end{verbatim}
1377\end{alltt}
1378
1379\subsubsection{vel2fb}
1380
1381The program vel2fb will convert TAO/PIRATA/RAMA currents files to feedback format. The program
1382is called in the following way:
1383
1384\begin{alltt}
1385\footnotesize
1386\begin{verbatim}
1387vel2fb.exe outputfile inputfile1 inputfile2 ...
1388\end{verbatim}
1389\end{alltt}
1390
1391\subsection{building the obstools}
1392
1393To build the obstools use in the tools directory use ./maketools -n OBSTOOLS -m [ARCH].
1394
1395\subsection{Dataplot}
1396
1397An IDL program called dataplot is included which uses a graphical interface to visualise
1398observations and feedback files. It is possible to zoom in, plot individual profiles and
1399calculate some basic statistics. To plot some data run IDL and then:
1400\begin{alltt}
1401\footnotesize
1402\begin{verbatim}
1403IDL> dataplot, "filename"
1404\end{verbatim}
1405\end{alltt}
1406
1407To read multiple files into dataplot, for example multiple feedback files from different
1408processors or from different days, the easiest method is to use the spawn command to generate
1409a list of files which can then be passed to dataplot.
1410\begin{alltt}
1411\footnotesize
1412\begin{verbatim}
1413IDL> spawn, 'ls profb*.nc', files
1414IDL> dataplot, files
1415\end{verbatim}
1416\end{alltt}
1417
1418Fig~\ref{fig:obsdataplotmain} shows the main window which is launched when dataplot starts.
1419This is split into three parts. At the top there is a menu bar which contains a variety of
1420drop down menus. Areas - zooms into prespecified regions; plot - plots the data as a
1421timeseries or a T-S diagram if appropriate; Find - allows data to be searched; Config - sets
1422various configuration options.
1423
1424The middle part is a plot of the geographical location of the observations. This will plot the
1425observation value, the model background value or observation minus background value depending
1426on the option selected in the radio button at the bottom of the window. The plotting colour
1427range can be changed by clicking on the colour bar. The title of the plot gives some basic
1428information about the date range and depth range shown, the extreme values, and the mean and
1429rms values. It is possible to zoom in using a drag-box. You may also zoom in or out using the
1430mouse wheel.
1431
1432The bottom part of the window controls what is visible in the plot above. There are two bars
1433which select the level range plotted (for profile data). The other bars below select the date
1434range shown. The bottom of the figure allows the option to plot the mean, root mean square,
1435standard deviation or mean square values. As mentioned above you can choose to plot the
1436observation value, the model background value or observation minus background value. The next
1437group of radio buttons selects the map projection. This can either be regular latitude
1438longitude grid, or north or south polar stereographic. The next group of radio buttons will
1439plot bad observations, switch to salinity and plot density for profile observations. The
1440rightmost group of buttons will print the plot window as a postscript, save it as png, or exit
1441from dataplot.
1442
1443%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1444\begin{figure}     \begin{center}
1445%\includegraphics[width=10cm,height=12cm,angle=-90.]{Fig_OBS_dataplot_main}
1446\includegraphics[width=9cm,angle=-90.]{Fig_OBS_dataplot_main}
1447\caption{      \label{fig:obsdataplotmain}
1448Main window of dataplot.}
1449\end{center}     \end{figure}
1450%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1451
1452If a profile point is clicked with the mouse button a plot of the observation and background
1453values as a function of depth (Fig~\ref{fig:obsdataplotprofile}).
1454
1455%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1456\begin{figure}     \begin{center}
1457%\includegraphics[width=10cm,height=12cm,angle=-90.]{Fig_OBS_dataplot_prof}
1458\includegraphics[width=7cm,angle=-90.]{Fig_OBS_dataplot_prof}
1459\caption{      \label{fig:obsdataplotprofile}
1460Profile plot from dataplot produced by right clicking on a point in the main window.}
1461\end{center}     \end{figure}
1462%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1463
1464
1465
1466
1467\end{document}
Note: See TracBrowser for help on using the repository browser.