source: branches/2017/dev_merge_2017/DOC/texfiles/chapters/chap_OBS.tex @ 9376

Last change on this file since 9376 was 9376, checked in by nicolasmartin, 2 years ago

Global reorganisation of DOC directory: refactoring & cleaning of TeX source code for readability

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