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 @ 8919

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

Updated the NEMO book documentation with information on the observation averaging operator.
File size: 62.7 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 SSS data)
566have a footprint which is similar 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 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 n2dint=0: Distance-weighted interpolation
574\item n2dint=1: Distance-weighted interpolation (small angle)
575\item n2dint=2: Bilinear interpolation (geographical grid)
576\item n2dint=3: Bilinear remapping interpolation (general grid)
577\item n2dint=4: Polynomial interpolation
578\item n2dint=5: Radial footprint averaging with radius specified in the namelist as rn_???_avglamscl in degrees or metres (set using ln_???_fp_indegs)
579\item n2dint=6: Rectangular footprint averaging with E/W and N/S size specified in the namelist as rn_???_avglamscl and rn_???_avgphiscl in degrees or metres (set using 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 n2dint default option can be overridden for surface observation types using namelist values 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
697%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
698\begin{figure}      \begin{center}
699\includegraphics[width=10cm,height=12cm]{Fig_OBS_avg_rec}
700\caption{      \label{fig:obsavgrec}
701Weights associated with each model grid box (blue lines and numbers) for an observation at -170.5E, 56.0N with a footprint of $1\deg \time 1\deg$} 
702\end{center}      \end{figure}
703%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
704
705
706%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
707\begin{figure}      \begin{center}
708\includegraphics[width=10cm,height=12cm]{Fig_OBS_avg_rad}
709\caption{      \label{fig:obsavgrad}
710As for figure \ref{obsavgrec} but for a radial footprint with diameter $1\deg$.
711\end{center}      \end{figure}
712%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
713
714
715\subsection{Grid search}
716
717For many grids used by the NEMO model, such as the ORCA family,
718the horizontal grid coordinates $i$ and $j$ are not simple functions
719of latitude and longitude. Therefore, it is not always straightforward
720to determine the grid points surrounding any given observational position.
721Before the interpolation can be performed, a search
722algorithm is then required to determine the corner points of
723the quadrilateral cell in which the observation is located.
724This is the most difficult and time consuming part of the
7252D interpolation procedure.
726A robust test for determining if an observation falls
727within a given quadrilateral cell is as follows. Let
728${\rm P}({\lambda_{}}_{\rm P} ,{\phi_{}}_{\rm P} )$ denote the observation point,
729and let ${\rm A}({\lambda_{}}_{\rm A} ,{\phi_{}}_{\rm A} )$,
730${\rm B}({\lambda_{}}_{\rm B} ,{\phi_{}}_{\rm B} )$,
731${\rm C}({\lambda_{}}_{\rm C} ,{\phi_{}}_{\rm C} )$ 
732and
733${\rm D}({\lambda_{}}_{\rm D} ,{\phi_{}}_{\rm D} )$ denote
734the bottom left, bottom right, top left and top right
735corner points of the cell, respectively.
736To determine if P is inside
737the cell, we verify that the cross-products
738\begin{eqnarray}
739\begin{array}{lllll}
740{{\bf r}_{}}_{\rm PA} \times {{\bf r}_{}}_{\rm PC}
741& = & [({\lambda_{}}_{\rm A}\; -\; {\lambda_{}}_{\rm P} )
742      ({\phi_{}}_{\rm C}   \; -\; {\phi_{}}_{\rm P} )
743    - ({\lambda_{}}_{\rm C}\; -\; {\lambda_{}}_{\rm P} )
744      ({\phi_{}}_{\rm A}   \; -\; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
745{{\bf r}_{}}_{\rm PB} \times {{\bf r}_{}}_{\rm PA}
746& = & [({\lambda_{}}_{\rm B}\; -\; {\lambda_{}}_{\rm P} )
747      ({\phi_{}}_{\rm A}   \; -\; {\phi_{}}_{\rm P} )
748    - ({\lambda_{}}_{\rm A}\; -\; {\lambda_{}}_{\rm P} )
749      ({\phi_{}}_{\rm B}   \; -\; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
750{{\bf r}_{}}_{\rm PC} \times {{\bf r}_{}}_{\rm PD}
751& = & [({\lambda_{}}_{\rm C}\; -\; {\lambda_{}}_{\rm P} )
752      ({\phi_{}}_{\rm D}   \; -\; {\phi_{}}_{\rm P} )
753    - ({\lambda_{}}_{\rm D}\; -\; {\lambda_{}}_{\rm P} )
754      ({\phi_{}}_{\rm C}   \; -\; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
755{{\bf r}_{}}_{\rm PD} \times {{\bf r}_{}}_{\rm PB}
756& = & [({\lambda_{}}_{\rm D}\; -\; {\lambda_{}}_{\rm P} )
757      ({\phi_{}}_{\rm B}   \; -\; {\phi_{}}_{\rm P} )
758    - ({\lambda_{}}_{\rm B}\; -\; {\lambda_{}}_{\rm P} )
759      ({\phi_{}}_{\rm D}  \;  - \; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\
760\end{array}
761\label{eq:cross}
762\end{eqnarray}
763point in the opposite direction to the unit normal
764$\widehat{\bf k}$ (i.e., that the coefficients of
765$\widehat{\bf k}$ are negative),
766where ${{\bf r}_{}}_{\rm PA}$, ${{\bf r}_{}}_{\rm PB}$,
767etc. correspond to the vectors between points P and A,
768P and B, etc.. The method used is
769similar to the method used in
770the SCRIP interpolation package \citep{Jones_1998}.
771
772In order to speed up the grid search, there is the possibility to construct
773a lookup table for a user specified resolution. This lookup
774table contains the lower and upper bounds on the $i$ and $j$ indices
775to be searched for on a regular grid. For each observation position,
776the closest point on the regular grid of this position is computed and
777the $i$ and $j$ ranges of this point searched to determine the precise
778four points surrounding the observation.
779
780\subsection{Parallel aspects of horizontal interpolation}
781\label{OBS_parallel}
782
783For horizontal interpolation, there is the basic problem that the
784observations are unevenly distributed on the globe. In numerical
785models, it is common to divide the model grid into subgrids (or
786domains) where each subgrid is executed on a single processing element
787with explicit message passing for exchange of information along the
788domain boundaries when running on a massively parallel processor (MPP)
789system. This approach is used by \NEMO.
790
791For observations there is no natural distribution since the
792observations are not equally distributed on the globe.
793Two options have been made available: 1) geographical distribution;
794and 2) round-robin.
795
796\subsubsection{Geographical distribution of observations among processors}
797
798%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
799\begin{figure}      \begin{center}
800\includegraphics[width=10cm,height=12cm,angle=-90.]{Fig_ASM_obsdist_local}
801\caption{      \label{fig:obslocal}
802Example of the distribution of observations with the geographical distribution of observational data.} 
803\end{center}      \end{figure}
804%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
805
806This is the simplest option in which the observations are distributed according
807to the domain of the grid-point parallelization. Figure~\ref{fig:obslocal}
808shows an example of the distribution of the {\em in situ} data on processors
809with a different colour for each observation
810on a given processor for a 4 $\times$ 2 decomposition with ORCA2.
811The grid-point domain decomposition is clearly visible on the plot.
812
813The advantage of this approach is that all
814information needed for horizontal interpolation is available without
815any MPP communication. Of course, this is under the assumption that
816we are only using a $2 \times 2$ grid-point stencil for the interpolation
817(e.g., bilinear interpolation). For higher order interpolation schemes this
818is no longer valid. A disadvantage with the above scheme is that the number of
819observations on each processor can be very different. If the cost of
820the actual interpolation is expensive relative to the communication of
821data needed for interpolation, this could lead to load imbalance.
822
823\subsubsection{Round-robin distribution of observations among processors}
824
825%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
826\begin{figure}     \begin{center}
827\includegraphics[width=10cm,height=12cm,angle=-90.]{Fig_ASM_obsdist_global}
828\caption{      \label{fig:obsglobal}
829Example of the distribution of observations with the round-robin distribution of observational data.}
830\end{center}     \end{figure}
831%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
832
833An alternative approach is to distribute the observations equally
834among processors and use message passing in order to retrieve
835the stencil for interpolation. The simplest distribution of the observations
836is to distribute them using a round-robin scheme. Figure~\ref{fig:obsglobal}
837shows the distribution of the {\em in situ} data on processors for the
838round-robin distribution of observations with a different colour for
839each observation on a given processor for a 4 $\times$ 2 decomposition
840with ORCA2 for the same input data as in Fig.~\ref{fig:obslocal}.
841The observations are now clearly randomly distributed on the globe.
842In order to be able to perform horizontal interpolation in this case,
843a subroutine has been developed that retrieves any grid points in the
844global space.
845
846\subsection{Vertical interpolation operator}
847
848Vertical interpolation is achieved using either a cubic spline or
849linear interpolation. For the cubic spline, the top and
850bottom boundary conditions for the second derivative of the
851interpolating polynomial in the spline are set to zero.
852At the bottom boundary, this is done using the land-ocean mask.
853
854\newpage
855
856% ================================================================
857% Offline observation operator documentation
858% ================================================================
859
860%\usepackage{framed}
861
862\section{Offline observation operator}
863\label{OBS_ooo}
864
865\subsection{Concept}
866
867The obs oper maps model variables to observation space. It is possible to apply this mapping
868without running the model. The software which performs this functionality is known as the
869\textbf{offline obs oper}. The obs oper is divided into three stages. An initialisation phase,
870an interpolation phase and an output phase. The implementation of which is outlined in the
871previous sections. During the interpolation phase the offline obs oper populates the model
872arrays by reading saved model fields from disk.
873
874There are two ways of exploiting this offline capacity. The first is to mimic the behaviour of
875the online system by supplying model fields at regular intervals between the start and the end
876of the run. This approach results in a single model counterpart per observation. This kind of
877usage produces feedback files the same file format as the online obs oper.
878The second is to take advantage of the offline setting in which multiple model counterparts can
879be calculated per observation. In this case it is possible to consider all forecasts verifying
880at the same time. By forecast, I mean any method which produces an estimate of physical reality
881which is not an observed value. In the case of class 4 files this means forecasts, analyses, persisted
882analyses and climatological values verifying at the same time. Although the class 4 file format
883doesn't account for multiple ensemble members or multiple experiments per observation, it is possible
884to include these components in the same or multiple files.
885
886%--------------------------------------------------------------------------------------------------------
887% offline_oper.exe
888%--------------------------------------------------------------------------------------------------------
889
890\subsection{Using the offline observation operator}
891
892\subsubsection{Building}
893
894In addition to \emph{OPA\_SRC} the offline obs oper requires the inclusion
895of the \emph{OOO\_SRC} directory. \emph{OOO\_SRC} contains a replacement \textbf{nemo.f90} and
896\textbf{nemogcm.F90} which overwrites the resultant \textbf{nemo.exe}. This is the approach taken
897by \emph{SAS\_SRC} and \emph{OFF\_SRC}.
898
899%--------------------------------------------------------------------------------------------------------
900% Running
901%--------------------------------------------------------------------------------------------------------
902\subsubsection{Running}
903
904The simplest way to use the executable is to edit and append the \textbf{ooo.nml} namelist to
905a full NEMO namelist and then to run the executable as if it were nemo.exe.
906
907\subsubsection{Quick script}
908
909A useful Python utility to control the namelist options can be found in \textbf{OBSTOOLS/OOO}. The
910functions which locate model fields and observation files can be manually specified. The package
911can be installed by appropriate use of the included setup.py script.
912
913Documentation can be auto-generated by Sphinx by running \emph{make html} in the \textbf{doc} directory.
914
915%--------------------------------------------------------------------------------------------------------
916% Configuration section
917%--------------------------------------------------------------------------------------------------------
918\subsection{Configuring the offline observation operator}
919The observation files and settings understood by \textbf{namobs} have been outlined in the online
920obs oper section. In addition there are two further namelists wich control the operation of the offline
921obs oper. \textbf{namooo} which controls the input model fields and \textbf{namcl4} which controls the
922production of class 4 files.
923
924\subsubsection{Single field}
925
926In offline mode model arrays are populated at appropriate time steps via input files.
927At present, \textbf{tsn} and \textbf{sshn} are populated by the default read routines.
928These routines will be expanded upon in future versions to allow the specification of any
929model variable. As such, input files must be global versions of the model domain with
930\textbf{votemper}, \textbf{vosaline} and optionally \textbf{sshn} present.
931
932For each field read there must be an entry in the \textbf{namooo} namelist specifying the
933name of the file to read and the index along the \emph{time\_counter}. For example, to
934read the second time counter from a single file the namelist would be.
935
936\begin{alltt}
937\tiny
938\begin{verbatim} 
939!----------------------------------------------------------------------
940!       namooo Offline obs_oper namelist
941!----------------------------------------------------------------------
942!   ooo_files    specifies the files containing the model counterpart
943!   nn_ooo_idx   specifies the time_counter index within the model file
944&namooo
945   ooo_files = "foo.nc"
946   nn_ooo_idx = 2
947/
948\end{verbatim} 
949\end{alltt}
950
951\subsubsection{Multiple fields per run}
952
953Model field iteration is controlled via \textbf{nn\_ooo\_freq} which specifies
954the number of model steps at which the next field gets read. For example, if
95512 hourly fields are to be interpolated in a setup where 288 steps equals 24 hours.
956
957\begin{alltt}
958\tiny
959\begin{verbatim} 
960!----------------------------------------------------------------------
961!       namooo Offline obs_oper namelist
962!----------------------------------------------------------------------
963!   ooo_files    specifies the files containing the model counterpart
964!   nn_ooo_idx   specifies the time_counter index within the model file
965!   nn_ooo_freq  specifies number of time steps between read operations
966&namooo
967   ooo_files = "foo.nc" "foo.nc"
968   nn_ooo_idx = 1 2
969   nn_ooo_freq = 144
970/
971\end{verbatim} 
972\end{alltt}
973
974The above namelist will result in feedback files whose first 12 hours contain
975the first field of foo.nc and the second 12 hours contain the second field.
976
977%\begin{framed}
978\textbf{Note} Missing files can be denoted as "nofile".
979%\end{framed}
980
981It is easy to see how a collection of fields taken fron a number of files
982at different indices can be combined at a particular frequency in time to
983generate a pseudo model evolution. As long as all that is needed is a single
984model counterpart at a regular interval then namooo is all that needs to
985be edited. However, a far more interesting approach can be taken in which
986multiple forecasts, analyses, persisted analyses and climatologies are
987considered against the same set of observations. For this a slightly more
988complicated approach is needed. It is referred to as \emph{Class 4} since
989it is the fourth metric defined by the GODAE intercomparison project.
990
991%--------------------------------------------------------------------------------------------------------
992% Class 4 file section
993%--------------------------------------------------------------------------------------------------------
994\subsubsection{Multiple model counterparts per observation a.k.a Class 4}
995
996A generalisation of feedback files to allow multiple model components per observation. For a single
997observation, as well as previous forecasts verifying at the same time there are also analyses, persisted
998analyses and climatologies.
999
1000
1001The above namelist performs two basic functions. It organises the fields
1002given in \textbf{namooo} into groups so that observations can be matched
1003up multiple times. It also controls the metadata and the output variable
1004of the class 4 file when a write routine is called.
1005
1006%\begin{framed}
1007\textbf{Note: ln\_cl4} must be set to \emph{.TRUE.} in \textbf{namobs} 
1008to use class 4 outputs.
1009%\end{framed}
1010
1011\subsubsection{Class 4 naming convention}
1012
1013The standard class 4 file naming convention is as follows.
1014
1015\noindent
1016\linebreak
1017\textbf{\$\{prefix\}\_\$\{yyyymmdd\}\_\$\{sys\}\_\$\{cfg\}\_\$\{vn\}\_\$\{kind\}\_\$\{nproc\}.nc}
1018
1019\noindent
1020\linebreak
1021Much of the namelist is devoted to specifying this convention. The
1022following namelist settings control the elements of the output
1023file names. Each should be specified as a single string of character data.
1024
1025\begin{description}
1026\item[cl4\_prefix]
1027Prefix for class 4 files e.g. class4
1028\item[cl4\_date]
1029YYYYMMDD validity date
1030\item[cl4\_sys]
1031The name of the class 4 model system e.g. FOAM
1032\item[cl4\_cfg]
1033The name of the class 4 model configuration e.g. orca025
1034\item[cl4\_vn]
1035The name of the class 4 model version e.g. 12.0
1036\end{description}
1037
1038\noindent
1039The kind is specified by the observation type internally to the obs oper. The processor
1040number is specified internally in NEMO.
1041
1042\subsubsection{Class 4 file global attributes}
1043
1044Global attributes necessary to fulfill the class 4 file definition. These
1045are also useful pieces of information when collaborating with external
1046partners.
1047
1048\begin{description}
1049\item[cl4\_contact]
1050Contact email for class 4 files.
1051\item[cl4\_inst]
1052The name of the producers institution.
1053\item[cl4\_cfg]
1054The name of the class 4 model configuration e.g. orca025
1055\item[cl4\_vn]
1056The name of the class 4 model version e.g. 12.0
1057\end{description}
1058
1059\noindent
1060The obs\_type,
1061creation date and validity time are specified internally to the obs oper.
1062
1063\subsubsection{Class 4 model counterpart configuration}
1064
1065As seen previously it is possible to perform a single sweep of the
1066obs oper and specify a collection of model fields equally spaced
1067along that sweep. In the class 4 case the single sweep is replaced
1068with multiple sweeps and a certain ammount of book keeping is
1069needed to ensure each model counterpart makes its way to the
1070correct piece of memory in the output files.
1071
1072\noindent
1073\linebreak
1074In terms of book keeping, the offline obs oper needs to know how many
1075full sweeps need to be performed. This is specified via the
1076\textbf{cl4\_match\_len} variable and is the total number of model
1077counterparts per observation. For example, a 3 forecasts plus 3 persistence
1078fields plus an analysis field would be 7 counterparts per observation.
1079
1080\begin{alltt}
1081\tiny
1082\begin{verbatim}
1083   cl4_match_len = 7
1084\end{verbatim}
1085\end{alltt}
1086
1087Then to correctly allocate a class 4 file the forecast axis must be defined. This
1088is controlled via \textbf{cl4\_fcst\_len}, which in out above example would be 3.
1089
1090\begin{alltt}
1091\tiny
1092\begin{verbatim}
1093   cl4_fcst_len = 3
1094\end{verbatim}
1095\end{alltt}
1096
1097Then for each model field it is necessary to designate what class 4 variable and
1098index along the forecast dimension the model counterpart should be stored in the
1099output file. As well as a value for that lead time in hours, this will be useful
1100when interpreting the data afterwards.
1101
1102\begin{alltt}
1103\tiny
1104\begin{verbatim}
1105   cl4_vars = "forecast" "forecast" "forecast" "persistence" "persistence"
1106              "persistence" "best_estimate"
1107   cl4_fcst_idx = 1 2 3 1 2 3 1
1108   cl4_leadtime = 12 36 60
1109\end{verbatim}
1110\end{alltt}
1111
1112In terms of files and indices of fields inside each file the class 4 approach
1113makes use of the \textbf{namooo} namelist. If our fields are in separate files
1114with a single field per file our example inputs will be specified.
1115
1116\begin{alltt}
1117\tiny
1118\begin{verbatim}
1119   ooo_files = "F.1.nc" "F.2.nc" "F.3.nc" "P.1.nc" "P.2.nc" "P.3.nc" "A.1.nc"
1120   nn_ooo_idx = 1 1 1 1 1 1 1
1121\end{verbatim}
1122\end{alltt}
1123
1124When we combine all of the naming conventions, global attributes and i/o instructions
1125the class 4 namelist becomes.
1126
1127\begin{alltt}
1128\tiny
1129\begin{verbatim}
1130!----------------------------------------------------------------------
1131!       namooo Offline obs_oper namelist
1132!----------------------------------------------------------------------
1133!   ooo_files    specifies the files containing the model counterpart
1134!   nn_ooo_idx   specifies the time_counter index within the model file
1135!   nn_ooo_freq  specifies number of time steps between read operations
1136&namooo
1137   ooo_files = "F.1.nc" "F.2.nc" "F.3.nc" "P.1.nc" "P.2.nc" "P.3.nc" "A.1.nc"
1138   nn_ooo_idx = 1 1 1 1 1 1 1
1139/
1140!----------------------------------------------------------------------
1141!       namcl4 Offline obs_oper class 4 namelist
1142!----------------------------------------------------------------------
1143!
1144!  Naming convention
1145!  -----------------
1146!  cl4_prefix    specifies the output file prefix
1147!  cl4_date      specifies the output file validity date
1148!  cl4_sys       specifies the model counterpart system
1149!  cl4_cfg       specifies the model counterpart configuration
1150!  cl4_vn        specifies the model counterpart version
1151!  cl4_inst      specifies the model counterpart institute
1152!  cl4_contact   specifies the file producers contact details
1153!
1154!  I/O specification
1155!  -----------------
1156!  cl4_vars      specifies the names of the output file netcdf variable
1157!  cl4_fcst_idx  specifies output file forecast index
1158!  cl4_fcst_len  specifies forecast axis length
1159!  cl4_match_len specifies number of unique matches per observation
1160!  cl4_leadtime  specifies the forecast axis lead time
1161!
1162&namcl4
1163   cl4_match_len = 7
1164   cl4_fcst_len = 3
1165   cl4_fcst_idx = 1 2 3 1 2 3 1
1166   cl4_vars = "forecast" "forecast" "forecast" "persistence" "persistence"
1167              "persistence" "best_estimate"
1168   cl4_leadtime = 12 36 60
1169   cl4_prefix = "class4"
1170   cl4_date = "20130101"
1171   cl4_vn = "12.0"
1172   cl4_sys = "FOAM"
1173   cl4_cfg = "AMM7"
1174   cl4_contact = "example@example.com"
1175   cl4_inst = "UK Met Office"
1176/
1177\end{verbatim}
1178\end{alltt}
1179
1180\subsubsection{Climatology interpolation}
1181
1182The climatological counterpart is generated at the start of the run by restarting
1183the model from climatology through appropriate use of \textbf{namtsd}. To override
1184the offline observation operator read routine and to take advantage of the restart
1185settings, specify the first entry in \textbf{cl4\_vars} as "climatology". This will then
1186pipe the restart from climatology into the output class 4 file. As in every other
1187class 4 matchup the input file, input index and output index must be specified.
1188These can be replaced with dummy data since they are not used but they must be
1189present to cycle through the matchups correctly.
1190
1191\subsection{Advanced usage}
1192
1193In certain cases it may be desirable to include both multiple model fields per
1194observation window with multiple match ups per observation. This can be achieved
1195by specifying \textbf{nn\_ooo\_freq} as well as the class 4 settings. Care must
1196be taken in generating the ooo\_files list such that the files are arranged into
1197consecutive blocks of single match ups. For example, 2 forecast fields
1198of 12 hourly data would result in 4 separate read operations but only 2 write
1199operations, 1 per forecast.
1200
1201\begin{alltt}
1202\tiny
1203\begin{verbatim}
1204   ooo_files = "F1.nc" "F1.nc" "F2.nc" "F2.nc"
1205...
1206   cl4_fcst_idx = 1 2
1207\end{verbatim}
1208\end{alltt}
1209
1210The above notation reveals the internal split between match up iterators and file
1211iterators. This technique has not been used before so experimentation is needed
1212before results can be trusted.
1213
1214
1215
1216
1217\newpage
1218
1219\section{Observation Utilities}
1220\label{OBS_obsutils}
1221
1222Some tools for viewing and processing of observation and feedback files are provided in the
1223NEMO repository for convenience. These include OBSTOOLS which are a collection of Fortran
1224programs which are helpful to deal with feedback files. They do such tasks as observation file
1225conversion, printing of file contents, some basic statistical analysis of feedback files. The
1226other tool is an IDL program called dataplot which uses a graphical interface to visualise
1227observations and feedback files. OBSTOOLS and dataplot are described in more detail below. 
1228
1229\subsection{Obstools}
1230
1231A series of Fortran utilities is provided with NEMO called OBSTOOLS. This are helpful in
1232handling observation files and the feedback file output from the NEMO observation operator.
1233The utilities are as follows
1234
1235\subsubsection{c4comb}
1236
1237The program c4comb combines multiple class 4 files produced by individual processors in an
1238MPI run of NEMO offline obs\_oper into a single class 4 file. The program is called in the following way:
1239
1240\begin{alltt}
1241\footnotesize
1242\begin{verbatim}
1243c4comb.exe outputfile inputfile1 inputfile2 ...
1244\end{verbatim}
1245\end{alltt}
1246
1247\subsubsection{corio2fb}
1248
1249The program corio2fb converts profile observation files from the Coriolis format to the
1250standard feedback format. The program is called in the following way:
1251
1252\begin{alltt}
1253\footnotesize
1254\begin{verbatim}
1255corio2fb.exe outputfile inputfile1 inputfile2 ...
1256\end{verbatim}
1257\end{alltt}
1258
1259\subsubsection{enact2fb}
1260
1261The program enact2fb converts profile observation files from the ENACT format to the standard
1262feedback format. The program is called in the following way:
1263
1264\begin{alltt}
1265\footnotesize
1266\begin{verbatim}
1267enact2fb.exe outputfile inputfile1 inputfile2 ...
1268\end{verbatim}
1269\end{alltt}
1270
1271\subsubsection{fbcomb}
1272
1273The program fbcomb combines multiple feedback files produced by individual processors in an
1274MPI run of NEMO into a single feedback file. The program is called in the following way:
1275
1276\begin{alltt}
1277\footnotesize
1278\begin{verbatim}
1279fbcomb.exe outputfile inputfile1 inputfile2 ...
1280\end{verbatim}
1281\end{alltt}
1282
1283\subsubsection{fbmatchup}
1284
1285The program fbmatchup will match observations from two feedback files. The program is called
1286in the following way:
1287
1288\begin{alltt}
1289\footnotesize
1290\begin{verbatim}
1291fbmatchup.exe outputfile inputfile1 varname1 inputfile2 varname2 ...
1292\end{verbatim}
1293\end{alltt}
1294
1295
1296\subsubsection{fbprint}
1297
1298The program fbprint will print the contents of a feedback file or files to standard output.
1299Selected information can be output using optional arguments. The program is called in the
1300following way:
1301
1302\begin{alltt}
1303\footnotesize
1304\begin{verbatim}
1305fbprint.exe [options] inputfile
1306
1307options:
1308     -b            shorter output
1309     -q            Select observations based on QC flags
1310     -Q            Select observations based on QC flags
1311     -B            Select observations based on QC flags
1312     -u            unsorted
1313     -s ID         select station ID 
1314     -t TYPE       select observation type
1315     -v NUM1-NUM2  select variable range to print by number
1316                      (default all)
1317     -a NUM1-NUM2  select additional variable range to print by number
1318                      (default all)
1319     -e NUM1-NUM2  select extra variable range to print by number
1320                      (default all)
1321     -d            output date range
1322     -D            print depths
1323     -z            use zipped files
1324\end{verbatim}
1325\end{alltt}
1326
1327\subsubsection{fbsel}
1328
1329The program fbsel will select or subsample observations. The program is called in the
1330following way:
1331
1332\begin{alltt}
1333\footnotesize
1334\begin{verbatim}
1335fbsel.exe <input filename> <output filename>
1336\end{verbatim}
1337\end{alltt}
1338
1339\subsubsection{fbstat}
1340
1341The program fbstat will output summary statistics in different global areas into a number of
1342files. The program is called in the following way:
1343
1344\begin{alltt}
1345\footnotesize
1346\begin{verbatim}
1347fbstat.exe [-nmlev] <filenames>
1348\end{verbatim}
1349\end{alltt}
1350
1351\subsubsection{fbthin}
1352
1353The program fbthin will thin the data to 1 degree resolution. The code could easily be
1354modified to thin to a different resolution. The program is called in the following way:
1355
1356\begin{alltt}
1357\footnotesize
1358\begin{verbatim}
1359fbthin.exe inputfile outputfile
1360\end{verbatim}
1361\end{alltt}
1362
1363\subsubsection{sla2fb}
1364
1365The program sla2fb will convert an AVISO SLA format file to feedback format. The program is
1366called in the following way:
1367
1368\begin{alltt}
1369\footnotesize
1370\begin{verbatim}
1371sla2fb.exe [-s type] outputfile inputfile1 inputfile2 ...
1372
1373Option:
1374     -s            Select altimeter data_source
1375\end{verbatim}
1376\end{alltt}
1377
1378\subsubsection{vel2fb}
1379
1380The program vel2fb will convert TAO/PIRATA/RAMA currents files to feedback format. The program
1381is called in the following way:
1382
1383\begin{alltt}
1384\footnotesize
1385\begin{verbatim}
1386vel2fb.exe outputfile inputfile1 inputfile2 ...
1387\end{verbatim}
1388\end{alltt}
1389
1390\subsection{building the obstools}
1391
1392To build the obstools use in the tools directory use ./maketools -n OBSTOOLS -m [ARCH].
1393
1394\subsection{Dataplot}
1395
1396An IDL program called dataplot is included which uses a graphical interface to visualise
1397observations and feedback files. It is possible to zoom in, plot individual profiles and
1398calculate some basic statistics. To plot some data run IDL and then:
1399\begin{alltt}
1400\footnotesize
1401\begin{verbatim}
1402IDL> dataplot, "filename"
1403\end{verbatim}
1404\end{alltt}
1405
1406To read multiple files into dataplot, for example multiple feedback files from different
1407processors or from different days, the easiest method is to use the spawn command to generate
1408a list of files which can then be passed to dataplot.
1409\begin{alltt}
1410\footnotesize
1411\begin{verbatim}
1412IDL> spawn, 'ls profb*.nc', files
1413IDL> dataplot, files
1414\end{verbatim}
1415\end{alltt}
1416
1417Fig~\ref{fig:obsdataplotmain} shows the main window which is launched when dataplot starts.
1418This is split into three parts. At the top there is a menu bar which contains a variety of
1419drop down menus. Areas - zooms into prespecified regions; plot - plots the data as a
1420timeseries or a T-S diagram if appropriate; Find - allows data to be searched; Config - sets
1421various configuration options.
1422
1423The middle part is a plot of the geographical location of the observations. This will plot the
1424observation value, the model background value or observation minus background value depending
1425on the option selected in the radio button at the bottom of the window. The plotting colour
1426range can be changed by clicking on the colour bar. The title of the plot gives some basic
1427information about the date range and depth range shown, the extreme values, and the mean and
1428rms values. It is possible to zoom in using a drag-box. You may also zoom in or out using the
1429mouse wheel.
1430
1431The bottom part of the window controls what is visible in the plot above. There are two bars
1432which select the level range plotted (for profile data). The other bars below select the date
1433range shown. The bottom of the figure allows the option to plot the mean, root mean square,
1434standard deviation or mean square values. As mentioned above you can choose to plot the
1435observation value, the model background value or observation minus background value. The next
1436group of radio buttons selects the map projection. This can either be regular latitude
1437longitude grid, or north or south polar stereographic. The next group of radio buttons will
1438plot bad observations, switch to salinity and plot density for profile observations. The
1439rightmost group of buttons will print the plot window as a postscript, save it as png, or exit
1440from dataplot.
1441
1442%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1443\begin{figure}     \begin{center}
1444%\includegraphics[width=10cm,height=12cm,angle=-90.]{Fig_OBS_dataplot_main}
1445\includegraphics[width=9cm,angle=-90.]{Fig_OBS_dataplot_main}
1446\caption{      \label{fig:obsdataplotmain}
1447Main window of dataplot.}
1448\end{center}     \end{figure}
1449%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1450
1451If a profile point is clicked with the mouse button a plot of the observation and background
1452values as a function of depth (Fig~\ref{fig:obsdataplotprofile}).
1453
1454%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1455\begin{figure}     \begin{center}
1456%\includegraphics[width=10cm,height=12cm,angle=-90.]{Fig_OBS_dataplot_prof}
1457\includegraphics[width=7cm,angle=-90.]{Fig_OBS_dataplot_prof}
1458\caption{      \label{fig:obsdataplotprofile}
1459Profile plot from dataplot produced by right clicking on a point in the main window.}
1460\end{center}     \end{figure}
1461%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1462
1463
1464
1465
1466\end{document}
Note: See TracBrowser for help on using the repository browser.