Newer
Older
\documentclass[../main/NEMO_manual]{subfiles}
\begin{document}
\chapter{Observation and Model Comparison (OBS)}
\label{chap:OBS}
%\subsubsection*{Changes record}
%\begin{tabular}{l||l|m{0.65\linewidth}}
% Release & Author & Modifications \\
% {\em 4.0} & {\em D. J. Lea} & {\em \NEMO\ 4.0 updates} \\
% {\em 3.6} & {\em M. Martin, A. Ryan} & {\em Add averaging operator, standalone obs oper} \\
% {\em 3.4} & {\em D. J. Lea, M. Martin, ...} & {\em Initial version} \\
% {\em --\texttt{"}--} & {\em ... K. Mogensen, A. Vidard, A. Weaver} & {\em ---\texttt{"}---} \\
%\end{tabular}
\chaptertoc
\paragraph{Changes record} ~\\
{\footnotesize
\begin{tabularx}{\textwidth}{l||X|X}
Release & Author(s) & Modifications \\
\hline
{\em 5.0} & {\em D. Ford} & {\em Update to accommodate interface changes} \\
{\em 4.0} & {\em ...} & {\em ...} \\
{\em 3.6} & {\em ...} & {\em ...} \\
{\em 3.4} & {\em ...} & {\em ...} \\
{\em <=3.4} & {\em ...} & {\em ...}
\end{tabularx}
}
\clearpage
The observation and model comparison code, the observation operator (OBS), reads in observation files
(temperature profiles, salinity profiles, velocity profiles, sea surface temperature, sea surface salinity, sea level anomaly, sea surface velocity, sea ice concentration, sea ice thickness)
and calculates an interpolated model equivalent value at the observation location and nearest model time step.
The resulting data are saved in a ``feedback'' file (or files).
The code was originally developed for use with the NEMOVAR data assimilation code,
but can be used for validation or verification of the model or with any other data assimilation system.
Its use can be extended to other variables, for instance biogeochemistry, with small changes to the code.
The OBS code is called from \mdl{nemogcm} for model initialisation and to calculate the model equivalent values for observations on the 0th time step.
The code is then called again after each time step from \mdl{step}.
The code is only activated if the \nam{obs}{obs} namelist logical \np{ln_diaobs}{ln\_diaobs} is set to true.
For all data types a 2D horizontal interpolator or averager is needed to
interpolate/average the model fields to the observation location.
For {\em in situ} profiles, a 1D vertical interpolator is needed in addition to
provide model fields at the observation depths.
This now works in a generalised vertical coordinate system.
Some profile observation types (\eg\ tropical moored buoys) are made available as daily averaged quantities.
The observation operator code can be set-up to calculate the equivalent daily average model temperature fields using
the \np{nn_profdavtypes}{nn\_profdavtypes} namelist array.
Some SST observations are equivalent to a night-time average value and
the observation operator code can calculate equivalent night-time average model fields by
setting the namelist value \np{ln_night}{ln\_night} to true.
Otherwise (by default) the model value from the nearest time step to the observation time is used.
The code is controlled by the namelist \nam{obs}{obs} and multiple instances of the namelist \nam{obs_dta}{obs\_dta}, one for each observation group.
See the following sections for more details on setting up the namelists.
In \autoref{sec:OBS_example} a test example of the observation operator code is introduced, including
where to obtain data and how to setup the namelists.
In \autoref{sec:OBS_advanced} some more technical details of the available options are introduced, and we
also show more complete namelists.
In \autoref{sec:OBS_headers} example ``feedback'' file NetCDF headers are shown.
In \autoref{sec:OBS_newvars} a description of how to add an observation operator for a new variable is given.
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
In \autoref{sec:OBS_theory} some of the theoretical aspects of the observation operator are described including
interpolation methods and running on multiple processors.
In \autoref{sec:OBS_sao} the standalone observation operator code is described.
In \autoref{sec:OBS_obsutils} we describe some utilities to help work with the files produced by the OBS code.
%% =================================================================================================
\section{Running the observation operator code example}
\label{sec:OBS_example}
In this section an example of running the observation operator code is described using
profile observation data which can be freely downloaded.
It shows how to adapt an existing run and build of \NEMO\ to run the observation operator. Note also the observation operator and the assimilation increments code are run in the ORCA2\_ICE\_OBS SETTE test.
\begin{enumerate}
\item Compile \NEMO.
\item Download some EN4 data from \href{http://www.metoffice.gov.uk/hadobs}{www.metoffice.gov.uk/hadobs}.
Choose observations which are valid for the period of your test run because
the observation operator compares the model and observations for a matching date and time.
\item Compile the OBSTOOLS code in the \path{tools} directory using:
\begin{cmds}
./maketools -n OBSTOOLS -m [ARCH]
\end{cmds}
replacing \texttt{[ARCH]} with the build architecture file for your machine. Note the tools are checked out from a separate location of the repository (under \path{/utils/tools}).
\item Convert the EN4 data into feedback format:
\begin{cmds}
enact2fb.exe profiles_01.nc EN.4.1.1.f.profiles.g10.YYYYMM.nc
\end{cmds}
\item Include the following in the \NEMO\ namelist to run the observation operator on this data:
\end{enumerate}
Top-level options are defined through the \nam{obs}{obs} namelist variables.
\begin{itemize}
\item \np{ln_diaobs}{ln\_diaobs} turns on the observation operator code.
\item \np{nn_obsgroups}{nn\_obsgroups} specifies the number of observation groups and therefore number of \nam{obs_dta}{obs\_dta} namelists required; in this example it should be 1, as the temperature and salinity data are grouped together in the same input file(s).
\item The model grid points for a particular observation latitude and longitude are found using
the grid searching part of the code.
This can be expensive, particularly for large numbers of observations,
setting \np{ln_grid_search_lookup}{ln\_grid\_search\_lookup} allows the use of a lookup table which
is saved into an \np{cn_gridsearchfile}{cn\_gridsearchfile} file (or files).
This will need to be generated the first time if it does not exist in the run directory.
However, once produced it will significantly speed up future grid searches.
\np{rn_gridsearchres}{rn\_gridsearchres} specifies the resolution (in degrees) of the grid search.
\item Setting \np{ln_grid_global}{ln\_grid\_global} means that the code distributes the observations evenly between processors.
Alternatively each processor will work with observations located within the model subdomain
(see \autoref{subsec:OBS_parallel}).
\item \np{rn_dobsini}{rn\_dobsini} and \np{rn_dobsend}{rn\_dobsend} can be used to restrict the period the observation operator is run for.
\end{itemize}
A \nam{obs_dta}{obs\_dta} namelist is then needed for each observation group.
Typically, an individual variable (\eg\ sea surface temperature) would be its own group, but if variables are collocated in the input files (\eg\ temperature and salinity profiles, or meridional and zonal velocities) they would be considered as a single group.
\begin{itemize}
\item \np{cn_groupname}{cn\_groupname} is used to name the output ``feedback'' file(s).
\item \np{ln_enabled}{ln\_enabled} must be .true. for the observation operator to be run on this group, setting it to .false. allows it to be turned off without needing to remove the instance of \nam{obs_dta}{obs\_dta}.
\item Setting \np{ln_prof}{ln\_prof} specifies the observations consist of vertical profiles, as in this example.
\np{ln_surf}{ln\_surf} indicates surface data.
One and only one of \np{ln_prof}{ln\_prof} and \np{ln_surf}{ln\_surf} must be .true..
\item \np{cn_obsfiles}{cn\_obsfiles} specifies the input filename or array of filenames.
\item \np{cn_obstypes}{cn\_obstypes} specifies the observation variable(s).
In this example it should be set to 'POTM','PSAL' for profile temperature and profile salinity.
A full list of \np{cn_obstypes}{cn\_obstypes} options is given in \autoref{sec:OBS_headers}.
\end{itemize}
Other \nam{obs_dta}{obs\_dta} options are available to control specific behaviour, as detailed in \autoref{sec:OBS_advanced}.
A number of utilities are now provided to plot the feedback files, convert and recombine the files.
These are explained in more detail in \autoref{sec:OBS_obsutils}.
Utilities to convert other input data formats into the feedback format are also described in
\autoref{sec:OBS_obsutils}.
%% =================================================================================================
\section{Technical details and full namelist options}
\label{sec:OBS_advanced}
Here we show more complete information about the available options in the \nam{obs_dta}{obs\_dta} namelists.
Each item in \nam{obs}{obs} has already been described in \autoref{sec:OBS_example}.
The first six entries in \nam{obs_dta}{obs\_dta}, which are the ones most likely to be modified by users,
have also been described in \autoref{sec:OBS_example}.
\begin{listing}
\nlst{namobs}
\caption{\forcode{&namobs}}
\label{lst:namobs}
\end{listing}
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
\begin{listing}
\nlst{namobs_dta}
\caption{\forcode{&namobs_dta}}
\label{lst:namobs_dta}
\end{listing}
The remaining options in \nam{obs_dta}{obs\_dta} offer more specialist control, with some only
relevant to specific observation types.
The following are relevant to all observation types:
\begin{itemize}
\item Setting \np{ln_nea}{ln\_nea} rejects observations within a grid cell that neighbours land.
\item Setting \np{ln_bound_reject}{ln\_bound\_reject} rejects observations within a grid cell that neighbours a boundary point.
\item Setting \np{ln_ignmis}{ln\_ignmis} merely issues a warning if a file listed in \np{cn_obsfiles}{cn\_obsfiles} is missing,
otherwise execution of the code ceases with an error.
\item \np{nn_2dint}{nn\_2dint} sets the type of horizontal interpolation method, see \autoref{sec:OBS_theory} for details.
\end{itemize}
The following are only relevant to profile observations:
\begin{itemize}
\item \np{nn_1dint}{nn\_1dint} sets the type of vertical interpolation method, see \autoref{sec:OBS_theory} for details.
\item \np{nn_profdavtypes}{nn\_profdavtypes} is an array of profile observation types corresponding to entries in the
STATION\_TYPE variable in the feedback file (see \autoref{sec:OBS_headers}) which should be treated as
daily averages rather than at a specific time.
Set to -1 if no daily averaged types.
\item Setting \np{ln_all_at_all}{ln\_all\_at\_all} means that if there is more than one variable in
\np{cn_obstypes}{cn\_obstypes}, then a model counterpart is calculated for all variables at any location
for which there is an observation of any of the variables.
\end{itemize}
The following are only relevant to surface observations:
\begin{itemize}
\item Setting \np{ln_fp_indegs}{ln\_fp\_indegs} means that the footprint of the observations (see options below)
is in units of degrees rather than metres.
\item \np{rn_avglamscl}{rn\_avglamscl} and \np{rn_avgphiscl}{rn\_avgphiscl} specify the east/west and north/south
footprint of the observations respectively.
This is only relevant if \np{nn_2dint}{nn\_2dint} is set to 5 or 6.
See \autoref{sec:OBS_theory} for details.
\item Setting \np{ln_night}{ln\_night} calculates a night-time average model counterpart rather than one at
the nearest model time step.
\item Setting \np{ln_time_mean_bkg}{ln\_time\_mean\_bkg} calculates a model counterpart averaged over
\np{rn_time_mean_period}{rn\_time\_mean\_period} (in hours) rather than one at the nearest model time step.
\item Setting \np{ln_obsbias}{ln\_obsbias} allows a NetCDF file containing an estimate of the observation bias
on the model grid to be read in, and this bias to be removed from the observation values.
\np{cn_obsbiasfiles}{cn\_obsbiasfiles} is a list of observation bias file names.
\np{cn_type_to_biascorrect}{cn\_type\_to\_biascorrect} is the variable in \np{cn_obstypes}{cn\_obstypes} to bias correct.
\np{cn_obsbiasfile_varname}{cn\_obsbiasfile\_varname} is the NetCDF variable to read from \np{cn_obsbiasfiles}{cn\_obsbiasfiles}.
\end{itemize}
The following are only relevant to sea level anomaly observations ('SLA' is in \np{cn_obstypes}{cn\_obstypes}),
which require use of a mean dynamic topography (MDT) (see \autoref{sec:OBS_headers}) which may need bias correcting:
\begin{itemize}
\item Setting \np{ln_altbias}{ln\_altbias} means that the MDT used is bias corrected using a bias file on the
model grid in the NetCDF file \np{cn_altbiasfile}{cn\_altbiasfile}.
\item \np{nn_msshc}{nn\_msshc} sets the mean sea surface height (MSSH) correction to apply to the MDT.
0 means no correction is applied.
1 means a correction based on the mean difference between the MDT and model SSH is applied.
2 means a correction of \np{rn_mdtcorr}{rn\_mdtcorr} is applied.
Either correction is only applied at latitudes between \np{rn_mdtcutoff}{rn\_mdtcutoff} degrees north and
\np{rn_mdtcutoff}{rn\_mdtcutoff} degrees south.
\end{itemize}
The following is only relevant to temperature and salinity observations ('POTM', 'PSAL', 'SST', or 'SSS' in
\np{cn_obstypes}{cn\_obstypes}):
\begin{itemize}
\item Setting \np{ln_output_clim}{ln\_output\_clim} means that as well as a model counterpart, an equivalent
climatological counterpart is output.
The climatology used is the same as used for tracer damping (see \autoref{sec:TRA_dmp}).
Therefore, this option also requires \np{ln_tradmp}{ln\_tradmp} to be set in \nam{tra_dmp}{tra\_dmp}.
\end{itemize}
The following is only relevant to sea ice freeboard observations ('FBD' is in \np{cn_obstypes}{cn\_obstypes}):
\begin{itemize}
\item \np{rn_radar_snow_penetr}{rn\_radar\_snow\_penetr} sets the snow depth penetration factor for the
radar ice freeboard conversion.
\end{itemize}
%% =================================================================================================
\section{Example feedback type observation file headers}
\label{sec:OBS_headers}
The observation operator code uses the feedback observation file format for all data types.
All the observation files must be in NetCDF format.
Some example headers (produced using \mbox{\textit{ncdump~-h}}) for temperature and salinity profile data,
sea level anomaly and sea surface temperature are in the following subsections.
The VARIABLES in the feedback files must match those specified in \np{cn_obstypes}{cn\_obstypes} in \nam{obs_dta}{obs\_dta}.
The following variable names can be used:
\begin{itemize}
\item POTM: temperature profiles
\item PSAL: salinity profiles
\item UVEL: zonal velocities (profile or surface)
\item VVEL: meridional velocities (profile or surface)
\item SST: sea surface temperature
\item SLA: sea level anomaly
\item SSS: sea surface salinity
\item ICECONC: sea ice concentration
\item SIT: sea ice thickness
\item FBD: sea ice freeboard
\end{itemize}
Other variables, such as biogeochemical variables from PISCES, are not currently available.
However, to define a new variable a user need only add an extra option to the relevant SELECT CASE statement in the routine dia\_obs, that matches a new observation type with the desired model variable (see \autoref{sec:OBS_newvars}).
%% =================================================================================================
\subsection{Temperature and salinity profile feedback file}
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
\begin{clines}
netcdf profiles_01 {
dimensions:
N_OBS = 603 ;
N_LEVELS = 150 ;
N_VARS = 2 ;
N_QCF = 2 ;
N_ENTRIES = 1 ;
N_EXTRA = 1 ;
STRINGNAM = 8 ;
STRINGGRID = 1 ;
STRINGWMO = 8 ;
STRINGTYP = 4 ;
STRINGJULD = 14 ;
variables:
char VARIABLES(N_VARS, STRINGNAM) ;
VARIABLES:long_name = "List of variables in feedback files" ;
char ENTRIES(N_ENTRIES, STRINGNAM) ;
ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
char EXTRA(N_EXTRA, STRINGNAM) ;
EXTRA:long_name = "List of extra variables" ;
char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
STATION_IDENTIFIER:long_name = "Station identifier" ;
char STATION_TYPE(N_OBS, STRINGTYP) ;
STATION_TYPE:long_name = "Code instrument type" ;
double LONGITUDE(N_OBS) ;
LONGITUDE:long_name = "Longitude" ;
LONGITUDE:units = "degrees_east" ;
LONGITUDE:_Fillvalue = 99999.f ;
double LATITUDE(N_OBS) ;
LATITUDE:long_name = "Latitude" ;
LATITUDE:units = "degrees_north" ;
LATITUDE:_Fillvalue = 99999.f ;
double DEPTH(N_OBS, N_LEVELS) ;
DEPTH:long_name = "Depth" ;
DEPTH:units = "metre" ;
DEPTH:_Fillvalue = 99999.f ;
int DEPTH_QC(N_OBS, N_LEVELS) ;
DEPTH_QC:long_name = "Quality on depth" ;
DEPTH_QC:Conventions = "q where q =[0,9]" ;
DEPTH_QC:_Fillvalue = 0 ;
int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
double JULD(N_OBS) ;
JULD:long_name = "Julian day" ;
JULD:units = "days since JULD_REFERENCE" ;
JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
JULD:_Fillvalue = 99999.f ;
char JULD_REFERENCE(STRINGJULD) ;
JULD_REFERENCE:long_name = "Date of reference for julian days" ;
JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
int OBSERVATION_QC(N_OBS) ;
OBSERVATION_QC:long_name = "Quality on observation" ;
OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
OBSERVATION_QC:_Fillvalue = 0 ;
int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
int POSITION_QC(N_OBS) ;
POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
POSITION_QC:Conventions = "q where q =[0,9]" ;
POSITION_QC:_Fillvalue = 0 ;
int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
POSITION_QC_FLAGS:_Fillvalue = 0 ;
int JULD_QC(N_OBS) ;
JULD_QC:long_name = "Quality on date and time" ;
JULD_QC:Conventions = "q where q =[0,9]" ;
JULD_QC:_Fillvalue = 0 ;
int JULD_QC_FLAGS(N_OBS, N_QCF) ;
JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
JULD_QC_FLAGS:_Fillvalue = 0 ;
int ORIGINAL_FILE_INDEX(N_OBS) ;
ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
float POTM_OBS(N_OBS, N_LEVELS) ;
POTM_OBS:long_name = "Potential temperature" ;
POTM_OBS:units = "Degrees Celsius" ;
POTM_OBS:_Fillvalue = 99999.f ;
float POTM_Hx(N_OBS, N_LEVELS) ;
POTM_Hx:long_name = "Model interpolated potential temperature" ;
POTM_Hx:units = "Degrees Celsius" ;
POTM_Hx:_Fillvalue = 99999.f ;
int POTM_QC(N_OBS) ;
POTM_QC:long_name = "Quality on potential temperature" ;
POTM_QC:Conventions = "q where q =[0,9]" ;
POTM_QC:_Fillvalue = 0 ;
int POTM_QC_FLAGS(N_OBS, N_QCF) ;
POTM_QC_FLAGS:long_name = "Quality flags on potential temperature" ;
POTM_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
POTM_QC_FLAGS:_Fillvalue = 0 ;
int POTM_LEVEL_QC(N_OBS, N_LEVELS) ;
POTM_LEVEL_QC:long_name = "Quality for each level on potential temperature" ;
POTM_LEVEL_QC:Conventions = "q where q =[0,9]" ;
POTM_LEVEL_QC:_Fillvalue = 0 ;
int POTM_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
POTM_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on potential temperature" ;
POTM_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
POTM_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
int POTM_IOBSI(N_OBS) ;
POTM_IOBSI:long_name = "ORCA grid search I coordinate" ;
int POTM_IOBSJ(N_OBS) ;
POTM_IOBSJ:long_name = "ORCA grid search J coordinate" ;
int POTM_IOBSK(N_OBS, N_LEVELS) ;
POTM_IOBSK:long_name = "ORCA grid search K coordinate" ;
char POTM_GRID(STRINGGRID) ;
POTM_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
float PSAL_OBS(N_OBS, N_LEVELS) ;
PSAL_OBS:long_name = "Practical salinity" ;
PSAL_OBS:units = "PSU" ;
PSAL_OBS:_Fillvalue = 99999.f ;
float PSAL_Hx(N_OBS, N_LEVELS) ;
PSAL_Hx:long_name = "Model interpolated practical salinity" ;
PSAL_Hx:units = "PSU" ;
PSAL_Hx:_Fillvalue = 99999.f ;
int PSAL_QC(N_OBS) ;
PSAL_QC:long_name = "Quality on practical salinity" ;
PSAL_QC:Conventions = "q where q =[0,9]" ;
PSAL_QC:_Fillvalue = 0 ;
int PSAL_QC_FLAGS(N_OBS, N_QCF) ;
PSAL_QC_FLAGS:long_name = "Quality flags on practical salinity" ;
PSAL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
PSAL_QC_FLAGS:_Fillvalue = 0 ;
int PSAL_LEVEL_QC(N_OBS, N_LEVELS) ;
PSAL_LEVEL_QC:long_name = "Quality for each level on practical salinity" ;
PSAL_LEVEL_QC:Conventions = "q where q =[0,9]" ;
PSAL_LEVEL_QC:_Fillvalue = 0 ;
int PSAL_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
PSAL_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on practical salinity" ;
PSAL_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
PSAL_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
int PSAL_IOBSI(N_OBS) ;
PSAL_IOBSI:long_name = "ORCA grid search I coordinate" ;
int PSAL_IOBSJ(N_OBS) ;
PSAL_IOBSJ:long_name = "ORCA grid search J coordinate" ;
int PSAL_IOBSK(N_OBS, N_LEVELS) ;
PSAL_IOBSK:long_name = "ORCA grid search K coordinate" ;
char PSAL_GRID(STRINGGRID) ;
PSAL_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
float TEMP(N_OBS, N_LEVELS) ;
TEMP:long_name = "Insitu temperature" ;
TEMP:units = "Degrees Celsius" ;
TEMP:_Fillvalue = 99999.f ;
// global attributes:
:title = "NEMO observation operator output" ;
:Convention = "NEMO unified observation operator output" ;
}
\end{clines}
%% =================================================================================================
\subsection{Sea level anomaly feedback file}
\begin{clines}
netcdf sla_01 {
dimensions:
N_OBS = 41301 ;
N_LEVELS = 1 ;
N_VARS = 1 ;
N_QCF = 2 ;
N_ENTRIES = 1 ;
N_EXTRA = 1 ;
STRINGNAM = 8 ;
STRINGGRID = 1 ;
STRINGWMO = 8 ;
STRINGTYP = 4 ;
STRINGJULD = 14 ;
variables:
char VARIABLES(N_VARS, STRINGNAM) ;
VARIABLES:long_name = "List of variables in feedback files" ;
char ENTRIES(N_ENTRIES, STRINGNAM) ;
ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
char EXTRA(N_EXTRA, STRINGNAM) ;
EXTRA:long_name = "List of extra variables" ;
char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
STATION_IDENTIFIER:long_name = "Station identifier" ;
char STATION_TYPE(N_OBS, STRINGTYP) ;
STATION_TYPE:long_name = "Code instrument type" ;
double LONGITUDE(N_OBS) ;
LONGITUDE:long_name = "Longitude" ;
LONGITUDE:units = "degrees_east" ;
LONGITUDE:_Fillvalue = 99999.f ;
double LATITUDE(N_OBS) ;
LATITUDE:long_name = "Latitude" ;
LATITUDE:units = "degrees_north" ;
LATITUDE:_Fillvalue = 99999.f ;
double DEPTH(N_OBS, N_LEVELS) ;
DEPTH:long_name = "Depth" ;
DEPTH:units = "metre" ;
DEPTH:_Fillvalue = 99999.f ;
int DEPTH_QC(N_OBS, N_LEVELS) ;
DEPTH_QC:long_name = "Quality on depth" ;
DEPTH_QC:Conventions = "q where q =[0,9]" ;
DEPTH_QC:_Fillvalue = 0 ;
int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
double JULD(N_OBS) ;
JULD:long_name = "Julian day" ;
JULD:units = "days since JULD_REFERENCE" ;
JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
JULD:_Fillvalue = 99999.f ;
char JULD_REFERENCE(STRINGJULD) ;
JULD_REFERENCE:long_name = "Date of reference for julian days" ;
JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
int OBSERVATION_QC(N_OBS) ;
OBSERVATION_QC:long_name = "Quality on observation" ;
OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
OBSERVATION_QC:_Fillvalue = 0 ;
int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
int POSITION_QC(N_OBS) ;
POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
POSITION_QC:Conventions = "q where q =[0,9]" ;
POSITION_QC:_Fillvalue = 0 ;
int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
POSITION_QC_FLAGS:_Fillvalue = 0 ;
int JULD_QC(N_OBS) ;
JULD_QC:long_name = "Quality on date and time" ;
JULD_QC:Conventions = "q where q =[0,9]" ;
JULD_QC:_Fillvalue = 0 ;
int JULD_QC_FLAGS(N_OBS, N_QCF) ;
JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
JULD_QC_FLAGS:_Fillvalue = 0 ;
int ORIGINAL_FILE_INDEX(N_OBS) ;
ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
float SLA_OBS(N_OBS, N_LEVELS) ;
SLA_OBS:long_name = "Sea level anomaly" ;
SLA_OBS:units = "metre" ;
SLA_OBS:_Fillvalue = 99999.f ;
float SLA_Hx(N_OBS, N_LEVELS) ;
SLA_Hx:long_name = "Model interpolated sea level anomaly" ;
SLA_Hx:units = "metre" ;
SLA_Hx:_Fillvalue = 99999.f ;
int SLA_QC(N_OBS) ;
SLA_QC:long_name = "Quality on sea level anomaly" ;
SLA_QC:Conventions = "q where q =[0,9]" ;
SLA_QC:_Fillvalue = 0 ;
int SLA_QC_FLAGS(N_OBS, N_QCF) ;
SLA_QC_FLAGS:long_name = "Quality flags on sea level anomaly" ;
SLA_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
SLA_QC_FLAGS:_Fillvalue = 0 ;
int SLA_LEVEL_QC(N_OBS, N_LEVELS) ;
SLA_LEVEL_QC:long_name = "Quality for each level on sea level anomaly" ;
SLA_LEVEL_QC:Conventions = "q where q =[0,9]" ;
SLA_LEVEL_QC:_Fillvalue = 0 ;
int SLA_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
SLA_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on sea level anomaly" ;
SLA_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
SLA_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
int SLA_IOBSI(N_OBS) ;
SLA_IOBSI:long_name = "ORCA grid search I coordinate" ;
int SLA_IOBSJ(N_OBS) ;
SLA_IOBSJ:long_name = "ORCA grid search J coordinate" ;
int SLA_IOBSK(N_OBS, N_LEVELS) ;
SLA_IOBSK:long_name = "ORCA grid search K coordinate" ;
char SLA_GRID(STRINGGRID) ;
SLA_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
float MDT(N_OBS, N_LEVELS) ;
MDT:long_name = "Mean Dynamic Topography" ;
MDT:units = "metre" ;
MDT:_Fillvalue = 99999.f ;
// global attributes:
:title = "NEMO observation operator output" ;
:Convention = "NEMO unified observation operator output" ;
}
\end{clines}
To use Sea Level Anomaly (SLA) data the mean dynamic topography (MDT) must be provided in a separate file defined on
the model grid called \textit{slaReferenceLevel.nc}.
The MDT is required in order to produce the model equivalent sea level anomaly from the model sea surface height.
Below is an example header for this file (on the ORCA025 grid).
\begin{clines}
dimensions:
x = 1442 ;
y = 1021 ;
variables:
float nav_lon(y, x) ;
nav_lon:units = "degrees_east" ;
float nav_lat(y, x) ;
nav_lat:units = "degrees_north" ;
float sossheig(y, x) ;
sossheig:_FillValue = -1.e+30f ;
sossheig:coordinates = "nav_lon nav_lat" ;
sossheig:long_name = "Mean Dynamic Topography" ;
sossheig:units = "metres" ;
sossheig:grid = "orca025T" ;
\end{clines}
%% =================================================================================================
\subsection{Sea surface temperature feedback file}
\begin{clines}
netcdf sst_01 {
dimensions:
N_OBS = 33099 ;
N_LEVELS = 1 ;
N_VARS = 1 ;
N_QCF = 2 ;
N_ENTRIES = 1 ;
STRINGNAM = 8 ;
STRINGGRID = 1 ;
STRINGWMO = 8 ;
STRINGTYP = 4 ;
STRINGJULD = 14 ;
variables:
char VARIABLES(N_VARS, STRINGNAM) ;
VARIABLES:long_name = "List of variables in feedback files" ;
char ENTRIES(N_ENTRIES, STRINGNAM) ;
ENTRIES:long_name = "List of additional entries for each variable in feedback files" ;
char STATION_IDENTIFIER(N_OBS, STRINGWMO) ;
STATION_IDENTIFIER:long_name = "Station identifier" ;
char STATION_TYPE(N_OBS, STRINGTYP) ;
STATION_TYPE:long_name = "Code instrument type" ;
double LONGITUDE(N_OBS) ;
LONGITUDE:long_name = "Longitude" ;
LONGITUDE:units = "degrees_east" ;
LONGITUDE:_Fillvalue = 99999.f ;
double LATITUDE(N_OBS) ;
LATITUDE:long_name = "Latitude" ;
LATITUDE:units = "degrees_north" ;
LATITUDE:_Fillvalue = 99999.f ;
double DEPTH(N_OBS, N_LEVELS) ;
DEPTH:long_name = "Depth" ;
DEPTH:units = "metre" ;
DEPTH:_Fillvalue = 99999.f ;
int DEPTH_QC(N_OBS, N_LEVELS) ;
DEPTH_QC:long_name = "Quality on depth" ;
DEPTH_QC:Conventions = "q where q =[0,9]" ;
DEPTH_QC:_Fillvalue = 0 ;
int DEPTH_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
DEPTH_QC_FLAGS:long_name = "Quality flags on depth" ;
DEPTH_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
double JULD(N_OBS) ;
JULD:long_name = "Julian day" ;
JULD:units = "days since JULD_REFERENCE" ;
JULD:Conventions = "relative julian days with decimal part (as parts of day)" ;
JULD:_Fillvalue = 99999.f ;
char JULD_REFERENCE(STRINGJULD) ;
JULD_REFERENCE:long_name = "Date of reference for julian days" ;
JULD_REFERENCE:Conventions = "YYYYMMDDHHMMSS" ;
int OBSERVATION_QC(N_OBS) ;
OBSERVATION_QC:long_name = "Quality on observation" ;
OBSERVATION_QC:Conventions = "q where q =[0,9]" ;
OBSERVATION_QC:_Fillvalue = 0 ;
int OBSERVATION_QC_FLAGS(N_OBS, N_QCF) ;
OBSERVATION_QC_FLAGS:long_name = "Quality flags on observation" ;
OBSERVATION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
OBSERVATION_QC_FLAGS:_Fillvalue = 0 ;
int POSITION_QC(N_OBS) ;
POSITION_QC:long_name = "Quality on position (latitude and longitude)" ;
POSITION_QC:Conventions = "q where q =[0,9]" ;
POSITION_QC:_Fillvalue = 0 ;
int POSITION_QC_FLAGS(N_OBS, N_QCF) ;
POSITION_QC_FLAGS:long_name = "Quality flags on position" ;
POSITION_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
POSITION_QC_FLAGS:_Fillvalue = 0 ;
int JULD_QC(N_OBS) ;
JULD_QC:long_name = "Quality on date and time" ;
JULD_QC:Conventions = "q where q =[0,9]" ;
JULD_QC:_Fillvalue = 0 ;
int JULD_QC_FLAGS(N_OBS, N_QCF) ;
JULD_QC_FLAGS:long_name = "Quality flags on date and time" ;
JULD_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
JULD_QC_FLAGS:_Fillvalue = 0 ;
int ORIGINAL_FILE_INDEX(N_OBS) ;
ORIGINAL_FILE_INDEX:long_name = "Index in original data file" ;
ORIGINAL_FILE_INDEX:_Fillvalue = -99999 ;
float SST_OBS(N_OBS, N_LEVELS) ;
SST_OBS:long_name = "Sea surface temperature" ;
SST_OBS:units = "Degree centigrade" ;
SST_OBS:_Fillvalue = 99999.f ;
float SST_Hx(N_OBS, N_LEVELS) ;
SST_Hx:long_name = "Model interpolated sea surface temperature" ;
SST_Hx:units = "Degree centigrade" ;
SST_Hx:_Fillvalue = 99999.f ;
int SST_QC(N_OBS) ;
SST_QC:long_name = "Quality on sea surface temperature" ;
SST_QC:Conventions = "q where q =[0,9]" ;
SST_QC:_Fillvalue = 0 ;
int SST_QC_FLAGS(N_OBS, N_QCF) ;
SST_QC_FLAGS:long_name = "Quality flags on sea surface temperature" ;
SST_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
SST_QC_FLAGS:_Fillvalue = 0 ;
int SST_LEVEL_QC(N_OBS, N_LEVELS) ;
SST_LEVEL_QC:long_name = "Quality for each level on sea surface temperature" ;
SST_LEVEL_QC:Conventions = "q where q =[0,9]" ;
SST_LEVEL_QC:_Fillvalue = 0 ;
int SST_LEVEL_QC_FLAGS(N_OBS, N_LEVELS, N_QCF) ;
SST_LEVEL_QC_FLAGS:long_name = "Quality flags for each level on sea surface temperature" ;
SST_LEVEL_QC_FLAGS:Conventions = "NEMOVAR flag conventions" ;
SST_LEVEL_QC_FLAGS:_Fillvalue = 0 ;
int SST_IOBSI(N_OBS) ;
SST_IOBSI:long_name = "ORCA grid search I coordinate" ;
int SST_IOBSJ(N_OBS) ;
SST_IOBSJ:long_name = "ORCA grid search J coordinate" ;
int SST_IOBSK(N_OBS, N_LEVELS) ;
SST_IOBSK:long_name = "ORCA grid search K coordinate" ;
char SST_GRID(STRINGGRID) ;
SST_GRID:long_name = "ORCA grid search grid (T,U,V)" ;
// global attributes:
:title = "NEMO observation operator output" ;
:Convention = "NEMO unified observation operator output" ;
}
\end{clines}
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
%% =================================================================================================
\section{Adding code for a new variable}
\label{sec:OBS_newvars}
As detailed in \autoref{sec:OBS_headers}, the OBS interface is only defined for specific physics variables.
However, the internal observation operator code is largely generic, and for observation types which can be
directly matched to model variables, is relatively straightforward to implement for additional variables.
This could include biogeochemical variables from PISCES or another coupled model.
In the routine dia\_obs in \mdl{diaobs} (\path{src/OCE/OBS/diaobs.F90}) there are two \textbf{SELECT CASE}
constructs, one for surface variables and one for profile variables.
To implement an observation operator for a new variable, an extra \textbf{CASE} need simply be added to the
relevant \textbf{SELECT CASE} construct.
To give a generic example for a profile variable:
\begin{forlines}
CASE('OBSNAME')
zprofvar(:,:,:) = MODELVARIABLE(:,:,:,Kmm)
\end{forlines}
where MODELVARIABLE is the model variable to be matched to the observation type, and OBSNAME is the name of
the observation variable used in the feedback file and in \np{cn_obstypes}{cn\_obstypes} in the
\nam{obs_dta}{obs\_dta} namelist.
OBSNAME must have a maximum of eight characters.
Then, a corresponding namelist and feedback file(s) just need to be supplied, following the method
described in \autoref{sec:OBS_example}.
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
%% =================================================================================================
\section{Theoretical details}
\label{sec:OBS_theory}
%% =================================================================================================
\subsection{Horizontal interpolation and averaging methods}
For most observation types, the horizontal extent of the observation is small compared to the model grid size and so
the model equivalent of the observation is calculated by interpolating from
the four surrounding grid points to the observation location.
Some satellite observations (\eg\ microwave satellite SST data, or satellite SSS data) have a footprint which
is similar in size or larger than the model grid size (particularly when the grid size is small).
In those cases the model counterpart should be calculated by averaging the model grid points over
the same size as the footprint.
\NEMO\ therefore has the capability to specify either an interpolation or an averaging
(for surface observation types only).
The main namelist option associated with the interpolation/averaging is \np{nn_2dint}{nn\_2dint}.
This default option can be set to values from 0 to 6.
Values between 0 to 4 are associated with interpolation while values 5 or 6 are associated with averaging.
\begin{itemize}
\item \np[=0]{nn_2dint}{nn\_2dint}: Distance-weighted interpolation
\item \np[=1]{nn_2dint}{nn\_2dint}: Distance-weighted interpolation (small angle)
\item \np[=2]{nn_2dint}{nn\_2dint}: Bilinear interpolation (geographical grid)
\item \np[=3]{nn_2dint}{nn\_2dint}: Bilinear remapping interpolation (general grid)
\item \np[=4]{nn_2dint}{nn\_2dint}: Polynomial interpolation
\item \np[=5]{nn_2dint}{nn\_2dint}: Radial footprint averaging with diameter specified in the namelist as
\texttt{rn\_avglamscl} in degrees or metres (set using \texttt{ln\_fp\_indegs})
\item \np[=6]{nn_2dint}{nn\_2dint}: Rectangular footprint averaging with E/W and N/S size specified in
the namelist as \texttt{rn\_avglamscl} and \texttt{rn\_avgphiscl} in degrees or metres
(set using \texttt{ln\_fp\_indegs})
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
\end{itemize}
Below is some more detail on the various options for interpolation and averaging available in \NEMO.
%% =================================================================================================
\subsubsection{Horizontal interpolation}
Consider an observation point ${\mathrm P}$ with longitude and latitude (${\lambda_{}}_{\mathrm P}$, $\phi_{\mathrm P}$) and
the four nearest neighbouring model grid points ${\mathrm A}$, ${\mathrm B}$, ${\mathrm C}$ and ${\mathrm D}$ with
longitude and latitude ($\lambda_{\mathrm A}$, $\phi_{\mathrm A}$),($\lambda_{\mathrm B}$, $\phi_{\mathrm B}$) etc.
All horizontal interpolation methods implemented in \NEMO\ estimate the value of a model variable $x$ at point $P$ as
a weighted linear combination of the values of the model variables at the grid points ${\mathrm A}$, ${\mathrm B}$ etc.:
\begin{align*}
{x_{}}_{\mathrm P} =
\frac{1}{w} \left( {w_{}}_{\mathrm A} {x_{}}_{\mathrm A} +
{w_{}}_{\mathrm B} {x_{}}_{\mathrm B} +
{w_{}}_{\mathrm C} {x_{}}_{\mathrm C} +
{w_{}}_{\mathrm D} {x_{}}_{\mathrm D} \right)
\end{align*}
where ${w_{}}_{\mathrm A}$, ${w_{}}_{\mathrm B}$ etc. are the respective weights for the model field at
points ${\mathrm A}$, ${\mathrm B}$ etc., and $w = {w_{}}_{\mathrm A} + {w_{}}_{\mathrm B} + {w_{}}_{\mathrm C} + {w_{}}_{\mathrm D}$.
Four different possibilities are available for computing the weights.
\begin{enumerate}
\item {\bfseries Great-Circle distance-weighted interpolation.}
The weights are computed as a function of the great-circle distance $s(P, \cdot)$ between $P$ and
the model grid points $A$, $B$ etc.
For example, the weight given to the field ${x_{}}_{\mathrm A}$ is specified as the product of the distances
from ${\mathrm P}$ to the other points:
\begin{alignat*}{2}
{w_{}}_{\mathrm A} = s({\mathrm P}, {\mathrm B}) \, s({\mathrm P}, {\mathrm C}) \, s({\mathrm P}, {\mathrm D})
\end{alignat*}
where
\begin{alignat*}{2}
s\left({\mathrm P}, {\mathrm M} \right) & = & \hspace{0.25em} \cos^{-1} \! \left\{
\sin {\phi_{}}_{\mathrm P} \sin {\phi_{}}_{\mathrm M}
+ \cos {\phi_{}}_{\mathrm P} \cos {\phi_{}}_{\mathrm M}
\cos ({\lambda_{}}_{\mathrm M} - {\lambda_{}}_{\mathrm P})
\right\}
\end{alignat*}
and $M$ corresponds to $B$, $C$ or $D$.
A more stable form of the great-circle distance formula for small distances ($x$ near 1)
involves the arcsine function (\eg\ see p.~101 of \citet{daley.barker_bk01}:
\begin{alignat*}{2}
s\left( {\mathrm P}, {\mathrm M} \right) = \sin^{-1} \! \left\{ \sqrt{ 1 - x^2 } \right\}
\end{alignat*}
where
\begin{alignat*}{2}
x = {a_{}}_{\mathrm M} {a_{}}_{\mathrm P} + {b_{}}_{\mathrm M} {b_{}}_{\mathrm P} + {c_{}}_{\mathrm M} {c_{}}_{\mathrm P}
\end{alignat*}
and
\begin{alignat*}{3}
& {a_{}}_{\mathrm M} & = && \quad \sin {\phi_{}}_{\mathrm M}, \\
& {a_{}}_{\mathrm P} & = && \quad \sin {\phi_{}}_{\mathrm P}, \\
& {b_{}}_{\mathrm M} & = && \quad \cos {\phi_{}}_{\mathrm M} \cos {\phi_{}}_{\mathrm M}, \\
& {b_{}}_{\mathrm P} & = && \quad \cos {\phi_{}}_{\mathrm P} \cos {\phi_{}}_{\mathrm P}, \\
& {c_{}}_{\mathrm M} & = && \quad \cos {\phi_{}}_{\mathrm M} \sin {\phi_{}}_{\mathrm M}, \\
& {c_{}}_{\mathrm P} & = && \quad \cos {\phi_{}}_{\mathrm P} \sin {\phi_{}}_{\mathrm P}.
\end{alignat*}
\item {\bfseries Great-Circle distance-weighted interpolation with small angle approximation.}
Similar to the previous interpolation but with the distance $s$ computed as
\begin{alignat*}{2}
s\left( {\mathrm P}, {\mathrm M} \right)
& = & \sqrt{ \left( {\phi_{}}_{\mathrm M} - {\phi_{}}_{\mathrm P} \right)^{2}
+ \left( {\lambda_{}}_{\mathrm M} - {\lambda_{}}_{\mathrm P} \right)^{2}
\cos^{2} {\phi_{}}_{\mathrm M} }
\end{alignat*}
where $M$ corresponds to $A$, $B$, $C$ or $D$.
\item {\bfseries Bilinear interpolation for a regular spaced grid.}
The interpolation is split into two 1D interpolations in the longitude and latitude directions, respectively.
\item {\bfseries Bilinear remapping interpolation for a general grid.}
An iterative scheme that involves first mapping a quadrilateral cell into
a cell with coordinates (0,0), (1,0), (0,1) and (1,1).
This method is based on the \href{https://github.com/SCRIP-Project/SCRIP}{SCRIP interpolation package}.
\end{enumerate}
%% =================================================================================================
\subsubsection{Horizontal averaging}
For each surface observation type:
\begin{itemize}
\item The standard grid-searching code is used to find the nearest model grid point to the observation location
(see next subsection).
\item The maximum number of grid points required for that observation in each local grid domain is calculated. Some of these points may later turn out to have zero weight depending on the shape of the footprint.
\item The longitudes and latitudes of the grid points surrounding the nearest model grid box are extracted using
existing MPI routines.
\item The weights for each grid point associated with each observation are calculated,
either for radial or rectangular footprints.
For grid points completely within the footprint, the weight is one;
for grid points completely outside the footprint, the weight is zero.
For grid points which are partly within the footprint the ratio between the area of the footprint within
the grid box and the total area of the grid box is used as the weight.
\item The weighted average of the model grid points associated with each observation is calculated,
and this is then given as the model counterpart of the observation.
\end{itemize}
Examples of the weights calculated for an observation with rectangular and radial footprints are shown in
\autoref{fig:OBS_avgrec} and~\autoref{fig:OBS_avgrad}.
\begin{figure}
\centering
\includegraphics[width=0.66\textwidth]{OBS_avg_rec}
\caption[Observational weights with a rectangular footprint]{
Weights associated with each model grid box (blue lines and numbers)
for an observation at -170.5\deg{E}, 56.0\deg{N} with a rectangular footprint of 1\deg\ x 1\deg.}
\label{fig:OBS_avgrec}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=0.66\textwidth]{OBS_avg_rad}
\caption[Observational weights with a radial footprint]{
Weights associated with each model grid box (blue lines and numbers)
for an observation at -170.5\deg{E}, 56.0\deg{N} with a radial footprint with diameter 1\deg.}
\label{fig:OBS_avgrad}
\end{figure}
%% =================================================================================================
\subsection{Grid search}
For many grids used by the \NEMO\ model, such as the ORCA family, the horizontal grid coordinates $i$ and $j$ are not simple functions of latitude and longitude.
Therefore, it is not always straightforward to determine the grid points surrounding any given observational position.
Before the interpolation can be performed, a search algorithm is then required to determine the corner points of
the quadrilateral cell in which the observation is located.
This is the most difficult and time consuming part of the 2D interpolation procedure.
A robust test for determining if an observation falls within a given quadrilateral cell is as follows.
Let ${\mathrm P}({\lambda_{}}_{\mathrm P} ,{\phi_{}}_{\mathrm P} )$ denote the observation point,
and let ${\mathrm A}({\lambda_{}}_{\mathrm A} ,{\phi_{}}_{\mathrm A} )$, ${\mathrm B}({\lambda_{}}_{\mathrm B} ,{\phi_{}}_{\mathrm B} )$,
${\mathrm C}({\lambda_{}}_{\mathrm C} ,{\phi_{}}_{\mathrm C} )$ and ${\mathrm D}({\lambda_{}}_{\mathrm D} ,{\phi_{}}_{\mathrm D} )$
denote the bottom left, bottom right, top left and top right corner points of the cell, respectively.
To determine if P is inside the cell, we verify that the cross-products
\begin{align*}
\begin{array}{lllll}
{{\mathbf r}_{}}_{\mathrm PA} \times {{\mathbf r}_{}}_{\mathrm PC}
& = & [({\lambda_{}}_{\mathrm A}\; -\; {\lambda_{}}_{\mathrm P} )
({\phi_{}}_{\mathrm C} \; -\; {\phi_{}}_{\mathrm P} )
- ({\lambda_{}}_{\mathrm C}\; -\; {\lambda_{}}_{\mathrm P} )
({\phi_{}}_{\mathrm A} \; -\; {\phi_{}}_{\mathrm P} )] \; \widehat{\mathbf k} \\
{{\mathbf r}_{}}_{\mathrm PB} \times {{\mathbf r}_{}}_{\mathrm PA}
& = & [({\lambda_{}}_{\mathrm B}\; -\; {\lambda_{}}_{\mathrm P} )
({\phi_{}}_{\mathrm A} \; -\; {\phi_{}}_{\mathrm P} )
- ({\lambda_{}}_{\mathrm A}\; -\; {\lambda_{}}_{\mathrm P} )
({\phi_{}}_{\mathrm B} \; -\; {\phi_{}}_{\mathrm P} )] \; \widehat{\mathbf k} \\
{{\mathbf r}_{}}_{\mathrm PC} \times {{\mathbf r}_{}}_{\mathrm PD}
& = & [({\lambda_{}}_{\mathrm C}\; -\; {\lambda_{}}_{\mathrm P} )
({\phi_{}}_{\mathrm D} \; -\; {\phi_{}}_{\mathrm P} )
- ({\lambda_{}}_{\mathrm D}\; -\; {\lambda_{}}_{\mathrm P} )
({\phi_{}}_{\mathrm C} \; -\; {\phi_{}}_{\mathrm P} )] \; \widehat{\mathbf k} \\
{{\mathbf r}_{}}_{\mathrm PD} \times {{\mathbf r}_{}}_{\mathrm PB}
& = & [({\lambda_{}}_{\mathrm D}\; -\; {\lambda_{}}_{\mathrm P} )
({\phi_{}}_{\mathrm B} \; -\; {\phi_{}}_{\mathrm P} )
- ({\lambda_{}}_{\mathrm B}\; -\; {\lambda_{}}_{\mathrm P} )
({\phi_{}}_{\mathrm D} \; - \; {\phi_{}}_{\mathrm P} )] \; \widehat{\mathbf k} \\
\end{array}
% \label{eq:OBS_cross}
\end{align*}
point in the opposite direction to the unit normal $\widehat{\mathbf k}$
(\ie\ that the coefficients of $\widehat{\mathbf k}$ are negative),
where ${{\mathbf r}_{}}_{\mathrm PA}$, ${{\mathbf r}_{}}_{\mathrm PB}$, etc. correspond to
the vectors between points P and A, P and B, etc..
The method used is similar to the method used in the \href{https://github.com/SCRIP-Project/SCRIP}{SCRIP interpolation package}.
In order to speed up the grid search, there is the possibility to construct a lookup table for a user specified resolution.
This lookup table contains the lower and upper bounds on the $i$ and $j$ indices to
be searched for on a regular grid.
For each observation position, the closest point on the regular grid of this position is computed and
the $i$ and $j$ ranges of this point searched to determine the precise four points surrounding the observation.
%% =================================================================================================
\subsection{Parallel aspects of horizontal interpolation}
\label{subsec:OBS_parallel}
For horizontal interpolation, there is the basic problem that
the observations are unevenly distributed on the globe.
In \NEMO\ the model grid is divided into subgrids (or domains) where
each subgrid is executed on a single processing element with explicit message passing for
exchange of information along the domain boundaries when running on a massively parallel processor (MPP) system.
For observations there is no natural distribution since the observations are not equally distributed on the globe.
Two options have been made available:
1) geographical distribution;
and 2) round-robin.
%% =================================================================================================
\subsubsection{Geographical distribution of observations among processors}
\begin{figure}
\centering
\includegraphics[width=0.66\textwidth]{OBS_obsdist_local}
\caption[Observations with the geographical distribution]{
Example of the distribution of observations with
the geographical distribution of observational data}
\label{fig:OBS_local}
\end{figure}
This is the simplest option in which the observations are distributed according to
the domain of the grid-point parallelization.
\autoref{fig:OBS_local} shows an example of the distribution of the {\em in situ} data on processors with
a different colour for each observation on a given processor for a 4 $\times$ 2 decomposition with ORCA2.
The grid-point domain decomposition is clearly visible on the plot.
The advantage of this approach is that all information needed for horizontal interpolation is available without
any MPP communication.
This is under the assumption that we are dealing with point observations and only using a $2 \times 2$ grid-point stencil for
the interpolation (\eg\ bilinear interpolation).
For higher order interpolation schemes this is no longer valid.
A disadvantage with the above scheme is that the number of observations on each processor can be very different.
If the cost of the actual interpolation is expensive relative to the communication of data needed for interpolation,
this could lead to load imbalance.
%% =================================================================================================
\subsubsection{Round-robin distribution of observations among processors}
\begin{figure}
\centering
\includegraphics[width=0.66\textwidth]{OBS_obsdist_global}
\caption[Observations with the round-robin distribution]{
Example of the distribution of observations with
the round-robin distribution of observational data.}
\label{fig:OBS_global}
\end{figure}
An alternative approach is to distribute the observations equally among processors and
use message passing in order to retrieve the stencil for interpolation.
The simplest distribution of the observations is to distribute them using a round-robin scheme.
\autoref{fig:OBS_global} shows the distribution of the {\em in situ} data on processors for
the round-robin distribution of observations with a different colour for each observation on a given processor for
a 4 $\times$ 2 decomposition with ORCA2 for the same input data as in \autoref{fig:OBS_local}.
The observations are now clearly randomly distributed on the globe.
In order to be able to perform horizontal interpolation in this case,
a subroutine has been developed that retrieves any grid points in the global space.
%% =================================================================================================
\subsection{Vertical interpolation operator}
Vertical interpolation is achieved using either a cubic spline or linear interpolation.
For the cubic spline, the top and bottom boundary conditions for the second derivative of
the interpolating polynomial in the spline are set to zero.
At the bottom boundary, this is done using the land-ocean mask.
For profile observation types we do both vertical and horizontal interpolation. \NEMO\ has a generalised vertical coordinate system this means the vertical level depths can vary with location. Therefore, it is necessary first to perform vertical interpolation of the model value to the observation depths for each of the four surrounding grid points. After this the model values, at these points, at the observation depth, are horizontally interpolated to the observation location.