diff --git a/pgrid/README.md b/pgrid/README.md
new file mode 100644
index 00000000..605f7cfb
--- /dev/null
+++ b/pgrid/README.md
@@ -0,0 +1,54 @@
+# README for pgrid
+
+### This collection of programs is designed to make gridfiles for ROMS.  It works for analytical and realistic cases, and handles complex tasks like mask editing and smoothing.  It also creates other files associated with rivers, nudging to climatology, and vertical grid parameters, all in the form expected by LO/forcing and ROMS.
+
+---
+
+First edit the gridname and (if needed) a few directory locations at the top of gfun.py.  This gridname will then be used by all subsequent code.
+
+In order to keep track of several choices typically made about a grid, we use "dch," a dict of "default choices":
+```
+dch =  gfun.default_choices(Gr)
+```
+These are initialized in gfun.default_choices, but you typically override some of them in start_grid.py when doing the initial grid specification.  The choices are saved in a pickle file:
+```
+pickle.dump(dch, open(Gr['gdir'] + 'choices.p', 'wb'))
+```
+You can also go back and change things in dch later using Z_edit_dch.py.
+
+Each time you run a piece of code it makes a new grid file with the name altered to indicate what happened.  The names start as: grid_m00_r00_s00_x00.nc with the letters and numbers indicating changes to: mask, river, smoothing, or extras.
+
+You can use plot_grid.py to look at any of the grids.  You can override the default grid by calling it in ipython like:
+- run plot_grid.py -g 'sal0'
+
+Suggested order to run the code, and which dch items are used at each step:
+
+`start_grid.py` (can take awhile for big grids; be patient; e.g. 15 min for cas6)
+- analytical
+- t_dir/t_list
+- use_z_offset/z_offset
+
+`make_mask.py`
+- z_land
+- unmask_coast
+- remove_islands
+
+`edit_mask.py` (to get rid of obvious issues like Lake Washington and river channels)
+
+`carve_rivers.py`
+
+`edit_mask.py` (for real this time, perhaps running many times) NOTE: you can use an optional command line argument -d ## to change "dval" the carving depth used for lines or points from its default value of 5 m.
+
+`smooth_grid.py`
+- use_min_depth/min_depth
+
+`carve_rivers.py` (not needed?)
+
+`smooth_grid.py` (not needed?)
+
+`make_extras.py`
+- min_depth (enforced for whole grid)
+
+`grid_to_LiveOcean.py` (should be run on remote machine for large grids) (also saves S-coordinate info, so be careful what file you choose in the code)
+- nudging_edges
+- nudging_days 
diff --git a/pgrid/README.txt b/pgrid/README.txt
deleted file mode 100644
index 0555fd57..00000000
--- a/pgrid/README.txt
+++ /dev/null
@@ -1,50 +0,0 @@
-This collection of programs is designed to make gridfiles for ROMS.  It works for analytical and realistic cases, and handles complex tasks like mask editing and smoothing.  It also creates other files associated with rivers, nudging to climatology, and vertical grid parameters, all in the form expected by LiveOcean/forcing and ROMS.
-
-First edit the gridname and (if needed) a few directory locations at the top of gfun.py.  This gridname will then be used by all subsequent code.
-
-In order to keep track of several choices typically made about a grid, we use "dch," a dict of "default choices":
-- dch =  gfun.default_choices(Gr)
-These are initialized in gfun.default_choices, but you typically override some of them in start_grid.py when doing the initial grid specification.  The choices are saved in a pickle file:
-- pickle.dump(dch, open(Gr['gdir'] + 'choices.p', 'wb'))
-You can also go back and change things in dch later using Z_edit_dch.py.
-
-Each time you run a piece of code it makes a new grid file with the name altered to indicate what happened.  The names start as: grid_m00_r00_s00_x00.nc with the letters and numbers indicating changes to: mask, river, smoothing, or extras.
-
-You can use plot_grid.py to look at any of the grids.  You can override the default grid by calling it in ipython like:
-- run plot_grid.py -g 'sal0'
-
-Suggested order to run the code, and which dch items are used at each step:
-
-* start_grid.py (can take awhile for big grids; be patient; e.g. 15 min for cas6)
-    analytical
-    t_dir/t_list
-    use_z_offset/z_offset
-
-* make_mask.py
-    z_land
-    unmask_coast
-    remove_islands
-
-* edit_mask.py (to get rid of obvious issues like Lake Washington and river channels)
-
-* carve_rivers.py
-
-* edit_mask.py (for real this time, perhaps running many times)
-	NOTE: you can use an optional command line argument -d ## to change "dval"
-	the carving depth used for lines or points from its default value of 5 m.
-
-* smooth_grid.py
-    use_min_depth/min_depth
-
-* carve_rivers.py (not needed)
-
-* smooth_grid.py (not needed)
-
-* make_extras.py
-    min_depth (enforced for whole grid)
-
-* grid_to_LiveOcean.py (should be run on fjord for large grids)
-    nudging_edges
-    nudging_days
-	(also saves S-coordinate info, so be careful what file you choose in the code)
-
diff --git a/pgrid/gfun.py b/pgrid/gfun.py
index 3181e4cf..3e1eb3d3 100644
--- a/pgrid/gfun.py
+++ b/pgrid/gfun.py
@@ -1,22 +1,26 @@
-# -*- coding: utf-8 -*-
 """
 Organizational functions for pgrid.
 """
 
 # **** USER EDIT ********
-#gridname = 'aestus3'
-gridname = 'cas6'
+gridname = 'sal0'; base_gridname = 'cas6'; base_tag = 'v3'
 # **** END USER EDIT ****
 
-import os; import sys
-sys.path.append(os.path.abspath('../../LiveOcean/alpha'))
-import Lfun
-Ldir = Lfun.Lstart()
+from pathlib import Path
+from lo_tools import Lfun
 
-sys.path.append(os.path.abspath('../../LiveOcean/plotting'))
+Ldir = Lfun.Lstart(gridname=base_gridname, tag=base_tag)
 
-dir0 = Ldir['parent']
-pgdir = dir0 + 'ptools_output/pgrid/'
+def gstart():
+    """
+    This returns a dict of Path objects that tell where various things are,
+    or where they should go.
+    """
+    pgdir = Ldir['LOo'] / 'pgrid'
+    gdir = pgdir / gridname # where grid.nc will end up
+    ri_dir = Ldir['LOo'] / 'pre' / 'river' / Ldir['gtag'] / 'tracks'
+    Gr ={'gridname': gridname,'pgdir': pgdir, 'gdir': gdir,'ri_dir': ri_dir}
+    return Gr
 
 def default_choices(Gr, wet_dry=False):
     # Default choices (can override in each case)    
@@ -34,12 +38,13 @@ def default_choices(Gr, wet_dry=False):
     dch['use_z_offset'] = True
     dch['z_offset'] = -1.06
     # specify topography files to use
-    dch['t_dir'] = Gr['dir0'] + 'ptools_data/topo/'    
+    t_dir = Ldir['data'] / 'topo'
+    dch['t_dir'] = Ldir['data'] / 'topo'
     # list of topo files: coarsest to finest
-    dch['t_list'] = ['srtm15/topo15.nc',
-              'cascadia/cascadia_gridded.nc',
-             'psdem/PS_183m.nc',
-             'ttp_patch/TTP_Regional_27m_patch.nc']
+    dch['t_list'] = [t_dir / 'srtm15' / 'topo15.nc',
+              t_dir / 'cascadia' / 'cascadia_gridded.nc',
+             t_dir / 'psdem' / 'PS_183m.nc',
+             t_dir / 'ttp_patch' / 'TTP_Regional_27m_patch.nc']
  
     # MASKING
     # list of existing masks to work from
@@ -67,46 +72,16 @@ def default_choices(Gr, wet_dry=False):
         
     return dch
 
-def gstart(gridname=gridname):
-    
-    if gridname in ['aestus1', 'aestus2']:
-        ri_dir = dir0 + 'ptools_output/river/analytical/'
-    else:
-        ri_dir = dir0 + 'ptools_output/river/pnw_all_2016_07/'
-    
-    gdir = pgdir + gridname + '/'
-    Gr ={'gridname': gridname, 'dir0': dir0, 'pgdir': pgdir, 'gdir': gdir,
-         'ri_dir': ri_dir}
-    return Gr
-
-def select_file(Gr, using_old_grid=False):
+def select_file(Gr):
     # interactive selection
-    if using_old_grid==True:
-        fn_list = []
-        dir0 = Ldir['parent'] + 'LiveOcean_data/grids/'
-        gn_list = ['cascadia1', 'cascadia2']
-        for gn in gn_list:
-            fn_list.append(dir0 + gn + '/grid.nc')
-    elif using_old_grid==False:
-        print('\n** %s in <<%s>> **\n' % ('Choose file to edit', Gr['gridname']))
-        fn_list_raw = os.listdir(Gr['gdir'])
-        fn_list = []
-        for item in fn_list_raw:
-            if item[-3:] == '.nc':
-                fn_list.append(item)
-        fn_list.sort()
-    Nfn = len(fn_list)
-    fn_dict = dict(zip(range(Nfn), fn_list))
-    for nfn in range(Nfn):
-        print(str(nfn) + ': ' + fn_list[nfn])
-    my_nfn = int(input('-- Input number -- '))
-    fn = fn_dict[my_nfn]
+    fn = choose_item(Gr['gdir'], tag='.nc', itext='** Choose grid from list **')
     return fn
 
 def increment_filename(fn, tag='_m'):
+    fns = str(fn)
     # create the new file name
-    gni = fn.find(tag)
-    new_num = ('00' + str(int(fn[gni+2: gni+4]) + 1))[-2:]
-    fn_new = fn.replace(fn[gni:gni+4], tag + new_num)
-
+    gni = fns.find(tag)
+    new_num = ('00' + str(int(fns[gni+2: gni+4]) + 1))[-2:]
+    fns_new = fns.replace(fns[gni:gni+4], tag + new_num)
+    fn_new = Path(fns_new)
     return fn_new
diff --git a/pgrid/gfun_utility.py b/pgrid/gfun_utility.py
index 99361950..7ef67a4d 100755
--- a/pgrid/gfun_utility.py
+++ b/pgrid/gfun_utility.py
@@ -1,26 +1,11 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
 """
-Created on Wed Mar 29 15:55:45 2017
-
-@author: PM5
-
 Utility functions for pgrid.
 """
 import numpy as np
-import zfun
-import h5py
-import matfun
-import netCDF4 as nc
+import xarray as xr
 import seawater as sw
 
-import os
-import sys
-pth = os.path.abspath('../../LiveOcean/alpha')
-if pth not in sys.path:
-    sys.path.append(pth)
-import zrfun
-import Lfun
+from lo_tools import Lfun, zfun, zrfun
 
 def simple_grid(aa, res):
     dlat = aa[3] - aa[2]
@@ -74,46 +59,18 @@ def stretched_grid(lon_list, x_res_list, lat_list, y_res_list):
     return np.array(plon_list), np.array(plat_list)
 
 def load_bathy_nc(t_fn):
-    ds = nc.Dataset(t_fn)
-    tlon_vec = ds['lon'][:]
-    tlat_vec = ds['lat'][:]
-    tz = ds['z'][:]
+    ds = xr.open_dataset(t_fn)
+    tlon_vec = ds['lon'].values
+    tlat_vec = ds['lat'].values
+    tz = ds['z'].values
     ds.close()
     return tlon_vec, tlat_vec, tz
-
-def load_bathy(t_fn):
-    try:
-        a = h5py.File(t_fn)
-        aa = dict()
-        print(' using h5py')
-        for item in a.keys():
-            aa[item] = a[item][:]
-            #print('    ' + item + ' ' + str(aa[item].shape))
-        # reshaping because of how things are packed
-        tlon_vec = a['lon'][:].flatten()
-        tlat_vec = a['lat'][:].flatten()
-        tz = a['z'][:].T
-        #tmask = a['mask'][:].T
-        # The field "mask" was created in the topo processing
-        # as a matrix of ints with 1=ocean, 0=land.
-        # We will supersede it here with our own masking.
-        a.close()
-    except:
-        # for some reason h5py does not work with, for example,
-        # psdem/PS_183m.mat
-        tmat = matfun.loadmat(t_fn)
-        print(' using matfun')
-        # reshaping because of how things are packed
-        tlon_vec = tmat['lon'][:].flatten()
-        tlat_vec = tmat['lat'][:].flatten()
-        tz = tmat['z'][:]
-    return tlon_vec, tlat_vec, tz
     
 def load_bathy2(t_fn, lon_vec, lat_vec):
     # load a section of the new NetCDF Smith-Sandwell bathy
-    ds = nc.Dataset(t_fn)
-    Lon = ds['lon'][:]
-    Lat = ds['lat'][:]
+    ds = xr.open_dataset(t_fn)
+    tlon_vec = ds['lon'].values
+    tlat_vec = ds['lat'].values
     i0 = zfun.find_nearest_ind(Lon, lon_vec[0])
     i1 = zfun.find_nearest_ind(Lon, lon_vec[-1])
     j0 = zfun.find_nearest_ind(Lat, lat_vec[0])
@@ -125,14 +82,10 @@ def load_bathy2(t_fn, lon_vec, lat_vec):
     return tlon_vec, tlat_vec, tz
 
 def make_nc(out_fn, plon, plat, lon, lat, z, m, dch):
-    # get rid of old version
-    try:
-        os.remove(out_fn)
-    except OSError:
-        pass # assume error was because the file did not exist
     
+    out_fn.unlink(missing_ok=True)
     # create new NetCDF file
-    foo = nc.Dataset(out_fn, 'w', format='NETCDF3_CLASSIC')
+    foo = xr.Dataset(out_fn)
     
     # create dimensions
     M, L = lon.shape # use ROMS teminology
diff --git a/pgrid/start_grid.py b/pgrid/start_grid.py
index dfce778c..321b78bf 100644
--- a/pgrid/start_grid.py
+++ b/pgrid/start_grid.py
@@ -1,9 +1,4 @@
-# -*- coding: utf-8 -*-
 """
-Created on Fri Apr 29 15:34:17 2016
-
-@author: PM5
-
 Code to initialize the creation of a ROMS grid file.
 
 NOTE: the gridname is set in gfun.gstart().
@@ -12,26 +7,24 @@
 when manipulating or plotting I refer to [lon,lat]_rho as [lon,lat],
 and [lon,lat]_psi_ex as [plon,plat].
 """
+import numpy as np
+import pickle
 
-from importlib import reload
 import gfun
-reload(gfun)
-Gr =gfun.gstart()
 import gfun_utility as gfu
+from importlib import reload
+reload(gfun)
 reload(gfu)
+from lo_tools import Lfun, zfun
 
-import numpy as np
-import pickle
-
-import Lfun
-import zfun
+Gr =gfun.gstart()
 
 Lfun.make_dir(Gr['gdir'], clean=True)
 
 fn = 'grid_m00_r00_s00_x00.nc'
-out_fn = Gr['gdir'] + fn
-print(50*'*')
-print(out_fn)
+out_fn = Gr['gdir'] / fn
+print(60*'*')
+print(str(out_fn).center(60,'-'))
 
 dch = gfun.default_choices(Gr)
 
@@ -43,96 +36,19 @@
 if Gr['gridname'] == 'sal0':
     # start of a salish nest grid
     aa = [-124, -122, 47, 49]
-    res = 300 # target resolution (m)
-    plon_vec, plat_vec = gfu.simple_grid(aa, res)
-    dch['t_list'] = ['cascadia/cascadia_gridded.nc',
-             'psdem/PS_183m.nc',
-             'ttp_patch/TTP_Regional_27m_patch.nc']
-    dch['nudging_edges'] = ['north', 'west']
-
-elif Gr['gridname'] == 'hc1':
-    # mid Hood Canal nest
-    aa = [-123, -122.55, 47.5, 47.9]
-    res = 100 # target resolution (m)
+    res = 600 # target resolution (m)
     plon_vec, plat_vec = gfu.simple_grid(aa, res)
-    dch['t_list'] = ['psdem/PS_27m.nc']
+    # t_dir = Ldir['data'] / 'topo'
+    # dch['t_list'] = [
+    #           t_dir / 'cascadia' / 'cascadia_gridded.nc',
+    #          t_dir / 'psdem' / 'PS_183m.nc',
+    #          t_dir / 'ttp_patch' / 'TTP_Regional_27m_patch.nc']
     dch['nudging_edges'] = ['north', 'west']
-    dch['nudging_days'] = (0.1, 1.0)
-
-elif Gr['gridname'] == 'cas3': # a stretched MoSSea-like grid, became cas4/5
-    maxres = 1500
-    minres = 500
-    lon_list = [-127.4, -124, -122]
-    x_res_list = [maxres, minres, minres]
-    lat_list = [42, 47, 49, 50.3]
-    y_res_list = [maxres, minres, minres, maxres]
-    plon_vec, plat_vec = gfu.stretched_grid(lon_list, x_res_list,
-                                        lat_list, y_res_list)
-    dch['nudging_edges'] = ['south', 'west']
-    # new: a list of good masks to work from
-    dch['maskfiles'] = ['cas2/grid_m05_r01_s01_x01.nc', 'sal0/grid_m06_r03_s05_x02.nc']
-
-elif Gr['gridname'] == 'cas6': # an extended version of the excellent cas4/5
-    maxres = 1500
-    minres = 500
-    extres = 3000
-    lon_list = [-130, -127.4, -124, -122]
-    x_res_list = [extres, maxres, minres, minres]
-    lat_list = [42, 47, 49, 50.3, 52]
-    y_res_list = [maxres, minres, minres, maxres, extres]
-    plon_vec, plat_vec = gfu.stretched_grid(lon_list, x_res_list,
-                                        lat_list, y_res_list)
-    dch['nudging_edges'] = ['north', 'south', 'west']
-    # new: a list of good masks to work from
-    dch['maskfiles'] = ['cas5/grid_m05_r01_s02_x02.nc']
-    
-elif Gr['gridname'] == 'sj0':
-    # San Juan Islands nest, first version
-    aa = [-123.3, -122.65, 48.3, 48.8]
-    res = 100 # target resolution (m)
-    plon_vec, plat_vec = gfu.simple_grid(aa, res)
-    dch['t_list'] = ['cascadia/cascadia_gridded.nc']
-    dch['nudging_edges'] = ['north', 'south', 'east', 'west']
-    dch['nudging_days'] = (0.1, 1.0)
-
-elif Gr['gridname'] == 'aestus1': # idealized model
-    lon_list = [-1, 0, 1, 2, 3]
-    x_res_list = [5000, 1000, 1000, 5000, 5000]
-    lat_list = [44, 44.9, 45.1, 46]
-    y_res_list = [5000, 1000, 1000, 5000]
-    plon_vec, plat_vec = gfu.stretched_grid(lon_list, x_res_list,
-                                        lat_list, y_res_list)
-    dch['analytical'] = True
-    dch['nudging_edges'] = ['north', 'south', 'west']
-    
-elif Gr['gridname'] == 'aestus2':
-    # idealized model, higher resolution and larger domain than aestus1
-    lon_list = [-2, 0, 1, 2, 3]
-    x_res_list = [2500, 500, 500, 2500, 2500]
-    lat_list = [43, 44.9, 45.1, 47]
-    y_res_list = [2500, 500, 500, 2500]
-    plon_vec, plat_vec = gfu.stretched_grid(lon_list, x_res_list,
-                                        lat_list, y_res_list)
-    dch['analytical'] = True
-    dch['nudging_edges'] = ['north', 'south', 'west']
-    
-elif Gr['gridname'] == 'aestus3':
-    # idealized model, like aestus2 but cut off just inside the estuary mouth,
-    # with the goal of more manipulative forcing of the exchange flow
-    lon_list = [-2, 0, .1]
-    x_res_list = [2500, 500, 500]
-    lat_list = [43, 44.9, 45.1, 47]
-    y_res_list = [2500, 500, 500, 2500]
-    plon_vec, plat_vec = gfu.stretched_grid(lon_list, x_res_list,
-                                        lat_list, y_res_list)
-    dch['analytical'] = True
-    dch['nudging_edges'] = ['north', 'south', 'west']
 
 # save the default choices for use by other code
-pickle.dump(dch, open(Gr['gdir'] + 'choices.p', 'wb'))
+pickle.dump(dch, open(Gr['gdir'] / 'choices.p', 'wb'))
 
 plon, plat = np.meshgrid(plon_vec, plat_vec)
-#ax_lims = (plon_vec[0], plon_vec[-1], plat_vec[0], plat_vec[-1])
 
 # make box centers
 lon_vec = plon_vec[:-1] + np.diff(plon_vec)/2
@@ -157,9 +73,8 @@
     # add bathymetry automatically from files
     # m is the start of a mask: 1=water, 0=land
     m = np.ones_like(lon)
-    for t_file in dch['t_list']:
-        t_fn = dch['t_dir'] + t_file
-        print('\nOPENING BATHY FILE: ' + t_file)
+    for t_fn in dch['t_list']:
+        print('\nOPENING BATHY FILE: ' + t_fn.name)
         tlon_vec, tlat_vec, tz = gfu.load_bathy_nc(t_fn)
         if isinstance(tz, np.ma.masked_array):
             tz1 = tz.data