diff --git a/pcmdi_metrics/variability_mode/lib/eof_analysis.py b/pcmdi_metrics/variability_mode/lib/eof_analysis.py
index d367716d3..e8864bbe3 100644
--- a/pcmdi_metrics/variability_mode/lib/eof_analysis.py
+++ b/pcmdi_metrics/variability_mode/lib/eof_analysis.py
@@ -298,13 +298,21 @@ def gain_pcs_fraction(full_field, eof_pattern, pcs, debug=False):
         print("from gain_pcs_fraction:")
         print("full_field.shape (before grower): ", full_field.shape)
         print("eof_pattern.shape (before grower): ", eof_pattern.shape)
+        print("pcs.shape: ", pcs.shape)
+        print("pcs[0:5]:", pcs[0:5])
+        print("full_field: max, min:", np.max(full_field), np.min(full_field))
+        print("eof_pattern: max, min:", np.max(eof_pattern), np.min(eof_pattern))
+        print("pcs: max, min:", np.max(pcs), np.min(pcs))
+
     # Extend eof_pattern (add 3rd dimension as time then copy same 2d value for all time step)
-    reconstructed_field = genutil.grower(full_field, eof_pattern)[
-        1
-    ]  # Matching dimension (add time axis)
+    # Matching dimension (add time axis)
+    reconstructed_field = genutil.grower(full_field, eof_pattern)[1]
+
+    # Reconstruct field
     for t in range(0, len(pcs)):
         reconstructed_field[t] = MV2.multiply(reconstructed_field[t], pcs[t])
-    # 2-2) Get variance of reconstructed field
+
+    # 2-2) Get variance of the reconstructed field
     variance_partial = genutil.statistics.variance(reconstructed_field, axis="t")
     variance_partial_area_ave = cdutil.averager(
         variance_partial, axis="xy", weights="weighted"
@@ -313,6 +321,17 @@ def gain_pcs_fraction(full_field, eof_pattern, pcs, debug=False):
     fraction = MV2.divide(variance_partial_area_ave, variance_total_area_ave)
     # debugging
     if debug:
+        print(
+            "reconstructed_field: max, min:",
+            np.max(reconstructed_field),
+            np.min(reconstructed_field),
+        )
+        print(
+            "variance_partial: max, min:",
+            np.max(variance_partial),
+            np.min(variance_partial),
+        )
+
         print("full_field.shape (after grower): ", full_field.shape)
         print("reconstructed_field.shape: ", reconstructed_field.shape)
         print("variance_partial_area_ave: ", variance_partial_area_ave)
@@ -381,9 +400,16 @@ def get_anomaly_timeseries(timeseries, season):
         if season != "monthly":
             # Get seasonal mean time series
             # each season chunk should have 100% of data
-            timeseries_ano = getattr(cdutil, season)(
-                timeseries_ano, criteriaarg=[1.0, None]
-            )
+            if season.upper() in ["DJF", "MAM", "JJA", "SON"]:
+                timeseries_ano = getattr(cdutil, season)(
+                    timeseries_ano, criteriaarg=[1.0, None]
+                )
+            else:
+                custom_season_accessor = cdutil.times.Seasons(season.upper())
+                timeseries_ano = custom_season_accessor(
+                    timeseries_ano, criteriaarg=[1.0, None]
+                )
+
     # return result
     return timeseries_ano