Merge branch 'feature/test_update' into develop

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,14 +529,14 @@ 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:
@@ -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 

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")
     

src/wrf/latlonutils.py

@@ -167,7 +167,7 @@ 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):
+        # Always use the 0th time for non-moving domains to avoid problems
         lat_timeidx = 0
             
         if is_multi_file(wrfin):

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:
@@ -1150,9 +1187,46 @@ def _set_line_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 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

src/wrf/util.py

@@ -3881,11 +3881,22 @@ 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
     
     
     

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,

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)
                 

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,24 +97,28 @@
     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 /)
     
@@ -113,12 +132,40 @@
     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)