diff --git a/grizli/combine.py b/grizli/combine.py
index df651426..e31e5397 100644
--- a/grizli/combine.py
+++ b/grizli/combine.py
@@ -173,7 +173,7 @@ def combine_flt(files=[], output='exposures_cmb.fits', grow=1,
             print('{0:3d} {1:s} {2:6.1f} {3:6.1f} {4:10.2f}'.format(i+1, file,
                                 x0[i], y0[i], im[0].header['EXPTIME']))
 
-        dq = utils.unset_dq_bits(im['DQ'].data, okbits=608,
+        dq = utils.mod_dq_bits(im['DQ'].data, okbits=608,
                                            verbose=False)
         wht = 1./im['ERR'].data**2
         wht[(im['ERR'].data == 0) | (dq > 0) | (~np.isfinite(wht))] = 0
diff --git a/grizli/jwst_utils.py b/grizli/jwst_utils.py
index 799ee5c1..fc22ac09 100644
--- a/grizli/jwst_utils.py
+++ b/grizli/jwst_utils.py
@@ -1036,7 +1036,7 @@ def exposure_oneoverf_correction(file, axis=None, thresholds=[5,4,3], erode_mask
         else:
             erode_mask = True
             
-    dq = utils.unset_dq_bits(im['DQ'].data, 4)
+    dq = utils.mod_dq_bits(im['DQ'].data, okbits=4)
     dqmask = dq == 0
     mask = dqmask
 
diff --git a/grizli/model.py b/grizli/model.py
index 8471f06f..8f10dcab 100644
--- a/grizli/model.py
+++ b/grizli/model.py
@@ -1687,7 +1687,7 @@ def __init__(self, sci=None, err=None, dq=None,
             self.MW_EBV = 0.
 
     def unset_dq(self):
-        """Flip OK data quality bits using utils.unset_dq_bits
+        """Flip OK data quality bits using utils.mod_dq_bits
 
         OK bits are defined as
 
@@ -1709,7 +1709,7 @@ def unset_dq(self):
         else:
             okbits = okbits_instrument[self.instrument]
 
-        self.data['DQ'] = utils.unset_dq_bits(self.data['DQ'], okbits=okbits)
+        self.data['DQ'] = utils.mod_dq_bits(self.data['DQ'], okbits=okbits)
 
 
     def flag_negative(self, sigma=-3):
diff --git a/grizli/prep.py b/grizli/prep.py
index 8210feab..49c6c980 100644
--- a/grizli/prep.py
+++ b/grizli/prep.py
@@ -6499,13 +6499,13 @@ def visit_grism_sky(grism={}, apply=True, column_average=True, verbose=True, ext
     if isACS:
         bits = 64+32
     elif isJWST:
-        bits = 1+4+6+32768 #+16777200+1049600+26232800+9438210+9438220
+        bits = 2+4+32768 #+16777200+1049600+26232800+9438210+9438220
     else:
         bits = 576
 
     for i in range(Nexp):
         flt = pyfits.open(grism['files'][i])
-        dq = utils.unset_dq_bits(flt['DQ', ext].data, okbits=bits)
+        dq = utils.mod_dq_bits(flt['DQ', ext].data, okbits=bits)
         dq_mask = dq == 0
 
         # Data
@@ -6895,7 +6895,7 @@ def fix_star_centers(root='macs1149.6+2223-rot-ca5-22-032.0-f105w',
 
                 sci = images[i]['SCI'].data[sly, slx]
                 dq = images[i]['DQ'].data[sly, slx]
-                dqm = dq - (dq & 2048)
+                dqm = utils.mod_dq_bits(dq,okbits=2048)
                 err = images[i]['ERR'].data[sly, slx]
                 mask = satpix[sly, slx]
 
diff --git a/grizli/utils.py b/grizli/utils.py
index f9609977..5f5fa353 100644
--- a/grizli/utils.py
+++ b/grizli/utils.py
@@ -1621,46 +1621,9 @@ def flt_to_dict(fobj, primary_keys=DEFAULT_PRIMARY_KEYS, extensions=[('SCI', i+1
 
     return flt_dict
 
-
-def get_set_bits(value):
-    """
-    Compute which binary bits are set for an integer
-    """
-
-    if hasattr(value, '__iter__'):
-        values = value
-        single = False
-    else:
-        values = [value]
-        single = True
-
-    result = []
-
-    for v in values:
-        try:
-            bitstr = np.binary_repr(v)[::-1]
-        except:
-            result.append([])
-
-        nset = bitstr.count('1')
-        setbits = []
-
-        j = -1
-        for i in range(nset):
-            j = bitstr.index('1', j+1)
-            setbits.append(j)
-
-        result.append(setbits)
-
-    if single:
-        return result[0]
-    else:
-        return result
-
-
-def unset_dq_bits(value, okbits=32+64+512, verbose=False):
+def mod_dq_bits(value, okbits=32+64+512, badbits=0, verbose=False):
     """
-    Unset bit flags from a DQ array
+    Modify bit flags from a DQ array
 
     For WFC3/IR, the following DQ bits can usually be unset:
 
@@ -1675,6 +1638,9 @@ def unset_dq_bits(value, okbits=32+64+512, verbose=False):
     okbits : int
         Bits to unset
 
+    badbits : int
+        Bits to set
+
     verbose : bool
         Print some information
 
@@ -1683,16 +1649,12 @@ def unset_dq_bits(value, okbits=32+64+512, verbose=False):
     new_value : int, `~numpy.ndarray`
 
     """
-    bin_bits = np.binary_repr(okbits)
-    n = len(bin_bits)
-    for i in range(n):
-        if bin_bits[-(i+1)] == '1':
-            if verbose:
-                print(2**i)
 
-            value -= (value & 2**i)
+    if verbose:
+         print(f'Unset bits: {np.binary_repr(okbits)}')
+         print(f'Set bits: {np.binary_repr(badbits)}')
 
-    return value
+    return (value & ~okbits) | badbits
 
 
 def detect_with_photutils(sci, err=None, dq=None, seg=None, detect_thresh=2.,
@@ -6312,11 +6274,11 @@ def drizzle_from_visit(visit, output=None, pixfrac=1., kernel='point',
                     dq = flt[('DQ', ext)].data & (1+1024+4096)
                     dq |= bpdata.astype(dq.dtype)
                     
-                    # dq0 = unset_dq_bits(flt[('DQ', ext)].data, bits) | bpdata
+                    # dq0 = mod_dq_bits(flt[('DQ', ext)].data, okbits=bits) | bpdata
                     # print('xxx', (dq > 0).sum(), (dq0 > 0).sum())
                     
                 else:
-                    dq = unset_dq_bits(flt[('DQ', ext)].data, bits) | bpdata
+                    dq = mod_dq_bits(flt[('DQ', ext)].data, okbits=bits) | bpdata
                     
                     
                 wht = 1/err**2