Namelist option for time axis position. (CICE-Consortium#839)

* Add option to change location in interval of time axis

* Only use hist_time_axis when hist_avg is true

* Add more comments and information in the documentation

* Add a check on hist_time_axis as well as a global attribute

* Abort if hist_time_axis is not set correctly.

cicecore/cicedyn/analysis/ice_history_shared.F90

@@ -132,6 +132,8 @@ module ice_history_shared
                                        time_end(max_nstrm), &
                                        time_bounds(2)
 
+      character (len=char_len), public :: hist_time_axis
+
       real (kind=dbl_kind), allocatable, public :: &
          a2D (:,:,:,:)    , & ! field accumulations/averages, 2D
          a3Dz(:,:,:,:,:)  , & ! field accumulations/averages, 3D vertical

cicecore/cicedyn/general/ice_init.F90

@@ -81,7 +81,7 @@ subroutine input_data
           runid, runtype, use_restart_time, restart_format, lcdf64
       use ice_history_shared, only: hist_avg, history_dir, history_file, &
                              incond_dir, incond_file, version_name, &
-                             history_precision, history_format
+                             history_precision, history_format, hist_time_axis
       use ice_flux, only: update_ocn_f, l_mpond_fresh
       use ice_flux, only: default_season
       use ice_flux_bgc, only: cpl_bgc
@@ -185,6 +185,7 @@ subroutine input_data
         restart_ext,    use_restart_time, restart_format, lcdf64,       &
         pointer_file,   dumpfreq,       dumpfreq_n,      dump_last,     &
         diagfreq,       diag_type,      diag_file,       history_format,&
+        hist_time_axis,                                                 &
         print_global,   print_points,   latpnt,          lonpnt,        &
         debug_forcing,  histfreq,       histfreq_n,      hist_avg,      &
         history_dir,    history_file,   history_precision, cpl_bgc,     &
@@ -324,6 +325,8 @@ subroutine input_data
       histfreq_base = 'zero' ! output frequency reference date
       hist_avg(:) = .true.   ! if true, write time-averages (not snapshots)
       history_format = 'default' ! history file format
+      hist_time_axis = 'end' ! History file time axis averaging interval position
+
       history_dir  = './'    ! write to executable dir for default
       history_file = 'iceh'  ! history file name prefix
       history_precision = 4  ! precision of history files
@@ -906,6 +909,7 @@ subroutine input_data
       call broadcast_scalar(history_file,         master_task)
       call broadcast_scalar(history_precision,    master_task)
       call broadcast_scalar(history_format,       master_task)
+      call broadcast_scalar(hist_time_axis,       master_task)
       call broadcast_scalar(write_ic,             master_task)
       call broadcast_scalar(cpl_bgc,              master_task)
       call broadcast_scalar(incond_dir,           master_task)
@@ -1570,6 +1574,11 @@ subroutine input_data
          abort_list = trim(abort_list)//":24"
       endif
 
+      if(trim(hist_time_axis) /= 'begin' .and. trim(hist_time_axis) /= 'middle' .and. trim(hist_time_axis) /= 'end') then
+         write (nu_diag,*) subname//' ERROR: hist_time_axis value not valid = '//trim(hist_time_axis)
+         abort_list = trim(abort_list)//":29"
+      endif
+
       if(dumpfreq_base /= 'init' .and. dumpfreq_base /= 'zero') then
          write (nu_diag,*) subname//' ERROR: bad value for dumpfreq_base, allowed values: init, zero'
          abort_list = trim(abort_list)//":25"
@@ -2316,6 +2325,7 @@ subroutine input_data
          write(nu_diag,1031) ' history_file     = ', trim(history_file)
          write(nu_diag,1021) ' history_precision= ', history_precision
          write(nu_diag,1031) ' history_format   = ', trim(history_format)
+         write(nu_diag,1031) ' hist_time_axis   = ', trim(hist_time_axis)
          if (write_ic) then
             write(nu_diag,1039) ' Initial condition will be written in ', &
                                trim(incond_dir)

cicecore/cicedyn/infrastructure/io/io_netcdf/ice_history_write.F90

@@ -21,7 +21,7 @@
 
       module ice_history_write
 
-      use ice_constants, only: c0, c360, spval, spval_dbl
+      use ice_constants, only: c0, c360, p5, spval, spval_dbl
       use ice_fileunits, only: nu_diag
       use ice_exit, only: abort_ice
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
@@ -137,8 +137,6 @@ subroutine ice_write_hist (ns)
 
       if (my_task == master_task) then
 
-        ltime2 = timesecs/secday
-
         call construct_filename(ncfile(ns),'nc',ns)
 
         ! add local directory path name to ncfile
@@ -718,6 +716,12 @@ subroutine ice_write_hist (ns)
                          'ERROR: global attribute time_period_freq')
         endif
 
+        if (hist_avg(ns)) then
+           status = nf90_put_att(ncid,nf90_global,'time_axis_position',trim(hist_time_axis))
+           if (status /= nf90_noerr) call abort_ice(subname// &
+                         'ERROR: global attribute time axis position')
+        endif
+
         title = 'CF-1.0'
         status =  &
              nf90_put_att(ncid,nf90_global,'conventions',title)
@@ -749,7 +753,16 @@ subroutine ice_write_hist (ns)
       !-----------------------------------------------------------------
       ! write time variable
       !-----------------------------------------------------------------
-
+
+        ltime2 = timesecs/secday ! hist_time_axis = 'end' (default)
+
+        ! Some coupled models require the time axis "stamp" to be in the middle
+        ! or even beginning of averaging interval.
+        if (hist_avg(ns)) then
+           if (trim(hist_time_axis) == "begin" ) ltime2 = time_beg(ns)
+           if (trim(hist_time_axis) == "middle") ltime2 = p5*(time_beg(ns)+time_end(ns))
+        endif
+
         status = nf90_inq_varid(ncid,'time',varid)
         if (status /= nf90_noerr) call abort_ice(subname// &
                       'ERROR: getting time varid')

cicecore/cicedyn/infrastructure/io/io_pio2/ice_history_write.F90

@@ -18,7 +18,7 @@
       module ice_history_write
 
       use ice_kinds_mod
-      use ice_constants, only: c0, c360, spval, spval_dbl
+      use ice_constants, only: c0, c360, p5, spval, spval_dbl
       use ice_fileunits, only: nu_diag
       use ice_exit, only: abort_ice
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
@@ -185,8 +185,6 @@ subroutine ice_write_hist (ns)
       call ice_pio_initdecomp(ndim3=nzslyr,    ndim4=ncat_hist, iodesc=iodesc4ds, precision=history_precision)
       call ice_pio_initdecomp(ndim3=nfsd_hist, ndim4=ncat_hist, iodesc=iodesc4df, precision=history_precision)
 
-      ltime2 = timesecs/secday
-
       ! option of turning on double precision history files
       lprecision = pio_real
       if (history_precision == 8) lprecision = pio_double
@@ -678,6 +676,9 @@ subroutine ice_write_hist (ns)
            status = pio_put_att(File,pio_global,'time_period_freq',trim(time_period_freq))
         endif
 
+        if (hist_avg(ns)) &
+           status = pio_put_att(File,pio_global,'time_axis_position',trim(hist_time_axis))
+
         title = 'CF-1.0'
         status =  &
              pio_put_att(File,pio_global,'conventions',trim(title))
@@ -706,6 +707,15 @@ subroutine ice_write_hist (ns)
       ! write time variable
       !-----------------------------------------------------------------
 
+        ltime2 = timesecs/secday ! hist_time_axis = 'end' (default)
+
+        ! Some coupled models require the time axis "stamp" to be in the middle
+        ! or even beginning of averaging interval.
+        if (hist_avg(ns)) then
+           if (trim(hist_time_axis) == "begin" ) ltime2 = time_beg(ns)
+           if (trim(hist_time_axis) == "middle") ltime2 = p5*(time_beg(ns)+time_end(ns))
+        endif
+
         status = pio_inq_varid(File,'time',varid)
         status = pio_put_var(File,varid,(/1/),ltime2)
 

configuration/scripts/ice_in

@@ -53,6 +53,7 @@
     history_file   = 'iceh'
     history_precision = 4
     history_format = 'default'
+    hist_time_axis = 'end'
     write_ic       = .true.
     incond_dir     = './history/'
     incond_file    = 'iceh_ic'

doc/source/cice_index.rst

@@ -322,6 +322,7 @@ either Celsius or Kelvin units).  Deprecated parameters are listed at the end.
    "history_file", "history output file prefix", ""
    "history_format", "history file format", ""
    "history_precision", "history output precision: 4 or 8 byte", "4"
+   "hist_time_axis", "history file time axis interval location: begin, middle, end", "end"
    "hm", "land/boundary mask, thickness (T-cell)", ""
    "hmix", "ocean mixed layer depth", "20. m"
    "hour", "hour of the year", ""

doc/source/user_guide/ug_case_settings.rst

@@ -191,6 +191,7 @@ setup_nml
    "``history_format``", "``default``", "read/write history files in default format", "``default``"
    "", "``pio_pnetcdf``", "read/write restart files with pnetcdf in pio", ""
    "``history_precision``", "integer", "history file precision: 4 or 8 byte", "4"
+   "``hist_time_axis``","character","history file time axis interval location: begin, middle, end","end"
    "``ice_ic``", "``default``", "equal to internal", "``default``"
    "", "``internal``", "initial conditions set based on ice\_data\_type,conc,dist inputs", ""
    "", "``none``", "no ice", ""

doc/source/user_guide/ug_implementation.rst

@@ -1197,8 +1197,11 @@ with a given ``histfreq`` value, or if an element of ``histfreq_n`` is 0, then
 no file will be written at that frequency. The output period can be
 discerned from the filenames.  All history streams will be either instantaneous
 or averaged as specified by the ``hist_avg`` namelist setting and the frequency
-will be relative to a reference date specified by ``histfreq_base``.  More
-information about how the frequency is computed is found in :ref:`timemanager`.
+will be relative to a reference date specified by ``histfreq_base``.  Also, some
+Earth Sytem Models require the history file time axis to be centered in the averaging
+interval. The flag ``hist_time_axis`` will allow the user to chose ``begin``, ``middle``,
+or ``end`` for the time stamp. More information about how the frequency is 
+computed is found in :ref:`timemanager`.
 
 For example, in the namelist: