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/UKMO/r5518_INGV1_WAVE-coupling/DOC/TexFiles/Chapters – NEMO

source: branches/UKMO/r5518_INGV1_WAVE-coupling/DOC/TexFiles/Chapters/Chap_OBS.tex @ 7152

Last change on this file since 7152 was 7152, checked in by jcastill, 7 years ago

Initial implementation of wave coupling branch - INGV wave branch + UKMO wave coupling branch

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