diff --git a/.readthedocs.yaml b/.readthedocs.yaml
new file mode 100644
index 000000000..f83760cce
--- /dev/null
+++ b/.readthedocs.yaml
@@ -0,0 +1,29 @@
+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Set the version of Python and other tools you might need
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.7"
+    # You can also specify other tool versions:
+    # nodejs: "19"
+    # rust: "1.64"
+    # golang: "1.19"
+
+# Build documentation in the docs/ directory with Sphinx
+sphinx:
+   configuration: doc/source/conf.py
+
+# If using Sphinx, optionally build your docs in additional formats such as PDF
+# formats:
+#    - pdf
+
+# Optionally declare the Python requirements required to build your docs
+python:
+   install:
+   - requirements: doc/requirements.txt
diff --git a/cicecore/cicedyn/analysis/ice_history.F90 b/cicecore/cicedyn/analysis/ice_history.F90
index 54b6ce934..598f05a61 100644
--- a/cicecore/cicedyn/analysis/ice_history.F90
+++ b/cicecore/cicedyn/analysis/ice_history.F90
@@ -2225,7 +2225,7 @@ subroutine accum_hist (dt)
       n4Dfcum = n4Dscum + num_avail_hist_fields_4Df ! should equal num_avail_hist_fields_tot
 
       do ns = 1,nstreams
-         if (.not. hist_avg) then  ! write snapshots
+         if (.not. hist_avg(ns)) then  ! write snapshots
            do n = 1,n2D
               if (avail_hist_fields(n)%vhistfreq == histfreq(ns)) &
                   a2D(:,:,n,:) = c0
diff --git a/cicecore/cicedyn/analysis/ice_history_shared.F90 b/cicecore/cicedyn/analysis/ice_history_shared.F90
index 70aa5e14c..f4e1f3ebf 100644
--- a/cicecore/cicedyn/analysis/ice_history_shared.F90
+++ b/cicecore/cicedyn/analysis/ice_history_shared.F90
@@ -38,7 +38,7 @@ module ice_history_shared
       integer (kind=int_kind), public :: history_precision
 
       logical (kind=log_kind), public :: &
-         hist_avg  ! if true, write averaged data instead of snapshots
+         hist_avg(max_nstrm)  ! if true, write averaged data instead of snapshots
 
       character (len=char_len_long), public :: &
          history_file  , & ! output file for history
@@ -743,7 +743,7 @@ subroutine construct_filename(ncfile,suffix,ns)
               imonth,'-',iday,'-',isec,'.',trim(suffix)
         else
 
-           if (hist_avg) then
+           if (hist_avg(ns)) then
               if (histfreq(ns) == '1' .or. histfreq(ns) == 'h'.or.histfreq(ns) == 'H') then
                  ! do nothing
               elseif (new_year) then
@@ -763,7 +763,7 @@ subroutine construct_filename(ncfile,suffix,ns)
 !echmod ! of other groups (including RASM which uses CESMCOUPLED)
 !echmod         if (ns > 1) write(cstream,'(i1.1)') ns-1
 
-           if (hist_avg) then    ! write averaged data
+           if (hist_avg(ns)) then    ! write averaged data
               if (histfreq(ns) == '1' .and. histfreq_n(ns) == 1)  then ! timestep
                  write(ncfile,'(a,a,i4.4,a,i2.2,a,i2.2,a,i5.5,a,a)')  &
                        history_file(1:lenstr(history_file))//trim(cstream),'_inst.', &
diff --git a/cicecore/cicedyn/general/ice_init.F90 b/cicecore/cicedyn/general/ice_init.F90
index 4c8fb1fee..2f2e5802b 100644
--- a/cicecore/cicedyn/general/ice_init.F90
+++ b/cicecore/cicedyn/general/ice_init.F90
@@ -322,7 +322,7 @@ subroutine input_data
       histfreq(5) = 'y'      ! output frequency option for different streams
       histfreq_n(:) = 1      ! output frequency
       histfreq_base = 'zero' ! output frequency reference date
-      hist_avg = .true.      ! if true, write time-averages (not snapshots)
+      hist_avg(:) = .true.   ! if true, write time-averages (not snapshots)
       history_format = 'default' ! history file format
       history_dir  = './'    ! write to executable dir for default
       history_file = 'iceh'  ! history file name prefix
@@ -901,7 +901,7 @@ subroutine input_data
       enddo
       call broadcast_array(histfreq_n,            master_task)
       call broadcast_scalar(histfreq_base,        master_task)
-      call broadcast_scalar(hist_avg,             master_task)
+      call broadcast_array(hist_avg,              master_task)
       call broadcast_scalar(history_dir,          master_task)
       call broadcast_scalar(history_file,         master_task)
       call broadcast_scalar(history_precision,    master_task)
@@ -2311,8 +2311,7 @@ subroutine input_data
          write(nu_diag,1033) ' histfreq         = ', histfreq(:)
          write(nu_diag,1023) ' histfreq_n       = ', histfreq_n(:)
          write(nu_diag,1031) ' histfreq_base    = ', trim(histfreq_base)
-         write(nu_diag,1011) ' hist_avg         = ', hist_avg
-         if (.not. hist_avg) write(nu_diag,1039) ' History data will be snapshots'
+         write(nu_diag,*)    ' hist_avg         = ', hist_avg(:)
          write(nu_diag,1031) ' history_dir      = ', trim(history_dir)
          write(nu_diag,1031) ' history_file     = ', trim(history_file)
          write(nu_diag,1021) ' history_precision= ', history_precision
diff --git a/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90 b/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90
index 9df51635d..526d0d96d 100644
--- a/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90
+++ b/cicecore/cicedyn/infrastructure/io/io_binary/ice_history_write.F90
@@ -157,7 +157,7 @@ subroutine ice_write_hist(ns)
             write (nu_hdr, 995) nrec,trim(avail_hist_fields(n)%vname), &
                trim(avail_hist_fields(n)%vcomment)
 
-            if (histfreq(ns) == '1' .or. .not. hist_avg         &
+            if (histfreq(ns) == '1' .or. .not. hist_avg(ns)     &
                 .or. write_ic                                   &
                 .or. n==n_divu(ns)      .or. n==n_shear(ns)     &  ! snapshots
                 .or. n==n_sig1(ns)      .or. n==n_sig2(ns)      &
@@ -187,7 +187,7 @@ subroutine ice_write_hist(ns)
             write (nu_hdr, 994) nrec,trim(avail_hist_fields(n)%vname), &
                trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn
 
-            if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then
+            if (histfreq(ns) == '1' .or. .not. hist_avg(ns) .or. write_ic) then
                write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), &
                   'time_rep','instantaneous'
             else
@@ -211,7 +211,7 @@ subroutine ice_write_hist(ns)
             write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), &
                trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),k
 
-            if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then
+            if (histfreq(ns) == '1' .or. .not. hist_avg(ns) .or. write_ic) then
                write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), &
                   'time_rep','instantaneous'
             else
@@ -235,7 +235,7 @@ subroutine ice_write_hist(ns)
             write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), &
                trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k
 
-            if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then
+            if (histfreq(ns) == '1' .or. .not. hist_avg(ns) .or. write_ic) then
                write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), &
                   'time_rep','instantaneous'
             else
@@ -259,7 +259,7 @@ subroutine ice_write_hist(ns)
             write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), &
                trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k
 
-            if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then
+            if (histfreq(ns) == '1' .or. .not. hist_avg(ns) .or. write_ic) then
                write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), &
                   'time_rep','instantaneous'
             else
@@ -283,7 +283,7 @@ subroutine ice_write_hist(ns)
             write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), &
                trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k
 
-            if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then
+            if (histfreq(ns) == '1' .or. .not. hist_avg(ns) .or. write_ic) then
                write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), &
                   'time_rep','instantaneous'
             else
@@ -308,7 +308,7 @@ subroutine ice_write_hist(ns)
             write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), &
                trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k
 
-            if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then
+            if (histfreq(ns) == '1' .or. .not. hist_avg(ns) .or. write_ic) then
                write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), &
                   'time_rep','instantaneous'
             else
@@ -334,7 +334,7 @@ subroutine ice_write_hist(ns)
             write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), &
                trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k
 
-            if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then
+            if (histfreq(ns) == '1' .or. .not. hist_avg(ns) .or. write_ic) then
                write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), &
                   'time_rep','instantaneous'
             else
@@ -360,7 +360,7 @@ subroutine ice_write_hist(ns)
             write (nu_hdr, 993) nrec,trim(avail_hist_fields(n)%vname), &
                trim(avail_hist_fields(n)%vdesc),trim(avail_hist_fields(n)%vunit),nn,k
 
-            if (histfreq(ns) == '1' .or. .not. hist_avg .or. write_ic) then
+            if (histfreq(ns) == '1' .or. .not. hist_avg(ns) .or. write_ic) then
                write (nu_hdr, 996) nrec,trim(avail_hist_fields(n)%vname), &
                   'time_rep','instantaneous'
             else
diff --git a/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90 b/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90
index 10d750300..25178ed6e 100644
--- a/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90
+++ b/cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90
@@ -159,10 +159,10 @@ subroutine ice_write_hist (ns)
       ! define dimensions
       !-----------------------------------------------------------------
 
-        if (hist_avg .and. .not. write_ic) then
-          status = nf90_def_dim(ncid,'d2',2,boundid)
+        if (hist_avg(ns) .and. .not. write_ic) then
+          status = nf90_def_dim(ncid,'nbnd',2,boundid)
           if (status /= nf90_noerr) call abort_ice(subname// &
-             'ERROR: defining dim d2')
+             'ERROR: defining dim nbnd')
         endif
 
         status = nf90_def_dim(ncid,'ni',nx_global,imtid)
@@ -213,7 +213,7 @@ subroutine ice_write_hist (ns)
         if (status /= nf90_noerr) call abort_ice(subname// &
                       'ERROR: defining var time')
 
-        status = nf90_put_att(ncid,varid,'long_name','model time')
+        status = nf90_put_att(ncid,varid,'long_name','time')
         if (status /= nf90_noerr) call abort_ice(subname// &
                       'ice Error: time long_name')
 
@@ -230,7 +230,7 @@ subroutine ice_write_hist (ns)
            if (status /= nf90_noerr) call abort_ice(subname// &
                          'ERROR: time calendar')
         elseif (days_per_year == 365 .and. .not.use_leap_years ) then
-           status = nf90_put_att(ncid,varid,'calendar','NoLeap')
+           status = nf90_put_att(ncid,varid,'calendar','noleap')
            if (status /= nf90_noerr) call abort_ice(subname// &
                          'ERROR: time calendar')
         elseif (use_leap_years) then
@@ -241,7 +241,7 @@ subroutine ice_write_hist (ns)
            call abort_ice(subname//'ERROR: invalid calendar settings')
         endif
 
-        if (hist_avg .and. .not. write_ic) then
+        if (hist_avg(ns) .and. .not. write_ic) then
           status = nf90_put_att(ncid,varid,'bounds','time_bounds')
           if (status /= nf90_noerr) call abort_ice(subname// &
                       'ERROR: time bounds')
@@ -251,14 +251,14 @@ subroutine ice_write_hist (ns)
       ! Define attributes for time bounds if hist_avg is true
       !-----------------------------------------------------------------
 
-        if (hist_avg .and. .not. write_ic) then
+        if (hist_avg(ns) .and. .not. write_ic) then
           dimid(1) = boundid
           dimid(2) = timid
           status = nf90_def_var(ncid,'time_bounds',lprecision,dimid(1:2),varid)
           if (status /= nf90_noerr) call abort_ice(subname// &
                         'ERROR: defining var time_bounds')
           status = nf90_put_att(ncid,varid,'long_name', &
-                                'boundaries for time-averaging interval')
+                                'time interval endpoints')
           if (status /= nf90_noerr) call abort_ice(subname// &
                         'ERROR: time_bounds long_name')
           write(cdate,'(i8.8)') idate0
@@ -268,6 +268,22 @@ subroutine ice_write_hist (ns)
           status = nf90_put_att(ncid,varid,'units',title)
           if (status /= nf90_noerr) call abort_ice(subname// &
                         'ERROR: time_bounds units')
+          if (days_per_year == 360) then
+             status = nf90_put_att(ncid,varid,'calendar','360_day')
+             if (status /= nf90_noerr) call abort_ice(subname// &
+                           'ERROR: time calendar')
+          elseif (days_per_year == 365 .and. .not.use_leap_years ) then
+             status = nf90_put_att(ncid,varid,'calendar','noleap')
+             if (status /= nf90_noerr) call abort_ice(subname// &
+                           'ERROR: time calendar')
+          elseif (use_leap_years) then
+             status = nf90_put_att(ncid,varid,'calendar','Gregorian')
+             if (status /= nf90_noerr) call abort_ice(subname// &
+                           'ERROR: time calendar')
+          else
+             call abort_ice(subname//'ERROR: invalid calendar settings')
+          endif
+
         endif
 
       !-----------------------------------------------------------------
@@ -745,7 +761,7 @@ subroutine ice_write_hist (ns)
       ! write time_bounds info
       !-----------------------------------------------------------------
 
-        if (hist_avg .and. .not. write_ic) then
+        if (hist_avg(ns) .and. .not. write_ic) then
           status = nf90_inq_varid(ncid,'time_bounds',varid)
           if (status /= nf90_noerr) call abort_ice(subname// &
                         'ERROR: getting time_bounds id')
@@ -1279,7 +1295,7 @@ subroutine ice_write_hist_attrs(ncid, varid, hfield, ns)
       call ice_write_hist_fill(ncid,varid,hfield%vname,history_precision)
 
       ! Add cell_methods attribute to variables if averaged
-      if (hist_avg .and. .not. write_ic) then
+      if (hist_avg(ns) .and. .not. write_ic) then
          if    (TRIM(hfield%vname(1:4))/='sig1' &
            .and.TRIM(hfield%vname(1:4))/='sig2' &
            .and.TRIM(hfield%vname(1:9))/='sistreave' &
@@ -1292,7 +1308,7 @@ subroutine ice_write_hist_attrs(ncid, varid, hfield, ns)
       endif
 
       if ((histfreq(ns) == '1' .and. histfreq_n(ns) == 1) &
-          .or..not. hist_avg                              &
+          .or..not. hist_avg(ns)                          &
           .or. write_ic                                   &
           .or.TRIM(hfield%vname(1:4))=='divu' &
           .or.TRIM(hfield%vname(1:5))=='shear' &
diff --git a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90 b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90
index 25f9850ce..35ec7bed2 100644
--- a/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90
+++ b/cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90
@@ -195,8 +195,8 @@ subroutine ice_write_hist (ns)
       ! define dimensions
       !-----------------------------------------------------------------
 
-        if (hist_avg .and. .not. write_ic) then
-          status = pio_def_dim(File,'d2',2,boundid)
+        if (hist_avg(ns) .and. .not. write_ic) then
+          status = pio_def_dim(File,'nbnd',2,boundid)
         endif
 
         status = pio_def_dim(File,'ni',nx_global,imtid)
@@ -215,7 +215,7 @@ subroutine ice_write_hist (ns)
       !-----------------------------------------------------------------
 
         status = pio_def_var(File,'time',pio_double,(/timid/),varid)
-        status = pio_put_att(File,varid,'long_name','model time')
+        status = pio_put_att(File,varid,'long_name','time')
 
         write(cdate,'(i8.8)') idate0
         write(title,'(a,a4,a1,a2,a1,a2,a1,i2.2,a1,i2.2,a1,i2.2)') 'days since ', &
@@ -226,24 +226,35 @@ subroutine ice_write_hist (ns)
         if (days_per_year == 360) then
            status = pio_put_att(File,varid,'calendar','360_day')
         elseif (days_per_year == 365 .and. .not.use_leap_years ) then
-           status = pio_put_att(File,varid,'calendar','NoLeap')
+           status = pio_put_att(File,varid,'calendar','noleap')
         elseif (use_leap_years) then
            status = pio_put_att(File,varid,'calendar','Gregorian')
         else
            call abort_ice(subname//'ERROR: invalid calendar settings')
         endif
 
-        if (hist_avg .and. .not. write_ic) then
+        if (hist_avg(ns) .and. .not. write_ic) then
           status = pio_put_att(File,varid,'bounds','time_bounds')
         endif
 
         ! Define attributes for time_bounds if hist_avg is true
-        if (hist_avg .and. .not. write_ic) then
+        if (hist_avg(ns) .and. .not. write_ic) then
           dimid2(1) = boundid
           dimid2(2) = timid
           status = pio_def_var(File,'time_bounds',pio_double,dimid2,varid)
           status = pio_put_att(File,varid,'long_name', &
-                                'boundaries for time-averaging interval')
+                                'time interval endpoints')
+
+          if (days_per_year == 360) then
+             status = pio_put_att(File,varid,'calendar','360_day')
+          elseif (days_per_year == 365 .and. .not.use_leap_years ) then
+             status = pio_put_att(File,varid,'calendar','noleap')
+          elseif (use_leap_years) then
+             status = pio_put_att(File,varid,'calendar','Gregorian')
+          else
+             call abort_ice(subname//'ERROR: invalid calendar settings')
+          endif
+
           write(cdate,'(i8.8)') idate0
           write(title,'(a,a4,a1,a2,a1,a2,a1,i2.2,a1,i2.2,a1,i2.2)') 'days since ', &
                 cdate(1:4),'-',cdate(5:6),'-',cdate(7:8),' ', &
@@ -702,7 +713,7 @@ subroutine ice_write_hist (ns)
       ! write time_bounds info
       !-----------------------------------------------------------------
 
-        if (hist_avg .and. .not. write_ic) then
+        if (hist_avg(ns) .and. .not. write_ic) then
           status = pio_inq_varid(File,'time_bounds',varid)
           time_bounds=(/time_beg(ns),time_end(ns)/)
           bnd_start  = (/1,1/)
@@ -1250,7 +1261,7 @@ subroutine ice_write_hist_attrs(File, varid, hfield, ns)
       call ice_write_hist_fill(File,varid,hfield%vname,history_precision)
 
       ! Add cell_methods attribute to variables if averaged
-      if (hist_avg .and. .not. write_ic) then
+      if (hist_avg(ns) .and. .not. write_ic) then
          if    (TRIM(hfield%vname(1:4))/='sig1' &
            .and.TRIM(hfield%vname(1:4))/='sig2' &
            .and.TRIM(hfield%vname(1:9))/='sistreave' &
@@ -1261,7 +1272,7 @@ subroutine ice_write_hist_attrs(File, varid, hfield, ns)
       endif
 
       if ((histfreq(ns) == '1' .and. histfreq_n(ns) == 1) &
-          .or..not. hist_avg                              &
+          .or..not. hist_avg(ns)                          &
           .or. write_ic                                   &
           .or.TRIM(hfield%vname(1:4))=='divu' &
           .or.TRIM(hfield%vname(1:5))=='shear' &
diff --git a/configuration/scripts/ice_in b/configuration/scripts/ice_in
index 32db0270b..8dc046da5 100644
--- a/configuration/scripts/ice_in
+++ b/configuration/scripts/ice_in
@@ -48,7 +48,7 @@
     histfreq       = 'm','x','x','x','x'
     histfreq_n     =  1 , 1 , 1 , 1 , 1
     histfreq_base  = 'zero'
-    hist_avg       = .true.
+    hist_avg       = .true.,.true.,.true.,.true.,.true.
     history_dir    = './history/'
     history_file   = 'iceh'
     history_precision = 4
diff --git a/configuration/scripts/options/set_nml.histinst b/configuration/scripts/options/set_nml.histinst
index f2f0995c8..31d566d76 100644
--- a/configuration/scripts/options/set_nml.histinst
+++ b/configuration/scripts/options/set_nml.histinst
@@ -1 +1 @@
-hist_avg        = .false.
+hist_avg        = .false.,.false.,.false.,.false.,.false.
diff --git a/configuration/scripts/options/set_nml.qc b/configuration/scripts/options/set_nml.qc
index feefb376d..5de4dd28e 100644
--- a/configuration/scripts/options/set_nml.qc
+++ b/configuration/scripts/options/set_nml.qc
@@ -12,5 +12,5 @@ dumpfreq_n     = 12
 diagfreq       = 24
 histfreq       = 'd','x','x','x','x'
 f_hi           = 'd'
-hist_avg       = .false.
+hist_avg       = .false.,.false.,.false.,.false.,.false.
 distribution_wght = 'blockall'
diff --git a/configuration/scripts/options/set_nml.run3dt b/configuration/scripts/options/set_nml.run3dt
index 4ff27ce22..11a8c0f85 100644
--- a/configuration/scripts/options/set_nml.run3dt
+++ b/configuration/scripts/options/set_nml.run3dt
@@ -2,6 +2,6 @@ npt_unit        = '1'
 npt             = 3
 dump_last       = .true.
 histfreq        = '1','x','x','x','x'
-hist_avg        = .false.
+hist_avg        = .false.,.false.,.false.,.false.,.false.
 f_uvel          = '1' 
 f_vvel          = '1' 
diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst
index 000004bb9..0c0ab6971 100644
--- a/doc/source/cice_index.rst
+++ b/doc/source/cice_index.rst
@@ -314,10 +314,10 @@ either Celsius or Kelvin units).  Deprecated parameters are listed at the end.
    "highfreq", "high-frequency atmo coupling", "F"
    "hin_old", "ice thickness prior to growth/melt", "m"
    "hin_max", "category thickness limits", "m"
-   "hist_avg", "if true, write averaged data instead of snapshots", "T"
-   "histfreq", "units of history output frequency: y, m, w, d or 1", ""
+   "hist_avg", "if true, write averaged data instead of snapshots", "T,T,T,T,T"
+   "histfreq", "units of history output frequency: y, m, w, d or 1", "m,x,x,x,x"
    "histfreq_base", "reference date for history output", ""
-   "histfreq_n", "integer output frequency in histfreq units", ""
+   "histfreq_n", "integer output frequency in histfreq units", "1,1,1,1,1"
    "history_dir", "path to history output files", ""
    "history_file", "history output file prefix", ""
    "history_format", "history file format", ""
diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst
index 9906fba87..2a7240c78 100644
--- a/doc/source/user_guide/ug_case_settings.rst
+++ b/doc/source/user_guide/ug_case_settings.rst
@@ -176,7 +176,7 @@ setup_nml
    "", "zero", "restart output frequency relative to year-month-day of 0000-01-01", ""
    "``dumpfreq_n``", "integer", "write restart frequency with ``dumpfreq``", "1"
    "``dump_last``", "logical", "write restart on last time step of simulation", "``.false.``"
-   "``hist_avg``", "logical", "write time-averaged data", "``.true.``"
+   "``hist_avg``", "logical", "write time-averaged data", "``.true.,.true.,.true.,.true.,.true.``"
    "``histfreq``", "``d``", "write history every ``histfreq_n`` days", "'1','h','d','m','y'"
    "", "``h``", "write history every ``histfreq_n`` hours", ""
    "", "``m``", "write history every ``histfreq_n`` months", ""
diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst
index d9ea07a02..acc75b3d8 100644
--- a/doc/source/user_guide/ug_implementation.rst
+++ b/doc/source/user_guide/ug_implementation.rst
@@ -1154,7 +1154,8 @@ io package.  The namelist variable ``history_format`` further refines the
 format approach or style for some io packages.
 
 Model output data can be written as instantaneous or average data as specified
-by the ``hist_avg`` namelist flag.  The data is written at the period(s) given by ``histfreq`` and
+by the ``hist_avg`` namelist array and is customizable by stream.  The data is 
+written at the period(s) given by ``histfreq`` and
 ``histfreq_n`` relative to a reference date specified by ``histfreq_base``.  
 The files are written to binary or netCDF files prepended by ``history_file``
 in **ice_in**. These settings for history files are set in the 
@@ -1206,7 +1207,7 @@ For example, in the namelist:
   histfreq = ’1’, ’h’, ’d’, ’m’, ’y’
   histfreq_n = 1, 6, 0, 1, 1
   histfreq_base = 'zero'
-  hist_avg = .true.
+  hist_avg = .true.,.true.,.true.,.true.,.true.
   f_hi = ’1’
   f_hs = ’h’
   f_Tsfc = ’d’