diff --git a/MANIFEST.in b/MANIFEST.in
index 0d5d0b6c..fc872927 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1 +1 @@
-recursive-include docs example scripts
+recursive-include docs example
diff --git a/docs/api.rst b/docs/api.rst
index fcdd6e68..a44a8435 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -35,9 +35,17 @@ xfaster_exec
   submits the job to a grid system
 
 spec_tools
--------------
+----------
 .. automodule:: xfaster.spec_tools
   :members:
   :undoc-members:
 
-  More generic functions used by methods in xfaster_class
+  More generic functions used by methods in :py:mod:`~xfaster.xfaster_class`
+
+gcorr_tools
+-----------
+.. automodule:: xfaster.gcorr_tools
+  :members:
+  :undoc-members:
+
+  Tools for iteratively computing corrections to the ``g_ell`` factor
diff --git a/docs/gcorr.rst b/docs/gcorr.rst
new file mode 100644
index 00000000..23638bb5
--- /dev/null
+++ b/docs/gcorr.rst
@@ -0,0 +1,79 @@
+G-Correction Calculation
+========================
+
+Code package for calculating correction factors to g_ell using pre-generated
+simulation maps for each single mask. This procedure assumes you are running
+this code on a cluster-- it is extremely slow to do otherwise.
+
+There is one main command-line utility that calls two functions from the
+:py:mod:`~xfaster.gcorr_tools` library:
+
+1. :py:func:`~xfaster.gcorr_tools.xfaster_gcorr_once` -- This function calls
+   XFaster to run or submit jobs for gcorr runs.
+2. :py:func:`~xfaster.gcorr_tools.process_gcorr` -- A function that computes the
+   gcorr correction from the ensemble of bandpowers, updates the running total,
+   and backs up the necessary files from each iteration.
+3. ``xfaster gcorr`` -- An xfaster command-line utility, that calls function 1
+   and 2 to get a new gcorr.npz file each time. This is the main code you'll
+   run.
+
+There is also a config file with options specific to computing gcorr.
+
+Procedure
+---------
+
+1. Edit the gcorr config file to suit your purposes. Examples are given for
+   signal and null runs in the examples directory. Required fields are:
+
+   * ``null`` -- must be true for null tests and false for signal runs
+   * ``map_tags`` -- comma-separated list of map tags
+   * ``data_subset`` -- the globabble data_subset argument to xfaster_run, but
+     without map tags. So, "full", "chunk*", etc.
+   * ``output_root`` -- the parent directory where your gcorr XFaster runs will
+     be written
+   * ``nsim`` -- the number of simulations to use to compute gcorr
+   * ``[xfaster_opts]`` -- this is where you'll put any options that will be
+     directly input to :py:func:`~xfaster.xfaster_exec.xfaster_run`
+   * ``[submit_opts]`` -- this is where you'll put any options that will be
+     directly input to :py:func:`~xfaster.xfaster_exec.xfaster_submit`, in
+     addition to those in ``[xfaster_opts]``
+
+2. Run ``xfaster gcorr`` once to get the full set of XFaster output files in the
+   output directory.  Since we haven't computed gcorr yet, this will set
+   ``apply_gcorr=False``. Make sure to use as many OMP threads as possible since
+   this is the step where the sims_xcorr file, which benefits the most from
+   extra threads, is computed.  Your command should look like this::
+
+      xfaster gcorr path-to-my-gcorr-config.ini
+
+3. Run ``xfaster gcorr`` until convergence is reached. In practice, you will run
+   the command above and wait for it to finish. If you include the
+   ``--max-iters`` option with a non-zero value, the code will try to determine
+   whether convergence or max_iters has been reached and stop on its own.
+   Otherwise, you can look at the correction-to-the-correction that it both
+   prints and plots (it should converge to 1s for TT, EE, BB), and rerun the
+   same command if it hasn't converged.  In much more detail, here's what the
+   code does:
+
+   1. Call :py:func:`~xfaster.gcorr_tools.xfaster_gcorr_once` for the 0th sim
+      seed while also reloading gcorr (if this is not the first iteration).
+      This does a couple things-- saves the new gcorr in the ``masks_xcorr``
+      file, so later seeds will use the right thing. And recompute the transfer
+      function, which doesn't depend on the ``sim_index``, so is only necessary
+      to do once.
+   2. After the transfer functions are all on disk, submit individual jobs for
+      all the other seeds, just doing the bandpowers step for those.
+   3. Once they're all done, run :py:func:`~xfaster.gcorr_tools.compute_gcal`,
+      and save a correction-to-gcorr as ``gcorr_corr_<tag>_iter<iter>.npz`` in the
+      rundir.
+   4. If not the first iteration, load up the correction-to-gcorr computed for
+      this iteration. Multiply it by the total gcorr, and save that to the
+      output directory as ``gcorr_total_<tag>_iter<iter>.npz``.
+   5. Plot gcorr total and the correction to gcorr total. Save in rundir/plots.
+   6. Clear out rundir bandpowers/transfer functions/logs.
+   7. Exit.
+
+4. After convergence is reached, copy the gcorr_total file for the last
+   iteration from the rundir to the mask directory, labeling it
+   ``mask_map_<tag>_gcorr.npz`` for signal or ``mask_map_<tag>_gcorr_null.npz``
+   for null.
diff --git a/docs/index.rst b/docs/index.rst
index 04a96292..3b10faf8 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -20,6 +20,7 @@ More details of the math can be found in :ref:`Algorithm<Algorithm>`, and a deta
   quickstart
   algorithm
   notebooks/XFaster_Tutorial
+  gcorr
   api
 
 Indices and tables
diff --git a/docs/notebooks/XFaster_Tutorial.ipynb b/docs/notebooks/XFaster_Tutorial.ipynb
index 614972e7..5b935af2 100644
--- a/docs/notebooks/XFaster_Tutorial.ipynb
+++ b/docs/notebooks/XFaster_Tutorial.ipynb
@@ -420,7 +420,7 @@
    "source": [
     "We have now finished all the set-up. Now it's time to start calculating things. First, we compute the cross spectra of our masks--which will be needed for computing the $K_{\\ell \\ell^{\\prime}}$ mode-coupling matrix--and the $g_\\ell$ mode-counting factor. This is done with [get_mask_weights()](../api.rst#xfaster.xfaster_class.XFaster.get_mask_weights).\n",
     "\n",
-    "For the example, we will not apply an empirical correction to $g_\\ell$, the calibration of which is discussed in Section 2.3.2 of the [XFaster paper](https://arxiv.org/abs/2104.01172).  Otherwise, we would set `apply_gcorr=True`, and the code would look in the masks directory for this correction file for each map. The `reload_gcorr` option is only useful when doing the empirical calibration; it reloads the file, while by-passing the checkpoint tree that is usually performed after `get_mask_weights`."
+    "For the example, we will not apply an empirical correction to $g_\\ell$, the calibration of which is discussed in Section 2.3.2 of the [XFaster paper](https://arxiv.org/abs/2104.01172).  Otherwise, we would set `apply_gcorr=True`, and the code would look in the masks directory for this correction file for each map. The `reload_gcorr` option is only useful when doing the empirical calibration; it reloads the file, while by-passing the checkpoint tree that is usually performed after `get_mask_weights`.  Instructions are [provided](../gcorr.rst) for constructing this calibration using the example dataset."
    ]
   },
   {
diff --git a/scripts/gcorr/gcorr_config_null.ini b/example/gcorr_config_null.ini
similarity index 77%
rename from scripts/gcorr/gcorr_config_null.ini
rename to example/gcorr_config_null.ini
index a1394296..af9d6c6d 100644
--- a/scripts/gcorr/gcorr_config_null.ini
+++ b/example/gcorr_config_null.ini
@@ -3,16 +3,16 @@
 null = true
 map_tags = 95,150
 data_subset = full
-output_root = ../../example/null/gcorr_run
+output_root = null/gcorr_run
 nsim = 100
 
 # Options we can directly pass to XFaster
 [xfaster_opts]
-config = ../../example/config_example.ini
+config = config_example.ini
 # Make null map set using the same set as for the signal run.
-data_root = ../../example/maps_example_half1
+data_root = maps_example_half1
 # Need second data root for null run
-data_root2 = ../../example/maps_example_half2
+data_root2 = maps_example_half2
 signal_type = synfast
 # Noise type needed for nulls
 noise_type = gaussian
diff --git a/scripts/gcorr/gcorr_config_signal.ini b/example/gcorr_config_signal.ini
similarity index 82%
rename from scripts/gcorr/gcorr_config_signal.ini
rename to example/gcorr_config_signal.ini
index 73025f64..7e69baad 100644
--- a/scripts/gcorr/gcorr_config_signal.ini
+++ b/example/gcorr_config_signal.ini
@@ -3,13 +3,13 @@
 null = false
 map_tags = 95,150
 data_subset = full
-output_root = ../../example/gcorr_run
+output_root = gcorr_run
 nsim = 100
 
 # Options we can directly pass to XFaster
 [xfaster_opts]
-config = ../../example/config_example.ini
-data_root = ../../example/maps_example
+config = config_example.ini
+data_root = maps_example
 # data_root2 = # change for nulls
 signal_type = synfast
 # noise type ignored for signal gcorr
diff --git a/scripts/gcorr/README b/scripts/gcorr/README
deleted file mode 100644
index 90077d28..00000000
--- a/scripts/gcorr/README
+++ /dev/null
@@ -1,72 +0,0 @@
-Code package for calculating correction factors to g_ell using pre-generated
-simulation maps for each single mask. This procedure assumes you are running
-this code on a cluster-- it is extremely slow to do otherwise.
-
-There is one main script that calls two functions from the xfaster.gcorr_tools
-library:
-
-1. xfaster_gcorr ---- This function calls XFaster to run or submit jobs for
-   gcorr runs.
-2. process_gcorr ---- A function that computes the gcorr correction from the
-   ensemble of bandpowers, updates the running total, and backs up the necessary
-   files from each iteration.
-3. iterate_gcorr.py ---- A iteration script, that calls function 1 and 2 to get
-   a new gcorr.npz file each time. This is the main script you'll run.
-
-There is also a config file with options specific to computing gcorr.
-
----------------------
-Gcorr run procedure:
-
-1. Edit the gcorr config file to suit your purposes. An example is given.
-   Required fields are:
-   * null -- must be true for null tests and false for signal runs
-   * map_tags -- comma-separated list of map tags
-   * data_subset -- the globabble data_subset argument to xfaster_run,
-       but without map tags. So, "full", "chunk*", etc.
-   * output_root -- the parent directory where your gcorr XFaster runs will
-       be written
-   * nsim -- the number of simulations to use to compute gcorr
-   * [xfaster_opts] -- this is where you'll put any options
-       that will be directly input to xfaster_run
-   * [submit_opts] -- this is where you'll put any options
-       that will be directly input to xfaster_submit,
-       in addition to those in [xfaster_opts]
-
-2. Run iterate_gcorr.py once to get the full set of XFaster output files in the
-   output directory.  Since we haven't computed gcorr yet, this will set
-   apply_gcorr=False. Make sure to use as many OMP threads as possible since
-   this is the step where the sims_xcorr file, which benefits the most from
-   extra threads, is computed.  Your command should look like this:
-
-      python iterate_gcorr.py path-to-my-gcorr-config.ini
-
-3. Run iterate_gcorr.py until convergence is reached. In practice, you will run
-   the command above and wait for it to finish. If you include the `--max-iters`
-   option with a non-zero value, the code will try to determine whether
-   convergence or max_iters has been reached and stop on its own.  Otherwise,
-   you can look at the correction-to-the-correction that it both prints and
-   plots (it should converge to 1s for TT, EE, BB), and if it hasn't converged,
-   up+enter (redo) the same command you just did.  In much more detail, here's
-   what the code does:
-
-   1. Call xfaster_gcorr for the 0th sim seed while also reloading gcorr (if
-      this is not the first iteration).  This does a couple things-- saves the
-      new gcorr in the masks_xcorr file, so later seeds will use the right
-      thing. And recompute the transfer function, which doesn't depend on the
-      sim_index, so is only necessary to do once.
-   2. After the transfer functions are all on disk, submit individual jobs for
-      all the other seeds, just doing the bandpowers step for those.
-   3. Once they're all done, run compute_gcal, and save a correction-to-gcorr
-      as gcorr_corr_<tag>_iter<iter>.npz in the rundir.
-   4. If not the first iteration, load up the correction-to-gcorr computed for
-      this iteration. Multiply it by the total gcorr, and save that to the
-      output directory as gcorr_total_<tag>_iter<iter>.npz.
-   5. Plot gcorr total and the correction to gcorr total. Save in rundir/plots.
-   6. Clear out rundir bandpowers/transfer functions/logs.
-   7. Exit.
-
-4. After convergence is reached, copy the gcorr_total file for the last
-   iteration from the rundir to the mask directory, labeling it
-   mask_map_<tag>_gcorr.npz for signal or mask_map_<tag>_gcorr_null.npz for
-   null.
diff --git a/scripts/gcorr/iterate_gcorr.py b/scripts/gcorr/iterate_gcorr.py
deleted file mode 100644
index 6034cf3a..00000000
--- a/scripts/gcorr/iterate_gcorr.py
+++ /dev/null
@@ -1,176 +0,0 @@
-"""
-A iteration script used to update g_corr factors. Can either be run in series
-or submit jobs to run in parallel.
-"""
-import os
-import argparse as ap
-import subprocess as sp
-from xfaster import gcorr_tools as gt
-from xfaster.batch_tools import batch_sub
-from matplotlib import use
-
-use("agg")
-
-P = ap.ArgumentParser()
-P.add_argument("config", help="The config file for gcorr computation")
-P.add_argument(
-    "-s",
-    "--submit",
-    action="store_true",
-    help="Submit jobs, instead of running serially on current session",
-)
-P.add_argument(
-    "-t",
-    "--map-tags",
-    nargs="+",
-    help="Subset of map tags to iterate over",
-)
-P.add_argument(
-    "-i",
-    "--iternums",
-    type=int,
-    nargs="+",
-    help="Iteration number to start with for each tag in --map-tags.  If one "
-    "number is given, use the same index for all tags.  All files for iterations "
-    "above this number will be removed.  If not supplied, iterations will increment "
-    "to the next index.",
-)
-P.add_argument(
-    "--keep-iters",
-    action="store_true",
-    help="Store outputs from each iteration in a separate directory",
-)
-P.add_argument(
-    "--max-iters",
-    default=0,
-    type=int,
-    help="Maximum number of iterations to run.  If 0, run once and exit.",
-)
-P.add_argument(
-    "-c",
-    "--converge-criteria",
-    type=float,
-    default=0.01,
-    help="Maximum fractional change in gcorr that indicates "
-    "convergence and stops iteration",
-)
-P.add_argument(
-    "--allow-extreme",
-    action="store_true",
-    help="Do not clip gcorr corrections that are too large at each "
-    "iteration.  Try this if iterations are not converging.",
-)
-P.add_argument(
-    "--gcorr-fit-hist",
-    action="store_true",
-    help="Fit bandpower histogram to a lognorm distribution to compute gcorr",
-)
-P.add_argument(
-    "-a",
-    "--analyze-only",
-    action="store_true",
-    help="Compute and store gcorr files for the current iteration.  Typically used "
-    "internally by the script, and should not be called by the user.",
-)
-
-args = P.parse_args()
-
-# load configuration
-g_cfg = gt.get_gcorr_config(args.config)
-
-tags = g_cfg["gcorr_opts"]["map_tags"]
-if args.map_tags:
-    tags = [t for t in args.map_tags if t in tags]
-
-iternums = {t: None for t in tags}
-if args.iternums:
-    if len(args.iternums) == 1:
-        args.iternums = args.iternums * len(tags)
-    for t, i in zip(tags, args.iternums):
-        iternums[t] = i
-
-null = g_cfg["gcorr_opts"].get("null", False)
-nsim = g_cfg["gcorr_opts"]["nsim"]
-rundir = g_cfg["gcorr_opts"]["output_root"]
-xf_submit_opts = g_cfg.get("submit_opts", {})
-submit_opts = xf_submit_opts.copy()
-submit_opts.pop("num_jobs")
-
-# sim ensemble options
-run_opts = dict(
-    data_subset=g_cfg["gcorr_opts"]["data_subset"],
-    data_subset2=g_cfg["gcorr_opts"].get("data_subset2", None),
-    output_root=rundir,
-    null=null,
-    num_sims=nsim,
-    **g_cfg["xfaster_opts"],
-)
-if args.submit:
-    run_opts.update(submit=True, **xf_submit_opts)
-
-# gcorr analysis options
-gcorr_opts = dict(
-    output_root=rundir,
-    null=null,
-    num_sims=nsim,
-    gcorr_fit_hist=args.gcorr_fit_hist,
-    allow_extreme=args.allow_extreme,
-    keep_iters=args.keep_iters,
-    converge_criteria=args.converge_criteria,
-    max_iters=args.max_iters,
-)
-
-if args.analyze_only:
-    assert len(tags) == 1, "Analyze one tag at a time"
-    if gt.process_gcorr(output_tag=tags[0], iternum=iternums[tags[0]], **gcorr_opts):
-        raise RuntimeError("Stopping iterations")
-    raise SystemExit
-
-# build command for this script
-cmd = ["python", os.path.abspath(__file__), os.path.abspath(args.config)]
-for k in [
-    "allow_extreme",
-    "gcorr_fit_hist",
-    "keep_iters",
-]:
-    if getattr(args, k):
-        cmd += ["--{}".format(k.replace("_", "-"))]
-for k in ["converge_criteria", "max_iters"]:
-    cmd += ["--{}".format(k.replace("_", "-")), str(getattr(args, k))]
-
-if args.max_iters == 0:
-    args.max_iters = 1
-
-# run
-for tag in tags:
-    # setup for next iteration
-    iternum0 = gt.get_next_iter(
-        output_root=rundir, output_tag=tag, iternum=iternums[tag]
-    )
-    if iternum0 > args.max_iters:
-        raise ValueError(
-            "Tag {} iteration {} > max {}".format(tag, iternum0, args.max_iters)
-        )
-    gcorr_job = None
-
-    for iternum in range(iternum0, args.max_iters):
-        print("Starting {} iteration {}".format(tag, iternum))
-
-        # compute ensemble bandpowers
-        if args.submit:
-            run_opts["dep_afterok"] = gcorr_job
-        bp_jobs = gt.xfaster_gcorr(output_tag=tag, iternum=iternum, **run_opts)
-
-        # compute gcorr
-        if args.submit:
-            # submit analysis job
-            gcorr_job = batch_sub(
-                cmd + ["-a", "-t", tag, "-i", str(iternum)],
-                name="gcorr_{}".format(tag),
-                workdir=os.path.abspath(os.path.join(rundir, tag, "logs")),
-                dep_afterok=bp_jobs,
-                **submit_opts,
-            )
-        else:
-            if gt.process_gcorr(output_tag=tag, iternum=iternum, **gcorr_opts):
-                break
diff --git a/xfaster/gcorr_tools.py b/xfaster/gcorr_tools.py
index 6427d9b4..3aabf0cf 100644
--- a/xfaster/gcorr_tools.py
+++ b/xfaster/gcorr_tools.py
@@ -105,7 +105,7 @@ def get_next_iter(output_root="xfaster_gcal", output_tag=None, iternum=None):
     return iternum
 
 
-def xfaster_gcorr(
+def xfaster_gcorr_once(
     output_root="xfaster_gcal",
     output_tag=None,
     data_subset="full",
diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py
index e639e27a..8a254e76 100644
--- a/xfaster/xfaster_exec.py
+++ b/xfaster/xfaster_exec.py
@@ -18,6 +18,7 @@
     "xfaster_submit",
     "xfaster_dump",
     "xfaster_diff",
+    "xfaster_gcorr",
     "xfaster_main",
     "XFasterJobGroup",
 ]
@@ -1073,6 +1074,7 @@ def add_arg(
         ("submit", "submit xfaster job"),
         ("dump", "dump archive from xfaster job to stdout"),
         ("diff", "compare two archive files"),
+        ("gcorr", "iterate over gcorr computation"),
     ]:
         PP = S.add_parser(mode, help=helptext, **parser_opts)
 
@@ -1126,6 +1128,79 @@ def add_arg(
             )
             continue
 
+        if mode == "gcorr":
+            add_arg(
+                PP, "config", positional=True, help="Config file for gcorr computation"
+            )
+            add_arg(
+                PP,
+                "map_tags",
+                short="-t",
+                nargs="+",
+                metavar="TAG",
+                help="Subset of map tags to iterate over",
+            )
+            add_arg(
+                PP,
+                "iternums",
+                short="-i",
+                nargs="+",
+                metavar="N",
+                argtype=int,
+                help="Iteration number to start with for each tag in --map-tags",
+            )
+            add_arg(
+                PP,
+                "keep_iters",
+                default=False,
+                action="store_true",
+                help="Store outputs from each iteration in a separate directory",
+            )
+            add_arg(
+                PP,
+                "max_iters",
+                argtype=int,
+                help="Maximum number of iterations to run (run once and exit if 0)",
+            )
+            add_arg(
+                PP,
+                "converge_criteria",
+                short="-c",
+                argtype=float,
+                help="Maximum fractional change in gcorr that indicates convergence and stops iteration",
+            )
+            add_arg(
+                PP,
+                "allow_extreme",
+                default=False,
+                action="store_true",
+                help="Do not clip large gcorr corrections at each iteration (useful if not converging)",
+            )
+            add_arg(
+                PP,
+                "gcorr_fit_hist",
+                default=False,
+                action="store_true",
+                help="Fit bandpower histogram to a lognorm distribution to compute gcorr",
+            )
+            add_arg(
+                PP,
+                "analyze_only",
+                short="-a",
+                default=False,
+                action="store_true",
+                help="Compute and store gcorr files for the current iteration (internal use only)",
+            )
+            add_arg(
+                PP,
+                "submit",
+                short="-s",
+                default=False,
+                action="store_true",
+                help="Submit jobs with dependencies rather than running serially",
+            )
+            continue
+
         # common options
         G = PP.add_argument_group("common options")
         add_arg(G, "config", required=True)
@@ -1399,7 +1474,7 @@ def add_arg(
             v = v[0]
 
     # test mode
-    if args.mode not in ["submit", "dump", "diff"]:
+    if args.mode not in ["submit", "dump", "diff", "gcorr"]:
         if args.test:
             msg = ",\n".join(
                 "{}={!r}".format(k, v) for k, v in sorted(vars(args).items())
@@ -1872,6 +1947,162 @@ def compare(d1, d2, prefix=""):
     compare(data1, data2)
 
 
+def xfaster_gcorr(
+    config,
+    map_tags=None,
+    iternums=None,
+    keep_iters=False,
+    max_iters=0,
+    converge_criteria=0.01,
+    allow_extreme=False,
+    gcorr_fit_hist=False,
+    analyze_only=False,
+    submit=False,
+):
+    """
+    Iteration driver for computing gcorr factors.
+
+    Arguments
+    ---------
+    config : str
+        The config file for gcorr computation.
+    map_tags : list of str
+        Subset of map tags to iterate over.
+    iternums : list of int
+        Iteration number to start with for each tag in `map_tags`.  If one
+        number is given, use the same index for all tags.  All files for
+        iterations above this number will be removed.  If not supplied,
+        iterations will increment to the next index.
+    keep_iters : bool
+        If True, store outputs from each iteration in a separate directory.
+    max_iters : int
+        Maximum number of iterations to run.  If 0, run once and exit.
+    converge_criteria : float
+        Maximum fractional change in gcorr that indicates convergence and stops
+        iteration.
+    allow_extreme : bool
+        If True, do not clip gcorr corrections that are too large at each
+        iteration.  Try this if iterations are not converging.
+    gcorr_fit_hist : bool
+        If True, fit bandpower histogram to a lognorm distribution to compute gcorr.
+    analyze_only : bool
+        If True, compute and store gcorr files for the current iteration.
+        Typically used internally by the script, and should not be called by the
+        user.
+    submit : bool
+        If True, submit a sequence of jobs with dependencies.  Otherwise, runs
+        serially on the hots.
+    """
+    from . import gcorr_tools as gt
+    from matplotlib import use
+
+    use("agg")
+
+    # load configuration
+    g_cfg = gt.get_gcorr_config(config)
+
+    tags = g_cfg["gcorr_opts"]["map_tags"]
+    if map_tags:
+        tags = [t for t in map_tags if t in tags]
+
+    if iternums:
+        if len(iternums) == 1:
+            iternums = iternums * len(tags)
+        iternums = dict(zip(tags, iternums))
+    else:
+        iternums = {t: None for t in tags}
+    for t in tags:
+        iternums.setdefault(t, None)
+
+    null = g_cfg["gcorr_opts"].get("null", False)
+    nsim = g_cfg["gcorr_opts"]["nsim"]
+    rundir = g_cfg["gcorr_opts"]["output_root"]
+    xf_submit_opts = g_cfg.get("submit_opts", {})
+    submit_opts = xf_submit_opts.copy()
+    submit_opts.pop("num_jobs")
+
+    # sim ensemble options
+    run_opts = dict(
+        data_subset=g_cfg["gcorr_opts"]["data_subset"],
+        data_subset2=g_cfg["gcorr_opts"].get("data_subset2", None),
+        output_root=rundir,
+        null=null,
+        num_sims=nsim,
+        **g_cfg["xfaster_opts"],
+    )
+    if submit:
+        run_opts.update(submit=True, **xf_submit_opts)
+
+    # gcorr analysis options
+    gcorr_opts = dict(
+        output_root=rundir,
+        null=null,
+        num_sims=nsim,
+        gcorr_fit_hist=gcorr_fit_hist,
+        allow_extreme=allow_extreme,
+        keep_iters=keep_iters,
+        converge_criteria=converge_criteria,
+        max_iters=max_iters,
+    )
+
+    if analyze_only:
+        assert len(tags) == 1, "Analyze one tag at a time"
+        if gt.process_gcorr(
+            output_tag=tags[0], iternum=iternums[tags[0]], **gcorr_opts
+        ):
+            raise RuntimeError("Stopping iterations")
+        raise SystemExit
+
+    # build command for this script
+    cmd = ["xfaster", "gcorr", os.path.abspath(config)]
+    for k in [
+        "allow_extreme",
+        "gcorr_fit_hist",
+        "keep_iters",
+    ]:
+        if locals()[k]:
+            cmd += ["--{}".format(k.replace("_", "-"))]
+    for k in ["converge_criteria", "max_iters"]:
+        cmd += ["--{}".format(k.replace("_", "-")), str(locals()[k])]
+
+    if max_iters == 0:
+        max_iters = 1
+
+    # run
+    for tag in tags:
+        # setup for next iteration
+        iternum0 = gt.get_next_iter(
+            output_root=rundir, output_tag=tag, iternum=iternums[tag]
+        )
+        if iternum0 > max_iters:
+            raise ValueError(
+                "Tag {} iteration {} > max {}".format(tag, iternum0, max_iters)
+            )
+        gcorr_job = None
+
+        for iternum in range(iternum0, max_iters):
+            print("Starting {} iteration {}".format(tag, iternum))
+
+            # compute ensemble bandpowers
+            if submit:
+                run_opts["dep_afterok"] = gcorr_job
+            bp_jobs = gt.xfaster_gcorr_once(output_tag=tag, iternum=iternum, **run_opts)
+
+            # compute gcorr
+            if submit:
+                # submit analysis job
+                gcorr_job = bt.batch_sub(
+                    cmd + ["-a", "-t", tag, "-i", str(iternum)],
+                    name="gcorr_{}".format(tag),
+                    workdir=os.path.abspath(os.path.join(rundir, tag, "logs")),
+                    dep_afterok=bp_jobs,
+                    **submit_opts,
+                )
+            else:
+                if gt.process_gcorr(output_tag=tag, iternum=iternum, **gcorr_opts):
+                    break
+
+
 def xfaster_main():
     """
     Main entry point for command-line interface.
@@ -1895,3 +2126,7 @@ def xfaster_main():
     elif mode == "diff":
         # diff archive files
         xfaster_diff(**args)
+
+    elif mode == "gcorr":
+        # run gcorr iteration
+        xfaster_gcorr(**args)