Fix time_average for instantaneous output

configure.ac

@@ -125,9 +125,15 @@ AM_CONDITIONAL([WITH_ASCIIDOC], [test -n "$ASCIIDOC"])
 
 AC_CHECK_PROGS(BATS, [bats])
 AC_CHECK_PROGS(PROVE, [prove])
-
 AM_CONDITIONAL([WITH_CHECK_PROGS], [test -n "$PROVE" -a -n "$BATS"])
 
+AC_PYTHON_MODULE(xarray, [], python3)
+if test ${HAVE_PYMOD_XARRAY} = yes; then
+  AM_CONDITIONAL([SKIP_XARRAY_TESTS], false )
+else
+  AM_CONDITIONAL([SKIP_XARRAY_TESTS], true )
+fi
+
 # Check if openacc.h exists
 if test "$enable_acc" = yes ; then
    AC_CHECK_HEADERS([openacc.h], [], [AC_MSG_ERROR(["Cannot find OpenACC header file"])] )

m4/ax_python_module.m4

@@ -0,0 +1,55 @@
+# ===========================================================================
+#     http://www.gnu.org/software/autoconf-archive/ax_python_module.html
+# ===========================================================================
+#
+# SYNOPSIS
+#
+#   AX_PYTHON_MODULE(modname[, fatal, python])
+#
+# DESCRIPTION
+#
+#   Checks for Python module.
+#
+#   If fatal is non-empty then absence of a module will trigger an error.
+#   The third parameter can either be "python" for Python 2 or "python3" for
+#   Python 3; defaults to Python 2.
+#
+# LICENSE
+#
+#   Copyright (c) 2008 Andrew Collier
+#
+#   Copying and distribution of this file, with or without modification, are
+#   permitted in any medium without royalty provided the copyright notice
+#   and this notice are preserved. This file is offered as-is, without any
+#   warranty.
+
+#serial 8
+
+AU_ALIAS([AC_PYTHON_MODULE], [AX_PYTHON_MODULE])
+AC_DEFUN([AX_PYTHON_MODULE],[
+
+    if test -z "$3";
+    then
+        PYTHON_TEST="python2"
+    else
+        PYTHON_TEST="$3"
+    fi
+
+    PYTHON_NAME=`basename $PYTHON_TEST`
+    AC_MSG_CHECKING($PYTHON_NAME module: $1)
+    $PYTHON_TEST -c "import $1" 2>/dev/null
+    if test $? -eq 0;
+    then
+        AC_MSG_RESULT(yes)
+        eval AS_TR_CPP(HAVE_PYMOD_$1)=yes
+    else
+        AC_MSG_RESULT(no)
+        eval AS_TR_CPP(HAVE_PYMOD_$1)=no
+        #
+        if test -n "$2"
+        then
+            AC_MSG_ERROR(failed to find required module $1)
+            exit 1
+        fi
+    fi
+])

postprocessing/timavg/time_average.f90

@@ -597,27 +597,26 @@ program time_average
      do itime = 0, numtime
 
       ! get time coordinate value
+        tavg = 1
         if (itime > 0) then
             istat =  NF90_GET_VAR (ncid_in, varid_recdim, time, start=(/itime/))
             if (istat /= NF90_NOERR) call error_handler ('getting time coord value', ncode=istat)
-
-            !--- read time bnds_info
-            istat = NF90_GET_VAR (ncid_in, time_bnds_id, tavg(1:2), (/1, itime/))
-            if (istat /= NF90_NOERR) call error_handler &
-              ('reading time bnds', ncode=istat)
-            tavg(3) = tavg(2) - tavg(1)
-        else
-            tavg = 1
-        endif
-
-        if (do_climo_avg .and. itime == numtime) then
-          ! use midpoint of last time interval for climatological time
-          if (change_bounds_to_climo) then
-            climo_time = 0.5*(tavg(1)+tavg(2))
-          else
-            ! use last time if already climo time
-            climo_time = time
-          endif
+            if (do_avg) then
+                !--- read time bnds_info
+                istat = NF90_GET_VAR (ncid_in, time_bnds_id, tavg(1:2), (/1, itime/))
+                if (istat /= NF90_NOERR) call error_handler &
+                   ('reading time bnds', ncode=istat)
+                tavg(3) = tavg(2) - tavg(1)
+                if (do_climo_avg .and. itime == numtime) then
+                    ! use midpoint of last time interval for climatological time
+                    if (change_bounds_to_climo) then
+                        climo_time = 0.5*(tavg(1)+tavg(2))
+                    else
+                        ! use last time if already climo time
+                        climo_time = time
+                    endif
+                endif
+           endif
         endif
 !-----------------------------------------------------------------------
 !-------------------- Loop through variables ---------------------------
@@ -1092,6 +1091,7 @@ subroutine check_for_climo_avg ( ncid, tname, tunits, ntimes,  &
 
  ! read the time bounds for the two time periods
 
+  tbnds_present = .false.
   if (tbnds_name(1:1) .ne. ' ') tbnds_present = .true.
   if (tbnds_present) then
       istat = NF90_INQ_VARID ( ncid, trim(tbnds_name), varid(1) )

site-configs/gfdl-ws/env.sh

@@ -42,9 +42,10 @@ mpi_version=4.1.1
 module remove-path MODULEPATH /app/spack/${env_version}/modulefiles/linux-rhel8-x86_64
 module prepend-path MODULEPATH /app/spack/${env_version}/modulefiles/linux-rhel8-x86_64
 
-# bats and nccmp are needed for tests
+# bats, nccmp, and python are needed for tests
 module load bats
 module load nccmp
+module load python
 
 # Load the GCC compilers
 module load gcc/$gcc_version

site-configs/gfdl/env.sh

@@ -42,9 +42,11 @@ mpi_version=4.1.1
 module remove-path MODULEPATH /app/spack/${env_version}/modulefiles/linux-rhel7-x86_64
 module prepend-path MODULEPATH /app/spack/${env_version}/modulefiles/linux-rhel7-x86_64
 
-# bats and nccmp are needed for tests
+# bats, nccmp, and python are needed for tests
 module load bats
 module load nccmp
+module load python
+
 # Need newer autoconf/automake than what pan has at the system level
 module load autoconf
 module load automake

site-configs/ncrc5/env.sh

@@ -32,6 +32,7 @@
 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 module rm PrgEnv-pgi PrgEnv-intel PrgEnv-gnu PrgEnv-cray
+module load python
 module load PrgEnv-gnu/8.3.3
 module load gcc/12.2.0
 module load cray-hdf5/1.12.2.3

site-configs/ncrc6/env.sh

@@ -32,6 +32,7 @@
 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 module rm PrgEnv-pgi PrgEnv-intel PrgEnv-gnu PrgEnv-cray
+module load python
 module load PrgEnv-gnu/8.5.0
 module load gcc-native/12.3
 module load cray-hdf5/1.12.2.11

t/Makefile.am

@@ -23,8 +23,15 @@ else
 skip_MPI="skip"
 endif
 
+if SKIP_XARRAY_TESTS
+skipflag="skip"
+else
+skipflag=""
+endif
+
 TESTS_ENVIRONMENT = top_srcdir=$(abs_top_srcdir); export top_srcdir;\
                     export skip_mpi=$(skip_MPI);                    \
+                    export skipflag=$(skipflag);                 \
     for d in $$( find $(abs_top_builddir)/postprocessing $(abs_top_builddir)/tools -mindepth 1 -not -path '*/\.*' -type d ); do \
       PATH="$$d"'$(PATH_SEPARATOR)'"$$PATH"; \
     done;
@@ -68,7 +75,8 @@ TESTS = Test01-check_programs_exist.sh \
         Test29-make_vgrid.sh \
         Test31-fregrid_stretched.sh \
 	Test32-fregrid_no_stretched.sh \
-	Test33-reference_make_hgrid.sh
+	Test33-reference_make_hgrid.sh \
+	Test34-timavg.sh
 
 EXTRA_DIST = $(TESTS) \
              Test02-input\

t/Test34-timavg.sh

@@ -0,0 +1,43 @@
+#!/usr/bin/env bats
+
+#***********************************************************************
+#                   GNU Lesser General Public License
+#
+# This file is part of the GFDL FRE NetCDF tools package (FRE-NCTools).
+#
+# FRE-NCTools is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or (at
+# your option) any later version.
+#
+# FRE-NCTools is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with FRE-NCTools.  If not, see
+# <http://www.gnu.org/licenses/>.
+#***********************************************************************
+
+load test_utils
+
+@test "Test timavg" {
+if [ -z "$skipflag" ]; then
+cat <<_EOF > instantaneous.nml
+ &input
+    file_names(1) =  'timavg_instantaneous_in.nc' ,
+    file_name_out =  'timavg_instantaneous_out.nc' ,
+ &end
+_EOF
+
+cat <<_EOF > standard.nml
+ &input
+    file_names(1) =  'timavg_standard_in.nc' ,
+    file_name_out =  'timavg_standard_out.nc' ,
+ &end
+_EOF
+
+    python3 $top_srcdir/t/test_time_avg.py
+fi
+}

t/test_time_avg.py

@@ -0,0 +1,101 @@
+#!/usr/bin/env python3
+
+import xarray as xr
+import numpy as np
+import os
+
+class TestTimeAvg:
+    def __init__(self, Test_case):
+        self.Test_case = Test_case
+        if self.Test_case == 1:
+            self.In_file = "timavg_instantaneous_in.nc"
+            self.Out_file = "timavg_instantaneous_out.nc"
+            self.Namelist_file = "instantaneous.nml"
+            self.Error_msg = "timavg with instantaneous input"
+        elif self.Test_case == 2:
+            self.In_file = "timavg_standard_in.nc"
+            self.Out_file = "timavg_standard_out.nc"
+            self.Namelist_file = "standard.nml"
+            self.Error_msg = "timavg with standard input"
+
+
+    def create_instantaneous(self):
+        Time_attrs = {"units": "days since 2000-01-01", "axis": "T", "calendar_type": "JULIAN", "calendar": "julian"}
+        ds1 = xr.Dataset(
+        data_vars={
+            ## This is backwards #FORTRAN4LYFE
+            "a": (("Time", "y", "x"), self.a),
+            "Time": (("Time", self.time_data, Time_attrs))
+        },
+        coords={
+            "x": np.arange(self.nx)*1.0,
+            "y": np.arange(self.ny)*1.0,
+        },
+        )
+        ds1.to_netcdf(self.In_file, unlimited_dims="Time")
+
+    def create_standard(self):
+        Time_attrs = {"units": "days since 2000-01-01", "axis": "T", "calendar_type": "JULIAN", "calendar": "julian",
+                      "bounds": "time_bnds"}
+        ds1 = xr.Dataset(
+        data_vars={
+            ## This is backwards #FORTRAN4LYFE
+            "a": (("Time", "y", "x"), self.a),
+            "Time": (("Time", self.time_data, Time_attrs)),
+            "time_bnds": (("Time", "nv"), self.time_bnds_data)
+        },
+        coords={
+            "x": np.arange(self.nx)*1.0,
+            "y": np.arange(self.ny)*1.0,
+            "nv": np.arange(2)*1.0+1,
+        },
+        )
+        ds1.to_netcdf(self.In_file, unlimited_dims="Time")
+
+
+    def create_input(self):
+        self.nx = 96
+        self.ny = 96
+        self.nt = 12
+        self.time_data = np.arange(self.nt)*1.0+1
+
+        self.a = np.zeros([self.nt, self.ny, self.nx])
+        for k in range(0,self.nt):
+            for j in range(0,self.ny):
+              for i in range(0,self.nx):
+                  self.a[k, j, i] = (i+1.)*1000 + (j+1.)*10 + (k+1.)/100.
+
+        if self.Test_case == 1:
+            self.create_instantaneous()
+        elif self.Test_case == 2:
+            self.time_bnds_data = np.zeros([self.nt, 2])
+            self.time_bnds_data[:,1] = self.time_data[0:]
+            self.time_bnds_data[1:,0] = self.time_data[0:-1]
+            self.create_standard()
+
+
+    def run_test(self):
+      exit_code = os.system("TAVG.exe < " + self.Namelist_file)
+      if exit_code != 0 :
+          raise Exception(self.Error_msg + ":: Failed to complete")
+
+
+    def check_output(self):
+      out_file = xr.open_dataset(self.Out_file)
+      if not (out_file.a.values == np.average(self.a, axis=0)).all():
+          raise Exception(self.Error_msg + ":: answers were not the expected result!")
+
+
+# Test time_average when the input is instantaneous (so no time bounds)
+Test_Instantaneous = 1
+test_class = TestTimeAvg(Test_Instantaneous)
+test_class.create_input()
+test_class.run_test()
+test_class.check_output()
+
+# Test time_average when the input is standard (with time_bounds, no average_* variables)
+Test_Standard = 2
+test_class = TestTimeAvg(Test_Standard)
+test_class.create_input()
+test_class.run_test()
+test_class.check_output()