diff --git a/autotest/t007_test.py b/autotest/t007_test.py
index a02324f7e..5ef574b3e 100644
--- a/autotest/t007_test.py
+++ b/autotest/t007_test.py
@@ -887,6 +887,17 @@ def test_read_usgs_model_reference():
     model_ws = os.path.join('temp', 't007')
     mrf = os.path.join(model_ws, 'usgs.model.reference')
     shutil.copy('../examples/data/usgs.model.reference', mrf)
+
+    xul, yul = 0, 0
+    with open(mrf) as foo:
+        for line in foo:
+            if 'xul' in line.lower():
+                xul = float(line.strip().split()[1])
+            elif "yul" in line.lower():
+                yul = float(line.strip().split()[1])
+            else:
+                continue
+
     fm = flopy.modflow
     m = fm.Modflow(modelname='junk', model_ws=model_ws)
     # feet and days
@@ -901,6 +912,11 @@ def test_read_usgs_model_reference():
     mg.read_usgs_model_reference_file(mrf)
     m2.modelgrid = mg
 
+    if abs(mg.xvertices[0, 0] - xul) > 0.01:
+        raise AssertionError()
+    if abs(mg.yvertices[0, 0] - yul) > 0.01:
+        raise AssertionError
+
     assert m2.modelgrid.xoffset == mg.xoffset
     assert m2.modelgrid.yoffset == mg.yoffset
     assert m2.modelgrid.angrot == mg.angrot
diff --git a/examples/data/mf6/test003_gwfs_disv/test003_gwfs_disv.dbf b/examples/data/mf6/test003_gwfs_disv/test003_gwfs_disv.dbf
index 09f199449..638500eea 100644
Binary files a/examples/data/mf6/test003_gwfs_disv/test003_gwfs_disv.dbf and b/examples/data/mf6/test003_gwfs_disv/test003_gwfs_disv.dbf differ
diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py
index 4f0ca183f..42d9ae7c6 100644
--- a/flopy/discretization/grid.py
+++ b/flopy/discretization/grid.py
@@ -420,6 +420,7 @@ def attribs_from_namfile_header(self, namefile):
         # check for reference info in the nam file header
         if namefile is None:
             return False
+        xul, yul = None, None
         header = []
         with open(namefile, 'r') as f:
             for line in f:
@@ -443,7 +444,6 @@ def attribs_from_namfile_header(self, namefile):
             elif "xul" in item.lower():
                 try:
                     xul = float(item.split(':')[1])
-                    self._xoff = self._xul_to_xll(xul)
                     warnings.warn(
                         'xul/yul have been deprecated. Use xll/yll instead.',
                         DeprecationWarning)
@@ -452,7 +452,6 @@ def attribs_from_namfile_header(self, namefile):
             elif "yul" in item.lower():
                 try:
                     yul = float(item.split(':')[1])
-                    self._yoff = self._yul_to_yll(yul)
                     warnings.warn(
                         'xul/yul have been deprecated. Use xll/yll instead.',
                         DeprecationWarning)
@@ -475,11 +474,21 @@ def attribs_from_namfile_header(self, namefile):
                     start_datetime = item.split(':')[1].strip()
                 except:
                     pass
+
+        # we need to rotate the modelgrid first, then we can
+        # calculate the xll and yll from xul and yul
+        if (xul, yul) != (None, None):
+            self.set_coord_info(xoff=self._xul_to_xll(xul),
+                                yoff=self._yul_to_yll(yul),
+                                angrot=self._angrot)
+
         return True
 
     def read_usgs_model_reference_file(self, reffile='usgs.model.reference'):
         """read spatial reference info from the usgs.model.reference file
         https://water.usgs.gov/ogw/policy/gw-model/modelers-setup.html"""
+        xul = None
+        yul = None
         if os.path.exists(reffile):
             with open(reffile) as input:
                 for line in input:
@@ -493,17 +502,9 @@ def read_usgs_model_reference_file(self, reffile='usgs.model.reference'):
                                 elif info[0] == 'yll':
                                     self._yoff = float(data)
                                 elif info[0] == 'xul':
-                                    self._xoff = self._xul_to_xll(
-                                        float(data))
-                                    warnings.warn(
-                                        'xul/yul have been deprecated. Use xll/yll instead.',
-                                        DeprecationWarning)
+                                    xul = float(data)
                                 elif info[0] == 'yul':
-                                    self._yoff = self._yul_to_yll(
-                                        float(data))
-                                    warnings.warn(
-                                        'xul/yul have been deprecated. Use xll/yll instead.',
-                                        DeprecationWarning)
+                                    yul = float(data)
                                 elif info[0] == 'rotation':
                                     self._angrot = float(data)
                                 elif info[0] == 'epsg':
@@ -512,6 +513,14 @@ def read_usgs_model_reference_file(self, reffile='usgs.model.reference'):
                                     self._proj4 = data
                                 elif info[0] == 'start_date':
                                     start_datetime = data
+
+            # model must be rotated first, before setting xoff and yoff
+            # when xul and yul are provided.
+            if (xul, yul) != (None, None):
+                self.set_coord_info(xoff=self._xul_to_xll(xul),
+                                    yoff=self._yul_to_yll(yul),
+                                    angrot=self._angrot)
+
             return True
         else:
             return False