diff --git a/AUTHORS.md b/AUTHORS.md
index d82ccb599fe9..c1a294c8c38b 100644
--- a/AUTHORS.md
+++ b/AUTHORS.md
@@ -120,6 +120,7 @@ Mor
 Nat-1
 Nicola Fonzi
 Nijso Beishuizen
+Nikita Ageev
 Ole Burghardt
 Patrick Mischke
 Paul Urbanczyk
diff --git a/SU2_PY/SU2/eval/functions.py b/SU2_PY/SU2/eval/functions.py
index fd6c1103a519..15803c688957 100644
--- a/SU2_PY/SU2/eval/functions.py
+++ b/SU2_PY/SU2/eval/functions.py
@@ -105,7 +105,7 @@ def function(func_name, config, state=None):
             multipoint(config, state)
 
         # Geometry
-        elif func_name in su2io.optnames_geo:
+        elif su2io.is_optname_geo(func_name):
             geometry(func_name, config, state)
 
         else:
@@ -890,8 +890,8 @@ def geometry(func_name, config, state=None):
 
     # return output
     funcs = su2util.ordered_bunch()
-    for key in su2io.optnames_geo:
-        if key in state["FUNCTIONS"]:
+    for key in state["FUNCTIONS"]:
+        if su2io.is_optname_geo(key):
             funcs[key] = state["FUNCTIONS"][key]
     return funcs
 
diff --git a/SU2_PY/SU2/eval/gradients.py b/SU2_PY/SU2/eval/gradients.py
index 5af57f6c33b0..bd7b778ab393 100644
--- a/SU2_PY/SU2/eval/gradients.py
+++ b/SU2_PY/SU2/eval/gradients.py
@@ -109,7 +109,7 @@ def gradient(func_name, method, config, state=None):
                 grads = multipoint(func_name, config, state)
 
             # Geometry (actually a finite difference)
-            elif func_output in su2io.optnames_geo:
+            elif su2io.is_optname_geo(func_output):
                 grads = geometry(func_name, config, state)
 
             else:
@@ -1080,8 +1080,8 @@ def geometry(func_name, config, state=None):
 
     # return output
     grads = su2util.ordered_bunch()
-    for key in su2io.optnames_geo:
-        if key in state["GRADIENTS"]:
+    for key in state["GRADIENTS"]:
+        if su2io.is_optname_geo(key):
             grads[key] = state["GRADIENTS"][key]
     return grads
 
diff --git a/SU2_PY/SU2/io/tools.py b/SU2_PY/SU2/io/tools.py
index aaeccf9433fc..504cfb513391 100755
--- a/SU2_PY/SU2/io/tools.py
+++ b/SU2_PY/SU2/io/tools.py
@@ -30,6 +30,7 @@
 # -------------------------------------------------------------------
 
 import os
+import re
 import shutil, glob
 from SU2.util import ordered_bunch
 from .historyMap import history_header_map as historyOutFields
@@ -271,20 +272,30 @@ def getTurboPerfIndex(nZones=1):
     "NACELLE_MAX_TWIST",
 ]
 
-PerStation = []
-for i in range(20):
-    PerStation.append("STATION" + str(i) + "_AREA")
-    PerStation.append("STATION" + str(i) + "_LENGTH")
-    PerStation.append("STATION" + str(i) + "_WIDTH")
-    PerStation.append("STATION" + str(i) + "_WATERLINE_WIDTH")
-    PerStation.append("STATION" + str(i) + "_HEIGHT")
-    PerStation.append("STATION" + str(i) + "_THICKNESS")
-    PerStation.append("STATION" + str(i) + "_CHORD")
-    PerStation.append("STATION" + str(i) + "_LE_RADIUS")
-    PerStation.append("STATION" + str(i) + "_TOC")
-    PerStation.append("STATION" + str(i) + "_TWIST")
-
-optnames_geo.extend(PerStation)
+# SU2_GEO numbers the stations of GEO_LOCATION_STATIONS from 1 (STATION1_*),
+# see SU2_GEO.cpp, and does not limit their number.
+optnames_geo_station = [
+    "AREA",
+    "LENGTH",
+    "WIDTH",
+    "WATERLINE_WIDTH",
+    "HEIGHT",
+    "THICKNESS",
+    "CHORD",
+    "LE_RADIUS",
+    "TOC",
+    "TWIST",
+]
+
+optnames_geo_station_regex = re.compile(
+    r"^STATION[1-9][0-9]*_(%s)$" % "|".join(optnames_geo_station)
+)
+
+
+def is_optname_geo(name):
+    """Check if name is a geometry function, including STATION#_* for any station."""
+    return name in optnames_geo or bool(optnames_geo_station_regex.match(name))
+
 
 #: optnames_geo
 
diff --git a/config_template.cfg b/config_template.cfg
index 5e611f22ed48..72f2d042017b 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -2212,7 +2212,7 @@ GEO_DESCRIPTION= AIRFOIL
 % Z location of the waterline
 GEO_WATERLINE_LOCATION= 0.0
 %
-% Coordinate of the stations to be analyzed
+% Coordinate of the stations to be analyzed (STATION1_* ... STATIONn_* in SU2_GEO output)
 GEO_LOCATION_STATIONS= (0.0, 0.5, 1.0)
 %
 % Definition of the nacelle location (higlite coordinates, tilt angle, toe angle)
@@ -2735,7 +2735,7 @@ REORIENT_ELEMENTS= YES
 %    AIRFOIL_AREA, AIRFOIL_THICKNESS, AIRFOIL_CHORD, AIRFOIL_TOC, AIRFOIL_AOA,
 %    WING_VOLUME, WING_MIN_THICKNESS, WING_MAX_THICKNESS, WING_MAX_CHORD, WING_MIN_TOC, WING_MAX_TWIST, WING_MAX_CURVATURE, WING_MAX_DIHEDRAL
 %    STATION#_WIDTH, STATION#_AREA, STATION#_THICKNESS, STATION#_CHORD, STATION#_TOC,
-%    STATION#_TWIST (where # is the index of the station defined in GEO_LOCATION_STATIONS)
+%    STATION#_TWIST (where # is the 1-based index of the station defined in GEO_LOCATION_STATIONS)
 %
 % Available design variables
 % 2D Design variables