diff --git a/src/wrf/interp.py b/src/wrf/interp.py
index 12c4366..23fb4ca 100755
--- a/src/wrf/interp.py
+++ b/src/wrf/interp.py
@@ -7,7 +7,8 @@ from .extension import (_interpz3d, _vertcross, _interpline, _smooth2d,
                         _monotonic, _vintrp, _interpz3d_lev2d)
 
 from .metadecorators import set_interp_metadata
-from .util import extract_vars, is_staggered, get_id, to_np, get_iterable
+from .util import (extract_vars, is_staggered, get_id, to_np, get_iterable,
+                   is_moving_domain, is_latlon_pair)
 from .py3compat import py3range
 from .interputils import get_xy, get_xy_z_params, to_xy_coords
 from .constants import Constants, default_fill, ConversionFactors
@@ -257,10 +258,6 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
         :class:`numpy.ndarray` object with no metadata.
         
     """
-    if timeidx is None:
-        raise ValueError("'timeidx' must be a positive or negative integer")
-    
-    
     # Some fields like uvmet have an extra left dimension for the product
     # type, we'll handle that iteration here.
     multi = True if field3d.ndim - vert.ndim == 1 else False
@@ -275,9 +272,45 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
         end_point_xy = None
         pivot_point_xy = None
         
+        if (latlon is True or is_latlon_pair(start_point) or 
+            is_latlon_pair(pivot_point)):
+            
+            if wrfin is not None:
+                is_moving = is_moving_domain(wrfin)
+            else:
+                is_moving = False
+        
+            if timeidx is None:
+                if wrfin is not None:
+                    # Moving nests aren't supported with ALL_TIMES because the 
+                    # domain could move outside of the line, which causes 
+                    # crashes or different line lengths.
+                    if is_moving:
+                        raise ValueError("Requesting all times with a moving nest "
+                                         "is not supported when using lat/lon "
+                                         "cross sections because the domain could "
+                                         "move outside of the cross section. "
+                                         "You must request each time "
+                                         "individually.")
+                    else:
+                        # Domain not moving, just use 0
+                        _timeidx = 0
+                
+                # If using grid coordinates, then don't care about lat/lon
+                # coordinates. Just use 0.
+                else:
+                    _timeidx = 0
+            else:
+                if is_moving:
+                    _timeidx = timeidx
+                else:
+                    # When using non-moving nests, set the time to 0 
+                    # to avoid problems downstream
+                    _timeidx = 0
+        
         if pivot_point is not None:
             if pivot_point.lat is not None and pivot_point.lon is not None:
-                xy_coords = to_xy_coords(pivot_point, wrfin, timeidx, 
+                xy_coords = to_xy_coords(pivot_point, wrfin, _timeidx, 
                                          stagger, projection, ll_point)
                 pivot_point_xy = (xy_coords.x, xy_coords.y)
             else:
@@ -285,14 +318,14 @@ def vertcross(field3d, vert, levels=None, missing=default_fill(np.float64),
                 
         if start_point is not None and end_point is not None:
             if start_point.lat is not None and start_point.lon is not None:
-                xy_coords = to_xy_coords(start_point, wrfin, timeidx, 
+                xy_coords = to_xy_coords(start_point, wrfin, _timeidx, 
                                          stagger, projection, ll_point)
                 start_point_xy = (xy_coords.x, xy_coords.y)
             else:
                 start_point_xy = (start_point.x, start_point.y)
                 
             if end_point.lat is not None and end_point.lon is not None:
-                xy_coords = to_xy_coords(end_point, wrfin, timeidx, 
+                xy_coords = to_xy_coords(end_point, wrfin, _timeidx, 
                                          stagger, projection, ll_point)
                 end_point_xy = (xy_coords.x, xy_coords.y)
             else:
@@ -442,8 +475,6 @@ def interpline(field2d, pivot_point=None,
         
         
     """
-    if timeidx is None:
-        raise ValueError("'timeidx' must be a positive or negative integer")
     
     try:
         xy = cache["xy"]
@@ -452,9 +483,45 @@ def interpline(field2d, pivot_point=None,
         end_point_xy = None
         pivot_point_xy = None
         
+        if (latlon is True or is_latlon_pair(start_point) or 
+            is_latlon_pair(pivot_point)):
+            
+            if wrfin is not None:
+                is_moving = is_moving_domain(wrfin)
+            else:
+                is_moving = False
+        
+            if timeidx is None:
+                if wrfin is not None:
+                    # Moving nests aren't supported with ALL_TIMES because the 
+                    # domain could move outside of the line, which causes 
+                    # crashes or different line lengths.
+                    if is_moving:
+                        raise ValueError("Requesting all times with a moving nest "
+                                         "is not supported when using a lat/lon "
+                                         "line because the domain could "
+                                         "move outside of line. "
+                                         "You must request each time "
+                                         "individually.")
+                    else:
+                        # Domain not moving, just use 0
+                        _timeidx = 0
+                
+                # If using grid coordinates, then don't care about lat/lon
+                # coordinates. Just use 0.
+                else:
+                    _timeidx = 0
+            else:
+                if is_moving:
+                    _timeidx = timeidx
+                else:
+                    # When using non-moving nests, set the time to 0 
+                    # to avoid problems downstream
+                    _timeidx = 0
+        
         if pivot_point is not None:
             if pivot_point.lat is not None and pivot_point.lon is not None:
-                xy_coords = to_xy_coords(pivot_point, wrfin, timeidx, 
+                xy_coords = to_xy_coords(pivot_point, wrfin, _timeidx, 
                                          stagger, projection, ll_point)
                 pivot_point_xy = (xy_coords.x, xy_coords.y)
             else:
@@ -462,22 +529,22 @@ def interpline(field2d, pivot_point=None,
                 
         if start_point is not None and end_point is not None:
             if start_point.lat is not None and start_point.lon is not None:
-                xy_coords = to_xy_coords(start_point, wrfin, timeidx, 
+                xy_coords = to_xy_coords(start_point, wrfin, _timeidx, 
                                          stagger, projection, ll_point)
                 start_point_xy = (xy_coords.x, xy_coords.y)
             else:
                 start_point_xy = (start_point.x, start_point.y)
                 
             if end_point.lat is not None and end_point.lon is not None:
-                xy_coords = to_xy_coords(end_point, wrfin, timeidx, 
+                xy_coords = to_xy_coords(end_point, wrfin, _timeidx, 
                                          stagger, projection, ll_point)
                 end_point_xy = (xy_coords.x, xy_coords.y)
             else:
                 end_point_xy = (end_point.x, end_point.y)
-                
+        
         xy = get_xy(field2d, pivot_point_xy, angle, start_point_xy, 
                     end_point_xy)
-        
+    
     return _interpline(field2d, xy)
 
 
@@ -514,10 +581,14 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
             above for the *vert_coord* parameter.
             
         extrapolate (:obj:`bool`, optional): Set to True to extrapolate 
-            values below ground.  This is only performed when the 
-            vertical coordinate type is pressure or height. For temperature 
-            vertical coordinate types, setting this to True will set 
-            values below ground to the lowest model level. Default is False.
+            values below ground.  This is only performed when *vert_coord* is 
+            a pressure or height type, and the *field_type* is a pressure type 
+            (with height vertical coordinate), a height type (with pressure as 
+            the vertical coordinate), or a temperature type (with either height 
+            or pressure as the vertical coordinate). If those conditions are  
+            not met, or *field_type* is None, then the lowest model level 
+            will be used. Extrapolation is performed using standard atmosphere.
+            Default is False. 
             
         field_type (:obj:`str`, optional): 
             The type of field.  Default is None.
@@ -533,8 +604,11 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
                 * 'theta', 'th': potential temperature [K]
                 * 'theta-e', 'thetae', 'eth': equivalent potential temperature
             
-        log_p (:obj:`bool`, optional): Use the log of the pressure for 
-            interpolation instead of pressure. Default is False.
+        log_p (:obj:`bool`, optional): Set to True to use the log of the 
+            vertical coordinate for interpolation. This is mainly intended 
+            for pressure vertical coordinate types, but note that the log 
+            will still be taken for any vertical coordinate type when 
+            this is set to True. Default is False.
             
         timeidx (:obj:`int`, optional):
             The time index to use when extracting auxiallary variables used in 
diff --git a/src/wrf/interputils.py b/src/wrf/interputils.py
index 04807d1..8e04c86 100644
--- a/src/wrf/interputils.py
+++ b/src/wrf/interputils.py
@@ -156,13 +156,6 @@ def _calc_xy(xdim, ydim, pivot_point=None, angle=None,
         if x1 >= xdim or y1 >= ydim:
             raise ValueError("end_point {} is outside of domain "
                              "with shape {}".format(end_point, (xdim, ydim)))
-        
-        # From the original NCL code, but the error above will be thrown
-        # instead.
-        if ( x1 > xdim-1 ): 
-            x1 = xdim - 1
-        if ( y1 > ydim-1): 
-            y1 = ydim - 1
     else:
         raise ValueError("invalid start/end or pivot/angle arguments")
     
diff --git a/src/wrf/latlonutils.py b/src/wrf/latlonutils.py
index adf27ab..e22c1a4 100644
--- a/src/wrf/latlonutils.py
+++ b/src/wrf/latlonutils.py
@@ -167,8 +167,8 @@ def _get_proj_params(wrfin, timeidx, stagger, method, squeeze, cache, _key):
     
     # Only need one file and one time if the domain is not moving
     if not is_moving:
-        if is_multi_time_req(timeidx):
-            lat_timeidx = 0
+        # Always use the 0th time for non-moving domains to avoid problems
+        lat_timeidx = 0
             
         if is_multi_file(wrfin):
             if not is_mapping(wrfin):
@@ -179,7 +179,7 @@ def _get_proj_params(wrfin, timeidx, stagger, method, squeeze, cache, _key):
                 key = _key[first_entry]
                 return _get_proj_params(wrfin, timeidx, stagger, 
                                         method, squeeze, cache, key)
-            
+    
     xlat = extract_vars(wrfin, lat_timeidx, (latvar,), method, squeeze, cache,
                            meta=False, _key=_key)[latvar]
     xlon = extract_vars(wrfin, lat_timeidx, (lonvar,), method, squeeze, cache,
diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py
index 6957713..ea00d75 100644
--- a/src/wrf/metadecorators.py
+++ b/src/wrf/metadecorators.py
@@ -9,7 +9,8 @@ import numpy.ma as ma
 from .extension import _interpline
 from .util import (extract_vars, either, from_args, arg_location,
                    is_coordvar, latlon_coordvars, to_np, 
-                   from_var, iter_left_indexes, is_mapping)
+                   from_var, iter_left_indexes, is_mapping,
+                   is_moving_domain, is_latlon_pair)
 from .coordpair import CoordPair
 from .py3compat import viewkeys, viewitems, py3range
 from .interputils import get_xy_z_params, get_xy, to_xy_coords
@@ -932,9 +933,45 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
     end_point_xy = None
     pivot_point_xy = None
     
+    if (inc_latlon is True or is_latlon_pair(start_point) or 
+        is_latlon_pair(pivot_point)):
+        
+        if wrfin is not None:
+            is_moving = is_moving_domain(wrfin)
+        else:
+            is_moving = False
+    
+        if timeidx is None:
+            if wrfin is not None:
+                # Moving nests aren't supported with ALL_TIMES because the 
+                # domain could move outside of the line, which causes 
+                # crashes or different line lengths.
+                if is_moving:
+                    raise ValueError("Requesting all times with a moving nest "
+                                     "is not supported when using lat/lon "
+                                     "cross sections because the domain could "
+                                     "move outside of the cross section. "
+                                     "You must request each time "
+                                     "individually.")
+                else:
+                    # Domain not moving, just use 0
+                    _timeidx = 0
+            
+            # If using grid coordinates, then don't care about lat/lon
+            # coordinates. Just use 0.
+            else:
+                _timeidx = 0
+        else:
+            if is_moving:
+                _timeidx = timeidx
+            else:
+                # When using non-moving nests, set the time to 0 
+                # to avoid problems downstream
+                _timeidx = 0
+    
     if pivot_point is not None:
         if pivot_point.lat is not None and pivot_point.lon is not None:
-            xy_coords = to_xy_coords(pivot_point, wrfin, timeidx, 
+            xy_coords = to_xy_coords(pivot_point, wrfin, _timeidx, 
                                      stagger, projection, ll_point)
             pivot_point_xy = (xy_coords.x, xy_coords.y)
         else:
@@ -942,14 +979,14 @@ def _set_cross_meta(wrapped, instance, args, kwargs):
             
     if start_point is not None and end_point is not None:
         if start_point.lat is not None and start_point.lon is not None:
-            xy_coords = to_xy_coords(start_point, wrfin, timeidx, 
+            xy_coords = to_xy_coords(start_point, wrfin, _timeidx, 
                                      stagger, projection, ll_point)
             start_point_xy = (xy_coords.x, xy_coords.y)
         else:
             start_point_xy = (start_point.x, start_point.y)
             
         if end_point.lat is not None and end_point.lon is not None:
-            xy_coords = to_xy_coords(end_point, wrfin, timeidx, 
+            xy_coords = to_xy_coords(end_point, wrfin, _timeidx, 
                                      stagger, projection, ll_point)
             end_point_xy = (xy_coords.x, xy_coords.y)
         else:
@@ -1149,10 +1186,47 @@ def _set_line_meta(wrapped, instance, args, kwargs):
     start_point_xy = None
     end_point_xy = None
     pivot_point_xy = None
+    
+    if (inc_latlon is True or is_latlon_pair(start_point) or 
+        is_latlon_pair(pivot_point)):
+        
+        if wrfin is not None:
+            is_moving = is_moving_domain(wrfin)
+        else:
+            is_moving = False
+    
+        if timeidx is None:
+            if wrfin is not None:
+                # Moving nests aren't supported with ALL_TIMES because the 
+                # domain could move outside of the line, which causes 
+                # crashes or different line lengths.
+                if is_moving:
+                    raise ValueError("Requesting all times with a moving nest "
+                                     "is not supported when using a lat/lon "
+                                     "line because the domain could "
+                                     "move outside of line. "
+                                     "You must request each time "
+                                     "individually.")
+                else:
+                    # Domain not moving, just use 0
+                    _timeidx = 0
+            
+            # If using grid coordinates, then don't care about lat/lon
+            # coordinates. Just use 0.
+            else:
+                _timeidx = 0
+        else:
+            if is_moving:
+                _timeidx = timeidx
+            else:
+                # When using non-moving nests, set the time to 0 
+                # to avoid problems downstream
+                _timeidx = 0
+        
         
     if pivot_point is not None:
         if pivot_point.lat is not None and pivot_point.lon is not None:
-            xy_coords = to_xy_coords(pivot_point, wrfin, timeidx, 
+            xy_coords = to_xy_coords(pivot_point, wrfin, _timeidx, 
                                      stagger, projection, ll_point)
             pivot_point_xy = (xy_coords.x, xy_coords.y)
         else:
@@ -1161,19 +1235,20 @@ def _set_line_meta(wrapped, instance, args, kwargs):
             
     if start_point is not None and end_point is not None:
         if start_point.lat is not None and start_point.lon is not None:
-            xy_coords = to_xy_coords(start_point, wrfin, timeidx, 
+            xy_coords = to_xy_coords(start_point, wrfin, _timeidx, 
                                      stagger, projection, ll_point)
             start_point_xy = (xy_coords.x, xy_coords.y)
         else:
             start_point_xy = (start_point.x, start_point.y)
             
         if end_point.lat is not None and end_point.lon is not None:
-            xy_coords = to_xy_coords(end_point, wrfin, timeidx, 
+            xy_coords = to_xy_coords(end_point, wrfin, _timeidx, 
                                      stagger, projection, ll_point)
             end_point_xy = (xy_coords.x, xy_coords.y)
         else:
             end_point_xy = (end_point.x, end_point.y)
     
+
     xy = get_xy(field2d, pivot_point_xy, angle, start_point_xy, end_point_xy)
     
     # Make a copy so we don't modify a user supplied cache
diff --git a/src/wrf/util.py b/src/wrf/util.py
index be174b4..67d30eb 100644
--- a/src/wrf/util.py
+++ b/src/wrf/util.py
@@ -3880,13 +3880,24 @@ def pairs_to_latlon(pairs):
         
         return lats, lons
         
-            
 
+def is_latlon_pair(pair):
+    """Return True if the :class:`wrf.CoordPair` is a lat/lon pair
+    
+    Args:
+    
+        pair (:class:`wrf.CoordPair`): A single :class:`wrf.CoordPair` object.
+        
+    Returns:
+    
+        :obj:`bool`: True if the pair is a lat/lon pair.
+    
+    """
+    if pair is not None:
+        return (pair.lat is not None and pair.lon is not None)
+    else:
+        return False
     
-
-
-
-
     
     
     
diff --git a/test/ci_tests/ci_result_file.nc b/test/ci_tests/ci_result_file.nc
index e15d100..de544e1 100644
Binary files a/test/ci_tests/ci_result_file.nc and b/test/ci_tests/ci_result_file.nc differ
diff --git a/test/ci_tests/ci_test_file.nc b/test/ci_tests/ci_test_file.nc
index b0f8945..5765888 100644
Binary files a/test/ci_tests/ci_test_file.nc and b/test/ci_tests/ci_test_file.nc differ
diff --git a/test/ci_tests/make_test_file.py b/test/ci_tests/make_test_file.py
index fcdfac3..0b1bfd0 100644
--- a/test/ci_tests/make_test_file.py
+++ b/test/ci_tests/make_test_file.py
@@ -146,17 +146,17 @@ def make_result_file(opts):
         for latlonmeth in LATLON_METHS:
             if latlonmeth == "xy":
                 # Hardcoded values from other unit tests
-                lats = [-55, -60, -65]
-                lons = [25, 30, 35]
+                lats = [22.0, 25.0, 27.0]
+                lons = [-90.0, -87.5, -83.75]
                 
                 xy = ll_to_xy(infile, lats[0], lons[0])
                 add_to_ncfile(outfile, xy, "xy")
             else:
                 # Hardcoded from other unit tests
-                i_s = np.asarray([10, 100, 150], int) - 1
-                j_s = np.asarray([10, 100, 150], int) - 1
+                x_s = np.asarray([10, 50, 90], int)
+                y_s = np.asarray([10, 50, 90], int)
                 
-                ll = xy_to_ll(infile, i_s[0], j_s[0])
+                ll = xy_to_ll(infile, x_s[0], y_s[0])
                 add_to_ncfile(outfile, ll, "ll")
 
 def main(opts):
@@ -166,7 +166,8 @@ def main(opts):
     
 if __name__ == "__main__":
     DEFAULT_FILE = ("/Users/ladwig/Documents/wrf_files/"
-                    "wrf_vortex_multi/wrfout_d02_2005-08-28_12:00:00")
+                    "wrf_vortex_multi/moving_nest/"
+                    "wrfout_d02_2005-08-28_12:00:00")
     parser = argparse.ArgumentParser(description="Generate conda test files "
                                      "for unit testing.")
     parser.add_argument("-f", "--filename", required=False,
diff --git a/test/ci_tests/utests.py b/test/ci_tests/utests.py
index a61cf5c..67affd5 100644
--- a/test/ci_tests/utests.py
+++ b/test/ci_tests/utests.py
@@ -193,8 +193,8 @@ def make_latlon_test(testid, wrf_in, referent, single, multi=False, repeat=3,
             # same whether there are multiple or single files
             ref_vals = refnc.variables["xy"][:]
             # Lats/Lons taken from NCL script, just hard-coding for now
-            lats = [-55, -60, -65]
-            lons = [25, 30, 35]
+            lats = [22.0, 25.0, 27.0]
+            lons = [-90.0, -87.5, -83.75]
             
             xy = ll_to_xy(in_wrfnc, lats[0], lons[0])
                 
@@ -208,10 +208,10 @@ def make_latlon_test(testid, wrf_in, referent, single, multi=False, repeat=3,
             
              # i_s, j_s taken from NCL script, just hard-coding for now
              # NCL uses 1-based indexing for this, so need to subtract 1
-            i_s = np.asarray([10, 100, 150], int) - 1
-            j_s = np.asarray([10, 100, 150], int) - 1
+            x_s = np.asarray([10, 50, 90], int)
+            y_s = np.asarray([10, 50, 90], int)
             
-            ll = xy_to_ll(in_wrfnc, i_s[0], j_s[0])
+            ll = xy_to_ll(in_wrfnc, x_s[0], y_s[0])
             
             nt.assert_allclose(to_np(ll), ref_vals)
                 
diff --git a/test/ncl_get_var.ncl b/test/ncl_get_var.ncl
index 4c05ac7..a7b8aa5 100644
--- a/test/ncl_get_var.ncl
+++ b/test/ncl_get_var.ncl
@@ -4,20 +4,29 @@
     
     ;system("printenv")
     
-    if (.not. isvar("in_file")) then
-        in_file = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00.nc"
-        ;in_file = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d02_2005-08-28_00:00:00.nc"
+    if (.not. isvar("dir")) then
+        dir = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest"
+        ;in_file = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00.nc"
+        ;in_file = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest/wrfout_d02_2005-08-28_00:00:00.nc"
     end if 
     
+    if (.not. isvar("pattern")) then
+        pattern = "wrfout_d02_*"
+    end if
+    
     if (.not. isvar("out_file")) then
         out_file = "/tmp/wrftest.nc"
     end if
-    input_file = addfile(in_file,"r")
+    
+    pat = dir + "/" + pattern
+    cmd = "ls " + pat
+    fils = systemfunc (cmd) ; file paths
+    input_file = addfiles(fils,"r")
     
     system("/bin/rm -f " + out_file) ; remove if exists
     fout = addfile(out_file, "c")
     
-    time = 0
+    time = -1
     
     wrf_vars = [/"avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", \
                 "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", \
@@ -48,12 +57,15 @@
     
     
     ; Do the interpolation routines manually
-    ; 3D vertical cross section
-    z := wrf_user_getvar(input_file, "z", 0)        ; grid point height
-    p := wrf_user_getvar(input_file, "pressure", 0) ; total pressure
+    ;;;;;;;;;;;;;;;;;;; 3D vertical cross section
+    time = -1
+    
+    z := wrf_user_getvar(input_file, "z", time)        ; grid point height
+    p := wrf_user_getvar(input_file, "pressure", time) ; total pressure
+    
     
     dimsz = dimsizes(z)
-    pivot = (/ dimsz(2)/2, dimsz(1)/2 /)    ; pivot point is center of domain
+    pivot = (/ dimsz(3)/2, dimsz(2)/2 /)    ; pivot point is center of domain
     
     ht_cross = wrf_user_intrp3d(z,p,"v",pivot,90.0,False)
     
@@ -63,13 +75,16 @@
     
     fout->ht_cross = ht_cross
     fout->p_cross = p_cross
-;    
+ 
+
+    time = 0
+    
     ; For the new cross section routine
     xopt = True
     xopt@use_pivot = True
     xopt@angle = 90.0
     xopt@file_handle = input_file
-    xopt@timeidx = 0
+    xopt@timeidx = time
     xopt@linecoords = True
     
     ht_vertcross1 = wrf_user_vertcross(z, p, pivot, xopt)
@@ -82,25 +97,29 @@
     xopt@angle = 90.0
     xopt@levels = (/1000., 850., 700., 500., 250./)
     xopt@file_handle = input_file
-    xopt@timeidx = 0
+    xopt@timeidx = time
     xopt@linecoords = True
     
     ht_vertcross2 = wrf_user_vertcross(z, p, pivot, xopt)
-    ht_vertcross2!0 = "vertical2"
-    ht_vertcross2!1 = "cross_line_idx2"
-    
+    ht_vertcross2!1 = "vertical2"
+    ht_vertcross2!2 = "cross_line_idx2"
     
     fout->ht_vertcross2 = ht_vertcross2
     
+    ; Can only use a single time for lat/lon version at this time
+    
+    vertdims = dimsizes(ht_vertcross2)
+    htdims = dimsizes(z)
+    
     lats = wrf_user_getvar(input_file, "lat", 0)
     lons = wrf_user_getvar(input_file, "lon", 0) 
-    
+        
     start_lat = min(lats) + .25d*(max(lats) - min(lats))
-    end_lat = min(lats) + .75d*(max(lats) - min(lats))
-            
+    end_lat = min(lats) + .65d*(max(lats) - min(lats))
+                
     start_lon = min(lons) + .25d*(max(lons) - min(lons))
-    end_lon = min(lons) + .75d*(max(lons) - min(lons))
-    
+    end_lon = min(lons) + .65d*(max(lons) - min(lons))
+        
     start_end = (/ start_lon, start_lat, end_lon, end_lat /)
     
     ; For the new cross section routine
@@ -111,14 +130,42 @@
     xopt@timeidx = 0
     xopt@linecoords = True
     xopt@autolevels = 1000
-    
+        
     ht_vertcross3 = wrf_user_vertcross(z, p, start_end, xopt)
-    ht_vertcross3!0 = "vertical3"
-    ht_vertcross3!1 = "cross_line_idx3"
+    
+    ht_vertcross3!0 = "Time"
+    ht_vertcross3!1 = "vertical3"
+    ht_vertcross3!2 = "cross_line_idx3"
     
     fout->ht_vertcross3 = ht_vertcross3
     
-    ; 3D horizontal interpolation
+    ; Test the moving nest with lat/lon over time
+    
+    times = wrf_user_getvar(input_file, "times", -1)
+    ntimes = dimsizes(times)
+    
+    do i=0,ntimes-1
+      xopt@timeidx = i
+      name = sprinti("ht_vertcross_t%i", i)
+      p_var := p(i,:,:,:)
+      z_var := z(i,:,:,:)
+
+      ht_vertcross := wrf_user_vertcross(z_var, p_var, start_end, xopt)
+      
+      dim0name = sprinti("vertical_t%i",i)
+      dim1name = sprinti("cross_line_idx_t%i",i)
+      ht_vertcross!0 = dim0name
+      ht_vertcross!1 = dim1name
+      
+      fout->$name$ = ht_vertcross
+    end do
+    
+    ;;;;;;;;;;;;;;;;;;;;;;;; 3D horizontal interpolation
+    
+    time = -1
+    
+    z := wrf_user_getvar(input_file, "z", time)        ; grid point height
+    p := wrf_user_getvar(input_file, "pressure", time) ; total pressure
     
     ; First, check backwards compat
     plev = 500.   ; 500 MB
@@ -158,7 +205,7 @@
     fout->z2_multi = z2_multi
     fout->p2_multi = p2_multi
     
-    pblh = wrf_user_getvar(input_file, "PBLH", 0)
+    pblh = wrf_user_getvar(input_file, "PBLH", time)
     opts := False
     opts@inc2dlevs = True
     p_lev2d = wrf_user_interplevel(p, z, pblh, opts)
@@ -166,19 +213,90 @@
     fout->p_lev2d = p_lev2d
     
     
-    ; 2D interpolation along line
-    t2 = wrf_user_getvar(input_file, "T2", 0)
+    ;;;;;;;;;;;;;;;;;;;;;;;; 2D interpolation along line
+    
+    time = -1
+    
+    t2 = wrf_user_getvar(input_file, "T2", time)
     dimst2 = dimsizes(t2)
-    pivot = (/ dimst2(1)/2, dimst2(0)/2 /)
+    pivot = (/ dimst2(2)/2, dimst2(1)/2 /)
     
     t2_line = wrf_user_intrp2d(t2, pivot, 90.0, False)
     
-    fout->t2_line = (/t2_line/)
+    fout->t2_line = t2_line
+    
+    ; For the new interplevel routine
+    xopt := True
+    xopt@use_pivot = True
+    xopt@angle = 90.0
+    xopt@latlon = False
+    xopt@file_handle = input_file
+    xopt@timeidx = 0
+    xopt@linecoords = True
+    
+    t2_line2 = wrf_user_interpline(t2, pivot, xopt)
     
+    fout->t2_line2 = t2_line2
+    
+    lats = wrf_user_getvar(input_file, "lat", 0)
+    lons = wrf_user_getvar(input_file, "lon", 0) 
+        
+    start_lat = min(lats) + .25d*(max(lats) - min(lats))
+    end_lat = min(lats) + .65d*(max(lats) - min(lats))
+                
+    start_lon = min(lons) + .25d*(max(lons) - min(lons))
+    end_lon = min(lons) + .65d*(max(lons) - min(lons))
+        
+    start_end = (/ start_lon, start_lat, end_lon, end_lat /)
+    
+    ; For the new line routine
+    xopt := True
+    xopt@use_pivot = False
+    xopt@latlon = True
+    xopt@file_handle = input_file
+    xopt@timeidx = 0
+    xopt@linecoords = True
+    
+    t2_line3 = wrf_user_interpline(t2, start_end, xopt)
+    t2_line3!1 = "line_idx_t2_line3"
+    
+    fout->t2_line3 = t2_line3
+    
+    times = wrf_user_getvar(input_file, "times", -1)
+    ntimes = dimsizes(times)
+    
+    do i=0,ntimes-1
+      xopt@timeidx = i
+      name = sprinti("t2_line_t%i", i)
+      dim0name = sprinti("lineidx_t%i",i)
+      var := t2(i,:,:)
+      t2_line := wrf_user_interpline(var, start_end, xopt)
+      t2_line!0 = dim0name
+      fout->$name$ = t2_line
+    end do
+    
+    ; Make sure the 1 time case still works
+    t2 := wrf_user_getvar(input_file, "T2", 0)
+    
+    ; For the new line routine
+    xopt := True
+    xopt@use_pivot = False
+    xopt@latlon = True
+    xopt@file_handle = input_file
+    xopt@timeidx = 0
+    xopt@linecoords = True
+    
+    t2_line4 = wrf_user_interpline(t2, start_end, xopt)
+    t2_line4!0 = "t2_line4_idx"
+    
+    fout->t2_line4 = t2_line4
+    
+    
+    ;;;;;;;;;;;;;;;;;;;;;;; 3D interpolation to new vertical coordinates
+    time = -1
     
-    ; 3D interpolation to new vertical coordinates
     ; interp t to theta
-    fld1 = wrf_user_getvar(input_file, "tk", -1)
+    fld1 = wrf_user_getvar(input_file, "tk", time)
     vert_coord       = "theta"
     interp_levels    = (/200,300,500,1000/)
  
@@ -231,7 +349,7 @@
     fout->fld_tk_ght_agl = fld5_intrp
     
     ; interp ht to pres
-    fld6 = wrf_user_getvar(input_file, "height", -1)
+    fld6 = wrf_user_getvar(input_file, "height", time)
     vert_coord := "pressure"
     opts@field_type  = "ght"
     interp_levels    := (/500,50/)
@@ -242,7 +360,7 @@
     
     
     ; interp pres to theta
-    fld7 = wrf_user_getvar(input_file, "pressure", -1)
+    fld7 = wrf_user_getvar(input_file, "pressure", time)
     vert_coord := "theta"
     opts@field_type  = "pressure"
     interp_levels    := (/200,300,500,1000/)
@@ -253,7 +371,7 @@
     
     
     ; interp theta-e to pressure
-    fld8 = wrf_user_getvar(input_file, "eth", -1)
+    fld8 = wrf_user_getvar(input_file, "eth", time)
     vert_coord := "pressure"
     opts@field_type  = "T"
     interp_levels    := (/850,500,5/)
@@ -263,16 +381,32 @@
     fout->fld_thetae_pres = fld8_intrp
     
     
-    ; lat/lon to x/y and x/y to lat/lon routines
-    lats := (/-55, -60, -65 /)
-    lons := (/25, 30, 35 /)
-    i_s = (/10, 100, 150 /)
-    j_s = (/10, 100, 150 /)
-    ij = wrf_user_ll_to_ij(input_file, lons, lats, True) 
-    ll = wrf_user_ij_to_ll(input_file, i_s, j_s, True)
+    ;;;;;;;;;;;;;;;;;;; lat/lon to x/y and x/y to lat/lon routines
+    
+    lats := (/22.0, 25.0, 27.0 /)
+    lons := (/-90.0, -87.5, -83.75 /)
+    x_s = (/10, 50, 90 /)
+    y_s = (/10, 50, 90 /)
+    
+    opt = True
+    opt@useTime = -1
+    opt@returnInt = False
+    
+    xy1 = wrf_user_ll_to_xy(input_file, lons, lats, opt) 
+    ll1 = wrf_user_xy_to_ll(input_file, x_s, y_s, opt)
+    
+    fout->xy1 = xy1
+    fout->ll1 = ll1
+    
+    opt = True
+    opt@useTime = 8
+    opt@returnInt = True
+    
+    xy2 = wrf_user_ll_to_xy(input_file, lons, lats, opt) 
+    ll2 = wrf_user_xy_to_ll(input_file, x_s, y_s, opt)
     
-    fout->ij = ij
-    fout->ll = ll
+    fout->xy2 = xy2
+    fout->ll2 = ll2
     
     delete(fout)
 
diff --git a/test/utests.py b/test/utests.py
index 0febdd4..613b13e 100644
--- a/test/utests.py
+++ b/test/utests.py
@@ -4,6 +4,7 @@ import numpy as np
 import numpy.ma as ma
 import os, sys
 import subprocess
+import glob
 
 from wrf import (getvar, interplevel, interpline, vertcross, vinterp,
                  disable_xarray, xarray_enabled, to_np,
@@ -11,44 +12,51 @@ from wrf import (getvar, interplevel, interpline, vertcross, vinterp,
                  extract_global_attrs, viewitems, CoordPair, ll_points)
 from wrf.util import is_multi_file
 
-NCL_EXE = "/Users/ladwig/nclbuild/6.3.0/bin/ncl"
-TEST_FILE = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00"
-OUT_NC_FILE = "/tmp/wrftest.nc"
+NCL_EXE = "/Users/ladwig/miniconda2/envs/ncl_build/bin/ncl"
+NCARG_ROOT = "/Users/ladwig/miniconda2/envs/ncl_build"
+#TEST_FILE = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00"
+DIRS = ["/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest",
+        "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/static_nest"]
+PATTERN = "wrfout_d02_*"
+REF_NC_FILES = ["/tmp/wrftest_moving.nc", "/tmp/wrftest_static.nc"]
+NEST = ["moving", "static"]
 
 # Python 3
 if sys.version_info > (3,):
     xrange = range
 
 def setUpModule():
-    ncarg_root = os.environ.get("NCARG_ROOT", None)
-    if ncarg_root is None:
-        raise RuntimeError("$NCARG_ROOT environment variable not set")
-        
+    #ncarg_root = os.environ.get("NCARG_ROOT", None)
+    #if ncarg_root is None:
+    #    raise RuntimeError("$NCARG_ROOT environment variable not set")
     
+    os.environ["NCARG_ROOT"] = NCARG_ROOT
+    os.environ["NCARG_NCARG"] = os.path.join(NCARG_ROOT, "lib", "ncarg")
+    os.environ["OMP_NUM_THREADS"] = "4"
+     
+        
     this_path = os.path.realpath(__file__)
     ncl_script = os.path.join(os.path.dirname(this_path),
                               "ncl_get_var.ncl")
-    ncfile = TEST_FILE + ".nc" # NCL requires extension
     
-    # This needs to be set when PyNIO is installed, since PyNIOs data does
-    # not contain the dat file for the CAPE calcluations
-    os.environ["NCARG_NCARG"] = os.path.join(os.environ["NCARG_ROOT"],
-                                             "lib", "ncarg")
-    cmd = "%s %s 'in_file=\"%s\"' 'out_file=\"%s\"'" % (NCL_EXE,
-                                                        ncl_script,
-                                                        ncfile,
-                                                        OUT_NC_FILE)
+    for dir,outfile in zip(DIRS, REF_NC_FILES):
+        cmd = "%s %s 'dir=\"%s\"' 'pattern=\"%s\"'  'out_file=\"%s\"'" % (
+            NCL_EXE,
+            ncl_script,
+            dir,
+            PATTERN,
+            outfile)
     
-    #print cmd
+        print cmd
     
-    if not os.path.exists(OUT_NC_FILE):
-        status = subprocess.call(cmd, shell=True)
-        if (status != 0):
-            raise RuntimeError("NCL script failed.  Could not set up test.")
+        if not os.path.exists(outfile):
+            status = subprocess.call(cmd, shell=True)
+            if (status != 0):
+                raise RuntimeError("NCL script failed. Could not set up test.")
 
 # Using helpful information at: 
 # http://eli.thegreenplace.net/2014/04/02/dynamically-generating-python-test-cases
-def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
+def make_test(varname, dir, pattern, referent, multi=False, pynio=False):
     def test(self):
 
         try:
@@ -61,40 +69,28 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
         except:
             pass
         
-        if not multi:
-            timeidx = 0
-            if not pynio:
-                in_wrfnc = NetCDF(wrf_in)
-                try:
-                    in_wrfnc.set_auto_mask(False)
-                except:
-                    pass
-            else:
-                # Note: Python doesn't like it if you reassign an outer scoped
-                # variable (wrf_in in this case)
-                if not wrf_in.endswith(".nc"):
-                    wrf_file = wrf_in + ".nc"
-                else:
-                    wrf_file = wrf_in
-                in_wrfnc = Nio.open_file(wrf_file)
-        else:
-            timeidx = None
+        timeidx = 0 if not multi else None
+        pat = os.path.join(dir, pattern)
+        wrf_files = glob.glob(pat)
+        wrf_files.sort()
+        
+        wrfin = []
+        for fname in wrf_files:
             if not pynio:
-                nc = NetCDF(wrf_in)
+                f = NetCDF(fname)
                 try:
-                    nc.set_auto_mask(False)
+                    f.set_always_mask(False)
                 except:
                     pass
-                in_wrfnc = [nc for i in xrange(repeat)]
+                wrfin.append(f)
             else:
-                if not wrf_in.endswith(".nc"):
-                    wrf_file = wrf_in + ".nc"
+                if not fname.endswith(".nc"):
+                    _fname = fname + ".nc"
                 else:
-                    wrf_file = wrf_in
-                nc = Nio.open_file(wrf_file)
-                in_wrfnc = [nc for i in xrange(repeat)]
-            
-        
+                    _fname = fname
+                f = Nio.open_file(_fname)
+                wrfin.append(f)
+
         refnc = NetCDF(referent)
         try:
             refnc.set_auto_mask(False)
@@ -104,69 +100,42 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
         # These have a left index that defines the product type
         multiproduct = varname in ("uvmet", "uvmet10", "cape_2d", "cape_3d", 
                                    "cfrac")
+        multi2d = ("uvmet10", "cape_2d", "cfrac")
+        multi3d = ("uvmet", "cape_3d")
         
         # These varnames don't have NCL functions to test against
         ignore_referent = ("zstag", "geopt_stag")
         
         if varname not in ignore_referent:
             if not multi:
-                ref_vals = refnc.variables[varname][:]
-            else:
-                data = refnc.variables[varname][:]
-                if (varname != "uvmet" and varname != "uvmet10" 
-                    and varname != "cape_2d" and varname != "cape_3d"):
-                    new_dims = [repeat] + [x for x in data.shape]
-                elif (varname == "uvmet" or varname == "uvmet10" 
-                      or varname == "cape_3d"):
-                    new_dims = [2] + [repeat] + [x for x in data.shape[1:]]
-                elif (varname == "cape_2d"):
-                    new_dims = [4] + [repeat] + [x for x in data.shape[1:]]
-                elif (varname == "cfrac"):
-                    new_dims = [3] + [repeat] + [x for x in data.shape[1:]]
-    
-                    
-                masked=False
-                if (isinstance(data, ma.core.MaskedArray)):
-                    masked=True
-                
-                if not masked:
-                    ref_vals = np.zeros(new_dims, data.dtype)
+                if varname in multi2d:
+                    ref_vals = refnc.variables[varname][...,0,:,:]
+                elif varname in multi3d:
+                    ref_vals = refnc.variables[varname][...,0,:,:,:]
                 else:
-                    ref_vals = ma.asarray(np.zeros(new_dims, data.dtype))
-                
-                for i in xrange(repeat):
-                    if not multiproduct:
-                        ref_vals[i,:] = data[:]
-                    
-                        if masked:
-                            ref_vals.mask[i,:] = data.mask[:]
-                            
-                    else:
-                        for prod in xrange(ref_vals.shape[0]):
-                            ref_vals[prod,i,:] = data[prod,:]
-                            
-                            if masked:
-                                ref_vals.mask[prod,i,:] = data.mask[prod,:]
+                    ref_vals = refnc.variables[varname][0,:]
+            else:
+                ref_vals = refnc.variables[varname][:]
         
         if (varname == "tc"):
-            my_vals = getvar(in_wrfnc, "temp", timeidx=timeidx, units="c")
+            my_vals = getvar(wrfin, "temp", timeidx=timeidx, units="c")
             tol = 1/100.
             atol = .1 # Note:  NCL uses 273.16 as conversion for some reason
             nt.assert_allclose(to_np(my_vals), ref_vals, tol, atol)
         elif (varname == "height_agl"):
             # Change the vert_type to height_agl when NCL gets updated.
-            my_vals = getvar(in_wrfnc, "z", timeidx=timeidx, msl=False)
+            my_vals = getvar(wrfin, "z", timeidx=timeidx, msl=False)
             tol = 1/100.
             atol = .1 # Note:  NCL uses 273.16 as conversion for some reason
             nt.assert_allclose(to_np(my_vals), ref_vals, tol, atol)
         elif (varname == "cfrac"):
             # Change the vert_type to height_agl when NCL gets updated.
-            my_vals = getvar(in_wrfnc, "cfrac", timeidx=timeidx)
+            my_vals = getvar(wrfin, "cfrac", timeidx=timeidx)
             tol = 1/100.
             atol = .1 # Note:  NCL uses 273.16 as conversion for some reason
             nt.assert_allclose(to_np(my_vals), ref_vals, tol, atol)
         elif (varname == "pw"):
-            my_vals = getvar(in_wrfnc, "pw", timeidx=timeidx)
+            my_vals = getvar(wrfin, "pw", timeidx=timeidx)
             tol = .5/100.0
             atol = 0 # NCL uses different constants and doesn't use same
                      # handrolled virtual temp in method
@@ -176,7 +145,7 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
                 print (np.amax(np.abs(to_np(my_vals) - ref_vals)))
                 raise
         elif (varname == "cape_2d"):
-            cape_2d = getvar(in_wrfnc, varname, timeidx=timeidx)
+            cape_2d = getvar(wrfin, varname, timeidx=timeidx)
             tol = 0/100.
             atol = 200.0
             # Let's only compare CAPE values until the F90 changes are 
@@ -185,7 +154,7 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
             # causing wildly different values in LCL
             nt.assert_allclose(to_np(cape_2d[0,:]), ref_vals[0,:], tol, atol)
         elif (varname == "cape_3d"):
-            cape_3d = getvar(in_wrfnc, varname, timeidx=timeidx)
+            cape_3d = getvar(wrfin, varname, timeidx=timeidx)
             # Changing the R and CP constants, while keeping TK within
             # 2%, can lead to some big changes in CAPE.  Tolerances 
             # have been set wide when comparing the with the original
@@ -197,11 +166,11 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
             #print np.amax(np.abs(to_np(cape_3d[0,:]) - ref_vals[0,:]))
             nt.assert_allclose(to_np(cape_3d), ref_vals, tol, atol)
         elif (varname == "zstag" or varname == "geopt_stag"):
-            v = getvar(in_wrfnc, varname, timeidx=timeidx)
+            v = getvar(wrfin, varname, timeidx=timeidx)
             # For now, only make sure it runs without crashing since no NCL
             # to compare with yet.
         else:
-            my_vals = getvar(in_wrfnc, varname, timeidx=timeidx)
+            my_vals = getvar(wrfin, varname, timeidx=timeidx)
             tol = 2/100.
             atol = 0.1
             #print (np.amax(np.abs(to_np(my_vals) - ref_vals)))
@@ -218,7 +187,7 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False):
     
     return test
 
-def _get_refvals(referent, varname, repeat, multi):
+def _get_refvals(referent, varname, multi):
     try:
         from netCDF4 import Dataset as NetCDF
     except:
@@ -226,34 +195,36 @@ def _get_refvals(referent, varname, repeat, multi):
         
     refnc = NetCDF(referent)
     try:
-        refnc.set_auto_mask(False)
+        pass
+        #refnc.set_auto_mask(False)
     except:
         pass
     
+    multi2d = ("uvmet10", "cape_2d", "cfrac")
+    multi3d = ("uvmet", "cape_3d")
+    interpline = ("t2_line","t2_line2", "t2_line3")
+    
     if not multi:
-        ref_vals = refnc.variables[varname][:]
-    else:
-        data = refnc.variables[varname][:]
-        new_dims = [repeat] + [x for x in data.shape]
-        masked=False
-        if (isinstance(data, ma.core.MaskedArray)):
-            masked=True
-          
-        if not masked:
-            ref_vals = np.zeros(new_dims, data.dtype)
+        if varname in multi2d:
+            ref_vals = refnc.variables[varname][...,0,:,:]
+        elif varname in multi3d:
+            ref_vals = refnc.variables[varname][...,0,:,:,:]
         else:
-            ref_vals = ma.asarray(np.zeros(new_dims, data.dtype))
-              
-        for i in xrange(repeat):
-            ref_vals[i,:] = data[:]
-              
-            if masked:
-                ref_vals.mask[i,:] = data.mask[:]
+            v = refnc.variables[varname][:]
+            if v.ndim == 2:
+                if varname in interpline:
+                    ref_vals = v[0,:]
+                else:
+                    ref_vals = v
+            else:
+                ref_vals = v[0,:]
+    else:
+        ref_vals = refnc.variables[varname][:]
                 
     return ref_vals
 
-def make_interp_test(varname, wrf_in, referent, multi=False, 
-                     repeat=3, pynio=False):
+def make_interp_test(varname, dir, pattern, referent, multi=False, 
+                     pynio=False):
     def test(self):
         try:
             from netCDF4 import Dataset as NetCDF
@@ -265,55 +236,44 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
         except:
             pass
         
-        if not multi:
-            timeidx = 0
-            if not pynio:
-                in_wrfnc = NetCDF(wrf_in)
-                try:
-                    in_wrfnc.set_auto_mask(False)
-                except:
-                    pass
-            else:
-                # Note: Python doesn't like it if you reassign an outer scoped
-                # variable (wrf_in in this case)
-                if not wrf_in.endswith(".nc"):
-                    wrf_file = wrf_in + ".nc"
-                else:
-                    wrf_file = wrf_in
-                in_wrfnc = Nio.open_file(wrf_file)
-        else:
-            timeidx = None
+        timeidx = 0 if not multi else None
+        pat = os.path.join(dir, pattern)
+        wrf_files = glob.glob(pat)
+        wrf_files.sort()
+        
+        wrfin = []
+        for fname in wrf_files:
             if not pynio:
-                nc = NetCDF(wrf_in)
+                f = NetCDF(fname)
                 try:
-                    nc.set_auto_mask(False)
+                    f.set_always_mask(False)
                 except:
                     pass
-                in_wrfnc = [nc for i in xrange(repeat)]
+                wrfin.append(f)
             else:
-                if not wrf_in.endswith(".nc"):
-                    wrf_file = wrf_in + ".nc"
+                if not fname.endswith(".nc"):
+                    _fname = fname + ".nc"
                 else:
-                    wrf_file = wrf_in
-                nc = Nio.open_file(wrf_file)
-                in_wrfnc = [nc for i in xrange(repeat)]
+                    _fname = fname
+                f = Nio.open_file(_fname)
+                wrfin.append(f)
         
         if (varname == "interplevel"):
-            ref_ht_500 = _get_refvals(referent, "z_500", repeat, multi)
-            ref_p_5000 = _get_refvals(referent, "p_5000", repeat, multi)
-            ref_ht_multi = _get_refvals(referent, "z_multi", repeat, multi)
-            ref_p_multi = _get_refvals(referent, "p_multi", repeat, multi)
+            ref_ht_500 = _get_refvals(referent, "z_500", multi)
+            ref_p_5000 = _get_refvals(referent, "p_5000", multi)
+            ref_ht_multi = _get_refvals(referent, "z_multi", multi)
+            ref_p_multi = _get_refvals(referent, "p_multi", multi)
             
-            ref_ht2_500 = _get_refvals(referent, "z2_500", repeat, multi)
-            ref_p2_5000 = _get_refvals(referent, "p2_5000", repeat, multi)
-            ref_ht2_multi = _get_refvals(referent, "z2_multi", repeat, multi)
-            ref_p2_multi = _get_refvals(referent, "p2_multi", repeat, multi)
+            ref_ht2_500 = _get_refvals(referent, "z2_500", multi)
+            ref_p2_5000 = _get_refvals(referent, "p2_5000", multi)
+            ref_ht2_multi = _get_refvals(referent, "z2_multi", multi)
+            ref_p2_multi = _get_refvals(referent, "p2_multi", multi)
             
-            ref_p_lev2d = _get_refvals(referent, "p_lev2d", repeat, multi)
+            ref_p_lev2d = _get_refvals(referent, "p_lev2d", multi)
             
-            hts = getvar(in_wrfnc, "z", timeidx=timeidx)
-            p = getvar(in_wrfnc, "pressure", timeidx=timeidx)
-            wspd_wdir = getvar(in_wrfnc, "wspd_wdir", timeidx=timeidx)
+            hts = getvar(wrfin, "z", timeidx=timeidx)
+            p = getvar(wrfin, "pressure", timeidx=timeidx)
+            wspd_wdir = getvar(wrfin, "wspd_wdir", timeidx=timeidx)
             
             # Make sure the numpy versions work first
             hts_500 = interplevel(to_np(hts), to_np(p), 500)
@@ -346,7 +306,7 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             nt.assert_allclose(to_np(p_multi), ref_p_multi)
             nt.assert_allclose(to_np(p_multi), ref_p2_multi)
             
-            pblh = getvar(in_wrfnc, "PBLH", timeidx=timeidx)
+            pblh = getvar(wrfin, "PBLH", timeidx=timeidx)
             p_lev2d = interplevel(to_np(p), to_np(hts), to_np(pblh))
             p_lev2d = interplevel(p, hts, pblh)
             nt.assert_allclose(to_np(p_lev2d), ref_p_lev2d)
@@ -382,17 +342,17 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
                 
             
         elif (varname == "vertcross"):
-            ref_ht_cross = _get_refvals(referent, "ht_cross", repeat, multi)
-            ref_p_cross = _get_refvals(referent, "p_cross", repeat, multi)
-            ref_ht_vertcross1 = _get_refvals(referent, "ht_vertcross1", repeat, 
+            ref_ht_cross = _get_refvals(referent, "ht_cross", multi)
+            ref_p_cross = _get_refvals(referent, "p_cross", multi)
+            ref_ht_vertcross1 = _get_refvals(referent, "ht_vertcross1", 
                                             multi)
-            ref_ht_vertcross2 = _get_refvals(referent, "ht_vertcross2", repeat, 
+            ref_ht_vertcross2 = _get_refvals(referent, "ht_vertcross2", 
                                             multi)
-            ref_ht_vertcross3 = _get_refvals(referent, "ht_vertcross3", repeat, 
+            ref_ht_vertcross3 = _get_refvals(referent, "ht_vertcross3", 
                                             multi)
             
-            hts = getvar(in_wrfnc, "z", timeidx=timeidx)
-            p = getvar(in_wrfnc, "pressure", timeidx=timeidx)
+            hts = getvar(wrfin, "z", timeidx=timeidx)
+            p = getvar(wrfin, "pressure", timeidx=timeidx)
             
             pivot_point = CoordPair(hts.shape[-1] / 2, hts.shape[-2] / 2) 
             
@@ -403,33 +363,43 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             p_min = np.ceil(np.amin(p)/10) * 10      # top value
             
             p_autolevels = int((p_max - p_min) /10)
-
+            
             # Make sure the numpy versions work first
             
             ht_cross = vertcross(to_np(hts), to_np(p), 
                                  pivot_point=pivot_point, angle=90.,
                                  autolevels=p_autolevels)
+            
             ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90.,
                                  autolevels=p_autolevels)
             
             nt.assert_allclose(to_np(ht_cross), ref_ht_cross, rtol=.01)
             
-            # Test the manual projection method with lat/lon
             lats = hts.coords["XLAT"]
             lons = hts.coords["XLONG"]
-            ll_point = ll_points(lats, lons)
-            pivot = CoordPair(lat=lats[int(lats.shape[-2]/2), 
+                
+            # Test the manual projection method with lat/lon
+            # Only do this for the non-multi case, since the domain 
+            # might be moving
+            if not multi:
+                if lats.ndim > 2: # moving nest
+                    lats = lats[0,:]
+                    lons = lons[0,:]
+                    
+                ll_point = ll_points(lats, lons)
+            
+                pivot = CoordPair(lat=lats[int(lats.shape[-2]/2), 
                                        int(lats.shape[-1]/2)],
-                              lon=lons[int(lons.shape[-2]/2), 
+                                  lon=lons[int(lons.shape[-2]/2), 
                                        int(lons.shape[-1]/2)])
             
-            v1 = vertcross(hts,p,wrfin=in_wrfnc,pivot_point=pivot_point,
+                v1 = vertcross(hts,p,wrfin=wrfin,pivot_point=pivot_point,
                            angle=90.0)
-            v2 = vertcross(hts,p,projection=hts.attrs["projection"], 
+                v2 = vertcross(hts,p,projection=hts.attrs["projection"], 
                           ll_point=ll_point,
                           pivot_point=pivot_point, angle=90.)
             
-            nt.assert_allclose(to_np(v1), to_np(v2), rtol=.01)
+                nt.assert_allclose(to_np(v1), to_np(v2), rtol=.01)
             
             # Test opposite
             
@@ -440,7 +410,7 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
                                rtol=.01)
             # Test point to point
             start_point = CoordPair(0, hts.shape[-2]/2)
-            end_point = CoordPair(-1,hts.shape[-2]/2)
+            end_point = CoordPair(-1, hts.shape[-2]/2)
           
              
             p_cross2 = vertcross(p,hts,start_point=start_point, 
@@ -468,17 +438,19 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             nt.assert_allclose(to_np(ht_cross), 
                                to_np(ref_ht_vertcross2), atol=.01)
             
+            idxs = (0, slice(None)) if lats.ndim > 2 else (slice(None),)
             
-            start_lat = np.amin(lats) + .25*(np.amax(lats) - np.amin(lats))
-            end_lat = np.amin(lats) + .75*(np.amax(lats) - np.amin(lats))
+            start_lat = np.amin(lats[idxs]) + .25*(np.amax(lats[idxs]) - np.amin(lats[idxs]))
+            end_lat = np.amin(lats[idxs]) + .65*(np.amax(lats[idxs]) - np.amin(lats[idxs]))
             
-            start_lon = np.amin(lons) + .25*(np.amax(lons) - np.amin(lons))
-            end_lon = np.amin(lons) + .75*(np.amax(lons) - np.amin(lons))
+            start_lon = np.amin(lons[idxs]) + .25*(np.amax(lons[idxs]) - np.amin(lons[idxs]))
+            end_lon = np.amin(lons[idxs]) + .65*(np.amax(lons[idxs]) - np.amin(lons[idxs]))
             
             start_point = CoordPair(lat=start_lat, lon=start_lon)
             end_point = CoordPair(lat=end_lat, lon=end_lon)
             
-            # ll_point and projection came from above 
+            ll_point = ll_points(lats[idxs], lons[idxs])
+            
             ht_cross = vertcross(hts, p, 
                                  start_point=start_point, 
                                  end_point=end_point, 
@@ -492,12 +464,35 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
                                to_np(ref_ht_vertcross3),
                                rtol=.01)
             
+            if multi:
+                ntimes = hts.shape[0]
+                
+                for t in range(ntimes):
+                    hts = getvar(wrfin, "z", timeidx=t)
+                    p = getvar(wrfin, "pressure", timeidx=t)
+                    
+                    ht_cross = vertcross(hts, p, 
+                                 start_point=start_point, 
+                                 end_point=end_point, 
+                                 wrfin=wrfin,
+                                 timeidx=t,
+                                 latlon=True,
+                                 autolevels=1000)
+                    
+                    refname = "ht_vertcross_t{}".format(t)
+                    ref_ht_vertcross = _get_refvals(referent, refname, False)
+                    
+                    nt.assert_allclose(to_np(ht_cross), 
+                               to_np(ref_ht_vertcross),rtol=.02)
+            
               
         elif (varname == "interpline"):
             
-            ref_t2_line = _get_refvals(referent, "t2_line", repeat, multi)
+            ref_t2_line = _get_refvals(referent, "t2_line", multi)
+            ref_t2_line2 = _get_refvals(referent, "t2_line2", multi)
+            ref_t2_line3 = _get_refvals(referent, "t2_line3", multi)
             
-            t2 = getvar(in_wrfnc, "T2", timeidx=timeidx)
+            t2 = getvar(wrfin, "T2", timeidx=timeidx)
             pivot_point = CoordPair(t2.shape[-1] / 2, t2.shape[-2] / 2)
             
             # Make sure the numpy version works
@@ -507,15 +502,24 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             
             nt.assert_allclose(to_np(t2_line1), ref_t2_line)
             
+            # Test the new NCL wrf_user_interplevel result
+            nt.assert_allclose(to_np(t2_line1), ref_t2_line2)
+            
             # Test the manual projection method with lat/lon
             lats = t2.coords["XLAT"]
             lons = t2.coords["XLONG"]
+            if multi:
+                if lats.ndim > 2: # moving nest
+                    lats = lats[0,:]
+                    lons = lons[0,:]
+            
             ll_point = ll_points(lats, lons)
+            
             pivot = CoordPair(lat=lats[int(lats.shape[-2]/2), 
                                        int(lats.shape[-1]/2)],
                               lon=lons[int(lons.shape[-2]/2), 
                                        int(lons.shape[-1]/2)])
-            l1 = interpline(t2,wrfin=in_wrfnc,pivot_point=pivot_point,
+            l1 = interpline(t2,wrfin,pivot_point=pivot_point,
                            angle=90.0)
             l2 = interpline(t2,projection=t2.attrs["projection"], 
                           ll_point=ll_point, 
@@ -530,17 +534,74 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
                                   end_point=end_point)
             
             nt.assert_allclose(to_np(t2_line1), to_np(t2_line2))
+            
+            # Now test the start/end with lat/lon points
+            
+            start_lat = float(np.amin(lats) + .25*(np.amax(lats) 
+                                                   - np.amin(lats)))
+            end_lat = float(np.amin(lats) + .65*(np.amax(lats) 
+                                                 - np.amin(lats)))
+            
+            start_lon = float(np.amin(lons) + .25*(np.amax(lons) 
+                                                   - np.amin(lons)))
+            end_lon = float(np.amin(lons) + .65*(np.amax(lons) 
+                                                 - np.amin(lons)))
+            
+            start_point = CoordPair(lat=start_lat, lon=start_lon)
+            end_point = CoordPair(lat=end_lat, lon=end_lon)
+            
+            t2_line3 = interpline(t2,wrfin=wrfin,timeidx=0,
+                                  start_point=start_point,
+                                  end_point=end_point,latlon=True)
+            
+            
+            nt.assert_allclose(to_np(t2_line3), ref_t2_line3, rtol=.01)
+            
+            # Test all time steps
+            if multi:
+                refnc = NetCDF(referent)
+                ntimes = t2.shape[0]
+                
+                for t in range(ntimes):
+                    t2 = getvar(wrfin, "T2", timeidx=t)
+                    
+                    line = interpline(t2,wrfin=wrfin,timeidx=t,
+                                  start_point=start_point,
+                                  end_point=end_point,latlon=True)
+                    
+                    refname = "t2_line_t{}".format(t)
+                    refline = refnc.variables[refname][:]
+                    
+                    nt.assert_allclose(to_np(line), 
+                               to_np(refline),rtol=.005)
+            
+                refnc.close()
+            
+            # Test NCLs single time case
+            if not multi:
+                refnc = NetCDF(referent)
+                ref_t2_line4 = refnc.variables["t2_line4"][:]
+                
+                t2 = getvar(wrfin, "T2", timeidx=0)
+                line = interpline(t2,wrfin=wrfin,timeidx=0,
+                                  start_point=start_point,
+                                  end_point=end_point,latlon=True)
+                
+                nt.assert_allclose(to_np(line), 
+                               to_np(ref_t2_line4),rtol=.005)
+                refnc.close()
+
         elif (varname == "vinterp"):
             # Tk to theta
-            fld_tk_theta = _get_refvals(referent, "fld_tk_theta", repeat, multi)
+            fld_tk_theta = _get_refvals(referent, "fld_tk_theta", multi)
             fld_tk_theta = np.squeeze(fld_tk_theta)
             
-            tk = getvar(in_wrfnc, "temp", timeidx=timeidx, units="k")
+            tk = getvar(wrfin, "temp", timeidx=timeidx, units="k")
             
             interp_levels = [200,300,500,1000]
             
             # Make sure the numpy version works
-            field = vinterp(in_wrfnc, 
+            field = vinterp(wrfin, 
                             field=to_np(tk), 
                             vert_coord="theta", 
                             interp_levels=interp_levels, 
@@ -549,7 +610,7 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
                             timeidx=timeidx, 
                             log_p=True)
             
-            field = vinterp(in_wrfnc, 
+            field = vinterp(wrfin, 
                             field=tk, 
                             vert_coord="theta", 
                             interp_levels=interp_levels, 
@@ -566,12 +627,12 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             nt.assert_allclose(to_np(field), fld_tk_theta, tol, atol)
             
             # Tk to theta-e
-            fld_tk_theta_e = _get_refvals(referent, "fld_tk_theta_e", repeat, multi)
+            fld_tk_theta_e = _get_refvals(referent, "fld_tk_theta_e", multi)
             fld_tk_theta_e = np.squeeze(fld_tk_theta_e)
             
             interp_levels = [200,300,500,1000]
             
-            field = vinterp(in_wrfnc, 
+            field = vinterp(wrfin, 
                             field=tk, 
                             vert_coord="theta-e", 
                             interp_levels=interp_levels, 
@@ -588,12 +649,12 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             nt.assert_allclose(to_np(field), fld_tk_theta_e, tol, atol)
             
             # Tk to pressure
-            fld_tk_pres = _get_refvals(referent, "fld_tk_pres", repeat, multi)
+            fld_tk_pres = _get_refvals(referent, "fld_tk_pres", multi)
             fld_tk_pres = np.squeeze(fld_tk_pres)
             
             interp_levels = [850,500]
             
-            field = vinterp(in_wrfnc, 
+            field = vinterp(wrfin, 
                             field=tk, 
                             vert_coord="pressure", 
                             interp_levels=interp_levels, 
@@ -608,11 +669,11 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             nt.assert_allclose(to_np(field), fld_tk_pres, tol, atol)
             
             # Tk to geoht_msl
-            fld_tk_ght_msl = _get_refvals(referent, "fld_tk_ght_msl", repeat, multi)
+            fld_tk_ght_msl = _get_refvals(referent, "fld_tk_ght_msl", multi)
             fld_tk_ght_msl = np.squeeze(fld_tk_ght_msl)
             interp_levels = [1,2]
             
-            field = vinterp(in_wrfnc, 
+            field = vinterp(wrfin, 
                             field=tk, 
                             vert_coord="ght_msl", 
                             interp_levels=interp_levels, 
@@ -626,11 +687,11 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             nt.assert_allclose(to_np(field), fld_tk_ght_msl, tol, atol)
             
             # Tk to geoht_agl
-            fld_tk_ght_agl = _get_refvals(referent, "fld_tk_ght_agl", repeat, multi)
+            fld_tk_ght_agl = _get_refvals(referent, "fld_tk_ght_agl", multi)
             fld_tk_ght_agl = np.squeeze(fld_tk_ght_agl)
             interp_levels = [1,2]
             
-            field = vinterp(in_wrfnc, 
+            field = vinterp(wrfin, 
                             field=tk, 
                             vert_coord="ght_agl", 
                             interp_levels=interp_levels, 
@@ -644,12 +705,12 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             nt.assert_allclose(to_np(field), fld_tk_ght_agl, tol, atol)
             
             # Hgt to pressure
-            fld_ht_pres = _get_refvals(referent, "fld_ht_pres", repeat, multi)
+            fld_ht_pres = _get_refvals(referent, "fld_ht_pres", multi)
             fld_ht_pres = np.squeeze(fld_ht_pres)
             
-            z = getvar(in_wrfnc, "height", timeidx=timeidx, units="m")
+            z = getvar(wrfin, "height", timeidx=timeidx, units="m")
             interp_levels = [500,50]
-            field = vinterp(in_wrfnc, 
+            field = vinterp(wrfin, 
                             field=z, 
                             vert_coord="pressure", 
                             interp_levels=interp_levels, 
@@ -663,12 +724,12 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             nt.assert_allclose(to_np(field), fld_ht_pres, tol, atol)
             
             # Pressure to theta
-            fld_pres_theta = _get_refvals(referent, "fld_pres_theta", repeat, multi)
+            fld_pres_theta = _get_refvals(referent, "fld_pres_theta", multi)
             fld_pres_theta = np.squeeze(fld_pres_theta)
             
-            p = getvar(in_wrfnc, "pressure", timeidx=timeidx)
+            p = getvar(wrfin, "pressure", timeidx=timeidx)
             interp_levels = [200,300,500,1000]
-            field = vinterp(in_wrfnc, 
+            field = vinterp(wrfin, 
                             field=p, 
                             vert_coord="theta", 
                             interp_levels=interp_levels, 
@@ -682,12 +743,12 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
             nt.assert_allclose(to_np(field), fld_pres_theta, tol, atol)
             
             # Theta-e to pres
-            fld_thetae_pres = _get_refvals(referent, "fld_thetae_pres", repeat, multi)
+            fld_thetae_pres = _get_refvals(referent, "fld_thetae_pres", multi)
             fld_thetae_pres = np.squeeze(fld_thetae_pres)
             
-            eth = getvar(in_wrfnc, "eth", timeidx=timeidx)
+            eth = getvar(wrfin, "eth", timeidx=timeidx)
             interp_levels = [850,500,5]
-            field = vinterp(in_wrfnc, 
+            field = vinterp(wrfin, 
                             field=eth, 
                             vert_coord="pressure", 
                             interp_levels=interp_levels, 
@@ -702,25 +763,31 @@ def make_interp_test(varname, wrf_in, referent, multi=False,
     
     return test
 
-def extract_proj_params(wrfnc):
+def extract_proj_params(wrfnc, timeidx=0):
     attrs = extract_global_attrs(wrfnc, ("MAP_PROJ", "TRUELAT1", "TRUELAT2",
                                          "STAND_LON", "POLE_LAT", "POLE_LON", 
                                          "DX", "DY"))
     
     result = {key.lower(): val for key,val in viewitems(attrs)}
     
+    _timeidx = timeidx
     if is_multi_file(wrfnc):
-        wrfnc = wrfnc[0]
+        wrfnc0 = wrfnc[0]
+        num_times_per_file = len(wrfnc0.dimensions["Time"])
+        file_idx = timeidx // num_times_per_file
+        _timeidx = timeidx % num_times_per_file
+        
+        wrfnc = wrfnc[file_idx]
     
     result["known_x"] = 0
     result["known_y"] = 0
-    result["ref_lat"] = wrfnc.variables["XLAT"][0,0,0]
-    result["ref_lon"] = wrfnc.variables["XLONG"][0,0,0]
+    result["ref_lat"] = wrfnc.variables["XLAT"][_timeidx,0,0]
+    result["ref_lon"] = wrfnc.variables["XLONG"][_timeidx,0,0]
     
     return result
 
-def make_latlon_test(testid, wrf_in, referent, single, multi=False, repeat=3, 
-                     pynio=False):
+def make_latlon_test(testid, dir, pattern, referent, single, 
+                     multi=False, pynio=False):
     def test(self):
         try:
             from netCDF4 import Dataset as NetCDF
@@ -732,113 +799,135 @@ def make_latlon_test(testid, wrf_in, referent, single, multi=False, repeat=3,
         except:
             pass
         
-        if not multi:
-            timeidx = 0
-            if not pynio:
-                in_wrfnc = NetCDF(wrf_in)
-                try:
-                    in_wrfnc.set_auto_mask(False)
-                except:
-                    pass
-            else:
-                # Note: Python doesn't like it if you reassign an outer scoped
-                # variable (wrf_in in this case)
-                if not wrf_in.endswith(".nc"):
-                    wrf_file = wrf_in + ".nc"
-                else:
-                    wrf_file = wrf_in
-                in_wrfnc = Nio.open_file(wrf_file)
-        else:
-            timeidx = None
+        timeidx = 0 if not multi else None
+        pat = os.path.join(dir, pattern)
+        wrf_files = glob.glob(pat)
+        wrf_files.sort()
+        
+        refnc = NetCDF(referent)
+        try:
+            refnc.set_always_mask(False)
+        except:
+            pass
+        
+        wrfin = []
+        for fname in wrf_files:
             if not pynio:
-                nc = NetCDF(wrf_in)
+                f = NetCDF(fname)
                 try:
-                    nc.set_auto_mask(False)
+                    f.set_auto_mask(False)
                 except:
                     pass
-                in_wrfnc = [nc for i in xrange(repeat)]
+                wrfin.append(f)
             else:
-                if not wrf_in.endswith(".nc"):
-                    wrf_file = wrf_in + ".nc"
+                if not fname.endswith(".nc"):
+                    _fname = fname + ".nc"
                 else:
-                    wrf_file = wrf_in
-                nc = Nio.open_file(wrf_file)
-                in_wrfnc = [nc for i in xrange(repeat)]
-                
-        refnc = NetCDF(referent)
-        try:
-            refnc.set_auto_mask(False)
-        except:
-            pass
+                    _fname = fname
+                f = Nio.open_file(_fname)
+                wrfin.append(f)
                 
         if testid == "xy":
-            # Since this domain is not moving, the reference values are the 
-            # same whether there are multiple or single files
-            ref_vals = refnc.variables["ij"][:]
+            
             # Lats/Lons taken from NCL script, just hard-coding for now
-            lats = [-55, -60, -65]
-            lons = [25, 30, 35]
+            lats = [22.0, 25.0, 27.0]
+            lons = [-90.0, -87.5, -83.75]
             
             # Just call with a single lat/lon
             if single:
-                xy = ll_to_xy(in_wrfnc, lats[0], lons[0])
-                xy = xy + 1 # NCL uses fortran indexing
+                timeidx = 8
+                ref_vals = refnc.variables["xy2"][:]
+                
+                xy = ll_to_xy(wrfin, lats[0], lons[0], timeidx=timeidx,
+                              as_int=True)
                 ref = ref_vals[:,0]
                 
                 nt.assert_allclose(to_np(xy), ref)
                 
                 # Next make sure the 'proj' version works
-                projparams = extract_proj_params(in_wrfnc)
-                xy_proj = ll_to_xy_proj(lats[0], lons[0], **projparams)
+                projparams = extract_proj_params(wrfin, timeidx=timeidx)
+                xy_proj = ll_to_xy_proj(lats[0], lons[0],
+                                        as_int=True, 
+                                        **projparams)
                 
-                nt.assert_allclose(to_np(xy_proj), to_np(xy-1))
+                nt.assert_allclose(to_np(xy_proj), to_np(xy))
                 
             
             else:
-                xy = ll_to_xy(in_wrfnc, lats, lons)
-                xy = xy + 1 # NCL uses fortran indexing
+                ref_vals = refnc.variables["xy1"][:]
+                xy = ll_to_xy(wrfin, lats, lons, timeidx=None, as_int=False)
+                
                 ref = ref_vals[:]
                 
                 nt.assert_allclose(to_np(xy), ref)
                 
-                # Next make sure the 'proj' version works
-                projparams = extract_proj_params(in_wrfnc)
-                xy_proj = ll_to_xy_proj(lats, lons, **projparams)
+                if xy.ndim > 2:
+                    # Moving nest
+                    is_moving = True
+                    numtimes = xy.shape[-2]
+                else:
+                    is_moving = False
+                    numtimes = 1
+                
+                for tidx in range(9):
                 
-                nt.assert_allclose(to_np(xy_proj), to_np(xy-1))
+                    # Next make sure the 'proj' version works
+                    projparams = extract_proj_params(wrfin, timeidx=tidx)
+                    xy_proj = ll_to_xy_proj(lats, lons, as_int=False,
+                                            **projparams)
+                    
+                    if is_moving:
+                        idxs = (slice(None), tidx, slice(None))
+                    else:
+                        idxs = (slice(None),)
+                
+                    nt.assert_allclose(to_np(xy_proj), to_np(xy[idxs]))
                 
         else:
-            # Since this domain is not moving, the reference values are the 
-            # same whether there are multiple or single files
-            ref_vals = refnc.variables["ll"][:]
-            
              # i_s, j_s taken from NCL script, just hard-coding for now
              # NCL uses 1-based indexing for this, so need to subtract 1
-            i_s = np.asarray([10, 100, 150], int) - 1
-            j_s = np.asarray([10, 100, 150], int) - 1
+            x_s = np.asarray([10, 50, 90], int)
+            y_s = np.asarray([10, 50, 90], int)
             
             if single:
-                ll = xy_to_ll(in_wrfnc, i_s[0], j_s[0])
+                timeidx=8
+                ref_vals = refnc.variables["ll2"][:]
+                ll = xy_to_ll(wrfin, x_s[0], y_s[0], timeidx=timeidx)
                 ref = ref_vals[::-1,0]
                 
                 nt.assert_allclose(to_np(ll), ref)
                 
                 # Next make sure the 'proj' version works
-                projparams = extract_proj_params(in_wrfnc)
-                ll_proj = xy_to_ll_proj(i_s[0], j_s[0], **projparams)
+                projparams = extract_proj_params(wrfin, timeidx=8)
+                ll_proj = xy_to_ll_proj(x_s[0], y_s[0], **projparams)
                 
                 nt.assert_allclose(to_np(ll_proj), to_np(ll))
             else:
-                ll = xy_to_ll(in_wrfnc, i_s, j_s)
+                ref_vals = refnc.variables["ll1"][:]
+                ll = xy_to_ll(wrfin, x_s, y_s, timeidx=None)
                 ref = ref_vals[::-1,:]
                 
                 nt.assert_allclose(to_np(ll), ref)
                 
-                # Next make sure the 'proj' version works
-                projparams = extract_proj_params(in_wrfnc)
-                ll_proj = xy_to_ll_proj(i_s, j_s, **projparams)
+                if ll.ndim > 2:
+                    # Moving nest
+                    is_moving = True
+                    numtimes = ll.shape[-2]
+                else:
+                    is_moving = False
+                    numtimes = 1
                 
-                nt.assert_allclose(to_np(ll_proj), to_np(ll))
+                for tidx in range(numtimes):
+                    # Next make sure the 'proj' version works
+                    projparams = extract_proj_params(wrfin, timeidx=tidx)
+                    ll_proj = xy_to_ll_proj(x_s, y_s, **projparams)
+                    
+                    if is_moving:
+                        idxs = (slice(None), tidx, slice(None))
+                    else:
+                        idxs = (slice(None),)
+                        
+                    nt.assert_allclose(to_np(ll_proj), to_np(ll[idxs]))
         
     return test
 
@@ -869,82 +958,85 @@ if __name__ == "__main__":
     interp_methods = ["interplevel", "vertcross", "interpline", "vinterp"]
     latlon_tests = ["xy", "ll"]
     
-    try:
-        import netCDF4
-    except ImportError:
-        pass
-    else:
-        for var in wrf_vars:
-            if var in ignore_vars:
-                continue
-            
-            test_func1 = make_test(var, TEST_FILE, OUT_NC_FILE)
-            test_func2 = make_test(var, TEST_FILE, OUT_NC_FILE, multi=True)
-            setattr(WRFVarsTest, 'test_{0}'.format(var), test_func1)
-            setattr(WRFVarsTest, 'test_multi_{0}'.format(var), test_func2)
-            
-        for method in interp_methods:
-            test_interp_func1 = make_interp_test(method, TEST_FILE, 
-                                                 OUT_NC_FILE)
-            test_interp_func2 = make_interp_test(method, TEST_FILE, 
-                                                 OUT_NC_FILE, multi=True)
-            setattr(WRFInterpTest, 'test_{0}'.format(method), 
-                    test_interp_func1)
-            setattr(WRFInterpTest, 'test_multi_{0}'.format(method), 
-                    test_interp_func2)
-        
-        for testid in latlon_tests:
-            for single in (True, False):
-                for multi in (True, False):
-                    test_ll_func = make_latlon_test(testid, TEST_FILE, 
-                                                    OUT_NC_FILE, 
-                                                    single=single, multi=multi, 
-                                                    repeat=3, pynio=False)
-                    multistr = "" if not multi else "_multi"
-                    singlestr = "_nosingle" if not single else "_single"
-                    test_name = "test_{}{}{}".format(testid, singlestr, 
-                                                       multistr)
-                    setattr(WRFLatLonTest, test_name, test_ll_func)
+    for dir, ref_nc_file, nest in zip(DIRS, REF_NC_FILES, NEST):
+        try:
+            import netCDF4
+        except ImportError:
+            pass
+        else:
+            for var in wrf_vars:
+                if var in ignore_vars:
+                    continue
                 
-    try:
-        import PyNIO
-    except ImportError:
-        pass
-    else:
-        for var in wrf_vars:
-            if var in ignore_vars:
-                continue
-            
-            test_func1 = make_test(var, TEST_FILE, OUT_NC_FILE, pynio=True)
-            test_func2 = make_test(var, TEST_FILE, OUT_NC_FILE, multi=True,
-                                   pynio=True)
-            setattr(WRFVarsTest, 'test_pynio_{0}'.format(var), test_func1)
-            setattr(WRFVarsTest, 'test_pynio_multi_{0}'.format(var), 
-                    test_func2)
-            
-        for method in interp_methods:
-            test_interp_func1 = make_interp_test(method, TEST_FILE, 
-                                                 OUT_NC_FILE)
-            test_interp_func2 = make_interp_test(method, TEST_FILE, 
-                                                 OUT_NC_FILE, multi=True)
-            setattr(WRFInterpTest, 'test_pynio_{0}'.format(method), 
-                    test_interp_func1)
-            setattr(WRFInterpTest, 'test_pynio_multi_{0}'.format(method), 
-                    test_interp_func2)
-            
-        for testid in latlon_tests:
-            for single in (True, False):
-                for multi in (True, False):
-                    test_ll_func = make_latlon_test(testid, TEST_FILE, 
-                                                    OUT_NC_FILE, 
-                                                    single=single, multi=multi, 
-                                                    repeat=3, pynio=False)
-                    multistr = "" if not multi else "_multi"
-                    singlestr = "_nosingle" if not single else "_single"
-                    test_name = "test_pynio_{}{}{}".format(testid, 
-                                                              singlestr, 
-                                                              multistr)
-                    setattr(WRFLatLonTest, test_name, test_ll_func)
+                test_func1 = make_test(var, dir, PATTERN, ref_nc_file)
+                test_func2 = make_test(var, dir, PATTERN, ref_nc_file, multi=True)
+                setattr(WRFVarsTest, 'test_{0}_{1}'.format(nest,var), test_func1)
+                setattr(WRFVarsTest, 'test_{0}_multi_{1}'.format(nest,var), test_func2)
+                
+            for method in interp_methods:
+                test_interp_func1 = make_interp_test(method, dir, PATTERN, 
+                                                     ref_nc_file)
+                test_interp_func2 = make_interp_test(method, dir, PATTERN, 
+                                                     ref_nc_file, multi=True)
+                setattr(WRFInterpTest, 'test_{0}_{1}'.format(nest,method), 
+                        test_interp_func1)
+                setattr(WRFInterpTest, 'test_{0}_multi_{1}'.format(nest,method), 
+                        test_interp_func2)
+            
+            for testid in latlon_tests:
+                for single in (True, False):
+                    for multi in (True, False):
+                        test_ll_func = make_latlon_test(testid, dir, PATTERN, 
+                                                        ref_nc_file, 
+                                                        single=single, 
+                                                        multi=multi, 
+                                                        pynio=False)
+                        multistr = "" if not multi else "_multi"
+                        singlestr = "_nosingle" if not single else "_single"
+                        test_name = "test_{}_{}{}{}".format(nest, testid, singlestr, 
+                                                           multistr)
+                        setattr(WRFLatLonTest, test_name, test_ll_func)
+                    
+        try:
+            import PyNIO
+        except ImportError:
+            pass
+        else:
+            for var in wrf_vars:
+                if var in ignore_vars:
+                    continue
+                
+                test_func1 = make_test(var, dir, PATTERN, ref_nc_file, pynio=True)
+                test_func2 = make_test(var, dir, PATTERN, ref_nc_file, multi=True,
+                                       pynio=True)
+                setattr(WRFVarsTest, 'test_pynio_{0}_{1}'.format(nest,var), test_func1)
+                setattr(WRFVarsTest, 'test_pynio_{0}_multi_{1}'.format(nest,var), 
+                        test_func2)
+                
+            for method in interp_methods:
+                test_interp_func1 = make_interp_test(method, dir, PATTERN, 
+                                                     ref_nc_file)
+                test_interp_func2 = make_interp_test(method, dir, PATTERN, 
+                                                     ref_nc_file, multi=True)
+                setattr(WRFInterpTest, 'test_pynio_{0}_{1}'.format(nest,method), 
+                        test_interp_func1)
+                setattr(WRFInterpTest, 'test_pynio_{0}_multi_{1}'.format(nest,method), 
+                        test_interp_func2)
+                
+            for testid in latlon_tests:
+                for single in (True, False):
+                    for multi in (True, False):
+                        test_ll_func = make_latlon_test(testid, dir, PATTERN, 
+                                                        ref_nc_file, 
+                                                        single=single, 
+                                                        multi=multi, 
+                                                        pynio=False)
+                        multistr = "" if not multi else "_multi"
+                        singlestr = "_nosingle" if not single else "_single"
+                        test_name = "test_pynio_{}_{}{}{}".format(nest, testid, 
+                                                                  singlestr, 
+                                                                  multistr)
+                        setattr(WRFLatLonTest, test_name, test_ll_func)
      
     ut.main()
     
\ No newline at end of file