Add functionality to change hist_avg for each stream (CICE-Consortium…

…#827)

* Add functionality to change hist_avg for each stream

* Fix some documentation

* Try to fix sphinx problem

* Fix hist_avg documentation

* Add some metadata changes to time and time_bounds

.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

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

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.', &

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

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

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' &

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' &

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

configuration/scripts/options/set_nml.histinst

@@ -1 +1 @@
-hist_avg        = .false.
+hist_avg        = .false.,.false.,.false.,.false.,.false.

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'

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'