diff --git a/docs/iris/gallery_code/meteorology/plot_lagged_ensemble.py b/docs/iris/gallery_code/meteorology/plot_lagged_ensemble.py
index 5f2ab724b3..cb82a663d4 100644
--- a/docs/iris/gallery_code/meteorology/plot_lagged_ensemble.py
+++ b/docs/iris/gallery_code/meteorology/plot_lagged_ensemble.py
@@ -41,7 +41,7 @@ def realization_metadata(cube, field, fname):
         import iris.coords
 
         realization_coord = iris.coords.AuxCoord(
-            np.int32(realization_number), "realization"
+            np.int32(realization_number), "realization", units="1"
         )
         cube.add_aux_coord(realization_coord)
 
diff --git a/docs/iris/src/userguide/navigating_a_cube.rst b/docs/iris/src/userguide/navigating_a_cube.rst
index 055617e047..581d1a67cf 100644
--- a/docs/iris/src/userguide/navigating_a_cube.rst
+++ b/docs/iris/src/userguide/navigating_a_cube.rst
@@ -229,7 +229,7 @@ by field basis *before* they are automatically merged together:
         # Add our own realization coordinate if it doesn't already exist.
         if not cube.coords('realization'):
             realization = np.int32(filename[-6:-3])
-            ensemble_coord = icoords.AuxCoord(realization, standard_name='realization')
+            ensemble_coord = icoords.AuxCoord(realization, standard_name='realization', units="1")
             cube.add_aux_coord(ensemble_coord)
 
     filename = iris.sample_data_path('GloSea4', '*.pp')
diff --git a/docs/iris/src/whatsnew/contributions_3.0.0/incompatiblechange_2020-May-15_change_default_unit_loading.txt b/docs/iris/src/whatsnew/contributions_3.0.0/incompatiblechange_2020-May-15_change_default_unit_loading.txt
new file mode 100644
index 0000000000..be048990f3
--- /dev/null
+++ b/docs/iris/src/whatsnew/contributions_3.0.0/incompatiblechange_2020-May-15_change_default_unit_loading.txt
@@ -0,0 +1 @@
+* When loading data from netcdf-CF files, where a variable has no "units" property, the corresponding Iris object will have "units='unknown'". Prior to Iris 3.0, these cases defaulted to "units='1'".
\ No newline at end of file
diff --git a/lib/iris/analysis/__init__.py b/lib/iris/analysis/__init__.py
index c94baf0a92..a1e56533fd 100644
--- a/lib/iris/analysis/__init__.py
+++ b/lib/iris/analysis/__init__.py
@@ -803,7 +803,9 @@ def post_process(self, collapsed_cube, data_result, coords, **kwargs):
         # order cube.
         for point in points:
             cube = collapsed_cube.copy()
-            coord = iris.coords.AuxCoord(point, long_name=coord_name)
+            coord = iris.coords.AuxCoord(
+                point, long_name=coord_name, units="percent"
+            )
             cube.add_aux_coord(coord)
             cubes.append(cube)
 
diff --git a/lib/iris/coords.py b/lib/iris/coords.py
index 09569f4ba5..d2dcd35c92 100644
--- a/lib/iris/coords.py
+++ b/lib/iris/coords.py
@@ -65,7 +65,7 @@ def __init__(
         standard_name=None,
         long_name=None,
         var_name=None,
-        units="no-unit",
+        units=None,
         attributes=None,
     ):
         """
@@ -687,7 +687,7 @@ def __init__(
         standard_name=None,
         long_name=None,
         var_name=None,
-        units="no-unit",
+        units=None,
         attributes=None,
     ):
         """
@@ -793,7 +793,7 @@ def __init__(
         standard_name=None,
         long_name=None,
         var_name=None,
-        units="1",
+        units=None,
         attributes=None,
         measure=None,
     ):
@@ -1283,7 +1283,7 @@ def __init__(
         standard_name=None,
         long_name=None,
         var_name=None,
-        units="1",
+        units=None,
         bounds=None,
         attributes=None,
         coord_system=None,
@@ -2262,7 +2262,7 @@ def from_regular(
         standard_name=None,
         long_name=None,
         var_name=None,
-        units="1",
+        units=None,
         attributes=None,
         coord_system=None,
         circular=False,
@@ -2325,7 +2325,7 @@ def __init__(
         standard_name=None,
         long_name=None,
         var_name=None,
-        units="1",
+        units=None,
         bounds=None,
         attributes=None,
         coord_system=None,
diff --git a/lib/iris/fileformats/_pyke_rules/fc_rules_cf.krb b/lib/iris/fileformats/_pyke_rules/fc_rules_cf.krb
index ad2c181b0b..2afc823795 100644
--- a/lib/iris/fileformats/_pyke_rules/fc_rules_cf.krb
+++ b/lib/iris/fileformats/_pyke_rules/fc_rules_cf.krb
@@ -1195,6 +1195,8 @@ fc_extras
     UD_UNITS_LON = ['degrees_east', 'degree_east', 'degree_e', 'degrees_e',
                     'degreee', 'degreese', 'degrees', 'degrees east',
                     'degree east', 'degree e', 'degrees e']
+    UNKNOWN_UNIT_STRING = "?"
+    NO_UNIT_STRING = "-"
 
     #
     # CF Dimensionless Vertical Coordinates
@@ -1651,9 +1653,9 @@ fc_extras
 
     ################################################################################
     def get_attr_units(cf_var, attributes):
-        attr_units = getattr(cf_var, CF_ATTR_UNITS, cf_units._UNIT_DIMENSIONLESS)
+        attr_units = getattr(cf_var, CF_ATTR_UNITS, UNKNOWN_UNIT_STRING)
         if not attr_units:
-            attr_units = '1'
+            attr_units = UNKNOWN_UNIT_STRING
 
         # Sanitise lat/lon units.
         if attr_units in UD_UNITS_LAT or attr_units in UD_UNITS_LON:
@@ -1668,10 +1670,10 @@ fc_extras
                 cf_var.cf_name, attr_units)
             warnings.warn(msg)
             attributes['invalid_units'] = attr_units
-            attr_units = cf_units._UNKNOWN_UNIT_STRING
+            attr_units = UNKNOWN_UNIT_STRING
 
         if np.issubdtype(cf_var.dtype, np.str_):
-            attr_units = cf_units._NO_UNIT_STRING
+            attr_units = NO_UNIT_STRING
 
         # Get any assoicated calendar for a time reference coordinate.
         if cf_units.as_unit(attr_units).is_time_reference():
diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py
index 867b0c9faf..f34dc45e72 100644
--- a/lib/iris/fileformats/netcdf.py
+++ b/lib/iris/fileformats/netcdf.py
@@ -22,6 +22,7 @@
 import warnings
 
 import dask.array as da
+import cf_units
 import netCDF4
 import numpy as np
 import numpy.ma as ma
@@ -1760,7 +1761,7 @@ def _inner_create_cf_cellmeasure_or_ancil_variable(
         # Add the data to the CF-netCDF variable.
         cf_var[:] = data
 
-        if dimensional_metadata.units != "unknown":
+        if dimensional_metadata.units.is_udunits():
             _setncattr(cf_var, "units", str(dimensional_metadata.units))
 
         if dimensional_metadata.standard_name is not None:
@@ -1926,7 +1927,7 @@ def _create_cf_coord_variable(self, cube, dimension_names, coord):
         # Deal with CF-netCDF units and standard name.
         standard_name, long_name, units = self._cf_coord_identity(coord)
 
-        if units != "unknown":
+        if cf_units.as_unit(units).is_udunits():
             _setncattr(cf_var, "units", units)
 
         if standard_name is not None:
@@ -2371,7 +2372,7 @@ def store(data, cf_var, fill_value):
         if cube.long_name:
             _setncattr(cf_var, "long_name", cube.long_name)
 
-        if cube.units != "unknown":
+        if cube.units.is_udunits():
             _setncattr(cf_var, "units", str(cube.units))
 
         # Add the CF-netCDF calendar attribute.
diff --git a/lib/iris/fileformats/pp_load_rules.py b/lib/iris/fileformats/pp_load_rules.py
index b6aacb382f..53d9f4dc35 100644
--- a/lib/iris/fileformats/pp_load_rules.py
+++ b/lib/iris/fileformats/pp_load_rules.py
@@ -147,6 +147,7 @@ def _convert_vertical_coords(
             model_level_number,
             standard_name="model_level_number",
             attributes={"positive": "down"},
+            units="1",
         )
         coords_and_dims.append((coord, dim))
 
@@ -197,6 +198,7 @@ def _convert_vertical_coords(
                 model_level_number,
                 long_name="soil_model_level_number",
                 attributes={"positive": "down"},
+                units="1",
             )
             coords_and_dims.append((coord, dim))
         elif np.any(brsvd1 != brlev):
@@ -235,6 +237,7 @@ def _convert_vertical_coords(
             model_level_number,
             standard_name="model_level_number",
             attributes={"positive": "up"},
+            units="1",
         )
         level_pressure = _dim_or_aux(
             bhlev,
@@ -243,7 +246,10 @@ def _convert_vertical_coords(
             bounds=np.vstack((bhrlev, brsvd2)).T,
         )
         sigma = AuxCoord(
-            blev, long_name="sigma", bounds=np.vstack((brlev, brsvd1)).T
+            blev,
+            long_name="sigma",
+            bounds=np.vstack((brlev, brsvd1)).T,
+            units="1",
         )
         coords_and_dims.extend(
             [(model_level_number, dim), (level_pressure, dim), (sigma, dim)]
@@ -265,6 +271,7 @@ def _convert_vertical_coords(
             model_level_number,
             standard_name="model_level_number",
             attributes={"positive": "up"},
+            units="1",
         )
         level_height = _dim_or_aux(
             blev,
@@ -274,7 +281,10 @@ def _convert_vertical_coords(
             attributes={"positive": "up"},
         )
         sigma = AuxCoord(
-            bhlev, long_name="sigma", bounds=np.vstack((bhrlev, brsvd2)).T
+            bhlev,
+            long_name="sigma",
+            bounds=np.vstack((bhrlev, brsvd2)).T,
+            units="1",
         )
         coords_and_dims.extend(
             [(model_level_number, dim), (level_height, dim), (sigma, dim)]
@@ -846,7 +856,7 @@ def _convert_scalar_realization_coords(lbrsvd4):
     coords_and_dims = []
     if lbrsvd4 != 0:
         coords_and_dims.append(
-            (DimCoord(lbrsvd4, standard_name="realization"), None)
+            (DimCoord(lbrsvd4, standard_name="realization", units="1"), None)
         )
     return coords_and_dims
 
@@ -1078,7 +1088,7 @@ def _all_other_rules(f):
         and f.lbmon == f.lbmond
     ):
         aux_coords_and_dims.append(
-            (AuxCoord(f.lbmon, long_name="month_number"), None)
+            (AuxCoord(f.lbmon, long_name="month_number", units="1"), None)
         )
         aux_coords_and_dims.append(
             (
diff --git a/lib/iris/tests/integration/test_netcdf.py b/lib/iris/tests/integration/test_netcdf.py
index 8c6e0f6659..267e5beb50 100644
--- a/lib/iris/tests/integration/test_netcdf.py
+++ b/lib/iris/tests/integration/test_netcdf.py
@@ -81,7 +81,9 @@ def test_hybrid_height_and_pressure(self):
                 1200.0, long_name="level_pressure", units="hPa"
             )
         )
-        cube.add_aux_coord(iris.coords.DimCoord(0.5, long_name="other sigma"))
+        cube.add_aux_coord(
+            iris.coords.DimCoord(0.5, long_name="other sigma", units="1")
+        )
         cube.add_aux_coord(
             iris.coords.DimCoord(
                 1000.0, long_name="surface_air_pressure", units="hPa"
diff --git a/lib/iris/tests/integration/test_pp.py b/lib/iris/tests/integration/test_pp.py
index 6fbf180ac5..b9b096d782 100644
--- a/lib/iris/tests/integration/test_pp.py
+++ b/lib/iris/tests/integration/test_pp.py
@@ -299,7 +299,7 @@ def test_hybrid_height_with_non_standard_coords(self):
         delta_lower, delta, delta_upper = 150, 200, 250
 
         cube = Cube(np.zeros((ny, nx)), "air_temperature")
-        level_coord = AuxCoord(0, "model_level_number")
+        level_coord = AuxCoord(0, "model_level_number", units="1")
         cube.add_aux_coord(level_coord)
         delta_coord = AuxCoord(
             delta,
@@ -308,7 +308,10 @@ def test_hybrid_height_with_non_standard_coords(self):
             units="m",
         )
         sigma_coord = AuxCoord(
-            sigma, bounds=[[sigma_lower, sigma_upper]], long_name="mavis"
+            sigma,
+            bounds=[[sigma_lower, sigma_upper]],
+            long_name="mavis",
+            units="1",
         )
         surface_altitude_coord = AuxCoord(
             np.zeros((ny, nx)), "surface_altitude", units="m"
@@ -343,7 +346,7 @@ def test_hybrid_pressure_with_non_standard_coords(self):
         delta_lower, delta, delta_upper = 0.15, 0.2, 0.25
 
         cube = Cube(np.zeros((ny, nx)), "air_temperature")
-        level_coord = AuxCoord(0, "model_level_number")
+        level_coord = AuxCoord(0, "model_level_number", units="1")
         cube.add_aux_coord(level_coord)
         delta_coord = AuxCoord(
             delta,
@@ -352,7 +355,10 @@ def test_hybrid_pressure_with_non_standard_coords(self):
             units="Pa",
         )
         sigma_coord = AuxCoord(
-            sigma, bounds=[[sigma_lower, sigma_upper]], long_name="mavis"
+            sigma,
+            bounds=[[sigma_lower, sigma_upper]],
+            long_name="mavis",
+            units="1",
         )
         surface_air_pressure_coord = AuxCoord(
             np.zeros((ny, nx)), "surface_air_pressure", units="Pa"
diff --git a/lib/iris/tests/results/analysis/first_quartile_foo_1d.cml b/lib/iris/tests/results/analysis/first_quartile_foo_1d.cml
index a9e69c291e..f027f2d9f8 100644
--- a/lib/iris/tests/results/analysis/first_quartile_foo_1d.cml
+++ b/lib/iris/tests/results/analysis/first_quartile_foo_1d.cml
@@ -6,7 +6,7 @@
         
       
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/analysis/first_quartile_foo_1d_fast_percentile.cml b/lib/iris/tests/results/analysis/first_quartile_foo_1d_fast_percentile.cml
index a9e69c291e..f027f2d9f8 100644
--- a/lib/iris/tests/results/analysis/first_quartile_foo_1d_fast_percentile.cml
+++ b/lib/iris/tests/results/analysis/first_quartile_foo_1d_fast_percentile.cml
@@ -6,7 +6,7 @@
         
       
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/analysis/first_quartile_foo_2d.cml b/lib/iris/tests/results/analysis/first_quartile_foo_2d.cml
index 34c9e746f6..1bc809ce63 100644
--- a/lib/iris/tests/results/analysis/first_quartile_foo_2d.cml
+++ b/lib/iris/tests/results/analysis/first_quartile_foo_2d.cml
@@ -11,7 +11,7 @@
         
       
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/analysis/first_quartile_foo_2d_fast_percentile.cml b/lib/iris/tests/results/analysis/first_quartile_foo_2d_fast_percentile.cml
index 34c9e746f6..1bc809ce63 100644
--- a/lib/iris/tests/results/analysis/first_quartile_foo_2d_fast_percentile.cml
+++ b/lib/iris/tests/results/analysis/first_quartile_foo_2d_fast_percentile.cml
@@ -11,7 +11,7 @@
         
       
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/analysis/first_quartile_foo_bar_2d.cml b/lib/iris/tests/results/analysis/first_quartile_foo_bar_2d.cml
index b3f135cede..cadd1e8b65 100644
--- a/lib/iris/tests/results/analysis/first_quartile_foo_bar_2d.cml
+++ b/lib/iris/tests/results/analysis/first_quartile_foo_bar_2d.cml
@@ -9,7 +9,7 @@
         
       
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/analysis/first_quartile_foo_bar_2d_fast_percentile.cml b/lib/iris/tests/results/analysis/first_quartile_foo_bar_2d_fast_percentile.cml
index b3f135cede..cadd1e8b65 100644
--- a/lib/iris/tests/results/analysis/first_quartile_foo_bar_2d_fast_percentile.cml
+++ b/lib/iris/tests/results/analysis/first_quartile_foo_bar_2d_fast_percentile.cml
@@ -9,7 +9,7 @@
         
       
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/analysis/last_quartile_foo_3d_masked.cml b/lib/iris/tests/results/analysis/last_quartile_foo_3d_masked.cml
index 80fab0e150..059541e208 100644
--- a/lib/iris/tests/results/analysis/last_quartile_foo_3d_masked.cml
+++ b/lib/iris/tests/results/analysis/last_quartile_foo_3d_masked.cml
@@ -9,7 +9,7 @@
         
       
       
-        
+        
       
       
         
diff --git a/lib/iris/tests/results/analysis/last_quartile_foo_3d_notmasked.cml b/lib/iris/tests/results/analysis/last_quartile_foo_3d_notmasked.cml
index 80fab0e150..059541e208 100644
--- a/lib/iris/tests/results/analysis/last_quartile_foo_3d_notmasked.cml
+++ b/lib/iris/tests/results/analysis/last_quartile_foo_3d_notmasked.cml
@@ -9,7 +9,7 @@
         
       
       
-        
+        
       
       
         
diff --git a/lib/iris/tests/results/analysis/last_quartile_foo_3d_notmasked_fast_percentile.cml b/lib/iris/tests/results/analysis/last_quartile_foo_3d_notmasked_fast_percentile.cml
index 80fab0e150..059541e208 100644
--- a/lib/iris/tests/results/analysis/last_quartile_foo_3d_notmasked_fast_percentile.cml
+++ b/lib/iris/tests/results/analysis/last_quartile_foo_3d_notmasked_fast_percentile.cml
@@ -9,7 +9,7 @@
         
       
       
-        
+        
       
       
         
diff --git a/lib/iris/tests/results/analysis/third_quartile_foo_1d.cml b/lib/iris/tests/results/analysis/third_quartile_foo_1d.cml
index b14c51cfb3..038e7c8668 100644
--- a/lib/iris/tests/results/analysis/third_quartile_foo_1d.cml
+++ b/lib/iris/tests/results/analysis/third_quartile_foo_1d.cml
@@ -6,7 +6,7 @@
         
       
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/analysis/third_quartile_foo_1d_fast_percentile.cml b/lib/iris/tests/results/analysis/third_quartile_foo_1d_fast_percentile.cml
index b14c51cfb3..038e7c8668 100644
--- a/lib/iris/tests/results/analysis/third_quartile_foo_1d_fast_percentile.cml
+++ b/lib/iris/tests/results/analysis/third_quartile_foo_1d_fast_percentile.cml
@@ -6,7 +6,7 @@
         
       
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/coord_api/minimal.xml b/lib/iris/tests/results/coord_api/minimal.xml
index a35c93dc68..8f93fb6376 100644
--- a/lib/iris/tests/results/coord_api/minimal.xml
+++ b/lib/iris/tests/results/coord_api/minimal.xml
@@ -1,2 +1,2 @@
 
-
+
diff --git a/lib/iris/tests/results/integration/climatology/TestClimatology/reference_simpledata.cdl b/lib/iris/tests/results/integration/climatology/TestClimatology/reference_simpledata.cdl
index 1740926645..1f6bc36832 100644
--- a/lib/iris/tests/results/integration/climatology/TestClimatology/reference_simpledata.cdl
+++ b/lib/iris/tests/results/integration/climatology/TestClimatology/reference_simpledata.cdl
@@ -17,11 +17,11 @@ variables:
 	double time_climatology(time, bnds) ;
 	double latitude(latitude) ;
 		latitude:axis = "Y" ;
-		latitude:units = "1" ;
+		latitude:units = "degrees_north" ;
 		latitude:standard_name = "latitude" ;
 	double longitude(longitude) ;
 		longitude:axis = "X" ;
-		longitude:units = "1" ;
+		longitude:units = "degrees_east" ;
 		longitude:standard_name = "longitude" ;
 
 // global attributes:
diff --git a/lib/iris/tests/results/netcdf/int64_auxiliary_coord_netcdf3.cml b/lib/iris/tests/results/netcdf/int64_auxiliary_coord_netcdf3.cml
index 39cb8f2950..e48cf41d2a 100644
--- a/lib/iris/tests/results/netcdf/int64_auxiliary_coord_netcdf3.cml
+++ b/lib/iris/tests/results/netcdf/int64_auxiliary_coord_netcdf3.cml
@@ -6,7 +6,7 @@
     
     
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/netcdf/int64_dimension_coord_netcdf3.cml b/lib/iris/tests/results/netcdf/int64_dimension_coord_netcdf3.cml
index 1c59fc947e..78fec459e9 100644
--- a/lib/iris/tests/results/netcdf/int64_dimension_coord_netcdf3.cml
+++ b/lib/iris/tests/results/netcdf/int64_dimension_coord_netcdf3.cml
@@ -6,7 +6,7 @@
     
     
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/netcdf/netcdf_cell_methods.cml b/lib/iris/tests/results/netcdf/netcdf_cell_methods.cml
index 8dd0e43b71..ca4a0eb017 100644
--- a/lib/iris/tests/results/netcdf/netcdf_cell_methods.cml
+++ b/lib/iris/tests/results/netcdf/netcdf_cell_methods.cml
@@ -1,6 +1,6 @@
 
 
-  
+  
     
       
         
@@ -20,7 +20,7 @@
     
     
   
-  
+  
     
       
         
@@ -41,7 +41,7 @@
     
     
   
-  
+  
     
       
         
@@ -66,7 +66,7 @@
     
     
   
-  
+  
     
       
         
@@ -89,7 +89,7 @@
     
     
   
-  
+  
     
       
         
@@ -112,7 +112,7 @@
     
     
   
-  
+  
     
       
         
@@ -131,7 +131,7 @@
     
     
   
-  
+  
     
       
         
@@ -150,7 +150,7 @@
     
     
   
-  
+  
     
       
         
@@ -170,7 +170,7 @@
     
     
   
-  
+  
     
       
         
@@ -190,7 +190,7 @@
     
     
   
-  
+  
     
       
         
@@ -213,7 +213,7 @@
     
     
   
-  
+  
     
       
         
@@ -232,7 +232,7 @@
     
     
   
-  
+  
     
       
         
@@ -252,7 +252,7 @@
     
     
   
-  
+  
     
       
         
@@ -272,7 +272,7 @@
     
     
   
-  
+  
     
       
         
@@ -295,7 +295,7 @@
     
     
   
-  
+  
     
       
         
@@ -320,7 +320,7 @@
     
     
   
-  
+  
     
       
         
@@ -333,7 +333,7 @@
     
     
   
-  
+  
     
       
         
@@ -346,7 +346,7 @@
     
     
   
-  
+  
     
       
         
@@ -359,7 +359,7 @@
     
     
   
-  
+  
     
       
         
@@ -372,7 +372,7 @@
     
     
   
-  
+  
     
       
         
@@ -385,7 +385,7 @@
     
     
   
-  
+  
     
       
         
@@ -404,7 +404,7 @@
     
     
   
-  
+  
     
       
         
@@ -424,7 +424,7 @@
     
     
   
-  
+  
     
       
         
@@ -447,7 +447,7 @@
     
     
   
-  
+  
     
       
         
@@ -460,7 +460,7 @@
     
     
   
-  
+  
     
       
         
@@ -473,7 +473,7 @@
     
     
   
-  
+  
     
       
         
@@ -486,7 +486,7 @@
     
     
   
-  
+  
     
       
         
@@ -499,7 +499,7 @@
     
     
   
-  
+  
     
       
         
diff --git a/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems.cml b/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems.cml
index 27d4569236..ac41f4a8b8 100644
--- a/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems.cml
+++ b/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems.cml
@@ -13,11 +13,11 @@
         
       
       
-        
+		51, 52, 53, 54, 55, 56, 57, 58, 59, 60]" shape="(60,)" units="Unit('unknown')" value_type="int32" var_name="levelist"/>
       
       
         
@@ -39,11 +39,11 @@
         
       
       
-        
+		51, 52, 53, 54, 55, 56, 57, 58, 59, 60]" shape="(60,)" units="Unit('unknown')" value_type="int32" var_name="levelist"/>
       
       
         
diff --git a/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems_iter_0.cml b/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems_iter_0.cml
index d677191beb..4234b5cc84 100644
--- a/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems_iter_0.cml
+++ b/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems_iter_0.cml
@@ -13,11 +13,11 @@
         
       
       
-        
+		51, 52, 53, 54, 55, 56, 57, 58, 59, 60]" shape="(60,)" units="Unit('unknown')" value_type="int32" var_name="levelist"/>
       
       
         
diff --git a/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems_iter_1.cml b/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems_iter_1.cml
index 775f480c66..17d87a0190 100644
--- a/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems_iter_1.cml
+++ b/lib/iris/tests/results/netcdf/netcdf_global_xyzt_gems_iter_1.cml
@@ -13,11 +13,11 @@
         
       
       
-        
+		51, 52, 53, 54, 55, 56, 57, 58, 59, 60]" shape="(60,)" units="Unit('unknown')" value_type="int32" var_name="levelist"/>
       
       
         
diff --git a/lib/iris/tests/results/netcdf/netcdf_save_no_name.cdl b/lib/iris/tests/results/netcdf/netcdf_save_no_name.cdl
index e67316b2f7..f1399e88b3 100644
--- a/lib/iris/tests/results/netcdf/netcdf_save_no_name.cdl
+++ b/lib/iris/tests/results/netcdf/netcdf_save_no_name.cdl
@@ -6,11 +6,9 @@ variables:
 	double unknown(dim0, dim1) ;
 		unknown:coordinates = "unknown_scalar" ;
 	double dim0(dim0) ;
-		dim0:units = "1" ;
 	double dim1(dim1) ;
 		dim1:units = "m" ;
 	char unknown_scalar(string6) ;
-		unknown_scalar:units = "no_unit" ;
 
 // global attributes:
 		:Conventions = "CF-1.7" ;
diff --git a/lib/iris/tests/results/netcdf/uint32_auxiliary_coord_netcdf3.cml b/lib/iris/tests/results/netcdf/uint32_auxiliary_coord_netcdf3.cml
index 39cb8f2950..e48cf41d2a 100644
--- a/lib/iris/tests/results/netcdf/uint32_auxiliary_coord_netcdf3.cml
+++ b/lib/iris/tests/results/netcdf/uint32_auxiliary_coord_netcdf3.cml
@@ -6,7 +6,7 @@
     
     
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/netcdf/uint32_dimension_coord_netcdf3.cml b/lib/iris/tests/results/netcdf/uint32_dimension_coord_netcdf3.cml
index 1c59fc947e..78fec459e9 100644
--- a/lib/iris/tests/results/netcdf/uint32_dimension_coord_netcdf3.cml
+++ b/lib/iris/tests/results/netcdf/uint32_dimension_coord_netcdf3.cml
@@ -6,7 +6,7 @@
     
     
       
-        
+        
       
     
     
diff --git a/lib/iris/tests/results/unit/fileformats/netcdf/Saver/write/with_climatology.cdl b/lib/iris/tests/results/unit/fileformats/netcdf/Saver/write/with_climatology.cdl
index 3646627746..3c1033c17e 100644
--- a/lib/iris/tests/results/unit/fileformats/netcdf/Saver/write/with_climatology.cdl
+++ b/lib/iris/tests/results/unit/fileformats/netcdf/Saver/write/with_climatology.cdl
@@ -17,10 +17,10 @@ variables:
 	double time_climatology(time, bnds) ;
 	double latitude(latitude) ;
 		latitude:axis = "Y" ;
-		latitude:units = "1" ;
+		latitude:units = "degrees_north" ;
 		latitude:standard_name = "latitude" ;
 	double longitude(longitude) ;
 		longitude:axis = "X" ;
-		longitude:units = "1" ;
+		longitude:units = "degrees_east" ;
 		longitude:standard_name = "longitude" ;
 }
diff --git a/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/002000000000.44.101.131200.1920.09.01.00.00.b_0.cml b/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/002000000000.44.101.131200.1920.09.01.00.00.b_0.cml
index 3ea688d1fa..0bf359e9c4 100644
--- a/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/002000000000.44.101.131200.1920.09.01.00.00.b_0.cml
+++ b/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/002000000000.44.101.131200.1920.09.01.00.00.b_0.cml
@@ -1,6 +1,6 @@
 
 
-  
+  
     
       
       
diff --git a/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/008000000000.44.101.000128.1890.09.01.00.00.b_0.cml b/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/008000000000.44.101.000128.1890.09.01.00.00.b_0.cml
index 829c7ce38e..e5cec55565 100644
--- a/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/008000000000.44.101.000128.1890.09.01.00.00.b_0.cml
+++ b/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/008000000000.44.101.000128.1890.09.01.00.00.b_0.cml
@@ -1,6 +1,6 @@
 
 
-  
+  
     
       
       
diff --git a/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/st0fc699.b_0.cml b/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/st0fc699.b_0.cml
index 4f84609832..b484ebb305 100644
--- a/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/st0fc699.b_0.cml
+++ b/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/st0fc699.b_0.cml
@@ -1,6 +1,6 @@
 
 
-  
+  
     
       
       
diff --git a/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/st0fc942.b_0.cml b/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/st0fc942.b_0.cml
index caafa5845c..c594c748cd 100644
--- a/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/st0fc942.b_0.cml
+++ b/lib/iris/tests/results/usecases/pp_to_cf_conversion/from_netcdf/st0fc942.b_0.cml
@@ -1,6 +1,6 @@
 
 
-  
+  
     
       
       
diff --git a/lib/iris/tests/stock/__init__.py b/lib/iris/tests/stock/__init__.py
index ea6d03a442..42ddfac161 100644
--- a/lib/iris/tests/stock/__init__.py
+++ b/lib/iris/tests/stock/__init__.py
@@ -828,8 +828,8 @@ def jan_offset(day, year):
         units="days since 1970-01-01 00:00:00-00",
         climatological=True,
     )
-    lon_dim = DimCoord(lon, standard_name="longitude")
-    lat_dim = DimCoord(lat, standard_name="latitude")
+    lon_dim = DimCoord(lon, standard_name="longitude", units="degrees")
+    lat_dim = DimCoord(lat, standard_name="latitude", units="degrees")
 
     data_shape = (len(time_points), len(lat), len(lon))
     values = np.zeros(shape=data_shape, dtype=np.int8)
diff --git a/lib/iris/tests/test_aggregate_by.py b/lib/iris/tests/test_aggregate_by.py
index b4e1bad640..bc759f251d 100644
--- a/lib/iris/tests/test_aggregate_by.py
+++ b/lib/iris/tests/test_aggregate_by.py
@@ -89,7 +89,9 @@ def setUp(self):
         )
 
         model_level = iris.coords.DimCoord(
-            np.arange(z_points.size), standard_name="model_level_number"
+            np.arange(z_points.size),
+            standard_name="model_level_number",
+            units="1",
         )
 
         self.cube_single.add_aux_coord(self.coord_z_single, 0)
@@ -124,7 +126,9 @@ def setUp(self):
         )
 
         model_level = iris.coords.DimCoord(
-            np.arange(z1_points.size), standard_name="model_level_number"
+            np.arange(z1_points.size),
+            standard_name="model_level_number",
+            units="1",
         )
 
         self.cube_multi.add_aux_coord(self.coord_z1_multi, 0)
diff --git a/lib/iris/tests/test_concatenate.py b/lib/iris/tests/test_concatenate.py
index bbe5f5eba2..d45a884a2f 100644
--- a/lib/iris/tests/test_concatenate.py
+++ b/lib/iris/tests/test_concatenate.py
@@ -66,52 +66,58 @@ def _make_cube(
     cube_data = np.empty((y_size, x_size), dtype=np.float32)
     cube_data[:] = data
     cube = iris.cube.Cube(cube_data)
-    coord = DimCoord(y_range, long_name="y")
+    coord = DimCoord(y_range, long_name="y", units="1")
     coord.guess_bounds()
     cube.add_dim_coord(coord, 0)
-    coord = DimCoord(x_range, long_name="x")
+    coord = DimCoord(x_range, long_name="x", units="1")
     coord.guess_bounds()
     cube.add_dim_coord(coord, 1)
 
     if aux is not None:
         aux = aux.split(",")
         if "y" in aux:
-            coord = AuxCoord(y_range * 10, long_name="y-aux")
+            coord = AuxCoord(y_range * 10, long_name="y-aux", units="1")
             cube.add_aux_coord(coord, (0,))
         if "x" in aux:
-            coord = AuxCoord(x_range * 10, long_name="x-aux")
+            coord = AuxCoord(x_range * 10, long_name="x-aux", units="1")
             cube.add_aux_coord(coord, (1,))
         if "xy" in aux:
             payload = np.arange(y_size * x_size, dtype=np.float32).reshape(
                 y_size, x_size
             )
-            coord = AuxCoord(payload * 100 + offset, long_name="xy-aux")
+            coord = AuxCoord(
+                payload * 100 + offset, long_name="xy-aux", units="1"
+            )
             cube.add_aux_coord(coord, (0, 1))
 
     if cell_measure is not None:
         cell_measure = cell_measure.split(",")
         if "y" in cell_measure:
-            cm = CellMeasure(y_range * 10, long_name="y-aux")
+            cm = CellMeasure(y_range * 10, long_name="y-aux", units="1")
             cube.add_cell_measure(cm, (0,))
         if "x" in cell_measure:
-            cm = CellMeasure(x_range * 10, long_name="x-aux")
+            cm = CellMeasure(x_range * 10, long_name="x-aux", units="1")
             cube.add_cell_measure(cm, (1,))
         if "xy" in cell_measure:
             payload = x_range + y_range[:, np.newaxis]
-            cm = CellMeasure(payload * 100 + offset, long_name="xy-aux")
+            cm = CellMeasure(
+                payload * 100 + offset, long_name="xy-aux", units="1"
+            )
             cube.add_cell_measure(cm, (0, 1))
 
     if ancil is not None:
         ancil = ancil.split(",")
         if "y" in ancil:
-            av = AncillaryVariable(y_range * 10, long_name="y-aux")
+            av = AncillaryVariable(y_range * 10, long_name="y-aux", units="1")
             cube.add_ancillary_variable(av, (0,))
         if "x" in ancil:
-            av = AncillaryVariable(x_range * 10, long_name="x-aux")
+            av = AncillaryVariable(x_range * 10, long_name="x-aux", units="1")
             cube.add_ancillary_variable(av, (1,))
         if "xy" in ancil:
             payload = x_range + y_range[:, np.newaxis]
-            av = AncillaryVariable(payload * 100 + offset, long_name="xy-aux")
+            av = AncillaryVariable(
+                payload * 100 + offset, long_name="xy-aux", units="1"
+            )
             cube.add_ancillary_variable(av, (0, 1))
 
     if scalar is not None:
@@ -169,50 +175,56 @@ def _make_cube_3d(x, y, z, data, aux=None, offset=0):
     cube_data = np.empty((x_size, y_size, z_size), dtype=np.float32)
     cube_data[:] = data
     cube = iris.cube.Cube(cube_data)
-    coord = DimCoord(z_range, long_name="z")
+    coord = DimCoord(z_range, long_name="z", units="1")
     coord.guess_bounds()
     cube.add_dim_coord(coord, 0)
-    coord = DimCoord(y_range, long_name="y")
+    coord = DimCoord(y_range, long_name="y", units="1")
     coord.guess_bounds()
     cube.add_dim_coord(coord, 1)
-    coord = DimCoord(x_range, long_name="x")
+    coord = DimCoord(x_range, long_name="x", units="1")
     coord.guess_bounds()
     cube.add_dim_coord(coord, 2)
 
     if aux is not None:
         aux = aux.split(",")
         if "z" in aux:
-            coord = AuxCoord(z_range * 10, long_name="z-aux")
+            coord = AuxCoord(z_range * 10, long_name="z-aux", units="1")
             cube.add_aux_coord(coord, (0,))
         if "y" in aux:
-            coord = AuxCoord(y_range * 10, long_name="y-aux")
+            coord = AuxCoord(y_range * 10, long_name="y-aux", units="1")
             cube.add_aux_coord(coord, (1,))
         if "x" in aux:
-            coord = AuxCoord(x_range * 10, long_name="x-aux")
+            coord = AuxCoord(x_range * 10, long_name="x-aux", units="1")
             cube.add_aux_coord(coord, (2,))
         if "xy" in aux:
             payload = np.arange(x_size * y_size, dtype=np.float32).reshape(
                 y_size, x_size
             )
-            coord = AuxCoord(payload + offset, long_name="xy-aux")
+            coord = AuxCoord(payload + offset, long_name="xy-aux", units="1")
             cube.add_aux_coord(coord, (1, 2))
         if "xz" in aux:
             payload = np.arange(x_size * z_size, dtype=np.float32).reshape(
                 z_size, x_size
             )
-            coord = AuxCoord(payload * 10 + offset, long_name="xz-aux")
+            coord = AuxCoord(
+                payload * 10 + offset, long_name="xz-aux", units="1"
+            )
             cube.add_aux_coord(coord, (0, 2))
         if "yz" in aux:
             payload = np.arange(y_size * z_size, dtype=np.float32).reshape(
                 z_size, y_size
             )
-            coord = AuxCoord(payload * 100 + offset, long_name="yz-aux")
+            coord = AuxCoord(
+                payload * 100 + offset, long_name="yz-aux", units="1"
+            )
             cube.add_aux_coord(coord, (0, 1))
         if "xyz" in aux:
             payload = np.arange(
                 x_size * y_size * z_size, dtype=np.float32
             ).reshape(z_size, y_size, x_size)
-            coord = AuxCoord(payload * 1000 + offset, long_name="xyz-aux")
+            coord = AuxCoord(
+                payload * 1000 + offset, long_name="xyz-aux", units="1"
+            )
             cube.add_aux_coord(coord, (0, 1, 2))
 
     return cube
diff --git a/lib/iris/tests/test_coord_api.py b/lib/iris/tests/test_coord_api.py
index bdc6fcc609..c9c686b78b 100644
--- a/lib/iris/tests/test_coord_api.py
+++ b/lib/iris/tests/test_coord_api.py
@@ -247,7 +247,7 @@ def test_basic(self):
             "AuxCoord("
             "array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]),"
             " standard_name=None,"
-            " units=Unit('1'),"
+            " units=Unit('unknown'),"
             " attributes={'monty': 'python'})"
         )
         self.assertEqual(result, str(b))
@@ -337,7 +337,7 @@ def test_basic(self):
             "DimCoord("
             "array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]),"
             " standard_name=None,"
-            " units=Unit('1'),"
+            " units=Unit('unknown'),"
             " attributes={'monty': 'python'})"
         )
         self.assertEqual(result, str(b))
diff --git a/lib/iris/tests/test_iterate.py b/lib/iris/tests/test_iterate.py
index 85f5943b8e..e53eede6f4 100644
--- a/lib/iris/tests/test_iterate.py
+++ b/lib/iris/tests/test_iterate.py
@@ -475,17 +475,24 @@ def test_izip_nd_non_ortho(self):
     def test_izip_nd_ortho(self):
         cube1 = iris.cube.Cube(np.zeros((5, 5, 5, 5, 5), dtype="f8"))
         cube1.add_dim_coord(
-            iris.coords.DimCoord(np.arange(5, dtype="i8"), long_name="z"), [0]
+            iris.coords.DimCoord(
+                np.arange(5, dtype="i8"), long_name="z", units="1"
+            ),
+            [0],
         )
         cube1.add_aux_coord(
             iris.coords.AuxCoord(
-                np.arange(25, dtype="i8").reshape(5, 5), long_name="y"
+                np.arange(25, dtype="i8").reshape(5, 5),
+                long_name="y",
+                units="1",
             ),
             [1, 2],
         )
         cube1.add_aux_coord(
             iris.coords.AuxCoord(
-                np.arange(25, dtype="i8").reshape(5, 5), long_name="x"
+                np.arange(25, dtype="i8").reshape(5, 5),
+                long_name="x",
+                units="1",
             ),
             [3, 4],
         )
diff --git a/lib/iris/tests/test_netcdf.py b/lib/iris/tests/test_netcdf.py
index a550e1ed4b..c69a83edd5 100644
--- a/lib/iris/tests/test_netcdf.py
+++ b/lib/iris/tests/test_netcdf.py
@@ -16,9 +16,11 @@
 import os.path
 import shutil
 import stat
+from subprocess import check_call
 import tempfile
 from unittest import mock
 
+from cf_units import as_unit
 import netCDF4 as nc
 import numpy as np
 import numpy.ma as ma
@@ -27,6 +29,7 @@
 import iris.analysis.trajectory
 import iris.fileformats._pyke_rules.compiled_krb.fc_rules_cf_fc as pyke_rules
 import iris.fileformats.netcdf
+from iris.fileformats.netcdf import load_cubes as nc_load_cubes
 import iris.std_names
 import iris.util
 import iris.coord_systems as icoord_systems
@@ -292,6 +295,42 @@ def test_deferred_loading(self):
             cube[0][(0, 2), (1, 3)], ("netcdf", "netcdf_deferred_mix_1.cml")
         )
 
+    def test_default_units(self):
+        # Note: using a CDL string as a test data reference, rather than a binary file.
+        ref_cdl = """
+            netcdf cm_attr {
+            dimensions:
+                axv = 3 ;
+                ayv = 2 ;
+            variables:
+                int64 qqv(ayv, axv) ;
+                    qqv:long_name = "qq" ;
+                int64 ayv(ayv) ;
+                    ayv:long_name = "y" ;
+                int64 axv(axv) ;
+                    axv:units = "1" ;
+                    axv:long_name = "x" ;
+            data:
+                axv = 11, 12, 13;
+                ayv = 21, 22;
+            }
+            """
+        self.tmpdir = tempfile.mkdtemp()
+        cdl_path = os.path.join(self.tmpdir, "tst.cdl")
+        nc_path = os.path.join(self.tmpdir, "tst.nc")
+        # Write CDL string into a temporary CDL file.
+        with open(cdl_path, "w") as f_out:
+            f_out.write(ref_cdl)
+        # Use ncgen to convert this into an actual (temporary) netCDF file.
+        command = "ncgen -o {} {}".format(nc_path, cdl_path)
+        check_call(command, shell=True)
+        # Load with iris.fileformats.netcdf.load_cubes, and check expected content.
+        cubes = list(nc_load_cubes(nc_path))
+        self.assertEqual(len(cubes), 1)
+        self.assertEqual(cubes[0].units, as_unit("unknown"))
+        self.assertEqual(cubes[0].coord("y").units, as_unit("unknown"))
+        self.assertEqual(cubes[0].coord("x").units, as_unit(1))
+
     def test_units(self):
         # Test exercising graceful cube and coordinate units loading.
         cube0, cube1 = sorted(
@@ -608,10 +647,10 @@ def test_netcdf_multi_with_coords(self):
 
     def test_netcdf_multi_wtih_samedimcoord(self):
         time1 = iris.coords.DimCoord(
-            np.arange(10), standard_name="time", var_name="time"
+            np.arange(10), standard_name="time", var_name="time", units="1"
         )
         time2 = iris.coords.DimCoord(
-            np.arange(20), standard_name="time", var_name="time"
+            np.arange(20), standard_name="time", var_name="time", units="1"
         )
 
         self.cube4.add_dim_coord(time1, 0)
@@ -630,11 +669,13 @@ def test_netcdf_multi_wtih_samedimcoord(self):
     def test_netcdf_multi_conflict_name_dup_coord(self):
         # Duplicate coordinates with modified variable names lookup.
         latitude1 = iris.coords.DimCoord(
-            np.arange(10), standard_name="latitude"
+            np.arange(10), standard_name="latitude", units="1"
+        )
+        time2 = iris.coords.DimCoord(
+            np.arange(2), standard_name="time", units="1"
         )
-        time2 = iris.coords.DimCoord(np.arange(2), standard_name="time")
         latitude2 = iris.coords.DimCoord(
-            np.arange(2), standard_name="latitude"
+            np.arange(2), standard_name="latitude", units="1"
         )
 
         self.cube6.add_dim_coord(latitude1, 0)
@@ -711,10 +752,10 @@ def test_netcdf_save_conflicting_aux(self):
         # Test saving CF-netCDF with multi-dimensional auxiliary coordinates,
         # with conflicts.
         self.cube4.add_aux_coord(
-            iris.coords.AuxCoord(np.arange(10), "time"), 0
+            iris.coords.AuxCoord(np.arange(10), "time", units="1"), 0
         )
         self.cube6.add_aux_coord(
-            iris.coords.AuxCoord(np.arange(10, 20), "time"), 0
+            iris.coords.AuxCoord(np.arange(10, 20), "time", units="1"), 0
         )
 
         cubes = iris.cube.CubeList([self.cube4, self.cube6])
@@ -811,9 +852,11 @@ def test_netcdf_save_conflicting_names(self):
         # Test saving CF-netCDF with a dimension name corresponding to
         # an existing variable name (conflict).
         self.cube4.add_dim_coord(
-            iris.coords.DimCoord(np.arange(10), "time"), 0
+            iris.coords.DimCoord(np.arange(10), "time", units="1"), 0
+        )
+        self.cube6.add_aux_coord(
+            iris.coords.AuxCoord(1, "time", units="1"), None
         )
-        self.cube6.add_aux_coord(iris.coords.AuxCoord(1, "time"), None)
 
         cubes = iris.cube.CubeList([self.cube4, self.cube6])
         with self.temp_filename(suffix=".nc") as file_out:
diff --git a/lib/iris/tests/test_plot.py b/lib/iris/tests/test_plot.py
index 04418d8d40..600801312f 100644
--- a/lib/iris/tests/test_plot.py
+++ b/lib/iris/tests/test_plot.py
@@ -875,6 +875,7 @@ def test_non_cube_coordinate(self):
             pts,
             standard_name="model_level_number",
             attributes={"positive": "up"},
+            units="1",
         )
         self.draw("contourf", cube, coords=["grid_latitude", x])
 
diff --git a/lib/iris/tests/unit/analysis/interpolation/test_RectilinearInterpolator.py b/lib/iris/tests/unit/analysis/interpolation/test_RectilinearInterpolator.py
index d1eaf71030..fa4ca8b608 100644
--- a/lib/iris/tests/unit/analysis/interpolation/test_RectilinearInterpolator.py
+++ b/lib/iris/tests/unit/analysis/interpolation/test_RectilinearInterpolator.py
@@ -34,9 +34,15 @@
 class ThreeDimCube(tests.IrisTest):
     def setUp(self):
         cube = stock.simple_3d_w_multidim_coords()
-        cube.add_aux_coord(iris.coords.DimCoord(np.arange(2), "height"), 0)
-        cube.add_dim_coord(iris.coords.DimCoord(np.arange(3), "latitude"), 1)
-        cube.add_dim_coord(iris.coords.DimCoord(np.arange(4), "longitude"), 2)
+        cube.add_aux_coord(
+            iris.coords.DimCoord(np.arange(2), "height", units="1"), 0
+        )
+        cube.add_dim_coord(
+            iris.coords.DimCoord(np.arange(3), "latitude", units="1"), 1
+        )
+        cube.add_dim_coord(
+            iris.coords.DimCoord(np.arange(4), "longitude", units="1"), 2
+        )
         self.data = np.arange(24).reshape(2, 3, 4).astype(np.float32)
         cube.data = self.data
         self.cube = cube
diff --git a/lib/iris/tests/unit/analysis/regrid/test_RectilinearRegridder.py b/lib/iris/tests/unit/analysis/regrid/test_RectilinearRegridder.py
index 55fb2f4829..492283f843 100644
--- a/lib/iris/tests/unit/analysis/regrid/test_RectilinearRegridder.py
+++ b/lib/iris/tests/unit/analysis/regrid/test_RectilinearRegridder.py
@@ -1253,7 +1253,7 @@ def setUp(self):
             units="m",
             attributes={"positive": "up"},
         )
-        sigma = AuxCoord(1, long_name="sigma")
+        sigma = AuxCoord(1, long_name="sigma", units="1")
         surface_altitude = AuxCoord(
             (src.data - src.data.min()) * 50, "surface_altitude", units="m"
         )
diff --git a/lib/iris/tests/unit/analysis/test_PercentileAggregator.py b/lib/iris/tests/unit/analysis/test_PercentileAggregator.py
index 2b2524795c..cffac86291 100644
--- a/lib/iris/tests/unit/analysis/test_PercentileAggregator.py
+++ b/lib/iris/tests/unit/analysis/test_PercentileAggregator.py
@@ -71,7 +71,7 @@ def test_simple_single_point(self):
         self.assertIs(actual.data, data)
         name = "percentile_over_time"
         coord = actual.coord(name)
-        expected = AuxCoord(percent, long_name=name)
+        expected = AuxCoord(percent, long_name=name, units="percent")
         self.assertEqual(coord, expected)
 
     def test_simple_multiple_points(self):
@@ -89,7 +89,7 @@ def test_simple_multiple_points(self):
         self.assertArrayEqual(actual.data, expected)
         name = "percentile_over_time"
         coord = actual.coord(name)
-        expected = AuxCoord(percent, long_name=name)
+        expected = AuxCoord(percent, long_name=name, units="percent")
         self.assertEqual(coord, expected)
 
     def test_multi_single_point(self):
@@ -105,7 +105,7 @@ def test_multi_single_point(self):
         self.assertIs(actual.data, data)
         name = "percentile_over_time"
         coord = actual.coord(name)
-        expected = AuxCoord(percent, long_name=name)
+        expected = AuxCoord(percent, long_name=name, units="percent")
         self.assertEqual(coord, expected)
 
     def test_multi_multiple_points(self):
@@ -123,7 +123,7 @@ def test_multi_multiple_points(self):
         self.assertArrayEqual(actual.data, expected)
         name = "percentile_over_time"
         coord = actual.coord(name)
-        expected = AuxCoord(percent, long_name=name)
+        expected = AuxCoord(percent, long_name=name, units="percent")
         self.assertEqual(coord, expected)
 
 
diff --git a/lib/iris/tests/unit/analysis/test_WeightedPercentileAggregator.py b/lib/iris/tests/unit/analysis/test_WeightedPercentileAggregator.py
index 1c59ded1fc..878708e48a 100644
--- a/lib/iris/tests/unit/analysis/test_WeightedPercentileAggregator.py
+++ b/lib/iris/tests/unit/analysis/test_WeightedPercentileAggregator.py
@@ -82,7 +82,7 @@ def test_simple_single_point(self):
         self.assertIs(actual.data, data)
         name = "weighted_percentile_over_time"
         coord = actual.coord(name)
-        expected = AuxCoord(percent, long_name=name)
+        expected = AuxCoord(percent, long_name=name, units="percent")
         self.assertEqual(coord, expected)
 
     def test_simple_multiple_points(self):
@@ -107,7 +107,7 @@ def test_simple_multiple_points(self):
         self.assertIs(actual[1], total_weights)
         name = "weighted_percentile_over_time"
         coord = actual[0].coord(name)
-        expected = AuxCoord(percent, long_name=name)
+        expected = AuxCoord(percent, long_name=name, units="percent")
         self.assertEqual(coord, expected)
 
     def test_multi_single_point(self):
@@ -123,7 +123,7 @@ def test_multi_single_point(self):
         self.assertIs(actual.data, data)
         name = "weighted_percentile_over_time"
         coord = actual.coord(name)
-        expected = AuxCoord(percent, long_name=name)
+        expected = AuxCoord(percent, long_name=name, units="percent")
         self.assertEqual(coord, expected)
 
     def test_multi_multiple_points(self):
@@ -141,7 +141,7 @@ def test_multi_multiple_points(self):
         self.assertArrayEqual(actual.data, expected)
         name = "weighted_percentile_over_time"
         coord = actual.coord(name)
-        expected = AuxCoord(percent, long_name=name)
+        expected = AuxCoord(percent, long_name=name, units="percent")
         self.assertEqual(coord, expected)
 
 
diff --git a/lib/iris/tests/unit/aux_factory/test_HybridPressureFactory.py b/lib/iris/tests/unit/aux_factory/test_HybridPressureFactory.py
index 789b1d61d5..14944891f2 100644
--- a/lib/iris/tests/unit/aux_factory/test_HybridPressureFactory.py
+++ b/lib/iris/tests/unit/aux_factory/test_HybridPressureFactory.py
@@ -144,7 +144,9 @@ def setUp(self):
         self.delta = iris.coords.DimCoord(
             [0.0, 1.0, 2.0], long_name="level_pressure", units="Pa"
         )
-        self.sigma = iris.coords.DimCoord([1.0, 0.9, 0.8], long_name="sigma")
+        self.sigma = iris.coords.DimCoord(
+            [1.0, 0.9, 0.8], long_name="sigma", units="1"
+        )
         self.surface_air_pressure = iris.coords.AuxCoord(
             np.arange(4).reshape(2, 2), "surface_air_pressure", units="Pa"
         )
diff --git a/lib/iris/tests/unit/coord_categorisation/test_add_hour.py b/lib/iris/tests/unit/coord_categorisation/test_add_hour.py
index 6965ea7a2f..9b101362a5 100644
--- a/lib/iris/tests/unit/coord_categorisation/test_add_hour.py
+++ b/lib/iris/tests/unit/coord_categorisation/test_add_hour.py
@@ -70,7 +70,7 @@ def test_basic(self):
         cube = self.cube
         time_coord = self.time_coord
         expected_coord = iris.coords.AuxCoord(
-            self.hour_numbers % 24, long_name=coord_name
+            self.hour_numbers % 24, long_name=coord_name, units="1"
         )
 
         ccat.add_hour(cube, time_coord, coord_name)
diff --git a/lib/iris/tests/unit/cube/test_Cube.py b/lib/iris/tests/unit/cube/test_Cube.py
index 9c03f0f4d4..3b98be6454 100644
--- a/lib/iris/tests/unit/cube/test_Cube.py
+++ b/lib/iris/tests/unit/cube/test_Cube.py
@@ -276,7 +276,7 @@ def test_byteorder_true(self):
     def test_cell_measures(self):
         cube = stock.simple_3d_w_multidim_coords()
         cm_a = iris.coords.CellMeasure(
-            np.zeros(cube.shape[-2:]), measure="area"
+            np.zeros(cube.shape[-2:]), measure="area", units="1"
         )
         cube.add_cell_measure(cm_a, (1, 2))
         cm_v = iris.coords.CellMeasure(
@@ -1077,7 +1077,10 @@ def create_cube(lon_min, lon_max, bounds=False):
         0,
     )
     cube.add_aux_coord(
-        iris.coords.AuxCoord([1.0, 0.9, 0.8, 0.6], long_name="sigma"), 0
+        iris.coords.AuxCoord(
+            [1.0, 0.9, 0.8, 0.6], long_name="sigma", units="1"
+        ),
+        0,
     )
     cube.add_dim_coord(
         iris.coords.DimCoord([-45, 0, 45], "latitude", units="degrees"), 1
diff --git a/lib/iris/tests/unit/cube/test_CubeList.py b/lib/iris/tests/unit/cube/test_CubeList.py
index 985c5b6576..2e7b110d60 100644
--- a/lib/iris/tests/unit/cube/test_CubeList.py
+++ b/lib/iris/tests/unit/cube/test_CubeList.py
@@ -151,16 +151,18 @@ class Test_merge__time_triple(tests.IrisTest):
     @staticmethod
     def _make_cube(fp, rt, t, realization=None):
         cube = Cube(np.arange(20).reshape(4, 5))
-        cube.add_dim_coord(DimCoord(np.arange(5), long_name="x"), 1)
-        cube.add_dim_coord(DimCoord(np.arange(4), long_name="y"), 0)
-        cube.add_aux_coord(DimCoord(fp, standard_name="forecast_period"))
+        cube.add_dim_coord(DimCoord(np.arange(5), long_name="x", units="1"), 1)
+        cube.add_dim_coord(DimCoord(np.arange(4), long_name="y", units="1"), 0)
         cube.add_aux_coord(
-            DimCoord(rt, standard_name="forecast_reference_time")
+            DimCoord(fp, standard_name="forecast_period", units="1")
         )
-        cube.add_aux_coord(DimCoord(t, standard_name="time"))
+        cube.add_aux_coord(
+            DimCoord(rt, standard_name="forecast_reference_time", units="1")
+        )
+        cube.add_aux_coord(DimCoord(t, standard_name="time", units="1"))
         if realization is not None:
             cube.add_aux_coord(
-                DimCoord(realization, standard_name="realization")
+                DimCoord(realization, standard_name="realization", units="1")
             )
         return cube
 
diff --git a/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py b/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py
index 08595ed3f3..acea552fdf 100644
--- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py
+++ b/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py
@@ -161,7 +161,7 @@ def _simple_cube(self, dtype):
         points = np.arange(3, dtype=dtype)
         bounds = np.arange(6, dtype=dtype).reshape(3, 2)
         cube = Cube(data, "air_pressure_anomaly")
-        coord = DimCoord(points, bounds=bounds)
+        coord = DimCoord(points, bounds=bounds, units="1")
         cube.add_dim_coord(coord, 0)
         return cube
 
diff --git a/lib/iris/tests/unit/fileformats/pp_load_rules/test__all_other_rules.py b/lib/iris/tests/unit/fileformats/pp_load_rules/test__all_other_rules.py
index d10c1218ab..d44b5a1d54 100644
--- a/lib/iris/tests/unit/fileformats/pp_load_rules/test__all_other_rules.py
+++ b/lib/iris/tests/unit/fileformats/pp_load_rules/test__all_other_rules.py
@@ -269,7 +269,7 @@ def test_month_coord(self):
         res = _all_other_rules(field)[AUX_COORDS_INDEX]
 
         expected = [
-            (AuxCoord(3, long_name="month_number"), None),
+            (AuxCoord(3, long_name="month_number", units="1"), None),
             (AuxCoord("Mar", long_name="month", units=Unit("no unit")), None),
             (
                 DimCoord(
diff --git a/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_scalar_pseudo_level_coords.py b/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_scalar_pseudo_level_coords.py
index 70807408d0..b7074f3c00 100644
--- a/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_scalar_pseudo_level_coords.py
+++ b/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_scalar_pseudo_level_coords.py
@@ -23,7 +23,8 @@ class Test(TestField):
     def test_valid(self):
         coords_and_dims = _convert_scalar_pseudo_level_coords(lbuser5=21)
         self.assertEqual(
-            coords_and_dims, [(DimCoord([21], long_name="pseudo_level"), None)]
+            coords_and_dims,
+            [(DimCoord([21], long_name="pseudo_level", units="1"), None)],
         )
 
     def test_missing_indicator(self):
diff --git a/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_scalar_realization_coords.py b/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_scalar_realization_coords.py
index 4a4649c978..929f65c921 100644
--- a/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_scalar_realization_coords.py
+++ b/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_scalar_realization_coords.py
@@ -24,7 +24,7 @@ def test_valid(self):
         coords_and_dims = _convert_scalar_realization_coords(lbrsvd4=21)
         self.assertEqual(
             coords_and_dims,
-            [(DimCoord([21], standard_name="realization"), None)],
+            [(DimCoord([21], standard_name="realization", units="1"), None)],
         )
 
     def test_missing_indicator(self):
diff --git a/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_vertical_coords.py b/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_vertical_coords.py
index b3a6e537ac..b9a652c397 100644
--- a/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_vertical_coords.py
+++ b/lib/iris/tests/unit/fileformats/pp_load_rules/test__convert_vertical_coords.py
@@ -210,6 +210,7 @@ def _check_depth(
                         lblev,
                         standard_name="model_level_number",
                         attributes={"positive": "down"},
+                        units="1",
                     ),
                     dim,
                 )
@@ -354,6 +355,7 @@ def _check_soil_level(
                 lblev,
                 long_name="soil_model_level_number",
                 attributes={"positive": "down"},
+                units="1",
             )
             expect_result = [(coord, dim)]
         self.assertCoordsAndDimsListsMatch(coords_and_dims, expect_result)
@@ -604,6 +606,7 @@ def _check(
                     lblev,
                     standard_name="model_level_number",
                     attributes={"positive": "up"},
+                    units="1",
                 ),
                 dim,
             )
@@ -630,6 +633,7 @@ def _check(
                     blev,
                     long_name="sigma",
                     bounds=np.vstack((brlev, brsvd1)).T,
+                    units="1",
                 ),
                 dim,
             )
@@ -706,6 +710,7 @@ def _check(
                     lblev,
                     standard_name="model_level_number",
                     attributes={"positive": "up"},
+                    units="1",
                 ),
                 dim,
             )
@@ -732,6 +737,7 @@ def _check(
                     bhlev,
                     long_name="sigma",
                     bounds=np.vstack((bhrlev, brsvd2)).T,
+                    units="1",
                 ),
                 dim,
             )
diff --git a/lib/iris/tests/unit/fileformats/pyke_rules/compiled_krb/fc_rules_cf_fc/test_get_attr_units.py b/lib/iris/tests/unit/fileformats/pyke_rules/compiled_krb/fc_rules_cf_fc/test_get_attr_units.py
index c5e36e8d8e..b752de2370 100644
--- a/lib/iris/tests/unit/fileformats/pyke_rules/compiled_krb/fc_rules_cf_fc/test_get_attr_units.py
+++ b/lib/iris/tests/unit/fileformats/pyke_rules/compiled_krb/fc_rules_cf_fc/test_get_attr_units.py
@@ -44,7 +44,7 @@ def test_unicode_character(self):
         expected_attributes = {'invalid_units': u'\u266b'}
         cf_var = self._make_cf_var()
         attr_units = get_attr_units(cf_var, attributes)
-        self.assertEqual(attr_units, 'unknown')
+        self.assertEqual(attr_units, '?')
         self.assertEqual(attributes, expected_attributes)
 
 
diff --git a/lib/iris/tests/unit/fileformats/um/fast_load/test__convert_collation.py b/lib/iris/tests/unit/fileformats/um/fast_load/test__convert_collation.py
index 3dc6f96d48..7ce0573d25 100644
--- a/lib/iris/tests/unit/fileformats/um/fast_load/test__convert_collation.py
+++ b/lib/iris/tests/unit/fileformats/um/fast_load/test__convert_collation.py
@@ -335,6 +335,7 @@ def test_soil_level(self):
             points,
             long_name="soil_model_level_number",
             attributes={"positive": "down"},
+            units="1",
         )
         coords_and_dims = [(LONGITUDE, 2), (LATITUDE, 1), (level, (0,))]
         self.assertEqual(metadata.dim_coords_and_dims, coords_and_dims)
@@ -416,6 +417,7 @@ def test_vertical_hybrid_height(self):
                     [1, 2, 3],
                     "model_level_number",
                     attributes={"positive": "up"},
+                    units="1",
                 ),
                 (0,),
             ),
@@ -437,6 +439,7 @@ def test_vertical_hybrid_height(self):
                     [0.9994, 0.9979, 0.9957],
                     long_name="sigma",
                     bounds=[[1, 0.9989], [0.9989, 0.9970], [0.9970, 0.9944]],
+                    units="1",
                 ),
                 (0,),
             ),