diff --git a/doc/changes.rst b/doc/changes.rst
index 7049c5f73..7b5d71378 100644
--- a/doc/changes.rst
+++ b/doc/changes.rst
@@ -7,8 +7,11 @@ desisim change log
 
 * Major refactor of newexp to add connection to upstream mocks, surveysims,
   and fiber assignment (`PR #250`_).
+* Support latest (>DR4) data model in the templates metadata table and also
+  scale simulated templates by 1e17 erg/s/cm2/Angstrom (`PR #252`_).
 
 .. _`PR #250`: https://github.com/desihub/desisim/pull/250
+.. _`PR #252`: https://github.com/desihub/desisim/pull/252
 
 0.20.0 (2017-07-12)
 -------------------
diff --git a/py/desisim/io.py b/py/desisim/io.py
index d3a9bd513..fa27e5c3f 100644
--- a/py/desisim/io.py
+++ b/py/desisim/io.py
@@ -863,7 +863,7 @@ def write_templates(outfile, flux, wave, meta, objtype=None,
         CDELT1 = (wave[1]-wave[0], 'Wavelength step [Angstrom]'),
         LOGLAM = (0, 'linear wavelength steps, not log10'),
         AIRORVAC = ('vac', 'wavelengths in vacuum (vac) or air'),
-        BUNIT = ('erg/s/cm2/A', 'spectrum flux units')
+        BUNIT = ('1e-17 erg/(s cm2 Angstrom)', 'spectrum flux units')
         )
     hdr = fitsheader(header)
 
@@ -916,9 +916,17 @@ def empty_metatable(nmodel=1, objtype='ELG', subtype='', add_SNeIa=None):
     meta.add_column(Column(name='REDSHIFT', length=nmodel, dtype='f4',
                            data=np.zeros(nmodel)))
     meta.add_column(Column(name='MAG', length=nmodel, dtype='f4',
-                           data=np.zeros(nmodel)-1))
-    meta.add_column(Column(name='DECAM_FLUX', shape=(6,), length=nmodel, dtype='f4'))
-    meta.add_column(Column(name='WISE_FLUX', shape=(2,), length=nmodel, dtype='f4'))
+                           data=np.zeros(nmodel)-1, unit='mag'))
+    meta.add_column(Column(name='FLUX_G', length=nmodel, dtype='f4',
+                           unit='nanomaggies'))
+    meta.add_column(Column(name='FLUX_R', length=nmodel, dtype='f4',
+                           unit='nanomaggies'))
+    meta.add_column(Column(name='FLUX_Z', length=nmodel, dtype='f4',
+                           unit='nanomaggies'))
+    meta.add_column(Column(name='FLUX_W1', length=nmodel, dtype='f4',
+                           unit='nanomaggies'))
+    meta.add_column(Column(name='FLUX_W2', length=nmodel, dtype='f4',
+                           unit='nanomaggies'))
 
     meta.add_column(Column(name='OIIFLUX', length=nmodel, dtype='f4',
                            data=np.zeros(nmodel)-1, unit='erg/(s*cm2)'))
diff --git a/py/desisim/lya_spectra.py b/py/desisim/lya_spectra.py
index 279221298..f2c6cf05b 100644
--- a/py/desisim/lya_spectra.py
+++ b/py/desisim/lya_spectra.py
@@ -95,6 +95,7 @@ def get_spectra(lyafile, nqso=None, wave=None, templateid=None, normfilter='sdss
         flux1, _, meta1 = qso.make_templates(nmodel=1, redshift=np.atleast_1d(zqso[ii]), 
                                              mag=np.atleast_1d(mag_g[ii]), seed=templateseed[ii],
                                              nocolorcuts=nocolorcuts)
+        flux1 *= 1e-17
         for col in meta1.colnames:
             meta[col][ii] = meta1[col][0]
 
@@ -114,8 +115,8 @@ def get_spectra(lyafile, nqso=None, wave=None, templateid=None, normfilter='sdss
                                                        mask_invalid=True)[normfilter])
         factor = 10**(-0.4 * mag_g[ii]) / normmaggies
         flux1 *= factor
-        meta['DECAM_FLUX'][ii] *= factor
-        meta['WISE_FLUX'][ii] *= factor
+        for key in ('FLUX_G', 'FLUX_R', 'FLUX_Z', 'FLUX_W1', 'FLUX_W2'):
+            meta[key][ii] *= factor
 
         flux[ii, :] = flux1[:]
 
diff --git a/py/desisim/quickcat.py b/py/desisim/quickcat.py
index 972a138f3..5780386c8 100644
--- a/py/desisim/quickcat.py
+++ b/py/desisim/quickcat.py
@@ -160,16 +160,26 @@ def get_redshift_efficiency(simtype, targets, truth, targets_in_tile, obsconditi
     targetid = targets['TARGETID']
     n = len(targetid)
 
+    try:
+        if 'DECAM_FLUX' in targets.colnames:
+            true_gflux = targets['DECAM_FLUX'][:, 1]
+            true_rflux = targets['DECAM_FLUX'][:, 2]
+        else:
+            true_gflux = targets['FLUX_G']
+            true_rflux = targets['FLUX_R']
+    except:
+        raise Exception('Missing photometry needed to estimate redshift efficiency!')
+
+    if (obsconditions is None) and (truth['OIIFLUX'] not in truth.colnames):
+        raise Exception('Missing obsconditions and flux information to estimate redshift efficiency')
+
     if (simtype == 'ELG'):
         # Read the model OII flux threshold (FDR fig 7.12 modified to fit redmonster efficiency on OAK)
         filename = resource_filename('desisim', 'data/quickcat_elg_oii_flux_threshold.txt')
         fdr_z, modified_fdr_oii_flux_threshold = np.loadtxt(filename, unpack=True)
 
         # Get OIIflux from truth
-        try:
-            true_oii_flux = truth['OIIFLUX']
-        except:
-            raise Exception('Missing OII flux information to estimate redshift efficiency for ELGs')
+        true_oii_flux = truth['OIIFLUX']
 
         # Compute OII flux thresholds for truez
         oii_flux_threshold = np.interp(truth['TRUEZ'],fdr_z,modified_fdr_oii_flux_threshold)
@@ -186,11 +196,6 @@ def get_redshift_efficiency(simtype, targets, truth, targets_in_tile, obsconditi
         magr, magr_eff = np.loadtxt(filename, unpack=True)
 
         # Get Rflux from truth
-        try:
-            true_rflux = targets['DECAM_FLUX'][:,2]
-        except:
-            raise Exception('Missing Rmag information to estimate redshift efficiency for LRGs')
-
         r_mag = 22.5 - 2.5*np.log10(true_rflux)
 
         mean_eff_mag=np.interp(r_mag,magr,magr_eff)
@@ -205,11 +210,8 @@ def get_redshift_efficiency(simtype, targets, truth, targets_in_tile, obsconditi
         zc, qso_gmag_threshold_vs_z = np.loadtxt(filename, unpack=True)
 
         # Get Gflux from truth
-        try:
-            true_gmag = targets['DECAM_FLUX'][:,1]
-        except:
-            raise Exception('Missing Gmag information to estimate redshift efficiency for QSOs')
-
+        true_gmag = 22.5 - 2.5 * np.log10(true_gflux) 
+            
         # Computes QSO mag thresholds for truez
         qso_gmag_threshold=np.interp(truth['TRUEZ'],zc,qso_gmag_threshold_vs_z)
         assert (qso_gmag_threshold.size == true_gmag.size),"qso_gmag_threshold and true_gmag should have the same size"
@@ -233,9 +235,6 @@ def get_redshift_efficiency(simtype, targets, truth, targets_in_tile, obsconditi
         log.warning('using default redshift efficiency of {} for {}'.format(default_zeff, simtype))
         simulated_eff = default_zeff * np.ones(n)
 
-    if (obsconditions is None) and (truth['OIIFLUX'] is None) and (targets['DECAM_FLUX'] is None):
-        raise Exception('Missing obsconditions and flux information to estimate redshift efficiency')
-
     #- Get the corrections for observing conditions per tile, then
     #- correct targets on those tiles.  Parameterize in terms of failure
     #- rate instead of success rate to handle bookkeeping of targets that
@@ -343,6 +342,16 @@ def get_observed_redshifts(targets, truth, targets_in_tile, obsconditions):
     truez = truth['TRUEZ']
     targetid = truth['TARGETID']
 
+    try:
+        if 'DECAM_FLUX' in targets.colnames:
+            true_gflux = targets['DECAM_FLUX'][:, 1]
+            true_rflux = targets['DECAM_FLUX'][:, 2]
+        else:
+            true_gflux = targets['FLUX_G']
+            true_rflux = targets['FLUX_R']
+    except:
+        raise Exception('Missing photometry needed to estimate redshift efficiency!')
+
     zout = truez.copy()
     zerr = np.zeros(len(truez), dtype=np.float32)
     zwarn = np.zeros(len(truez), dtype=np.int32)
@@ -375,10 +384,8 @@ def get_observed_redshifts(targets, truth, targets_in_tile, obsconditions):
             elif (objtype == 'LRG'):
                 redbins=np.linspace(0.55,1.1,5)
                 mag = np.linspace(20.5,23.,50)
-                try:
-                    true_magr=targets['DECAM_FLUX'][ii,2]
-                except:
-                    raise Exception('Missing DECAM r flux information to estimate redshift error for LRGs')
+
+                true_magr = 22.5 - 2.5 * np.log10(true_rflux)
 
                 coefs = [[9.46882282e-05,-1.87022383e-03],[6.14601021e-05,-1.17406643e-03],\
                              [8.85342362e-05,-1.76079966e-03],[6.96202042e-05,-1.38632104e-03]]
@@ -419,10 +426,7 @@ def get_observed_redshifts(targets, truth, targets_in_tile, obsconditions):
                 redbins = np.linspace(0.5,3.5,7)
                 mag = np.linspace(21.,23.,50)
 
-                try:
-                    true_magg=targets['DECAM_FLUX'][ii,1]
-                except:
-                    raise Exception('Missing DECAM g flux information to estimate redshift error for QSOs')
+                true_magg = 22.5 - 2.5 * np.log10(true_gflux)
 
                 coefs = [[0.000156950059747,-0.00320719603886],[0.000461779391179,-0.00924485142818],\
                              [0.000458672517009,-0.0091254038977],[0.000461427968475,-0.00923812594293],\
diff --git a/py/desisim/scripts/quickgen.py b/py/desisim/scripts/quickgen.py
index f4a2972d7..43c3143e0 100644
--- a/py/desisim/scripts/quickgen.py
+++ b/py/desisim/scripts/quickgen.py
@@ -386,7 +386,7 @@ def main(args):
             jj.append(ii.tolist())
 
             # Sanity check on units; templates currently return ergs, not 1e-17 ergs...
-            assert (thisobj == 'SKY') or (np.max(truth['FLUX']) < 1e-6)
+            # assert (thisobj == 'SKY') or (np.max(truth['FLUX']) < 1e-6)
 
         # Sort the metadata table.
         jj = sum(jj,[])
diff --git a/py/desisim/targets.py b/py/desisim/targets.py
index 364c38e85..0a0b3c94a 100644
--- a/py/desisim/targets.py
+++ b/py/desisim/targets.py
@@ -335,17 +335,15 @@ def get_targets(nspec, program, tileid=None, seed=None, specmin=0):
         else:
             raise ValueError('Unable to simulate OBJTYPE={}'.format(objtype))
 
-        flux[ii] = 1e17 * simflux
+        flux[ii] = simflux
         meta[ii] = meta1
 
-        ugrizy = 22.5-2.5*np.log10(meta1['DECAM_FLUX'].data)
-        wise = 22.5-2.5*np.log10(meta1['WISE_FLUX'].data)
         fibermap['FILTER'][ii, :6] = ['DECAM_G', 'DECAM_R', 'DECAM_Z', 'WISE_W1', 'WISE_W2']
-        fibermap['MAG'][ii, 0] = ugrizy[:, 1]
-        fibermap['MAG'][ii, 1] = ugrizy[:, 2]
-        fibermap['MAG'][ii, 2] = ugrizy[:, 4]
-        fibermap['MAG'][ii, 3] = wise[:, 0]
-        fibermap['MAG'][ii, 4] = wise[:, 1]
+        fibermap['MAG'][ii, 0] = 22.5 - 2.5 * np.log10(meta['FLUX_G'][ii])
+        fibermap['MAG'][ii, 1] = 22.5 - 2.5 * np.log10(meta['FLUX_R'][ii])
+        fibermap['MAG'][ii, 2] = 22.5 - 2.5 * np.log10(meta['FLUX_Z'][ii])
+        fibermap['MAG'][ii, 3] = 22.5 - 2.5 * np.log10(meta['FLUX_W1'][ii])
+        fibermap['MAG'][ii, 4] = 22.5 - 2.5 * np.log10(meta['FLUX_W2'][ii])
 
     #- Load fiber -> positioner mapping and tile information
     fiberpos = desimodel.io.load_fiberpos()
diff --git a/py/desisim/templates.py b/py/desisim/templates.py
index fc3bd61cf..d9982a95a 100755
--- a/py/desisim/templates.py
+++ b/py/desisim/templates.py
@@ -11,6 +11,8 @@
 import sys
 import numpy as np
 
+import astropy.units as u
+
 from desisim.io import empty_metatable
 
 LIGHT = 2.99792458E5  #- speed of light in km/s
@@ -283,9 +285,10 @@ def __init__(self, objtype='ELG', minwave=3600.0, maxwave=10000.0, cdelt=0.2,
             corresponding to BASEFLUX (Angstrom).
           basemeta (astropy.Table): Table of meta-data [nbase] for each base template.
           pixbound (numpy.ndarray): Pixel boundaries of BASEWAVE (Angstrom).
-          rfilt (speclite.filters instance): DECam2014 r-band FilterSequence
-          normfilt (speclite.filters instance): FilterSequence of self.normfilter
-          decamwise (speclite.filters instance): DECam2014-* and WISE2010-* FilterSequence
+          rfilt (speclite.filters instance): DECam2014 r-band FilterSequence.
+          normfilt (speclite.filters instance): FilterSequence of self.normfilter.
+          decamwise (speclite.filters instance): DECam2014-[g,r,z] and WISE2010-[W1,W2]
+            FilterSequence.
 
         Optional Attributes:
           sne_baseflux (numpy.ndarray): Array [sne_nbase,sne_npix] of the base
@@ -341,7 +344,8 @@ def __init__(self, objtype='ELG', minwave=3600.0, maxwave=10000.0, cdelt=0.2,
         # Initialize the filter profiles.
         self.rfilt = filters.load_filters('decam2014-r')
         self.normfilt = filters.load_filters(self.normfilter)
-        self.decamwise = filters.load_filters('decam2014-*', 'wise2010-W1', 'wise2010-W2')
+        self.decamwise = filters.load_filters('decam2014-g', 'decam2014-r', 'decam2014-z',
+                                              'wise2010-W1', 'wise2010-W2')
 
     def _blurmatrix(self, vdisp, log=None):
         """Pre-compute the blur_matrix as a dictionary keyed by each unique value of
@@ -479,10 +483,12 @@ def make_galaxy_templates(self, nmodel=100, zrange=(0.6, 1.6), magrange=(21.0, 2
             templates instead of resampled observer frame.
           verbose (bool, optional): Be verbose!
 
-        Returns:
-          outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame spectra (erg/s/cm2/A).
-          wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
-          meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
+        Returns (outflux, wave, meta) tuple where:
+
+          * outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame
+            spectra (1e-17 erg/s/cm2/A).
+          * wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
+          * meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
 
         Raises:
           ValueError
@@ -695,20 +701,20 @@ def make_galaxy_templates(self, nmodel=100, zrange=(0.6, 1.6), magrange=(21.0, 2
                             restflux, zwave, mask_invalid=True)[self.normfilter])
                     magnorm = 10**(-0.4*mag[ii]) / normmaggies
 
-                synthnano = np.zeros((nbasechunk, len(self.decamwise)))
-                for ff, key in enumerate(maggies.columns):
-                    synthnano[:, ff] = 1E9 * maggies[key] * magnorm # nanomaggies
+                synthnano = dict()
+                for key in maggies.columns:
+                    synthnano[key] = 1E9 * maggies[key] * magnorm # nanomaggies
                 zlineflux = normlineflux[templateid] * magnorm
 
                 if nocolorcuts or self.colorcuts_function is None:
                     colormask = np.repeat(1, nbasechunk)
                 else:
                     colormask = self.colorcuts_function(
-                        gflux=synthnano[:, 1],
-                        rflux=synthnano[:, 2],
-                        zflux=synthnano[:, 4],
-                        w1flux=synthnano[:, 6],
-                        w2flux=synthnano[:, 7])
+                        gflux=synthnano['decam2014-g'],
+                        rflux=synthnano['decam2014-r'],
+                        zflux=synthnano['decam2014-z'],
+                        w1flux=synthnano['wise2010-W1'],
+                        w2flux=synthnano['wise2010-W2'])
 
                 # If the color-cuts pass then populate the output flux vector
                 # (suitably normalized) and metadata table, convolve with the
@@ -733,8 +739,11 @@ def make_galaxy_templates(self, nmodel=100, zrange=(0.6, 1.6), magrange=(21.0, 2
 
                     meta['TEMPLATEID'][ii] = tempid
                     meta['D4000'][ii] = d4000[tempid]
-                    meta['DECAM_FLUX'][ii] = synthnano[this, :6]
-                    meta['WISE_FLUX'][ii] = synthnano[this, 6:8]
+                    meta['FLUX_G'][ii] = synthnano['decam2014-g'][this]
+                    meta['FLUX_R'][ii] = synthnano['decam2014-r'][this]
+                    meta['FLUX_Z'][ii] = synthnano['decam2014-z'][this]
+                    meta['FLUX_W1'][ii] = synthnano['wise2010-W1'][this]
+                    meta['FLUX_W2'][ii] = synthnano['wise2010-W2'][this]
 
                     if self.normline is not None:
                         if self.normline == 'OII':
@@ -753,9 +762,9 @@ def make_galaxy_templates(self, nmodel=100, zrange=(0.6, 1.6), magrange=(21.0, 2
                         format(np.sum(success == 0)))
 
         if restframe:
-            return outflux, self.basewave, meta
+            return 1e17 * outflux, self.basewave, meta
         else:
-            return outflux, self.wave, meta
+            return 1e17 * outflux, self.wave, meta
 
 class ELG(GALAXY):
     """Generate Monte Carlo spectra of emission-line galaxies (ELGs)."""
@@ -815,10 +824,12 @@ def make_templates(self, nmodel=100, zrange=(0.6, 1.6), rmagrange=(21.0, 23.4),
           minoiiflux (float, optional): Minimum [OII] 3727 flux (default 0.0
             erg/s/cm2).
 
-        Returns:
-          outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame spectra (erg/s/cm2/A).
-          wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
-          meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
+        Returns (outflux, wave, meta) tuple where:
+
+          * outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame
+            spectra (1e-17 erg/s/cm2/A).
+          * wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
+          * meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
 
         Raises:
 
@@ -896,15 +907,13 @@ def make_templates(self, nmodel=100, zrange=(0.01, 0.4), rmagrange=(15.0, 19.5),
           minhbetaflux (float, optional): Minimum H-beta flux (default 0.0
             erg/s/cm2).
 
-        Returns:
-          outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame spectra (erg/s/cm2/A).
-          wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
-          meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
-
-        Raises:
+        Returns (outflux, wave, meta) tuple where:
 
+          * outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame
+            spectra (1e-17 erg/s/cm2/A).
+          * wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
+          * meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
         """
-
         outflux, wave, meta = self.make_galaxy_templates(nmodel=nmodel, zrange=zrange, magrange=rmagrange,
                                                          oiiihbrange=oiiihbrange, logvdisp_meansig=logvdisp_meansig,
                                                          minlineflux=minhbetaflux, redshift=redshift, vdisp=vdisp,
@@ -965,10 +974,12 @@ def make_templates(self, nmodel=100, zrange=(0.5, 1.0), zmagrange=(19.0, 20.4),
           agnlike (bool, optional): adopt AGN-like emission-line ratios (not yet
             supported; defaults False).
 
-        Returns:
-          outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame spectra (erg/s/cm2/A).
-          wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
-          meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
+        Returns (outflux, wave, meta) tuple where:
+
+          * outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame
+            spectra (1e-17 erg/s/cm2/A).
+          * wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
+          * meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
 
         Raises:
 
@@ -1026,8 +1037,9 @@ def __init__(self, objtype='STAR', subtype='', minwave=3600.0, maxwave=10000.0,
           basewave (numpy.ndarray): Array [npix] of rest-frame wavelengths
             corresponding to BASEFLUX (Angstrom).
           basemeta (astropy.Table): Table of meta-data [nbase] for each base template.
-          normfilt (speclite.filters instance): FilterSequence of self.normfilter
-          decamwise (speclite.filters instance): DECam2014-* and WISE2010-* FilterSequence
+          normfilt (speclite.filters instance): FilterSequence of self.normfilter.
+          decamwise (speclite.filters instance): DECam2014-[g,r,z] and WISE2010-[W1,W2]
+            FilterSequence.
 
         """
         from speclite import filters
@@ -1055,7 +1067,8 @@ def __init__(self, objtype='STAR', subtype='', minwave=3600.0, maxwave=10000.0,
 
         # Initialize the filter profiles.
         self.normfilt = filters.load_filters(self.normfilter)
-        self.decamwise = filters.load_filters('decam2014-*', 'wise2010-W1', 'wise2010-W2')
+        self.decamwise = filters.load_filters('decam2014-g', 'decam2014-r', 'decam2014-z',
+                                              'wise2010-W1', 'wise2010-W2')
 
     def make_star_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0),
                             magrange=(18.0, 23.5), seed=None, redshift=None,
@@ -1125,10 +1138,12 @@ def make_star_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0),
             templates instead of resampled observer frame.
           verbose (bool, optional): Be verbose!
 
-        Returns:
-          outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame spectra (erg/s/cm2/A).
-          wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
-          meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
+        Returns (outflux, wave, meta) tuple where:
+
+          * outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame
+            spectra (1e-17 erg/s/cm2/A).
+          * wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
+          * meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
 
         Raises:
           ValueError
@@ -1262,19 +1277,19 @@ def make_star_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0),
                         restflux, zwave, mask_invalid=True)[self.normfilter])
                 magnorm = 10**(-0.4*mag[ii]) / normmaggies
 
-                synthnano = np.zeros((nbasechunk, len(self.decamwise)))
-                for ff, key in enumerate(maggies.columns):
-                    synthnano[:, ff] = 1E9 * maggies[key] * magnorm
+                synthnano = dict()
+                for key in maggies.columns:
+                    synthnano[key] = 1E9 * maggies[key] * magnorm
 
                 if nocolorcuts or self.colorcuts_function is None:
                     colormask = np.repeat(1, nbasechunk)
                 else:
                     colormask = self.colorcuts_function(
-                        gflux=synthnano[:, 1],
-                        rflux=synthnano[:, 2],
-                        zflux=synthnano[:, 4],
-                        w1flux=synthnano[:, 6],
-                        w2flux=synthnano[:, 7])
+                        gflux=synthnano['decam2014-g'],
+                        rflux=synthnano['decam2014-r'],
+                        zflux=synthnano['decam2014-z'],
+                        w1flux=synthnano['wise2010-W1'],
+                        w2flux=synthnano['wise2010-W2'])
 
                 # If the color-cuts pass then populate the output flux vector
                 # (suitably normalized) and metadata table and finish up.
@@ -1291,8 +1306,11 @@ def make_star_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0),
                                                        extrapolate=True) * magnorm[this]
 
                     meta['TEMPLATEID'][ii] = tempid
-                    meta['DECAM_FLUX'][ii] = synthnano[this, :6]
-                    meta['WISE_FLUX'][ii] = synthnano[this, 6:8]
+                    meta['FLUX_G'][ii] = synthnano['decam2014-g'][this]
+                    meta['FLUX_R'][ii] = synthnano['decam2014-r'][this]
+                    meta['FLUX_Z'][ii] = synthnano['decam2014-z'][this]
+                    meta['FLUX_W1'][ii] = synthnano['wise2010-W1'][this]
+                    meta['FLUX_W2'][ii] = synthnano['wise2010-W2'][this]
 
                     if star_properties is None:
                         meta['TEFF'][ii] = self.basemeta['TEFF'][tempid]
@@ -1314,9 +1332,9 @@ def make_star_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0),
                         format(np.sum(success == 0)))
 
         if restframe:
-            return outflux, self.basewave, meta
+            return 1e17 * outflux, self.basewave, meta
         else:
-            return outflux, self.wave, meta
+            return 1e17 * outflux, self.wave, meta
 
 class STAR(SUPERSTAR):
     """Generate Monte Carlo spectra of generic stars."""
@@ -1358,10 +1376,12 @@ def make_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0),
             magnitude range.  Defaults to a uniform distribution between (18,
             23.5).
 
-        Returns:
-          outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame spectra (erg/s/cm2/A).
-          wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
-          meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
+        Returns (outflux, wave, meta) tuple where:
+
+          * outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame
+            spectra (1e-17 erg/s/cm2/A).
+          * wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
+          * meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
 
         Raises:
 
@@ -1417,10 +1437,12 @@ def make_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0), rmagrange=(16.0,
             magnitude range.  Defaults to a uniform distribution between (16,
             19).
 
-        Returns:
-          outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame spectra (erg/s/cm2/A).
-          wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
-          meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
+        Returns (outflux, wave, meta) tuple where:
+
+          * outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame
+            spectra (1e-17 erg/s/cm2/A).
+          * wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
+          * meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
 
         Raises:
 
@@ -1476,10 +1498,12 @@ def make_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0), rmagrange=(16.0,
             magnitude range.  Defaults to a uniform distribution between (16,
             20).
 
-        Returns:
-          outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame spectra (erg/s/cm2/A).
-          wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
-          meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
+        Returns (outflux, wave, meta) tuple where:
+
+          * outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame
+            spectra (1e-17 erg/s/cm2/A).
+          * wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
+          * meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
 
         Raises:
 
@@ -1532,10 +1556,12 @@ def make_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0), gmagrange=(16.0,
             magnitude range.  Defaults to a uniform distribution between (16,
             19).
 
-        Returns:
-          outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame spectra (erg/s/cm2/A).
-          wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
-          meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
+        Returns (outflux, wave, meta) tuple where:
+
+          * outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame
+            spectra (1e-17 erg/s/cm2/A).
+          * wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
+          * meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
 
         Raises:
           ValueError: If the INPUT_META or STAR_PROPERTIES table contains
@@ -1595,8 +1621,9 @@ def __init__(self, minwave=3600.0, maxwave=10000.0, cdelt=0.2, wave=None,
           wave (numpy.ndarray): Output wavelength array [Angstrom].
           cosmo (astropy.cosmology): Default cosmology object (currently
             hard-coded to FlatLCDM with H0=70, Omega0=0.3).
-          normfilt (speclite.filters instance): FilterSequence of self.normfilter
-          decamwise (speclite.filters instance): DECam2014-* and WISE2010-* FilterSequence
+          normfilt (speclite.filters instance): FilterSequence of self.normfilter.
+          decamwise (speclite.filters instance): DECam2014-[g,r,z] and WISE2010-[W1,W2]
+            FilterSequence.
 
         """
         from astropy.io import fits
@@ -1652,7 +1679,8 @@ def __init__(self, minwave=3600.0, maxwave=10000.0, cdelt=0.2, wave=None,
 
         # Initialize the filter profiles.
         self.normfilt = filters.load_filters(self.normfilter)
-        self.decamwise = filters.load_filters('decam2014-*', 'wise2010-W1', 'wise2010-W2')
+        self.decamwise = filters.load_filters('decam2014-g', 'decam2014-r', 'decam2014-z',
+                                              'wise2010-W1', 'wise2010-W2')
 
     def _sample_pcacoeff(self, nsample, coeff, rand):
         """Draw from the distribution of PCA coefficients."""
@@ -1713,10 +1741,13 @@ def make_templates(self, nmodel=100, zrange=(0.5, 4.0), rmagrange=(20.0, 22.5),
           nocolorcuts (bool, optional): Do not apply the fiducial rzW1W2 color-cuts
             cuts (default False).
           verbose (bool, optional): Be verbose!
-        Returns:
-          outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame spectra (erg/s/cm2/A).
-          wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
-          meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
+
+        Returns (outflux, wave, meta) tuple where:
+
+          * outflux (numpy.ndarray): Array [nmodel, npix] of observed-frame
+            spectra (1e-17 erg/s/cm2/A).
+          * wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
+          * meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
 
         Raises:
           ValueError
@@ -1876,19 +1907,20 @@ def make_templates(self, nmodel=100, zrange=(0.5, 4.0), rmagrange=(20.0, 22.5),
                         padflux, padzwave, mask_invalid=True)[self.normfilter])
                 magnorm = 10**(-0.4*mag[ii]) / normmaggies
 
-                synthnano = np.zeros((N_perz, len(self.decamwise)))
-                for ff, key in enumerate(maggies.columns):
-                    synthnano[:, ff] = 1E9 * maggies[key] * magnorm
+                synthnano = dict()
+                for key in maggies.columns:
+                    synthnano[key] = 1E9 * maggies[key] * magnorm
 
                 if nocolorcuts or self.colorcuts_function is None:
                     colormask = np.repeat(1, N_perz)
                 else:
                     colormask = self.colorcuts_function(
-                        gflux=synthnano[:, 1],
-                        rflux=synthnano[:, 2],
-                        zflux=synthnano[:, 4],
-                        w1flux=synthnano[:, 6],
-                        w2flux=synthnano[:, 7], optical=True)
+                        gflux=synthnano['decam2014-g'],
+                        rflux=synthnano['decam2014-r'],
+                        zflux=synthnano['decam2014-z'],
+                        w1flux=synthnano['wise2010-W1'],
+                        w2flux=synthnano['wise2010-W2'],
+                        optical=True)
 
                 # If the color-cuts pass then populate the output flux vector
                 # (suitably normalized) and metadata table and finish up.
@@ -1897,8 +1929,11 @@ def make_templates(self, nmodel=100, zrange=(0.5, 4.0), rmagrange=(20.0, 22.5),
                     outflux[ii, :] = resample_flux(self.wave, zwave[:, ii], flux[this, :],
                                                    extrapolate=True) * magnorm[this]
 
-                    meta['DECAM_FLUX'][ii] = synthnano[this, :6]
-                    meta['WISE_FLUX'][ii] = synthnano[this, 6:8]
+                    meta['FLUX_G'][ii] = synthnano['decam2014-g'][this]
+                    meta['FLUX_R'][ii] = synthnano['decam2014-r'][this]
+                    meta['FLUX_Z'][ii] = synthnano['decam2014-z'][this]
+                    meta['FLUX_W1'][ii] = synthnano['wise2010-W1'][this]
+                    meta['FLUX_W2'][ii] = synthnano['wise2010-W2'][this]
 
                     makemore = False
 
@@ -1913,4 +1948,4 @@ def make_templates(self, nmodel=100, zrange=(0.5, 4.0), rmagrange=(20.0, 22.5),
             log.warning('{} spectra could not be computed given the input priors!'.\
                         format(np.sum(success == 0)))
 
-        return outflux, self.wave, meta
+        return 1e17 * outflux, self.wave, meta
diff --git a/py/desisim/test/test_templates.py b/py/desisim/test/test_templates.py
index 96eb280d2..052344925 100644
--- a/py/desisim/test/test_templates.py
+++ b/py/desisim/test/test_templates.py
@@ -78,7 +78,7 @@ def test_OII(self):
         for i in range(len(meta)):
             z = meta['REDSHIFT'][i]
             ii = (3722*(1+z) < wave) & (wave < 3736*(1+z))
-            OIIflux = np.sum(flux[i,ii]*np.gradient(wave[ii]))
+            OIIflux = 1e-17 * np.sum(flux[i,ii] * np.gradient(wave[ii]))
             self.assertAlmostEqual(OIIflux, meta['OIIFLUX'][i], 2)
     
     @unittest.skipUnless(desi_basis_templates_available, '$DESI_BASIS_TEMPLATES was not detected.')
@@ -96,11 +96,11 @@ def test_HBETA(self):
             nmodel=10, zrange=(0.05, 0.4),
             logvdisp_meansig=[np.log10(75),0.0], 
             nocolorcuts=True, nocontinuum=True)
-    
+
         for i in range(len(meta)):
             z = meta['REDSHIFT'][i]
             ii = (4854*(1+z) < wave) & (wave < 4868*(1+z))
-            hbetaflux = np.sum(flux[i,ii]*np.gradient(wave[ii]))
+            hbetaflux = 1e-17 * np.sum(flux[i,ii] * np.gradient(wave[ii]))
             self.assertAlmostEqual(hbetaflux, meta['HBETAFLUX'][i], 2)
     
     @unittest.skipUnless(desi_basis_templates_available, '$DESI_BASIS_TEMPLATES was not detected.')
@@ -147,7 +147,7 @@ def test_input_meta(self):
             flux2, wave2, meta2 = Tx.make_templates(input_meta=meta1)
             badkeys = list()
             for key in meta1.colnames:
-                if key in ('DECAM_FLUX', 'WISE_FLUX', 'OIIFLUX', 'HBETAFLUX'):
+                if key in ('FLUX_G', 'FLUX_R', 'FLUX_Z', 'FLUX_W1', 'FLUX_W2', 'OIIFLUX', 'HBETAFLUX'):
                     #- not sure why the tolerances aren't closer
                     if not np.allclose(meta1[key], meta2[key], atol=5e-5):
                         print(meta1['OBJTYPE'][0], key, meta1[key], meta2[key])