diff --git a/tests/global_time_series/test_coupled_global.py b/tests/global_time_series/test_coupled_global.py
new file mode 100644
index 00000000..ff8891d6
--- /dev/null
+++ b/tests/global_time_series/test_coupled_global.py
@@ -0,0 +1,355 @@
+import unittest
+from typing import Any, Dict, List
+
+from zppy.templates.coupled_global import (
+    Metric,
+    Parameters,
+    Variable,
+    VariableGroup,
+    construct_generic_variables,
+    get_data_dir,
+    get_exps,
+    get_region,
+    get_variable_groups,
+    get_vars_original,
+    get_ylim,
+    land_csv_row_to_var,
+    param_get_list,
+)
+
+
+# Helper function
+def get_var_names(vars: List[Variable]):
+    return list(map(lambda v: v.variable_name, vars))
+
+
+# Run this test suite in the environment the global_time_series task runs in.
+# I.e., whatever `environment_commands` is set to for `[global_time_series]`
+# NOT the zppy dev environment.
+# Run: python -u -m unittest tests/global_time_series/test_coupled_global.py
+class TestCoupledGlobal(unittest.TestCase):
+
+    # Useful classes and their helper functions ###############################
+    def test_param_get_list(self):
+        self.assertEqual(param_get_list("None"), [])
+
+        self.assertEqual(param_get_list("a"), ["a"])
+        self.assertEqual(param_get_list("a,b,c"), ["a", "b", "c"])
+
+        self.assertEqual(param_get_list(""), [""])
+        self.assertEqual(param_get_list("a,"), ["a", ""])
+        self.assertEqual(param_get_list("a,b,c,"), ["a", "b", "c", ""])
+
+    def test_get_region(self):
+        self.assertEqual(get_region("glb"), "glb")
+        self.assertEqual(get_region("global"), "glb")
+        self.assertEqual(get_region("GLB"), "glb")
+        self.assertEqual(get_region("Global"), "glb")
+
+        self.assertEqual(get_region("n"), "n")
+        self.assertEqual(get_region("north"), "n")
+        self.assertEqual(get_region("northern"), "n")
+        self.assertEqual(get_region("N"), "n")
+        self.assertEqual(get_region("North"), "n")
+        self.assertEqual(get_region("Northern"), "n")
+
+        self.assertEqual(get_region("s"), "s")
+        self.assertEqual(get_region("south"), "s")
+        self.assertEqual(get_region("southern"), "s")
+        self.assertEqual(get_region("S"), "s")
+        self.assertEqual(get_region("South"), "s")
+        self.assertEqual(get_region("Southern"), "s")
+
+        self.assertRaises(ValueError, get_region, "not-a-region")
+
+    def test_Parameters_and_related_functions(self):
+        # Consider the following parameters given by a user.
+        args = [
+            "coupled_global.py",
+            "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051",
+            "v3.LR.historical_0051",
+            "v3.LR.historical_0051",
+            "1985",
+            "1989",
+            "Blue",
+            "5",
+            "None",
+            "false",
+            "TREFHT",
+            "None",
+            "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",
+            "None",
+            "1",
+            "1",
+            "glb,n,s",
+        ]
+        # Then:
+        parameters: Parameters = Parameters(args)
+        self.assertEqual(
+            parameters.case_dir,
+            "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051",
+        )
+        self.assertEqual(parameters.experiment_name, "v3.LR.historical_0051")
+        self.assertEqual(parameters.figstr, "v3.LR.historical_0051")
+        self.assertEqual(parameters.year1, 1985)
+        self.assertEqual(parameters.year2, 1989)
+        self.assertEqual(parameters.color, "Blue")
+        self.assertEqual(parameters.ts_num_years_str, "5")
+        self.assertEqual(parameters.plots_original, [])
+        self.assertEqual(parameters.atmosphere_only, False)
+        self.assertEqual(parameters.plots_atm, ["TREFHT"])
+        self.assertEqual(parameters.plots_ice, [])
+        self.assertEqual(
+            parameters.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",
+            ],
+        )
+        self.assertEqual(parameters.plots_ocn, [])
+        self.assertEqual(parameters.nrows, 1)
+        self.assertEqual(parameters.ncols, 1)
+        self.assertEqual(parameters.regions, ["glb", "n", "s"])
+
+        # test_get_data_dir
+        self.assertEqual(
+            get_data_dir(parameters, "atm", True),
+            "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051/post/atm/glb/ts/monthly/5yr/",
+        )
+        self.assertEqual(get_data_dir(parameters, "atm", False), "")
+        self.assertEqual(
+            get_data_dir(parameters, "ice", True),
+            "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051/post/ice/glb/ts/monthly/5yr/",
+        )
+        self.assertEqual(get_data_dir(parameters, "ice", False), "")
+        self.assertEqual(
+            get_data_dir(parameters, "lnd", True),
+            "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051/post/lnd/glb/ts/monthly/5yr/",
+        )
+        self.assertEqual(get_data_dir(parameters, "lnd", False), "")
+        self.assertEqual(
+            get_data_dir(parameters, "ocn", True),
+            "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051/post/ocn/glb/ts/monthly/5yr/",
+        )
+        self.assertEqual(get_data_dir(parameters, "ocn", False), "")
+
+        # test_get_exps
+        self.maxDiff = None
+        exps: List[Dict[str, Any]] = get_exps(parameters)
+        self.assertEqual(len(exps), 1)
+        expected = {
+            "atmos": "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051/post/atm/glb/ts/monthly/5yr/",
+            "ice": "",
+            "land": "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051/post/lnd/glb/ts/monthly/5yr/",
+            "ocean": "",
+            "moc": "",
+            "vol": "",
+            "name": "v3.LR.historical_0051",
+            "yoffset": 0.0,
+            "yr": ([1985, 1989],),
+            "color": "Blue",
+        }
+        self.assertEqual(exps[0], expected)
+        # Change up parameters
+        parameters.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".split(
+            ","
+        )
+        parameters.plots_atm = []
+        parameters.plots_lnd = []
+        exps = get_exps(parameters)
+        self.assertEqual(len(exps), 1)
+        expected = {
+            "atmos": "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051/post/atm/glb/ts/monthly/5yr/",
+            "ice": "",
+            "land": "",
+            "ocean": "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051/post/ocn/glb/ts/monthly/5yr/",
+            "moc": "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051/post/ocn/glb/ts/monthly/5yr/",
+            "vol": "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051/post/ocn/glb/ts/monthly/5yr/",
+            "name": "v3.LR.historical_0051",
+            "yoffset": 0.0,
+            "yr": ([1985, 1989],),
+            "color": "Blue",
+        }
+        self.assertEqual(exps[0], expected)
+        # Change up parameters
+        parameters.atmosphere_only = True
+        exps = get_exps(parameters)
+        self.assertEqual(len(exps), 1)
+        expected = {
+            "atmos": "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/test-616-20240930/v3.LR.historical_0051/post/atm/glb/ts/monthly/5yr/",
+            "ice": "",
+            "land": "",
+            "ocean": "",
+            "moc": "",
+            "vol": "",
+            "name": "v3.LR.historical_0051",
+            "yoffset": 0.0,
+            "yr": ([1985, 1989],),
+            "color": "Blue",
+        }
+        self.assertEqual(exps[0], expected)
+
+    # Metric
+
+    def test_Variable(self):
+        v = Variable(
+            "var_name",
+            original_units="units1",
+            final_units="units2",
+            group="group_name",
+            long_name="long name",
+        )
+        self.assertEqual(v.variable_name, "var_name")
+        self.assertEqual(v.metric, Metric.AVERAGE)
+        self.assertEqual(v.scale_factor, 1.0)
+        self.assertEqual(v.original_units, "units1")
+        self.assertEqual(v.final_units, "units2")
+        self.assertEqual(v.group, "group_name")
+        self.assertEqual(v.long_name, "long name")
+
+    def test_get_vars_original(self):
+        self.assertEqual(
+            get_var_names(get_vars_original(["net_toa_flux_restom"])), ["RESTOM"]
+        )
+        self.assertEqual(
+            get_var_names(get_vars_original(["net_atm_energy_imbalance"])),
+            ["RESTOM", "RESSURF"],
+        )
+        self.assertEqual(
+            get_var_names(get_vars_original(["global_surface_air_temperature"])),
+            ["TREFHT"],
+        )
+        self.assertEqual(
+            get_var_names(get_vars_original(["toa_radiation"])), ["FSNTOA", "FLUT"]
+        )
+        self.assertEqual(
+            get_var_names(get_vars_original(["net_atm_water_imbalance"])),
+            ["PRECC", "PRECL", "QFLX"],
+        )
+        self.assertEqual(
+            get_var_names(
+                get_vars_original(
+                    [
+                        "net_toa_flux_restom",
+                        "net_atm_energy_imbalance",
+                        "global_surface_air_temperature",
+                        "toa_radiation",
+                        "net_atm_water_imbalance",
+                    ]
+                )
+            ),
+            ["RESTOM", "RESSURF", "TREFHT", "FSNTOA", "FLUT", "PRECC", "PRECL", "QFLX"],
+        )
+        self.assertEqual(get_var_names(get_vars_original(["invalid_plot"])), [])
+
+    def test_land_csv_row_to_var(self):
+        # Test with first row of land csv, whitespace stripped
+        csv_row = "BCDEP,A,1.00000E+00,kg/m^2/s,kg/m^2/s,Aerosol Flux,total black carbon deposition (dry+wet) from atmosphere".split(
+            ","
+        )
+        v: Variable = land_csv_row_to_var(csv_row)
+        self.assertEqual(v.variable_name, "BCDEP")
+        self.assertEqual(v.metric, Metric.AVERAGE)
+        self.assertEqual(v.scale_factor, 1.0)
+        self.assertEqual(v.original_units, "kg/m^2/s")
+        self.assertEqual(v.final_units, "kg/m^2/s")
+        self.assertEqual(v.group, "Aerosol Flux")
+        self.assertEqual(
+            v.long_name, "total black carbon deposition (dry+wet) from atmosphere"
+        )
+
+    # construct_land_variables -- requires IO
+
+    def test_construct_generic_variables(self):
+        vars: List[str] = ["a", "b", "c"]
+        self.assertEqual(get_var_names(construct_generic_variables(vars)), vars)
+
+    # RequestedVariables -- requires IO
+
+    def test_VariableGroup(self):
+        var_str_list: List[str] = ["a", "b", "c"]
+        vars: List[Variable] = construct_generic_variables(var_str_list)
+        g: VariableGroup = VariableGroup("MyGroup", vars)
+        self.assertEqual(g.group_name, "MyGroup")
+        self.assertEqual(get_var_names(g.variables), var_str_list)
+
+    # TS -- requires IO
+    # OutputViewer -- requires IO
+
+    # Setup ###################################################################
+
+    # set_var -- requires IO
+    # process_data -- requires IO
+
+    # Plotting ####################################################################
+
+    def test_get_variable_groups(self):
+        a: Variable = Variable(variable_name="a", group="GroupA")
+        b: Variable = Variable(variable_name="b", group="GroupA")
+        x: Variable = Variable(variable_name="x", group="GroupX")
+        y: Variable = Variable(variable_name="y", group="GroupX")
+
+        def get_group_names(groups: List[VariableGroup]) -> List[str]:
+            return list(map(lambda g: g.group_name, groups))
+
+        self.assertEqual(
+            get_group_names(get_variable_groups([a, b, x, y])), ["GroupA", "GroupX"]
+        )
+
+    # getmoc -- requires IO
+    # add_line -- requires IO
+    # add_trend -- requires IO
+
+    def test_get_ylim(self):
+        # Min is equal, max is equal
+        self.assertEqual(get_ylim([-1, 1], [-1, 1]), [-1, 1])
+        # Min is lower, max is equal
+        self.assertEqual(get_ylim([-1, 1], [-2, 1]), [-2, 1])
+        # Min is equal, max is higher
+        self.assertEqual(get_ylim([-1, 1], [-1, 2]), [-1, 2])
+        # Min is lower, max is higher
+        self.assertEqual(get_ylim([-1, 1], [-2, 2]), [-2, 2])
+        # Min is lower, max is higher, multiple extreme_values
+        self.assertEqual(get_ylim([-1, 1], [-2, -0.5, 0.5, 2]), [-2, 2])
+        # No standard range
+        self.assertEqual(get_ylim([], [-2, 2]), [-2, 2])
+        # No extreme range
+        self.assertEqual(get_ylim([-1, 1], []), [-1, 1])
+
+    # Plotting functions -- require IO
+    # make_plot_pdfs -- requires IO
+
+    # Run coupled_global ######################################################
+
+    # run -- requires IO
+    # get_vars -- requires IO
+    # create_viewer -- requires IO
+    # create_viewer_index -- requires IO
+    # run_by_region -- requires IO
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/integration/generated/test_debug_chrysalis.cfg b/tests/integration/generated/test_debug_chrysalis.cfg
new file mode 100644
index 00000000..e30cc9f8
--- /dev/null
+++ b/tests/integration/generated/test_debug_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_debug` 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_complete_run.py
+output = "/lcrc/group/e3sm/ac.forsyth2/zppy_test_debug_output/test-601-output-viewer-v3/v2.LR.historical_0201"
+partition = "debug"
+qos = "regular"
+www = "/lcrc/group/e3sm/public_html/diagnostic_output/ac.forsyth2/zppy_debug_www/test-601-output-viewer-v3"
+
+[ts]
+active = True
+walltime = "00:30:00"
+years = "1850:1854:2",
+
+  [[ atm_monthly_180x360_aave ]]
+  frequency = "monthly"
+  input_files = "eam.h0"
+  input_subdir = "archive/atm/hist"
+  ts_fmt = "cmip"
+
+  [[ atm_daily_180x360_aave ]]
+  frequency = "daily"
+  input_files = "eam.h1"
+  input_subdir = "archive/atm/hist"
+  vars = "PRECT"
+
+  [[ land_monthly ]]
+  e3sm_to_cmip_environment_commands = ""
+  extra_vars = "landfrac"
+  frequency = "monthly"
+  input_files = "elm.h0"
+  input_subdir = "archive/lnd/hist"
+  vars = "LAISHA,LAISUN"
+  #vars = "FSH,RH2M"
+  ts_fmt = "cmip"
+
+[ilamb]
+active = True
+grid = '180x360_aave'
+nodes = 8
+partition = "compute"
+short_name = 'v2.LR.historical_0201'
+ts_num_years = 2
+years = "1850:1854:2",
+
+  [[ land_monthly ]]
diff --git a/tests/integration/generated/test_min_case_global_time_series_custom_chrysalis.cfg b/tests/integration/generated/test_min_case_global_time_series_custom_chrysalis.cfg
index 3a648737..3204bde9 100644
--- a/tests/integration/generated/test_min_case_global_time_series_custom_chrysalis.cfg
+++ b/tests/integration/generated/test_min_case_global_time_series_custom_chrysalis.cfg
@@ -23,7 +23,7 @@ years = "1985:1995:5",
   input_files = "eam.h0"
   input_subdir = "archive/atm/hist"
   mapping_file = "glb"
-  vars = "TREFHT,AODDUST"
+  vars = "TREFHT"
 
   [[ lnd_monthly_glb ]]
   frequency = "monthly"
@@ -37,7 +37,7 @@ active = True
 experiment_name = "v3.LR.historical_0051"
 figstr = "v3.LR.historical_0051"
 plots_original=""
-plots_atm = "TREFHT,AODDUST"
+plots_atm = "TREFHT"
 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 = 5
 walltime = "00:30:00"
diff --git a/tests/integration/generated/test_min_case_global_time_series_single_plots_chrysalis.cfg b/tests/integration/generated/test_min_case_global_time_series_single_plots_chrysalis.cfg
new file mode 100644
index 00000000..0516fc09
--- /dev/null
+++ b/tests/integration/generated/test_min_case_global_time_series_single_plots_chrysalis.cfg
@@ -0,0 +1,46 @@
+[default]
+case = "v3.LR.historical_0051"
+constraint = ""
+dry_run = "False"
+environment_commands = ""
+input = /lcrc/group/e3sm2/ac.wlin/E3SMv3/v3.LR.historical_0051
+input_subdir = archive/atm/hist
+mapping_file = "map_ne30pg2_to_cmip6_180x360_aave.20200201.nc"
+output = "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_single_plots_output/unique_id/v3.LR.historical_0051"
+partition = "debug"
+qos = "regular"
+www = "/lcrc/group/e3sm/public_html/diagnostic_output/ac.forsyth2/zppy_min_case_global_time_series_single_plots_www/unique_id"
+
+[ts]
+active = True
+e3sm_to_cmip_environment_commands = ""
+walltime = "00:30:00"
+years = "1985:1995:5",
+
+  [[ 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"
+  vars = "TREFHT"
+
+  [[ lnd_monthly_glb ]]
+  frequency = "monthly"
+  input_files = "elm.h0"
+  input_subdir = "archive/lnd/hist"
+  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.historical_0051"
+figstr = "v3.LR.historical_0051"
+ncols = 1
+nrows = 1
+plots_original=""
+plots_atm = "TREFHT"
+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 = 5
+walltime = "00:30:00"
+years = "1985-1995",
diff --git a/tests/integration/generated/test_min_case_global_time_series_viewers_chrysalis.cfg b/tests/integration/generated/test_min_case_global_time_series_viewers_chrysalis.cfg
new file mode 100644
index 00000000..0bd495ae
--- /dev/null
+++ b/tests/integration/generated/test_min_case_global_time_series_viewers_chrysalis.cfg
@@ -0,0 +1,46 @@
+[default]
+case = "v3.LR.historical_0051"
+constraint = ""
+dry_run = "False"
+environment_commands = ""
+input = /lcrc/group/e3sm2/ac.wlin/E3SMv3/v3.LR.historical_0051
+input_subdir = archive/atm/hist
+mapping_file = "map_ne30pg2_to_cmip6_180x360_aave.20200201.nc"
+output = "/lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_viewers_output/unique_id/v3.LR.historical_0051"
+partition = "debug"
+qos = "regular"
+www = "/lcrc/group/e3sm/public_html/diagnostic_output/ac.forsyth2/zppy_min_case_global_time_series_viewers_www/unique_id"
+
+[ts]
+active = True
+e3sm_to_cmip_environment_commands = ""
+walltime = "00:30:00"
+years = "1985:1995:5",
+
+  [[ 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"
+  vars = "TREFHT"
+
+  [[ lnd_monthly_glb ]]
+  frequency = "monthly"
+  input_files = "elm.h0"
+  input_subdir = "archive/lnd/hist"
+  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.historical_0051"
+figstr = "v3.LR.historical_0051"
+ncols = 1
+nrows = 1
+plots_original=""
+plots_atm = "TREFHT"
+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 = 5
+walltime = "00:30:00"
+years = "1985-1995",
diff --git a/tests/integration/generated/update_bundles_expected_files_chrysalis.sh b/tests/integration/generated/update_bundles_expected_files_chrysalis.sh
new file mode 100755
index 00000000..72b1cf50
--- /dev/null
+++ b/tests/integration/generated/update_bundles_expected_files_chrysalis.sh
@@ -0,0 +1,20 @@
+# Run this script to update expected files for test_bundles.py
+# Run from the top level of the zppy repo
+# Run as `./tests/integration/generated/update_bundles_expected_files.sh`
+
+# Remove old expected files.
+rm -rf /lcrc/group/e3sm/public_html/zppy_test_resources/expected_bundles
+# Your output will now become the new expectation.
+# Copy output so you don't have to rerun zppy to generate the output.
+cp -r /lcrc/group/e3sm/public_html/diagnostic_output/ac.forsyth2/zppy_test_bundles_www/test-601-output-viewer-v3/v2.LR.historical_0201 /lcrc/group/e3sm/public_html/zppy_test_resources/expected_bundles
+mkdir -p /lcrc/group/e3sm/public_html/zppy_test_resources/expected_bundles/bundle_files
+cp -r /lcrc/group/e3sm/ac.forsyth2/zppy_test_bundles_output/test-601-output-viewer-v3/v2.LR.historical_0201/post/scripts/bundle*.bash /lcrc/group/e3sm/public_html/zppy_test_resources/expected_bundles/bundle_files
+zppy_top_level=$(pwd)
+cd /lcrc/group/e3sm/public_html/zppy_test_resources/expected_bundles
+# Remove the image check failures, so they don't end up in the expected files.
+rm -rf image_check_failures_bundles
+# This file will list all the expected images.
+find . -type f -name '*.png' > ../image_list_expected_bundles.txt
+cd ${zppy_top_level}
+# Rerun test
+python -u -m unittest tests/integration/test_bundles.py
diff --git a/tests/integration/generated/update_complete_run_expected_files_chrysalis.sh b/tests/integration/generated/update_complete_run_expected_files_chrysalis.sh
new file mode 100755
index 00000000..3d692376
--- /dev/null
+++ b/tests/integration/generated/update_complete_run_expected_files_chrysalis.sh
@@ -0,0 +1,18 @@
+# Run this script to update expected files for test_complete_run.py
+# Run from the top level of the zppy repo
+# Run as `./tests/integration/generated/update_complete_run_expected_files.sh`
+
+# Remove old expected files.
+rm -rf /lcrc/group/e3sm/public_html/zppy_test_resources/expected_complete_run
+# Your output will now become the new expectation.
+# Copy output so you don't have to rerun zppy to generate the output.
+cp -r /lcrc/group/e3sm/public_html/diagnostic_output/ac.forsyth2/zppy_test_complete_run_www/test-601-output-viewer-v3/v2.LR.historical_0201 /lcrc/group/e3sm/public_html/zppy_test_resources/expected_complete_run
+zppy_top_level=$(pwd)
+cd /lcrc/group/e3sm/public_html/zppy_test_resources/expected_complete_run
+# Remove the image check failures, so they don't end up in the expected files.
+rm -rf image_check_failures_complete_run
+# This file will list all the expected images.
+find . -type f -name '*.png' > ../image_list_expected_complete_run.txt
+cd ${zppy_top_level}
+# Rerun test
+python -u -m unittest tests/integration/test_complete_run.py
diff --git a/tests/integration/template_min_case_global_time_series_custom.cfg b/tests/integration/template_min_case_global_time_series_custom.cfg
index 538ac7ac..f35eb364 100644
--- a/tests/integration/template_min_case_global_time_series_custom.cfg
+++ b/tests/integration/template_min_case_global_time_series_custom.cfg
@@ -23,7 +23,7 @@ years = "1985:1995:5",
   input_files = "eam.h0"
   input_subdir = "archive/atm/hist"
   mapping_file = "glb"
-  vars = "TREFHT,AODDUST"
+  vars = "TREFHT"
 
   [[ lnd_monthly_glb ]]
   frequency = "monthly"
@@ -37,7 +37,7 @@ active = True
 experiment_name = "#expand case_name#"
 figstr = "#expand case_name#"
 plots_original=""
-plots_atm = "TREFHT,AODDUST"
+plots_atm = "TREFHT"
 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 = 5
 walltime = "00:30:00"
diff --git a/tests/integration/template_min_case_global_time_series_viewers.cfg b/tests/integration/template_min_case_global_time_series_viewers.cfg
new file mode 100644
index 00000000..2401bc42
--- /dev/null
+++ b/tests/integration/template_min_case_global_time_series_viewers.cfg
@@ -0,0 +1,46 @@
+[default]
+case = "#expand case_name#"
+constraint = "#expand constraint#"
+dry_run = "#expand dry_run#"
+environment_commands = "#expand environment_commands#"
+input = /lcrc/group/e3sm2/ac.wlin/E3SMv3/#expand case_name#
+input_subdir = archive/atm/hist
+mapping_file = "map_ne30pg2_to_cmip6_180x360_aave.20200201.nc"
+output = "#expand user_output#zppy_min_case_global_time_series_viewers_output/#expand unique_id#/#expand case_name#"
+partition = "#expand partition_short#"
+qos = "#expand qos_short#"
+www = "#expand user_www#zppy_min_case_global_time_series_viewers_www/#expand unique_id#"
+
+[ts]
+active = True
+e3sm_to_cmip_environment_commands = "#expand e3sm_to_cmip_environment_commands#"
+walltime = "00:30:00"
+years = "1985:1995:5",
+
+  [[ 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"
+  vars = "TREFHT"
+
+  [[ lnd_monthly_glb ]]
+  frequency = "monthly"
+  input_files = "elm.h0"
+  input_subdir = "archive/lnd/hist"
+  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 = "#expand case_name#"
+figstr = "#expand case_name#"
+ncols = 1
+nrows = 1
+plots_original=""
+plots_atm = "TREFHT"
+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 = 5
+walltime = "00:30:00"
+years = "1985-1995",
diff --git a/tests/integration/utils.py b/tests/integration/utils.py
index 2390d60d..a52ca578 100644
--- a/tests/integration/utils.py
+++ b/tests/integration/utils.py
@@ -297,6 +297,7 @@ def generate_cfgs(unified_testing=False, dry_run=False):
         "min_case_global_time_series_custom",
         "min_case_global_time_series_original_8_no_ocn",
         "min_case_global_time_series_original_8",
+        "min_case_global_time_series_viewers",
         "min_case_ilamb_land_only",
         "min_case_ilamb",
         "min_case_mpas_analysis",
diff --git a/zppy/global_time_series.py b/zppy/global_time_series.py
index 674c8187..dd74b58a 100644
--- a/zppy/global_time_series.py
+++ b/zppy/global_time_series.py
@@ -1,4 +1,5 @@
 import os
+import shutil
 from typing import Any, Dict, List
 
 from zppy.bundle import handle_bundles
@@ -49,13 +50,18 @@ def global_time_series(config, script_dir, existing_bundles, job_ids_file):
             c["global_time_series_dir"] = os.path.join(script_dir, f"{prefix}_dir")
             if not os.path.exists(c["global_time_series_dir"]):
                 os.mkdir(c["global_time_series_dir"])
-            scripts = ["coupled_global.py", "readTS.py", "ocean_month.py"]
+            scripts = ["coupled_global.py", "ocean_month.py"]
             for script in scripts:
                 script_template = template_env.get_template(script)
                 script_file = os.path.join(c["global_time_series_dir"], script)
                 with open(script_file, "w") as f:
                     f.write(script_template.render(**c))
                 make_executable(script_file)
+            resources = ["zppy_land_fields.csv"]
+            for resource in resources:
+                resource_file = os.path.join(c["global_time_series_dir"], resource)
+                shutil.copy(f"zppy/templates/{resource}", resource_file)
+
             # Create script
             with open(bash_file, "w") as f:
                 f.write(template.render(**c))
diff --git a/zppy/templates/coupled_global.py b/zppy/templates/coupled_global.py
index 9fc41401..79e9e518 100644
--- a/zppy/templates/coupled_global.py
+++ b/zppy/templates/coupled_global.py
@@ -1,9 +1,14 @@
-# Script to plot some global atmosphere and ocean time series
+# Script to generate global time series plots
+import csv
+import distutils.dir_util
 import glob
 import math
+import os
+import stat
 import sys
 import traceback
-from typing import Any, Dict, List, Optional, Tuple
+from enum import Enum
+from typing import Any, Dict, List, Tuple
 
 import cftime
 import matplotlib as mpl
@@ -11,12 +16,596 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import xarray
+import xcdat
+from bs4 import BeautifulSoup
 from netCDF4 import Dataset
-from readTS import TS
+from output_viewer.build import build_page, build_viewer
+from output_viewer.index import (
+    OutputFile,
+    OutputGroup,
+    OutputIndex,
+    OutputPage,
+    OutputRow,
+)
+from output_viewer.utils import rechmod
 
 mpl.use("Agg")
 
 
+# Useful classes and their helper functions ###################################
+def param_get_list(param_value: str) -> List[str]:
+    if param_value == "None":
+        return []
+    else:
+        return param_value.split(",")
+
+
+def get_region(rgn: str) -> str:
+    if rgn.lower() in ["glb", "global"]:
+        rgn = "glb"
+    elif rgn.lower() in ["n", "north", "northern"]:
+        rgn = "n"
+    elif rgn.lower() in ["s", "south", "southern"]:
+        rgn = "s"
+    else:
+        raise ValueError(f"Invalid rgn={rgn}")
+    return rgn
+
+
+class Parameters(object):
+    def __init__(self, parameters):
+        self.case_dir: str = parameters[1]
+        self.experiment_name: str = parameters[2]
+        self.figstr: str = parameters[3]
+        self.year1: int = int(parameters[4])
+        self.year2: int = int(parameters[5])
+        self.color: str = parameters[6]
+        self.ts_num_years_str: str = parameters[7]
+        self.plots_original: List[str] = param_get_list(parameters[8])
+        self.atmosphere_only: bool = (
+            False if (parameters[9].lower() == "false") else True
+        )
+        self.plots_atm: List[str] = param_get_list(parameters[10])
+        self.plots_ice: List[str] = param_get_list(parameters[11])
+        self.plots_lnd: List[str] = param_get_list(parameters[12])
+        self.plots_ocn: List[str] = param_get_list(parameters[13])
+        self.nrows: int = int(parameters[14])
+        self.ncols: int = int(parameters[15])
+        # These regions are used often as strings,
+        # so making an Enum Region={GLOBAL, NORTH, SOUTH} would be limiting.
+        self.regions: List[str] = list(
+            map(lambda rgn: get_region(rgn), parameters[16].split(","))
+        )
+
+
+class Metric(Enum):
+    AVERAGE = 1
+    TOTAL = 2
+
+
+class Variable(object):
+    def __init__(
+        self,
+        variable_name,
+        metric=Metric.AVERAGE,
+        scale_factor=1.0,
+        original_units="",
+        final_units="",
+        group="",
+        long_name="",
+    ):
+        # The name of the EAM/ELM/etc. variable on the monthly h0 history file
+        self.variable_name: str = variable_name
+
+        # These fields are used for computation
+        # Global average over land area or global total
+        self.metric: Metric = metric
+        # The factor that should convert from original_units to final_units, after standard processing with nco
+        self.scale_factor: float = scale_factor
+        # Test string for the units as given on the history file (included here for possible testing)
+        self.original_units: str = original_units
+        # The units that should be reported in time series plots, based on metric and scale_factor
+        self.final_units: str = final_units
+
+        # These fields are used for plotting
+        # A name used to cluster variables together, to be separated in groups within the output web pages
+        self.group: str = group
+        # Descriptive text to add to the plot page to help users identify the variable
+        self.long_name: str = long_name
+
+
+def get_vars_original(plots_original: List[str]) -> List[Variable]:
+    vars_original: List[Variable] = []
+    if ("net_toa_flux_restom" in plots_original) or (
+        "net_atm_energy_imbalance" in plots_original
+    ):
+        vars_original.append(Variable("RESTOM"))
+    if "net_atm_energy_imbalance" in plots_original:
+        vars_original.append(Variable("RESSURF"))
+    if "global_surface_air_temperature" in plots_original:
+        vars_original.append(Variable("TREFHT"))
+    if "toa_radiation" in plots_original:
+        vars_original.append(Variable("FSNTOA"))
+        vars_original.append(Variable("FLUT"))
+    if "net_atm_water_imbalance" in plots_original:
+        vars_original.append(Variable("PRECC"))
+        vars_original.append(Variable("PRECL"))
+        vars_original.append(Variable("QFLX"))
+    return vars_original
+
+
+def land_csv_row_to_var(csv_row: List[str]) -> Variable:
+    # “A” or “T” for global average over land area or global total, respectively
+    metric: Metric
+    if csv_row[1] == "A":
+        metric = Metric.AVERAGE
+    elif csv_row[1] == "T":
+        metric = Metric.TOTAL
+    else:
+        raise ValueError(f"Invalid metric={csv_row[1]}")
+    return Variable(
+        variable_name=csv_row[0],
+        metric=metric,
+        scale_factor=float(csv_row[2]),
+        original_units=csv_row[3],
+        final_units=csv_row[4],
+        group=csv_row[5],
+        long_name=csv_row[6],
+    )
+
+
+def construct_land_variables(requested_vars: List[str]) -> List[Variable]:
+    var_list: List[Variable] = []
+    header = True
+    # If this file is being run stand-alone, then
+    # it will search the directory above the git directory
+    with open("zppy_land_fields.csv", newline="") as csv_file:
+        print("In File")
+        var_reader = csv.reader(csv_file)
+        for row in var_reader:
+            print(f"row={row}")
+            # Skip the header row
+            if header:
+                header = False
+            else:
+                # If set to "all" then we want all variables.
+                # Design note: we can't simply run all variables if requested_vars is empty because
+                # that would actually mean the user doesn't want to make *any* land plots.
+                if (requested_vars == ["all"]) or (row[0] in requested_vars):
+                    row_elements_strip_whitespace: List[str] = list(
+                        map(lambda x: x.strip(), row)
+                    )
+                    var_list.append(land_csv_row_to_var(row_elements_strip_whitespace))
+    return var_list
+
+
+def construct_generic_variables(requested_vars: List[str]) -> List[Variable]:
+    var_list: List[Variable] = []
+    for var_name in requested_vars:
+        var_list.append(Variable(var_name))
+    return var_list
+
+
+class RequestedVariables(object):
+    def __init__(self, parameters: Parameters):
+        self.vars_original: List[Variable] = get_vars_original(
+            parameters.plots_original
+        )
+        self.vars_land: List[Variable] = construct_land_variables(parameters.plots_lnd)
+        self.vars_atm: List[Variable] = construct_generic_variables(
+            parameters.plots_atm
+        )
+        self.vars_ice: List[Variable] = construct_generic_variables(
+            parameters.plots_ice
+        )
+        self.vars_ocn: List[Variable] = construct_generic_variables(
+            parameters.plots_ocn
+        )
+
+
+class VariableGroup(object):
+    def __init__(self, name: str, variables: List[Variable]):
+        self.group_name = name
+        self.variables = variables
+
+
+class TS(object):
+    def __init__(self, directory):
+
+        self.directory: str = directory
+
+        # `directory` will be of the form `{case_dir}/post/<component>/glb/ts/monthly/{ts_num_years_str}yr/`
+        self.f: xarray.core.dataset.Dataset = xcdat.open_mfdataset(
+            f"{directory}*.nc", center_times=True
+        )
+
+    def __del__(self):
+
+        self.f.close()
+
+    def globalAnnualHelper(
+        self,
+        var: str,
+        metric: Metric,
+        scale_factor: float,
+        original_units: str,
+        final_units: str,
+    ) -> Tuple[xarray.core.dataarray.DataArray, str]:
+
+        data_array: xarray.core.dataarray.DataArray
+        units: str = ""
+
+        # Constants, from AMWG diagnostics
+        Lv = 2.501e6
+        Lf = 3.337e5
+
+        # Is this a derived variable?
+        if var == "RESTOM":
+            FSNT, _ = self.globalAnnualHelper(
+                "FSNT", metric, scale_factor, original_units, final_units
+            )
+            FLNT, _ = self.globalAnnualHelper(
+                "FLNT", metric, scale_factor, original_units, final_units
+            )
+            data_array = FSNT - FLNT
+        elif var == "RESTOA":
+            print("NOT READY")
+            FSNTOA, _ = self.globalAnnualHelper(
+                "FSNTOA", metric, scale_factor, original_units, final_units
+            )
+            FLUT, _ = self.globalAnnualHelper(
+                "FLUT", metric, scale_factor, original_units, final_units
+            )
+            data_array = FSNTOA - FLUT
+        elif var == "LHFLX":
+            QFLX, _ = self.globalAnnualHelper(
+                "QFLX", metric, scale_factor, original_units, final_units
+            )
+            PRECC, _ = self.globalAnnualHelper(
+                "PRECC", metric, scale_factor, original_units, final_units
+            )
+            PRECL, _ = self.globalAnnualHelper(
+                "PRECL", metric, scale_factor, original_units, final_units
+            )
+            PRECSC, _ = self.globalAnnualHelper(
+                "PRECSC", metric, scale_factor, original_units, final_units
+            )
+            PRECSL, _ = self.globalAnnualHelper(
+                "PRECSL", metric, scale_factor, original_units, final_units
+            )
+            data_array = (Lv + Lf) * QFLX - Lf * 1.0e3 * (
+                PRECC + PRECL - PRECSC - PRECSL
+            )
+        elif var == "RESSURF":
+            FSNS, _ = self.globalAnnualHelper(
+                "FSNS", metric, scale_factor, original_units, final_units
+            )
+            FLNS, _ = self.globalAnnualHelper(
+                "FLNS", metric, scale_factor, original_units, final_units
+            )
+            SHFLX, _ = self.globalAnnualHelper(
+                "SHFLX", metric, scale_factor, original_units, final_units
+            )
+            LHFLX, _ = self.globalAnnualHelper(
+                "LHFLX", metric, scale_factor, original_units, final_units
+            )
+            data_array = FSNS - FLNS - SHFLX - LHFLX
+        elif var == "PREC":
+            PRECC, _ = self.globalAnnualHelper(
+                "PRECC", metric, scale_factor, original_units, final_units
+            )
+            PRECL, _ = self.globalAnnualHelper(
+                "PRECL", metric, scale_factor, original_units, final_units
+            )
+            data_array = 1.0e3 * (PRECC + PRECL)
+        else:
+            # Non-derived variables
+            if (metric == Metric.AVERAGE) or (metric == Metric.TOTAL):
+                annual_average_dataset_for_var: xarray.core.dataset.Dataset = (
+                    self.f.temporal.group_average(var, "year")
+                )
+                data_array = annual_average_dataset_for_var.data_vars[var]
+            # elif metric == Metric.TOTAL:
+            #     # TODO: Implement this!
+            #     raise NotImplementedError()
+            else:
+                # This shouldn't be possible
+                raise ValueError(f"Invalid Enum option for metric={metric}")
+            units = data_array.units
+            # `units` will be "1" if it's a dimensionless quantity
+            if (units != "1") and (original_units != "") and original_units != units:
+                raise ValueError(
+                    f"Units don't match up: Have {units} but expected {original_units}. This renders the supplied scale_factor ({scale_factor}) unusable."
+                )
+            data_array *= scale_factor
+            units = final_units
+        return data_array, units
+
+    def globalAnnual(
+        self, var: Variable
+    ) -> Tuple[xarray.core.dataarray.DataArray, str]:
+        return self.globalAnnualHelper(
+            var.variable_name,
+            var.metric,
+            var.scale_factor,
+            var.original_units,
+            var.final_units,
+        )
+
+
+# Copied from e3sm_diags
+class OutputViewer(object):
+    def __init__(self, path=".", index_name="Results"):
+        self.path = os.path.abspath(path)
+        self.index = OutputIndex(index_name)
+        self.cache = {}  # dict of { OutputPage: { OutputGroup: [OutputRow] } }
+        self.page = None
+        self.group = None
+        self.row = None
+
+    def add_page(self, page_title, *args, **kwargs):
+        """Add a page to the viewer's index"""
+        self.page = OutputPage(page_title, *args, **kwargs)
+        self.cache[self.page] = {}
+        self.index.addPage(self.page)
+
+    def set_page(self, page_title):
+        """Sets the page with the title name as the current page"""
+        for output_page in self.cache:
+            if page_title == output_page.title:
+                self.page = output_page
+                return
+        raise RuntimeError("There is no page titled: %s" % page_title)
+
+    def add_group(self, group_name):
+        """Add a group to the current page"""
+        if self.page is None:
+            raise RuntimeError("You must first insert a page with add_page()")
+        self.group = OutputGroup(group_name)
+        if self.group not in self.cache[self.page]:
+            self.cache[self.page][self.group] = []  # group doesn't have any rows yet
+        self.page.addGroup(self.group)
+
+    def set_group(self, group_name):
+        """Sets the group with the title name as the current group"""
+        for output_group in self.cache[self.page]:
+            if group_name == output_group.title:
+                self.group = output_group
+                return
+        raise RuntimeError("There is no group titled: %s" % group_name)
+
+    def add_row(self, row_name):
+        """Add a row with the title name to the current group"""
+        if self.group is None:
+            raise RuntimeError("You must first insert a group with add_group()")
+        self.row = OutputRow(row_name, [])
+        if self.row not in self.cache[self.page][self.group]:
+            self.cache[self.page][self.group].append(self.row)
+        self.page.addRow(self.row, len(self.page.groups) - 1)  # type: ignore
+
+    def set_row(self, row_name):
+        """Sets the row with the title name as the current row"""
+        for output_row in self.cache[self.page][self.group]:
+            if row_name == output_row.title:
+                self.row = output_row
+                return
+        raise RuntimeError("There is no row titled: %s" % row_name)
+
+    def add_cols(self, cols):
+        """Add multiple string cols to the current row"""
+        self.row.columns.append(cols)  # type: ignore
+
+    def add_col(self, col, is_file=False, **kwargs):
+        """Add a single col to the current row. Set is_file to True if the col is a file path."""
+        if is_file:
+            self.row.columns.append(OutputFile(col, **kwargs))  # type: ignore
+        else:
+            self.row.columns.append(col)  # type: ignore
+
+    def generate_page(self) -> str:
+        """
+        Generate and return the location of the current HTML page.
+        """
+        self.index.toJSON(os.path.join(self.path, "index.json"))
+
+        default_mask = stat.S_IMODE(os.stat(self.path).st_mode)
+        rechmod(self.path, default_mask)
+
+        if os.access(self.path, os.W_OK):
+            default_mask = stat.S_IMODE(
+                os.stat(self.path).st_mode
+            )  # mode of files to be included
+            url = build_page(
+                self.page,
+                os.path.join(self.path, "index.json"),
+                default_mask=default_mask,
+            )
+            return url
+
+        raise RuntimeError("Error geneating the page.")
+
+    def generate_viewer(self):
+        """Generate the webpage"""
+        self.index.toJSON(os.path.join(self.path, "index.json"))
+
+        default_mask = stat.S_IMODE(os.stat(self.path).st_mode)
+        rechmod(self.path, default_mask)
+
+        if os.access(self.path, os.W_OK):
+            default_mask = stat.S_IMODE(
+                os.stat(self.path).st_mode
+            )  # mode of files to be included
+            build_viewer(
+                os.path.join(self.path, "index.json"),
+                diag_name="Global Time Series",
+                default_mask=default_mask,
+            )
+
+
+# Setup #######################################################################
+
+
+def get_data_dir(parameters: Parameters, component: str, conditional: bool) -> str:
+    return (
+        f"{parameters.case_dir}/post/{component}/glb/ts/monthly/{parameters.ts_num_years_str}yr/"
+        if conditional
+        else ""
+    )
+
+
+def get_exps(parameters: Parameters) -> List[Dict[str, Any]]:
+    # Experiments
+    use_atmos: bool = (parameters.plots_atm != []) or (parameters.plots_original != [])
+    # Use set intersection: check if any of these 3 plots were requested
+    set_intersection: set = set(["change_ohc", "max_moc", "change_sea_level"]) & set(
+        parameters.plots_original
+    )
+    has_original_ocn_plots: bool = set_intersection != set()
+    use_ocn: bool = (parameters.plots_ocn != []) or (
+        (not parameters.atmosphere_only) and has_original_ocn_plots
+    )
+    ocean_dir = get_data_dir(parameters, "ocn", use_ocn)
+    exps: List[Dict[str, Any]] = [
+        {
+            "atmos": get_data_dir(parameters, "atm", use_atmos),
+            "ice": get_data_dir(parameters, "ice", parameters.plots_ice != []),
+            "land": get_data_dir(parameters, "lnd", parameters.plots_lnd != []),
+            "ocean": ocean_dir,
+            "moc": ocean_dir,
+            "vol": ocean_dir,
+            "name": parameters.experiment_name,
+            "yoffset": 0.0,
+            "yr": ([parameters.year1, parameters.year2],),
+            "color": f"{parameters.color}",
+        }
+    ]
+    return exps
+
+
+def set_var(
+    exp: Dict[str, Any],
+    exp_key: str,
+    var_list: List[Variable],
+    valid_vars: List[str],
+    invalid_vars: List[str],
+    rgn: str,
+) -> None:
+    if exp[exp_key] != "":
+        ts: TS = TS(exp[exp_key])
+        for var in var_list:
+            var_str: str = var.variable_name
+            try:
+                data_array: xarray.core.dataarray.DataArray
+                units: str
+                data_array, units = ts.globalAnnual(var)
+                valid_vars.append(str(var_str))
+            except Exception as e:
+                print(e)
+                print(f"globalAnnual failed for {var_str}")
+                invalid_vars.append(str(var_str))
+                continue
+            if data_array.sizes["rgn"] > 1:
+                # number of years x 3 regions = data_array.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}")
+                data_array = data_array.isel(rgn=n)  # Just use nth region
+            elif rgn != "glb":
+                # data_array only has one dimension -- glb.
+                # Therefore it is not possible to get n or s plots.
+                raise RuntimeError(
+                    f"var={var_str} only has global data. Cannot process rgn={rgn}"
+                )
+            exp["annual"][var_str] = (data_array, units)
+            if "year" not in exp["annual"]:
+                years: np.ndarray[cftime.DatetimeNoLeap] = data_array.coords[
+                    "time"
+                ].values
+                exp["annual"]["year"] = [x.year for x in years]
+        del ts
+
+
+def process_data(
+    parameters: Parameters, requested_variables: RequestedVariables, rgn: str
+) -> List[Dict[str, Any]]:
+    exps: List[Dict[str, Any]] = get_exps(parameters)
+    valid_vars: List[str] = []
+    invalid_vars: List[str] = []
+    exp: Dict[str, Any]
+    for exp in exps:
+        exp["annual"] = {}
+
+        set_var(
+            exp,
+            "atmos",
+            requested_variables.vars_original,
+            valid_vars,
+            invalid_vars,
+            rgn,
+        )
+        set_var(
+            exp, "atmos", requested_variables.vars_atm, valid_vars, invalid_vars, rgn
+        )
+        set_var(exp, "ice", requested_variables.vars_ice, valid_vars, invalid_vars, rgn)
+        set_var(
+            exp, "land", requested_variables.vars_land, valid_vars, invalid_vars, rgn
+        )
+        set_var(
+            exp, "ocean", requested_variables.vars_ocn, valid_vars, invalid_vars, rgn
+        )
+
+        # Optionally read ohc
+        if exp["ocean"] != "":
+            ts = TS(exp["ocean"])
+            exp["annual"]["ohc"], _ = ts.globalAnnual(Variable("ohc"))
+            # anomalies with respect to first year
+            exp["annual"]["ohc"][:] = exp["annual"]["ohc"][:] - exp["annual"]["ohc"][0]
+
+        if exp["vol"] != "":
+            ts = TS(exp["vol"])
+            exp["annual"]["volume"], _ = ts.globalAnnual(Variable("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}"
+    )
+    return exps
+
+
+# Plotting ####################################################################
+
+
+def get_variable_groups(variables: List[Variable]) -> List[VariableGroup]:
+    group_names: List[str] = []
+    groups: List[VariableGroup] = []
+    for v in variables:
+        g: str = v.group
+        if g not in group_names:
+            # A new group!
+            group_names.append(g)
+            groups.append(VariableGroup(g, [v]))
+        else:
+            # Add a new variable to this existing group
+            for group in groups:
+                if g == group.group_name:
+                    group.variables.append(v)
+    return groups
+
+
 # ---additional function to get moc time series
 def getmoc(dir_in):
     files = sorted(glob.glob(dir_in + "mocTimeSeries*.nc"))
@@ -61,8 +650,6 @@ def add_line(year, var, year1, year2, ax, format="%4.2f", lw=1, color="b"):
     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
-
 
 # -----------------------------------------------------------------------------
 # Function to add line showing linear trend over a specified period
@@ -450,13 +1037,13 @@ def plot(ax, xlim, exps, param_dict, rgn):  # noqa: C901
     extreme_values = []
     for exp in exps:
         # Relevant to "Plot 5: plot_change_ohc"
-        if param_dict["check_exp_ocean"] and (exp["ocean"] is None):
+        if param_dict["check_exp_ocean"] and (exp["ocean"] == ""):
             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):
+        if param_dict["check_exp_vol"] and (exp["vol"] == ""):
             continue
         # Relevant to "Plot 6: plot_max_moc"
         if param_dict["use_getmoc"]:
@@ -545,83 +1132,39 @@ def plot(ax, xlim, exps, param_dict, rgn):  # noqa: C901
 }
 
 
-def param_get_list(param_value):
-    if param_value == "None":
-        return []
-    else:
-        return param_value.split(",")
-
-
-def set_var(
-    exp: Dict[str, Any],
-    exp_key: str,
-    var_list: List[str],
-    valid_vars: List[str],
-    invalid_vars: List[str],
-    rgn: str,
-) -> None:
-    if exp[exp_key] is not None:
-        ts = TS(exp[exp_key])
-        for var in var_list:
-            try:
-                v: xarray.core.dataarray.DataArray
-                units: Optional[str]
-                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 v.sizes["rgn"] > 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.isel(rgn=n)  # Just use nth region
-            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, units)
-            if "year" not in exp["annual"]:
-                years: np.ndarray[cftime.DatetimeNoLeap] = v.coords["time"].values
-                exp["annual"]["year"] = [x.year for x in years]
-        del ts
-
-
-def make_plot_pdfs(
-    figstr, rgn, component, xlim, exps, plot_list, valid_plots, invalid_plots
+# FIXME: C901 'make_plot_pdfs' is too complex (20)
+def make_plot_pdfs(  # noqa: C901
+    parameters: Parameters,
+    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
+    plots_per_page = parameters.nrows * parameters.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")
+    pdf = matplotlib.backends.backend_pdf.PdfPages(
+        f"{parameters.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}")
+        if plots_per_page == 1:
+            fig = plt.figure(1, figsize=[13.5 / 2, 16.5 / 4])
+        else:
+            fig = plt.figure(1, figsize=[13.5, 16.5])
+        fig.suptitle(f"{parameters.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)
+            ax = plt.subplot(parameters.nrows, parameters.ncols, j + 1)
             if component == "original":
                 try:
                     plot_function = PLOT_DICT[plot_list[counter]]
@@ -662,210 +1205,194 @@ def make_plot_pdfs(
 
         fig.tight_layout()
         pdf.savefig(1)
-        if num_pages > 1:
-            fig.savefig(f"{figstr}_{rgn}_{component}_{page}.png", dpi=150)
+        if plots_per_page == 1:
+            fig.savefig(
+                f"{parameters.figstr}_{rgn}_{component}_{plot_name}.png", dpi=150
+            )
+        elif num_pages > 1:
+            fig.savefig(f"{parameters.figstr}_{rgn}_{component}_{page}.png", dpi=150)
         else:
-            fig.savefig(f"{figstr}_{rgn}_{component}.png", dpi=150)
+            fig.savefig(f"{parameters.figstr}_{rgn}_{component}.png", dpi=150)
         plt.clf()
     pdf.close()
 
 
+# Run coupled_global ##########################################################
 # -----------------------------------------------------------------------------
-# FIXME: C901 'run' is too complex (19)
-def run(parameters, rgn):  # noqa: C901
-    # These are the "Tableau 20" colors as RGB.
-    t20: List[Tuple[float, float, float]] = [
-        (31, 119, 180),
-        (174, 199, 232),
-        (255, 127, 14),
-        (255, 187, 120),
-        (44, 160, 44),
-        (152, 223, 138),
-        (214, 39, 40),
-        (255, 152, 150),
-        (148, 103, 189),
-        (197, 176, 213),
-        (140, 86, 75),
-        (196, 156, 148),
-        (227, 119, 194),
-        (247, 182, 210),
-        (127, 127, 127),
-        (199, 199, 199),
-        (188, 189, 34),
-        (219, 219, 141),
-        (23, 190, 207),
-        (158, 218, 229),
-    ]
-    # Scale the RGB values to the [0, 1] range, which is the format matplotlib accepts.
-    for i in range(len(t20)):
-        r, g, b = t20[i]
-        t20[i] = (r / 255.0, g / 255.0, b / 255.0)
-
-    # "Tableau 10" uses every other color
-    t10 = []
-    for i in range(0, len(t20), 2):
-        t10.append(t20[i])
-
-    # -----------------------------------------------------------------------------
-    # --- Atmos data ---
-
+def run(parameters: Parameters, requested_variables: RequestedVariables, rgn: str):
     # Experiments
-    case_dir = parameters[1]
-    experiment_name = parameters[2]
-    figstr = parameters[3]
-    year1 = int(parameters[4])
-    year2 = int(parameters[5])
-    color = parameters[6]
-    ts_num_years = parameters[7]
-    plots_original = param_get_list(parameters[8])
-    if parameters[9].lower() == "false":
-        atmosphere_only = False
-    else:
-        atmosphere_only = True
-    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: List[Dict[str, Any]] = [
-        {
-            "atmos": (
-                f"{case_dir}/post/atm/glb/ts/monthly/{ts_num_years}yr/"
-                if use_atmos
-                else None
-            ),
-            "ice": (
-                f"{case_dir}/post/ice/glb/ts/monthly/{ts_num_years}yr/"
-                if plots_ice
-                else None
-            ),
-            "land": (
-                f"{case_dir}/post/lnd/glb/ts/monthly/{ts_num_years}yr/"
-                if plots_lnd
-                else None
-            ),
-            "ocean": (
-                f"{case_dir}/post/ocn/glb/ts/monthly/{ts_num_years}yr/"
-                if use_ocn
-                else None
-            ),
-            "moc": (
-                f"{case_dir}/post/ocn/glb/ts/monthly/{ts_num_years}yr/"
-                if use_ocn
-                else None
-            ),
-            "vol": (
-                f"{case_dir}/post/ocn/glb/ts/monthly/{ts_num_years}yr/"
-                if use_ocn
-                else None
-            ),
-            "name": experiment_name,
-            "yoffset": 0.0,
-            "yr": ([year1, year2],),
-            "color": f"{color}",
-        }
-    ]
-
-    valid_vars: List[str] = []
-    invalid_vars: List[str] = []
-
-    # Read data
-    exp: Dict[str, Any]
-    for exp in exps:
-        exp["annual"] = {}
+    exps: List[Dict[str, Any]] = process_data(parameters, requested_variables, rgn)
 
-        # 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")
-            # 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")
-            # 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)]
+    xlim: List[float] = [float(parameters.year1), float(parameters.year2)]
 
     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
-    )
-
+    # Use list of tuples rather than a dict, to keep order
+    mapping: List[Tuple[str, List[str]]] = [
+        ("original", parameters.plots_original),
+        ("atm", parameters.plots_atm),
+        ("ice", parameters.plots_ice),
+        ("lnd", parameters.plots_lnd),
+        ("ocn", parameters.plots_ocn),
+    ]
+    for component, plot_list in mapping:
+        make_plot_pdfs(
+            parameters,
+            rgn,
+            component,
+            xlim,
+            exps,
+            plot_list,
+            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}"
+        f"These {rgn} regions plots could not be generated successfully: {invalid_plots}"
     )
 
 
-def run_by_region(parameters):
-    regions = parameters[14].split(",")
-    for rgn in regions:
-        if rgn.lower() in ["glb", "global"]:
-            rgn = "glb"
-        elif rgn.lower() in ["n", "north", "northern"]:
-            rgn = "n"
-        elif rgn.lower() in ["s", "south", "southern"]:
-            rgn = "s"
+def get_vars(requested_variables: RequestedVariables, component: str) -> List[Variable]:
+    vars: List[Variable]
+    if component == "original":
+        vars = requested_variables.vars_original
+    elif component == "atm":
+        vars = requested_variables.vars_atm
+    elif component == "ice":
+        vars = requested_variables.vars_ice
+    elif component == "lnd":
+        vars = requested_variables.vars_land
+    elif component == "ocn":
+        vars = requested_variables.vars_ocn
+    else:
+        raise ValueError(f"Invalid component={component}")
+    return vars
+
+
+def create_viewer(parameters: Parameters, vars: List[Variable], component: str) -> str:
+    viewer = OutputViewer(path=".")
+    viewer.add_page("Table", parameters.regions)
+    groups: List[VariableGroup] = get_variable_groups(vars)
+    for group in groups:
+        # Only groups that have at least one variable will be returned by `get_variable_groups`
+        # So, we know this group will be non-empty and should therefore be added to the viewer.
+        viewer.add_group(group.group_name)
+        for var in group.variables:
+            plot_name: str = var.variable_name
+            row_title: str
+            if var.long_name != "":
+                row_title = f"{plot_name}: {var.long_name}"
+            else:
+                row_title = plot_name
+            viewer.add_row(row_title)
+            for rgn in parameters.regions:
+                # v3.LR.historical_0051_glb_lnd_SOIL4C.png
+                # viewer/c-state/glb_lnd_soil4c.html
+                viewer.add_col(
+                    f"{parameters.figstr}_{rgn}_{component}_{plot_name}.png",
+                    is_file=True,
+                    title=f"{rgn}_{component}_{plot_name}",
+                )
+
+    url = viewer.generate_page()
+    viewer.generate_viewer()
+    # Copy the contents of `table` into the `viewer` directory
+    # (which initially only has `css` and `js` subdirectories)
+    # Because `viewer` already exists,
+    # `shutil.copytree` will not work.
+    distutils.dir_util.copy_tree("table", "viewer")
+    print(
+        os.getcwd()
+    )  # /lcrc/group/e3sm/ac.forsyth2/zppy_min_case_global_time_series_viewers_output/test-pr616-20241022v2/v3.LR.historical_0051/post/scripts/global_time_series_1985-1995_dir
+    # shutil.rmtree("table")
+    # new_url = f"viewer_{component}"
+    # # shutil.move("viewer", new_url)
+    # distutils.dir_util.copy_tree("viewer", new_url)
+    return url
+
+
+# Copied from E3SM Diags and modified
+def create_viewer_index(
+    root_dir: str, title_and_url_list: List[Tuple[str, str]]
+) -> str:
+    """
+    Creates the index page in root_dir which
+    joins the individual viewers.
+
+    Each tuple is on its own row.
+    """
+
+    def insert_data_in_row(row_obj, name, url):
+        """
+        Given a row object, insert the name and url.
+        """
+        td = soup.new_tag("td")
+        a = soup.new_tag("a")
+        a["href"] = url
+        a.string = name
+        td.append(a)
+        row_obj.append(td)
+
+    install_path = ""  # TODO: figure this out
+    path = os.path.join(install_path, "viewer", "index_template.html")
+    output = os.path.join(root_dir, "index.html")
+
+    soup = BeautifulSoup(open(path), "lxml")
+
+    # If no one changes it, the template only has
+    # one element in the find command below.
+    table = soup.find_all("table", {"class": "table"})[0]
+
+    # Adding the title.
+    tr = soup.new_tag("tr")
+    th = soup.new_tag("th")
+    th.string = "Output Sets"
+    tr.append(th)
+
+    # Adding each of the rows.
+    for row in title_and_url_list:
+        tr = soup.new_tag("tr")
+
+        if isinstance(row, list):
+            for elt in row:
+                name, url = elt
+                insert_data_in_row(tr, name, url)
         else:
-            raise RuntimeError(f"Invalid rgn={rgn}")
-        run(parameters, rgn)
+            name, url = row
+            insert_data_in_row(tr, name, url)
+
+        table.append(tr)
+
+    html = soup.prettify("utf-8")
+
+    with open(output, "wb") as f:
+        f.write(html)
+
+    return output
+
+
+def run_by_region(command_line_arguments):
+    parameters: Parameters = Parameters(command_line_arguments)
+    requested_variables = RequestedVariables(parameters)
+    for rgn in parameters.regions:
+        run(parameters, requested_variables, rgn)
+    plots_per_page = parameters.nrows * parameters.ncols
+    # TODO: Is this how we want to determine when to make a viewer or should we have a `make_viewer` parameter in the cfg?
+    if plots_per_page == 1:
+        # In this case, we don't want the summary PDF.
+        # Rather, we want to construct a viewer similar to E3SM Diags.
+        # TODO: determine directory paths for each viewer
+        # TODO: include "original"?
+        # for component in ["original", "atm", "ice", "lnd", "ocn"]:
+        title_and_url_list: List[Tuple[str, str]] = []
+        for component in ["lnd"]:
+            vars = get_vars(requested_variables, component)
+            url = create_viewer(parameters, vars, component)
+            print(url)
+            title_and_url_list.append((component, url))
+        # index_url: str = create_viewer_index(parameters.case_dir, title_and_url_list)
+        # print(f"Viewer index URL: {index_url}")
 
 
 if __name__ == "__main__":
diff --git a/zppy/templates/default.ini b/zppy/templates/default.ini
index 9bf5bdce..33043bd9 100644
--- a/zppy/templates/default.ini
+++ b/zppy/templates/default.ini
@@ -301,6 +301,10 @@ figstr = string(default="")
 # NOTE: always overrides value in [default]
 input_subdir = string(default="archive/ocn/hist")
 moc_file = string(default="")
+# Number of columns per page
+ncols = integer(default=2)
+# Number of rows per page
+nrows = integer(default=4)
 # 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.
@@ -319,6 +323,7 @@ plots_original = string(default="net_toa_flux_restom,global_surface_air_temperat
 # These should be a subset of the `vars` generated by the `ts` `glb` subtasks.
 plots_atm = string(default="")
 plots_ice = string(default="")
+# Set `plots_lnd = "all"` to run every variable in the land csv file.
 plots_lnd = string(default="")
 plots_ocn = string(default="")
 # regions to plot: glb, n, s (global, northern hemisphere, southern hemisphere)
diff --git a/zppy/templates/global_time_series.bash b/zppy/templates/global_time_series.bash
index 22dee4d7..5b0bdcc8 100644
--- a/zppy/templates/global_time_series.bash
+++ b/zppy/templates/global_time_series.bash
@@ -68,7 +68,7 @@ fi
 echo 'Update time series figures'
 cd ${global_ts_dir}
 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 }}
+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 }} {{ nrows }} {{ ncols }} {{ regions }}
 if [ $? != 0 ]; then
   cd {{ scriptDir }}
   echo 'ERROR (6)' > {{ prefix }}.status
@@ -82,6 +82,7 @@ results_dir_absolute_path={{ scriptDir }}/${results_dir}
 mkdir -p ${results_dir_absolute_path}
 cp *.pdf ${results_dir_absolute_path}
 cp *.png ${results_dir_absolute_path}
+cp -r viewer ${results_dir_absolute_path}
 
 ################################################################################
 # Copy output to web server
diff --git a/zppy/templates/readTS.py b/zppy/templates/readTS.py
deleted file mode 100644
index 16a9031e..00000000
--- a/zppy/templates/readTS.py
+++ /dev/null
@@ -1,78 +0,0 @@
-from typing import Optional, Tuple
-
-import xarray
-import xcdat  # noqa: F401
-
-
-class TS(object):
-    def __init__(self, directory):
-
-        self.directory: str = directory
-
-        # `directory` will be of the form `{case_dir}/post/<componen>/glb/ts/monthly/{ts_num_years}yr/`
-        self.f: xarray.core.dataset.Dataset = xcdat.open_mfdataset(
-            f"{directory}*.nc", center_times=True
-        )
-
-    def __del__(self):
-
-        self.f.close()
-
-    def globalAnnual(
-        self, var: str
-    ) -> Tuple[xarray.core.dataarray.DataArray, Optional[str]]:
-
-        v: xarray.core.dataarray.DataArray
-        units: Optional[str] = None
-
-        # Constants, from AMWG diagnostics
-        Lv = 2.501e6
-        Lf = 3.337e5
-
-        # Is this a derived variable?
-        if var == "RESTOM":
-
-            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")
-            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")
-            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")
-            v = FSNS - FLNS - SHFLX - LHFLX
-
-        elif var == "PREC":
-
-            PRECC, _ = self.globalAnnual("PRECC")
-            PRECL, _ = self.globalAnnual("PRECL")
-            v = 1.0e3 * (PRECC + PRECL)
-
-        else:
-            # Non-derived variables
-
-            annual_average_dataset_for_var: xarray.core.dataset.Dataset = (
-                self.f.temporal.group_average(var, "year")
-            )
-            v = annual_average_dataset_for_var.data_vars[var]
-            units = v.units
-
-        return v, units
diff --git a/zppy/templates/zppy_land_fields.csv b/zppy/templates/zppy_land_fields.csv
new file mode 100644
index 00000000..c220e70f
--- /dev/null
+++ b/zppy/templates/zppy_land_fields.csv
@@ -0,0 +1,354 @@
+Variable,Average or Total (A/T),Scale Factor,Original Units,Final Units,Group,Long name
+BCDEP,A,1.00000E+00, kg/m^2/s, kg/m^2/s,Aerosol Flux, total black carbon deposition (dry+wet) from atmosphere
+DSTDEP,A,1.00000E+00, kg/m^2/s, kg/m^2/s,Aerosol Flux, total dust deposition (dry+wet) from atmosphere
+DSTFLXT,A,1.00000E+00, kg/m2/s, kg/m^2/s,Aerosol Flux, total surface dust emission
+OCDEP,A,1.00000E+00, kg/m^2/s, kg/m^2/s,Aerosol Flux, total OC deposition (dry+wet) from atmosphere
+SNOBCMCL,A,1.00000E+00, kg/m2, kg/m^2,Aerosol State, mass of black carbon in snow column
+SNOBCMSL,A,1.00000E+00, kg/m2, kg/m^2,Aerosol State, mass of black carbon in top snow layer
+SNODSTMCL,A,1.00000E+00, kg/m2, kg/m^2,Aerosol State, mass of dust in snow column
+SNODSTMSL,A,1.00000E+00, kg/m2, kg/m^2,Aerosol State, mass of dust in top snow layer
+SNOOCMCL,A,1.00000E+00, kg/m2, kg/m^2,Aerosol State, mass of OC in snow column
+SNOOCMSL,A,1.00000E+00, kg/m2, kg/m^2,Aerosol State, mass of OC in top snow layer
+AGNPP,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, aboveground NPP
+AR,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, autotrophic respiration (maintenance + growth)
+BGNPP,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, belowground NPP
+CH4PROD,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, gridcell total production of CH4
+DWT_CONV_CFLUX_GRC,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, land conversion C flux
+DWT_SLASH_CFLUX,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, slash C flux to litter and CWD due to land use
+ER,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, total ecosystem respiration (autotrophic + heterotrophic)
+FCH4,T,3.15360E-02, kgC/m2/s, PgC/yr,C Flux, gridcell surface CH4 flux to atmosphere (+ to atm)
+FCH4TOCO2,T,3.15360E-02, gC/m2/s, PgC/yr,C Flux, gridcell oxidation of CH4 to CO2
+FPSN,A,1.00000E+00, umol/m2s, umol/m2s,C Flux, photosynthesis
+FPSN_WC,A,1.00000E+00, umol/m2s, umol/m2s,C Flux, Rubisco-limited photosynthesis
+FPSN_WJ,A,1.00000E+00, umol/m2s, umol/m2s,C Flux, RuBP-limited photosynthesis
+FPSN_WP,A,1.00000E+00, umol/m2s, umol/m2s,C Flux, product-limited photosynthesis
+FROOTC_ALLOC,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, fine root C allocation
+GPP,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, gross primary production
+GR,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, total growth respiration
+HR,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, total heterotrophic respiration
+LAND_USE_FLUX,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, total C emitted from land cover conversion and wood product pools
+LEAFC_ALLOC,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, leaf C allocation
+LITFALL,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, litterfall (leaves and fine roots)
+LITHR,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, litter heterotrophic respiration
+LITTERC_HR,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, litter C heterotrophic respiration
+LITTERC_LOSS,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, litter C loss
+MR,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, maintenance respiration
+NBP,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, net biome production (positive for sink)
+NEE,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, net ecosystem exchange of carbon (positive for source)
+NEP,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, net ecosystem production (no fire/landuse/harvest flux: positive for sink
+NPP,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, net primary production
+PSNSHA,A,1.00000E+00, umolCO2/m^2/s, umolCO2/m^2/s,C Flux, shaded leaf photosynthesis
+PSNSHADE_TO_CPOOL,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, C fixation from shaded canopy
+PSNSUN,A,1.00000E+00, umolCO2/m^2/s, umolCO2/m^2/s,C Flux, sunlit leaf photosynthesis
+PSNSUN_TO_CPOOL,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, C fixation from sunlit canopy
+RR,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, root respiration (fine root MR + total root GR)
+SOILC_HR,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, soil C heterotrophic respiration
+SOILC_LOSS,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, soil C loss
+SOMHR,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, soil organic matter heterotrophic respiration
+SOM_C_LEACHED,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, total flux of C from SOM pools due to leaching
+SR,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, total soil respiration (HR + root resp)
+WOODC_ALLOC,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, wood C allocation
+WOOD_HARVESTC,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, wood harvest C (to product pools)
+XR,T,3.15360E-02, gC/m^2/s, PgC/yr,C Flux, total excess respiration
+CPOOL,T,1.00000E-09, gC/m^2, PgC,C State, temporary photosynthate C pool
+CWDC,T,1.00000E-09, gC/m^2, PgC,C State, coarse woody debris C
+DEADCROOTC,T,1.00000E-09, gC/m^2, PgC,C State, dead coarse root C
+DEADSTEMC,T,1.00000E-09, gC/m^2, PgC,C State, dead stem C
+DISPVEGC,T,1.00000E-09, gC/m^2, PgC,C State, displayed vegetation carbon (excluding storage and cpool)
+FROOTC,T,1.00000E-09, gC/m^2, PgC,C State, fine root C
+LEAFC,T,1.00000E-09, gC/m^2, PgC,C State, leaf C
+LITR1C,T,1.00000E-09, gC/m^2, PgC,C State, LITR1 C
+LITR2C,T,1.00000E-09, gC/m^2, PgC,C State, LITR2 C
+LITR3C,T,1.00000E-09, gC/m^2, PgC,C State, LITR3 C
+LITTERC,T,1.00000E-09, gC/m^2, PgC,C State, litter C
+LIVECROOTC,T,1.00000E-09, gC/m^2, PgC,C State, live coarse root C
+LIVESTEMC,T,1.00000E-09, gC/m^2, PgC,C State, live stem C
+SEEDC_GRC,T,1.00000E-09, gC/m^2, PgC,C State, pool for seeding new PFTs via dynamic landcover
+SOIL1C,T,1.00000E-09, gC/m^2, PgC,C State, SOIL1 C
+SOIL2C,T,1.00000E-09, gC/m^2, PgC,C State, SOIL2 C
+SOIL3C,T,1.00000E-09, gC/m^2, PgC,C State, SOIL3 C
+SOIL4C,T,1.00000E-09, gC/m^2, PgC,C State, SOIL4 C
+SOILC,T,1.00000E-09, gC/m^2, PgC,C State, soil C
+STORVEGC,T,1.00000E-09, gC/m^2, PgC,C State, stored vegetation carbon (excluding cpool)
+TOTCOLC,T,1.00000E-09, gC/m^2, PgC,C State, total column-level C (no product pools)
+TOTECOSYSC,T,1.00000E-09, gC/m^2, PgC,C State, total ecosystem C (no product pools)
+TOTLITC,T,1.00000E-09, gC/m^2, PgC,C State, total litter carbon
+TOTLITC_1m,T,1.00000E-09, gC/m^2, PgC,C State, total litter C to 1 meter
+TOTPFTC,T,1.00000E-09, gC/m^2, PgC,C State, total PFT-level C (with cpool)
+TOTPRODC,T,1.00000E-09, gC/m^2, PgC,C State, total wood product C
+TOTSOMC,T,1.00000E-09, gC/m^2, PgC,C State, total soil organic matter C
+TOTSOMC_1m,T,1.00000E-09, gC/m^2, PgC,C State, total soil organic matter C to 1 meter
+TOTVEGC,T,1.00000E-09, gC/m^2, PgC,C State, total vegetation C (no cpool)
+TOTVEGC_ABG,T,1.00000E-09, gC/m^2, PgC,C State, total aboveground vegetation C (no cpool)
+WOODC,T,1.00000E-09, gC/m^2, PgC,C State, wood C
+XSMRPOOL,T,1.00000E-09, gC/m^2, PgC,C State, temporary photosynthate C pool
+EFLX_DYNBAL,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, dynamic land cover change conversion energy flux
+EFLX_GRND_LAKE,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, net heat flux into lake/snow surface (excluding light transmission)
+EFLX_LH_TOT,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, total latent heat flux [+ to atm]
+EFLX_LH_TOT_R,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, rural total evaporation
+ERRSEB,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, surface energy conservation error
+ERRSOI,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, soil/lake energy conservation error
+ERRSOL,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, solar radiation conservation error
+FCEV,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, canopy evaporation
+FCTR,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, canopy transpiration
+FGEV,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, ground evaporation
+FGR,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, heat flux into soil/snow including snow melt
+FGR12,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, heat flux between soil layers 1 and 2
+FGR_R,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, rural heat flux into soil/snow including snow melt
+FIRA,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, net longwave radiation
+FIRA_R,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, rural net longwave radiation
+FIRE,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, emitted longwave radiation
+FIRE_R,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, rural emitted longwave radiation
+FLDS,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, atmospheric longwave radiation
+FSA,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, absorbed solar radiation
+FSA_R,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, rural absorbed solar radiation
+FSDS,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, atmospheric incident solar radiation
+FSDSND,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, direct nir incident solar radiation
+FSDSNDLN,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, direct nir incident solar radiation at local noon
+FSDSNI,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, diffuse nir incident solar radiation
+FSDSVD,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, direct vis incident solar radiation
+FSDSVDLN,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, direct vis incident solar radiation at local noon
+FSDSVI,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, diffuse vis incident solar radiation
+FSDSVILN,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, diffuse vis incident solar radiation at local noon
+FSH,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, sensible heat
+FSH_G,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, sensible heat from ground
+FSH_NODYNLNDUSE,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, sensible heat not including correction for land use change
+FSH_R,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, rural sensible heat
+FSH_V,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, sensible heat from veg
+FSM,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, snow melt heat flux
+FSM_R,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, rural snow melt heat flux
+FSR,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, reflected solar radiation
+FSRND,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, direct nir reflected solar radiation
+FSRNDLN,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, direct nir reflected solar radiation at local noon
+FSRNI,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, diffuse nir reflected solar radiation
+FSRVD,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, direct vis reflected solar radiation
+FSRVDLN,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, direct vis reflected solar radiation at local noon
+FSRVI,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, diffuse vis reflected solar radiation
+PARVEGLN,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, absorbed par by vegetation at local noon
+SABG,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, solar rad absorbed by ground
+SABG_PEN,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, rural solar rad penetrating top soil or snow layer
+SABV,A,1.00000E+00, W/m^2, W/m^2,Energy Flux, solar rad absorbed by veg
+GC_HEAT1,A,1.00000E+00, J/m^2, J/m^2,Energy State, initial gridcell total heat content
+HC,A,1.00000E+00, MJ/m2, MJ/m2,Energy State, heat content of soil/snow/lake
+HCSOI,A,1.00000E+00, MJ/m2, MJ/m2,Energy State, soil heat content
+TBOT,A,1.00000E+00, K, K,Energy State, atmospheric air temperature
+TG,A,1.00000E+00, K, K,Energy State, ground temperature
+TG_R,A,1.00000E+00, K, K,Energy State, rural ground temperature
+TH2OSFC,A,1.00000E+00, K, K,Energy State, surface water temperature
+THBOT,A,1.00000E+00, K, K,Energy State, atmospheric air potential temperature
+TREFMNAV,A,1.00000E+00, K, K,Energy State, daily minimum of average 2-m temperature
+TREFMNAV_R,A,1.00000E+00, K, K,Energy State, rural daily minimum of average 2-m temperature
+TREFMXAV,A,1.00000E+00, K, K,Energy State, daily maximum of average 2-m temperature
+TREFMXAV_R,A,1.00000E+00, K, K,Energy State, rural daily maximum of average 2-m temperature
+TSA,A,1.00000E+00, K, K,Energy State, 2m air temperature
+TSA_R,A,1.00000E+00, K, K,Energy State, rural 2m air temperature
+TSOI_10CM,A,1.00000E+00, K, K,Energy State, soil temperature in top 10cm of soil
+TV,A,1.00000E+00, K, K,Energy State, vegetation temperature
+COL_FIRE_CLOSS,T,3.15360E-02, gC/m^2/s, PgC/yr,Fire, total column-level fire C loss (non-peat and non-conversion)
+FAREA_BURNED,A,1.00000E+00, proportion, proportion,Fire, timestep fractional area burned
+NFIRE,A,1.00000E+00, counts/km2/sec, counts/km2/sec,Fire, fire counts valid only in Reg.C
+PFT_FIRE_CLOSS,T,3.15360E-02, gC/m^2/s, PgC/yr,Fire, total PFT-level fire C loss (non-peat and non-conversion)
+DEFICIT,A,8.64000E+04, mm/s, mm/day,H2O Flux, runoff supply deficit
+DWB,A,8.64000E+04, mm/s, mm/day,H2O Flux, net change in total water mass
+QCHARGE,A,8.64000E+04, mm/s, mm/day,H2O Flux, aquifer recharge rate (vegetated landunits only)
+QDRAI,A,8.64000E+04, mm/s, mm/day,H2O Flux, sub-surface drainage
+QDRAI_PERCH,A,8.64000E+04, mm/s, mm/day,H2O Flux, perched water table drainage
+QDRAI_XS,A,8.64000E+04, mm/s, mm/day,H2O Flux, saturation excess drainage
+QDRIP,A,8.64000E+04, mm/s, mm/day,H2O Flux, throughfall
+QFLOOD,A,8.64000E+04, mm/s, mm/day,H2O Flux, runoff from river flooding
+QFLX_ICE_DYNBAL,A,8.64000E+04, mm/s, mm/day,H2O Flux, ice dynamic land cover change conversion runoff flux
+QFLX_LIQ_DYNBAL,A,8.64000E+04, mm/s, mm/day,H2O Flux, liq dynamic land cover change conversion runoff flux
+QH2OSFC,A,8.64000E+04, mm/s, mm/day,H2O Flux, surface water runoff
+QINFL,A,8.64000E+04, mm/s, mm/day,H2O Flux, infiltration
+QINTR,A,8.64000E+04, mm/s, mm/day,H2O Flux, interception
+QIRRIG_GRND,A,8.64000E+04, mm/s, mm/day,H2O Flux, groundwater irrigation
+QIRRIG_ORIG,A,8.64000E+04, mm/s, mm/day,H2O Flux, original total irrigation water demand (surface + ground)
+QIRRIG_REAL,A,8.64000E+04, mm/s, mm/day,H2O Flux, actual water added through irrigation (surface + ground)
+QIRRIG_SURF,A,8.64000E+04, mm/s, mm/day,H2O Flux, surface water irrigation
+QIRRIG_WM,A,8.64000E+04, mm/s, mm/day,H2O Flux, surface water irrigation demand sent to MOSART/WM
+QOVER,A,8.64000E+04, mm/s, mm/day,H2O Flux, surface runoff
+QRGWL,A,8.64000E+04, mm/s, mm/day,H2O Flux," surface runoff at glaciers (liquid only), wetlands, lakes"
+QRUNOFF,A,8.64000E+04, mm/s, mm/day,H2O Flux, total liquid runoff (does not include QSNWCPICE)
+QRUNOFF_NODYNLNDUSE,A,8.64000E+04, mm/s, mm/day,H2O Flux, total liquid runoff (does not include QSNWCPICE) not including correction for land use change
+QRUNOFF_R,A,8.64000E+04, mm/s, mm/day,H2O Flux, rural total runoff
+QSNOMELT,A,8.64000E+04, mm/s, mm/day,H2O Flux, snow melt
+QSNWCPICE,A,8.64000E+04, mm/s, mm/day,H2O Flux, excess snowfall due to snow capping
+QSNWCPICE_NODYNLNDUSE,A,8.64000E+04, mm/s, mm/day,H2O Flux, excess snowfall due to snow capping not including correction for land use change
+QSOIL,A,8.64000E+04, mm/s, mm/day,H2O Flux, ground evaporation
+QVEGE,A,8.64000E+04, mm/s, mm/day,H2O Flux, canopy evaporation
+QVEGT,A,8.64000E+04, mm/s, mm/day,H2O Flux, canopy transpiration
+RAIN,A,8.64000E+04, mm/s, mm/day,H2O Flux, atmospheric rain
+SNOW,A,8.64000E+04, mm/s, mm/day,H2O Flux, atmospheric snow
+SNOW_SINKS,A,8.64000E+04, mm/s, mm/day,H2O Flux, snow sinks (liquid water)
+SNOW_SOURCES,A,8.64000E+04, mm/s, mm/day,H2O Flux, snow sources (liquid water)
+SUPPLY,A,8.64000E+04, mm/s, mm/day,H2O Flux, runoff supply for land use
+ERRH2O,A,1.00000E+00, mm, mm,H2O State, total water conservation error
+ERRH2OSNO,A,1.00000E+00, mm, mm,H2O State, imbalance in snow depth (liquid water)
+GC_ICE1,A,1.00000E+00, mm, mm,H2O State, initial gridcell total ice content
+GC_LIQ1,A,1.00000E+00, mm, mm,H2O State, initial gridcell total liq content
+H2OCAN,A,1.00000E+00, mm, mm,H2O State, intercepted water
+H2OSFC,A,1.00000E+00, mm, mm,H2O State, surface water depth
+H2OSNO,A,1.00000E+00, mm, mm,H2O State, snow depth (liquid water)
+H2OSNO_TOP,A,1.00000E+00, kg/m2, mm,H2O State, mass of snow in top snow layer
+INT_SNOW,A,1.00000E+00, mm, mm,H2O State, accumulated SWE (vegetated landunits only)
+Q2M,A,1.00000E+00, kg/kg, kg/kg,H2O State, 2m specific humidity
+QBOT,A,1.00000E+00, kg/kg, kg/kg,H2O State, atmospheric specific humidity
+RH2M,A,1.00000E+00, %, kg/kg,H2O State, 2m relative humidity
+RH2M_R,A,1.00000E+00, %, kg/kg,H2O State, rural 2m relative humidity
+SNOWDP,A,1.00000E+00, m, m,H2O State, gridcell mean snow height
+SNOWICE,A,1.00000E+00, kg/m2, mm,H2O State, snow ice
+SNOWLIQ,A,1.00000E+00, kg/m2, mm,H2O State, snow liquid water
+SNOW_DEPTH,A,1.00000E+00, m, m,H2O State, snow height of snow covered area
+SOILWATER_10CM,A,1.00000E+00, kg/m2, mm,H2O State, soil liquid water + ice in top 10cm of soil (veg landunits only)
+TWS,A,1.00000E+00, mm, mm,H2O State, total water storage
+TWS_MONTH_BEGIN,A,1.00000E+00, mm, mm,H2O State, total water storage at the beginning of a month
+TWS_MONTH_END,A,1.00000E+00, mm, mm,H2O State, total water storage at the end of a month
+VOLR,A,1.00000E+00, m3, m3,H2O State, river channel total water storage
+VOLRMCH,A,1.00000E+00, m3, m3,H2O State, river channel main channel water storage
+WA,A,1.00000E+00, mm, mm,H2O State, water in the unconfined aquifer (vegetated landunits only)
+ACTUAL_IMMOB,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, actual N immobilization
+DENIT,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, total rate of denitrification
+DWT_CONV_NFLUX_GRC,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, land conversion N flux
+DWT_SLASH_NFLUX,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, slash N flux to litter and CWD due to land use
+F_DENIT,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, denitrification flux
+F_NIT,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, nitrification flux
+GROSS_NMIN,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, gross rate of N mineralization
+NDEP_TO_SMINN,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, atmospheric N deposition to soil mineral N
+NET_NMIN,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, net rate of N mineralization
+NFIX_TO_SMINN,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, symbiotic/asymbiotic N fixation to soil mineral N
+PFT_FIRE_NLOSS,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, total pft-level fire N loss
+SMINN_TO_NPOOL,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, deployment of soil mineral N uptake
+SMINN_TO_PLANT,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, plant uptake of soil mineral N
+SMIN_NO3_LEACHED,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, soil NO3 pool loss to leaching
+SMIN_NO3_RUNOFF,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, soil NO3 pool loss to runoff
+SUPPLEMENT_TO_SMINN,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, supplemental N supply
+WOOD_HARVESTN,T,3.15360E+01, gN/m^2/s, TgN/yr,N Flux, wood harvest N (to product pools)
+RETRANSN,T,1.00000E-06, gN/m^2, TgN,N State, plant pool of retranslocated N
+CWDN,T,1.00000E-06, gN/m^2, TgN,N State, coarse woody debris N
+DEADCROOTN,T,1.00000E-06, gN/m^2, TgN,N State, dead coarse root N
+DEADSTEMN,T,1.00000E-06, gN/m^2, TgN,N State, dead stem N
+DISPVEGN,T,1.00000E-06, gN/m^2, TgN,N State, displayed vegetation nitrogen
+FROOTN,T,1.00000E-06, gN/m^2, TgN,N State, fine root N
+LEAFN,T,1.00000E-06, gN/m^2, TgN,N State, leaf N
+LITR1N,T,1.00000E-06, gN/m^2, TgN,N State, LITR1 N
+LITR2N,T,1.00000E-06, gN/m^2, TgN,N State, LITR2 N
+LITR3N,T,1.00000E-06, gN/m^2, TgN,N State, LITR3 N
+LIVECROOTN,T,1.00000E-06, gN/m^2, TgN,N State, live coarse root N
+LIVESTEMN,T,1.00000E-06, gN/m^2, TgN,N State, live stem N
+NPOOL,T,1.00000E-06, gN/m^2, TgN,N State, temporary plant N pool
+SMINN,T,1.00000E-06, gN/m^2, TgN,N State, soil mineral N
+SMIN_NH4,T,1.00000E-06, gN/m^2, TgN,N State, soil mineral NH4
+SMIN_NO3,T,1.00000E-06, gN/m^2, TgN,N State, soil mineral NO3
+SOIL1N,T,1.00000E-06, gN/m^2, TgN,N State, SOIL1 N
+SOIL2N,T,1.00000E-06, gN/m^2, TgN,N State, SOIL2 N
+SOIL3N,T,1.00000E-06, gN/m^2, TgN,N State, SOIL3 N
+SOIL4N,T,1.00000E-06, gN/m^2, TgN,N State, SOIL4 N
+STORVEGN,T,1.00000E-06, gN/m^2, TgN,N State, stored vegetation nitrogen
+TOTCOLN,T,1.00000E-06, gN/m^2, TgN,N State, total column-level N (no product pools)
+TOTECOSYSN,T,1.00000E-06, gN/m^2, TgN,N State, total ecosystem N (no product pools)
+TOTLITN,T,1.00000E-06, gN/m^2, TgN,N State, total litter N
+TOTPFTN,T,1.00000E-06, gN/m^2, TgN,N State, total PFT-level N
+TOTSOMN,T,1.00000E-06, gN/m^2, TgN,N State, total soil organic matter N
+TOTVEGN,T,1.00000E-06, gN/m^2, TgN,N State, total vegetation N
+ACTUAL_IMMOB_P,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, actual P immobilization
+ADSORBTION_P,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, adsorb P flux
+BIOCHEM_PMIN,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, biochemical rate of P mineralization
+DESORPTION_P,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, desorp P flux
+DWT_CONV_PFLUX_GRC,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, land conversion P flux
+DWT_SLASH_PFLUX,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, slash P flux to litter and CWD due to land use
+GROSS_PMIN,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, gross rate of P mineralization
+NET_PMIN,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, net rate of P mineralization
+PDEP_TO_SMINP,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, atmospheric P deposition to soil mineral P
+RETRANSP,T,3.15360E+01, gP/m^2, TgP/yr,P Flux, plant pool of retranslocated P
+SMINP_LEACHED,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, soil mineral P pool loss to leaching
+SMINP_TO_PLANT,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, plant uptake of soil mineral P
+SMINP_TO_PPOOL,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, deployment of soil mineral P uptake
+SUPPLEMENT_TO_SMINP,T,3.15360E+01, gP/m^2/s, TgP/yr,P Flux, supplemental P supply
+CWDP,T,1.00000E-06, gP/m^2, TgP,P State, coarse woody debris P
+DEADCROOTP,T,1.00000E-06, gP/m^2, TgP,P State, dead coarse root P
+DEADSTEMP,T,1.00000E-06, gP/m^2, TgP,P State, dead stem P
+DISPVEGP,T,1.00000E-06, gP/m^2, TgP,P State, displayed vegetation phosphorus
+FROOTP,T,1.00000E-06, gP/m^2, TgP,P State, fine root P
+LABILEP,T,1.00000E-06, gP/m^2, TgP,P State, soil labile P
+LEAFP,T,1.00000E-06, gP/m^2, TgP,P State, leaf P
+LITR1P,T,1.00000E-06, gP/m^2, TgP,P State, LITR1 P
+LITR2P,T,1.00000E-06, gP/m^2, TgP,P State, LITR2 P
+LITR3P,T,1.00000E-06, gP/m^2, TgP,P State, LITR3 P
+LIVECROOTP,T,1.00000E-06, gP/m^2, TgP,P State, live coarse root P
+LIVESTEMP,T,1.00000E-06, gP/m^2, TgP,P State, live stem P
+OCCLP,T,1.00000E-06, gP/m^2, TgP,P State, soil occluded P
+PPOOL,T,1.00000E-06, gP/m^2, TgP,P State, temporary plant P pool
+PRIMP,T,1.00000E-06, gP/m^2, TgP,P State, soil primary P
+SECONDP,T,1.00000E-06, gP/m^2, TgP,P State, soil secondary P
+SMINP,T,1.00000E-06, gP/m^2, TgP,P State, soil mineral P
+SOIL1P,T,1.00000E-06, gP/m^2, TgP,P State, SOIL1 P
+SOIL2P,T,1.00000E-06, gP/m^2, TgP,P State, SOIL2 P
+SOIL3P,T,1.00000E-06, gP/m^2, TgP,P State, SOIL3 P
+SOIL4P,T,1.00000E-06, gP/m^2, TgP,P State, SOIL4 P
+SOLUTIONP,T,1.00000E-06, gP/m^2, TgP,P State, soil solution P
+STORVEGP,T,1.00000E-06, gP/m^2, TgP,P State, stored vegetation phosphorus
+TOTCOLP,T,1.00000E-06, gP/m^2, TgP,P State, total column-level P (no product pools)
+TOTECOSYSP,T,1.00000E-06, gP/m^2, TgP,P State, total ecosystem P (no product pools)
+TOTLITP,T,1.00000E-06, gP/m^2, TgP,P State, total litter P
+TOTLITP_1m,T,1.00000E-06, gP/m^2, TgP,P State, total litter P to 1 meter
+TOTPFTP,T,1.00000E-06, gP/m^2, TgP,P State, total PFT-level P
+TOTSOMP,T,1.00000E-06, gP/m^2, TgP,P State, total soil organic matter P
+TOTSOMP_1m,T,1.00000E-06, gP/m^2, TgP,P State, total soil organic matter P to 1 meter
+TOTVEGP,T,1.00000E-06, gP/m^2, TgP,P State, total vegetation P
+ALT,A,1.00000E+00, m, m,Physical State, current active layer thickness
+ALTMAX,A,1.00000E+00, m, m,Physical State, maximum annual active layer thickness
+FCOV,A,1.00000E+00, proportion, proportion,Physical State, fractional impermeable area
+FH2OSFC,A,1.00000E+00, proportion, proportion,Physical State, fraction of ground covered by surface water
+FINUNDATED,A,1.00000E+00, proportion, proportion,Physical State, fractional inundated area of vegetated columns
+FROST_TABLE,A,1.00000E+00, m, m,Physical State, frost table depth (vegetated landunits only)
+FSAT,A,1.00000E+00, unitless, unitless,Physical State, fractional area with water table at surface
+FSNO,A,1.00000E+00, proportion, proportion,Physical State, fraction of ground covered by snow
+FSNO_EFF,A,1.00000E+00, proportion, proportion,Physical State, effective fraction of ground covered by snow
+LAKEICETHICK,A,1.00000E+00, m, m,Physical State, thickness of lake ice (including physical expansion on freezing)
+SNOINTABS,A,1.00000E+00, %, %,Physical State, percent of incoming solar absorbed by lower snow layers
+SoilAlpha,A,1.00000E+00, unitless, unitless,Physical State, factor limiting ground evap
+TKE1,A,1.00000E+00, W/(mK), W/(mK),Physical State, top lake level eddy thermal conductivity
+ZBOT,A,1.00000E+00, m, m,Physical State, atmospheric reference height
+ZWT,A,1.00000E+00, m, m,Physical State, water table depth (vegetated landunits only)
+ZWT_CH4_UNSAT,A,1.00000E+00, m, m,Physical State, depth of water table for methane production used in non-inundated area
+ZWT_PERCH,A,1.00000E+00, m, m,Physical State, perched water table depth (vegetated landunits only)
+PBOT,A,1.00000E+00, Pa, Pa,Pressure, atmospheric pressure
+PCH4,A,1.00000E+00, Pa, Pa,Pressure, atmospheric partial pressure of CH4
+PCO2,A,1.00000E+00, Pa, Pa,Pressure, atmospheric partial pressure of CO2
+BUILDHEAT,A,1.00000E+00, W/m^2, W/m^2,Urban, heat flux from urban building interior to walls and roof
+EFLX_LH_TOT_U,A,1.00000E+00, W/m^2, W/m^2,Urban, urban total evaporation
+FGR_U,A,1.00000E+00, W/m^2, W/m^2,Urban, urban heat flux into soil/snow including snow melt
+FIRA_U,A,1.00000E+00, W/m^2, W/m^2,Urban, urban net longwave radiation
+FIRE_U,A,1.00000E+00, W/m^2, W/m^2,Urban, urban emitted longwave radiation
+FSA_U,A,1.00000E+00, W/m^2, W/m^2,Urban, urban absorbed solar radiation
+FSH_U,A,1.00000E+00, W/m^2, W/m^2,Urban, urban sensible heat
+FSM_U,A,1.00000E+00, W/m^2, W/m^2,Urban, urban snow melt heat flux
+HEAT_FROM_AC,A,1.00000E+00, W/m^2, W/m^2,Urban, sensible heat flux from air conditioning
+QRUNOFF_U,A,8.64000E+04, mm/s, mm/day,Urban, urban total runoff
+RH2M_U,A,1.00000E+00, %, %,Urban, urban 2m relative humidity
+SoilAlpha_U,A,1.00000E+00, unitless, unitless,Urban, urban factor limiting ground evap
+TBUILD,A,1.00000E+00, K, K,Urban, internal urban building temperature
+TG_U,A,1.00000E+00, K, K,Urban, urban ground temperature
+TREFMNAV_U,A,1.00000E+00, K, K,Urban, urban daily minimum of average 2-m temperature
+TREFMXAV_U,A,1.00000E+00, K, K,Urban, urban daily maximum of average 2-m temperature
+TSA_U,A,1.00000E+00, K, K,Urban, urban 2m air temperature
+URBAN_AC,A,1.00000E+00, W/m^2, W/m^2,Urban, urban air conditioning flux
+URBAN_HEAT,A,1.00000E+00, W/m^2, W/m^2,Urban, urban heating flux
+WASTEHEAT,A,1.00000E+00, W/m^2, W/m^2,Urban, sensible heat flux from heating/cooling sources of urban waste heat
+FPG_P,A,1.00000E+00, proportion, proportion,Veg State, fraction of potential gpp due to P limitation
+FPI_P,A,1.00000E+00, proportion, proportion,Veg State, fraction of potential immobilization of phosphorus
+FPG,A,1.00000E+00, proportion, proportion,Veg State, fraction of potential gpp due to N limitation
+FPI,A,1.00000E+00, proportion, proportion,Veg State, fraction of potential immobilization of nitrogen
+BTRAN,A,1.00000E+00, proportion, proportion,Veg State, transpiration beta factor
+ELAI,A,1.00000E+00, m^2/m^2, m^2/m^2,Veg State, exposed one-sided leaf area index
+ESAI,A,1.00000E+00, m^2/m^2, m^2/m^2,Veg State, exposed one-sided stem area index
+HTOP,A,1.00000E+00, m, m,Veg State, height of canopy top
+LAISHA,A,1.00000E+00,m^2/m^2,m^2/m^2,Veg State, shaded projected leaf area index
+LAISUN,A,1.00000E+00,m^2/m^2,m^2/m^2,Veg State, sunlit projected leaf area index
+TLAI,A,1.00000E+00,m^2/m^2,m^2/m^2,Veg State, total projected leaf area index
+TSAI,A,1.00000E+00,m^2/m^2,m^2/m^2,Veg State, total projected stem area index
+TAUX,A,1.00000E+00, kg/m/s^2, kg/m/s^2,Winds, zonal surface stress
+TAUY,A,1.00000E+00, kg/m/s^2, kg/m/s^2,Winds, meridional surface stress
+U10,A,1.00000E+00, m/s, m/s,Winds, 10-m wind
+U10WITHGUSTS,A,1.00000E+00, m/s, m/s,Winds, 10-m wind with gustiness enhancement included
+WIND,A,1.00000E+00, m/s, m/s,Winds, atmospheric wind velocity magnitude