diff --git a/numpy/f2py/auxfuncs.py b/numpy/f2py/auxfuncs.py
index 3f9b0ceafa21..75c20513fe80 100644
--- a/numpy/f2py/auxfuncs.py
+++ b/numpy/f2py/auxfuncs.py
@@ -16,6 +16,7 @@
 """
 import pprint
 import sys
+import re
 import types
 from functools import reduce
 
@@ -26,7 +27,7 @@
     'applyrules', 'debugcapi', 'dictappend', 'errmess', 'gentitle',
     'getargs2', 'getcallprotoargument', 'getcallstatement',
     'getfortranname', 'getpymethoddef', 'getrestdoc', 'getusercode',
-    'getusercode1', 'hasbody', 'hascallstatement', 'hascommon',
+    'getusercode1', 'getdimension', 'hasbody', 'hascallstatement', 'hascommon',
     'hasexternals', 'hasinitvalue', 'hasnote', 'hasresultnote',
     'isallocatable', 'isarray', 'isarrayofstrings',
     'ischaracter', 'ischaracterarray', 'ischaracter_or_characterarray',
@@ -417,6 +418,13 @@ def isexternal(var):
     return 'attrspec' in var and 'external' in var['attrspec']
 
 
+def getdimension(var):
+    dimpattern = r"\((.*?)\)"
+    if 'attrspec' in var.keys():
+        if any('dimension' in s for s in var['attrspec']):
+            return [re.findall(dimpattern, v) for v in var['attrspec']][0]
+
+
 def isrequired(var):
     return not isoptional(var) and isintent_nothide(var)
 
diff --git a/numpy/f2py/crackfortran.py b/numpy/f2py/crackfortran.py
index d2dbb56afbd2..92c112ad7e15 100755
--- a/numpy/f2py/crackfortran.py
+++ b/numpy/f2py/crackfortran.py
@@ -1417,10 +1417,10 @@ def analyzeline(m, case, line):
                 outmess(
                     'analyzeline: implied-DO list "%s" is not supported. Skipping.\n' % l[0])
                 continue
-            i = 0
-            j = 0
             llen = len(l[1])
-            for v in rmbadname([x.strip() for x in markoutercomma(l[0]).split('@,@')]):
+            for idx, v in enumerate(rmbadname(
+                    [x.strip() for x in markoutercomma(l[0]).split('@,@')])
+                                    ):
                 if v[0] == '(':
                     outmess(
                         'analyzeline: implied-DO list "%s" is not supported. Skipping.\n' % v)
@@ -1429,18 +1429,26 @@ def analyzeline(m, case, line):
                     # wrapping.
                     continue
                 fc = 0
-                while (i < llen) and (fc or not l[1][i] == ','):
-                    if l[1][i] == "'":
-                        fc = not fc
-                    i = i + 1
-                i = i + 1
+                vtype = vars[v].get('typespec')
+                vdim = getdimension(vars[v])
+
+                if (vtype == 'complex'):
+                    cmplxpat = r"\(.*?\)"
+                    matches = re.findall(cmplxpat, l[1])
+                else:
+                    matches = l[1].split(',')
+
                 if v not in vars:
                     vars[v] = {}
-                if '=' in vars[v] and not vars[v]['='] == l[1][j:i - 1]:
+                if '=' in vars[v] and not vars[v]['='] == matches[idx]:
                     outmess('analyzeline: changing init expression of "%s" ("%s") to "%s"\n' % (
-                        v, vars[v]['='], l[1][j:i - 1]))
-                vars[v]['='] = l[1][j:i - 1]
-                j = i
+                        v, vars[v]['='], matches[idx]))
+
+                if vdim is not None:
+                    # Need to assign multiple values to one variable
+                    vars[v]['='] = "(/{}/)".format(", ".join(matches))
+                else:
+                    vars[v]['='] = matches[idx]
                 last_name = v
         groupcache[groupcounter]['vars'] = vars
         if last_name is not None:
diff --git a/numpy/f2py/tests/src/crackfortran/data_stmts.f90 b/numpy/f2py/tests/src/crackfortran/data_stmts.f90
new file mode 100644
index 000000000000..0eb97a563899
--- /dev/null
+++ b/numpy/f2py/tests/src/crackfortran/data_stmts.f90
@@ -0,0 +1,18 @@
+! gh-23276
+module cmplxdat
+  implicit none
+  integer :: i, j
+  real :: x, y
+  real, dimension(2) :: z
+  complex(kind=8), target :: medium_ref_index
+  complex(kind=8), target :: ref_index_one, ref_index_two
+  complex(kind=8), dimension(2) :: my_array
+  real(kind=8), dimension(3) :: my_real_array = (/1.0d0, 2.0d0, 3.0d0/)
+
+  data i, j / 2, 3 /
+  data x, y / 1.5, 2.0 /
+  data z / 3.5, 7.0 /
+  data medium_ref_index / (1.d0, 0.d0) /
+  data ref_index_one, ref_index_two / (13.0d0, 21.0d0), (-30.0d0, 43.0d0) /
+  data my_array / (1.0d0, 2.0d0), (-3.0d0, 4.0d0) /
+end module cmplxdat
diff --git a/numpy/f2py/tests/test_data.py b/numpy/f2py/tests/test_data.py
new file mode 100644
index 000000000000..3b8ca544ec74
--- /dev/null
+++ b/numpy/f2py/tests/test_data.py
@@ -0,0 +1,33 @@
+import os
+import pytest
+import numpy as np
+
+from . import util
+from numpy.f2py.crackfortran import crackfortran
+
+
+class TestData(util.F2PyTest):
+    sources = [util.getpath("tests", "src", "crackfortran", "data_stmts.f90")]
+
+    # For gh-23276
+    def test_data_stmts(self):
+        assert self.module.cmplxdat.i == 2
+        assert self.module.cmplxdat.j == 3
+        assert self.module.cmplxdat.x == 1.5
+        assert self.module.cmplxdat.y == 2.0
+        assert self.module.cmplxdat.medium_ref_index == np.array(1.+0.j)
+        assert np.all(self.module.cmplxdat.z == np.array([3.5, 7.0]))
+        assert np.all(self.module.cmplxdat.my_array == np.array([ 1.+2.j, -3.+4.j]))
+        assert np.all(self.module.cmplxdat.my_real_array == np.array([ 1., 2., 3.]))
+        assert np.all(self.module.cmplxdat.ref_index_one == np.array([13.0 + 21.0j]))
+        assert np.all(self.module.cmplxdat.ref_index_two == np.array([-30.0 + 43.0j]))
+
+    def test_crackedlines(self):
+        mod = crackfortran(self.sources)
+        assert mod[0]['vars']['x']['='] == '1.5'
+        assert mod[0]['vars']['y']['='] == '2.0'
+        assert mod[0]['vars']['my_real_array']['='] == '(/1.0d0, 2.0d0, 3.0d0/)'
+        assert mod[0]['vars']['ref_index_one']['='] == '(13.0d0, 21.0d0)'
+        assert mod[0]['vars']['ref_index_two']['='] == '(-30.0d0, 43.0d0)'
+        assert mod[0]['vars']['my_array']['='] == '(/(1.0d0, 2.0d0), (-3.0d0, 4.0d0)/)'
+        assert mod[0]['vars']['z']['='] == '(/3.5,  7.0/)'