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/2013/dev_r3987_UKMO4_OBS/DOC/TexFiles/Chapters – NEMO

source: branches/2013/dev_r3987_UKMO4_OBS/DOC/TexFiles/Chapters/Chap_OBS.tex @ 4115

Last change on this file since 4115 was 4115, checked in by andrewryan, 11 years ago

clarified what produces class 4 files in c4comb documentation

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