diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90
index 3bb73e4c57..24e1c3b947 100644
--- a/src/framework/MOM_diag_mediator.F90
+++ b/src/framework/MOM_diag_mediator.F90
@@ -972,7 +972,7 @@ subroutine register_cell_measure(G, diag, Time)
   ! Local variables
   integer :: id
   id = register_diag_field('ocean_model', 'volcello', diag%axesTL, &
-                           Time, 'Ocean grid-cell volume', 'm3', &
+                           Time, 'Ocean grid-cell volume', units='m3', conversion=1.0, &
                            standard_name='ocean_volume', v_extensive=.true., &
                            x_cell_method='sum', y_cell_method='sum')
   call diag_associate_volume_cell_measure(diag, id)
@@ -3153,10 +3153,13 @@ function ocean_register_diag(var_desc, G, diag_CS, day)
   character(len=48) :: units            ! A variable's units.
   character(len=240) :: longname        ! A variable's longname.
   character(len=8) :: hor_grid, z_grid  ! Variable grid info.
+  real :: conversion ! A multiplicative factor for unit conversions for output,
+                     ! as might be needed to convert from intensive to extensive
+                     ! or for dimensional consistency testing [various] or [a A-1 ~> 1]
   type(axes_grp), pointer :: axes => NULL()
 
   call query_vardesc(var_desc, units=units, longname=longname, hor_grid=hor_grid, &
-                     z_grid=z_grid, caller="ocean_register_diag")
+                     z_grid=z_grid, conversion=conversion, caller="ocean_register_diag")
 
   ! Use the hor_grid and z_grid components of vardesc to determine the
   ! desired axes to register the diagnostic field for.
@@ -3211,8 +3214,8 @@ function ocean_register_diag(var_desc, G, diag_CS, day)
         "ocean_register_diag: unknown z_grid component "//trim(z_grid))
   end select
 
-  ocean_register_diag = register_diag_field("ocean_model", trim(var_name), &
-          axes, day, trim(longname), trim(units), missing_value=-1.0e+34)
+  ocean_register_diag = register_diag_field("ocean_model", trim(var_name), axes, day, &
+          trim(longname), units=trim(units), conversion=conversion, missing_value=-1.0e+34)
 
 end function ocean_register_diag
 
diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90
index 8ee192323a..141a57de44 100644
--- a/src/framework/MOM_io.F90
+++ b/src/framework/MOM_io.F90
@@ -555,10 +555,10 @@ subroutine create_MOM_file(IO_handle, filename, vars, novars, fields, &
     pack = 1
     if (present(checksums)) then
       fields(k) = IO_handle%register_field(axes(1:numaxes), vars(k)%name, vars(k)%units, &
-                          vars(k)%longname, pack=pack, checksum=checksums(k,:))
+                          vars(k)%longname, pack=pack, checksum=checksums(k,:), conversion=vars(k)%conversion)
     else
       fields(k) = IO_handle%register_field(axes(1:numaxes), vars(k)%name, vars(k)%units, &
-                          vars(k)%longname, pack=pack)
+                          vars(k)%longname, pack=pack, conversion=vars(k)%conversion)
     endif
   enddo
 
@@ -1880,6 +1880,8 @@ subroutine modify_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, &
   if (present(cmor_longname))   call safe_string_copy(cmor_longname, vd%cmor_longname, &
                                      "vd%cmor_longname of "//trim(vd%name), cllr)
 
+  if (present(conversion)) vd%conversion = conversion
+
   if (present(dim_names)) then
     do n=1,min(5,size(dim_names)) ; if (len_trim(dim_names(n)) > 0) then
       call safe_string_copy(dim_names(n), vd%dim_names(n), "vd%dim_names of "//trim(vd%name), cllr)
@@ -2084,6 +2086,9 @@ subroutine query_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, &
                                      "vd%cmor_units of "//trim(vd%name), cllr)
   if (present(cmor_longname))   call safe_string_copy(vd%cmor_longname, cmor_longname, &
                                      "vd%cmor_longname of "//trim(vd%name), cllr)
+
+  if (present(conversion)) conversion = vd%conversion
+
   if (present(position)) then
     position = vd%position
     if (position == -1) position = position_from_horgrid(vd%hor_grid)
diff --git a/src/framework/MOM_io_file.F90 b/src/framework/MOM_io_file.F90
index 261d4b628d..682f967099 100644
--- a/src/framework/MOM_io_file.F90
+++ b/src/framework/MOM_io_file.F90
@@ -128,6 +128,8 @@ module MOM_io_file
 type :: MOM_field
   character(len=:), allocatable :: label
     !< Identifier for the field in the handle's list
+  real :: conversion
+    !< A factor to use to rescale the field before output [a A-1 ~> 1]
 end type MOM_field
 
 
@@ -454,7 +456,7 @@ end function i_register_axis
 
   !> Interface to register a field to a netCDF file
   function i_register_field(handle, axes, label, units, longname, &
-      pack, standard_name, checksum) result(field)
+      pack, standard_name, checksum, conversion) result(field)
     import :: MOM_file, MOM_axis, MOM_field, int64
     class(MOM_file), intent(inout) :: handle
         !< Handle for a file that is open for writing
@@ -473,6 +475,8 @@ function i_register_field(handle, axes, label, units, longname, &
       !< The standard (e.g., CMOR) name for this variable
     integer(kind=int64), dimension(:), optional, intent(in) :: checksum
       !< Checksum values that can be used to verify reads.
+    real, optional, intent(in) :: conversion
+      !< A factor to use to rescale the field before output [a A-1 ~> 1]
     type(MOM_field) :: field
       !< IO handle for field in MOM_file
   end function i_register_field
@@ -1011,7 +1015,7 @@ end function register_axis_infra
 
 !> Register a field to the MOM framework file
 function register_field_infra(handle, axes, label, units, longname, pack, &
-    standard_name, checksum) result(field)
+    standard_name, checksum, conversion) result(field)
   class(MOM_infra_file), intent(inout) :: handle
     !< Handle for a file that is open for writing
   type(MOM_axis), dimension(:), intent(in) :: axes
@@ -1029,6 +1033,8 @@ function register_field_infra(handle, axes, label, units, longname, pack, &
     !< The standard (e.g., CMOR) name for this variable
   integer(kind=int64), dimension(:), optional, intent(in) :: checksum
     !< Checksum values that can be used to verify reads.
+  real, optional, intent(in) :: conversion
+    !< A factor to use to rescale the field before output [a A-1 ~> 1]
   type(MOM_field) :: field
     !< The field type where this information is stored
 
@@ -1047,6 +1053,7 @@ function register_field_infra(handle, axes, label, units, longname, pack, &
 
   call handle%fields%append(field_infra, label)
   field%label = label
+  field%conversion = 1.0 ; if (present(conversion)) field%conversion = conversion
 end function register_field_infra
 
 
@@ -1069,10 +1076,19 @@ subroutine write_field_4d_infra(handle, field_md, MOM_domain, field, tstamp, &
     !< Missing data fill value
 
   type(fieldtype) :: field_infra
+  real, allocatable :: unscaled_field(:,:,:,:) ! An unscaled version of field for output [a]
 
   field_infra = handle%fields%get(field_md%label)
-  call write_field(handle%handle_infra, field_infra, MOM_domain, field, &
-      tstamp=tstamp, tile_count=tile_count, fill_value=fill_value)
+  if (field_md%conversion == 1.0) then
+    call write_field(handle%handle_infra, field_infra, MOM_domain, field, &
+        tstamp=tstamp, tile_count=tile_count, fill_value=fill_value)
+  else
+    allocate(unscaled_field, source=field)
+    unscaled_field(:,:,:,:) = field_md%conversion * field(:,:,:,:)
+    call write_field(handle%handle_infra, field_infra, MOM_domain, unscaled_field, &
+        tstamp=tstamp, tile_count=tile_count, fill_value=fill_value)
+    deallocate(unscaled_field)
+  endif
 end subroutine write_field_4d_infra
 
 
@@ -1086,7 +1102,7 @@ subroutine write_field_3d_infra(handle, field_md, MOM_domain, field, tstamp, &
   type(MOM_domain_type), intent(in) :: MOM_domain
     !< The MOM_Domain that describes the decomposition
   real, intent(inout) :: field(:,:,:)
-    !< Field to write
+    !< Field to write, perhaps in arbitrary rescaled units [A ~> a]
   real, optional, intent(in) :: tstamp
     !< Model time of this field
   integer, optional, intent(in) :: tile_count
@@ -1095,10 +1111,20 @@ subroutine write_field_3d_infra(handle, field_md, MOM_domain, field, tstamp, &
     !< Missing data fill value
 
   type(fieldtype) :: field_infra
+  real, allocatable :: unscaled_field(:,:,:) ! An unscaled version of field for output [a]
 
   field_infra = handle%fields%get(field_md%label)
-  call write_field(handle%handle_infra, field_infra, MOM_domain, field, &
-      tstamp=tstamp, tile_count=tile_count, fill_value=fill_value)
+  if (field_md%conversion == 1.0) then
+    call write_field(handle%handle_infra, field_infra, MOM_domain, field, &
+        tstamp=tstamp, tile_count=tile_count, fill_value=fill_value)
+  else
+    allocate(unscaled_field, source=field)
+    unscaled_field(:,:,:) = field_md%conversion * field(:,:,:)
+    call write_field(handle%handle_infra, field_infra, MOM_domain, unscaled_field, &
+        tstamp=tstamp, tile_count=tile_count, fill_value=fill_value)
+    deallocate(unscaled_field)
+  endif
+
 end subroutine write_field_3d_infra
 
 
@@ -1121,10 +1147,19 @@ subroutine write_field_2d_infra(handle, field_md, MOM_domain, field, tstamp, &
     !< Missing data fill value
 
   type(fieldtype) :: field_infra
+  real, allocatable :: unscaled_field(:,:) ! An unscaled version of field for output [a]
 
   field_infra = handle%fields%get(field_md%label)
-  call write_field(handle%handle_infra, field_infra, MOM_domain, field, &
-      tstamp=tstamp, tile_count=tile_count, fill_value=fill_value)
+  if (field_md%conversion == 1.0) then
+    call write_field(handle%handle_infra, field_infra, MOM_domain, field, &
+        tstamp=tstamp, tile_count=tile_count, fill_value=fill_value)
+  else
+    allocate(unscaled_field, source=field)
+    unscaled_field(:,:) = field_md%conversion * field(:,:)
+    call write_field(handle%handle_infra, field_infra, MOM_domain, unscaled_field, &
+        tstamp=tstamp, tile_count=tile_count, fill_value=fill_value)
+    deallocate(unscaled_field)
+  endif
 end subroutine write_field_2d_infra
 
 
@@ -1140,9 +1175,17 @@ subroutine write_field_1d_infra(handle, field_md, field, tstamp)
     !< Model time of this field
 
   type(fieldtype) :: field_infra
+  real, allocatable :: unscaled_field(:) ! An unscaled version of field for output [a]
 
   field_infra = handle%fields%get(field_md%label)
-  call write_field(handle%handle_infra, field_infra, field, tstamp=tstamp)
+  if (field_md%conversion == 1.0) then
+    call write_field(handle%handle_infra, field_infra, field, tstamp=tstamp)
+  else
+    allocate(unscaled_field, source=field)
+    unscaled_field(:) = field_md%conversion * field(:)
+    call write_field(handle%handle_infra, field_infra, unscaled_field, tstamp=tstamp)
+    deallocate(unscaled_field)
+  endif
 end subroutine write_field_1d_infra
 
 
@@ -1158,9 +1201,11 @@ subroutine write_field_0d_infra(handle, field_md, field, tstamp)
     !< Model time of this field
 
   type(fieldtype) :: field_infra
+  real :: unscaled_field ! An unscaled version of field for output [a]
 
   field_infra = handle%fields%get(field_md%label)
-  call write_field(handle%handle_infra, field_infra, field, tstamp=tstamp)
+  unscaled_field = field_md%conversion*field
+  call write_field(handle%handle_infra, field_infra, unscaled_field, tstamp=tstamp)
 end subroutine write_field_0d_infra
 
 
@@ -1403,7 +1448,7 @@ end function register_axis_nc
 
 !> Register a field to the MOM netcdf file
 function register_field_nc(handle, axes, label, units, longname, pack, &
-    standard_name, checksum) result(field)
+    standard_name, checksum, conversion) result(field)
   class(MOM_netcdf_file), intent(inout) :: handle
     !< Handle for a file that is open for writing
   type(MOM_axis), intent(in) :: axes(:)
@@ -1421,6 +1466,8 @@ function register_field_nc(handle, axes, label, units, longname, pack, &
     !< The standard (e.g., CMOR) name for this variable
   integer(kind=int64), dimension(:), optional, intent(in) :: checksum
     !< Checksum values that can be used to verify reads.
+  real, optional, intent(in) :: conversion
+    !< A factor to use to rescale the field before output [a A-1 ~> 1]
   type(MOM_field) :: field
 
   type(netcdf_field) :: field_nc
@@ -1438,6 +1485,7 @@ function register_field_nc(handle, axes, label, units, longname, pack, &
     call handle%fields%append(field_nc, label)
   endif
   field%label = label
+  field%conversion = 1.0 ; if (present(conversion)) field%conversion = conversion
 end function register_field_nc
 
 
@@ -1475,11 +1523,19 @@ subroutine write_field_4d_nc(handle, field_md, MOM_domain, field, tstamp, &
     !< Missing data fill value
 
   type(netcdf_field) :: field_nc
+  real, allocatable :: unscaled_field(:,:,:,:) ! An unscaled version of field for output [a]
 
   if (.not. is_root_PE()) return
 
   field_nc = handle%fields%get(field_md%label)
-  call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp)
+  if (field_md%conversion == 1.0) then
+    call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp)
+  else
+    allocate(unscaled_field, source=field)
+    unscaled_field(:,:,:,:) = field_md%conversion * field(:,:,:,:)
+    call write_netcdf_field(handle%handle_nc, field_nc, unscaled_field, time=tstamp)
+    deallocate(unscaled_field)
+  endif
 end subroutine write_field_4d_nc
 
 
@@ -1502,11 +1558,19 @@ subroutine write_field_3d_nc(handle, field_md, MOM_domain, field, tstamp, &
     !< Missing data fill value
 
   type(netcdf_field) :: field_nc
+  real, allocatable :: unscaled_field(:,:,:) ! An unscaled version of field for output [a]
 
   if (.not. is_root_PE()) return
 
   field_nc = handle%fields%get(field_md%label)
-  call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp)
+  if (field_md%conversion == 1.0) then
+    call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp)
+  else
+    allocate(unscaled_field, source=field)
+    unscaled_field(:,:,:) = field_md%conversion * field(:,:,:)
+    call write_netcdf_field(handle%handle_nc, field_nc, unscaled_field, time=tstamp)
+    deallocate(unscaled_field)
+  endif
 end subroutine write_field_3d_nc
 
 
@@ -1529,11 +1593,19 @@ subroutine write_field_2d_nc(handle, field_md, MOM_domain, field, tstamp, &
     !< Missing data fill value
 
   type(netcdf_field) :: field_nc
+  real, allocatable :: unscaled_field(:,:) ! An unscaled version of field for output [a]
 
   if (.not. is_root_PE()) return
 
   field_nc = handle%fields%get(field_md%label)
-  call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp)
+  if (field_md%conversion == 1.0) then
+    call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp)
+  else
+    allocate(unscaled_field, source=field)
+    unscaled_field(:,:) = field_md%conversion * field(:,:)
+    call write_netcdf_field(handle%handle_nc, field_nc, unscaled_field, time=tstamp)
+    deallocate(unscaled_field)
+  endif
 end subroutine write_field_2d_nc
 
 
@@ -1549,11 +1621,19 @@ subroutine write_field_1d_nc(handle, field_md, field, tstamp)
     !< Model time of this field
 
   type(netcdf_field) :: field_nc
+  real, allocatable :: unscaled_field(:) ! An unscaled version of field for output [a]
 
   if (.not. is_root_PE()) return
 
   field_nc = handle%fields%get(field_md%label)
-  call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp)
+  if (field_md%conversion == 1.0) then
+    call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp)
+  else
+    allocate(unscaled_field, source=field)
+    unscaled_field(:) = field_md%conversion * field(:)
+    call write_netcdf_field(handle%handle_nc, field_nc, unscaled_field, time=tstamp)
+    deallocate(unscaled_field)
+  endif
 end subroutine write_field_1d_nc
 
 
@@ -1569,11 +1649,13 @@ subroutine write_field_0d_nc(handle, field_md, field, tstamp)
     !< Model time of this field
 
   type(netcdf_field) :: field_nc
+  real :: unscaled_field ! An unscaled version of field for output [a]
 
   if (.not. is_root_PE()) return
 
   field_nc = handle%fields%get(field_md%label)
-  call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp)
+  unscaled_field = field_md%conversion * field
+  call write_netcdf_field(handle%handle_nc, field_nc, unscaled_field, time=tstamp)
 end subroutine write_field_0d_nc