diff --git a/CHANGES.rst b/CHANGES.rst
index e8e3b3eeca..1306a4a6b7 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -58,8 +58,7 @@ RELEASE FILE IN doc/releases
       files that have frametype None (this prevent ``run_pypeit`` to crash)
     - Added a function ``check_spectrograph()`` (currently only defined for LRIS),
       that checks (during ``pypeit_setup``) if the selected spectrograph is the
-      corrected one for the data used.
-
+      corrected one for the data used. 
 
 1.13.0 (2 June 2023)
 --------------------
diff --git a/doc/releases/1.14.1dev.rst b/doc/releases/1.14.1dev.rst
index 9023193489..0133908452 100644
--- a/doc/releases/1.14.1dev.rst
+++ b/doc/releases/1.14.1dev.rst
@@ -43,6 +43,7 @@ Functionality/Performance Improvements and Additions
   to reflect the *type* of polynomial being fit.  Both names point to the same
   code, but the name ``'polynomial'`` is deprecated and will be removed in the
   future.
+- Introduced PCA method for telluric corrections
 - Added a new GUI for creating and editing PypeIt input files: ``pypeit_setup_gui``
 
 Instrument-specific Updates
diff --git a/pypeit/core/telluric.py b/pypeit/core/telluric.py
index 378e57ce81..f3bbc96147 100644
--- a/pypeit/core/telluric.py
+++ b/pypeit/core/telluric.py
@@ -39,7 +39,7 @@
 # TODO These codes should probably be in a separate qso_pca module. Also pickle functionality needs to be removed.
 # The prior is not used (that was the reason for pickling), so the components could be stored in fits format.
 # npca is not actually required here.
-def init_pca(filename,wave_grid,redshift, npca):
+def qso_init_pca(filename,wave_grid,redshift,npca):
     """
     This routine reads in the pickle file created by  coarse_pca.create_coarse_pca.
     The relevant pieces are the wavelengths (wave_pca_c), the PCA components (pca_comp_c), and the Gaussian mixture
@@ -72,19 +72,18 @@ def init_pca(filename,wave_grid,redshift, npca):
     pca_table = table.Table.read(file_with_path)
     wave_pca_c = pca_table['WAVE_PCA'][0].flatten()
     pca_comp_c = pca_table['PCA_COMP'][0][0,:,:]
-    coeffs_c =pca_table['PCA_COEFFS'][0][0,:,:]
-    #wave_pca_c, cont_all_c, pca_comp_c, coeffs_c, mean_pca, covar_pca, diff_pca, mix_fit, chi2, dof = pickle.load(open(filename,'rb'))
+    coeffs_c = pca_table['PCA_COEFFS'][0][0,:,:]
     num_comp = pca_comp_c.shape[0] # number of PCA components
     # Interpolate PCA components onto wave_grid
-    pca_interp = scipy.interpolate.interp1d(wave_pca_c*(1.0 + redshift),pca_comp_c, bounds_error=False, fill_value=0.0, axis=1)
+    pca_interp = scipy.interpolate.interp1d(wave_pca_c*(1.0+redshift), pca_comp_c, bounds_error=False, fill_value=0.0, axis=1)
     pca_comp_new = pca_interp(wave_grid)
     # Generate a mixture model for the coefficients prior, what should ngauss be?
     #prior = mixture.GaussianMixture(n_components = npca-1).fit(coeffs_c[:, 1:npca])
     # Construct the PCA dict
-    pca_dict = {'npca': npca, 'components': pca_comp_new, 'coeffs': coeffs_c, 'z_fid': redshift, 'dloglam': dloglam}
-    return pca_dict
+    qso_pca_dict = {'npca': npca, 'components': pca_comp_new, 'coeffs': coeffs_c, 'z_fid': redshift, 'dloglam': dloglam}
+    return qso_pca_dict
 
-def pca_eval(theta,pca_dict):
+def qso_pca_eval(theta,qso_pca_dict):
     """
     Function for evaluating the quasar PCA
 
@@ -92,18 +91,18 @@ def pca_eval(theta,pca_dict):
         theta (`numpy.ndarray`_):
           Parameter vector, where theta_pca[0] is redshift, theta_pca[1] is the normalization, and
           theta_pca[2:npca+1] are the PCA coefficients, where npca is the PCA dimensionality
-        pca_dict (dict):
-          Dictionary continaing the PCA information generated by init_pca
+        qso_pca_dict (dict):
+          Dictionary continaing the PCA information generated by qso_init_pca
 
     Returns:
         `numpy.ndarray`_: Evaluated PCA for the QSO
 
     """
 
-    C = pca_dict['components']
-    z_fid = pca_dict['z_fid']
-    dloglam = pca_dict['dloglam']
-    npca = pca_dict['npca']  # Size of the PCA currently being used, original PCA in the dict could be larger
+    C = qso_pca_dict['components']
+    z_fid = qso_pca_dict['z_fid']
+    dloglam = qso_pca_dict['dloglam']
+    npca = qso_pca_dict['npca']  # Size of the PCA currently being used, original PCA in the dict could be larger
     z_qso = theta[0]
     norm = theta[1]
     A = theta[2:]
@@ -112,7 +111,7 @@ def pca_eval(theta,pca_dict):
     return norm*np.exp(np.dot(np.append(1.0,A),C_now))
 
 # TODO The prior is not currently used, but this is left in here anyway.
-#def pca_lnprior(theta,pca_dict):
+#def qso_pca_lnprior(theta,qso_pca_dict):
 #    """
 #    Routine to evaluate the the ln of the prior probability lnPrior,
 #    from the Gaussian mixture model fit to the distriution of PCA
@@ -124,9 +123,9 @@ def pca_eval(theta,pca_dict):
 #            ``theta_pca[1]`` is the normalization, and
 #            ``theta_pca[2:npca+1]`` are the PCA coefficients, where npca
 #            is the PCA dimensionality
-#        pca_dict (dict):
+#        qso_pca_dict (dict):
 #            Dictionary continaing the PCA information generated by
-#            ``init_pca``
+#            ``qso_init_pca``
 #
 #
 #    Returns:
@@ -135,30 +134,82 @@ def pca_eval(theta,pca_dict):
 #        theta_pca[2:npca+1]``
 #
 #    """
-#    gaussian_mixture_model = pca_dict['prior']
+#    gaussian_mixture_model = qso_pca_dict['prior']
 #    A = theta[2:]
 #    return gaussian_mixture_model.score_samples(A.reshape(1,-1))
 
 
+def read_telluric_pca(filename, wave_min=None, wave_max=None, pad_frac=0.10):
+    """
+    Reads in the telluric PCA components from a file.
+
+    Optionally, this method also trims wavelength to be in within ``wave_min``
+    and ``wave_max`` and pads the data (see ``pad_frac``).
+
+    Args:
+        filename (:obj:`str`):
+            Telluric PCA filename
+        wave_min (:obj:`float`, optional):
+            Minimum wavelength at which the grid is desired
+        wave_max (:obj:`float`, optional):
+            Maximum wavelength at which the grid is desired.
+        pad_frac (:obj:`float`, optional):
+            Percentage padding to be added to the grid boundaries if
+            ``wave_min`` or ``wave_max`` are input; ignored otherwise. The
+            resulting grid will extend from ``(1.0 - pad_frac)*wave_min`` to
+            ``(1.0 + pad_frac)*wave_max``.
+
+    Returns:
+        :obj:`dict`: Dictionary containing the telluric PCA components.
+            - wave_grid: Telluric model wavelength grid
+            - dloglam: Wavelength sampling of telluric model
+            - tell_pad_pix: Number of pixels to pad at edges for convolution
+            - ncomp_tell_pca: Number of PCA components
+            - tell_pca: PCA component vectors
+            - bounds_tell_pca: Maximum/minimum coefficient
+            - coefs_tell_pca: Set of model coefficient values (for prior in future)
+            - teltype: Type of telluric model, i.e. 'pca'
+    """
+    # load_telluric_grid() takes care of path and existance check
+    hdul = data.load_telluric_grid(filename)
+    wave_grid_full = hdul[1].data
+    pca_comp_full = hdul[0].data
+    nspec_full = wave_grid_full.size
+    ncomp = hdul[0].header['NCOMP']
+    bounds = hdul[2].data
+    model_coefs = hdul[3].data
+
+    ind_lower = np.argmin(np.abs(wave_grid_full - (1.0 - pad_frac)*wave_min)) \
+                    if wave_min is not None else 0
+    ind_upper = np.argmin(np.abs(wave_grid_full - (1.0 + pad_frac)*wave_max)) \
+                    if wave_max is not None else nspec_full
+    wave_grid = wave_grid_full[ind_lower:ind_upper]
+    pca_comp_grid = pca_comp_full[:,ind_lower:ind_upper]
+
+    dwave, dloglam, resln_guess, pix_per_sigma = wvutils.get_sampling(wave_grid)
+    tell_pad_pix = int(np.ceil(10.0 * pix_per_sigma))
+
+    return dict(wave_grid=wave_grid, dloglam=dloglam,
+                tell_pad_pix=tell_pad_pix, ncomp_tell_pca=ncomp,
+                tell_pca=pca_comp_grid, bounds_tell_pca=bounds,
+                coefs_tell_pca=model_coefs, teltype='pca')
+            
 def read_telluric_grid(filename, wave_min=None, wave_max=None, pad_frac=0.10):
     """
-    Reads in the telluric grid from a file.
+    Reads in the telluric grid from a file. This method is no longer the
+    preferred approach; see "read_telluric_pca" for the PCA mode.
 
     Optionally, this method also trims the grid to be in within ``wave_min``
     and ``wave_max`` and pads the data (see ``pad_frac``).
 
-    .. todo::
-        List and describe the contents of the dictionary in the return
-        description.
-
     Args:
         filename (:obj:`str`):
            Telluric grid filename
-        wave_min (:obj:`float`):
+        wave_min (:obj:`float`, optional):
            Minimum wavelength at which the grid is desired
-        wave_max (:obj:`float`):
+        wave_max (:obj:`float`, optional):
            Maximum wavelength at which the grid is desired.
-        pad_frac (:obj:`float`):
+        pad_frac (:obj:`float`, optional):
            Percentage padding to be added to the grid boundaries if
            ``wave_min`` or ``wave_max`` are input; ignored otherwise. The
            resulting grid will extend from ``(1.0 - pad_frac)*wave_min`` to
@@ -166,16 +217,15 @@ def read_telluric_grid(filename, wave_min=None, wave_max=None, pad_frac=0.10):
 
     Returns:
         :obj:`dict`:  Dictionary containing the telluric grid
-            - wave_grid=
-            - dloglam=
-            - resln_guess=
-            - pix_per_sigma=
-            - tell_pad_pix=
-            - pressure_grid=
-            - temp_grid=
-            - h2o_grid=
-            - airmass_grid=
-            - tell_grid=
+            - wave_grid: Telluric model wavelength grid
+            - dloglam: Wavelength sampling of telluric model
+            - tell_pad_pix: Number of pixels to pad at edges for convolution
+            - pressure_grid: Atmospheric pressure values in telluric grid [mbar]
+            - temp_grid: Temperature values in telluric grid [degrees C]
+            - h2o_grid: Humidity values in telluric grid [%]
+            - airmass_grid: Airmass values in telluric grid
+            - tell_grid: Grid of telluric models
+            - teltype: Type of telluric model, i.e. 'grid'
     """
     # load_telluric_grid() takes care of path and existance check
     hdul = data.load_telluric_grid(filename)
@@ -201,17 +251,9 @@ def read_telluric_grid(filename, wave_min=None, wave_max=None, pad_frac=0.10):
     dwave, dloglam, resln_guess, pix_per_sigma = wvutils.get_sampling(wave_grid)
     tell_pad_pix = int(np.ceil(10.0 * pix_per_sigma))
 
-
-    return dict(wave_grid=wave_grid,
-                dloglam=dloglam, 
-                resln_guess=resln_guess,
-                pix_per_sigma=pix_per_sigma, 
-                tell_pad_pix=tell_pad_pix, 
-                pressure_grid=pg,
-                temp_grid=tg, 
-                h2o_grid=hg, 
-                airmass_grid=ag, 
-                tell_grid=model_grid)
+    return dict(wave_grid=wave_grid, dloglam=dloglam, tell_pad_pix=tell_pad_pix,
+                pressure_grid=pg, temp_grid=tg, h2o_grid=hg, airmass_grid=ag,
+                tell_grid=model_grid, teltype='grid')
 
 
 def interp_telluric_grid(theta, tell_dict):
@@ -259,8 +301,8 @@ def conv_telluric(tell_model, dloglam, res):
     Args:
         tell_model (`numpy.ndarray`_):
             Input telluric model at the native resolution of the telluric model grid. The shape of this input is in
-            general  different from the size of the telluric grid (read in by read_telluric_grid above) because it is
-            trimmed to relevant wavelenghts using ind_lower, ind_upper. See eval_telluric below.
+            general different from the size of the raw telluric model (read in by read_telluric_* above) because it is
+            trimmed to relevant wavelengths using ind_lower, ind_upper. See eval_telluric below.
         dloglam (float):
             Wavelength spacing of the telluric grid expressed as a dlog10(lambda), i.e. stored in the
             tell_dict as tell_dict['dloglam']
@@ -294,12 +336,13 @@ def conv_telluric(tell_model, dloglam, res):
 def shift_telluric(tell_model, loglam, dloglam, shift, stretch):
     """
     Routine to apply a shift to the telluric model. Note that the shift can be sub-pixel, i.e this routine interpolates.
+    Note that the shift units are pixels *of the telluric model*.
 
     Args:
         tell_model (`numpy.ndarray`_):
-            Input telluric model. The shape of this input is in general  different from the size of the telluric grid
-            (read in by read_telluric_grid above) because it is trimmed to relevant wavelenghts using ind_lower, ind_upper.
-            See eval_telluric below.
+            Input telluric model. The shape of this input is in general different from the size of the telluric models
+            (read in by read_telluric_* above) because it is trimmed to relevant wavelengths using ind_lower, ind_upper.
+            See eval_telluric_* below.
 
         loglam (`numpy.ndarray`_):
             The log10 of the wavelength grid on which the tell_model is evaluated.
@@ -320,43 +363,43 @@ def shift_telluric(tell_model, loglam, dloglam, shift, stretch):
     tell_model_shift = np.interp(loglam_shift, loglam, tell_model)
     return tell_model_shift
 
-
 def eval_telluric(theta_tell, tell_dict, ind_lower=None, ind_upper=None):
     """
-    Evaluate the telluric model at an arbitrary location in parameter space.
-    
-    The full atmosphere model parameter space is either 5 or 7 dimensional,
-    which is the size of the ``theta_tell`` input parameter vector.
+    Evaluate the telluric model.
 
-    The parameters provided by ``theta_tell`` must be: pressure, temperature,
-    humidity, airmass, spectral resolution, shift, and stretch. The latter two
-    can be omitted if a shift and stretch of the telluric model are not
-    included.
+    Parameters from ``theta_tell`` are: (if teltype == 'pca') the PCA coefficients
+    or (if teltype == 'grid') the telluric grid parameters pressure, temperature,
+    humidity, and airmass, in both cases followed by spectral resolution, shift,
+    and stretch.
 
     This routine performs the following steps:
 
-       1. nearest grid point interpolation of the telluric model onto a
-          new location in the 4-d space (pressure, temperature,
-          humidity, airmass)
+       1. summation of telluric PCA components multiplied by coefficients
+       
+       2. transformation of telluric PCA model from arsinh(tau) to transmission
 
-       2. convolution of the atmosphere model to the resolution set by
+       3. convolution of the atmosphere model to the resolution set by
           the spectral resolution.
 
-       3. (Optional) shift and stretch the telluric model.
+       4. shift and stretch the telluric model.
 
     Args:
         theta_tell (`numpy.ndarray`_):
-            Vector with the telluric model parameters. Must be 5 or 7
-            elements long. See method description.
+            Vector with tell_npca PCA coefficients (if teltype='pca')
+            or pressure, temperature, humidity, and airmass (if teltype='grid'),
+            followed by spectral resolution, shift, and stretch.
+            Final length is then tell_npca+3 or 7.
         tell_dict (:obj:`dict`):
-            Dictionary containing the telluric grid data. See
-            :func:`read_telluric_grid`.
+            Dictionary containing the telluric data. See
+            :func:`read_telluric_pca` if teltype=='pca'.
+            or
+            :func:`read_telluric_grid` if teltype=='grid'.
         ind_lower (:obj:`int`, optional):
             The index of the first pixel to include in the model. Selecting a
             wavelength region for the modeling makes things faster because we
             only need to convolve the portion that is needed for the current
             model fit.
-        ind_upper:
+        ind_upper (:obj:`int`, optional):
             The index (inclusive) of the last pixel to include in the model.
             Selecting a wavelength region for the modeling makes things
             faster because we only need to convolve the portion that is
@@ -364,13 +407,19 @@ def eval_telluric(theta_tell, tell_dict, ind_lower=None, ind_upper=None):
 
     Returns:
         `numpy.ndarray`_: Telluric model evaluated at the desired location
-        theta_tell in model atmosphere parameter space.
+        theta_tell in model atmosphere parameter space. Shape is given by
+        the size of `wave_grid' plus `tell_pad_pix' padding from the input
+        tell_dict.
+        
     """
     ntheta = len(theta_tell)
-    if ntheta not in [5, 7]:
-        msgs.error('Input model atmosphere parameters must have length 5 or 7.')
-
-    tellmodel_hires = interp_telluric_grid(theta_tell[:4], tell_dict)
+    teltype = tell_dict['teltype']
+    # FD: Currently assumes that shift and stretch are on.
+    # TODO: Make this work without shift and stretch.
+    if teltype == 'pca':
+        ntell = ntheta-3
+    elif teltype == 'grid':
+        ntell = 4
 
     # Set the wavelength range if not provided
     ind_lower = 0 if ind_lower is None else ind_lower
@@ -382,16 +431,32 @@ def eval_telluric(theta_tell, tell_dict, ind_lower=None, ind_upper=None):
     ## FW: There is an extreme case with ind_upper == ind_upper_pad, the previous -0 won't work
     ind_lower_final = ind_lower_pad if ind_lower_pad == ind_lower else ind_lower - ind_lower_pad
     ind_upper_final = ind_upper_pad if ind_upper_pad == ind_upper else ind_upper - ind_upper_pad
-    tellmodel_conv = conv_telluric(tellmodel_hires[ind_lower_pad:ind_upper_pad+1],
-                                   tell_dict['dloglam'], theta_tell[4])
-    if ntheta == 5:
-        return tellmodel_conv[ind_lower_final:ind_upper_final]
 
+    if teltype == 'pca':
+        # Evaluate PCA model after truncating the wavelength range
+        tellmodel_hires = np.zeros_like(tell_dict['tell_pca'][0])
+        tellmodel_hires[ind_lower_pad:ind_upper_pad+1] = np.dot(np.append(1,theta_tell[:ntell]),
+                                 tell_dict['tell_pca'][:ntell+1][:,ind_lower_pad:ind_upper_pad+1])
+                                 
+        # PCA model is inverse sinh of the optical depth, convert to transmission here
+        tellmodel_hires[ind_lower_pad:ind_upper_pad+1] = np.sinh(tellmodel_hires[ind_lower_pad:ind_upper_pad+1])
+        # It should generally be very rare, but trim negative optical depths here just in case.
+        clip = tellmodel_hires < 0
+        tellmodel_hires[clip] = 0
+        tellmodel_hires[ind_lower_pad:ind_upper_pad+1] = np.exp(-tellmodel_hires[ind_lower_pad:ind_upper_pad+1])
+    elif teltype == 'grid':
+        # Interpolate within the telluric grid
+        tellmodel_hires = interp_telluric_grid(theta_tell[:ntell], tell_dict)
+
+    tellmodel_conv = conv_telluric(tellmodel_hires[ind_lower_pad:ind_upper_pad+1],
+                                   tell_dict['dloglam'], theta_tell[-3])
+                                   
+    # Stretch and shift telluric wavelength grid
     tellmodel_out = shift_telluric(tellmodel_conv,
                                    np.log10(tell_dict['wave_grid'][ind_lower_pad:ind_upper_pad+1]),
-                                   tell_dict['dloglam'], theta_tell[5], theta_tell[6])
-    return tellmodel_out[ind_lower_final:ind_upper_final]
+                                   tell_dict['dloglam'], theta_tell[-2], theta_tell[-1])
 
+    return tellmodel_out[ind_lower_final:ind_upper_final]
 
 ############################
 #  Fitting routines        #
@@ -405,7 +470,45 @@ def tellfit_chi2(theta, flux, thismask, arg_dict):
 
     Args:
         theta (`numpy.ndarray`_):
-           Parameter vector for the object + telluric model. See documentation of tellfit for a detailed description.
+           
+            Parameter vector for the object + telluric model.
+
+            This is actually two concatenated parameter vectors, one for
+            the object and one for the telluric, i.e.:
+            
+            (in PCA mode)
+                theta_obj = theta[:-(tell_npca+3)]
+                theta_tell = theta[-(tell_npca+3):]
+                
+            (in grid mode)
+                theta_obj = theta[:-7]
+                theta_tell = theta[-7:]
+                
+            The telluric model theta_tell includes a either user-specified
+            number of PCA coefficients (in PCA mode) or ambient pressure,
+            temperature, humidity, and airmass (in grid mode) followed by
+            spectral resolution, shift, and stretch.
+            
+            That is, in PCA mode,
+            
+                pca_coeffs = theta_tell[:tell_npca]
+            
+            while in grid mode,
+            
+                pressure    = theta_tell[0]
+                temperature = theta_tell[1]
+                humidity    = theta_tell[2]
+                airmass     = theta_tell[3]
+                
+            with the last three indices of the array corresponding to
+            
+                resolution = theta_tell[-3]
+                shift      = theta_tell[-2]
+                stretch    = theta_tell[-1]
+
+            The object model theta_obj can have an arbitrary size and is
+            provided as an argument to obj_model_func
+
         flux (`numpy.ndarray`_):
            The flux of the object being fit
         thismask (`numpy.ndarray`_, boolean):
@@ -419,14 +522,23 @@ def tellfit_chi2(theta, flux, thismask, arg_dict):
            that is minimized to perform the fit.
 
     """
-
+    
     obj_model_func = arg_dict['obj_model_func']
     flux_ivar = arg_dict['ivar']
+    teltype = arg_dict['tell_dict']['teltype']
+
+    # TODO: make this work without shift and stretch?
+    # Number of telluric model parameters, plus shift, stretch, and resolution
+    if teltype == 'pca':
+        nfit = arg_dict['tell_npca']+3
+    elif teltype == 'grid':
+        nfit = 4+3
+
+    theta_obj = theta[:-nfit]
+    theta_tell = theta[-nfit:]
 
-    theta_obj = theta[:-7]
-    theta_tell = theta[-7:]
     tell_model = eval_telluric(theta_tell, arg_dict['tell_dict'],
-                               ind_lower=arg_dict['ind_lower'], ind_upper=arg_dict['ind_upper'])
+                                 ind_lower=arg_dict['ind_lower'], ind_upper=arg_dict['ind_upper'])
     obj_model, model_gpm = obj_model_func(theta_obj, arg_dict['obj_dict'])
 
     totalmask = thismask & model_gpm
@@ -446,24 +558,6 @@ def tellfit(flux, thismask, arg_dict, init_from_last=None):
     fits for any object model that the user provides.
 
     Args:
-
-        theta (`numpy.ndarray`_):
-
-            Parameter vector for the object + telluric model.
-
-            This is actually two concatenated paramter vectors, one for
-            the object and one for the telluric, i.e::
-
-                theta_obj = theta[:-7]
-                theta_tell = theta[-7:]
-
-            The telluric model theta_tell is currently hard wired to be six dimensional::
-
-                pressure, temperature, humidity, airmass, resln, shift = theta_tell
-
-            The object model theta_obj can have an arbitrary size and is
-            provided as an argument to obj_model_func
-
         flux (`numpy.ndarray`_):
             The flux of the object being fit
         thismask (`numpy.ndarray`_, boolean):
@@ -477,8 +571,8 @@ def tellfit(flux, thismask, arg_dict, init_from_last=None):
                 - ``arg_dict['flux_ivar']``:  Inverse variance for the
                   flux array
                 - ``arg_dict['tell_dict']``: Dictionary containing the
-                  telluric grid and its parameters read in by
-                  read_telluric_grid
+                  telluric model and its parameters read in by
+                  read_telluric_pca or read_telluric_grid
                 - ``arg_dict['ind_lower']``: Lower index into the
                   telluric model wave_grid to trim down the telluric
                   model.
@@ -491,14 +585,11 @@ def tellfit(flux, thismask, arg_dict, init_from_last=None):
                   object model arguments which is passed to the
                   obj_model_func
 
-        init_from_last (object):
+        init_from_last (object, optional):
              Optional. Result object returned by the differential
              evolution optimizer for the last iteration. If this is passed the code
              will initialize from the previous best-fit for faster convergence.
 
-        **kwargs_opt (dict):
-            Optional arguments for the differential evolution
-            optimization
 
     Returns:
         tuple:  Returns three objects:
@@ -525,6 +616,13 @@ def tellfit(flux, thismask, arg_dict, init_from_last=None):
     nparams = len(bounds) # Number of parameters in the model
     popsize = arg_dict['popsize'] # Note this does nothing if the init is done from a previous iteration or optimum
     nsamples = arg_dict['popsize']*nparams
+    teltype = arg_dict['tell_dict']['teltype']
+    # FD: Assumes shift and stretch are turned on.
+    if teltype == 'pca':
+        ntheta_tell = arg_dict['tell_npca']+3 # Total number of telluric model parameters in PCA mode
+    elif teltype == 'grid':
+        ntheta_tell = 4+3 # Total number of telluric model parameters in grid mode
+
     # Decide how to initialize
     if init_from_last is not None:
         # Use a Gaussian ball about the optimum from a previous iteration
@@ -538,9 +636,9 @@ def tellfit(flux, thismask, arg_dict, init_from_last=None):
         init_obj = np.array([[np.clip(param + ballsize*(bounds_obj[i][1] - bounds_obj[i][0]) * rng.standard_normal(1)[0],
                                       bounds_obj[i][0], bounds_obj[i][1]) for i, param in enumerate(arg_dict['obj_dict']['init_obj_opt_theta'])]
                              for jsamp in range(nsamples)])
-        tell_lhs = utils.lhs(7, samples=nsamples, seed_or_rng=rng)
+        tell_lhs = utils.lhs(ntheta_tell, samples=nsamples)
         init_tell = np.array([[bounds[-idim][0] + tell_lhs[isamp, idim] * (bounds[-idim][1] - bounds[-idim][0])
-                               for idim in range(7)] for isamp in range(nsamples)])
+                               for idim in range(ntheta_tell)] for isamp in range(nsamples)])
         init = np.hstack((init_obj, init_tell))
     else:
         # If this is the first iteration and no object model optimum is presented, use a latin hypercube which is the default
@@ -550,10 +648,11 @@ def tellfit(flux, thismask, arg_dict, init_from_last=None):
                                                    init = init, updating='immediate', popsize=popsize,
                                                    recombination=arg_dict['recombination'], maxiter=arg_dict['diff_evol_maxiter'],
                                                    polish=arg_dict['polish'], disp=arg_dict['disp'])
-    theta_obj  = result.x[:-7]
-    theta_tell = result.x[-7:]
+                                        
+    theta_obj  = result.x[:-ntheta_tell]
+    theta_tell = result.x[-ntheta_tell:]
     tell_model = eval_telluric(theta_tell, arg_dict['tell_dict'],
-                               ind_lower=arg_dict['ind_lower'], ind_upper=arg_dict['ind_upper'])
+                                 ind_lower=arg_dict['ind_lower'], ind_upper=arg_dict['ind_upper'])
     obj_model, modelmask = obj_model_func(theta_obj, arg_dict['obj_dict'])
     totalmask = thismask & modelmask
     chi_vec = totalmask*(flux - tell_model*obj_model)*np.sqrt(flux_ivar)
@@ -571,7 +670,7 @@ def tellfit(flux, thismask, arg_dict, init_from_last=None):
 
 # TODO This should be a general reader once we get our act together with the data model.
 #  For echelle:  read in all the orders into a (nspec, nporders) array
-#  FOr longslit: read in the stanard into a (nspec, 1) array
+#  For longslit: read in the standard into a (nspec, 1) array
 def unpack_orders(sobjs, ret_flam=False):
     """
     Utility function to unpack the sobjs object and return the 
@@ -888,8 +987,8 @@ def init_qso_model(obj_params, iord, wave, flux, ivar, mask, tellmodel):
 
     tellmodel : array shape (nspec,)
         This is a telluric model computed on the wave wavelength grid. Initialization usually requires some initial
-        best guess for the telluric absorption, which is computed from the midpoint of the telluric model grid parameter
-        space using the resolution of the spectrograph and the airmass of the observations.
+        best guess for the telluric absorption, which is computed from the mean of the telluric model grid using
+        the resolution of the spectrograph.
 
     Returns
     -------
@@ -905,17 +1004,17 @@ def init_qso_model(obj_params, iord, wave, flux, ivar, mask, tellmodel):
 
     """
 
-    pca_dict = init_pca(obj_params['pca_file'], wave, obj_params['z_qso'], obj_params['npca'])
-    pca_mean = np.exp(pca_dict['components'][0, :])
+    qso_pca_dict = qso_init_pca(obj_params['pca_file'], wave, obj_params['z_qso'], obj_params['npca'])
+    qso_pca_mean = np.exp(qso_pca_dict['components'][0, :])
     tell_mask = tellmodel > obj_params['tell_norm_thresh']
     # Create a reference model and bogus noise
-    flux_ref = pca_mean * tellmodel
-    ivar_ref = utils.inverse((pca_mean/100.0) ** 2)
+    flux_ref = qso_pca_mean * tellmodel
+    ivar_ref = utils.inverse((qso_pca_mean/100.0) ** 2)
     flam_norm_inv = coadd.robust_median_ratio(flux, ivar, flux_ref, ivar_ref, mask=mask, mask_ref=tell_mask)
     flam_norm = 1.0/flam_norm_inv
 
     # Set the bounds for the PCA and truncate to the right dimension
-    coeffs = pca_dict['coeffs'][:,1:obj_params['npca']]
+    coeffs = qso_pca_dict['coeffs'][:,1:obj_params['npca']]
     # Compute the min and max arrays of the coefficients which are not the norm, i.e. grab the coeffs that aren't the first one
     coeff_min = np.amin(coeffs, axis=0)  # only
     coeff_max = np.amax(coeffs, axis=0)
@@ -925,7 +1024,7 @@ def init_qso_model(obj_params, iord, wave, flux, ivar, mask, tellmodel):
     bounds_pca = [(i, j) for i, j in zip(coeff_min, coeff_max)]        # Coefficients:  determined from PCA model
     bounds_obj = bounds_z + bounds_flam + bounds_pca
     # Create the obj_dict
-    obj_dict = dict(npca=obj_params['npca'], pca_dict=pca_dict)
+    obj_dict = dict(npca=obj_params['npca'], pca_dict=qso_pca_dict)
 
     return obj_dict, bounds_obj
 
@@ -946,23 +1045,21 @@ def eval_qso_model(theta, obj_dict):
 
     Returns
     -------
-    pca_model : `numpy.ndarray`_
-        PCA vectors were already interpolated onto the telluric model grid by init_qso_model.
-        array with same shape as the PCA vectors (tored in the obj_dict['pca_dict'])
+    qso_pca_model : array with same shape as the PCA vectors (stored in the obj_dict['pca_dict'])
+       PCA vectors were already interpolated onto the telluric model grid by init_qso_model
 
-    gpm : `numpy.ndarray`_ 
-        Good pixel mask indicating where the model is valid.
-        array with same shape as the pca_model
+    gpm : `numpy.ndarray`_ : array with same shape as the qso_pca_model
+       Good pixel mask indicating where the model is valid
 
     """
 
-    pca_model = pca_eval(theta, obj_dict['pca_dict'])
+    qso_pca_model = qso_pca_eval(theta, obj_dict['pca_dict'])
     # TODO Is the prior evaluation slowing things down??
     # TODO Disablingthe prior for now as I think it slows things down for no big gain
     #ln_pca_pri = qso_pca.pca_lnprior(theta_PCA, arg_dict['pca_dict'])
     #ln_pca_pri = 0.0
     #flux_model, tell_model, spec_model, modelmask
-    return pca_model, (pca_model > 0.0)
+    return qso_pca_model, (qso_pca_model > 0.0)
 
 
 ##############
@@ -1202,8 +1299,10 @@ def eval_poly_model(theta, obj_dict):
 
 
 def sensfunc_telluric(wave, counts, counts_ivar, counts_mask, exptime, airmass, std_dict,
-                      telgridfile, log10_blaze_function=None, ech_orders=None, polyorder=8, mask_hydrogen_lines=True,
-                      mask_helium_lines=False, hydrogen_mask_wid=10., resln_guess=None, resln_frac_bounds=(0.5, 1.5),
+                      telgridfile, log10_blaze_function=None, ech_orders=None, polyorder=8,
+                      tell_npca=4, teltype='pca',
+                      mask_hydrogen_lines=True, mask_helium_lines=False, hydrogen_mask_wid=10.,
+                      resln_guess=None, resln_frac_bounds=(0.6, 1.4), pix_shift_bounds=(-5.0, 5.0),
                       delta_coeff_bounds=(-20.0, 20.0), minmax_coeff_bounds=(-5.0, 5.0),
                       sn_clip=30.0, ballsize=5e-4, only_orders=None, maxiter=3, lower=3.0,
                       upper=3.0, tol=1e-3, popsize=30, recombination=0.7, polish=True, disp=False,
@@ -1252,6 +1351,10 @@ def sensfunc_telluric(wave, counts, counts_ivar, counts_mask, exptime, airmass,
         If passed, provides the true order numbers for the spectra provided.
     polyorder : :obj:`int`, optional, default = 8
         Polynomial order for the sensitivity function fit.
+    teltype : :obj:`str`, optional, default = 'pca'
+        Method for evaluating telluric models, either `pca` or `grid`.
+    tell_npca : :obj:`int`, optional, default = 4
+        Number of telluric PCA components used, must be <= 10
     mask_hydrogen_lines : :obj:`bool`, optional
         If True, mask stellar hydrogen absorption lines before fitting sensitivity function. Default = True
     mask_helium_lines : :obj:`bool`, optional
@@ -1365,7 +1468,9 @@ def sensfunc_telluric(wave, counts, counts_ivar, counts_mask, exptime, airmass,
 
     # Since we are fitting a sensitivity function, first compute counts per second per angstrom.
     TelObj = Telluric(wave, counts, counts_ivar, mask_tot, telgridfile, obj_params,
-                      init_sensfunc_model, eval_sensfunc_model, log10_blaze_function=log10_blaze_function, ech_orders=ech_orders,
+                      init_sensfunc_model, eval_sensfunc_model, log10_blaze_function=log10_blaze_function,
+                      teltype=teltype, tell_npca=tell_npca,
+                      ech_orders=ech_orders, pix_shift_bounds=pix_shift_bounds,
                       resln_guess=resln_guess, resln_frac_bounds=resln_frac_bounds, sn_clip=sn_clip,
                       maxiter=maxiter,  lower=lower, upper=upper, tol=tol,
                       popsize=popsize, recombination=recombination, polish=polish, disp=disp,
@@ -1408,11 +1513,12 @@ def create_bal_mask(wave, bal_wv_min_max):
 
 
 
-def qso_telluric(spec1dfile, telgridfile, pca_file, z_qso, telloutfile, outfile, npca=8,
+def qso_telluric(spec1dfile, telgridfile,  pca_file, z_qso, telloutfile, outfile, npca=8,
                  pca_lower=1220.0, pca_upper=3100.0, bal_wv_min_max=None, delta_zqso=0.1,
+                 teltype='pca', tell_npca=4,
                  bounds_norm=(0.1, 3.0), tell_norm_thresh=0.9, sn_clip=30.0, only_orders=None,
                  maxiter=3, tol=1e-3, popsize=30, recombination=0.7, polish=True, disp=False,
-                 debug_init=False, debug=False, show=False):
+                 pix_shift_bounds=(-5.0,5.0), debug_init=False, debug=False, show=False):
     """
     Fit and correct a QSO spectrum for telluric absorption.
 
@@ -1445,6 +1551,10 @@ def qso_telluric(spec1dfile, telgridfile, pca_file, z_qso, telloutfile, outfile,
     delta_zqso : :obj:`float`, optional
         During the fit, the QSO redshift is allowed to vary within
         ``+/-delta_zqso``.
+    teltype : :obj:`str`, optional, default = 'pca'
+        Method for evaluating telluric models, either `pca` or `grid`.
+    tell_npca : :obj:`int`, optional, default = 4
+        Number of telluric PCA components used, must be <=10
     bounds_norm : :obj:`tuple`, optional
         A two-tuple with the lower and upper bounds on the fractional adjustment
         of the flux in the QSO model spectrum.  For example, a value of ``(0.1,
@@ -1527,14 +1637,16 @@ def qso_telluric(spec1dfile, telgridfile, pca_file, z_qso, telloutfile, outfile,
 
     # parameters lowered for testing
     TelObj = Telluric(wave, flux, ivar, mask_tot, telgridfile, obj_params, init_qso_model,
-                      eval_qso_model, sn_clip=sn_clip, maxiter=maxiter, tol=tol, popsize=popsize,
-                      recombination=recombination, polish=polish, disp=disp, debug=debug)
+                      eval_qso_model, pix_shift_bounds=pix_shift_bounds,
+                      sn_clip=sn_clip, maxiter=maxiter, tol=tol, popsize=popsize, teltype=teltype,
+                      tell_npca=tell_npca, recombination=recombination, polish=polish,
+                      disp=disp, debug=debug)
     TelObj.run(only_orders=only_orders)
     TelObj.to_file(telloutfile, overwrite=True)
 
     # Apply the telluric correction
     telluric = TelObj.model['TELLURIC'][0,:]
-    pca_model = TelObj.model['OBJ_MODEL'][0,:]
+    qso_pca_model = TelObj.model['OBJ_MODEL'][0,:]
 
     # Plot the telluric corrected and rescaled orders
     flux_corr = flux/(telluric + (telluric == 0.0))
@@ -1554,11 +1666,11 @@ def qso_telluric(spec1dfile, telgridfile, pca_file, z_qso, telloutfile, outfile,
                  alpha=0.7, zorder=5)
         plt.plot(wave, sig_corr, drawstyle='steps-mid', color='r', label='noise', alpha=0.3,
                  zorder=1)
-        plt.plot(wave, pca_model, color='cornflowerblue', linewidth=1.0, label='PCA model',
+        plt.plot(wave, qso_pca_model, color='cornflowerblue', linewidth=1.0, label='PCA model',
                  zorder=7, alpha=0.7)
-        plt.plot(wave, pca_model.max()*0.9*telluric, color='magenta', drawstyle='steps-mid',
+        plt.plot(wave, qso_pca_model.max()*0.9*telluric, color='magenta', drawstyle='steps-mid',
                  label='telluric', alpha=0.4)
-        plt.ylim(-0.1*pca_model.max(), 1.5*pca_model.max())
+        plt.ylim(-0.1*qso_pca_model.max(), 1.5*qso_pca_model.max())
         plt.legend()
         plt.xlabel('Wavelength')
         plt.ylabel('Flux')
@@ -1566,18 +1678,18 @@ def qso_telluric(spec1dfile, telgridfile, pca_file, z_qso, telloutfile, outfile,
 
     # save the telluric corrected spectrum
     save_coadd1d_tofits(outfile, wave, flux_corr, ivar_corr, mask_corr, wave_grid_mid=wave_grid_mid,
-                        spectrograph=header['PYP_SPEC'], telluric=telluric, obj_model=pca_model,
+                        spectrograph=header['PYP_SPEC'], telluric=telluric, obj_model=qso_pca_model,
                         header=header, ex_value='OPT', overwrite=True)
 
     return TelObj
 
-def star_telluric(spec1dfile, telgridfile, telloutfile, outfile, star_type=None, star_mag=None,
-                  star_ra=None, star_dec=None, func='legendre', model='exp', polyorder=5,
-                  mask_hydrogen_lines=True, mask_helium_lines=False, hydrogen_mask_wid=10.,
-                  delta_coeff_bounds=(-20.0, 20.0),
+def star_telluric(spec1dfile, telgridfile, telloutfile, outfile, star_type=None,
+                  star_mag=None, star_ra=None, star_dec=None, func='legendre', model='exp',
+                  polyorder=5, teltype='pca', tell_npca=4, mask_hydrogen_lines=True,
+                  mask_helium_lines=False, hydrogen_mask_wid=10., delta_coeff_bounds=(-20.0, 20.0),
                   minmax_coeff_bounds=(-5.0, 5.0), only_orders=None, sn_clip=30.0, maxiter=3,
                   tol=1e-3, popsize=30, recombination=0.7, polish=True, disp=False,
-                  debug_init=False, debug=False, show=False):
+                  pix_shift_bounds=(-5.0,5.0), debug_init=False, debug=False, show=False):
     """
     This needs a doc string.
 
@@ -1634,7 +1746,9 @@ def star_telluric(spec1dfile, telgridfile, telloutfile, outfile, star_type=None,
 
     # parameters lowered for testing
     TelObj = Telluric(wave, flux, ivar, mask_tot, telgridfile, obj_params, init_star_model,
-                      eval_star_model,  sn_clip=sn_clip, tol=tol, popsize=popsize,
+                      eval_star_model, pix_shift_bounds=pix_shift_bounds,
+                      teltype=teltype, tell_npca=tell_npca,
+                      sn_clip=sn_clip, tol=tol, popsize=popsize,
                       recombination=recombination, polish=polish, disp=disp, debug=debug)
     TelObj.run(only_orders=only_orders)
     TelObj.to_file(telloutfile, overwrite=True)
@@ -1679,11 +1793,11 @@ def star_telluric(spec1dfile, telgridfile, telloutfile, outfile, star_type=None,
     return TelObj
 
 def poly_telluric(spec1dfile, telgridfile, telloutfile, outfile, z_obj=0.0, func='legendre',
-                  model='exp', polyorder=3, fit_wv_min_max=None, mask_lyman_a=True,
-                  delta_coeff_bounds=(-20.0, 20.0), minmax_coeff_bounds=(-5.0, 5.0),
+                  model='exp', polyorder=3, fit_wv_min_max=None, mask_lyman_a=True, teltype='pca',
+                  tell_npca=4, delta_coeff_bounds=(-20.0, 20.0), minmax_coeff_bounds=(-5.0, 5.0),
                   only_orders=None, sn_clip=30.0, maxiter=3, tol=1e-3, popsize=30,
-                  recombination=0.7, polish=True, disp=False, debug_init=False, debug=False,
-                  show=False):
+                  recombination=0.7, polish=True, disp=False, pix_shift_bounds=(-5.0,5.0),
+                  debug_init=False, debug=False, show=False):
     """
     This needs a doc string.
 
@@ -1732,8 +1846,9 @@ def poly_telluric(spec1dfile, telgridfile, telloutfile, outfile, z_obj=0.0, func
 
     # parameters lowered for testing
     TelObj = Telluric(wave, flux, ivar, mask_tot, telgridfile, obj_params, init_poly_model,
-                      eval_poly_model, sn_clip=sn_clip, maxiter=maxiter, tol=tol, popsize=popsize,
-                      recombination=recombination, polish=polish, disp=disp, debug=debug)
+                      eval_poly_model, pix_shift_bounds=pix_shift_bounds,
+                      sn_clip=sn_clip, maxiter=maxiter, tol=tol, popsize=popsize, teltype=teltype,
+                      tell_npca=tell_npca, recombination=recombination, polish=polish, disp=disp, debug=debug)
     TelObj.run(only_orders=only_orders)
     TelObj.to_file(telloutfile, overwrite=True)
 
@@ -1783,7 +1898,12 @@ class Telluric(datamodel.DataContainer):
 
     The model telluric absorption spectra for the atmosphere are generated by
     HITRAN atmosphere models with a baseline atmospheric profile model of the
-    observatory in question. The modeling can be applied to both multi-slit
+    observatory in question. One can interpolate over grids of these absorption
+    spectra directly, which have been created for each observatory, or instead
+    (by default) use a PCA decomposition of the telluric models across all
+    observatories, which is more flexible and much lighter in file size.
+    
+    The modeling can be applied to both multi-slit
     or echelle data. The code performs a differential evolution optimization
     using `scipy.optimize.differential_evolution`_. Several iterations are
     performed with rejection of outliers (governed by the maxiter parameter
@@ -1814,26 +1934,29 @@ class Telluric(datamodel.DataContainer):
     object model, the bounds depend on the nature of the object model, which
     is why ``init_obj_model`` must be provided.
 
-    The telluric model is governed by seven parameters --- ``pressure``,
-    ``temperature``, ``humidity``, ``airmass``, ``resln``, ``shift``,
+    The telluric model is governed by seven parameters --- for the `grid`
+    approach: ``pressure``, ``temperature``, ``humidity``, ``airmass``; for
+    the `pca` approach 4 PCA coefficients; plus ``resln``, ``shift``, and
     ``stretch`` --- where ``resln`` is the resolution of the spectrograph and
     ``shift`` is a shift in pixels. The ``stretch`` is a stretch of the pixel
     scale. The airmass of the object will be used to initalize the fit (this
     helps with initalizing the object model), but the models are sufficiently
     flexible that often the best-fit airmass actually differs from the
-    airmass of the spectrum.
+    airmass of the spectrum. In the `pca` approach, the number of PCA
+    components used can be adjusted between 1 and 10, with a tradeoff between
+    speed and accuracy.
 
     This code can be run on stacked spectra covering a large range of
-    airmasses and will still provide good results. The resulting airmass will
-    be an effective value, and as per above may not have much relation to the
-    true airmass of the observation. The exception to this rule is extremely
-    high signal-to-noise ratio data S/N > 30, for which it can be difficult
-    to obtain a good fit within the noise of the data. In such cases, the
-    user should split up the data into chunks taken around the same airmass,
-    fit those individually with this class, and then combine the resulting
-    telluric corrected spectra. This will, in general, result in better fits,
-    and will also average down the residuals from the telluric model fit in
-    the final averaged spectrum.
+    airmasses and will still provide good results. The resulting airmass
+    (in the `grid` approach) will be an effective value, and as per above may
+    not have much relation to the true airmass of the observation. The
+    exception to this rule is extremely high signal-to-noise ratio data
+    (S/N > 30), for which it can be difficult to obtain a good fit within the
+    noise of the data. In such cases, the user should split up the data into
+    chunks taken around the same airmass, fit those individually with this
+    class, and then combine the resulting telluric corrected spectra. This
+    will, in general, result in better fits, and will also average down the
+    residuals from the telluric model fit in the final averaged spectrum.
 
     The datamodel attributes are:
 
@@ -1856,8 +1979,19 @@ class Telluric(datamodel.DataContainer):
             Good pixel gpm for the object in question. Same shape as
             ``wave``.
         telgridfile (:obj:`str`):
-            File containing grid of HITRAN atmosphere models; see
+            File containing grid of HITRAN atmosphere models or PCA
+            decomposition of such models, see
             :class:`~pypeit.par.pypeitpar.TelluricPar`.
+        teltype (:obj:`str`, optional):
+            Method for evaluating telluric models, either `pca` or `grid`.
+            The `grid` method directly interpolates a pre-computed grid of
+            HITRAN atmospheric models with nearest grid-point interpolation,
+            while the `pca` method instead fits for the coefficients of a
+            PCA decomposition of HITRAN atmospheric models, enabling a more
+            flexible and far more file-size efficient interpolation
+            of the telluric absorption model space.
+        tell_npca (:obj:`int`, optional):
+            Number of telluric PCA components used, must be <=10
         obj_params (:obj:`dict`):
             Dictionary of parameters for initializing the object model.
         init_obj_model (callable):
@@ -2043,7 +2177,12 @@ class Telluric(datamodel.DataContainer):
     """Datamodel version."""
 
     datamodel = {'telgrid': dict(otype=str,
-                                 descr='File containing grid of HITRAN atmosphere models'),
+                                 descr='File containing PCA components or grid of '
+                                       'HITRAN atmosphere models'),
+                 'teltype': dict(otype=str,
+                                 descr='Type of telluric model, `pca` or `grid`'),
+                 'tell_npca': dict(otype=int,
+                                    descr='Number of telluric PCA components used'),
                  'std_src': dict(otype=str, descr='Name of the standard source'),
                  'std_name': dict(otype=str, descr='Type of standard source'),
                  'std_cal': dict(otype=str,
@@ -2141,7 +2280,7 @@ class Telluric(datamodel.DataContainer):
                 ]
 
     @staticmethod
-    def empty_model_table(norders, nspec, n_obj_par=0):
+    def empty_model_table(norders, nspec, tell_npca=4, n_obj_par=0):
         """
         Construct an empty `astropy.table.Table`_ for the telluric model
         results.
@@ -2151,6 +2290,8 @@ def empty_model_table(norders, nspec, n_obj_par=0):
                 The number of slits/orders on the detector.
             nspec (:obj:`int`):
                 The number of spectral pixels on the detector.
+            tell_npca (:obj:`int`):
+                Number of telluric model parameters
             n_obj_par (:obj:`int`, optional):
                 The number of parameters used to construct the object model
                 spectrum.
@@ -2166,16 +2307,10 @@ def empty_model_table(norders, nspec, n_obj_par=0):
             table.Column(name='OBJ_MODEL', dtype=float, length=norders, shape=(nspec,),
                          description='Best-fitting object model spectrum'),
             # TODO: Why do we need both TELL_THETA and all the individual parameters...
-            table.Column(name='TELL_THETA', dtype=float, length=norders, shape=(7,),
+            table.Column(name='TELL_THETA', dtype=float, length=norders, shape=(tell_npca+3,),
                          description='Best-fitting telluric model parameters'),
-            table.Column(name='TELL_PRESS', dtype=float, length=norders,
-                         description='Best-fitting telluric model pressure'),
-            table.Column(name='TELL_TEMP', dtype=float, length=norders,
-                         description='Best-fitting telluric model temperature'),
-            table.Column(name='TELL_H2O', dtype=float, length=norders,
-                         description='Best-fitting telluric model humidity'),
-            table.Column(name='TELL_AIRMASS', dtype=float, length=norders,
-                         description='Best-fitting telluric model airmass'),
+            table.Column(name='TELL_PARAM', dtype=float, length=norders, shape=(tell_npca,),
+                         description='Best-fitting telluric atmospheric parameters or PCA coefficients'),
             table.Column(name='TELL_RESLN', dtype=float, length=norders,
                          description='Best-fitting telluric model spectral resolution'),
             table.Column(name='TELL_SHIFT', dtype=float, length=norders,
@@ -2203,9 +2338,9 @@ def empty_model_table(norders, nspec, n_obj_par=0):
             table.Column(name='WAVE_MAX', dtype=float, length=norders,
                          description='Maximum wavelength included in the fit')])
 
-    def __init__(self, wave, flux, ivar, gpm, telgridfile, obj_params, init_obj_model,
-                 eval_obj_model, log10_blaze_function=None, ech_orders=None, sn_clip=30.0, airmass_guess=1.5,
-                 resln_guess=None, resln_frac_bounds=(0.5, 1.5), pix_shift_bounds=(-5.0, 5.0),
+    def __init__(self, wave, flux, ivar, gpm, telgridfile, obj_params, init_obj_model, eval_obj_model,
+                 log10_blaze_function=None, ech_orders=None, sn_clip=30.0, teltype='pca', tell_npca=4,
+                 airmass_guess=1.5, resln_guess=None, resln_frac_bounds=(0.3, 1.5), pix_shift_bounds=(-5.0, 5.0),
                  pix_stretch_bounds=(0.9,1.1), maxiter=2, sticky=True, lower=3.0, upper=3.0,
                  seed=777, ballsize = 5e-4, tol=1e-3, diff_evol_maxiter=1000,  popsize=30,
                  recombination=0.7, polish=True, disp=False, sensfunc=False, debug=False):
@@ -2223,8 +2358,10 @@ def __init__(self, wave, flux, ivar, gpm, telgridfile, obj_params, init_obj_mode
 
         # 1) Assign arguments
         self.telgrid = telgridfile
+        self.teltype = teltype
         self.obj_params = obj_params
         self.init_obj_model = init_obj_model
+        self.tell_npca = tell_npca
         self.airmass_guess = airmass_guess
         self.eval_obj_model = eval_obj_model
         self.ech_orders = ech_orders
@@ -2263,8 +2400,15 @@ def __init__(self, wave, flux, ivar, gpm, telgridfile, obj_params, init_obj_mode
                 wave, flux, ivar, gpm)
         # 3) Read the telluric grid and initalize associated parameters
         wv_gpm = self.wave_in_arr > 1.0
-        self.tell_dict = read_telluric_grid(self.telgrid, wave_min=self.wave_in_arr[wv_gpm].min(),
-                                            wave_max=self.wave_in_arr[wv_gpm].max())
+        if self.teltype == 'pca':
+            self.tell_dict = read_telluric_pca(self.telgrid, wave_min=self.wave_in_arr[wv_gpm].min(),
+                                               wave_max=self.wave_in_arr[wv_gpm].max())
+        elif self.teltype == 'grid':
+            self.tell_npca = 4
+            self.tell_dict = read_telluric_grid(self.telgrid, wave_min=self.wave_in_arr[wv_gpm].min(),
+                                                wave_max=self.wave_in_arr[wv_gpm].max())
+
+            
         self.wave_grid = self.tell_dict['wave_grid']
         self.ngrid = self.wave_grid.size
         self.resln_guess = wvutils.get_sampling(self.wave_in_arr)[2] \
@@ -2313,8 +2457,8 @@ def __init__(self, wave, flux, ivar, gpm, telgridfile, obj_params, init_obj_mode
             msgs.info(f'Initializing object model for order: {iord}, {counter}/{self.norders}'
                       + f' with user supplied function: {self.init_obj_model.__name__}')
             tellmodel = eval_telluric(self.tell_guess, self.tell_dict,
-                                      ind_lower=self.ind_lower[iord],
-                                      ind_upper=self.ind_upper[iord])
+                                        ind_lower=self.ind_lower[iord],
+                                        ind_upper=self.ind_upper[iord])
             # TODO This is a pretty ugly way to pass in the blaze function. Particularly since now all the other models
             #  (star, qso, poly) are going to have this parameter set in their obj_params dictionary.
             #  Is there something more elegant that can be done with e.g. functools.partial?
@@ -2336,6 +2480,7 @@ def __init__(self, wave, flux, ivar, gpm, telgridfile, obj_params, init_obj_mode
             arg_dict_iord = dict(ivar=self.ivar_arr[self.ind_lower[iord]:self.ind_upper[iord]+1,iord],
                                  tell_dict=self.tell_dict, ind_lower=self.ind_lower[iord],
                                  ind_upper=self.ind_upper[iord],
+                                 tell_npca=self.tell_npca,
                                  obj_model_func=self.eval_obj_model, obj_dict=obj_dict,
                                  ballsize=self.ballsize, bounds=bounds_iord, rng=self.rng,
                                  diff_evol_maxiter=self.diff_evol_maxiter, tol=self.tol,
@@ -2380,8 +2525,12 @@ def run(self, only_orders=None):
                                               inmask=self.mask_arr[self.ind_lower[iord]:self.ind_upper[iord]+1,iord],
                                               maxiter=self.maxiter, lower=self.lower,
                                               upper=self.upper, sticky=self.sticky)
-            self.theta_obj_list[iord] = self.result_list[iord].x[:-7]
-            self.theta_tell_list[iord] = self.result_list[iord].x[-7:]
+            if self.teltype == 'pca':
+                self.theta_obj_list[iord] = self.result_list[iord].x[:-(self.tell_npca+3)]
+                self.theta_tell_list[iord] = self.result_list[iord].x[-(self.tell_npca+3):]
+            elif self.teltype == 'grid':
+                self.theta_obj_list[iord] = self.result_list[iord].x[:-(4+3)]
+                self.theta_tell_list[iord] = self.result_list[iord].x[-(4+3):]
             self.obj_model_list[iord], modelmask \
                     = self.eval_obj_model(self.theta_obj_list[iord], self.obj_dict_list[iord])
             self.tellmodel_list[iord] = eval_telluric(self.theta_tell_list[iord], self.tell_dict,
@@ -2432,7 +2581,7 @@ def init_output(self):
         """
 
         self.model = self.empty_model_table(self.norders, self.nspec_in,
-                                            n_obj_par=self.max_ntheta_obj)
+                                            tell_npca=self.tell_npca, n_obj_par=self.max_ntheta_obj)
         if 'output_meta_keys' in self.obj_params:
             for key in self.obj_params['output_meta_keys']:
                 if key.lower() in self.datamodel.keys():
@@ -2460,6 +2609,11 @@ def assign_output(self, iord):
         gdwave = self.wave_in_arr[:,iord] > 1.0
         wave_in_gd = self.wave_in_arr[gdwave,iord]
         wave_grid_now = self.wave_grid[self.ind_lower[iord]:self.ind_upper[iord]+1]
+        
+        if self.teltype == 'pca':
+            ntell = self.tell_npca
+        elif self.teltype == 'grid':
+            ntell = 4
 
         self.model['WAVE'][iord] = self.wave_in_arr[:,iord]
         self.model['TELLURIC'][iord][gdwave] \
@@ -2471,13 +2625,10 @@ def assign_output(self, iord):
                                              kind='linear', bounds_error=False,
                                              fill_value=0.0)(wave_in_gd)
         self.model['TELL_THETA'][iord] = self.theta_tell_list[iord]
-        self.model['TELL_PRESS'][iord] = self.theta_tell_list[iord][0]
-        self.model['TELL_TEMP'][iord] = self.theta_tell_list[iord][1]
-        self.model['TELL_H2O'][iord] = self.theta_tell_list[iord][2]
-        self.model['TELL_AIRMASS'][iord] = self.theta_tell_list[iord][3]
-        self.model['TELL_RESLN'][iord] = self.theta_tell_list[iord][4]
-        self.model['TELL_SHIFT'][iord] = self.theta_tell_list[iord][5]
-        self.model['TELL_STRETCH'][iord] = self.theta_tell_list[iord][6]
+        self.model['TELL_PARAM'][iord] = self.theta_tell_list[iord][:ntell]
+        self.model['TELL_RESLN'][iord] = self.theta_tell_list[iord][ntell]
+        self.model['TELL_SHIFT'][iord] = self.theta_tell_list[iord][ntell+1]
+        self.model['TELL_STRETCH'][iord] = self.theta_tell_list[iord][ntell+2]
         ntheta_iord = len(self.theta_obj_list[iord])
         self.model['OBJ_THETA'][iord][:ntheta_iord] = self.theta_obj_list[iord]
         self.model['CHI2'][iord] = self.result_list[iord].fun
@@ -2552,14 +2703,20 @@ def get_tell_guess(self):
         function.
 
         Returns:
-            :obj:`tuple`: The guess pressure, temperature, humidity,
-            airmass, resolution, shift, and stretch parameters. The first
-            three are the median of the parameters covered by the telluric
-            model grid.
+            :obj:`tuple`: The guess telluric model parameters,
+            resolution, shift, and stretch parameters.
         """
-        return np.median(self.tell_dict['pressure_grid']), np.median(self.tell_dict['temp_grid']), \
-               np.median(self.tell_dict['h2o_grid']), self.airmass_guess, self.resln_guess, \
-               0.0, 1.0
+        if self.teltype == 'grid':
+            guess = [np.median(self.tell_dict['pressure_grid'])]
+            guess.append(np.median(self.tell_dict['temp_grid']))
+            guess.append(np.median(self.tell_dict['h2o_grid']))
+            guess.append(self.airmass_guess)
+        else: 
+            guess = list(np.zeros(self.tell_npca))
+        guess.append(self.resln_guess)
+        guess.append(0.0)
+        guess.append(1.0)
+        return tuple(guess)
 
     def get_bounds_tell(self):
         """
@@ -2567,18 +2724,30 @@ def get_bounds_tell(self):
 
         Returns:
             :obj:`list`: A list of two-tuples, where each two-tuple proives
-            the minimum and maximum allowed model values for the pressure,
-            temperature, humidity, airmass, resolution, shift, and
-            stretch parameters.
+            the minimum and maximum allowed model values for the PCA coefficients,
+            (or for the `grid` approach: pressure, temperature, humidity, airmass)
+            as well as the resolution, shift, and stretch parameters.
         """
         # Set the bounds for the optimization
-        return [(self.tell_dict['pressure_grid'].min(), self.tell_dict['pressure_grid'].max()),
-                (self.tell_dict['temp_grid'].min(), self.tell_dict['temp_grid'].max()),
-                (self.tell_dict['h2o_grid'].min(), self.tell_dict['h2o_grid'].max()),
-                (self.tell_dict['airmass_grid'].min(), self.tell_dict['airmass_grid'].max()),
-                (self.resln_guess * self.resln_frac_bounds[0],
-                 self.resln_guess * self.resln_frac_bounds[1]),
-                self.pix_shift_bounds, self.pix_stretch_bounds]
+        bounds = []
+        if self.teltype == 'grid':
+            bounds.append((self.tell_dict['pressure_grid'].min(),
+                           self.tell_dict['pressure_grid'].max()))
+            bounds.append((self.tell_dict['temp_grid'].min(),
+                           self.tell_dict['temp_grid'].max()))
+            bounds.append((self.tell_dict['h2o_grid'].min(),
+                           self.tell_dict['h2o_grid'].max()))
+            bounds.append((self.tell_dict['airmass_grid'].min(),
+                           self.tell_dict['airmass_grid'].max()))
+        else: # i.e. for teltype == 'pca'
+            for ii in range(self.tell_npca):
+                bounds.append((self.tell_dict['bounds_tell_pca'][0][ii+1],
+                               self.tell_dict['bounds_tell_pca'][1][ii+1]))
+        bounds.append((self.resln_guess * self.resln_frac_bounds[0],
+                       self.resln_guess * self.resln_frac_bounds[1]))
+        bounds.append(self.pix_shift_bounds)
+        bounds.append(self.pix_stretch_bounds)
+        return bounds
 
     def sort_telluric(self):
         """
@@ -2600,10 +2769,14 @@ def sort_telluric(self):
         # Do a quick loop over all the orders to sort them in order of strongest
         # to weakest telluric absorption
         for iord in range(self.norders):
-            tm_grid = self.tell_dict['tell_grid'][...,self.ind_lower[iord]:self.ind_upper[iord]+1]
-            tell_model_mid = tm_grid[tm_grid.shape[0]//2, tm_grid.shape[1]//2,
-                                     tm_grid.shape[2]//2, tm_grid.shape[3]//2, :]
-            tell_med[iord] = np.mean(tell_model_mid)
+            if self.teltype == 'grid':
+                tm_grid = self.tell_dict['tell_grid'][...,self.ind_lower[iord]:self.ind_upper[iord]+1]
+                tell_model_mid = tm_grid[tm_grid.shape[0]//2, tm_grid.shape[1]//2,
+                                         tm_grid.shape[2]//2, tm_grid.shape[3]//2, :]
+                tell_med[iord] = np.mean(tell_model_mid)
+            else:
+                tell_model_mean = self.tell_dict['tell_pca'][0,self.ind_lower[iord]:self.ind_upper[iord]+1]
+                tell_med[iord] = np.mean(np.exp(-np.sinh(tell_model_mean)))
 
         # Perform fits in order of telluric strength
         return tell_med.argsort()
diff --git a/pypeit/par/pypeitpar.py b/pypeit/par/pypeitpar.py
index 7fd2fadbce..7a9d681f08 100644
--- a/pypeit/par/pypeitpar.py
+++ b/pypeit/par/pypeitpar.py
@@ -2209,7 +2209,7 @@ class TelluricPar(ParSet):
     """
 
     def __init__(self, telgridfile=None, sn_clip=None, resln_guess=None, resln_frac_bounds=None, pix_shift_bounds=None,
-                 delta_coeff_bounds=None, minmax_coeff_bounds=None, maxiter=None,
+                 delta_coeff_bounds=None, minmax_coeff_bounds=None, maxiter=None, tell_npca=None, teltype=None,
                  sticky=None, lower=None, upper=None, seed=None, tol=None, popsize=None, recombination=None, polish=None,
                  disp=None, objmodel=None, redshift=None, delta_redshift=None, pca_file=None, npca=None,
                  bal_wv_min_max=None, bounds_norm=None, tell_norm_thresh=None, only_orders=None, pca_lower=None,
@@ -2224,6 +2224,7 @@ def __init__(self, telgridfile=None, sn_clip=None, resln_guess=None, resln_frac_
         # Initialize the other used specifications for this parameter
         # set
         defaults = OrderedDict.fromkeys(pars.keys())
+        options = OrderedDict.fromkeys(pars.keys())
         dtypes = OrderedDict.fromkeys(pars.keys())
         descr = OrderedDict.fromkeys(pars.keys())
 
@@ -2234,6 +2235,18 @@ def __init__(self, telgridfile=None, sn_clip=None, resln_guess=None, resln_frac_
                                'must be downloaded from the GoogleDrive and installed in your PypeIt installation via ' \
                                'the pypeit_install_telluric script. NOTE: This parameter no longer includes the full ' \
                                'pathname to the Telluric Grid file, but is just the filename of the grid itself.'
+        
+        defaults['tell_npca'] = 5
+        dtypes['tell_npca'] = int
+        descr['tell_npca'] = 'Number of telluric PCA components used. Can be set to any number from 1 to 10.'
+        
+        defaults['teltype'] = 'pca'
+        options['teltype'] = TelluricPar.valid_teltype()
+        dtypes['teltype'] = str
+        descr['teltype'] = 'Method used to evaluate telluric models, either pca or grid. The grid option uses a ' \
+                           'fixed grid of pre-computed HITRAN+LBLRTM atmospheric transmission models for each ' \
+                           'observatory, whereas the pca option uses principal components of a larger model grid ' \
+                           'to compute an accurate pseudo-telluric model with a much lighter telgridfile.'
 
         defaults['sn_clip'] = 30.0
         dtypes['sn_clip'] = [int, float]
@@ -2254,12 +2267,12 @@ def __init__(self, telgridfile=None, sn_clip=None, resln_guess=None, resln_frac_
 
 
         pars['resln_frac_bounds'] = tuple_force(pars['resln_frac_bounds'])
-        defaults['resln_frac_bounds'] = (0.5,1.5)
+        defaults['resln_frac_bounds'] = (0.6,1.4)
         dtypes['resln_frac_bounds'] = tuple
         descr['resln_frac_bounds'] = 'Bounds for the resolution fit optimization which is part of the telluric model. ' \
-                                     'This range is in units of the resln_guess, so the (0.5, 1.5) would bound the ' \
+                                     'This range is in units of the resln_guess, so the (0.6, 1.4) would bound the ' \
                                      'spectral resolution fit to be within the range ' \
-                                     'bounds_resln = (0.5*resln_guess, 1.5*resln_guess)'
+                                     'bounds_resln = (0.6*resln_guess, 1.4*resln_guess)'
 
         pars['pix_shift_bounds'] = tuple_force(pars['pix_shift_bounds'])
         defaults['pix_shift_bounds'] = (-5.0,5.0)
@@ -2462,7 +2475,7 @@ def __init__(self, telgridfile=None, sn_clip=None, resln_guess=None, resln_frac_
     @classmethod
     def from_dict(cls, cfg):
         k = np.array([*cfg.keys()])
-        parkeys = ['telgridfile', 'sn_clip', 'resln_guess', 'resln_frac_bounds',
+        parkeys = ['telgridfile', 'teltype', 'sn_clip', 'resln_guess', 'resln_frac_bounds', 'tell_npca',
                    'pix_shift_bounds', 'delta_coeff_bounds', 'minmax_coeff_bounds',
                    'maxiter', 'sticky', 'lower', 'upper', 'seed', 'tol',
                    'popsize', 'recombination', 'polish', 'disp', 'objmodel','redshift', 'delta_redshift',
@@ -2479,12 +2492,27 @@ def from_dict(cls, cfg):
         for pk in parkeys:
             kwargs[pk] = cfg[pk] if pk in k else None
         return cls(**kwargs)
+        
+    @staticmethod
+    def valid_teltype():
+        """
+        Return the valid telluric methods.
+        """
+        return ['pca', 'grid']
 
     def validate(self):
         """
         Check the parameters are valid for the provided method.
         """
-        pass
+        if self.data['tell_npca'] < 1 or self.data['tell_npca'] > 10:
+            raise ValueError('Invalid value {:d} for tell_npca '.format(self.data['tell_npca'])+
+                             '(must be between 1 and 10).')
+                             
+        self.data['teltype'] = self.data['teltype'].lower()
+        if self.data['teltype'] not in TelluricPar.valid_teltype():
+            raise ValueError('Invalid teltype "{}"'.format(self.data['teltype'])+
+                             ', valid options are: {}.'.format(TelluricPar.valid_teltype()))
+        
         # JFH add something in here which checks that the recombination value provided is bewteen 0 and 1, although
         # scipy.optimize.differential_evoluiton probalby checks this.
 
diff --git a/pypeit/scripts/tellfit.py b/pypeit/scripts/tellfit.py
index 20ee30ae75..d0e0626f2f 100644
--- a/pypeit/scripts/tellfit.py
+++ b/pypeit/scripts/tellfit.py
@@ -155,13 +155,17 @@ def main(args):
                                            par['telluric']['pca_file'],
                                            par['telluric']['redshift'], modelfile, outfile,
                                            npca=par['telluric']['npca'],
+                                           teltype=par['telluric']['teltype'], tell_npca=par['telluric']['tell_npca'],
                                            pca_lower=par['telluric']['pca_lower'],
                                            pca_upper=par['telluric']['pca_upper'],
                                            bounds_norm=par['telluric']['bounds_norm'],
                                            tell_norm_thresh=par['telluric']['tell_norm_thresh'],
                                            only_orders=par['telluric']['only_orders'],
                                            bal_wv_min_max=par['telluric']['bal_wv_min_max'],
+                                           pix_shift_bounds=par['telluric']['pix_shift_bounds'],
                                            maxiter=par['telluric']['maxiter'],
+                                           popsize=par['telluric']['popsize'],
+                                           tol=par['telluric']['tol'],
                                            debug_init=args.debug, disp=args.debug,
                                            debug=args.debug, show=args.plot)
         elif par['telluric']['objmodel']=='star':
@@ -175,12 +179,16 @@ def main(args):
                                              model=par['telluric']['model'],
                                              polyorder=par['telluric']['polyorder'],
                                              only_orders=par['telluric']['only_orders'],
+                                             teltype=par['telluric']['teltype'], tell_npca=par['telluric']['tell_npca'],
                                              mask_hydrogen_lines=par['sensfunc']['mask_hydrogen_lines'],
                                              mask_helium_lines=par['sensfunc']['mask_helium_lines'],
                                              hydrogen_mask_wid=par['sensfunc']['hydrogen_mask_wid'],
                                              delta_coeff_bounds=par['telluric']['delta_coeff_bounds'],
                                              minmax_coeff_bounds=par['telluric']['minmax_coeff_bounds'],
+                                             pix_shift_bounds=par['telluric']['pix_shift_bounds'],
                                              maxiter=par['telluric']['maxiter'],
+                                             popsize=par['telluric']['popsize'],
+                                             tol=par['telluric']['tol'],
                                              debug_init=args.debug, disp=args.debug,
                                              debug=args.debug, show=args.plot)
         elif par['telluric']['objmodel']=='poly':
@@ -190,12 +198,16 @@ def main(args):
                                              func=par['telluric']['func'],
                                              model=par['telluric']['model'],
                                              polyorder=par['telluric']['polyorder'],
+                                             teltype=par['telluric']['teltype'], tell_npca=par['telluric']['tell_npca'],
                                              fit_wv_min_max=par['telluric']['fit_wv_min_max'],
                                              mask_lyman_a=par['telluric']['mask_lyman_a'],
                                              delta_coeff_bounds=par['telluric']['delta_coeff_bounds'],
                                              minmax_coeff_bounds=par['telluric']['minmax_coeff_bounds'],
                                              only_orders=par['telluric']['only_orders'],
+                                             pix_shift_bounds=par['telluric']['pix_shift_bounds'],
                                              maxiter=par['telluric']['maxiter'],
+                                             popsize=par['telluric']['popsize'],
+                                             tol=par['telluric']['tol'],
                                              debug_init=args.debug, disp=args.debug,
                                              debug=args.debug, show=args.plot)
         else:
diff --git a/pypeit/sensfunc.py b/pypeit/sensfunc.py
index 33937efba5..df2a2985b3 100644
--- a/pypeit/sensfunc.py
+++ b/pypeit/sensfunc.py
@@ -885,7 +885,10 @@ def compute_zeropoint(self):
                                                    ech_orders=self.meta_spec['ECH_ORDERS'],
                                                    resln_guess=self.par['IR']['resln_guess'],
                                                    resln_frac_bounds=self.par['IR']['resln_frac_bounds'],
+                                                   pix_shift_bounds=self.par['IR']['pix_shift_bounds'],
                                                    sn_clip=self.par['IR']['sn_clip'],
+                                                   teltype=self.par['IR']['teltype'],
+                                                   tell_npca=self.par['IR']['tell_npca'],
                                                    mask_hydrogen_lines=self.par['mask_hydrogen_lines'],
                                                    maxiter=self.par['IR']['maxiter'],
                                                    lower=self.par['IR']['lower'],
diff --git a/pypeit/spectrographs/gemini_flamingos.py b/pypeit/spectrographs/gemini_flamingos.py
index 6873fb4631..9e9d459404 100644
--- a/pypeit/spectrographs/gemini_flamingos.py
+++ b/pypeit/spectrographs/gemini_flamingos.py
@@ -143,7 +143,7 @@ def default_pypeit_par(cls):
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8
         # TODO: replace the telluric grid file for Gemini-S site.
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_LasCampanas_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         return par
 
diff --git a/pypeit/spectrographs/gemini_gmos.py b/pypeit/spectrographs/gemini_gmos.py
index b16d279c50..1015bff45a 100644
--- a/pypeit/spectrographs/gemini_gmos.py
+++ b/pypeit/spectrographs/gemini_gmos.py
@@ -873,7 +873,7 @@ def default_pypeit_par(cls):
         """
         par = super().default_pypeit_par()
         par['sensfunc']['algorithm'] = 'IR'
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_LasCampanas_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
         # Bound the detector with slit edges if no edges are found. These data are often trimmed
         # so we implement this here as the default.
         par['calibrations']['slitedges']['bound_detector'] = True
diff --git a/pypeit/spectrographs/gemini_gnirs.py b/pypeit/spectrographs/gemini_gnirs.py
index 921ecccaea..148464d7d1 100644
--- a/pypeit/spectrographs/gemini_gnirs.py
+++ b/pypeit/spectrographs/gemini_gnirs.py
@@ -297,7 +297,7 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 6
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
         return par
 
     def config_specific_par(self, scifile, inp_par=None):
diff --git a/pypeit/spectrographs/gtc_osiris.py b/pypeit/spectrographs/gtc_osiris.py
index 1a6864f8fa..f622ad99de 100644
--- a/pypeit/spectrographs/gtc_osiris.py
+++ b/pypeit/spectrographs/gtc_osiris.py
@@ -384,7 +384,7 @@ def config_specific_par(self, scifile, inp_par=None):
             par['calibrations']['wavelengths']['lamps'] = ['ArI','XeI','NeI']
             par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500I.fits'
             par['sensfunc']['algorithm'] = 'IR'
-            par['sensfunc']['IR']['telgridfile'] = "TelFit_MaunaKea_3100_26100_R20000.fits"
+            par['sensfunc']['IR']['telgridfile'] = "TellPCA_3000_26000_R10000.fits"
         else:
             msgs.warn('gtc_osiris.py: template arc missing for this grism! Trying holy-grail...')
             par['calibrations']['wavelengths']['method'] = 'holy-grail'
@@ -955,7 +955,7 @@ def config_specific_par(self, scifile, inp_par=None):
             par['calibrations']['wavelengths']['lamps'] = ['ArI','XeI','NeI']
             par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500I.fits'
             par['sensfunc']['algorithm'] = 'IR'
-            par['sensfunc']['IR']['telgridfile'] = "TelFit_MaunaKea_3100_26100_R20000.fits"
+            par['sensfunc']['IR']['telgridfile'] = "TellPCA_3000_26000_R10000.fits"
         else:
             msgs.warn('gtc_osiris.py: template arc missing for this grism! Trying holy-grail...')
             par['calibrations']['wavelengths']['method'] = 'holy-grail'
diff --git a/pypeit/spectrographs/keck_deimos.py b/pypeit/spectrographs/keck_deimos.py
index b7c711af6c..02d8758304 100644
--- a/pypeit/spectrographs/keck_deimos.py
+++ b/pypeit/spectrographs/keck_deimos.py
@@ -301,7 +301,7 @@ def default_pypeit_par(cls):
         par['scienceframe']['process']['objlim'] = 1.5
 
         # If telluric is triggered
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R15000.fits'
         return par
 
     def config_specific_par(self, scifile, inp_par=None):
diff --git a/pypeit/spectrographs/keck_hires.py b/pypeit/spectrographs/keck_hires.py
index 2757803548..b03fd2359a 100644
--- a/pypeit/spectrographs/keck_hires.py
+++ b/pypeit/spectrographs/keck_hires.py
@@ -157,7 +157,14 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 5 #[9, 11, 11, 9, 9, 8, 8, 7, 7, 7, 7, 7, 7, 7, 7]
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_10500_R120000.fits'
+        par['sensfunc']['IR']['pix_shift_bounds'] = (-40.0,40.0)
+        
+        # Telluric parameters
+        # HIRES is usually oversampled, so the helio shift can be large
+        par['telluric']['pix_shift_bounds'] = (-40.0,40.0)
+        # Similarly, the resolution guess is higher than it should be
+        par['telluric']['resln_frac_bounds'] = (0.25,1.25)
 
         # Coadding
         par['coadd1d']['wave_method'] = 'log10'
diff --git a/pypeit/spectrographs/keck_lris.py b/pypeit/spectrographs/keck_lris.py
index 5172bd92c8..43b5f53aa5 100644
--- a/pypeit/spectrographs/keck_lris.py
+++ b/pypeit/spectrographs/keck_lris.py
@@ -101,7 +101,7 @@ def default_pypeit_par(cls):
 
 
         # If telluric is triggered
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         return par
 
diff --git a/pypeit/spectrographs/keck_mosfire.py b/pypeit/spectrographs/keck_mosfire.py
index 9059901347..e874baea83 100644
--- a/pypeit/spectrographs/keck_mosfire.py
+++ b/pypeit/spectrographs/keck_mosfire.py
@@ -127,7 +127,7 @@ def default_pypeit_par(cls):
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 13
         par['sensfunc']['IR']['maxiter'] = 2
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
         return par
 
     # NOTE: This function is used by the dev-suite
diff --git a/pypeit/spectrographs/keck_nires.py b/pypeit/spectrographs/keck_nires.py
index bd4b3b1138..761fcb7253 100644
--- a/pypeit/spectrographs/keck_nires.py
+++ b/pypeit/spectrographs/keck_nires.py
@@ -140,7 +140,7 @@ def default_pypeit_par(cls):
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8
         par['sensfunc']['IR']['maxiter'] = 2
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         # Coadding
         par['coadd1d']['wave_method'] = 'log10'
diff --git a/pypeit/spectrographs/keck_nirspec.py b/pypeit/spectrographs/keck_nirspec.py
index eab7cafc91..60e84e0df7 100644
--- a/pypeit/spectrographs/keck_nirspec.py
+++ b/pypeit/spectrographs/keck_nirspec.py
@@ -129,7 +129,12 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R25000.fits'
+        par['sensfunc']['IR']['pix_shift_bounds'] = (-8.0,8.0)
+        
+        # Telluric parameters
+        par['telluric']['pix_shift_bounds'] = (-8.0,8.0)
+        
         return par
 
     def init_meta(self):
diff --git a/pypeit/spectrographs/magellan_fire.py b/pypeit/spectrographs/magellan_fire.py
index 2762b10d86..95d4dcb617 100644
--- a/pypeit/spectrographs/magellan_fire.py
+++ b/pypeit/spectrographs/magellan_fire.py
@@ -201,7 +201,7 @@ def default_pypeit_par(cls):
         par['sensfunc']['polyorder'] = 5
         par['sensfunc']['IR']['maxiter'] = 2
         # place holder for telgrid file
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_LasCampanas_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R15000.fits'
 
         # Coadding. I'm not sure what this should be for PRISM mode?
         par['coadd1d']['wave_method'] = 'log10'
@@ -418,7 +418,7 @@ def default_pypeit_par(cls):
         par['reduce']['findobj']['find_trim_edge'] = [50,50]
         par['flexure']['spec_method'] = 'skip'
 
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_LasCampanas_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         # Set the default exposure time ranges for the frame typing
         par['calibrations']['standardframe']['exprng'] = [None, 60]
diff --git a/pypeit/spectrographs/mmt_binospec.py b/pypeit/spectrographs/mmt_binospec.py
index 1b3998f1f4..90e178b658 100644
--- a/pypeit/spectrographs/mmt_binospec.py
+++ b/pypeit/spectrographs/mmt_binospec.py
@@ -214,7 +214,7 @@ def default_pypeit_par(cls):
 
         # Sensitivity function parameters
         par['sensfunc']['polyorder'] = 7
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         return par
 
diff --git a/pypeit/spectrographs/mmt_mmirs.py b/pypeit/spectrographs/mmt_mmirs.py
index d9341cd176..4425543c29 100644
--- a/pypeit/spectrographs/mmt_mmirs.py
+++ b/pypeit/spectrographs/mmt_mmirs.py
@@ -202,7 +202,7 @@ def default_pypeit_par(cls):
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8
         # ToDo: replace the telluric grid file for MMT site.
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         return par
 
diff --git a/pypeit/spectrographs/p200_dbsp.py b/pypeit/spectrographs/p200_dbsp.py
index caecea72f0..0ea3032236 100644
--- a/pypeit/spectrographs/p200_dbsp.py
+++ b/pypeit/spectrographs/p200_dbsp.py
@@ -525,7 +525,7 @@ def default_pypeit_par(cls):
 
         par['sensfunc']['algorithm'] = 'UVIS'
         par['sensfunc']['UVIS']['polycorrect'] = False
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_Lick_3100_11100_R10000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
         return par
 
     def config_specific_par(self, scifile, inp_par=None):
diff --git a/pypeit/spectrographs/p200_tspec.py b/pypeit/spectrographs/p200_tspec.py
index 1d4000d047..641dc2d60e 100644
--- a/pypeit/spectrographs/p200_tspec.py
+++ b/pypeit/spectrographs/p200_tspec.py
@@ -215,7 +215,7 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_MaunaKea_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         # Coadding
         par['coadd1d']['wave_method'] = 'log10'
diff --git a/pypeit/spectrographs/shane_kast.py b/pypeit/spectrographs/shane_kast.py
index e58691757f..0f667a44c0 100644
--- a/pypeit/spectrographs/shane_kast.py
+++ b/pypeit/spectrographs/shane_kast.py
@@ -57,7 +57,7 @@ def default_pypeit_par(cls):
         par['calibrations']['standardframe']['exprng'] = [1, 61]
         #
         par['scienceframe']['exprng'] = [61, None]
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_Lick_3100_11100_R10000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
         return par
 
     def init_meta(self):
@@ -460,7 +460,7 @@ def default_pypeit_par(cls):
 
         # TODO In case someone wants to use the IR algorithm for shane kast this is the telluric file. Note the IR
         # algorithm is not the default.
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_Lick_3100_11100_R10000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
         return par
 
     def init_meta(self):
@@ -684,7 +684,7 @@ def default_pypeit_par(cls):
         par['calibrations']['wavelengths']['rms_thresh_frac_fwhm'] = 0.09
         par['calibrations']['wavelengths']['sigdetect'] = 5.
         par['calibrations']['wavelengths']['use_instr_flag'] = True
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_Lick_3100_11100_R10000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
 
         return par
diff --git a/pypeit/spectrographs/soar_goodman.py b/pypeit/spectrographs/soar_goodman.py
index 1f1f9c5bc4..293bcaf025 100644
--- a/pypeit/spectrographs/soar_goodman.py
+++ b/pypeit/spectrographs/soar_goodman.py
@@ -330,7 +330,7 @@ def default_pypeit_par(cls):
         par['scienceframe']['exprng'] = [90, None]
 
         #par['sensfunc']['algorithm'] = 'IR'
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_LasCampanas_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R15000.fits'
 
         # TODO: Temporary fix for failure mode.  Remove once Ryan provides a
         # fix.
@@ -527,7 +527,7 @@ def default_pypeit_par(cls):
         par['scienceframe']['exprng'] = [90, None]
 
         # par['sensfunc']['algorithm'] = 'IR'
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_LasCampanas_3100_26100_R20000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R15000.fits'
 
         return par
 
diff --git a/pypeit/spectrographs/vlt_fors.py b/pypeit/spectrographs/vlt_fors.py
index 254245ff44..aff71e4f86 100644
--- a/pypeit/spectrographs/vlt_fors.py
+++ b/pypeit/spectrographs/vlt_fors.py
@@ -69,7 +69,7 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 5
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_Paranal_VIS_9800_25000_R25000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
 
 
@@ -457,4 +457,4 @@ def parse_dither_pattern(self, file_list, ext=None):
 
 #            dither_id.append(hdr['FRAMEID'])
 #            offset_arcsec[ifile] = hdr['YOFFSET']
-        return dither_pattern, dither_id, offset_arcsec
\ No newline at end of file
+        return dither_pattern, dither_id, offset_arcsec
diff --git a/pypeit/spectrographs/vlt_sinfoni.py b/pypeit/spectrographs/vlt_sinfoni.py
index 68f40af6b3..df29ac8653 100644
--- a/pypeit/spectrographs/vlt_sinfoni.py
+++ b/pypeit/spectrographs/vlt_sinfoni.py
@@ -157,7 +157,7 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 7
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_Paranal_NIR_9800_25000_R25000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
 
         return par
 
diff --git a/pypeit/spectrographs/vlt_xshooter.py b/pypeit/spectrographs/vlt_xshooter.py
index b2ffa58ceb..01cc8fc088 100644
--- a/pypeit/spectrographs/vlt_xshooter.py
+++ b/pypeit/spectrographs/vlt_xshooter.py
@@ -344,7 +344,12 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_Paranal_NIR_9800_25000_R25000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R25000.fits'
+        par['sensfunc']['IR']['pix_shift_bounds'] = (-10.0,10.0)
+        
+        # Telluric parameters
+        par['telluric']['pix_shift_bounds'] = (-10.0,10.0)
+        par['telluric']['resln_frac_bounds'] = (0.4,2.0)
 
         # Coadding
         par['coadd1d']['wave_method'] = 'log10'
@@ -719,7 +724,12 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] = 8 #[9, 11, 11, 9, 9, 8, 8, 7, 7, 7, 7, 7, 7, 7, 7]
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_Paranal_VIS_4900_11100_R25000.fits'
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R25000.fits'
+        par['sensfunc']['IR']['pix_shift_bounds'] = (-10.0,10.0)
+        
+        # Telluric parameters
+        par['telluric']['pix_shift_bounds'] = (-10.0,10.0)
+        par['telluric']['resln_frac_bounds'] = (0.4,2.0)
 
         # Coadding
         par['coadd1d']['wave_method'] = 'log10'
@@ -1007,8 +1017,11 @@ def default_pypeit_par(cls):
         # Sensitivity function parameters
         par['sensfunc']['algorithm'] = 'IR'
         par['sensfunc']['polyorder'] =  8
-        par['sensfunc']['IR']['telgridfile'] = 'TelFit_LasCampanas_3100_26100_R20000.fits'
-        # This is a hack until we we have a Paranal file generated that covers the UVB wavelength range.
+        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R25000.fits'
+        par['sensfunc']['IR']['pix_shift_bounds'] = (-8.0,8.0)
+        
+        # Telluric parameters
+        par['telluric']['pix_shift_bounds'] = (-8.0,8.0)
 
         # Coadding
         par['coadd1d']['wave_method'] = 'log10'
diff --git a/pypeit/tests/tstutils.py b/pypeit/tests/tstutils.py
index a6393a2c15..c120e2f388 100644
--- a/pypeit/tests/tstutils.py
+++ b/pypeit/tests/tstutils.py
@@ -36,7 +36,8 @@
 
 # Tests require the Telluric file (Mauna Kea)
 par = Spectrograph.default_pypeit_par()
-tell_test_grid = data.Paths.telgrid / 'TelFit_MaunaKea_3100_26100_R20000.fits'
+tell_test_grid = data.get_telgrid_filepath('TellPCA_3000_26000_R25000.fits')
+#tell_test_grid = data.Paths.telgrid / 'TelFit_MaunaKea_3100_26100_R20000.fits'
 telluric_required = pytest.mark.skipif(not tell_test_grid.is_file(),
                                        reason='no Mauna Kea telluric file')