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

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

source: branches/NERC/dev_r5549_BDY_ZEROGRAD/DOC/TexFiles/Chapters/Chap_OBS.tex @ 6808

Last change on this file since 6808 was 6808, checked in by jamesharle, 8 years ago

merge with trunk@6232 for consistency with SSB code

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