diff --git a/tests/integration/generated/test_debug2_chrysalis.cfg b/tests/integration/generated/test_debug2_chrysalis.cfg
new file mode 100644
index 00000000..4ea749ce
--- /dev/null
+++ b/tests/integration/generated/test_debug2_chrysalis.cfg
@@ -0,0 +1,58 @@
+# Use this cfg to test a pull request before merging.
+# 1. Copy over the contents of template_complete_run.cfg.
+# 2. Delete parts that aren't needed. `output` and `www` should include `zppy_test_debug2` rather than `zppy_test_complete_run`.
+# 3. Iteratively debug/develop with this cfg.
+# 4. Once it's working as expected, copy important changes to template_complete_run.cfg.
+
+[default]
+case = v2.LR.historical_0201
+constraint = ""
+dry_run = "False"
+environment_commands = ""
+input = "/lcrc/group/e3sm/ac.forsyth2//E3SMv2/v2.LR.historical_0201"
+input_subdir = archive/atm/hist
+mapping_file = "map_ne30pg2_to_cmip6_180x360_aave.20200201.nc"
+# To run this test, edit `output` and `www` in this file, along with `actual_images_dir` in test_debug2.py
+output = "/lcrc/group/e3sm/ac.forsyth2/zppy_test_debug2_output/unique_id/v2.LR.historical_0201"
+partition = "debug"
+qos = "regular"
+www = "/lcrc/group/e3sm/public_html/diagnostic_output/ac.forsyth2/zppy_test_debug2_www/unique_id"
+
+[ts]
+active = True
+e3sm_to_cmip_environment_commands = ""
+environment_commands = "source /home/ac.forsyth2/miniconda3/etc/profile.d/conda.sh; conda activate zppy_dev_with_nco_20240405"
+walltime = "00:30:00"
+years = "1850:1854:2",
+
+  [[ atm_monthly_glb ]]
+  # Note global average won't work for 3D variables.
+  frequency = "monthly"
+  input_files = "eam.h0"
+  input_subdir = "archive/atm/hist"
+  mapping_file = "glb"
+  years = "1850:1860:5",
+
+[mpas_analysis]
+active = True
+anomalyRefYear = 1850
+climo_years ="1850-1854", "1855-1860",
+enso_years = "1850-1854", "1855-1860",
+mesh = "EC30to60E2r2"
+parallelTaskCount = 6
+partition = "compute"
+qos = "regular"
+ts_years = "1850-1854", "1850-1860",
+walltime = "00:30:00"
+
+[global_time_series]
+active = True
+climo_years ="1850-1854", "1855-1860",
+experiment_name = "v2.LR.historical_0201"
+figstr = "v2_historical_0201"
+moc_file=mocTimeSeries_1850-1860.nc
+regions="glb"
+ts_num_years = 5
+ts_years = "1850-1854", "1850-1860",
+walltime = "00:30:00"
+years = "1850-1860",
diff --git a/tests/integration/generated/test_debug3_chrysalis.cfg b/tests/integration/generated/test_debug3_chrysalis.cfg
new file mode 100644
index 00000000..16e92b9a
--- /dev/null
+++ b/tests/integration/generated/test_debug3_chrysalis.cfg
@@ -0,0 +1,57 @@
+# Use this cfg to test a pull request before merging.
+# 1. Copy over the contents of template_complete_run.cfg.
+# 2. Delete parts that aren't needed. `output` and `www` should include `zppy_test_debug3` rather than `zppy_test_complete_run`.
+# 3. Iteratively debug/develop with this cfg.
+# 4. Once it's working as expected, copy important changes to template_complete_run.cfg.
+
+[default]
+case = v2.LR.BGC-LNDATM.CONTRL
+constraint = ""
+dry_run = "False"
+environment_commands = ""
+#input = "/lcrc/group/e3sm/ac.forsyth2//E3SMv2/v2.LR.historical_0201"
+input = "/lcrc/group/e3sm/ac.sfeng1/E3SM_Simulations/20221127.v2.LR.BGC-LNDATM.CONTRL.ne30pg2_r05_EC30to60E2r2.chrysalis"
+input_subdir = archive/atm/hist
+mapping_file = "map_r05_to_cmip6_180x360_aave.20231110.nc"
+# To run this test, edit `output` and `www` in this file, along with `actual_images_dir` in test_complete_run.py
+output = "/lcrc/group/e3sm/ac.forsyth2/zppy_test_debug3_output/unique_id/20221127.v2.LR.BGC-LNDATM.CONTRL.ne30pg2_r05_EC30to60E2r2.chrysalis"
+partition = "debug"
+qos = "regular"
+#vars = "TREFHT,FSNTOA,FLUT,PRECC,PRECL,QFLX,TS,FSNT,FLNT,TOTSOMC,TOTECOSYSC,NEE,TOTVEGC"
+www = "/lcrc/group/e3sm/public_html/diagnostic_output/ac.forsyth2/zppy_test_debug3_www/unique_id"
+
+[ts]
+active = True
+e3sm_to_cmip_environment_commands = ""
+environment_commands = "source /home/ac.forsyth2/miniconda3/etc/profile.d/conda.sh; conda activate zppy_dev_with_nco_20240405"
+walltime = "00:30:00"
+
+  [[ atm_monthly_glb ]]
+  # Note global average won't work for 3D variables.
+  frequency = "monthly"
+  input_files = "eam.h0"
+  input_subdir = "archive/atm/hist"
+  mapping_file = "glb"
+  # First 6 are for the original plots.
+  vars = "TREFHT,FSNTOA,FLUT,PRECC,PRECL,QFLX,TS,FSNT,FLNT"
+  years = "1880:1890:5",
+
+  [[ lnd_monthly_glb ]]
+  frequency = "monthly"
+  input_files = "elm.h0"
+  input_subdir = "archive/lnd/hist"
+  mapping_file = "glb"
+  vars = "TOTSOMC,TOTECOSYSC,NBP,TOTVEGC"
+  years = "1880:1890:5",
+
+[global_time_series]
+active = True
+atmosphere_only = "True"
+experiment_name = "v2.LR.BGC-LNDATM.CONTRL"
+figstr = "v2.LR.BGC-LNDATM.CONTRL"
+plots_original = "net_toa_flux_restom,global_surface_air_temperature,toa_radiation,net_atm_energy_imbalance,change_ohc,max_moc,change_sea_level,net_atm_water_imbalance"
+plots_atm = "RESSURF,TS,FSNT,FLNT"
+plots_lnd = "TOTSOMC,TOTECOSYSC,NBP,TOTVEGC"
+ts_num_years = 5
+walltime = "00:30:00"
+years = "1880-1890",
diff --git a/tests/integration/post.v3.LR.piControl.cfg b/tests/integration/post.v3.LR.piControl.cfg
new file mode 100644
index 00000000..5167f450
--- /dev/null
+++ b/tests/integration/post.v3.LR.piControl.cfg
@@ -0,0 +1,41 @@
+[default]
+input = /lcrc/group/e3sm2/ac.golaz/E3SMv3/v3.LR.piControl
+output = /lcrc/group/e3sm/ac.forsyth2/zppy_test_v3_output/v3/v3.LR.piControl
+case = v3.LR.piControl
+www = /lcrc/group/e3sm/public_html/diagnostic_output/ac.forsyth2/zppy_test_v3_www/v3
+partition = compute
+environment_commands = "source /lcrc/soft/climate/e3sm-unified/load_latest_e3sm_unified_chrysalis.sh"
+campaign = "water_cycle"
+
+[ts]
+active = True
+environment_commands = "source /home/ac.forsyth2/miniconda3/etc/profile.d/conda.sh; conda activate zppy_dev_with_nco_20240405"
+years = "0001:0050:10",
+walltime = "00:50:00"
+
+  [[ atm_monthly_glb ]]
+  input_subdir = "archive/atm/hist"
+  input_files = "eam.h0"
+  frequency = "monthly"
+  mapping_file = "glb"
+
+  [[ lnd_monthly_glb ]]
+  input_subdir = "archive/lnd/hist"
+  input_files = "elm.h0"
+  frequency = "monthly"
+  mapping_file = "glb"
+  vars = "FSH,RH2M,LAISHA,LAISUN,QINTR,QOVER,QRUNOFF,QSOIL,QVEGE,QVEGT,SOILWATER_10CM,TSA,H2OSNO,TOTLITC,CWDC,SOIL1C,SOIL2C,SOIL3C,SOIL4C,WOOD_HARVESTC,TOTVEGC,NBP,GPP,AR,HR"
+
+
+[global_time_series]
+active = True
+experiment_name = "v3.LR.piControl"
+figstr = "v3.LR.piControl"
+atmosphere_only = True
+#plots_original = "net_toa_flux_restom,global_surface_air_temperature,toa_radiation,net_atm_energy_imbalance,change_ohc,max_moc,change_sea_level,net_atm_water_imbalance"
+plots_lnd = "FSH,RH2M,LAISHA,LAISUN,QINTR,QOVER,QRUNOFF,QSOIL,QVEGE,QVEGT,SOILWATER_10CM,TSA,H2OSNO,TOTLITC,CWDC,SOIL1C,SOIL2C,SOIL3C,SOIL4C,WOOD_HARVESTC,TOTVEGC,NBP,GPP,AR,HR"
+ts_num_years = 10
+walltime = "00:30:00"
+years = "1-50",
+climo_years ="1-50",
+ts_years ="1-50",
diff --git a/tests/integration/template_debug2.cfg b/tests/integration/template_debug2.cfg
new file mode 100644
index 00000000..e5c2cd14
--- /dev/null
+++ b/tests/integration/template_debug2.cfg
@@ -0,0 +1,58 @@
+# Use this cfg to test a pull request before merging.
+# 1. Copy over the contents of template_complete_run.cfg.
+# 2. Delete parts that aren't needed. `output` and `www` should include `zppy_test_debug2` rather than `zppy_test_complete_run`.
+# 3. Iteratively debug/develop with this cfg.
+# 4. Once it's working as expected, copy important changes to template_complete_run.cfg.
+
+[default]
+case = v2.LR.historical_0201
+constraint = "#expand constraint#"
+dry_run = "#expand dry_run#"
+environment_commands = "#expand environment_commands#"
+input = "#expand user_input#/E3SMv2/v2.LR.historical_0201"
+input_subdir = archive/atm/hist
+mapping_file = "map_ne30pg2_to_cmip6_180x360_aave.20200201.nc"
+# To run this test, edit `output` and `www` in this file, along with `actual_images_dir` in test_debug2.py
+output = "#expand user_output#zppy_test_debug2_output/#expand unique_id#/v2.LR.historical_0201"
+partition = "#expand partition_short#"
+qos = "#expand qos_short#"
+www = "#expand user_www#zppy_test_debug2_www/#expand unique_id#"
+
+[ts]
+active = True
+e3sm_to_cmip_environment_commands = "#expand e3sm_to_cmip_environment_commands#"
+environment_commands = "source /home/ac.forsyth2/miniconda3/etc/profile.d/conda.sh; conda activate zppy_dev_with_nco_20240405"
+walltime = "00:30:00"
+years = "1850:1854:2",
+
+  [[ atm_monthly_glb ]]
+  # Note global average won't work for 3D variables.
+  frequency = "monthly"
+  input_files = "eam.h0"
+  input_subdir = "archive/atm/hist"
+  mapping_file = "glb"
+  years = "1850:1860:5",
+
+[mpas_analysis]
+active = True
+anomalyRefYear = 1850
+climo_years ="1850-1854", "1855-1860",
+enso_years = "1850-1854", "1855-1860",
+mesh = "EC30to60E2r2"
+parallelTaskCount = 6
+partition = "#expand partition_long#"
+qos = "#expand qos_long#"
+ts_years = "1850-1854", "1850-1860",
+walltime = "#expand mpas_analysis_walltime#"
+
+[global_time_series]
+active = True
+climo_years ="1850-1854", "1855-1860",
+experiment_name = "v2.LR.historical_0201"
+figstr = "v2_historical_0201"
+moc_file=mocTimeSeries_1850-1860.nc
+regions="glb"
+ts_num_years = 5
+ts_years = "1850-1854", "1850-1860",
+walltime = "00:30:00"
+years = "1850-1860",
diff --git a/tests/integration/template_debug3.cfg b/tests/integration/template_debug3.cfg
new file mode 100644
index 00000000..9defc3e9
--- /dev/null
+++ b/tests/integration/template_debug3.cfg
@@ -0,0 +1,57 @@
+# Use this cfg to test a pull request before merging.
+# 1. Copy over the contents of template_complete_run.cfg.
+# 2. Delete parts that aren't needed. `output` and `www` should include `zppy_test_debug3` rather than `zppy_test_complete_run`.
+# 3. Iteratively debug/develop with this cfg.
+# 4. Once it's working as expected, copy important changes to template_complete_run.cfg.
+
+[default]
+case = v2.LR.BGC-LNDATM.CONTRL
+constraint = "#expand constraint#"
+dry_run = "#expand dry_run#"
+environment_commands = "#expand environment_commands#"
+#input = "#expand user_input#/E3SMv2/v2.LR.historical_0201"
+input = "/lcrc/group/e3sm/ac.sfeng1/E3SM_Simulations/20221127.v2.LR.BGC-LNDATM.CONTRL.ne30pg2_r05_EC30to60E2r2.chrysalis"
+input_subdir = archive/atm/hist
+mapping_file = "map_r05_to_cmip6_180x360_aave.20231110.nc"
+# To run this test, edit `output` and `www` in this file, along with `actual_images_dir` in test_complete_run.py
+output = "#expand user_output#zppy_test_debug3_output/#expand unique_id#/20221127.v2.LR.BGC-LNDATM.CONTRL.ne30pg2_r05_EC30to60E2r2.chrysalis"
+partition = "#expand partition_short#"
+qos = "#expand qos_short#"
+#vars = "TREFHT,FSNTOA,FLUT,PRECC,PRECL,QFLX,TS,FSNT,FLNT,TOTSOMC,TOTECOSYSC,NEE,TOTVEGC"
+www = "#expand user_www#zppy_test_debug3_www/#expand unique_id#"
+
+[ts]
+active = True
+e3sm_to_cmip_environment_commands = "#expand e3sm_to_cmip_environment_commands#"
+environment_commands = "source /home/ac.forsyth2/miniconda3/etc/profile.d/conda.sh; conda activate zppy_dev_with_nco_20240405"
+walltime = "00:30:00"
+
+  [[ atm_monthly_glb ]]
+  # Note global average won't work for 3D variables.
+  frequency = "monthly"
+  input_files = "eam.h0"
+  input_subdir = "archive/atm/hist"
+  mapping_file = "glb"
+  # First 6 are for the original plots.
+  vars = "TREFHT,FSNTOA,FLUT,PRECC,PRECL,QFLX,TS,FSNT,FLNT"
+  years = "1880:1890:5",
+
+  [[ lnd_monthly_glb ]]
+  frequency = "monthly"
+  input_files = "elm.h0"
+  input_subdir = "archive/lnd/hist"
+  mapping_file = "glb"
+  vars = "TOTSOMC,TOTECOSYSC,NBP,TOTVEGC"
+  years = "1880:1890:5",
+
+[global_time_series]
+active = True
+atmosphere_only = "True"
+experiment_name = "v2.LR.BGC-LNDATM.CONTRL"
+figstr = "v2.LR.BGC-LNDATM.CONTRL"
+plots_original = "net_toa_flux_restom,global_surface_air_temperature,toa_radiation,net_atm_energy_imbalance,change_ohc,max_moc,change_sea_level,net_atm_water_imbalance"
+plots_atm = "RESSURF,TS,FSNT,FLNT"
+plots_lnd = "TOTSOMC,TOTECOSYSC,NBP,TOTVEGC"
+ts_num_years = 5
+walltime = "00:30:00"
+years = "1880-1890",
diff --git a/tests/integration/utils.py b/tests/integration/utils.py
index 581d978c..635136dc 100644
--- a/tests/integration/utils.py
+++ b/tests/integration/utils.py
@@ -260,7 +260,7 @@ def generate_cfgs(unified_testing=False, dry_run=False):
     else:
         expansions["dry_run"] = "False"
 
-    cfg_names = ["bundles", "complete_run", "debug"]
+    cfg_names = ["bundles", "complete_run", "debug", "debug2", "debug3"]
     for cfg_name in cfg_names:
         cfg_template = f"{git_top_level}/tests/integration/template_{cfg_name}.cfg"
         cfg_generated = (
diff --git a/zppy/global_time_series.py b/zppy/global_time_series.py
index c46fe992..e3cb55a6 100644
--- a/zppy/global_time_series.py
+++ b/zppy/global_time_series.py
@@ -15,7 +15,8 @@
 
 
 # -----------------------------------------------------------------------------
-def global_time_series(config, scriptDir, existing_bundles, job_ids_file):
+# FIXME: C901 'run' is too complex (19)
+def global_time_series(config, scriptDir, existing_bundles, job_ids_file):  # noqa: C901
 
     # Initialize jinja2 template engine
     templateLoader = jinja2.FileSystemLoader(
@@ -52,6 +53,59 @@ def global_time_series(config, scriptDir, existing_bundles, job_ids_file):
             if skip:
                 continue
 
+            # Handle legacy parameter
+            if c["plot_names"]:
+                print("warning: plot_names for global_time_series is deprecated.")
+                print(
+                    "Setting plot_names will override the new parameter, plots_original."
+                )
+                c["plots_original"] = c["plot_names"]
+
+            # Determine which components are needed
+            c["use_atm"] = False
+            c["use_ice"] = False
+            c["use_lnd"] = False
+            c["use_ocn"] = False
+            if c["plots_original"]:
+                c["use_atm"] = True
+                if c["atmosphere_only"]:
+                    print(
+                        "warning: atmosphere_only for global_time_series is deprecated."
+                    )
+                    print(
+                        "preferred method: remove the 3 ocean plots (change_ohc,max_moc,change_sea_level) from plots_original."
+                    )
+                has_original_ocn_plots = (
+                    ("change_ohc" in c["plots_original"])
+                    or ("max_moc" in c["plots_original"])
+                    or ("change_sea_level" in c["plots_original"])
+                )
+                if (not c["atmosphere_only"]) and has_original_ocn_plots:
+                    c["use_ocn"] = True
+            else:
+                # For better string processing in global_time_series.bash
+                c["plots_original"] = "None"
+            if c["plots_atm"]:
+                c["use_atm"] = True
+            else:
+                # For better string processing in global_time_series.bash
+                c["plots_atm"] = "None"
+            if c["plots_ice"]:
+                c["use_ice"] = True
+            else:
+                # For better string processing in global_time_series.bash
+                c["plots_ice"] = "None"
+            if c["plots_lnd"]:
+                c["use_lnd"] = True
+            else:
+                # For better string processing in global_time_series.bash
+                c["plots_lnd"] = "None"
+            if c["plots_ocn"]:
+                c["use_ocn"] = True
+            else:
+                # For better string processing in global_time_series.bash
+                c["plots_ocn"] = "None"
+
             # Load useful scripts
             c["global_time_series_dir"] = os.path.join(
                 scriptDir, "{}_dir".format(prefix)
@@ -74,21 +128,37 @@ def global_time_series(config, scriptDir, existing_bundles, job_ids_file):
             # List of dependencies
             dependencies = []
             # Add Time Series dependencies
-            # Iterate from year1 to year2 incrementing by the number of years per time series file.
-            for yr in range(c["year1"], c["year2"], c["ts_num_years"]):
-                start_yr = yr
-                end_yr = yr + c["ts_num_years"] - 1
-                dependencies.append(
-                    os.path.join(
-                        scriptDir,
-                        "ts_%s_%04d-%04d-%04d.status"
-                        % ("atm_monthly_glb", start_yr, end_yr, c["ts_num_years"]),
+            if c["use_atm"]:
+                # Iterate from year1 to year2 incrementing by the number of years per time series file.
+                for yr in range(c["year1"], c["year2"], c["ts_num_years"]):
+                    start_yr = yr
+                    end_yr = yr + c["ts_num_years"] - 1
+                    dependencies.append(
+                        os.path.join(
+                            scriptDir,
+                            "ts_%s_%04d-%04d-%04d.status"
+                            % ("atm_monthly_glb", start_yr, end_yr, c["ts_num_years"]),
+                        )
                     )
-                )
-            if not c["atmosphere_only"]:
+            if c["use_lnd"]:
+                for yr in range(c["year1"], c["year2"], c["ts_num_years"]):
+                    start_yr = yr
+                    end_yr = yr + c["ts_num_years"] - 1
+                    dependencies.append(
+                        os.path.join(
+                            scriptDir,
+                            "ts_%s_%04d-%04d-%04d.status"
+                            % ("lnd_monthly_glb", start_yr, end_yr, c["ts_num_years"]),
+                        )
+                    )
+            if c["use_ocn"]:
                 # Add MPAS Analysis dependencies
                 ts_year_sets = getYears(c["ts_years"])
                 climo_year_sets = getYears(c["climo_years"])
+                if (not ts_year_sets) or (not climo_year_sets):
+                    raise Exception(
+                        "ts_years and climo_years must both be set for ocn plots."
+                    )
                 for ts_year_set, climo_year_set in zip(ts_year_sets, climo_year_sets):
                     c["ts_year1"] = ts_year_set[0]
                     c["ts_year2"] = ts_year_set[1]
diff --git a/zppy/templates/coupled_global.py b/zppy/templates/coupled_global.py
index 6fcaba47..4957ba55 100644
--- a/zppy/templates/coupled_global.py
+++ b/zppy/templates/coupled_global.py
@@ -2,6 +2,7 @@
 import glob
 import math
 import sys
+import traceback
 from typing import Any, List, Tuple
 
 import matplotlib as mpl
@@ -55,7 +56,7 @@ def add_line(year, var, year1, year2, ax, format="%4.2f", lw=1, color="b"):
     i2 = (np.abs(year - year2)).argmin()
 
     tmp = np.average(var[i1 : i2 + 1])
-    ax.plot((year[i1], year[i2]), (tmp, tmp), lw=lw, color=color)
+    ax.plot((year[i1], year[i2]), (tmp, tmp), lw=lw, color=color, label="average")
     ax.text(ax.get_xlim()[1] + 1, tmp, format % tmp, va="center", color=color)
 
     return
@@ -86,7 +87,7 @@ def add_trend(
     if verbose:
         print(fit)
     fit_fn = np.poly1d(fit)
-    ax.plot(x, fit_fn(x), lw=lw, ls="--", c=color)
+    ax.plot(x, fit_fn(x), lw=lw, ls="--", c=color, label="trend")
     if ohc:
         # Earth radius 6371229. from MPAS-O output files
         heat_uptake = fit[0] / (4.0 * math.pi * (6371229.0) ** 2 * 365.0 * 86400.0)
@@ -112,20 +113,41 @@ def add_trend(
 # -----------------------------------------------------------------------------
 # Function to get ylim
 def get_ylim(standard_range, extreme_values):
-    standard_min = standard_range[0]
-    standard_max = standard_range[1]
-    if extreme_values == []:
-        return [standard_min, standard_max]
-    extreme_min = np.amin(extreme_values)
-    extreme_max = np.amax(extreme_values)
-    if standard_min <= extreme_min:
-        ylim_min = standard_min
+    if len(extreme_values) > 0:
+        has_extreme_values = True
+        extreme_min = np.amin(extreme_values)
+        extreme_max = np.amax(extreme_values)
+    else:
+        has_extreme_values = False
+        extreme_min = None
+        extreme_max = None
+    if len(standard_range) == 2:
+        has_standard_range = True
+        standard_min = standard_range[0]
+        standard_max = standard_range[1]
     else:
+        has_standard_range = False
+        standard_min = None
+        standard_max = None
+    if has_extreme_values and has_standard_range:
+        # Use at least the standard range,
+        # perhaps a wider window to include extremes
+        if standard_min <= extreme_min:
+            ylim_min = standard_min
+        else:
+            ylim_min = extreme_min
+        if standard_max >= extreme_max:
+            ylim_max = standard_max
+        else:
+            ylim_max = extreme_max
+    elif has_extreme_values and not has_standard_range:
         ylim_min = extreme_min
-    if standard_max >= extreme_max:
+        ylim_max = extreme_max
+    elif has_standard_range and not has_extreme_values:
+        ylim_min = standard_min
         ylim_max = standard_max
     else:
-        ylim_max = extreme_max
+        raise ValueError("Not enough range information supplied")
     return [ylim_min, ylim_max]
 
 
@@ -134,6 +156,7 @@ def get_ylim(standard_range, extreme_values):
 
 # 1
 def plot_net_toa_flux_restom(ax, xlim, exps, rgn):
+    print("Plot 1: plot_net_toa_flux_restom")
     param_dict = {
         "2nd_var": False,
         "axhline_y": 0,
@@ -152,7 +175,7 @@ def plot_net_toa_flux_restom(ax, xlim, exps, rgn):
         "shorten_year": False,
         "title": "Net TOA flux (restom)",
         "use_getmoc": False,
-        "var": lambda exp: np.array(exp["annual"]["RESTOM"]),
+        "var": lambda exp: np.array(exp["annual"]["RESTOM"][0]),
         "verbose": False,
         "vol": False,
         "ylabel": "W m-2",
@@ -162,6 +185,7 @@ def plot_net_toa_flux_restom(ax, xlim, exps, rgn):
 
 # 2
 def plot_global_surface_air_temperature(ax, xlim, exps, rgn):
+    print("Plot 2: plot_global_surface_air_temperature")
     if rgn == "glb":
         region_title = "Global"
     elif rgn == "n":
@@ -184,11 +208,11 @@ def plot_global_surface_air_temperature(ax, xlim, exps, rgn):
         "lw": 1.0,
         "ohc": False,
         "set_axhline": False,
-        "set_legend": False,
+        "set_legend": True,
         "shorten_year": False,
         "title": f"{region_title} surface air temperature",
         "use_getmoc": False,
-        "var": lambda exp: np.array(exp["annual"]["TREFHT"]) - 273.15,
+        "var": lambda exp: np.array(exp["annual"]["TREFHT"][0]) - 273.15,
         "verbose": False,
         "vol": False,
         "ylabel": "degC",
@@ -198,6 +222,7 @@ def plot_global_surface_air_temperature(ax, xlim, exps, rgn):
 
 # 3
 def plot_toa_radiation(ax, xlim, exps, rgn):
+    print("Plot 3: plot_toa_radiation")
     param_dict = {
         "2nd_var": True,
         "axhline_y": None,
@@ -216,7 +241,7 @@ def plot_toa_radiation(ax, xlim, exps, rgn):
         "shorten_year": False,
         "title": "TOA radiation: SW (solid), LW (dashed)",
         "use_getmoc": False,
-        "var": lambda exp: np.array(exp["annual"]["FSNTOA"]),
+        "var": lambda exp: np.array(exp["annual"]["FSNTOA"][0]),
         "verbose": None,
         "vol": None,
         "ylabel": "W m-2",
@@ -226,6 +251,7 @@ def plot_toa_radiation(ax, xlim, exps, rgn):
 
 # 4
 def plot_net_atm_energy_imbalance(ax, xlim, exps, rgn):
+    print("Plot 4: plot_net_atm_energy_imbalance")
     param_dict = {
         "2nd_var": False,
         "axhline_y": None,
@@ -240,12 +266,12 @@ def plot_net_atm_energy_imbalance(ax, xlim, exps, rgn):
         "lw": 1.0,
         "ohc": False,
         "set_axhline": False,
-        "set_legend": False,
+        "set_legend": True,
         "shorten_year": False,
         "title": "Net atm energy imbalance (restom-ressurf)",
         "use_getmoc": False,
-        "var": lambda exp: np.array(exp["annual"]["RESTOM"])
-        - np.array(exp["annual"]["RESSURF"]),
+        "var": lambda exp: np.array(exp["annual"]["RESTOM"][0])
+        - np.array(exp["annual"]["RESSURF"][0]),
         "verbose": False,
         "vol": False,
         "ylabel": "W m-2",
@@ -255,6 +281,7 @@ def plot_net_atm_energy_imbalance(ax, xlim, exps, rgn):
 
 # 5
 def plot_change_ohc(ax, xlim, exps, rgn):
+    print("Plot 5: plot_change_ohc")
     param_dict = {
         "2nd_var": False,
         "axhline_y": 0,
@@ -283,6 +310,7 @@ def plot_change_ohc(ax, xlim, exps, rgn):
 
 # 6
 def plot_max_moc(ax, xlim, exps, rgn):
+    print("Plot 6: plot_max_moc")
     param_dict = {
         "2nd_var": False,
         "axhline_y": 10,
@@ -311,6 +339,7 @@ def plot_max_moc(ax, xlim, exps, rgn):
 
 # 7
 def plot_change_sea_level(ax, xlim, exps, rgn):
+    print("Plot 7: plot_change_sea_level")
     param_dict = {
         "2nd_var": False,
         "axhline_y": None,
@@ -343,6 +372,7 @@ def plot_change_sea_level(ax, xlim, exps, rgn):
 
 # 8
 def plot_net_atm_water_imbalance(ax, xlim, exps, rgn):
+    print("Plot 8: plot_net_atm_water_imbalance")
     param_dict = {
         "2nd_var": False,
         "axhline_y": None,
@@ -357,7 +387,7 @@ def plot_net_atm_water_imbalance(ax, xlim, exps, rgn):
         "lw": 1.0,
         "ohc": False,
         "set_axhline": False,
-        "set_legend": False,
+        "set_legend": True,
         "shorten_year": False,
         "title": "Net atm water imbalance (evap-prec)",
         "use_getmoc": False,
@@ -365,9 +395,12 @@ def plot_net_atm_water_imbalance(ax, xlim, exps, rgn):
             365
             * 86400
             * (
-                np.array(exp["annual"]["QFLX"])
+                np.array(exp["annual"]["QFLX"][0])
                 - 1e3
-                * (np.array(exp["annual"]["PRECC"]) + np.array(exp["annual"]["PRECL"]))
+                * (
+                    np.array(exp["annual"]["PRECC"][0])
+                    + np.array(exp["annual"]["PRECL"][0])
+                )
             )
         ),
         "verbose": False,
@@ -377,16 +410,52 @@ def plot_net_atm_water_imbalance(ax, xlim, exps, rgn):
     plot(ax, xlim, exps, param_dict, rgn)
 
 
-def plot(ax, xlim, exps, param_dict, rgn):
+# Generic plot function
+def plot_generic(ax, xlim, exps, var_name, rgn):
+    print(f"plot_generic for {var_name}")
+    param_dict = {
+        "2nd_var": False,
+        "axhline_y": 0,
+        "check_exp_ocean": False,
+        "check_exp_vol": False,
+        "check_exp_year": True,
+        "default_ylim": [],
+        "do_add_line": True,
+        "do_add_trend": True,
+        "format": "%4.2f",
+        "glb_only": False,
+        "lw": 1.0,
+        "ohc": False,
+        "set_axhline": False,
+        "set_legend": True,
+        "shorten_year": False,
+        "title": var_name,
+        "use_getmoc": False,
+        "var": lambda exp: np.array(exp["annual"][var_name][0]),
+        "verbose": False,
+        "vol": False,
+        "ylabel": lambda exp: np.array(exp["annual"][var_name][1]),
+    }
+    plot(ax, xlim, exps, param_dict, rgn)
+
+
+# FIXME: C901 'plot' is too complex (19)
+def plot(ax, xlim, exps, param_dict, rgn):  # noqa: C901
     if param_dict["glb_only"] and (rgn != "glb"):
         return
     ax.set_xlim(xlim)
     extreme_values = []
     for exp in exps:
+        # Relevant to "Plot 5: plot_change_ohc"
         if param_dict["check_exp_ocean"] and (exp["ocean"] is None):
             continue
+        # Relevant to "Plot 7: plot_change_sea_level"
+        # This must be checked before plot 6,
+        # otherwise, `param_dict["var"]` will be run,
+        # but `exp["annual"]["volume"]` won't exist.
         if param_dict["check_exp_vol"] and (exp["vol"] is None):
             continue
+        # Relevant to "Plot 6: plot_max_moc"
         if param_dict["use_getmoc"]:
             if exp["moc"]:
                 [year, var] = getmoc(exp["moc"])
@@ -399,25 +468,27 @@ def plot(ax, xlim, exps, param_dict, rgn):
         extreme_values.append(np.amin(var))
         if param_dict["shorten_year"]:
             year = year[: len(var)]
-        ax.plot(
-            year,
-            var,
-            lw=param_dict["lw"],
-            marker=None,
-            c=exp["color"],
-            label=exp["name"],
-        )
+        try:
+            ax.plot(
+                year,
+                var,
+                lw=param_dict["lw"],
+                marker=None,
+                c=exp["color"],
+                label=exp["name"],
+            )
+        except Exception:
+            raise RuntimeError(f"{param_dict['title']} could not be plotted.")
         if param_dict["2nd_var"]:
             # Specifically for plot_toa_radiation
             # TODO: if more plots require a 2nd variable, we can change `var` to be a list,
             # but that will be a more significant refactoring.
-            var = np.array(exp["annual"]["FLUT"])
+            var = np.array(exp["annual"]["FLUT"][0])
             ax.plot(year, var, lw=1.0, marker=None, ls=":", c=exp["color"])
             continue
         if param_dict["check_exp_year"] and exp["yr"] is None:
             continue
         elif param_dict["do_add_line"] or param_dict["do_add_trend"]:
-            print(exp["name"])
             for yrs in exp["yr"]:
                 if param_dict["do_add_line"]:
                     add_line(
@@ -444,13 +515,17 @@ def plot(ax, xlim, exps, param_dict, rgn):
                         verbose=param_dict["verbose"],
                         vol=param_dict["vol"],
                     )
-
-    ax.set_ylim(get_ylim(param_dict["default_ylim"], extreme_values))
+    ylim = get_ylim(param_dict["default_ylim"], extreme_values)
+    ax.set_ylim(ylim)
     if param_dict["set_axhline"]:
         ax.axhline(y=param_dict["axhline_y"], lw=1, c="0.5")
     ax.set_title(param_dict["title"])
     ax.set_xlabel("Year")
-    ax.set_ylabel(param_dict["ylabel"])
+    units = param_dict["ylabel"]
+    c = callable(units)
+    if c:
+        units = units(exps[0])
+    ax.set_ylabel(units)
     if param_dict["set_legend"]:
         ax.legend(loc="best")
 
@@ -467,6 +542,123 @@ def plot(ax, xlim, exps, param_dict, rgn):
 }
 
 
+def param_get_list(param_value):
+    if param_value == "None":
+        return []
+    else:
+        return param_value.split(",")
+
+
+def set_var(exp, exp_key, var_list, valid_vars, invalid_vars, rgn):
+    if exp[exp_key] is not None:
+        ts = TS(exp[exp_key])
+        for var in var_list:
+            try:
+                v, units = ts.globalAnnual(var)
+                valid_vars.append(str(var))
+            except Exception as e:
+                print(e)
+                print(f"globalAnnual failed. Invalid var = {var}")
+                invalid_vars.append(str(var))
+                continue
+            if len(v.shape) > 1:
+                # number of years x 3 regions = v.shape
+                # 3 regions = global, northern hemisphere, southern hemisphere
+                # We get here if we used the updated `ts` task
+                # (using `rgn_avg` rather than `glb_avg`).
+                if rgn == "glb":
+                    n = 0
+                elif rgn == "n":
+                    n = 1
+                elif rgn == "s":
+                    n = 2
+                else:
+                    raise RuntimeError(f"Invalid rgn={rgn}")
+                v = v[:, n]  # Just use nth column
+            elif rgn != "glb":
+                # v only has one dimension -- glb.
+                # Therefore it is not possible to get n or s plots.
+                raise RuntimeError(
+                    f"var={var} only has global data. Cannot process rgn={rgn}"
+                )
+            exp["annual"][var] = v
+            exp["annual"][var] = (v, units)
+            if "year" not in exp["annual"]:
+                time = v.getTime()
+                exp["annual"]["year"] = [x.year for x in time.asComponentTime()]
+        del ts
+
+
+def make_plot_pdfs(
+    figstr, rgn, component, xlim, exps, plot_list, valid_plots, invalid_plots
+):
+    num_plots = len(plot_list)
+    if num_plots == 0:
+        return
+    nrows = 4
+    ncols = 2
+    plots_per_page = nrows * ncols
+    num_pages = math.ceil(num_plots / plots_per_page)
+
+    counter = 0
+    # https://stackoverflow.com/questions/58738992/save-multiple-figures-with-subplots-into-a-pdf-with-multiple-pages
+    pdf = matplotlib.backends.backend_pdf.PdfPages(f"{figstr}_{rgn}_{component}.pdf")
+    for page in range(num_pages):
+        fig = plt.figure(1, figsize=[13.5, 16.5])
+        fig.suptitle(f"{figstr}_{rgn}_{component}")
+        for j in range(plots_per_page):
+            # The final page doesn't need to be filled out with plots.
+            if counter >= num_plots:
+                break
+            ax = plt.subplot(nrows, ncols, j + 1)
+            if component == "original":
+                try:
+                    plot_function = PLOT_DICT[plot_list[counter]]
+                except KeyError:
+                    raise KeyError(f"Invalid plot name: {plot_list[counter]}")
+                try:
+                    plot_function(ax, xlim, exps, rgn)
+                    valid_plots.append(plot_list[counter])
+                except Exception:
+                    traceback.print_exc()
+                    plot_name = plot_list[counter]
+                    required_vars = []
+                    if plot_name == "net_toa_flux_restom":
+                        required_vars = ["RESTOM"]
+                    elif plot_name == "net_atm_energy_imbalance":
+                        required_vars = ["RESTOM", "RESSURF"]
+                    elif plot_name == "global_surface_air_temperature":
+                        required_vars = ["TREFHT"]
+                    elif plot_name == "toa_radiation":
+                        required_vars = ["FSNTOA", "FLUT"]
+                    elif plot_name == "net_atm_water_imbalance":
+                        required_vars = ["PRECC", "PRECL", "QFLX"]
+                    print(
+                        f"Failed plot_function for {plot_name}. Check that {required_vars} are available."
+                    )
+                    invalid_plots.append(plot_name)
+                counter += 1
+            else:
+                try:
+                    plot_name = plot_list[counter]
+                    plot_generic(ax, xlim, exps, plot_name, rgn)
+                    valid_plots.append(plot_name)
+                except Exception:
+                    traceback.print_exc()
+                    print(f"plot_generic failed. Invalid plot={plot_name}")
+                    invalid_plots.append(plot_name)
+                counter += 1
+
+        fig.tight_layout()
+        pdf.savefig(1)
+        if num_pages > 1:
+            fig.savefig(f"{figstr}_{rgn}_{component}_{page}.png", dpi=150)
+        else:
+            fig.savefig(f"{figstr}_{rgn}_{component}.png", dpi=150)
+        plt.clf()
+    pdf.close()
+
+
 # -----------------------------------------------------------------------------
 # FIXME: C901 'run' is too complex (19)
 def run(parameters, rgn):  # noqa: C901
@@ -514,26 +706,63 @@ def run(parameters, rgn):  # noqa: C901
     year2 = int(parameters[5])
     color = parameters[6]
     ts_num_years = parameters[7]
-    if parameters[8].lower() == "false":
+    plots_original = param_get_list(parameters[8])
+    if parameters[9].lower() == "false":
         atmosphere_only = False
     else:
         atmosphere_only = True
-    plot_list = parameters[9].split(",")
+    plots_atm = param_get_list(parameters[10])
+    plots_ice = param_get_list(parameters[11])
+    plots_lnd = param_get_list(parameters[12])
+    plots_ocn = param_get_list(parameters[13])
+    vars_original = []
+    if "net_toa_flux_restom" or "net_atm_energy_imbalance" in plots_original:
+        vars_original.append("RESTOM")
+    if "net_atm_energy_imbalance" in plots_original:
+        vars_original.append("RESSURF")
+    if "global_surface_air_temperature" in plots_original:
+        vars_original.append("TREFHT")
+    if "toa_radiation" in plots_original:
+        vars_original.append("FSNTOA")
+        vars_original.append("FLUT")
+    if "net_atm_water_imbalance" in plots_original:
+        vars_original.append("PRECC")
+        vars_original.append("PRECL")
+        vars_original.append("QFLX")
+    use_atmos = plots_atm or plots_original
+    has_original_ocn_plots = (
+        ("change_ohc" in plots_original)
+        or ("max_moc" in plots_original)
+        or ("change_sea_level" in plots_original)
+    )
+    use_ocn = plots_ocn or (not atmosphere_only and has_original_ocn_plots)
     exps = [
         {
-            "atmos": "{}/post/atm/glb/ts/monthly/{}yr/glb.xml".format(
+            "atmos": None
+            if not use_atmos
+            else "{}/post/atm/glb/ts/monthly/{}yr/glb.xml".format(
+                case_dir, ts_num_years
+            ),
+            "ice": None
+            if not plots_ice
+            else "{}/post/ice/glb/ts/monthly/{}yr/glb.xml".format(
+                case_dir, ts_num_years
+            ),
+            "land": None
+            if not plots_lnd
+            else "{}/post/lnd/glb/ts/monthly/{}yr/glb.xml".format(
                 case_dir, ts_num_years
             ),
             "ocean": None
-            if atmosphere_only
+            if not use_ocn
             else "{}/post/ocn/glb/ts/monthly/{}yr/glb.xml".format(
                 case_dir, ts_num_years
             ),
             "moc": None
-            if atmosphere_only
+            if not use_ocn
             else "{}/post/ocn/glb/ts/monthly/{}yr/".format(case_dir, ts_num_years),
             "vol": None
-            if atmosphere_only
+            if not use_ocn
             else "{}/post/ocn/glb/ts/monthly/{}yr/glb.xml".format(
                 case_dir, ts_num_years
             ),
@@ -544,97 +773,76 @@ def run(parameters, rgn):  # noqa: C901
         }
     ]
 
-    # Variables to extract
-    vars = ["RESTOM", "RESSURF", "TREFHT", "FSNTOA", "FLUT", "PRECC", "PRECL", "QFLX"]
+    valid_vars: List[str] = []
+    invalid_vars: List[str] = []
 
     # Read data
     exp: Any
     for exp in exps:
-        print(exp["atmos"])
-        ts = TS(exp["atmos"])
         exp["annual"] = {}
-        for var in vars:
-            print(var)
-            v = ts.globalAnnual(var)
-            if len(v.shape) > 1:
-                # number of years x 3 regions = v.shape
-                # 3 regions = global, northern hemisphere, southern hemisphere
-                # We get here if we used the updated `ts` task
-                # (using `rgn_avg` rather than `glb_avg`).
-                if rgn == "glb":
-                    n = 0
-                elif rgn == "n":
-                    n = 1
-                elif rgn == "s":
-                    n = 2
-                else:
-                    raise RuntimeError(f"Invalid rgn={rgn}")
-                v = v[:, n]  # Just use nth column
-            elif rgn != "glb":
-                # v only has one dimension -- glb.
-                # Therefore it is not possible to get n or s plots.
-                raise RuntimeError(
-                    f"var={var} only has global data. Cannot process rgn={rgn}"
-                )
-            exp["annual"][var] = v
-            if "year" not in exp["annual"]:
-                time = v.getTime()
-                exp["annual"]["year"] = [x.year for x in time.asComponentTime()]
-        del ts
+
+        # Use vars_original rather than plots_original,
+        # since the plots have different names than the variables
+        set_var(exp, "atmos", vars_original, valid_vars, invalid_vars, rgn)
+        set_var(exp, "atmos", plots_atm, valid_vars, invalid_vars, rgn)
+        set_var(exp, "ice", plots_ice, valid_vars, invalid_vars, rgn)
+        set_var(exp, "land", plots_lnd, valid_vars, invalid_vars, rgn)
+        set_var(exp, "ocean", plots_ocn, valid_vars, invalid_vars, rgn)
 
         # Optionally read ohc
         if exp["ocean"] is not None:
             ts = TS(exp["ocean"])
-            exp["annual"]["ohc"] = ts.globalAnnual("ohc")
+            exp["annual"]["ohc"], _ = ts.globalAnnual("ohc")
             # annomalies with respect to first year
             exp["annual"]["ohc"][:] = exp["annual"]["ohc"][:] - exp["annual"]["ohc"][0]
 
         if exp["vol"] is not None:
             ts = TS(exp["vol"])
-            exp["annual"]["volume"] = ts.globalAnnual("volume")
+            exp["annual"]["volume"], _ = ts.globalAnnual("volume")
             # annomalies with respect to first year
             exp["annual"]["volume"][:] = (
                 exp["annual"]["volume"][:] - exp["annual"]["volume"][0]
             )
 
+    print(
+        f"{rgn} region globalAnnual was computed successfully for these variables: {valid_vars}"
+    )
+    print(
+        f"{rgn} region globalAnnual could not be computed for these variables: {invalid_vars}"
+    )
+
     # -----------------------------------------------------------------------------
     # --- Generate plots ---
 
     xlim = [float(year1), float(year2)]
 
-    num_plots = len(plot_list)
-    nrows = 4
-    ncols = 2
-    plots_per_page = nrows * ncols
-    num_pages = math.ceil(num_plots / plots_per_page)
-
-    i = 0
-    # https://stackoverflow.com/questions/58738992/save-multiple-figures-with-subplots-into-a-pdf-with-multiple-pages
-    pdf = matplotlib.backends.backend_pdf.PdfPages(f"{figstr}_{rgn}.pdf")
-    for page in range(num_pages):
-        fig = plt.figure(1, figsize=[13.5, 16.5])
-        fig.suptitle(f"{figstr}_{rgn}")
-        for j in range(plots_per_page):
-            if i < num_plots:
-                ax = plt.subplot(nrows, ncols, j + 1)
-                try:
-                    PLOT_DICT[plot_list[i]](ax, xlim, exps, rgn)
-                except KeyError:
-                    raise KeyError(f"Invalid plot name: {plot_list[i]}")
-                i += 1
-
-        fig.tight_layout()
-        pdf.savefig(1)
-        if num_pages > 1:
-            fig.savefig(f"{figstr}_{rgn}_{page}.png", dpi=150)
-        else:
-            fig.savefig(f"{figstr}_{rgn}.png", dpi=150)
-        plt.clf()
-    pdf.close()
+    valid_plots: List[str] = []
+    invalid_plots: List[str] = []
+
+    make_plot_pdfs(
+        figstr, rgn, "original", xlim, exps, plots_original, valid_plots, invalid_plots
+    )
+    make_plot_pdfs(
+        figstr, rgn, "atm", xlim, exps, plots_atm, valid_plots, invalid_plots
+    )
+    make_plot_pdfs(
+        figstr, rgn, "ice", xlim, exps, plots_ice, valid_plots, invalid_plots
+    )
+    make_plot_pdfs(
+        figstr, rgn, "lnd", xlim, exps, plots_lnd, valid_plots, invalid_plots
+    )
+    make_plot_pdfs(
+        figstr, rgn, "ocn", xlim, exps, plots_ocn, valid_plots, invalid_plots
+    )
+
+    print(f"These {rgn} region plots generated successfully: {valid_plots}")
+    print(
+        f"These {rgn} region plots could not be generated successfully: {invalid_plots}"
+    )
 
 
 def run_by_region(parameters):
-    regions = parameters[10].split(",")
+    regions = parameters[14].split(",")
     for rgn in regions:
         if rgn.lower() in ["glb", "global"]:
             rgn = "glb"
diff --git a/zppy/templates/default.ini b/zppy/templates/default.ini
index ad56c93b..0f85e7ff 100644
--- a/zppy/templates/default.ini
+++ b/zppy/templates/default.ini
@@ -123,10 +123,10 @@ scratch = string(default="")
 # See https://e3sm-project.github.io/e3sm_diags/_build/html/master/available-parameters.html
 backend = string(default="mpl")
 cfg = string(default="")
-# Name of the subsection of `[climo]` to use for "diurnal_cycle" runs
+# Name of the frequency from `[climo]` to use for "diurnal_cycle" runs
 climo_diurnal_frequency = string(default="")
 # climo_diurnal_input_files -- NOTE: this parameter is created internally
-# Name of the frequency from `[climo]` to use for "diurnal_cycle" runs
+# Name of the subsection of `[climo]` to use for "diurnal_cycle" runs
 climo_diurnal_subsection = string(default="")
 # Name of the `[climo]` subtask for "diurnal_cycle" runs
 climo_subsection = string(default="")
@@ -284,7 +284,8 @@ stream_ocn = string(default="streams.ocean")
 walltime = string(default="06:00:00")
 
 [global_time_series]
-# Set to True to skip figures requiring ocean data
+# Deprecated; kept for backwards compatibility
+# Set to True to skip figures requiring ocean data; only affects original plots
 atmosphere_only = boolean(default=False)
 climo_years = string_list(default=list(""))
 # The color to be used for the graphs.
@@ -296,8 +297,26 @@ figstr = string(default="")
 # NOTE: always overrides value in [default]
 input_subdir = string(default="archive/ocn/hist")
 moc_file = string(default="")
+# Deprecated; legacy name for plots_original; kept for backwards compatibility
+# plots_original replaces it with the same default.
+# So, if a cfg used the default value before, the behavior will remain the same.
+# And if a cfg set plot_names explicitally, then global_time_series.py will override plots_original.
+plot_names = string(default="")
 # The names of the plots you want displayed
-plot_names = string(default="net_toa_flux_restom,global_surface_air_temperature,toa_radiation,net_atm_energy_imbalance,change_ohc,max_moc,change_sea_level,net_atm_water_imbalance")
+# Variable requirements:
+# net_toa_flux_restom requires RESTOM
+# net_atm_energy_imbalance requires RESTOM, RESSURF
+# global_surface_air_temperature requires TREFHT
+# toa_radiation requires FSNTOA, FLUT
+# net_atm_water_imbalance requires PRECC, PRECL, QFLX
+# Remove the 3 ocean plots (change_ohc,max_moc,change_sea_level) if you don't have ocean data.
+plots_original = string(default="net_toa_flux_restom,global_surface_air_temperature,toa_radiation,net_atm_energy_imbalance,change_ohc,max_moc,change_sea_level,net_atm_water_imbalance")
+# The names of the extra plots you want displayed (i.e., define specific variables)
+# These should be a subset of the `vars` generated by the `ts` `glb` subtasks.
+plots_atm = string(default="")
+plots_ice = string(default="")
+plots_lnd = string(default="")
+plots_ocn = string(default="")
 # regions to plot: glb, n, s (global, northern hemisphere, southern hemisphere)
 regions = string(default="glb,n,s")
 # The number of years in a time-series file
diff --git a/zppy/templates/global_time_series.bash b/zppy/templates/global_time_series.bash
index 8b69c7b0..2154e66b 100644
--- a/zppy/templates/global_time_series.bash
+++ b/zppy/templates/global_time_series.bash
@@ -34,24 +34,36 @@ www={{ www }}
 case_dir={{ output }}
 global_ts_dir={{ global_time_series_dir }}
 
-# Options
-atmosphere_only={{ atmosphere_only }}
-atmosphere_only=${atmosphere_only,,}
-
 ################################################################################
 
-echo 'Create xml files for atm'
 export CDMS_NO_MPI=true
-cd ${case_dir}/post/atm/glb/ts/monthly/${ts_num_years}yr
-cdscan -x glb.xml *.nc
-if [ $? != 0 ]; then
-  cd {{ scriptDir }}
-  echo 'ERROR (1)' > {{ prefix }}.status
-  exit 1
+
+use_atm={{ use_atm }}
+if [[ ${use_atm,,} == "true" ]]; then
+    echo 'Create xml files for atm'
+    cd ${case_dir}/post/atm/glb/ts/monthly/${ts_num_years}yr
+    cdscan -x glb.xml *.nc
+    if [ $? != 0 ]; then
+      cd {{ scriptDir }}
+      echo 'ERROR (1)' > {{ prefix }}.status
+      exit 1
+    fi
 fi
 
-if [[ ${atmosphere_only} == "false" ]]; then
+use_lnd={{ use_lnd }}
+if [[ ${use_lnd,,} == "true" ]]; then
+    echo 'Create xml files for lnd'
+    cd ${case_dir}/post/lnd/glb/ts/monthly/${ts_num_years}yr
+    cdscan -x glb.xml *.nc
+    if [ $? != 0 ]; then
+      cd {{ scriptDir }}
+      echo 'ERROR (2)' > {{ prefix }}.status
+      exit 2
+    fi
+fi
 
+use_ocn={{ use_ocn }}
+if [[ ${use_ocn,,} == "true" ]]; then
     echo 'Create ocean time series'
     cd ${global_ts_dir}
     mkdir -p ${case_dir}/post/ocn/glb/ts/monthly/${ts_num_years}yr
@@ -59,19 +71,18 @@ if [[ ${atmosphere_only} == "false" ]]; then
     python ocean_month.py ${input} ${case_dir} ${start_yr} ${end_yr} ${ts_num_years}
     if [ $? != 0 ]; then
       cd {{ scriptDir }}
-      echo 'ERROR (2)' > {{ prefix }}.status
-      exit 2
+      echo 'ERROR (3)' > {{ prefix }}.status
+      exit 3
     fi
 
-
     echo 'Create xml for for ocn'
     export CDMS_NO_MPI=true
     cd ${case_dir}/post/ocn/glb/ts/monthly/${ts_num_years}yr
     cdscan -x glb.xml mpaso.glb*.nc
     if [ $? != 0 ]; then
       cd {{ scriptDir }}
-      echo 'ERROR (3)' > {{ prefix }}.status
-      exit 3
+      echo 'ERROR (4)' > {{ prefix }}.status
+      exit 4
     fi
 
     echo 'Copy moc file'
@@ -79,19 +90,20 @@ if [[ ${atmosphere_only} == "false" ]]; then
     cp ${moc_file} ../../../../../ocn/glb/ts/monthly/${ts_num_years}yr/
     if [ $? != 0 ]; then
       cd {{ scriptDir }}
-      echo 'ERROR (4)' > {{ prefix }}.status
-      exit 4
+      echo 'ERROR (5)' > {{ prefix }}.status
+      exit 5
     fi
 
 fi
 
 echo 'Update time series figures'
 cd ${global_ts_dir}
-python coupled_global.py ${case_dir} ${experiment_name} ${figstr} ${start_yr} ${end_yr} {{ color }} ${ts_num_years} ${atmosphere_only} "{{ plot_names }}" {{ regions }}
+atmosphere_only={{ atmosphere_only }}
+python coupled_global.py ${case_dir} ${experiment_name} ${figstr} ${start_yr} ${end_yr} {{ color }} ${ts_num_years} {{ plots_original }} ${atmosphere_only,,} {{ plots_atm }} {{ plots_ice }} {{ plots_lnd }} {{ plots_ocn }} {{ regions }}
 if [ $? != 0 ]; then
   cd {{ scriptDir }}
-  echo 'ERROR (5)' > {{ prefix }}.status
-  exit 5
+  echo 'ERROR (6)' > {{ prefix }}.status
+  exit 6
 fi
 
 
@@ -113,8 +125,8 @@ f=${www}/${case}/global_time_series
 mkdir -p ${f}
 if [ $? != 0 ]; then
   cd {{ scriptDir }}
-  echo 'ERROR (6)' > {{ prefix }}.status
-  exit 6
+  echo 'ERROR (7)' > {{ prefix }}.status
+  exit 7
 fi
 
 {% if machine in ['pm-cpu', 'pm-gpu'] %}
@@ -135,8 +147,8 @@ done
 rsync -a --delete ${results_dir_absolute_path} ${www}/${case}/global_time_series
 if [ $? != 0 ]; then
   cd {{ scriptDir }}
-  echo 'ERROR (7)' > {{ prefix }}.status
-  exit 7
+  echo 'ERROR (8)' > {{ prefix }}.status
+  exit 8
 fi
 
 {% if machine in ['pm-cpu', 'pm-gpu'] %}
diff --git a/zppy/templates/readTS.py b/zppy/templates/readTS.py
index 5527e3e3..ccd77ad0 100644
--- a/zppy/templates/readTS.py
+++ b/zppy/templates/readTS.py
@@ -15,6 +15,8 @@ def __del__(self):
 
     def globalAnnual(self, var):
 
+        units = None
+
         # Constants, from AMWG diagnostics
         Lv = 2.501e6
         Lf = 3.337e5
@@ -22,46 +24,49 @@ def globalAnnual(self, var):
         # Is this a derived variable?
         if var == "RESTOM":
 
-            FSNT = self.globalAnnual("FSNT")
-            FLNT = self.globalAnnual("FLNT")
+            FSNT, _ = self.globalAnnual("FSNT")
+            FLNT, _ = self.globalAnnual("FLNT")
             v = FSNT - FLNT
 
         elif var == "RESTOA":
 
             print("NOT READY")
-            FSNTOA = self.globalAnnual("FSNTOA")
-            FLUT = self.globalAnnual("FLUT")
+            FSNTOA, _ = self.globalAnnual("FSNTOA")
+            FLUT, _ = self.globalAnnual("FLUT")
             v = FSNTOA - FLUT
 
         elif var == "LHFLX":
 
-            QFLX = self.globalAnnual("QFLX")
-            PRECC = self.globalAnnual("PRECC")
-            PRECL = self.globalAnnual("PRECL")
-            PRECSC = self.globalAnnual("PRECSC")
-            PRECSL = self.globalAnnual("PRECSL")
+            QFLX, _ = self.globalAnnual("QFLX")
+            PRECC, _ = self.globalAnnual("PRECC")
+            PRECL, _ = self.globalAnnual("PRECL")
+            PRECSC, _ = self.globalAnnual("PRECSC")
+            PRECSL, _ = self.globalAnnual("PRECSL")
             v = (Lv + Lf) * QFLX - Lf * 1.0e3 * (PRECC + PRECL - PRECSC - PRECSL)
 
         elif var == "RESSURF":
 
-            FSNS = self.globalAnnual("FSNS")
-            FLNS = self.globalAnnual("FLNS")
-            SHFLX = self.globalAnnual("SHFLX")
-            LHFLX = self.globalAnnual("LHFLX")
+            FSNS, _ = self.globalAnnual("FSNS")
+            FLNS, _ = self.globalAnnual("FLNS")
+            SHFLX, _ = self.globalAnnual("SHFLX")
+            LHFLX, _ = self.globalAnnual("LHFLX")
             v = FSNS - FLNS - SHFLX - LHFLX
 
         elif var == "PREC":
 
-            PRECC = self.globalAnnual("PRECC")
-            PRECL = self.globalAnnual("PRECL")
+            PRECC, _ = self.globalAnnual("PRECC")
+            PRECL, _ = self.globalAnnual("PRECL")
             v = 1.0e3 * (PRECC + PRECL)
 
         else:
+            # Non-derived variables
 
             # Read variable
             v = self.f(var)
+            units = v.units
 
             # Annual average
+            cdutil.setTimeBoundsMonthly(v)
             v = cdutil.YEAR(v)
 
-        return v
+        return v, units
diff --git a/zppy/templates/ts.bash b/zppy/templates/ts.bash
index 0b03cedd..1e811976 100644
--- a/zppy/templates/ts.bash
+++ b/zppy/templates/ts.bash
@@ -83,7 +83,6 @@ cat input.txt | ncclimo \
 {%- if vars != '' %}
 -v ${vars} \
 {%- endif %}
---mem_mb=0 \
 --split \
 {%- if extra_vars != '' %}
 --var_xtr={{extra_vars}} \