diff --git a/ioda/create_ioda_obs.m b/ioda/create_ioda_obs.m
index 32b35c1..afea2fa 100644
--- a/ioda/create_ioda_obs.m
+++ b/ioda/create_ioda_obs.m
@@ -13,13 +13,21 @@ function create_ioda_obs(S, file)
 %   group: MetaData {
 %     variables:
 %       int64 dateTime(Location) ;
-%       string date_time(Location) ;
+%       int64 dateTimeAverageBegin(timeWindow) ;
+%       int64 dateTimeAverageEnd(timeWindow) ;
 %       float depth(Location) ;
 %       float latitude(Location) ;
 %       float longitude(Location) ;
-%       int provenance(Location) ;  
+%       int provenance(Location) ;
 %       int sequenceNumber(Location) ;
+%       float spatialAverage ;
+%       int stateID(nvars);
+%       int surveyIndex(survey)  ;
+%       int64 surveyTime(survey)  ;
 %       string variables_name(nvars) ;
+%       float x_grid(Location) ;
+%       float y_grid(Location) ;
+%       float z_grid(Location) ;
 %   }
 %
 %   group: ObsError {
@@ -36,21 +44,36 @@ function create_ioda_obs(S, file)
 %     variables:
 %       float myObsVariableName(Location) ;
 %   }
-%  
+%
 % On Input:
 %
 %    S        Observations file creation parameters (struct):
 %
 %               S.ncfile            NetCDF file name (string)
+%               S.ncvname(:)        IODA NetCDF4 variable name
 %               S.nlocs             number of observations
+%               S.nsurvey           number of survey times
 %               S.nvars             number of variables
-%               S.ncvname(:)        IODA NetCDF4 variable name
-%               S.variable_names(:) IODA variable standard name
-%               S.units(:)          variables units
+%               S.timeWindow        mumber of time-aeraged windows
+%               S.dateTime          observation time
 %               S.depth             depth of observation
+%               S.latitude          observation latitude
+%               S.longitude         observation latitude
 %               S.provenance        observation provenance
-%               S.TimeIODA          IODA reference time YYYYMMDDHH
+%               S.sequenceNumber    observation sequence number
+%               S.areaAvgRadius     area-averaged radius
+%               S.timeAvgBegin      time-averaged window start
+%               S.timeAvgEnd        time-averaged window end
+%               S.stateId           state variable identifier
+%               S.surveyIndex       indices as it appear in dateTime
+%               S.surveyTime        observation survey times
+%               S.units(:)          variables units
+%               S.variable_names(:) IODA variable standard name
+%               S.x_grid            fractional x-grid location
+%               S.y_grid            fractional y-grid location
+%               S.z_grid            fractional z-grid location
 %               S.DateIODA          IODA reference date string
+%               S.TimeIODA          IODA reference time YYYYMMDDHH
 %
 %    file     Output observation file name (string, OPTIONAL). If provided,
 %               it creates this file instead the one in S.ncfile.
@@ -63,14 +86,10 @@ function create_ioda_obs(S, file)
 %
 % group: MetaData {
 %
-%   string date_time(Location) ;
-%          date_time:long_name = "ISO 8601 UTC date and time string" ;
 %   string variables_name(nvars) ;
 %          variables_name:long_name = "observation UFO/IODA standard name" ;
 %
-% data:  
-%
-%   date_time = "2004-01-03T14:21:00Z", "2004-01-03T14:21:00Z", ... ;
+% data:
 %
 %   variables_name = "absolute_dynamic_topography" ;
 % }
@@ -105,9 +124,10 @@ function create_ioda_obs(S, file)
 
 % Set indices for local dimension structure vector element.
 
-nlocs = 1;           % "Location" dimension index
-nvars = 2;           % "nvars"    dimension index
-nstr  = 5;           % "nstring"  dimension index
+nlocs   = 1;         % "Location"   dimension index
+nvars   = 2;         % "nvars"      dimension index
+nsurv   = 3;         % "survey"     dimension index
+nwindow = 4;         % "timeWindow" dimension index
 
 % Set indices for NetCDF/HDF5 Group local structure vector element.
 
@@ -130,21 +150,70 @@ function create_ioda_obs(S, file)
   if (~isempty(S.depth))
     do_depth = true;           % sub-surface observations
   end
-end  
-  
+end
+
+do_fractional = false;         % do not include fractional locations
+if (isfield(S, 'x_grid'))
+  if (~isempty(S.x_grid))
+    do_fractional = true;      % include fractional locations
+  end
+end
+
 do_provenance = false;         % do not include observation origin
 if (isfield(S, 'provenance'))
   if (~isempty(S.provenance))
     do_provenance = true;      % include observation origin
-  end 
+  end
+end
+
+do_stateID = false;            % do not include state variable ID
+surface_obs = false;
+if (isfield(S, 'stateID'))
+  if (~isempty(S.stateID))
+    do_stateID = true;         % include state variable ID
+  end
+  for i = 1:S.nvars
+    Vname = char(S.ncvname(i));
+    if (contains(Vname, 'Surface'))
+      surface_obs = true;
+    end
+  end
+end
+
+do_areaAvg = false;
+if (isfield(S, 'areaAvgRadius'))
+  if (~isnan(S.areaAvgRadius))
+    do_areaAvg = true;
+  end
+end
+
+do_timeAvgBegin = false;
+if (isfield(S, 'timeAvgBegin'))
+  if (~isnan(S.timeAvgBegin))
+    do_timeAvgBegin = true;
+  end
 end
-  
+
+do_timeAvgEnd = false;
+if (isfield(S, 'timeAvgEnd'))
+  if (~isnan(S.timeAvgEnd))
+    do_timeAvgEnd = true;
+  end
+end
+
 %--------------------------------------------------------------------------
 % Set NetCDF4 file dimensions.
 %--------------------------------------------------------------------------
 
-D(nlocs).name = 'Location';     D(nlocs).size = S.nlocs;
-D(nvars).name = 'nvars';        D(nvars).size = S.nvars;
+D(nlocs).name = 'Location';       D(nlocs).size = S.nlocs;
+D(nvars).name = 'nvars';          D(nvars).size = S.nvars;
+D(nsurv).name = 'survey';         D(nsurv).size = S.nsurvey;
+
+if (isfield(S, 'nwindow'))
+  if (~isnan(S.nwindow))
+    D(nwindow).name = 'timeWindow'; D(nwindow).size = S.nwindow;
+  end
+end
 
 %--------------------------------------------------------------------------
 % Create NetCDF4 observation file.
@@ -153,7 +222,7 @@ function create_ioda_obs(S, file)
 disp(' ');
 disp(['*** Creating observations file:  ', ncfile]);
 
-mode = netcdf.getConstant('NETCDF4'); 
+mode = netcdf.getConstant('NETCDF4');
 mode = bitor(mode, netcdf.getConstant('NC_CLOBBER'));
 
 ncid = netcdf.create(ncfile, mode);
@@ -198,11 +267,19 @@ function create_ioda_obs(S, file)
     case 'Location'
       D(i).vid = netcdf.defVar(ncid, Vname, nc_int, D(i).did);
       netcdf.putAtt(ncid, D(i).vid, 'suggested_chunck_dim',             ...
-                    int32(D(i).size));
+                    int32(512));
     case 'nvars'
       D(i).vid = netcdf.defVar(ncid, Vname, nc_int, D(i).did);
       netcdf.putAtt(ncid, D(i).vid, 'suggested_chunck_dim',             ...
                     int32(100));
+    case 'survey'
+      D(i).vid = netcdf.defVar(ncid, Vname, nc_int, D(i).did);
+      netcdf.putAtt(ncid, D(i).vid, 'suggested_chunck_dim',             ...
+                    int32(100));
+    case 'timeWindow'
+      D(i).vid = netcdf.defVar(ncid, Vname, nc_int, D(i).did);
+      netcdf.putAtt(ncid, D(i).vid, 'suggested_chunck_dim',             ...
+                    int32(100));
   end
 end
 
@@ -221,14 +298,33 @@ function create_ioda_obs(S, file)
 %--------------------------------------------------------------------------
 
 MetaVars = {'dateTime', 'latitude', 'longitude', 'sequenceNumber',      ...
-            'variables_name'};
-MetaVars = [MetaVars, 'date_time'];
+            'surveyTime' 'surveyIndex' 'variables_name'};
 if (do_depth)
   MetaVars = [MetaVars, 'depth'];
 end
+if (do_fractional)
+  MetaVars = [MetaVars, 'x_grid', 'y_grid'];
+  if (do_depth)
+    MetaVars = [MetaVars, 'z_grid'];
+  end
+end
 if (do_provenance)
   MetaVars = [MetaVars, 'provenance'];
 end
+if (do_stateID)
+  MetaVars = [MetaVars, 'stateID'];
+end
+if (do_areaAvg)
+  MetaVars = [MetaVars, 'spatialAverage'];
+end
+if (do_timeAvgBegin)
+  MetaVars = [MetaVars, 'dateTimeAverageBegin'];
+end
+if (do_timeAvgEnd)
+  MetaVars = [MetaVars, 'dateTimeAverageEnd'];
+end
+
+
 G(Meta).vars = sort(MetaVars);
 
 for i = 1:length(G(Meta).vars)
@@ -242,14 +338,22 @@ function create_ioda_obs(S, file)
       epoch  = datenum(num2str(S.datetime_ref),'yyyymmddHH');
       string = ['seconds since ' datestr(epoch, 'yyyy-mm-ddTHH:MM:SSZ')];
       netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'units', string);
-    case 'date_time'
-      G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_string,     ...
-                                     D(nlocs).did);
-      string = 'ISO 8601 UTC date and time string';
+    case 'dateTimeAverageBegin'
+      G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_int64,      ...
+                                     D(nwindow).did);
+      string = 'start of time averaging filter';
       netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'long_name', string);
-      string = blanks(0);
-      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), '_FillValue', string,  ...
-                    nc_string);
+      epoch  = datenum(num2str(S.datetime_ref),'yyyymmddHH');
+      string = ['seconds since ' datestr(epoch, 'yyyy-mm-ddTHH:MM:SSZ')];
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'units', string);
+    case 'dateTimeAverageEnd'
+      G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_int64,      ...
+                                     D(nwindow).did);
+      string = 'end of time averaging filter';
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'long_name', string);
+      epoch  = datenum(num2str(S.datetime_ref),'yyyymmddHH');
+      string = ['seconds since ' datestr(epoch, 'yyyy-mm-ddTHH:MM:SSZ')];
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'units', string);
     case 'depth'
       G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_real,       ...
                                      D(nlocs).did);
@@ -274,16 +378,53 @@ function create_ioda_obs(S, file)
                                      D(nlocs).did);
       string = 'observation origin identifier';
       netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'long_name', string);
+      if (any(S.stateID == 1))
+        string = 'repetitive observation at different times and error';
+        netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'negative', string);
+      end
       if (isfield(S, 'flag_values'))
         netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'flag_values',       ...
-	              S.flag_values);
+                      S.flag_values);
         netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'flag_meanings',     ...
-	              S.flag_meanings);
-      end	
+                      S.flag_meanings);
+      end
+    case 'stateID'
+      G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_int,        ...
+                                     D(nvars).did);
+      string = 'state variable index';
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'long_name', string);
+      if (surface_obs)
+        netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'surface_level',     ...
+                    int32(S.N));
+      end
+%     values = 1:7;
+%     netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'flag_values',         ...
+%                   int32(values));
+%     values = 'zeta ubar vbar u v temperature salinity';
+%     netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'flag_meanings',       ...
+%                   values);
     case 'sequenceNumber'
       G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_int,        ...
                                      D(nlocs).did);
-      string = 'observations sequence number';
+      string = 'observation sequence number';
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'long_name', string);
+    case 'spatialAverage'
+      G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_real, []);
+      string = 'half-length of spatial averaging filter';
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'long_name', string);
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'units', 'meter');
+    case 'surveyTime'
+      G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_int64,      ...
+                                     D(nsurv).did);
+      string = 'observation survey time';
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'long_name', string);
+      epoch  = datenum(num2str(S.datetime_ref),'yyyymmddHH');
+      string = ['seconds since ' datestr(epoch, 'yyyy-mm-ddTHH:MM:SSZ')];
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'units', string);
+    case 'surveyIndex'
+      G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_int,        ...
+                                     D(nsurv).did);
+      string = 'observation survey time indices as they appear in dateTime';
       netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'long_name', string);
     case 'variables_name'
       G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_string,     ...
@@ -293,6 +434,21 @@ function create_ioda_obs(S, file)
       string = blanks(0);
       netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), '_FillValue', string,  ...
                     nc_string);
+    case 'x_grid'
+      G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_real,       ...
+                                     D(nlocs).did);
+      string = 'observation fractional x-grid location';
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'long_name', string);
+    case 'y_grid'
+      G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_real,       ...
+                                     D(nlocs).did);
+      string = 'observation fractional y-grid location';
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'long_name', string);
+    case 'z_grid'
+      G(Meta).vid(i) = netcdf.defVar(G(Meta).gid, Vname, nc_real,       ...
+                                     D(nlocs).did);
+      string = 'observation fractional z-grid location';
+      netcdf.putAtt(G(Meta).gid, G(Meta).vid(i), 'long_name', string);
   end
 end
 
diff --git a/ioda/ioda_read.m b/ioda/ioda_read.m
index 2bb2dee..d2a21af 100644
--- a/ioda/ioda_read.m
+++ b/ioda/ioda_read.m
@@ -22,7 +22,7 @@
 %              S.nvars            number of variables
 %              S.epoch            IODA time reference YYYYMMDDHH
 %              S.datenum          epoch date number
-%              S.dateTimeRef      IODA time reference string     
+%              S.dateTimeRef      IODA time reference string
 %              S.variable_names   UFO/IODA standard name
 %              S.dateTime         seconds since yyyy-mm-ddTHH:MM:SSZ
 %              S.date_time        date and time ISO 8601 UTC string
@@ -42,7 +42,7 @@
 %    Licensed under a MIT/X style license                                 %
 %    See License_ROMS.md                            Hernan G. Arango      %
 %=========================================================================%
-  
+
 % Initialize.
 
 S = struct('ncfile'           , [],                                     ...
@@ -106,7 +106,9 @@
 
 S.dateTime = double(ncread(ncfile, '/MetaData/dateTime'));
 
-S.date_time = ncread(ncfile, '/MetaData/date_time');
+if (any(strcmp({G.Variables.Name}, 'date_time')))
+  S.date_time = ncread(ncfile, '/MetaData/date_time');
+end
 
 % Set IODA NetCDF variables.
 
@@ -196,6 +198,30 @@
   end
 end
 
+if (any(strcmp({I.Groups.Name}, 'hofx0_1')))
+  for i = 1:S.nvars
+    Vname = strcat('/hofx0_1/', S.iodaVarName{i});
+    field = double(ncread(ncfile, Vname));
+    S.hofx0_1{i} = field;
+  end
+end
+
+if (any(strcmp({I.Groups.Name}, 'hofx0_2')))
+  for i = 1:S.nvars
+    Vname = strcat('/hofx0_2/', S.iodaVarName{i});
+    field = double(ncread(ncfile, Vname));
+    S.hofx0_2{i} = field;
+  end
+end
+
+if (any(strcmp({I.Groups.Name}, 'hofx0_3')))
+  for i = 1:S.nvars
+    Vname = strcat('/hofx0_3/', S.iodaVarName{i});
+    field = double(ncread(ncfile, Vname));
+    S.hofx0_3{i} = field;
+  end
+end
+
 % Read in 'hofx1' Group: Final H(x).
 
 if (any(strcmp({I.Groups.Name}, 'hofx1')))
@@ -206,6 +232,30 @@
   end
 end
 
+if (any(strcmp({I.Groups.Name}, 'hofx1_1')))
+  for i = 1:S.nvars
+    Vname = strcat('/hofx0_1/', S.iodaVarName{i});
+    field = double(ncread(ncfile, Vname));
+    S.hofx1_1{i} = field;
+  end
+end
+
+if (any(strcmp({I.Groups.Name}, 'hofx1_2')))
+  for i = 1:S.nvars
+    Vname = strcat('/hofx0_2/', S.iodaVarName{i});
+    field = double(ncread(ncfile, Vname));
+    S.hofx1_2{i} = field;
+  end
+end
+
+if (any(strcmp({I.Groups.Name}, 'hofx1_3')))
+  for i = 1:S.nvars
+    Vname = strcat('/hofx0_3/', S.iodaVarName{i});
+    field = double(ncread(ncfile, Vname));
+    S.hofx1_3{i} = field;
+  end
+end
+
 % Read in 'oman' Group: Observation minus analysis.
 
 if (any(strcmp({I.Groups.Name}, 'oman')))
diff --git a/ioda/plot_ioda.m b/ioda/plot_ioda.m
index f1d2a4c..de9f148 100644
--- a/ioda/plot_ioda.m
+++ b/ioda/plot_ioda.m
@@ -1,8 +1,8 @@
 function S = plot_ioda (ncname, varargin)
-  
+
 % PLOT_IODA:  Plots field from an IODA/UFO files.
 %
-% plot_ioda(ncname, ptype, wrtPNG)
+% plot_ioda(ncname, ptype, PNGprefix)
 %
 % On Input:
 %
@@ -15,9 +15,10 @@
 %                    'ObsValue'    ~> observation value
 %                    'oman'        ~> observation minus analysis
 %                    'ombg'        ~> observation minus background
+%                    []            ~> plots selected groups
 %
-%    wrtPNG        Switch to write out PNG file (true or false)
-%                    (Optional, default: false)
+%    PNGprefix     PNG filename prefix (string; OPTIONAL)
+%                    (default: PNG files are not generated)
 %
 % On Output:
 %
@@ -41,29 +42,29 @@
 
 switch numel(varargin)
   case 0
-    group   = "ObsValue";                        % default
-    got_grp = false;
-    wrtPNG  = false;
+    group     = "ObsValue";                        % default
+    got_grp   = false;
+    PNGprefix = [];
   case 1
-    group   = varargin{1}
+    group     = varargin{1};
     if (isempty(group))
       group   = "ObsValue";
       got_grp = false;
     else
-      got_grp = true; 
+      got_grp = true;
     end
-    wrtPNG  = false;
+    PNGprefix = [];
   case 2
-    group   = varargin{1};
+    group     = varargin{1};
     if (isempty(group))
       group   = "ObsValue";
       got_grp = false;
     else
-      got_grp = true; 
+      got_grp = true;
     end
-    wrtPNG  = varargin{2};
+    PNGprefix = varargin{2};
 end
-  
+
 % Check if requested group is available
 
 if (got_grp)
@@ -96,7 +97,7 @@
         for j=1:addit
           s = [ s ' '];
         end
-      end 
+      end
       if (i < Ngroups - 1)
         stri = int2str(i+3);
         if (length(stri) == 1)
@@ -104,13 +105,13 @@
         end
         grpnam = I.Groups(i+2).Name;
         s = [ s '  ' stri ') ' grpnam ];
-      end 
+      end
       disp(s);
     end
     disp(blanks(1));
-    error(['PLOT_IODA: cannot find NetCDF4 Group: ', group]);  
-  
-  end  
+    error(['PLOT_IODA: cannot find NetCDF4 Group: ', group]);
+
+  end
 end
 
 % Read input IODA/UFO NetCDF4 file.
@@ -138,7 +139,7 @@
   shortname = {'vsur'};
 elseif (strfind(ncname, 'uv_sur'))
   shortname = {'u_sur', 'v_sur'};
-end 
+end
 S.shortname = shortname;
 
 % Set Date for time axis labels.
@@ -172,7 +173,7 @@
   Vname = S.iodaVarName{n};
 
 %--------------------------------------------------------------------------
-% Plot "ObsValue" Group.  
+% Plot "ObsValue" Group.
 %--------------------------------------------------------------------------
 
   if ( (~got_grp && any(strcmp({I.Groups.Name}, 'ObsValue'))) ||        ...
@@ -185,9 +186,9 @@
       ylabel('depth');
       title(['Group:  ObsValue,   File:  ', untexlabel(MyFile)]);
       grid on;
-      if (wrtPNG)
-        png_file = strcat(shortname{n}, '_obsZ.png');
-        print(png_file, '-dpng', '-r300')
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_obsZ.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
       end
     else
       figure;
@@ -197,16 +198,34 @@
       zlabel(Vname);
       grid on;
       title(['Group:  ObsValue,   File:  ', untexlabel(MyFile)]);
-      if (wrtPNG)
-        png_file = strcat(shortname{n}, '_obsXY.png');
-        print(png_file, '-dpng', '-r300')
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_obsXY.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
+      end
+
+      if (strcmp(Vname, 'absoluteDynamicTopography'))
+        ind = find(S.provenance > 0);
+        if ~isempty(ind)                % removes repetitive observations
+          figure;                       % with negative provenance
+          Value = S.ObsValue{n};
+          h2 = plot3(S.longitude(ind), S.latitude(ind), Value(ind), 'r+');
+          xlabel('Longitude');
+          ylabel('latitude');
+          zlabel(Vname);
+          grid on;
+          title(['Group:  UniqueObs,   File:  ', untexlabel(MyFile)]);
+          if (~isempty(PNGprefix))
+            png_file = strcat(PNGprefix, shortname{n}, '_noreps_obsXY.png');
+            exportgraphics(gcf, png_file, 'resolution', 300);
+          end
+        end
       end
-    end  
-  
+    end
+
   end
-    
+
 %--------------------------------------------------------------------------
-% Plot "hofx" Group.  
+% Plot "hofx" Group.
 %--------------------------------------------------------------------------
 
   if ( (~got_grp && any(strcmp({I.Groups.Name}, 'hofx'))) ||            ...
@@ -220,11 +239,11 @@
     ylabel(Vname);
     grid on;
     title(['Group:  hofx,   File:  ', untexlabel(MyFile)]);
-    if (wrtPNG)
-      png_file = strcat(shortname{n}, '_enum.png');
-      print(png_file, '-dpng', '-r300')
+    if (~isempty(PNGprefix))
+      png_file = strcat(PNGprefix, shortname{n}, '_enum.png');
+      exportgraphics(gcf, png_file, 'resolution', 300);
     end
-       
+
     figure;
     h2 = plot(t, S.ObsValue{n}, 'r+',                                   ...
               t, S.hofx{n}, 'bo');
@@ -234,11 +253,11 @@
     ylabel(Vname);
     grid on;
     title(['Group:  hofx,   File:  ', untexlabel(MyFile)]);
-    if (wrtPNG)
-      png_file = strcat(shortname{n}, '_date.png');
-      print(png_file, '-dpng', '-r300')
+    if (~isempty(PNGprefix))
+      png_file = strcat(PNGprefix, shortname{n}, '_date.png');
+      exportgraphics(gcf, png_file, 'resolution', 300);
     end
-       
+
     figure;
     h3 = scatter(Nlocs, S.ObsValue{n} - S.hofx{n}, 'bo',                ...
                  'MarkerFaceColor', 'b');
@@ -248,11 +267,11 @@
     xlabel('Observation Number')
     ylabel([Vname, ' Innovation Vector']);
     title(['d = OBS - HofX for File:  ', untexlabel(MyFile)]);
-    if (wrtPNG)
-      png_file = strcat(shortname{n}, '_innov.png');
-      print(png_file, '-dpng', '-r300')
+    if (~isempty(PNGprefix))
+      png_file = strcat(PNGprefix, shortname{n}, '_innov.png');
+      exportgraphics(gcf, png_file, 'resolution', 300);
     end
-       
+
     if (any(strcmp({M.Variables.Name}, 'depth')))
       figure;
       h4 = plot(S.ObsValue{n}, S.depth, 'r+',                           ...
@@ -262,9 +281,9 @@
       ylabel('depth');
       grid on;
       title(['HofX File:  ', untexlabel(MyFile)]);
-      if (wrtPNG)
-        png_file = strcat(shortname{n}, '_depth.png');
-        print(png_file, '-dpng', '-r300')
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_depth.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
       end
     else
       figure;
@@ -273,21 +292,21 @@
       legend('OBS', 'HofX', 'Location', 'best');
       xlabel('Longitude');
       ylabel('latitude');
-      zlabel(Vname); 
+      zlabel(Vname);
       grid on;
       title(['Group:  hofx,   File:  ', untexlabel(MyFile)]);
-      if (wrtPNG)
-        png_file = strcat(shortname{n}, '_lonlat.png');
-        print(png_file, '-dpng', '-r300')
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_lonlat.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
       end
     end
-  
+
   end
 
 %--------------------------------------------------------------------------
 %  Plot "hofx0" Group: initial HofX.
 %--------------------------------------------------------------------------
-  
+
   if ( (~got_grp && any(strcmp({I.Groups.Name}, 'hofx0'))) ||           ...
        (strcmp(group, 'hofx0')) )
 
@@ -299,46 +318,95 @@
     ylabel(Vname);
     grid on;
     title(['Group:  hofx0,   File:  ', untexlabel(MyFile)]);
-    if (wrtPNG)
-      png_file = strcat(shortname{n}, '_hofx0.png');
-      print(png_file, '-dpng', '-r300')
+    if (~isempty(PNGprefix))
+      png_file = strcat(PNGprefix, shortname{n}, '_hofx0.png');
+      exportgraphics(gcf, png_file, 'resolution', 300);
     end
-  
+
     if (any(strcmp({M.Variables.Name}, 'depth')))
       figure;
       h4 = plot(S.ObsValue{n}, S.depth, 'r+',                           ...
-                S.hofx1{n}, S.depth, 'bo');
+                S.hofx0{n}, S.depth, 'bo');
       legend('OBS', 'Initial HofX', 'Location', 'best');
       xlabel(Vname);
       ylabel('depth');
       grid on;
       title(['Group:  hofx0,   File:  ', untexlabel(MyFile)]);
-      if (wrtPNG)
-        png_file = strcat(shortname{n}, '_hofx0_depth.png');
-        print(png_file, '-dpng', '-r300')
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_hofx0_depth.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
       end
     else
       figure;
       h4 = plot3(S.longitude, S.latitude, S.ObsValue{n}, 'r+',         ...
-                 S.longitude, S.latitude, S.hofx1{n}, 'bo');
+                 S.longitude, S.latitude, S.hofx0{n}, 'bo');
       legend('OBS', 'Initial HofX', 'Location', 'best');
       xlabel('Longitude');
       ylabel('latitude');
-      zlabel(Vname); 
+      zlabel(Vname);
       grid on;
       title(['Group:  hofx0,   File:  ', untexlabel(MyFile)]);
-      if (wrtPNG)
-        png_file = strcat(shortname{n}, '_hofx0_lonlat.png');
-        print(png_file, '-dpng', '-r300')
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_hofx0_lonlat.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
+      end
+    end
+
+  end
+
+% Ensemble HofX.
+
+  if ( (~got_grp && any(strcmp({I.Groups.Name}, 'hofx0_1'))) ||         ...
+       (strcmp(group, 'hofx0_1')) )
+
+    figure;
+    h1 = plot(Nlocs, S.ObsValue{n}, 'r+',                               ...
+              Nlocs, S.hofx0_1{n}, 'bo');
+    legend('OBS', 'Initial HofX', 'Location', 'best');
+    xlabel('Observation Number')
+    ylabel(Vname);
+    grid on;
+    title(['Group:  hofx0 member 1,   File:  ', untexlabel(MyFile)]);
+    if (~isempty(PNGprefix))
+      png_file = strcat(PNGprefix, shortname{n}, '_hofx0_1.png');
+      exportgraphics(gcf, png_file, 'resolution', 300);
+    end
+
+    if (any(strcmp({M.Variables.Name}, 'depth')))
+      figure;
+      h4 = plot(S.ObsValue{n}, S.depth, 'r+',                           ...
+                S.hofx0_1{n}, S.depth, 'bo');
+      legend('OBS', 'Initial HofX', 'Location', 'best');
+      xlabel(Vname);
+      ylabel('depth');
+      grid on;
+      title(['Group:  hofx0 member 1,   File:  ', untexlabel(MyFile)]);
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_hofx0_1_depth.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
+      end
+    else
+      figure;
+      h4 = plot3(S.longitude, S.latitude, S.ObsValue{n}, 'r+',         ...
+                 S.longitude, S.latitude, S.hofx0_1{n}, 'bo');
+      legend('OBS', 'Initial HofX', 'Location', 'best');
+      xlabel('Longitude');
+      ylabel('latitude');
+      zlabel(Vname);
+      grid on;
+      title(['Group:  hofx0 member 1,   File:  ', untexlabel(MyFile)]);
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_hofx0_1_lonlat.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
       end
     end
-  
+
   end
 
 %--------------------------------------------------------------------------
 %  Plot "hofx1" Group: Final HofX.
 %--------------------------------------------------------------------------
-  
+
   if ( (~got_grp && any(strcmp({I.Groups.Name}, 'hofx1'))) ||           ...
        (strcmp(group, 'hofx1')) )
 
@@ -350,11 +418,11 @@
     ylabel(Vname);
     grid on;
     title(['Group:  hofx1,   File:  ', untexlabel(MyFile)]);
-    if (wrtPNG)
-      png_file = strcat(shortname{n}, '_hofx1.png');
-      print(png_file, '-dpng', '-r300')
+    if (~isempty(PNGprefix))
+      png_file = strcat(PNGprefix, shortname{n}, '_hofx1.png');
+      exportgraphics(gcf, png_file, 'resolution', 300);
     end
- 
+
     if (any(strcmp({M.Variables.Name}, 'depth')))
       figure;
       h4 = plot(S.ObsValue{n}, S.depth, 'r+',                           ...
@@ -364,9 +432,9 @@
       ylabel('depth');
       grid on;
       title(['Group:  hofx1,   File:  ', untexlabel(MyFile)]);
-      if (wrtPNG)
-        png_file = strcat(shortname{n}, '_hofx1_depth.png');
-        print(png_file, '-dpng', '-r300')
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_hofx1_depth.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
       end
     else
       figure;
@@ -375,21 +443,70 @@
       legend('OBS', 'Final HofX', 'Location', 'best');
       xlabel('Longitude');
       ylabel('latitude');
-      zlabel(Vname); 
+      zlabel(Vname);
       grid on;
       title(['Group:  hofx1,   File:  ', untexlabel(MyFile)]);
-      if (wrtPNG)
-        png_file = strcat(shortname{n}, '_hofx1_lonlat.png');
-        print(png_file, '-dpng', '-r300')
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_hofx1_lonlat.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
+      end
+    end
+
+  end
+
+% Ensemble HofX.
+
+  if ( (~got_grp && any(strcmp({I.Groups.Name}, 'hofx1_1'))) ||         ...
+       (strcmp(group, 'hofx1_1')) )
+
+    figure;
+    h1 = plot(Nlocs, S.ObsValue{n}, 'r+',                               ...
+              Nlocs, S.hofx1_1{n}, 'bo');
+    legend('OBS', 'Final HofX', 'Location', 'best');
+    xlabel('Observation Number')
+    ylabel(Vname);
+    grid on;
+    title(['Group:  hofx1 member 1,   File:  ', untexlabel(MyFile)]);
+    if (~isempty(PNGprefix))
+      png_file = strcat(PNGprefix, shortname{n}, '_hofx1_1.png');
+      exportgraphics(gcf, png_file, 'resolution', 300);
+    end
+
+    if (any(strcmp({M.Variables.Name}, 'depth')))
+      figure;
+      h4 = plot(S.ObsValue{n}, S.depth, 'r+',                           ...
+                S.hofx1_1{n}, S.depth, 'bo');
+      legend('OBS', 'Final HofX', 'Location', 'best');
+      xlabel(Vname);
+      ylabel('depth');
+      grid on;
+      title(['Group:  hofx1 member 1,   File:  ', untexlabel(MyFile)]);
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_hofx1_1_depth.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
+      end
+    else
+      figure;
+      h4 = plot3(S.longitude, S.latitude, S.ObsValue{n}, 'r+',         ...
+                 S.longitude, S.latitude, S.hofx1_1{n}, 'bo');
+      legend('OBS', 'Final HofX', 'Location', 'best');
+      xlabel('Longitude');
+      ylabel('latitude');
+      zlabel(Vname);
+      grid on;
+      title(['Group:  hofx1 member 1,   File:  ', untexlabel(MyFile)]);
+      if (~isempty(PNGprefix))
+        png_file = strcat(PNGprefix, shortname{n}, '_hofx1_1_lonlat.png');
+        exportgraphics(gcf, png_file, 'resolution', 300);
       end
     end
-  
+
   end
 
 %--------------------------------------------------------------------------
 %  Plot observation minus analysis
 %--------------------------------------------------------------------------
-  
+
   if ( (~got_grp && any(strcmp({I.Groups.Name}, 'oman'))) ||            ...
        (strcmp(group, 'oman')) )
 
@@ -398,18 +515,18 @@
     xlabel('Observation Number')
     ylabel(Vname);
     grid on;
-    title(['oman = OBS - Analysis,   File:  ', untexlabel(MyFile)]);    
-    if (wrtPNG)
-      png_file = strcat(shortname{n}, '_oman.png');
-      print(png_file, '-dpng', '-r300')
+    title(['oman = OBS - Analysis,   File:  ', untexlabel(MyFile)]);
+    if (~isempty(PNGprefix))
+      png_file = strcat(PNGprefix, shortname{n}, '_oman.png');
+      exportgraphics(gcf, png_file, 'resolution', 300);
     end
-  
+
   end
 
 %--------------------------------------------------------------------------
 %  Plot observation minus background.
 %--------------------------------------------------------------------------
-  
+
   if ( (~got_grp && any(strcmp({I.Groups.Name}, 'oman'))) ||            ...
        (strcmp(group, 'ombg')) )
 
@@ -418,16 +535,16 @@
     xlabel('Observation Number')
     ylabel(Vname);
     grid on;
-    title(['ombg = OBS - Background,   File:  ', untexlabel(MyFile)]);    
-    if (wrtPNG)
-      png_file = strcat(shortname{n}, '_ombg.png');
-      print(png_file, '-dpng', '-r300')
+    title(['ombg = OBS - Background,   File:  ', untexlabel(MyFile)]);
+    if (~isempty(PNGprefix))
+      png_file = strcat(PNGprefix, shortname{n}, '_ombg.png');
+      exportgraphics(gcf, png_file, 'resolution', 300);
     end
-  
+
   end
-  
+
 % plot next state variable, if any.
-  
+
 end
 
 return
diff --git a/ioda/roms2ioda.m b/ioda/roms2ioda.m
index 3ef0a55..e100fd2 100644
--- a/ioda/roms2ioda.m
+++ b/ioda/roms2ioda.m
@@ -1,9 +1,10 @@
-function roms2ioda(ObsData, HisName, prefix, suffix)
+function roms2ioda(ObsData, HisName, prefix, suffix, varargin)
 
 %
 % ROMS2IODA:  Converts ROMS observation NetCDF to IODA NetCDF4 files
 %
-% roms2ioda(ObsData, HisName)
+% roms2ioda(ObsData, HisName, prefix, suffix, ...
+%           SSHareaAvg, SSHtimeAvg, UVtimeAvg)
 %
 % This function converts a ROMS 4D-Var observation NetCDF file into several
 % IODA NetCDF files. One file per each observation type is usually the way
@@ -13,6 +14,9 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 %
 % For example:           wc13_sst_20040103.nc4
 %
+% The area-averaged and time-averaged values will be used in ROMS H(x)
+% computations with area-averaged and/or time-averaged filters.
+%
 % On Input:
 %
 %    ObsData     ROMS standard 4D-Var observation NetCDF filename (string)
@@ -21,7 +25,7 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 %    HisName     ROMS application history NetCDF filename (string)
 %
 %    prefix      Output file prefix associated with application (string)
-%   
+%
 %    suffix      Output file suffix associated with date and time, use
 %                  YYYYMMDD or YYYYMMDDhh (numeric or string).
 %
@@ -31,12 +35,41 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 %                  For example, if suffix =  20040103  or
 %                                  suffix = '20040103' then
 %
-%            int64 dateTime(Location);
+%                  int64 dateTime(Location);
 %                  dateTime:units = "seconds since 2004-01-03T00:00:00Z";
 %
 %                  Therefore, input and output time of the observations
 %                  may have different time references!
 %
+%    SSHareaAvg  SSH observations area-averaged radius (km, OPTIONAL)
+%                  Use unique SSH observations (remove repetitive)
+%
+%                  float spatialAverage
+%
+%    SSHtimeAvg  SSH observations time-averaged window since reference
+%                  time (hours; OPTIONAL).
+%
+%                  For example 24 hours, 2004-01-03 to 2004-01-04
+%                                        2004-01-04 to 2004-01-05
+%                                        2004-01-05 to 2004-01-06
+%                                        and so on ...
+%
+%                  int64 dateTimeAverageStart(timeWindow);
+%
+%                  int64 dateTimeAverageEnd(timeWindow);
+%
+%    UVtimeAvg   CODAR observations time-averaged window since reference
+%                  time (hours; OPTIONAL)
+%
+%                  For example 24 hours, 2004-01-03 to 2004-01-04
+%                                        2004-01-04 to 2004-01-05
+%                                        2004-01-05 to 2004-01-06
+%                                        and so on ...
+%
+%                  int64 dateTimeAverageStart(timeWindow);
+%
+%                  int64 dateTimeAverageEnd(timeWindow);
+%
 % Usage: Example to convert WC13 4D-Var observation file to create the
 %        following IODA-2 files:
 %                                  wc13_adt_20040103.nc4
@@ -65,6 +98,25 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 
 % Initialize.
 
+switch numel(varargin)
+  case 0
+    SSHareaAvg = NaN;
+    SSHtimeAvg = NaN;
+    UVtimeAvg  = NaN;
+  case 1
+    SSHareaAvg = varargin{1};
+    SSHtimeAvg = NaN;
+    UVtimeAvg  = NaN;
+  case 2
+    SSHareaAvg = varargin{1};
+    SSHtimeAvg = varargin{2};
+    UVtimeAvg  = NaN;
+  case 3
+    SSHareaAvg = varargin{1};
+    SSHtimeAvg = varargin{2};
+    UVtimeAvg  = varargin{3};
+end
+
 if (ischar(ObsData))
   S = obs_read(ObsData);
 else
@@ -108,6 +160,10 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 S.Zgrid = S.depth;
 S = obs_k2z(S, HisName);
 
+% Estimate data assimilation cycle time-window.
+
+days_window = floor((max(S.time)-min(S.time))+0.5);
+
 %--------------------------------------------------------------------------
 % Inquire observation structure about type meassurent types.
 %--------------------------------------------------------------------------
@@ -130,6 +186,28 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 itemp = find(S.type == 6);
 isalt = find(S.type == 7);
 
+% Currently, one of the strategies to improve the impact of altimetry
+% track data is to repeat their observations several times over the
+% assimilation cycle, but with different errors.
+%
+% Identify such repetitive observations by setting their provenance
+% to negative values.
+
+if (~isempty(issh))
+  [~,IA,~] = unique(complex(S.lat(issh), S.lon(issh)));
+  if (~isempty(IA))
+    sshProv = -abs(S.provenance(issh));      % Set negative SSH provenance
+    for n = 1:length(IA)
+      sshProv(IA(n))= -sshProv(IA(n));       % Turn positive unique values
+    end
+    S.provenance(issh) = sshProv;            % overwrite SSH provenance
+  end
+
+  if (~isnan(SSHareaAvg) || ~isnan(SShtimeAbg))
+    issh = find(S.provenance(issh) > 0);
+  end
+end
+
 % Get provenance indices for each state variable.
 
 P = struct('ssh', [], 'uvel', [], 'vvel', [], 'temp', [], 'salt', []);
@@ -175,15 +253,31 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
     has_depth  = false;
     Obs = extract_observations(S, issh, DateTimeIODA, has_depth);
     Obs.ncfile         = [prefix '_adt_' suffix '.nc4'];
+    Obs.N              = G.N;
     Obs.nvars          = 1;
     Obs.units          = {'meter'};
     Obs.ncvname        = {'absoluteDynamicTopography'};
+    Obs.stateID        = 1;
+    Obs.areaAvgRadius  = SSHareaAvg;
+    Obs.timeAvgWindow  = SSHtimeAvg;
     Obs.variables_name = {'absolute_dynamic_topography'};
     Obs.datetime_ref   = DateTimeIODA;
     if (got_flagAtt)
       Obs.flag_values   = int32(P.flag_values(P.ssh));
       Obs.flag_meanings = string(join(P.flag_meanings(P.ssh)));
     end
+    if (~isnan(SSHareaAvg))
+      Obs.areaAvgRadius = SSHareaAvg * 1000;        % to meters
+    end
+    if (~isnan(SSHtimeAvg))
+      delta  = SSHtimeAvg * 3600;                    % to secods
+      window = days_window * 86400;
+      Tstr = 0:delta:window-delta;
+      Tend = delta:delta:window;
+      Obs.timeAvgBegin = Tstr;
+      Obs.timeAvgEnd   = Tend;
+      Obs.nwindow      = length(Tstr);
+    end
     create_ioda_obs(Obs);
     write_observations(Obs);
   end
@@ -191,15 +285,17 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 
 % Process sea surface temperature.
 
-if (got_temp)     
+if (got_temp)
   isst = find(S.type == 6 & S.Zgrid == G.N);
   if (~isempty(isst))
     has_depth = false;
     Obs = extract_observations(S, isst, DateTimeIODA, has_depth);
     Obs.ncfile         = [prefix '_sst_' suffix '.nc4'];
+    Obs.N              = G.N;
     Obs.nvars          = 1;
     Obs.units          = {'C'};
     Obs.ncvname        = {'seaSurfaceTemperature'};
+    Obs.stateID        = 6;
     Obs.variables_name = {'sea_surface_temperature'};
     Obs.datetime_ref   = DateTimeIODA;
     if (got_flagAtt)
@@ -219,14 +315,16 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 
 % Process sub-surface insitu temperature.
 
-if (got_temp)     
+if (got_temp)
   ktemp = find(S.type == 6 & S.Zgrid ~= G.N);
   if (~isempty(ktemp))
     has_depth = true;
     Obs = extract_observations(S, ktemp, DateTimeIODA, has_depth);
     Obs.ncfile         = [prefix '_temp_' suffix '.nc4'];
+    Obs.N              = G.N;
     Obs.nvars          = 1;
     Obs.units          = {'C'};
+    Obs.stateID        = 6;
     Obs.ncvname        = {'waterTemperature'};
     Obs.variables_name = {'sea_water_temperature'};
     Obs.datetime_ref   = DateTimeIODA;
@@ -234,12 +332,12 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
       sst_ind = contains(upper(P.flag_meanings(P.temp)), 'SST');
       if (any(sst_ind))                             % remove SST provenance
         P.temp(sst_ind)   = [];
-	Obs.flag_values   = int32(P.flag_values(P.temp));
+        Obs.flag_values   = int32(P.flag_values(P.temp));
         Obs.flag_meanings = string(join(P.flag_meanings(P.temp)));
-      else 
-	Obs.flag_values   = int32(P.flag_values(P.temp));
+      else
+        Obs.flag_values   = int32(P.flag_values(P.temp));
         Obs.flag_meanings = string(join(P.flag_meanings(P.temp)));
-      end      
+      end
     end
     create_ioda_obs(Obs);
     write_observations(Obs);
@@ -250,34 +348,36 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 % observations in the input NetCDF are insitu temperature. Below, they
 % are converted to potential temperature.
 
-if (got_temp)     
+if (got_temp)
   ktemp = find(S.type == 6 & S.Zgrid ~= G.N);
   if (~isempty(ktemp))
     has_depth = true;
     Obs = extract_observations(S, ktemp, DateTimeIODA, has_depth);
-    
+
     p = sw_pres(abs(Obs.depth), Obs.latitude);      % decibars
     s = ones(size(Obs.ObsValue)).*35.0;
     pr = zeros(size(Obs.ObsValue));
     ptemp = sw_ptmp(s, Obs.ObsValue, p, pr);
     Obs.ObsValue = ptemp;
-    
+
     Obs.ncfile         = [prefix '_ptemp_' suffix '.nc4'];
+    Obs.N              = G.N;
     Obs.nvars          = 1;
     Obs.units          = {'C'};
     Obs.ncvname        = {'waterPotentialTemperature'};
+    Obs.stateID        = 6;
     Obs.variables_name = {'sea_water_potential_temperature'};
     Obs.datetime_ref   = DateTimeIODA;
     if (got_flagAtt)
       sst_ind = contains(upper(P.flag_meanings(P.temp)), 'SST');
       if (any(sst_ind))                             % remove SST provenance
         P.temp(sst_ind)   = [];
-	Obs.flag_values   = int32(P.flag_values(P.temp));
+        Obs.flag_values   = int32(P.flag_values(P.temp));
         Obs.flag_meanings = string(join(P.flag_meanings(P.temp)));
-      else 
-	Obs.flag_values   = int32(P.flag_values(P.temp));
+      else
+        Obs.flag_values   = int32(P.flag_values(P.temp));
         Obs.flag_meanings = string(join(P.flag_meanings(P.temp)));
-      end      
+      end
     end
     create_ioda_obs(Obs);
     write_observations(Obs);
@@ -286,14 +386,16 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 
 % Process salinity.
 
-if (got_salt)     
+if (got_salt)
   if (~isempty(isalt))
     has_depth = true;
     Obs = extract_observations(S, isalt, DateTimeIODA, has_depth);
     Obs.ncfile         = [prefix '_salt_' suffix '.nc4'];
+    Obs.N              = G.N;
     Obs.nvars          = 1;
     Obs.units          = {'dimensionless'};
     Obs.ncvname        = {'salinity'};
+    Obs.stateID        = 7;
     Obs.variables_name = {'sea_water_salinity'};
     Obs.datetime_ref   = DateTimeIODA;
     if (got_flagAtt)
@@ -305,35 +407,98 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
   end
 end
 
+% Process gridded HF radar near surface (z=-2 m) velocities. The
+% CODAR velocity components in the native ROMS observation file
+% are rotated to grid curvilinear coordinates, we need to rotate
+% back to geographycal EAST and NORTH directions before they are
+% written to the IODA-type file.
+
+if (got_uvel && got_vvel)
+  if ~(isempty(iuvel) || isempty(iuvel))
+    has_depth = true;
+    Uobs = extract_observations(S, iuvel, DateTimeIODA, has_depth);
+    Vobs = extract_observations(S, ivvel, DateTimeIODA, has_depth);
+    Obs  = Uobs;
+    Obs.ObsError       = [];
+    Obs.ObsError{1}    = Uobs.ObsError;
+    Obs.ObsError{2}    = Vobs.ObsError;
+    Obs.ObsValue       = [];
+    Obs.ObsValue{1}    = Uobs.ObsValue;
+    Obs.ObsValue{2}    = Vobs.ObsValue;
+    Obs.PreQC          = [];
+    Obs.PreQC{1}       = Uobs.PreQC;
+    Obs.PreQC{2}       = Uobs.PreQC;
+    Obs.ncfile         = [prefix '_uv_codar_' suffix '.nc4'];
+    Obs.N              = G.N;
+    Obs.nvars          = 2;
+    Obs.units          = {'m s-1', 'm s-1'};
+    Obs.ncvname        = {'waterZonalVelocity',                      ...
+                          'waterMeridionalVelocity'};
+    Obs.stateID        = [4, 5];
+    Obs.timeAvgWindow  = UVtimeAvg;
+    Obs.variables_name = {'eastward_sea_water_velocity',             ...
+                          'meridional_sea_water_velocity'};
+    Obs.datetime_ref   = DateTimeIODA;
+    if (got_flagAtt)
+      Obs.flag_values   = int32(P.flag_values(P.uvel));
+      Obs.flag_meanings = string(join(P.flag_meanings(P.uvel)));
+    end
+    if (~isnan(UVtimeAvg))
+      delta  = UVtimeAvg * 3600;                    % to secods
+      window = days_window * 86400;
+      Tstr = 0:delta:window-delta;
+      Tend = delta:delta:window;
+      Obs.timeAvgBegin = Tstr;
+      Obs.timeAvgEnd   = Tend;
+      Obs.nwindow      = length(Tstr);
+    end
+    Obs = rotate_codar(Obs, G);
+    create_ioda_obs(Obs);
+    write_observations(Obs);
+  end
+end
+
 return
 
 %--------------------------------------------------------------------------
 % Extract requested observation type from structure.
 
-function [Obs] = extract_observations(S, ind, DateTimeIODA, has_depth);
+function [Obs] = extract_observations(S, ind, DateTimeIODA, has_depth)
 
 % Initialize structure.
 
 Obs = struct('ncfile'        , [],                                      ...
              'source'        , [],                                      ...
+             'N'             , [],                                      ...
              'nlocs'         , [],                                      ...
              'nobs'          , [],                                      ...
+             'nsurvey'       , [],                                      ...
              'nvars'         , [],                                      ...
+             'nwindow'       , NaN,                                     ...
              'units'         , [],                                      ...
              'ncvname'       , [],                                      ...
              'variable_names', [],                                      ...
              'TimeIODA'      , [],                                      ...
-	     'DateIODA'      , [],                                      ...
+             'DateIODA'      , [],                                      ...
              'dateTime'      , [],                                      ...
+             'surveyTime'    , [],                                      ...
+             'surveyIndex'   , [],                                      ...
              'date_time'     , [],                                      ...
              'depth'         , [],                                      ...
              'latitude'      , [],                                      ...
              'longitude'     , [],                                      ...
              'provenance'    , [],                                      ...
+             'stateID'       , [],                                      ...
              'sequenceNumber', [],                                      ...
+             'areaAvgRadius' , NaN,                                     ...
+             'timeAvgBegin'  , NaN,                                     ...
+             'timeAvgEnd'    , NaN,                                     ...
+             'x_grid'        , [],                                      ...
+             'y_grid'        , [],                                      ...
+             'z_grid'        , [],                                      ...
              'ObsError'      , [],                                      ...
-	     'ObsValue'      , [],                                      ...
-	     'PreQC'         , []);
+             'ObsValue'      , [],                                      ...
+             'PreQC'         , []);
 
 % Fill some values in the structure.
 
@@ -349,20 +514,30 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 Obs.longitude      = S.lon(ind);
 Obs.provenance     = S.provenance(ind);
 Obs.sequenceNumber = int32(1:Obs.nlocs);
+Obs.x_grid         = S.Xgrid(ind);
+Obs.y_grid         = S.Ygrid(ind);
+if (has_depth)
+  Obs.z_grid       = S.Zgrid(ind);
+end
 Obs.ObsError       = sqrt(S.error(ind));              % standard deviation
 Obs.ObsValue       = S.value(ind);
 Obs.PreQC          = int32(zeros(size(Obs.longitude)));
 
 % Compute dateTime seconds from IODA file reference time.
 
-ioda_epoch     = datenum(num2str(DateTimeIODA), 'yyyymmddHH');
+ioda_epoch      = datenum(num2str(DateTimeIODA), 'yyyymmddHH');
+
+time_days       = S.reference_time + S.time(ind);
+ioda_time_days  = time_days - ioda_epoch;
 
-time_days      = S.reference_time + S.time(ind);
-ioda_time_days = time_days - ioda_epoch;        
+Obs.dateTime    = int64(ioda_time_days * 86400);      % seconds
 
-Obs.dateTime   = int64(ioda_time_days * 86400);        % seconds
+[survey,IA,IC]  = unique(Obs.dateTime, 'stable');
+Obs.surveyTime  = int64(survey);
+Obs.surveyIndex = int32(IA);
+Obs.nsurvey     = length(Obs.surveyTime);
 
-Obs.date_time  = cellstr(datestr(ioda_epoch + ioda_time_days,          ...
+Obs.date_time   = cellstr(datestr(ioda_epoch + ioda_time_days,          ...
                                  'yyyy-mm-ddTHH:MM:SSZ'));
 
 % Set IODA file time reference global attributes.
@@ -372,17 +547,57 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 
 return
 
+%--------------------------------------------------------------------------
+% Extract requested observation type from structure.
+
+function [Obs] = rotate_codar (S, G)
+
+% Interpolate ROMS curvilinear angle to CODAR observation locations.
+% Use scatteredInterpolant since horizontal grid coordinates are not
+% plaid.
+
+F = scatteredInterpolant(G.lon_rho(:), G.lat_rho(:), G.angle(:));
+
+angle = F(S.longitude, S.latitude);
+
+% Rotate velocity component to geographical EAST and NORTH.
+
+Urot = S.ObsValue{1};
+Vrot = S.ObsValue{2};
+
+cos_angle = cos(angle);
+sin_angle = sin(angle);
+
+Ugeo = Urot .* cos_angle - Vrot .* sin_angle;
+Vgeo = Vrot .* cos_angle + Urot .* sin_angle;
+
+% Replace velocity components.
+
+Obs = S;
+
+Obs.ObsValue{1} = Ugeo;
+Obs.ObsValue{2} = Vgeo;
+
+return
+
 %--------------------------------------------------------------------------
 % Writes requested observation type data to ouput NetCDF file.
 
-function write_observations(Obs);
+function write_observations(Obs)
 
 disp(['*** Writing  observations file:  ', Obs.ncfile]);
 
 % Write out 'MetaData' Group variables.
 
 ncwrite(Obs.ncfile, 'MetaData/dateTime', int64(Obs.dateTime));
-ncwrite(Obs.ncfile, 'MetaData/date_time', Obs.date_time);
+if (~isnan(Obs.timeAvgBegin))
+  ncwrite(Obs.ncfile, 'MetaData/dateTimeAverageBegin',                  ...
+          int64(Obs.timeAvgBegin));
+end
+if (~isnan(Obs.timeAvgEnd))
+  ncwrite(Obs.ncfile, 'MetaData/dateTimeAverageEnd',                    ...
+          int64(Obs.timeAvgEnd));
+end
 if (~isempty(Obs.depth))
   ncwrite(Obs.ncfile, 'MetaData/depth', Obs.depth);
 end
@@ -391,27 +606,53 @@ function roms2ioda(ObsData, HisName, prefix, suffix)
 if (~isempty(Obs.provenance))
   ncwrite(Obs.ncfile, 'MetaData/provenance', Obs.provenance);
 end
+ncwrite(Obs.ncfile, 'MetaData/stateID', int32(Obs.stateID));
 ncwrite(Obs.ncfile, 'MetaData/sequenceNumber', Obs.sequenceNumber);
+if (~isnan(Obs.areaAvgRadius))
+  ncwrite(Obs.ncfile, 'MetaData/spatialAverage', Obs.areaAvgRadius);
+end
+ncwrite(Obs.ncfile, 'MetaData/surveyIndex', int32(Obs.surveyIndex));
+ncwrite(Obs.ncfile, 'MetaData/surveyTime', int64(Obs.surveyTime));
 ncwrite(Obs.ncfile, 'MetaData/variables_name', Obs.variables_name);
+ncwrite(Obs.ncfile, 'MetaData/x_grid', Obs.x_grid);
+ncwrite(Obs.ncfile, 'MetaData/y_grid', Obs.y_grid);
+if (~isempty(Obs.z_grid))
+  ncwrite(Obs.ncfile, 'MetaData/z_grid', Obs.z_grid);
+end
 
 % Write out 'ObsError' Group variables.
 
-for i = 1:Obs.nvars
-  Vname = ['/ObsError/' Obs.ncvname{i}];
-  ncwrite(Obs.ncfile, Vname, Obs.ObsError{i});
+if (iscell(Obs.ObsError))
+  for i = 1:Obs.nvars
+    Vname = ['/ObsError/' Obs.ncvname{i}];
+    ncwrite(Obs.ncfile, Vname, Obs.ObsError{i});
+  end
+else
+  Vname = ['/ObsError/' char(Obs.ncvname)];
+  ncwrite(Obs.ncfile, Vname, Obs.ObsError);
 end
 
 % Write out 'ObsValue' Group variables.
 
-for i = 1:Obs.nvars
-  Vname = ['/ObsValue/' Obs.ncvname{i}];
+if (iscell(Obs.ObsValue))
+  for i = 1:Obs.nvars
+    Vname = ['/ObsValue/' Obs.ncvname{i}];
+    ncwrite(Obs.ncfile, Vname, Obs.ObsValue{i});
+  end
+else
+  Vname = ['/ObsValue/' char(Obs.ncvname)];
   ncwrite(Obs.ncfile, Vname, Obs.ObsValue);
 end
 
 % Write out 'ObsQC' Group variables.
 
-for i = 1:Obs.nvars
-  Vname = ['/PreQC/' Obs.ncvname{i}];
+if (iscell(Obs.PreQC))
+  for i = 1:Obs.nvars
+    Vname = ['/PreQC/' Obs.ncvname{i}];
+    ncwrite(Obs.ncfile, Vname, Obs.PreQC{i});
+  end
+else
+  Vname = ['/PreQC/' char(Obs.ncvname)];
   ncwrite(Obs.ncfile, Vname, Obs.PreQC);
 end
 
diff --git a/utility/plot_section.m b/utility/plot_section.m
index af41736..4b73715 100644
--- a/utility/plot_section.m
+++ b/utility/plot_section.m
@@ -57,13 +57,13 @@
 %    Licensed under a MIT/X style license                                 %
 %    See License_ROMS.md                            Hernan G. Arango      %
 %=========================================================================%
-  
+
 % Initialize.
 
 if (orient == 'c')
   F = struct('ncname'     , [], 'Vname'     , [],                       ...
              'Tindex'     , [], 'Tname'     , [], 'Tstring'   , [],     ...
-	     'column'     , [],                                         ...
+             'column'     , [],                                         ...
              'value'      , [], 'min'       , [], 'max'       , [],     ...
              'X'          , [], 'Z'         , [],                       ...
              'Imin'       , [], 'Jmin'      , [],                       ...
@@ -72,7 +72,7 @@
 else
   F = struct('ncname'     , [], 'Vname'     , [],                       ...
              'Tindex'     , [], 'Tname'     , [], 'Tstring'   , [],     ...
-	     'row'        , [],                                         ...
+             'row'        , [],                                         ...
              'value'      , [], 'min'       , [], 'max'       , [],     ...
              'X'          , [], 'Z'         , [],                       ...
              'Imin'       , [], 'Jmin'      , [],                       ...
@@ -180,7 +180,7 @@
 
 if (getdata)
 
-  V = nc_vnames(Hname);  
+  V = nc_vnames(Hname);
 
   switch (Vname)
     case {'Hz', 'hz'}
@@ -189,7 +189,7 @@
       else
         I = nc_vinfo(Hname, 'z_w');
       end
-    otherwise 
+    otherwise
       I = nc_vinfo(Hname, Vname);
   end
 
@@ -336,7 +336,7 @@
     N = size(Z,3);
   else
     error([' PLOT_SECTION - cannot plot a section for a 2D field: ''',  ...
-	   Vname, ''])
+           Vname, ''])
   end
 
   if (isfield(G,Mname))
@@ -352,8 +352,10 @@
   end
 
   if (~G.spherical)
-    X = 0.001 .* X;
-    Y = 0.001 .* Y;
+    if (max(X(:)) > 5000 || max(Y(:)) > 5000)
+      X = 0.001 .* X;                           % scale to kilometers
+      Y = 0.001 .* Y;
+    end
   end
 end
 
@@ -401,7 +403,7 @@
       Z = G.z_r;
    otherwise
       field = nc_read(Hname,Vname,Tindex,ReplaceValue,PreserveType);
-  end  
+  end
 
 end
 
@@ -435,7 +437,7 @@
       x = squeeze(G.x_rho(:,index));
     end
     z = -squeeze(G.h(:,index));
-end    
+end
 
 V = V .* scale;
 
@@ -493,7 +495,7 @@
   hold on;
   area(x, z, min(z), 'FaceColor', Land, 'EdgeColor', Land);
   hold off;
-  
+
   if (~isempty(Tname))
     ht = title([untexlabel(Vname), ':', blanks(4),                      ...
                 'Record = ', num2str(Tindex), ',', blanks(4),           ...
diff --git a/utility/scoord.m b/utility/scoord.m
index 4e1b326..fd273be 100644
--- a/utility/scoord.m
+++ b/utility/scoord.m
@@ -97,7 +97,7 @@
 if (Vtransform < 1 | Vtransform > 2),
   disp(' ');
   disp([setstr(7),'*** Error:  SCOORD - Illegal parameter Vtransform = ' ...
-        num2str(Vtransfrom), setstr(7)]);
+        num2str(Vtransform), setstr(7)]);
   return
 end,