source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/chap_OBS.tex @ 10419

Last change on this file since 10419 was 10419, checked in by smasson, 21 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10418, see #2133

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