diff --git a/demo/app/plot-raw-histograms.gnu b/demo/app/histogram-plot.gnu
similarity index 95%
rename from demo/app/plot-raw-histograms.gnu
rename to demo/app/histogram-plot.gnu
index fe081ff59..00808f81f 100644
--- a/demo/app/plot-raw-histograms.gnu
+++ b/demo/app/histogram-plot.gnu
@@ -4,11 +4,11 @@ set term png size 2600, 1200
 set output "input-tensor-histograms.png"
 set multiplot layout 3,3
 do for [name in "potential_temperature pressure qc qr qs qv temperature"] {
-  filename = name . ".plt"
+  filename = name . "-inputs.plt"
   set title sprintf("Histogram for %s",name)
   plot filename using 1:2 with linespoints title columnheader
 }
-print "\n Output-tensor histograms have been written to the graphics file 'output-tensor-histograms.png'."
+print " Input-tensor histograms have been written to the graphics file 'input-tensor-histograms.png'.\n"
 
 unset multiplot
 
@@ -19,4 +19,4 @@ do for [name in "dpotential_temperature_dt dqc_dt dqr_dt dqs_dt dqv_dt"] {
   set title sprintf("Histogram for %s",name)
   plot filename using 1:2 with linespoints title columnheader
 }
-print " Input-tensor histograms have been written to the graphics file 'input-tensor-histograms.png'.\n"
+print "\n Output-tensor histograms have been written to the graphics file 'output-tensor-histograms.png'."
diff --git a/demo/app/tensor-statistics.F90 b/demo/app/tensor-statistics.F90
index 4819e23a3..dcec19177 100644
--- a/demo/app/tensor-statistics.F90
+++ b/demo/app/tensor-statistics.F90
@@ -11,235 +11,195 @@ program tensor_statistics
 
   ! External dependencies:
   use julienne_m, only : command_line_t, file_t, string_t, call_julienne_assert_
-  use ieee_arithmetic, only : ieee_is_nan
   use iso_fortran_env, only : int64, real64
     
   ! Internal dependencies:
   use NetCDF_file_m, only: NetCDF_file_t
-  use NetCDF_variable_m, only: NetCDF_variable_t, time_derivative_t
-  use histogram_m, only: histogram_t, to_file
+  use NetCDF_variable_m, only: NetCDF_variable_t, time_derivative_t, histogram
+  use histogram_m, only: to_file
   use time_data_m, only: time_data_t, icar_output_file_t
   use fiats_m, only : training_configuration_t, training_data_files_t
   implicit none
 
-  character(len=*), parameter :: usage =                new_line('') // new_line('') // & 
-    'Usage: '                                        // new_line('') // new_line('') // &
-    './build/run-fpm.sh run tensor-statistics -- --bins <integer> \' // new_line('') // &
-    '  [--raw] [--start <integer>] [--end <integer>] [--stride <integer>]'           // &
-                                                        new_line('') // new_line('') // &
-    'where angular brackets denote user-provided values and square brackets denote optional arguments.' // new_line('')
-
+  character(len=*), parameter :: usage = new_line('') &
+    // new_line('') // 'Usage: ' &
+    // new_line('') &
+    // new_line('') // './build/run-fpm.sh run tensor-statistics -- --bins <integer>' &
+    // new_line('') &
+    // new_line('') // 'where angular brackets (<>) denote a user-provided value.' &
+    // new_line('') &
+    // new_line('') // '--bins  number of histogram bins in which to collect variable values' &
+    // new_line('')
+
+  type command_line_arguments_t
+    integer :: num_bins = 0
+  end type
   integer(int64) t_start, t_finish, clock_rate
-  integer num_bins, start_step, stride
-  integer, allocatable :: end_step
-  logical raw
-  type(string_t), allocatable :: input_file_names(:), output_file_names(:)
-
-  call system_clock(t_start, clock_rate)
-  call get_command_line_arguments(num_bins, start_step, end_step, stride, raw)
-
-  associate( &
-     training_configuration => training_configuration_t(file_t("training_configuration.json")) &
-    ,training_data_files => training_data_files_t(file_t("training_data_files.json")) &
-  )
-
-    call compute_histograms( &
-       input_tensor_file_names  = training_data_files%fully_qualified_inputs_files() &
-      ,output_tensor_file_names = training_data_files%fully_qualified_outputs_files() &
-      ,time_data_file_name      = training_data_files%fully_qualified_time_file() &
-      ,input_component_names    = training_configuration%input_variable_names() &
-      ,output_component_names   = training_configuration%output_variable_names() &
-      ,raw = raw &
-    )
-  end associate
-  call system_clock(t_finish)
-
-  print '(a,g0)',"System clock time: ", real(t_finish - t_start, real64)/real(clock_rate, real64)
-  print *
-  print '(a)',"______________________________________________________"
-  print '(a)',"The *.plt files contain tensor ranges and histograms."
-  print '(a)',"If you have gnuplot installed, please execute the"
-  print '(a)',"following command to produce histogram visualizations:"
-  print *
-  associate(file_name => trim(merge("plot-raw-histograms       ", "plot-normalized-histograms", raw)) // ".gnu")
-   print*,"  gnuplot app/" // file_name
+
+  associate(args => command_line_arguments())
+    if (args%num_bins < 1) error stop usage
+    call system_clock(t_start, clock_rate)
+    call calc_write_stats(args, training_data_files_t(file_t("training_data_files.json")), training_configuration_t(file_t("training_configuration.json")))
+    call system_clock(t_finish)
+    call print_trailer(args, elapsed_time = real(t_finish - t_start, real64)/real(clock_rate, real64))
   end associate
-  print *
-  print '(a)',"_______________ tensor_statistics done________________"
 
 contains
 
-  subroutine get_command_line_arguments(num_bins, start_step, end_step, stride, raw)
-    integer, intent(out) :: num_bins, start_step, stride
-    integer, intent(out), allocatable :: end_step
-    logical, intent(out) :: raw
-
-    ! local variables
+  type(command_line_arguments_t) function command_line_arguments()
     type(command_line_t) command_line
-    character(len=:), allocatable :: stride_string, bins_string, start_string, end_string
-
+    character(len=:), allocatable :: bins_string
     bins_string = command_line%flag_value("--bins")
-    start_string = command_line%flag_value("--start")
-    end_string = command_line%flag_value("--end")
-    stride_string = command_line%flag_value("--stride")
-    raw = command_line%argument_present(["--raw"])
-
-    associate(required_arguments => len(bins_string)/=0)
-       if (.not. required_arguments) error stop usage 
-    end associate
-
-    read(bins_string,*) num_bins
-
-    if (len(stride_string)==0) then
-      stride = 1
-    else
-      read(stride_string,*) stride
-    end if
-
-    if (len(start_string)==0) then
-      start_step = 1
-    else
-      read(start_string,*) start_step
-    end if
-
-    if (len(end_string)/=0) then
-      if (.not. allocated(end_step)) allocate(end_step)
-      read(end_string,*) end_step
-    end if
- 
-  end subroutine get_command_line_arguments
-
-  subroutine compute_histograms( &
-    input_tensor_file_names, output_tensor_file_names, input_component_names, output_component_names, time_data_file_name, raw &
-  )
-    type(string_t), intent(in) :: input_tensor_file_names(:), output_tensor_file_names(:)
-    type(string_t), intent(in) :: input_component_names(:),   output_component_names(:), time_data_file_name
-    logical, intent(in) :: raw
-
-    type(time_derivative_t), allocatable :: derivative(:)
-    type(NetCDF_variable_t), allocatable :: input_variable(:), output_variable(:)
-    type(NetCDF_variable_t) input_time, output_time
-    double precision, allocatable, dimension(:) :: dt
-    double precision, parameter :: tolerance = 1.E-07
-    integer(int64) t_histo_start, t_histo_finish 
-    integer t, t_end, v
-
-    allocate(input_variable(size(input_component_names)))
-
-    print '(a)',"Reading physics-based model inputs from " // input_tensor_file_names(1)%string()
-
-    input_file: &
-    associate(NetCDF_file => netCDF_file_t(input_tensor_file_names(1)))
-
-      do v=1, size(input_variable) 
-        print '(a)',"- reading " // input_component_names(v)%string()
-        call input_variable(v)%input(input_component_names(v), NetCDF_file, rank=4)
-      end do
-
-      do v = 2, size(input_variable)
-        call_julienne_assert(input_variable(v)%conformable_with(input_variable(1)))
-      end do
-
-      print '(a)',"- reading time from NetCDF file"
-      call input_time%input("time", NetCDF_file, rank=1)
-
-    end associate input_file
-
-    print '(a)',"Calculating input tensor histograms."
-    call system_clock(t_histo_start)
-    associate(histograms => [(input_variable(v)%histogram(num_bins, raw), v = 1, size(input_variable))])
-      call system_clock(t_histo_finish)
-      print '(i0,a,g0,a)',size(histograms), " input histograms done in ", real(t_histo_finish - t_histo_start, real64)/real(clock_rate, real64), " sec."
-
-      print '(a)',"Writing input tensor histograms file"
-      block
-        type(file_t) histograms_file
-        integer h
-
-        if (raw) then
-          do h = 1, size(histograms)
-            histograms_file = to_file(histograms(h))
-            associate(gnuplot_file_name => histograms(h)%variable_name() // ".plt")
-              print '(a)',"- writing " // gnuplot_file_name
-              call histograms_file%write_lines(string_t(gnuplot_file_name))
-            end associate
-          end do
-        else
-          histograms_file = to_file(histograms)
-          associate(gnuplot_file_name => "inputs_stats.plt")
-            print '(a)',"- writing " // gnuplot_file_name
-            call histograms_file%write_lines(string_t(gnuplot_file_name))
-          end associate
-        end if
-      end block
-    end associate
-
-    allocate(output_variable(size(output_component_names)))
-
-    print '(a)',"Reading physics-based model outputs from " // output_tensor_file_names(1)%string()
-
-    output_file: &
-    associate(NetCDF_file => netCDF_file_t(output_tensor_file_names(1)))
-
-      do v=1, size(output_variable)
-        print '(a)', "- reading " // output_component_names(v)%string()
-        call output_variable(v)%input(output_component_names(v), NetCDF_file, rank=4)
-      end do
-
-      do v = 2, size(output_variable)
-        call_julienne_assert(output_variable(v)%conformable_with(output_variable(1)))
-      end do
-
-      print '(a)',"- reading time"
-      call output_time%input("time", NetCDF_file, rank=1)
-
-    end associate output_file
-
-    print '(a)',"Calculating the desired neural-network model outputs: time derivatives of a subset of the inputs"
-
-    allocate(derivative(size(output_variable)))
-
-    print '(a)',"- reading time from JSON file"
-    associate(time_data => time_data_t(file_t(time_data_file_name)))
-      do v = 1, size(derivative)
-        derivative_name: &
-        associate(derivative_name => "d" // output_component_names(v)%string() // "_dt")
-          print '(a)',"- " // derivative_name
-          derivative(v) = time_derivative_t(old = input_variable(v), new = output_variable(v), dt=time_data%dt())
-          call_julienne_assert(.not. derivative(v)%any_nan())
-        end associate derivative_name
-      end do
-    end associate
+    if (len(bins_string) /= 0) read(bins_string,*) command_line_arguments%num_bins
+  end function
+
+  subroutine print_trailer(args, elapsed_time)
+    type(command_line_arguments_t), intent(in) :: args
+    real(real64), intent(in) :: elapsed_time
+
+    print '(a,g0)',"System clock time: ", elapsed_time
+    print *
+    print '(a)',"______________________________________________________"
+    print '(a)',"The *.plt files contain tensor ranges and histograms."
+    print '(a)',"With gnuplot installed, please execute the following"
+    print '(a)',"following command to produce histogram visualizations:"
+    print *
+    print '(a)',"  gnuplot app/histogram-plot.gnu"
+    print *
+    print '(a)',"_______________ tensor_statistics done________________"
+  end subroutine
 
-    print '(a)',"Calculating histograms for the desired neural-network outputs."
-    call system_clock(t_histo_start)
-    associate(histograms => [(derivative(v)%histogram(num_bins, raw), v = 1, size(derivative))])
-      call system_clock(t_histo_finish)
-      print '(i0,a,g0,a)',size(histograms), " output histograms done in ", real(t_histo_finish - t_histo_start, real64)/real(clock_rate, real64), " sec."
-
-      print '(a)',"Writing desired-output tensor histograms file"
-      block
-        type(file_t) histograms_file
-        integer h
-
-        if (raw) then
-          do h = 1, size(histograms)
-            histograms_file = to_file(histograms(h))
-            associate(gnuplot_file_name => histograms(h)%variable_name() // ".plt")
-              print '(a)',"- writing " // gnuplot_file_name
-              call histograms_file%write_lines(string_t(gnuplot_file_name))
-            end associate
-          end do
-        else
-          histograms_file = to_file(histograms)
-          associate(gnuplot_file_name => "outputs_stats.plt")
-            print '(a)',"- writing " // gnuplot_file_name
-            call histograms_file%write_lines(string_t(gnuplot_file_name))
-          end associate
-        end if
-      end block
+  subroutine calc_write_stats(args, training_data_files, training_configuration)
+    type(command_line_arguments_t), intent(in) :: args
+    type(training_data_files_t)   , intent(in) :: training_data_files
+    type(training_configuration_t), intent(in) :: training_configuration
+
+    type(NetCDF_file_t)    , allocatable, dimension(:)   :: NetCDF_input_file, NetCDF_output_file
+    type(time_derivative_t), allocatable, dimension(:,:) :: derivative
+    type(NetCDF_variable_t), allocatable, dimension(:,:) :: input_variable, output_variable
+    integer(int64) t_histo_start, t_histo_finish
+    integer v, f
+
+    associate( &
+       input_tensor_file_names  => training_data_files%fully_qualified_inputs_files() &
+      ,output_tensor_file_names => training_data_files%fully_qualified_outputs_files() &
+      ,time_data_file_name      => training_data_files%fully_qualified_time_file() &
+      ,input_component_names    => training_configuration%input_variable_names() &
+      ,output_component_names   => training_configuration%output_variable_names() &
+    )
+      allocate(NetCDF_input_file(size(input_tensor_file_names)))
+      allocate(input_variable(size(input_component_names), size(NetCDF_input_file)))
+
+      input_variable_histograms: &
+      associate(num_input_files => size(NetCDF_input_file), num_variables => size(input_variable,1))
+
+        read_files: &
+        do f = 1, num_input_files
+
+          print '(a)',"Reading physics-based model inputs from " // input_tensor_file_names(f)%string()
+          NetCDF_input_file(f) = netCDF_file_t(input_tensor_file_names(f))
+
+          read_variables: &
+          do v = 1, num_variables
+            print '(a)',"- reading " // input_component_names(v)%string() // " from " // input_tensor_file_names(f)%string()
+            call input_variable(v,f)%input(input_component_names(v), NetCDF_input_file(f), rank=4)
+            call_julienne_assert(input_variable(v,f)%conformable_with(input_variable(1,f)))
+          end do read_variables
+
+        end do read_files
+
+        print '(a)',"Calculating input tensor histograms."
+        call system_clock(t_histo_start, clock_rate)
+
+        compute_histograms: &
+        associate(histograms => [(histogram(input_variable(v,:), args%num_bins), v = 1, num_variables)])
+          call system_clock(t_histo_finish)
+          print '(a,i0,a,g0,a)',"computed ", size(histograms), " histograms in " &
+            ,real(t_histo_finish - t_histo_start, real64)/real(clock_rate, real64), " sec."
+
+          print '(a)',"Writing input tensor histogram file(s)"
+
+          write_histograms: &
+          associate(input_file_names => training_data_files%inputs_files())
+            do v = 1, num_variables
+              associate(histograms_file => to_file(histograms(v)))
+                associate(gnuplot_file_name => histograms(v)%variable_name() // "-" // "inputs.plt")
+                  print '(a)',"- writing " // gnuplot_file_name
+                  call histograms_file%write_lines(gnuplot_file_name)
+                end associate
+              end associate
+            end do
+          end associate write_histograms
+        end associate compute_histograms
+      end associate input_variable_histograms
+
+      ! Currently unused because time stamps were not recorded correctly in the training-data files
+      ! type(NetCDF_variable_t) input_time, output_time
+      ! print '(a)',"- reading time from NetCDF file"
+      ! call input_time%input("time", NetCDF_input_file(1), rank=1)
+      ! print '(a)',"- reading time"
+      ! call output_time%input("time", NetCDF_output_file(1), rank=1)
+
+      allocate(NetCDF_output_file(size(output_tensor_file_names)))
+      allocate(output_variable(size(output_component_names), size(NetCDF_output_file)))
+
+      derivative_histograms: &
+      associate(num_output_files => size(NetCDF_output_file), num_output_variables => size(output_variable,1))
+
+        print '(a)',"- reading time from JSON file"
+        associate(time_data => time_data_t(file_t(time_data_file_name)))
+
+          print '(a)',"Calculating the desired neural-network model outputs: time derivatives of the outputs"
+          allocate(derivative(num_output_variables, num_output_files))
+
+          read_files: &
+          do f = 1, num_output_files
+
+            print '(a)',"Reading physics-based model outputs from " // output_tensor_file_names(f)%string()
+            NetCDF_output_file(f) = netCDF_file_t(output_tensor_file_names(f))
+
+            read_variables: &
+            do v = 1, num_output_variables
+
+              print '(a)',"- reading " // output_component_names(v)%string() // " from " // output_tensor_file_names(f)%string()
+              call output_variable(v,f)%input(output_component_names(v), NetCDF_output_file(f), rank=4)
+              call_julienne_assert(output_variable(v,f)%conformable_with(output_variable(1,f)))
+
+              derivative_name: &
+              associate(derivative_name => "d" // output_component_names(v)%string() // "_dt")
+                print '(a)',"- calculating " // derivative_name
+                derivative(v,f) = time_derivative_t(old = input_variable(v,1), new = output_variable(v,1), dt=time_data%dt())
+                call_julienne_assert(.not. derivative(v,f)%any_nan())
+              end associate derivative_name
+            end do read_variables
+          end do read_files
+        end associate
+
+        print '(a)',"Calculating histograms for the desired neural-network outputs."
+        call system_clock(t_histo_start, clock_rate)
+
+        compute_histograms: &
+        associate(histograms => [(histogram(derivative(v,:), args%num_bins), v = 1, num_output_variables)])
+          call system_clock(t_histo_finish)
+          print '(a,i0,a,g0,a)',"computed ", size(histograms), " histograms in " &
+            ,real(t_histo_finish - t_histo_start, real64)/real(clock_rate, real64), " sec."
+
+          print '(a)',"Writing desired-output tensor histograms file"
+          block
+            type(file_t) histograms_file
+            integer h
+
+            do h = 1, size(histograms)
+              histograms_file = to_file(histograms(h))
+              associate(gnuplot_file_name => histograms(h)%variable_name() // ".plt")
+                print '(a)',"- writing " // gnuplot_file_name
+                call histograms_file%write_lines(string_t(gnuplot_file_name))
+              end associate
+            end do
+          end block
+        end associate compute_histograms
+      end associate derivative_histograms
     end associate
-
   end subroutine
 
 end program tensor_statistics
diff --git a/demo/app/plot-normalized-histograms.gnu b/demo/app/unified-histogram-plot.gnu
similarity index 100%
rename from demo/app/plot-normalized-histograms.gnu
rename to demo/app/unified-histogram-plot.gnu
diff --git a/demo/src/NetCDF_variable_m.f90 b/demo/src/NetCDF_variable_m.f90
index ce455a42d..f3cf4d022 100644
--- a/demo/src/NetCDF_variable_m.f90
+++ b/demo/src/NetCDF_variable_m.f90
@@ -12,6 +12,7 @@ module NetCDF_variable_m
   public :: NetCDF_variable_t
   public :: tensors
   public :: time_derivative_t
+  public :: histogram
 
   type NetCDF_variable_t(k)
     integer, kind :: k = default_real
@@ -35,8 +36,6 @@ module NetCDF_variable_m
     procedure, private, non_overridable :: default_real_maximum         , double_precision_maximum
     generic :: operator(-)              => default_real_subtract        , double_precision_subtract
     procedure, private, non_overridable :: default_real_subtract        , double_precision_subtract
-    generic :: histogram                => default_real_histogram       , double_precision_histogram
-    procedure, private, non_overridable :: default_real_histogram       , double_precision_histogram
     generic :: compare                  => default_real_compare         , double_precision_compare
     procedure, private, non_overridable :: default_real_compare         , double_precision_compare
   end type
@@ -148,13 +147,13 @@ elemental module function double_precision_rank(self) result(my_rank)
 
     elemental module function default_real_end_step(self) result(end_step)
       implicit none
-      class(NetCDF_variable_t), intent(inout) :: self
+      class(NetCDF_variable_t), intent(in) :: self
       integer end_step
     end function
 
     elemental module function double_precision_end_step(self) result(end_step)
       implicit none
-      class(NetCDF_variable_t(double_precision)), intent(inout) :: self
+      class(NetCDF_variable_t(double_precision)), intent(in) :: self
       integer end_step
     end function
 
@@ -215,32 +214,38 @@ module function tensors(NetCDF_variables, step_start, step_end, step_stride)
 
     elemental module function default_real_end_time(self) result(end_time)
       implicit none
-      class(NetCDF_variable_t), intent(inout) :: self
+      class(NetCDF_variable_t), intent(in) :: self
       integer end_time
     end function
 
     elemental module function double_precision_end_time(self) result(end_time)
       implicit none
-      class(NetCDF_variable_t), intent(inout) :: self
+      class(NetCDF_variable_t), intent(in) :: self
       integer end_time
     end function
 
-    elemental module function default_real_histogram(self, num_bins, raw) result(histogram)
+  end interface
+
+  interface histogram
+
+    module function default_real_histogram(variables, num_bins) result(aggregate)
       implicit none
-      class(NetCDF_variable_t), intent(inout) :: self
+      class(NetCDF_variable_t), intent(in), target :: variables(:)
       integer, intent(in) :: num_bins
-      logical, intent(in) :: raw
-      type(histogram_t) histogram
+      type(histogram_t) aggregate
     end function
 
-    elemental module function double_precision_histogram(self, num_bins, raw) result(histogram)
+    module function double_precision_histogram(variables, num_bins) result(aggregate)
       implicit none
-      class(NetCDF_variable_t(double_precision)), intent(inout) :: self
+      class(NetCDF_variable_t(double_precision)), intent(in) :: variables
       integer, intent(in) :: num_bins
-      logical, intent(in) :: raw
-      type(histogram_t) histogram
+      type(histogram_t) aggregate
     end function
 
+  end interface
+
+  interface
+
     module subroutine default_real_compare(self, values)
       implicit none
       class(NetCDF_variable_t), intent(in) :: self
diff --git a/demo/src/NetCDF_variable_s.F90 b/demo/src/NetCDF_variable_s.F90
index 21ed17627..760286a8b 100644
--- a/demo/src/NetCDF_variable_s.F90
+++ b/demo/src/NetCDF_variable_s.F90
@@ -6,7 +6,10 @@
 submodule(NetCDF_variable_m) NetCDF_variable_s
   use ieee_arithmetic, only : ieee_is_nan
   use kind_parameters_m, only : default_real
-  use julienne_m, only : call_julienne_assert_, operator(.equalsExpected.)
+  use julienne_m, only : &
+     call_julienne_assert_ &
+    ,operator(.all.) &
+    ,operator(.equalsExpected.)
   use default_m, only : default_or_present_value
   implicit none
 
@@ -66,16 +69,22 @@
    end procedure
 
   module procedure default_real_histogram
-    select case(self%rank())
+
+    integer v
+
+    call_julienne_assert(.all. (variables(1)%rank() .equalsExpected.   variables%rank()))
+    call_julienne_assert(.all. (variables(1)%name_  .equalsExpected. [(variables(v)%name_, v = 1, size(variables))]))
+
+    select case(variables(1)%rank())
     case (4)
-      histogram = histogram_t(self%values_4D_ , self%name_, num_bins, raw)
+      aggregate = histogram_t([(variables(v)%values_4D_, v = 1, size(variables))], variables(1)%name_, num_bins)
     case default
       error stop 'NetCDF_variable_s(default_real_histogram): unsupported rank'
     end select 
   end procedure 
 
   module procedure double_precision_histogram
-      error stop 'NetCDF_variable_s(double_precision_histogram): unsupported rank'
+      error stop 'NetCDF_variable_s(double_precision_histogram): unsupported kind'
   end procedure 
 
   module procedure default_real_copy
diff --git a/demo/src/histogram_m.f90 b/demo/src/histogram_m.f90
index 21aaba85f..b4e779474 100644
--- a/demo/src/histogram_m.f90
+++ b/demo/src/histogram_m.f90
@@ -12,26 +12,30 @@ module histogram_m
     !! encapsulate the primary data associated with histograms
     private
     character(len=:), allocatable :: variable_name_
-    real unmapped_min_, unmapped_max_
-    real, allocatable :: bin_frequency_(:), bin_value_(:)
+    real variable_min_, variable_max_
+    real, allocatable :: bin_value_(:)
+    integer, allocatable :: bin_count_(:)
   contains
     procedure variable_name
-    procedure unmapped_range
-    procedure unmapped_min
-    procedure unmapped_max
-    procedure num_bins
-    procedure bin_value
     procedure bin_frequency
+    procedure bin_count
   end type
 
   interface histogram_t
 
-    pure module function construct(v, variable_name, num_bins, raw) result(histogram)
+    module function construct(v, variable_name, num_bins) result(histogram)
       implicit none
-      real, intent(in) :: v(:,:,:,:)
+      real, intent(in) :: v(:)
+      character(len=*), intent(in) :: variable_name
+      integer, intent(in) :: num_bins
+      type(histogram_t) histogram
+    end function
+
+    pure module function construct_in_range(v, variable_name, num_bins, v_min, v_max) result(histogram)
+      implicit none
+      real, intent(in) :: v(:), v_min, v_max
       character(len=*), intent(in) :: variable_name
       integer, intent(in) :: num_bins
-      logical, intent(in) :: raw
       type(histogram_t) histogram
     end function
 
@@ -55,49 +59,23 @@ pure module function to_common_file(histograms) result(file)
 
   interface
 
-    pure module function num_bins(self) result(bins)
-      implicit none
-      class(histogram_t), intent(in) :: self
-      integer bins
-    end function
-
     pure module function variable_name(self) result(name)
       implicit none
       class(histogram_t), intent(in) :: self
       character(len=:), allocatable :: name
     end function
 
-    pure module function unmapped_range(self) result(raw_range)
-      implicit none
-      class(histogram_t), intent(in) :: self
-      integer, parameter :: num_end_points = 2
-      real raw_range(num_end_points)
-    end function
-
-    elemental module function unmapped_min(self) result(range_minimum)
-      implicit none
-      class(histogram_t), intent(in) :: self
-      real range_minimum
-    end function
-
-    elemental module function unmapped_max(self) result(range_maximum)
-      implicit none
-      class(histogram_t), intent(in) :: self
-      real range_maximum
-    end function
-
-    elemental module function bin_value(self, bin) result(v)
+    pure module function bin_frequency(self, bin) result(frequency)
       implicit none
       class(histogram_t), intent(in) :: self
       integer, intent(in) :: bin
-      real v
+      real frequency
     end function
 
-    pure module function bin_frequency(self, bin) result(frequency)
+    pure module function bin_count(self) result(counts)
       implicit none
       class(histogram_t), intent(in) :: self
-      integer, intent(in) :: bin
-      real frequency
+      integer, allocatable :: counts(:)
     end function
 
   end interface
diff --git a/demo/src/histogram_s.F90 b/demo/src/histogram_s.F90
index 789bcdff0..c779c3046 100644
--- a/demo/src/histogram_s.F90
+++ b/demo/src/histogram_s.F90
@@ -19,28 +19,12 @@
     name = self%variable_name_
   end procedure
 
-  module procedure unmapped_range
-    raw_range = [self%unmapped_min_, self%unmapped_max_]
-  end procedure
-
-  module procedure unmapped_max
-    range_maximum = self%unmapped_max_
-  end procedure
-
-  module procedure unmapped_min
-    range_minimum = self%unmapped_min_
-  end procedure
-
-  module procedure num_bins
-    bins = size(self%bin_value_)
-  end procedure
-
-  module procedure bin_value
-    v = self%bin_value_(bin)
+  module procedure bin_count
+    counts = self%bin_count_
   end procedure
 
   module procedure bin_frequency
-    frequency = self%bin_frequency_(bin)
+    frequency = self%bin_count_(bin)/real(sum(self%bin_count_))
   end procedure
 
   module procedure to_separate_file
@@ -56,12 +40,14 @@
       allocate(comments(num_histograms))
 
       block
-        integer line
+        integer line, b
         do line  = 1, size(comments)
           associate( &
-            mode_frequency => string_t(int(100*maxval(histograms(line)%bin_frequency_))), &
-            range_min => string_t(histograms(line)%unmapped_min()), &
-            range_max =>  string_t(histograms(line)%unmapped_max()) &
+            mode_frequency => string_t(int(100*maxval( &
+              [( histograms(line)%bin_frequency(b), b = 1, size(histograms(line)%bin_count_) )] &
+            ))), &
+            range_min => string_t(histograms(line)%variable_min_), &
+            range_max =>  string_t(histograms(line)%variable_max_) &
            )
              comments(line) = "# " &
                // trim(histograms(line)%variable_name_) &
@@ -71,18 +57,14 @@
         end do
       end block
 
-      block
-        integer h
-        column_headings = "    bin    " // .cat. [("    " // string_t(histograms(h)%variable_name()) // "    ", h=1,num_histograms)]
-      end block
-
-      associate(num_bins =>  histograms(1)%num_bins())
-
+      associate(num_bins => size(histograms(1)%bin_value_))
         block 
           integer h, b ! histogram number, bin number
 
+          column_headings = "    bin    " // .cat. [("    " // string_t(histograms(h)%variable_name_) // "    ", h=1,num_histograms)]
+
           call_julienne_assert(num_bins .isAtLeast. 1)
-          call_julienne_assert(.all. (histograms(1)%num_bins() .equalsExpected. [(histograms(h)%num_bins() , h=1,size(histograms))]))
+          call_julienne_assert(.all. (num_bins .equalsExpected. [(size(histograms(h)%bin_value_) , h=1,size(histograms))]))
             
           allocate(columns(num_bins))
           do b = 1, num_bins
@@ -90,7 +72,6 @@
               .cat. [("  " // string_t(histograms(h)%bin_frequency(b)), h=1,num_histograms)]
           end do
         end block
-
       end associate
 
       file = file_t([comments, column_headings, columns])
@@ -99,67 +80,45 @@
  
   end procedure
 
-  pure function normalize(x, x_min, x_max) result(x_normalized)
-    real, intent(in) :: x(:,:,:,:), x_min, x_max
-    real, allocatable :: x_normalized(:,:,:,:)
-    call_julienne_assert(x_min/=x_max)
-    x_normalized = (x - x_min)/(x_max - x_min)
-  end function
-
   module procedure construct
+    histogram = construct_in_range(v, variable_name, num_bins, minval(v), maxval(v))
+  end procedure
 
-    integer i, j, k, n
-    integer, allocatable :: bin_count(:)
+  module procedure construct_in_range
+
+    integer b
     integer, parameter :: performance_threshold = 80
     real, parameter :: capture_maxval = 1.0001 ! ensure maxval(v_max) falls within the highest bin
-    real, allocatable :: v_mapped(:,:,:,:)
 
     histogram%variable_name_ = variable_name
-    histogram%unmapped_min_ = minval(v)
-    histogram%unmapped_max_ = maxval(v)
+    histogram%variable_min_  = v_min
+    histogram%variable_max_  = v_max
 
-    allocate(histogram%bin_frequency_(num_bins))
     allocate(histogram%bin_value_(num_bins))
-    allocate(bin_count(num_bins))
-
-    associate(v_min => (histogram%unmapped_min_), v_max => (histogram%unmapped_max_), cardinality => size(v))
-      if (raw)  then
-        v_mapped = v
-      else
-        v_mapped = normalize(v, v_min, v_max)
-      end if
-      associate(v_mapped_min => merge(v_min, 0., raw), v_mapped_max => capture_maxval*merge(v_max, 1., raw))
-        associate(dv => (v_mapped_max - v_mapped_min)/real(num_bins))
-          associate(v_bin_min => [(v_mapped_min + (i-1)*dv, i=1,num_bins)])
-            associate(v_bin_max => [v_bin_min(2:), v_mapped_max])
-              histogram%bin_value_ = 0.5*[v_bin_min + v_bin_max] ! switching to average yields problems likely related to roundoff
-              if (num_bins < performance_threshold) then
-                do concurrent(i = 1:num_bins)
-                  bin_count(i) = count(v_mapped >= v_bin_min(i) .and. v_mapped < v_bin_max(i))
-                end do
-              else
-                bin_count = 0
-                do i = 1,size(v_mapped,1)
-                  do j = 1,size(v_mapped,2)
-                    do k = 1,size(v_mapped,3)
-                      do n = 1,size(v_mapped,4)
-                        associate(bin => floor((v_mapped(i,j,k,n) - v_mapped_min)/dv) + 1)
-                          bin_count(bin) = bin_count(bin) + 1
-                        end associate
-                      end do
-                    end do
-                  end do
-                end do
-              end if
-              histogram%bin_frequency_ = real(bin_count) / real(cardinality)
-            end associate
+    allocate(histogram%bin_count_(num_bins))
+
+    associate(v_max_expanded => capture_maxval*v_max)
+      associate(dv => (v_max_expanded - v_min)/real(num_bins))
+        associate(v_bin_min => [(v_min + (b-1)*dv, b=1,num_bins)])
+          associate(v_bin_max => [v_bin_min(2:), v_max_expanded])
+            histogram%bin_value_ = 0.5*[v_bin_min + v_bin_max]
+            if (num_bins < performance_threshold) then
+              do concurrent(b = 1:num_bins) default(none) shared(histogram, v, v_bin_min, v_bin_max)
+                histogram%bin_count_(b) = count(v >= v_bin_min(b) .and. v < v_bin_max(b))
+              end do
+            else
+              histogram%bin_count_ = 0
+              do b = 1,size(v)
+                associate(bin => floor((v(b) - v_min)/dv) + 1)
+                  histogram%bin_count_(bin) = histogram%bin_count_(bin) + 1
+                end associate
+              end do
+            end if
           end associate
         end associate
       end associate
 #ifdef ASSERTIONS
-      associate(binned => sum(bin_count))
-        call_julienne_assert(cardinality .equalsExpected. binned)
-      end associate
+      call_julienne_assert(sum(histogram%bin_count_) .equalsExpected. size(v))
 #endif
     end associate
 
diff --git a/demo/test/driver.f90 b/demo/test/driver.f90
index 7eb3ab153..87f12e1a8 100644
--- a/demo/test/driver.f90
+++ b/demo/test/driver.f90
@@ -5,11 +5,13 @@ program test_suite_driver
   use julienne_m, only : test_fixture_t, test_harness_t
   use netCDF_file_test_m, only : netCDF_file_test_t
   use time_data_test_m, only : time_data_test_t
+  use histogram_test_m, only : histogram_test_t
   implicit none
 
   associate(test_harness => test_harness_t([ &
      test_fixture_t(netCDF_file_test_t()) &
     ,test_fixture_t(time_data_test_t()) &
+    ,test_fixture_t(histogram_test_t()) &
   ]))
     call test_harness%report_results
   end associate
diff --git a/demo/test/histogram_test_m.f90 b/demo/test/histogram_test_m.f90
new file mode 100644
index 000000000..2f98fd3f5
--- /dev/null
+++ b/demo/test/histogram_test_m.f90
@@ -0,0 +1,71 @@
+! Copyright (c) 2022-2025, The Regents of the University of California and Sourcery Institute
+! Terms of use are as specified in LICENSE.txt
+
+module histogram_test_m
+  !! Unit test for the histogram subroutine
+  use julienne_m, only : &
+     operator(.all.) &
+    ,operator(.also.) &
+    ,operator(.approximates.) &
+    ,operator(.equalsExpected.) &
+    ,operator(.within.) &
+    ,operator(//) &
+    ,string_t &
+    ,test_t &
+    ,test_result_t &
+    ,test_description_t &
+    ,test_diagnosis_t
+  use histogram_m, only : histogram_t, to_file
+  implicit none
+
+  private
+  public :: histogram_test_t
+
+  type, extends(test_t) :: histogram_test_t
+  contains
+    procedure, nopass :: subject
+    procedure, nopass :: results
+  end type
+
+contains
+
+  pure function subject() result(specimen)
+    character(len=:), allocatable :: specimen
+    specimen = "The histogram_t type" 
+  end function
+
+  function results() result(test_results)
+    type(histogram_test_t) histogram_test
+    type(test_result_t), allocatable :: test_results(:)
+
+    test_results = histogram_test%run([ &
+       test_description_t("constructing a uniform distribution from uniformly distributed points", check_uniform_distribution) &
+    ])
+  end function
+
+  type(test_diagnosis_t) function check_uniform_distribution() result(test_diagnosis)
+      
+    real u(1,2,3,4), w(  size(u,1), size(u,2), size(u,3), size(u,4))
+    real, target ::  v(2*size(u,1), size(u,2), size(u,3), size(u,4))
+    real, pointer :: v_1D(:)
+    integer, parameter :: num_points = size(v), num_bins = 7, remainder = mod(num_points,num_bins)
+    real, parameter :: v_min = -10., v_max = 10., dv = (v_max - v_min)/num_bins
+    integer b, p
+
+    v_1D(1:num_points) => v
+    v_1D = [( [(v_min + dv/2 + (b-1)*dv, p = 1, num_points/num_bins + merge(0, 1, b > remainder))], b = 1, num_bins )]
+    u = v(1:1,:,:,:)
+    w = v(2:2,:,:,:)
+
+    associate(block_distribution => num_points/num_bins + [(merge(0, 1, p > remainder), p = 1, num_bins)])
+      associate(histogram => histogram_t([v], "uniform", num_bins))
+        test_diagnosis = .all. (histogram%bin_count() .equalsExpected. block_distribution)
+      end associate
+      associate(histogram => histogram_t([u,v,w], "uniform", num_bins))
+        test_diagnosis = test_diagnosis .also. (.all. (histogram%bin_count() .equalsExpected. 2*block_distribution))
+      end associate
+    end associate
+
+  end function
+
+end module histogram_test_m
diff --git a/demo/test/netCDF_file_test_m.F90 b/demo/test/netCDF_file_test_m.f90
similarity index 100%
rename from demo/test/netCDF_file_test_m.F90
rename to demo/test/netCDF_file_test_m.f90
diff --git a/demo/test/time_data_test_m.F90 b/demo/test/time_data_test_m.f90
similarity index 98%
rename from demo/test/time_data_test_m.F90
rename to demo/test/time_data_test_m.f90
index 8ef65074e..a651149df 100644
--- a/demo/test/time_data_test_m.F90
+++ b/demo/test/time_data_test_m.f90
@@ -1,8 +1,6 @@
 ! Copyright (c) 2022-2025, The Regents of the University of California and Sourcery Institute
 ! Terms of use are as specified in LICENSE.txt
 
-#include "language-support.F90"
-
 module time_data_test_m
   !! Unit test for the time_data subroutine
   use julienne_m, only : &
diff --git a/src/fiats/training_data_files_m.f90 b/src/fiats/training_data_files_m.f90
index 5f1c65a0f..e9e12ef30 100644
--- a/src/fiats/training_data_files_m.f90
+++ b/src/fiats/training_data_files_m.f90
@@ -15,6 +15,8 @@ module training_data_files_m
     procedure :: to_json
     procedure :: fully_qualified_inputs_files
     procedure :: fully_qualified_outputs_files
+    procedure :: inputs_files
+    procedure :: outputs_files
     procedure :: fully_qualified_time_file
     procedure :: path
     generic :: operator(==) => equals
@@ -64,12 +66,24 @@ pure module function fully_qualified_inputs_files(self) result(names)
       type(string_t), allocatable :: names(:)
     end function
 
+    pure module function inputs_files(self) result(names)
+      implicit none
+      class(training_data_files_t), intent(in) :: self
+      type(string_t), allocatable :: names(:)
+    end function
+
     pure module function fully_qualified_outputs_files(self) result(names)
       implicit none
       class(training_data_files_t), intent(in) :: self
       type(string_t), allocatable :: names(:)
     end function
 
+    pure module function outputs_files(self) result(names)
+      implicit none
+      class(training_data_files_t), intent(in) :: self
+      type(string_t), allocatable :: names(:)
+    end function
+
     pure module function fully_qualified_time_file(self) result(name)
       implicit none
       class(training_data_files_t), intent(in) :: self
diff --git a/src/fiats/training_data_files_s.F90 b/src/fiats/training_data_files_s.F90
index 900323e12..efab6d99f 100644
--- a/src/fiats/training_data_files_s.F90
+++ b/src/fiats/training_data_files_s.F90
@@ -84,4 +84,14 @@
     names = [(self%path_ // "/" // self%outputs_prefix_ //self%infixes_(i) // suffix, i=1, size(self%infixes_))]
   end procedure
 
+  module procedure inputs_files
+    integer i
+    names = [(self%inputs_prefix_ //self%infixes_(i) // suffix, i=1, size(self%infixes_))]
+  end procedure
+
+  module procedure outputs_files
+    integer i
+    names = [(self%outputs_prefix_ //self%infixes_(i) // suffix, i=1, size(self%infixes_))]
+  end procedure
+
 end submodule training_data_files_s