diff --git a/docs/api.rst b/docs/api.rst
index e7cb2c2..9a04428 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -8,7 +8,7 @@ API
    :toctree: generated/
    :recursive:
 
-   xroms.roms_dataset
+   xroms
    derived
    interp
    roms_seawater
diff --git a/docs/calc.md b/docs/calc.md
index 2f49bfb..6632df1 100644
--- a/docs/calc.md
+++ b/docs/calc.md
@@ -394,26 +394,26 @@ ds.xroms.vertical_shear
 ds.xroms.vort
 ```
 
-### Horizontal divergence
+### Horizontal convergence
 
-Horizontal component of the currents divergence.
+Horizontal component of the currents convergence.
 
-    ds.xroms.div
+    ds.xroms.convergence
 
-    xroms.divergence(ds.u, ds.v, xgrid)
+    xroms.convergence(ds.u, ds.v, xgrid)
 
 ```{code-cell} ipython3
-ds.xroms.div
+ds.xroms.convergence
 ```
 
-### Normalized surface divergence
+### Normalized surface convergence
 
-Horizontal component of the currents divergence at the surface, normalized by $f$. This is only available through the accessor.
+Horizontal component of the currents convergence at the surface, normalized by $f$. This is only available through the accessor.
 
-    ds.xroms.div_norm
+    ds.xroms.convergence_norm
 
 ```{code-cell} ipython3
-ds.xroms.div_norm
+ds.xroms.convergence_norm
 ```
 
 ### Ertel potential vorticity
diff --git a/docs/conf.py b/docs/conf.py
index c4a179a..102afa1 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -21,6 +21,15 @@
 from importlib.metadata import version as imversion
 
 
+# to fix issue with xESMF
+# https://github.com/conda-forge/esmf-feedstock/issues/91#issuecomment-1387279692
+if "ESMFMKFILE" not in os.environ:
+    # RTD doesn't activate the env, and esmpy depends on a env var set there
+    # We assume the `os` package is in {ENV}/lib/pythonX.X/os.py
+    # See conda-forge/esmf-feedstock#91 and readthedocs/readthedocs.org#4067
+    os.environ["ESMFMKFILE"] = str(pathlib.Path(os.__file__).parent.parent / "esmf.mk")
+
+
 print("python exec:", sys.executable)
 print("sys.path:", sys.path)
 root = pathlib.Path(__file__).parent.parent.absolute()
@@ -32,7 +41,7 @@
 # -- Project information -----------------------------------------------------
 
 project = "xroms"
-copyright = "2020-2023"
+copyright = "2020-2024"
 author = "Rob Hetland, Kristen Thyng, Veronica Ruiz Xomchuk"
 
 release = imversion("xroms")
diff --git a/docs/whats_new.md b/docs/whats_new.md
index 913124b..ce78b28 100644
--- a/docs/whats_new.md
+++ b/docs/whats_new.md
@@ -1,5 +1,9 @@
 # What's New
 
+## v0.6.0 (February 9, 2023)
+* fixed error in `derived.py`'s `uv_geostrophic` function after being pointed out by @ak11283
+* updated docs so mostly well-formatted and working
+
 ## v0.5.3 (October 11, 2023)
 * change to `roms_dataset()` so that input flag `include_3D_metrics` also controls if `ds["3d"] = True`.
 
diff --git a/xroms/accessor.py b/xroms/accessor.py
index 9099dcf..f634019 100644
--- a/xroms/accessor.py
+++ b/xroms/accessor.py
@@ -75,6 +75,7 @@ def set_grid(self, xgrid):
 
         Examples
         --------
+
         >>> ds.xroms.set_grid(xgrid)
         """
         self._xgrid = xgrid
@@ -100,6 +101,7 @@ def speed(self):
 
         Examples
         --------
+
         >>> ds.xroms.speed
         """
 
@@ -120,6 +122,7 @@ def KE(self):
 
         Examples
         --------
+
         >>> ds.xroms.KE
         """
 
@@ -140,6 +143,7 @@ def ug(self):
 
         Examples
         --------
+
         >>> ds.xroms.ug
         """
 
diff --git a/xroms/derived.py b/xroms/derived.py
index 33a8fc2..a37e9cc 100644
--- a/xroms/derived.py
+++ b/xroms/derived.py
@@ -28,12 +28,14 @@ def speed(u, v, xgrid, hboundary="extend", hfill_value=None):
         for moving to rho grid.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary fill value
         selection for moving to rho grid.
@@ -122,12 +124,14 @@ def uv_geostrophic(zeta, f, xgrid, hboundary="extend", hfill_value=None, which="
         for moving f to rho grid.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         for moving f to rho grid.
@@ -135,6 +139,7 @@ def uv_geostrophic(zeta, f, xgrid, hboundary="extend", hfill_value=None, which="
         The value to use in the boundary condition with `boundary='fill'`.
     which: string, optional
         Which components of geostrophic velocity to return.
+
         * 'both': return both components of hgrad
         * 'xi': return only xi-direction.
         * 'eta': return only eta-direction.
@@ -147,10 +152,16 @@ def uv_geostrophic(zeta, f, xgrid, hboundary="extend", hfill_value=None, which="
 
     Notes
     -----
-    vg = g * zeta_eta / (d eta * f)  # on v grid
-    ug = -g * zeta_xi / (d xi * f)  # on u grid
+
+    ug = -g * zeta_eta / (d eta * f)  # on u grid
+
+    vg = g * zeta_xi / (d xi * f)  # on v grid
+
     Translation to Python of Matlab copy of surf_geostr_vel of IRD Roms_Tools.
 
+    Good resourcefor more information:
+    https://uw.pressbooks.pub/ocean285/chapter/geostrophic-balance/
+
     Examples
     --------
     >>> xroms.uv_geostrophic(ds.zeta, ds.f, xgrid)
@@ -162,11 +173,13 @@ def uv_geostrophic(zeta, f, xgrid, hboundary="extend", hfill_value=None, which="
     if which in ["both", "xi"]:
 
         # calculate derivatives of zeta
-        dzetadxi = hgrad(zeta, xgrid, which="xi")
+        dzetadeta = hgrad(zeta, xgrid, which="eta", hcoord="u")
 
         # calculate geostrophic velocities
         ug = (
-            -g * dzetadxi / to_u(f, xgrid, hboundary=hboundary, hfill_value=hfill_value)
+            -g
+            * dzetadeta
+            / to_u(f, xgrid, hboundary=hboundary, hfill_value=hfill_value)
         )
 
         ug.attrs["name"] = "ug"
@@ -175,13 +188,12 @@ def uv_geostrophic(zeta, f, xgrid, hboundary="extend", hfill_value=None, which="
         ug.name = ug.attrs["name"]
 
     if which in ["both", "eta"]:
+
         # calculate derivatives of zeta
-        dzetadeta = hgrad(zeta, xgrid, which="eta")
+        dzetadxi = hgrad(zeta, xgrid, which="xi", hcoord="v")
 
         # calculate geostrophic velocities
-        vg = (
-            g * dzetadeta / to_v(f, xgrid, hboundary=hboundary, hfill_value=hfill_value)
-        )
+        vg = g * dzetadxi / to_v(f, xgrid, hboundary=hboundary, hfill_value=hfill_value)
 
         vg.attrs["name"] = "vg"
         vg.attrs["long_name"] = "geostrophic v velocity"
@@ -214,12 +226,14 @@ def EKE(ug, vg, xgrid, hboundary="extend", hfill_value=None):
         for moving to rho grid.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         for moving to rho grid.
@@ -272,12 +286,14 @@ def dudz(u, xgrid, sboundary="extend", sfill_value=None):
         calculating z derivative.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary fill value
         associated with sboundary input.
@@ -321,12 +337,14 @@ def dvdz(v, xgrid, sboundary="extend", sfill_value=None):
         calculating z derivative.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary fill value
         associated with sboundary input.
@@ -372,12 +390,14 @@ def vertical_shear(dudz, dvdz, xgrid, hboundary="extend", hfill_value=None):
         for moving dudz and dvdz to rho grid.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         for moving to rho grid.
@@ -439,12 +459,14 @@ def relative_vorticity(
         for calculating horizontal derivatives of u and v.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -455,12 +477,14 @@ def relative_vorticity(
         for calculating horizontal derivatives of u and v.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary selection
         fill value.
@@ -537,12 +561,14 @@ def convergence(
         for calculating horizontal derivatives of u and v.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -553,12 +579,14 @@ def convergence(
         for calculating horizontal derivatives of u and v.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary selection
         fill value.
@@ -660,12 +688,14 @@ def ertel(
         horizontal grid changes too.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -678,12 +708,14 @@ def ertel(
         all vertical grid changes too.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary selection
         fill value.
diff --git a/xroms/interp.py b/xroms/interp.py
index 82834ee..d25a6b8 100644
--- a/xroms/interp.py
+++ b/xroms/interp.py
@@ -33,10 +33,12 @@ def interpll(var, lons, lats, which="pairs", **kwargs):
         Latitudes to interpolate to. Will be flattened upon input.
     which: str, optional
         Which type of interpolation to do:
+
         * "pairs": lons/lats as unstructured coordinate pairs
           (in xESMF language, LocStream).
         * "grid": 2D array of points with 1 dimension the lons and
           the other dimension the lats.
+
     **kwargs:
         passed on to xESMF Regridder class
 
@@ -56,9 +58,13 @@ def interpll(var, lons, lats, which="pairs", **kwargs):
 
     Examples
     --------
+
     To return 1D pairs of points, in this case 3 points:
+
     >>> xroms.interpll(var, [-96, -97, -96.5], [26.5, 27, 26.5], which='pairs')
+
     To return 2D pairs of points, in this case a 3x3 array of points:
+
     >>> xroms.interpll(var, [-96, -97, -96.5], [26.5, 27, 26.5], which='grid')
     """
 
@@ -178,41 +184,54 @@ def isoslice(var, iso_values, xgrid, iso_array=None, axis="Z"):
 
     Examples
     --------
+
     To calculate temperature onto fixed depths:
+
     >>> xroms.isoslice(ds.temp, np.linspace(0, -30, 50))
 
     To calculate temperature onto salinity:
+
     >>> xroms.isoslice(ds.temp, np.arange(0, 36), iso_array=ds.salt, axis='Z')
 
     Calculate lat-z slice of salinity along a constant longitude value (-91.5):
+
     >>> xroms.isoslice(ds.salt, -91.5, iso_array=ds.lon_rho, axis='X')
 
     Calculate slice of salt at 28 deg latitude
+
     >>> xroms.isoslice(ds.salt, 28, iso_array=ds.lat_rho, axis='Y')
 
     Interpolate temp to salinity values between 0 and 36 in the X direction
+
     >>> xroms.isoslice(ds.temp, np.linspace(0, 36, 50), iso_array=ds.salt, axis='X')
 
     Interpolate temp to salinity values between 0 and 36 in the Z direction
+
     >>> xroms.isoslice(ds.temp, np.linspace(0, 36, 50), iso_array=ds.salt, axis='Z')
 
     Calculate the depth of a specific isohaline (33):
+
     >>> xroms.isoslice(ds.salt, 33, iso_array=ds.z_rho, axis='Z')
 
     Calculate dye 10 meters above seabed. Either do this on the vertical
     rho grid, or first change to the w grid and then use `isoslice`. You may prefer
     to do the latter if there is a possibility that the distance above the seabed you are
     interpolating to (10 m) could be below the deepest rho grid depth.
+
     * on rho grid directly:
+
     >>> height_from_seabed = ds.z_rho + ds.h
     >>> height_from_seabed.name = 'z_rho'
     >>> xroms.isoslice(ds.dye_01, 10, iso_array=height_from_seabed, axis='Z')
+
     * on w grid:
+
     >>> var_w = ds.dye_01.xroms.to_grid(scoord='w').chunk({'s_w': -1})
     >>> ds['dye_01_w'] = var_w  # currently this is the easiest way to reattached coords xgcm variables
     >>> height_from_seabed = ds.z_w + ds.h
     >>> height_from_seabed.name = 'z_w'
     >>> xroms.isoslice(ds['dye_01_w'], 10, iso_array=height_from_seabed, axis='Z')
+
     """
 
     words = "Grid should be input."
diff --git a/xroms/roms_seawater.py b/xroms/roms_seawater.py
index b387a47..e6965cb 100644
--- a/xroms/roms_seawater.py
+++ b/xroms/roms_seawater.py
@@ -242,12 +242,14 @@ def N2(rho, xgrid, rho0=1025.0, sboundary="fill", sfill_value=np.nan):
         calculating z derivative.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary fill value
         associated with sboundary input.
@@ -306,12 +308,14 @@ def M2(
         for calculating horizontal derivatives of rho.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -322,12 +326,14 @@ def M2(
         calculating horizontal derivatives of rho.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary fill value
         associated with sboundary input.
diff --git a/xroms/tests/test_derived.py b/xroms/tests/test_derived.py
index 1ad83bb..80377f4 100644
--- a/xroms/tests/test_derived.py
+++ b/xroms/tests/test_derived.py
@@ -58,26 +58,28 @@ def test_KE():
 
 
 def test_uv_geostrophic():
-    vg = 0  # zeta only varies in xi direction
+    ug = 0  # zeta only varies in xi direction
     # test one corner of domain
-    f = ds.f[0, 0].values
+    f = ds.f[1, 1].values
     # correct dx value (eta=0)
     dzetadxi = (zeta[2] - zeta[0]) / (2 * dx)
-    ug = -g * dzetadxi / f
+    vg = g * dzetadxi / f
 
     assert np.allclose(
         xroms.uv_geostrophic(ds.zeta, ds.f, grid, which="xi")[0, 0, 0], ug
     )
-    assert np.allclose(xroms.uv_geostrophic(ds.zeta, ds.f, grid, which="eta"), vg)
+    assert np.allclose(
+        xroms.uv_geostrophic(ds.zeta, ds.f, grid, which="eta")[0, 0, 0], vg
+    )
 
 
 def test_EKE():
-    vg = 0  # zeta only varies in xi direction
+    ug = 0  # zeta only varies in xi direction
     # test one corner of domain
-    f = ds.f[0, 0].values
+    f = ds.f[1, 1].values
     # correct dx value (eta=0)
     dzetadxi = (zeta[2] - zeta[0]) / (2 * dx)
-    ug = -g * dzetadxi / f
+    vg = g * dzetadxi / f
     EKE = 0.5 * (ug**2 + vg**2)
 
     xug, xvg = xroms.uv_geostrophic(ds.zeta, ds.f, grid, which="both")
diff --git a/xroms/utilities.py b/xroms/utilities.py
index d2c4e37..ab217e8 100644
--- a/xroms/utilities.py
+++ b/xroms/utilities.py
@@ -118,9 +118,11 @@ def hgrad(
         Grid object associated with q.
     which: string, optional
         Which components of gradient to return.
+
         * 'both': return both components of hgrad.
         * 'xi': return only xi-direction.
         * 'eta': return only eta-direction.
+
     z: DataArray, ndarray, optional
         Depth [m]. If None, use z coordinate attached to q.
     hcoord: string, optional
@@ -135,12 +137,14 @@ def hgrad(
         will be used for all horizontal grid changes too.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -152,12 +156,14 @@ def hgrad(
         be used for all vertical grid changes too.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary selection
         fill value.
@@ -357,12 +363,14 @@ def ddxi(
         will be used for all horizontal grid changes too.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -374,12 +382,14 @@ def ddxi(
         be used for all vertical grid changes too.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary selection
         fill value.
@@ -463,12 +473,14 @@ def ddeta(
         will be used for grid changes too.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -480,12 +492,14 @@ def ddeta(
         be used for vertical grid changes too.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary selection
         fill value.
@@ -564,12 +578,14 @@ def ddz(
         will be used for grid changes too.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -581,12 +597,14 @@ def ddz(
         will be used for grid changes too.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary fill value
         associated with sboundary input.
@@ -743,12 +761,14 @@ def to_rho(var, xgrid, hboundary="extend", hfill_value=None):
         for grid changes.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -803,12 +823,14 @@ def to_psi(var, xgrid, hboundary="extend", hfill_value=None):
         for grid changes.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -864,12 +886,14 @@ def to_u(var, xgrid, hboundary="extend", hfill_value=None):
         for grid changes.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -924,12 +948,14 @@ def to_v(var, xgrid, hboundary="extend", hfill_value=None):
         for grid changes.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -984,12 +1010,14 @@ def to_s_rho(var, xgrid, sboundary="extend", sfill_value=None):
         for grid changes.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary selection
         fill value.
@@ -1034,12 +1062,14 @@ def to_s_w(var, xgrid, sboundary="extend", sfill_value=None):
         for grid changes.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary selection
         fill value.
@@ -1099,6 +1129,7 @@ def to_grid(
         Passed to `grid` method calls; horizontal boundary selection
         for grid changes.
         From xgcm documentation:
+
         A flag indicating how to handle boundaries:
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
@@ -1106,6 +1137,7 @@ def to_grid(
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     hfill_value: float, optional
         Passed to `grid` method calls; horizontal boundary selection
         fill value.
@@ -1116,12 +1148,14 @@ def to_grid(
         for grid changes.
         From xgcm documentation:
         A flag indicating how to handle boundaries:
+
         * None:  Do not apply any boundary conditions. Raise an error if
           boundary conditions are required for the operation.
         * 'fill':  Set values outside the array boundary to fill_value
           (i.e. a Neumann boundary condition.)
         * 'extend': Set values outside the array to the nearest array
           value. (i.e. a limited form of Dirichlet boundary condition.
+
     sfill_value: float, optional
         Passed to `grid` method calls; vertical boundary selection
         fill value.
@@ -1202,7 +1236,9 @@ def gridmean(var, xgrid, dim):
 
     Examples
     --------
+
     Note that the following two approaches are equivalent:
+
     >>> app1 = xroms.gridmean(ds.u, xgrid, ('Y','X'))
     >>> app2 = (ds.u*ds.dy_u*ds.dx_u).sum(('eta_rho','xi_u'))/(ds.dy_u*ds.dx_u).sum(('eta_rho','xi_u'))
     >>> np.allclose(app1, app2)
@@ -1251,7 +1287,9 @@ def gridsum(var, xgrid, dim):
 
     Examples
     --------
+
     Note that the following two approaches are equivalent:
+
     >>> app1 = xroms.gridsum(ds.u, xgrid, ('Z','X'))
     >>> app2 = (ds.u*ds.dz_u * ds.dx_u).sum(('s_rho','xi_u'))
     >>> np.allclose(app1, app2)
@@ -1289,6 +1327,7 @@ def xisoslice(iso_array, iso_value, projected_array, coord):
     This function has been possibly superseded by isoslice
     that wraps `xgcm.grid.transform` for the following reasons,
     but more testing is needed:
+
     * The implementation of `xgcm.grid.transform` is more robust
       than `xisoslice` which has extra code for in case iso_value
       is exactly in iso_array.
@@ -1330,16 +1369,20 @@ def xisoslice(iso_array, iso_value, projected_array, coord):
     --------
 
     Calculate lat-z slice of salinity along a constant longitude value (-91.5):
+
     >>> sl = xroms.utilities.xisoslice(ds.lon_rho, -91.5, ds.salt, 'xi_rho')
 
     Calculate a lon-lat slice at a constant z value (-10):
+
     >>> sl = xroms.utilities.xisoslice(ds.z_rho, -10, ds.temp, 's_rho')
 
     Calculate a lon-lat slice at a constant z value (-10) but without zeta changing in time:
+
     (use ds.z_rho0 which is relative to mean sea level and does not vary in time)
     >>> sl = xroms.utilities.xisoslice(ds.z_rho0, -10, ds.temp, 's_rho')
 
     Calculate the depth of a specific isohaline (33):
+
     >>> sl = xroms.utilities.xisoslice(ds.salt, 33, ds.z_rho, 's_rho')
 
     Calculate the salt 10 meters above the seabed. Either do this on the vertical
@@ -1360,19 +1403,24 @@ def xisoslice(iso_array, iso_value, projected_array, coord):
     related coordinates for plotting. For example, to accompany the lat-z slice,
     you may want the following:
 
-        # calculate z values (s_rho)
+    calculate z values (s_rho)
+
     >>> slz = xroms.utilities.xisoslice(ds.lon_rho, -91.5, ds.z_rho, 'xi_rho')
 
-        # calculate latitude values (eta_rho)
+    calculate latitude values (eta_rho)
+
     >>> sllat = xroms.utilities.xisoslice(ds.lon_rho, -91.5, ds.lat_rho, 'xi_rho')
 
-        # assign these as coords to be used in plot
+    assign these as coords to be used in plot
+
     >>> sl = sl.assign_coords(z=slz, lat=sllat)
 
-        # points that should be masked
+    points that should be masked
+
     >>> slmask = xroms.utilities.xisoslice(ds.lon_rho, -91.5, ds.mask_rho, 'xi_rho')
 
-        # drop masked values
+    drop masked values
+
     >>> sl = sl.where(slmask==1, drop=True)
     """
 
@@ -1442,11 +1490,15 @@ def subset(ds, X=None, Y=None):
         have been updated.
     X: slice, optional
         Slice in X dimension using form `X=slice(start, stop, step)`. For example,
+
         >>> X=slice(20,40,2)
+
         Indices are used for rho grid, and psi grid is reduced accordingly.
     Y: slice, optional
         Slice in Y dimension using form `Y=slice(start, stop, step)`. For example,
+
         >>> Y=slice(20,40,2)
+
         Indices are used for rho grid, and psi grid is reduced accordingly.
 
     Returns
@@ -1461,9 +1513,13 @@ def subset(ds, X=None, Y=None):
 
     Examples
     --------
+
     Subset only in Y direction:
+
     >>> xroms.subset(ds, Y=slice(50,100))
+
     Subset in X and Y:
+
     >>> xroms.subset(ds, X=slice(20,40), Y=slice(50,100))
     """
 
@@ -1506,6 +1562,7 @@ def order(var):
 
     Examples
     --------
+
     >>> xroms.order(var)
     """
 
diff --git a/xroms/xroms.py b/xroms/xroms.py
index 64ea7c3..486ec4c 100644
--- a/xroms/xroms.py
+++ b/xroms/xroms.py
@@ -813,6 +813,8 @@ def open_netcdf(
 ):
     """Return Dataset based on a single thredds or physical location.
 
+    This function is deprecated; use `xroms.open_netcdf` or `xroms.open_zarr` directly instead.
+
     Parameters
     ----------
     file: str
@@ -878,6 +880,8 @@ def open_mfnetcdf(
 ):
     """Return Dataset based on a list of netCDF files.
 
+    This function is deprecated; use `xroms.open_netcdf` or `xroms.open_zarr` directly instead.
+
     Parameters
     ----------
     files: list of strings