refac(histogram): store in sourcery file_t object

cloud-microphysics/app/tensor-statistics.f90

@@ -7,15 +7,15 @@ program tensor_statistics
   !! 2. Saves the resulting statistics to text files with space-separated columns and column labels.
 
   ! External dependencies:
-  use sourcery_m, only : command_line_t
+  use sourcery_m, only : command_line_t, file_t, string_t
   use assert_m, only : assert, intrinsic_array_t
   use inference_engine_m, only : rkind, ubounds_t
   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 histogram_m, only : histogram_t
+  use histogram_m, only : histogram_t, to_file
   implicit none
 
   character(len=*), parameter :: usage = &
@@ -42,7 +42,7 @@ program tensor_statistics
   print *,"---"
   print *
   print *,"---> Please see the *.plt files for the tensor ranges and histograms. <---"
-  print *,"---> Execute `gnuplot app/gnuplot.in` to graph the histograms.        <---"
+  print *,"---> Execute `gnuplot app/gnuplot.inp` to graph the histograms.       <---"
   print *
   print *,"______ tensor_statistics done_______"
 
@@ -63,7 +63,6 @@ subroutine get_command_line_arguments(base_name, num_bins, start_step, end_step,
     end_string = command_line%flag_value("--end")
     stride_string = command_line%flag_value("--stride")
 
-
     associate(required_arguments => len(base_name)/=0 .and. len(bins_string)/=0)
        if (.not. required_arguments) error stop usage 
     end associate
@@ -126,49 +125,34 @@ subroutine compute_histograms(base_name)
            ubounds_t(ubound(pressure_in)), ubounds_t(ubound(temperature_in)) &
           ]
 
-        block
-          integer bin, h, file_unit
-          integer, parameter :: num_inputs = 7
-          type(histogram_t) :: histogram(num_inputs)
-
-          print *,"Calculating input tensor histograms"
-
-          histogram(1) = histogram_t(pressure_in, "pressure", num_bins)
-          histogram(2) = histogram_t(potential_temperature_in, '"potential temperature"', num_bins)
-          histogram(3) = histogram_t(temperature_in, "temperature", num_bins)
-          histogram(4) = histogram_t(qv_in, "qv", num_bins)
-          histogram(5) = histogram_t(qc_in, "qc", num_bins)
-          histogram(6) = histogram_t(qr_in, "qr", num_bins)
-          histogram(7) = histogram_t(qs_in, "qs", num_bins)
-
-          associate(input_tensor_stats_file_name => base_name // "_inputs_stats.plt")
-            print *,"Writing input tensor statistics to " // input_tensor_stats_file_name
-            open(newunit=file_unit, file=input_tensor_stats_file_name, status="unknown")
-          end associate
-
-          do h = 1, size(histogram)
-            write(file_unit,*), &
-             "# unmapped range for ", histogram(h)%variable_name(),":", histogram(h)%unmapped_range()
-          end do
-
-          write(file_unit,'(5x,a,8(10x,a))'),"bin", (histogram(h)%variable_name(), h=1,size(histogram)) ! column headings
-
-          do bin = 1, histogram(1)%num_bins()
-            write(file_unit, *), histogram(1)%bin_midpoint(bin), (histogram(h)%bin_frequency(bin), h=1,size(histogram))
-          end do
-
-          close(file_unit)
-        end block
+        print *,"Calculating input tensor histograms and write to file"
+
+        associate(histograms => [ &
+          histogram_t(pressure_in, "pressure", num_bins) &
+         ,histogram_t(potential_temperature_in, '"potential temperature"', num_bins) &
+         ,histogram_t(temperature_in, "temperature", num_bins) &
+         ,histogram_t(qv_in, "qv", num_bins) &
+         ,histogram_t(qc_in, "qc", num_bins) &
+         ,histogram_t(qr_in, "qr", num_bins) &
+         ,histogram_t(qs_in, "qs", num_bins) &
+        ])
+          block
+            type(file_t) histograms_file
+
+            histograms_file = to_file(histograms)
+            call histograms_file%write_lines(string_t(base_name // "_inputs_stats.plt"))
+          end block
+        end associate
       end associate
     end associate
 
     associate(network_output_file_name => base_name // "_output.nc")
 
-      print *,"Reading network outputs from " // network_output_file_name
+      print *, "Reading network outputs from " // network_output_file_name
 
       associate(network_output_file => netCDF_file_t(network_output_file_name))
-        call network_output_file%input("potential_temperature", potential_temperature_out)
         ! Skipping the following unnecessary outputs: pressure, temperature, precipitation, snowfall
+        call network_output_file%input("potential_temperature", potential_temperature_out)
         call network_output_file%input("qv", qv_out)
         call network_output_file%input("qc", qc_out)
         call network_output_file%input("qr", qr_out)
@@ -206,39 +190,23 @@ subroutine compute_histograms(base_name)
       call assert(.not. any(ieee_is_nan(dqr_dt)), ".not. any(ieee_is_nan(dqr_dt)")
       call assert(.not. any(ieee_is_nan(dqs_dt)), ".not. any(ieee_is_nan(dqs_dt)")
 
-      block
-        integer bin, h, file_unit
-        integer, parameter :: num_outputs = 5
-        type(histogram_t) :: histogram(num_outputs)
+      print *,"Calculating output tensor histograms"
 
-        print *,"Calculating output tensor histograms"
-
-        histogram(1) = histogram_t(dpt_dt, "d(pt)/dt", num_bins)
-        histogram(2) = histogram_t(dqv_dt, "d(qv)/dt", num_bins)
-        histogram(3) = histogram_t(dqc_dt, "d(qc)/dt", num_bins)
-        histogram(4) = histogram_t(dqr_dt, "d(qr)/dt", num_bins)
-        histogram(5) = histogram_t(dqs_dt, "d(qs)/dt", num_bins)
-
-        associate(output_tensor_stats_file_name => base_name // "_outputs_stats.plt")
-          print *,"Writing output tensor statistics to " // output_tensor_stats_file_name
-          open(newunit=file_unit, file=output_tensor_stats_file_name, status="unknown")
-        end associate
-
-        do h = 1, size(histogram)
-          write(file_unit,*), &
-           "# unmapped range for ", histogram(h)%variable_name(),":", histogram(h)%unmapped_range()
-        end do
-
-        write(file_unit,'(5x,a,5(10x,a))'),"bin", (histogram(h)%variable_name(), h=1,size(histogram)) ! column headings
-
-        do bin = 1, histogram(1)%num_bins()
-          write(file_unit, *), histogram(1)%bin_midpoint(bin), (histogram(h)%bin_frequency(bin), h=1,size(histogram))
-        end do
+      associate(histograms => [ &
+         histogram_t(dpt_dt, "d(pt)/dt", num_bins) &
+        ,histogram_t(dqv_dt, "d(qv)/dt", num_bins) &
+        ,histogram_t(dqc_dt, "d(qc)/dt", num_bins) &
+        ,histogram_t(dqr_dt, "d(qr)/dt", num_bins) &
+        ,histogram_t(dqs_dt, "d(qs)/dt", num_bins) &
+      ])
+        block
+          type(file_t) histograms_file
 
-        close(file_unit)
-      end block
+          histograms_file = to_file(histograms)
+          call histograms_file%write_lines(string_t(base_name // "_outputs_stats.plt"))
+        end block
+      end associate
     end associate
-
   end subroutine
 
 end program tensor_statistics

cloud-microphysics/fpm.toml

@@ -6,6 +6,6 @@ maintainer = "rouson@lbl.gov"
 
 [dependencies]
 assert = {git = "https://github.com/sourceryinstitute/assert", tag = "1.5.0"}
-sourcery = {git = "https://github.com/sourceryinstitute/sourcery", tag = "4.5.0"}
+sourcery = {git = "https://github.com/sourceryinstitute/sourcery", tag = "4.5.3"}
 inference-engine = {path = "../"}
 netcdf-interfaces = {git = "https://github.com/LKedward/netcdf-interfaces.git", rev = "d2bbb71ac52b4e346b62572b1ca1620134481096"}

cloud-microphysics/src/histogram_m.f90

@@ -1,19 +1,24 @@
 ! Copyright (c), The Regents of the University of California
 ! Terms of use are as specified in LICENSE.txt
 module histogram_m
+  !! Generate and represent histograms
+  use sourcery_m, only : file_t
   implicit none
 
   private
-  public :: histogram_t
+  public :: histogram_t, to_file
 
   type histogram_t
+    !! encapsulate the primary data associated with histograms
     private
     character(len=:), allocatable :: variable_name_
     real unmapped_min_, unmapped_max_
     real, allocatable :: frequency_(:), bin_midpoint_(:)
   contains
     procedure variable_name
     procedure unmapped_range
+    procedure unmapped_min
+    procedure unmapped_max
     procedure num_bins
     procedure bin_midpoint
     procedure bin_frequency
@@ -33,11 +38,18 @@ pure module function histogram_on_unit_interval(v, variable_name, num_bins) resu
 
   interface
 
-  pure module function num_bins(self) result(bins)
-    implicit none
-    class(histogram_t), intent(in) :: self
-    integer bins
-  end function
+    !pure module function to_file(histograms) result(file)
+    module function to_file(histograms) result(file)
+      implicit none
+      type(histogram_t), intent(in) :: histograms(:)
+      type(file_t) file
+    end function
+
+    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
@@ -52,7 +64,19 @@ pure module function unmapped_range(self) result(raw_range)
       real raw_range(num_end_points)
     end function
 
-    pure module function bin_midpoint(self, bin) result(midpoint)
+    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_midpoint(self, bin) result(midpoint)
       implicit none
       class(histogram_t), intent(in) :: self
       integer, intent(in) :: bin

cloud-microphysics/src/histogram_s.f90

@@ -3,6 +3,7 @@
 submodule(histogram_m) histogram_s
   use assert_m, only : assert, intrinsic_array_t
   use kind_parameters_m, only : rkind
+  use sourcery_m, only : string_t, operator(.cat.)
   implicit none
 
 contains
@@ -15,6 +16,14 @@
     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_midpoint_)
   end procedure
@@ -27,6 +36,51 @@
     frequency = self%frequency_(bin)
   end procedure
 
+  module procedure to_file
+    type(string_t), allocatable :: comments(:), columns(:)
+    type(string_t) column_headings
+
+    associate(num_histograms => size(histograms))
+
+      allocate(comments(num_histograms))
+
+      block 
+        integer line
+        do line  = 1, size(comments)
+          comments(line) = "# " // histograms(line)%variable_name() // " range: " // &
+            string_t(histograms(line)%unmapped_min()) // "  " // string_t(histograms(line)%unmapped_max())
+        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())
+
+        block 
+          integer h, b ! histogram number, bin number
+
+          call assert(num_bins > 0, "histogram_s(to_file): num_bins > 0")
+          call assert(all(histograms(1)%num_bins() == [(histograms(h)%num_bins() , h=1,size(histograms))]), &
+            "histogram_s(to_file): uniform number of bins")
+
+          allocate(columns(num_bins))
+          do b = 1, num_bins
+            columns(b) =  string_t(histograms(1)%bin_midpoint(b)) // &
+              .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])
+
+    end associate
+
+  end procedure
+
   pure function normalize(x, x_min, x_max) result(x_normalized)
     implicit none
     real(rkind), intent(in) :: x(:,:,:,:), x_min, x_max