diff --git a/.github/windows_arm64_steps/action.yml b/.github/windows_arm64_steps/action.yml
new file mode 100644
index 0000000..65f6530
--- /dev/null
+++ b/.github/windows_arm64_steps/action.yml
@@ -0,0 +1,28 @@
+name: Build Dependencies(Win-ARM64)
+description: "Common setup steps for Win-ARM64 CI"
+runs:
+  using: "composite"
+  steps:
+    - name: Install LLVM
+      shell: pwsh
+      run: |
+        Invoke-WebRequest https://github.com/llvm/llvm-project/releases/download/llvmorg-20.1.8/LLVM-20.1.8-woa64.exe -UseBasicParsing -OutFile LLVM-woa64.exe
+        $expectedHash = "7c4ac97eb2ae6b960ca5f9caf3ff6124c8d2a18cc07a7840a4d2ea15537bad8e"
+        $fileHash = (Get-FileHash -Path "LLVM-woa64.exe" -Algorithm SHA256).Hash
+        if ($fileHash -ne $expectedHash) {
+            Write-Error "Checksum verification failed. The downloaded file may be corrupted or tampered with."
+            exit 1
+        }
+        Start-Process -FilePath ".\LLVM-woa64.exe" -ArgumentList "/S" -Wait
+        echo "C:\Program Files\LLVM\bin" | Out-File -FilePath $env:GITHUB_PATH -Encoding utf8 -Append
+
+    - name: Install pkgconf via vcpkg
+      shell: pwsh
+      run: |
+        & "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvarsall.bat" arm64
+        $env:VCPKG_ROOT = "C:\vcpkg"
+        Set-Location $env:VCPKG_ROOT
+        ./vcpkg install pkgconf:arm64-windows
+        $pkgconfPath = "$env:VCPKG_ROOT\installed\arm64-windows\tools\pkgconf"
+        Copy-Item "$pkgconfPath\pkgconf.exe" "$pkgconfPath\pkg-config.exe" -Force
+        echo "$pkgconfPath" | Out-File -FilePath $env:GITHUB_PATH -Encoding utf8 -Append
\ No newline at end of file
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 4153b60..6510c86 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -1,46 +1,46 @@
-name: MkDocs
+# name: MkDocs
 
-on:
-  push:
-    branches:
-      - main
-  pull_request:
-    branches:
-      - main
-      - dev
+# on:
+#   push:
+#     branches:
+#       - main
+#   pull_request:
+#     branches:
+#       - main
+#       - dev
 
-jobs:
-  deploy:
-    runs-on: ubuntu-latest
-    steps:
-    - name: Checkout Repository
-      uses: actions/checkout@v2
+# jobs:
+#   deploy:
+#     runs-on: ubuntu-latest
+#     steps:
+#     - name: Checkout Repository
+#       uses: actions/checkout@v2
 
-    - name: Set up Python
-      uses: actions/setup-python@v2
-      with:
-        python-version: '3.x'
+#     - name: Set up Python
+#       uses: actions/setup-python@v2
+#       with:
+#         python-version: '3.x'
 
-    - name: Install Dependencies
-      run: |
-        python -m pip install --upgrade pip
-        pip install mkdocs-material
-        pip install markdown
-        pip install pymdown-extensions
-        pip install Pygments 
-        pip install mkdocstrings
-        pip install 'mkdocstrings[python]'
-        pip install mkdocs-gen-files
-        pip install mkdocs-literate-nav
-        pip install mkdocs-include-markdown-plugin
+#     - name: Install Dependencies
+#       run: |
+#         python -m pip install --upgrade pip
+#         pip install mkdocs-material
+#         pip install markdown
+#         pip install pymdown-extensions
+#         pip install Pygments 
+#         pip install mkdocstrings
+#         pip install 'mkdocstrings[python]'
+#         pip install mkdocs-gen-files
+#         pip install mkdocs-literate-nav
+#         pip install mkdocs-include-markdown-plugin
         
-    - name: Test Build
-      run: |
-        # mkdocs build --strict # can't use strict mode because of some typ annotations
-        mkdocs build
+#     - name: Test Build
+#       run: |
+#         # mkdocs build --strict # can't use strict mode because of some typ annotations
+#         mkdocs build
         
 
-    - name: Build and Deploy
-      if: github.event_name == 'push'
-      run: |
-        mkdocs gh-deploy --force
+#     - name: Build and Deploy
+#       if: github.event_name == 'push'
+#       run: |
+#         mkdocs gh-deploy --force
diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml
index 85850f9..e2a2fd1 100644
--- a/.github/workflows/wheels.yml
+++ b/.github/workflows/wheels.yml
@@ -26,14 +26,16 @@ jobs:
         # - [macos-14, macosx, arm64, accelerate, "14.0"]
         
         # working
-        - [windows-2025, win, AMD64]
-        - [ubuntu-22.04, manylinux, x86_64]
-        - [macos-14, macosx, arm64, openblas, "12.3"]
-        - [macos-15-intel, macosx, x86_64, openblas, "10.14"]
+        # - [windows-2025, win, AMD64]
+        # - [ubuntu-22.04, manylinux, x86_64]
+        # - [macos-14, macosx, arm64, openblas, "12.3"]
+        - [windows-11-arm, win, ARM64, "", ""]
 
+        # removed 
+        # - [macos-15-intel, macosx, x86_64, openblas, "10.14"]
         
-        python: [["cp39", "3.9"], ["cp310", "3.10"], ["cp311", "3.11"], ["cp312", "3.12"], ["cp313", "3.13"]]
-        # python: [["cp312", "3.12"]]
+        # python: [["cp39", "3.9"], ["cp310", "3.10"], ["cp311", "3.11"], ["cp312", "3.12"], ["cp313", "3.13"]]
+        python: [["cp312", "3.12"]]
 
 
 
@@ -55,6 +57,19 @@ jobs:
         if: ${{ runner.os == 'Windows' }}
 
 
+      - name: Set environment variables for ARM64
+        if: matrix.buildplat[1] == 'win' && matrix.buildplat[2] == 'ARM64'
+        run: |
+            echo "CC=clang-cl" >> $env:GITHUB_ENV
+            echo "CXX=clang-cl" >> $env:GITHUB_ENV
+            echo "FC=flang" >> $env:GITHUB_ENV
+            echo "TARGET_ARCH=${{ matrix.buildplat[2] }}" >> $env:GITHUB_ENV
+   
+      - name: Set up Flang and pkgconf for ARM64
+        if: matrix.buildplat[1] == 'win' && matrix.buildplat[2] == 'ARM64'
+        uses: ./.github/windows_arm64_steps
+
+
       - name: Setup macOS
         if: startsWith( matrix.buildplat[0], 'macos-' )
         run: |
diff --git a/.gitignore b/.gitignore
index 01894e0..aa0b628 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,6 +18,9 @@ dev_make_install/optvl
 
 tests/test_om_wrapper*_out
 examples/run_opt_om*_out
+examples/*.avl
+examples/*.ps
+examples/opt_output_sweep/
 # ==============================================================================
 # Github python gitignore template
 # ==============================================================================
@@ -52,7 +55,6 @@ MANIFEST
 src/ad_src/forward_tmp/
 src/ad_src/reverse_tmp/
 src/ad_src/preprocessed_files/
-src/includes/AVL.INC
 examples/*.html
 examples/*.db
 examples/reports
diff --git a/meson.build b/meson.build
index 89d6410..2082aeb 100644
--- a/meson.build
+++ b/meson.build
@@ -8,7 +8,6 @@ project(
   default_options: [
     'buildtype=debugoptimized',
     'c_std=c99',
-    'fortran_std=legacy',
   ],
 )
 
@@ -111,13 +110,28 @@ ff_args = ff.get_supported_arguments(
 
 ff_args += [
   '-ffixed-line-length-80', 
-  '-std=legacy',
   '-fdefault-real-8',
   '-fdefault-double-8',
-  '-fPIC',
   '-O2'
 ]
 
+if host_machine.system() != 'windows'
+  ff_args += ['-fPIC']
+endif
+
+# --- add flags for fortran standard ---
+if ff.get_id() == 'gcc'
+  # -std=legacy is not supported by all Fortran compilers, but very useful with
+  # gfortran since it avoids a ton of warnings that we don't care about.
+  # Needs fixing in Meson, see https://github.com/mesonbuild/meson/issues/11633.
+  ff_args += ['-std=legacy']
+elif ff.get_id() == 'llvm-flang'
+
+  ff_args += ['-fno-init-global-zero', '-ffixed-form', '-w']
+  link_args += ['-fno-init-global-zero']
+  
+endif
+
 numpy_nodepr_api = '-DNPY_NO_DEPRECATED_API=NPY_1_9_API_VERSION'
 
 cc_args = [numpy_nodepr_api]
@@ -179,4 +193,5 @@ py3.extension_module('libavl',
                    
 )
 
+
 install_subdir('optvl', install_dir: py3.get_install_dir())
diff --git a/optvl/optvl_class.py b/optvl/optvl_class.py
index 8ab2a53..ca8f6ae 100644
--- a/optvl/optvl_class.py
+++ b/optvl/optvl_class.py
@@ -286,6 +286,7 @@ def __init__(
 
         # # get just the file name
         avl_lib_so_file = os.path.basename(avl_lib_so_file)
+        
 
         self.avl = MExt.MExt("libavl", module_name, "optvl", lib_so_file=avl_lib_so_file, debug=debug)._module
 
diff --git a/pyproject.toml b/pyproject.toml
index 9876685..f31724f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -33,7 +33,7 @@ version = "2.4.0" # this automatically updates __init__.py and setup_deprecated.
 
 
 [tool.cibuildwheel]
-skip = "cp36-* cp37-* pp* *_ppc64le *_i686 *_s390x"
+skip = "cp38-* pp* *_ppc64le *_i686 *_s390x"
 build-verbosity = "3"
 test-command = "bash {project}/tools/wheels/cibw_test_command.sh {project}"
 
@@ -62,6 +62,11 @@ select = "*-win_amd64"
 # environment = { CMAKE_PREFIX_PATH="c:/opt/64" }
 environment = { PKG_CONFIG_PATH="c:/opt/64/lib/pkgconfig" }
 
+# pyproject.toml
+[tool.cibuildwheel.config-settings]
+setup-args = ["-Dbuildtype=debug"]   # optional, but helps
+compile-args = ["-v"]                # this is the key one — makes ninja verbose
+
 [project.urls]
 homepage = "https://github.com/joanibal/OptVL"
 # documentation =
diff --git a/src/ad_src/forward_ad_src/amake_d.f b/src/ad_src/forward_ad_src/amake_d.f
index 17bf2b2..1d67176 100644
--- a/src/ad_src/forward_ad_src/amake_d.f
+++ b/src/ad_src/forward_ad_src/amake_d.f
@@ -18,6 +18,7 @@ SUBROUTINE UPDATE_SURFACES_D()
       use avl_heap_inc
       use avl_heap_diff_inc
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       INCLUDE 'AVL_ad_seeds.inc'
       INTEGER ii
       INTEGER isurf
@@ -201,6 +202,8 @@ SUBROUTINE UPDATE_SURFACES_D()
 C
       SUBROUTINE MAKESURF_D(isurf)
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
+      INCLUDE 'AVL_oml.INC'
       INCLUDE 'AVL_ad_seeds.inc'
 C
 C
@@ -1254,41 +1257,41 @@ SUBROUTINE MAKESURF_D(isurf)
 C---------- TE control point used only if surface sheds a wake
                 lvnc(idx_vor) = lfwake(isurf)
 C
-Cc#ifdef USE_CPOML
-C...        nodal grid associated with vortex strip (aft-panel nodes)
-C...        NOTE: airfoil in plane of wing, but not rotated perpendicular to dihedral;
-C...        retained in (x,z) plane at this point
-                CALL AKIMA(xlasec(1, isec, isurf), zlasec(1, isec, isurf
-     +                     ), nsl, xpt(ivc+1), zl_l, dsdx)
-                CALL AKIMA(xuasec(1, isec, isurf), zuasec(1, isec, isurf
-     +                     ), nsl, xpt(ivc+1), zu_l, dsdx)
-C
-                CALL AKIMA(xlasec(1, isec+1, isurf), zlasec(1, isec+1, 
-     +                     isurf), nsr, xpt(ivc+1), zl_r, dsdx)
-                CALL AKIMA(xuasec(1, isec+1, isurf), zuasec(1, isec+1, 
-     +                     isurf), nsr, xpt(ivc+1), zu_r, dsdx)
-C
-                xyn1(1, idx_vor) = rle1(1, idx_strip) + xpt(ivc+1)*
-     +            chord1(idx_strip)
-                xyn1(2, idx_vor) = rle1(2, idx_strip)
-                zl = (1.-f1)*zl_l + f1*zl_r
-                zu = (1.-f1)*zu_l + f1*zu_r
-                zlon1(idx_vor) = rle1(3, idx_strip) + zl*chord1(
-     +            idx_strip)
-                zupn1(idx_vor) = rle1(3, idx_strip) + zu*chord1(
-     +            idx_strip)
-C
-                xyn2(1, idx_vor) = rle2(1, idx_strip) + xpt(ivc+1)*
-     +            chord2(idx_strip)
-                xyn2(2, idx_vor) = rle2(2, idx_strip)
-                zl = (1.-f2)*zl_l + f2*zl_r
-                zu = (1.-f2)*zu_l + f2*zu_r
-                zlon2(idx_vor) = rle2(3, idx_strip) + zl*chord2(
-     +            idx_strip)
-                zupn2(idx_vor) = rle2(3, idx_strip) + zu*chord2(
-     +            idx_strip)
-C
-Cc#endif
+! Cc#ifdef USE_CPOML
+! C...        nodal grid associated with vortex strip (aft-panel nodes)
+! C...        NOTE: airfoil in plane of wing, but not rotated perpendicular to dihedral;
+! C...        retained in (x,z) plane at this point
+!                 CALL AKIMA(xlasec(1, isec, isurf), zlasec(1, isec, isurf
+!      +                     ), nsl, xpt(ivc+1), zl_l, dsdx)
+!                 CALL AKIMA(xuasec(1, isec, isurf), zuasec(1, isec, isurf
+!      +                     ), nsl, xpt(ivc+1), zu_l, dsdx)
+! C
+!                 CALL AKIMA(xlasec(1, isec+1, isurf), zlasec(1, isec+1, 
+!      +                     isurf), nsr, xpt(ivc+1), zl_r, dsdx)
+!                 CALL AKIMA(xuasec(1, isec+1, isurf), zuasec(1, isec+1, 
+!      +                     isurf), nsr, xpt(ivc+1), zu_r, dsdx)
+! C
+!                 xyn1(1, idx_vor) = rle1(1, idx_strip) + xpt(ivc+1)*
+!      +            chord1(idx_strip)
+!                 xyn1(2, idx_vor) = rle1(2, idx_strip)
+!                 zl = (1.-f1)*zl_l + f1*zl_r
+!                 zu = (1.-f1)*zu_l + f1*zu_r
+!                 zlon1(idx_vor) = rle1(3, idx_strip) + zl*chord1(
+!      +            idx_strip)
+!                 zupn1(idx_vor) = rle1(3, idx_strip) + zu*chord1(
+!      +            idx_strip)
+! C
+!                 xyn2(1, idx_vor) = rle2(1, idx_strip) + xpt(ivc+1)*
+!      +            chord2(idx_strip)
+!                 xyn2(2, idx_vor) = rle2(2, idx_strip)
+!                 zl = (1.-f2)*zl_l + f2*zl_r
+!                 zu = (1.-f2)*zu_l + f2*zu_r
+!                 zlon2(idx_vor) = rle2(3, idx_strip) + zl*chord2(
+!      +            idx_strip)
+!                 zupn2(idx_vor) = rle2(3, idx_strip) + zu*chord2(
+!      +            idx_strip)
+! C
+! Cc#endif
                 idx_vor = idx_vor + 1
               ENDDO
 C           
@@ -1341,6 +1344,8 @@ SUBROUTINE MAKESURF_D(isurf)
 C
       SUBROUTINE MAKESURF_MESH_D(isurf)
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
+      INCLUDE 'AVL_oml.INC'
       INCLUDE 'AVL_ad_seeds.inc'
 C working variables (AVL original)
       INTEGER isurf
@@ -2540,49 +2545,49 @@ SUBROUTINE MAKESURF_MESH_D(isurf)
      +        idx_node_yp1+1))/2-rle(1, idx_strip))/chord(idx_strip)
 C
 C
-            CALL AKIMA(xlasec(1, iptl, isurf), zlasec(1, iptl, isurf), 
-     +                 nsl, xptxind1, zl_l, dsdx)
-            CALL AKIMA(xuasec(1, iptl, isurf), zuasec(1, iptl, isurf), 
-     +                 nsl, xptxind1, zu_l, dsdx)
-C Interpolate cross section on right side
-C
-            CALL AKIMA(xlasec(1, iptr, isurf), zlasec(1, iptr, isurf), 
-     +                 nsr, xptxind1, zl_r, dsdx)
-            CALL AKIMA(xuasec(1, iptr, isurf), zuasec(1, iptr, isurf), 
-     +                 nsr, xptxind1, zu_r, dsdx)
-C Compute the left aft node of panel 
-C X-point
-C
-C
-C Y-point
-            xyn1(1, idx_vor) = rle1(1, idx_strip) + xptxind1*chord1(
-     +        idx_strip)
-C
-C Interpolate z from sections to left aft node of panel
-            xyn1(2, idx_vor) = rle1(2, idx_strip)
-C
-            zl = (1.-f1)*zl_l + f1*zl_r
-C Store left aft z-point
-            zu = (1.-f1)*zu_l + f1*zu_r
-C
-            zlon1(idx_vor) = rle1(3, idx_strip) + zl*chord1(idx_strip)
-C Compute the right aft node of panel 
-C X-point
-            zupn1(idx_vor) = rle1(3, idx_strip) + zu*chord1(idx_strip)
-C
-C Y-point
-            xyn2(1, idx_vor) = rle2(1, idx_strip) + xptxind1*chord2(
-     +        idx_strip)
-C
-C Interpolate z from sections to right aft node of panel
-            xyn2(2, idx_vor) = rle2(2, idx_strip)
-            zl = (1.-f2)*zl_l + f2*zl_r
-C Store right aft z-point
-            zu = (1.-f2)*zu_l + f2*zu_r
-C
-            zlon2(idx_vor) = rle2(3, idx_strip) + zl*chord2(idx_strip)
-            zupn2(idx_vor) = rle2(3, idx_strip) + zu*chord2(idx_strip)
-C
+!             CALL AKIMA(xlasec(1, iptl, isurf), zlasec(1, iptl, isurf), 
+!      +                 nsl, xptxind1, zl_l, dsdx)
+!             CALL AKIMA(xuasec(1, iptl, isurf), zuasec(1, iptl, isurf), 
+!      +                 nsl, xptxind1, zu_l, dsdx)
+! C Interpolate cross section on right side
+! C
+!             CALL AKIMA(xlasec(1, iptr, isurf), zlasec(1, iptr, isurf), 
+!      +                 nsr, xptxind1, zl_r, dsdx)
+!             CALL AKIMA(xuasec(1, iptr, isurf), zuasec(1, iptr, isurf), 
+!      +                 nsr, xptxind1, zu_r, dsdx)
+! C Compute the left aft node of panel 
+! C X-point
+! C
+! C
+! C Y-point
+!             xyn1(1, idx_vor) = rle1(1, idx_strip) + xptxind1*chord1(
+!      +        idx_strip)
+! C
+! C Interpolate z from sections to left aft node of panel
+!             xyn1(2, idx_vor) = rle1(2, idx_strip)
+! C
+!             zl = (1.-f1)*zl_l + f1*zl_r
+! C Store left aft z-point
+!             zu = (1.-f1)*zu_l + f1*zu_r
+! C
+!             zlon1(idx_vor) = rle1(3, idx_strip) + zl*chord1(idx_strip)
+! C Compute the right aft node of panel 
+! C X-point
+!             zupn1(idx_vor) = rle1(3, idx_strip) + zu*chord1(idx_strip)
+! C
+! C Y-point
+!             xyn2(1, idx_vor) = rle2(1, idx_strip) + xptxind1*chord2(
+!      +        idx_strip)
+! C
+! C Interpolate z from sections to right aft node of panel
+!             xyn2(2, idx_vor) = rle2(2, idx_strip)
+!             zl = (1.-f2)*zl_l + f2*zl_r
+! C Store right aft z-point
+!             zu = (1.-f2)*zu_l + f2*zu_r
+! C
+!             zlon2(idx_vor) = rle2(3, idx_strip) + zl*chord2(idx_strip)
+!             zupn2(idx_vor) = rle2(3, idx_strip) + zu*chord2(idx_strip)
+! C
             idx_vor = idx_vor + 1
           ENDDO
 C End vortex loop
@@ -2625,6 +2630,8 @@ SUBROUTINE MAKESURF_MESH_D(isurf)
 C
       SUBROUTINE SDUPL_D(nn, ypt, msg)
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
+      INCLUDE 'AVL_oml.INC'
       INCLUDE 'AVL_ad_seeds.inc'
       CHARACTER*(*) msg
       INTEGER idx_vor
@@ -2923,6 +2930,7 @@ SUBROUTINE SDUPL_D(nn, ypt, msg)
 C and uses the real mesh to compute normals if it is
       SUBROUTINE ENCALC_D()
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       INCLUDE 'AVL_ad_seeds.inc'
 C
       REAL ep(3), eq(3), es(3), eb(3), ec(3), ecxb(3)
diff --git a/src/ad_src/forward_ad_src/asetup_d.f b/src/ad_src/forward_ad_src/asetup_d.f
index 94ef1ec..efaedbe 100644
--- a/src/ad_src/forward_ad_src/asetup_d.f
+++ b/src/ad_src/forward_ad_src/asetup_d.f
@@ -170,23 +170,23 @@ SUBROUTINE VELSUM_D()
 C
       DO i=1,nvor
         DO k=1,3
-          vc(k, i) = 0.0
+          ! vc(k, i) = 0.0
           vv_diff(k, i) = 0.D0
           vv(k, i) = 0.0
 C--- h.v. velocity at control points and vortex midpoints
           DO j=1,nvor
-            vc(k, i) = vc(k, i) + wc_gam(k, i, j)*gam(j)
+            ! vc(k, i) = vc(k, i) + wc_gam(k, i, j)*gam(j)
             vv_diff(k, i) = vv_diff(k, i) + gam(j)*wv_gam_diff(k, i, j) 
      +        + wv_gam(k, i, j)*gam_diff(j)
             vv(k, i) = vv(k, i) + wv_gam(k, i, j)*gam(j)
           ENDDO
           DO n=1,numax
-            vc_u(k, i, n) = 0.0
+            ! vc_u(k, i, n) = 0.0
             vv_u_diff(k, i, n) = 0.D0
             vv_u(k, i, n) = 0.0
             DO j=1,nvor
-              vc_u(k, i, n) = vc_u(k, i, n) + wc_gam(k, i, j)*gam_u(j, n
-     +          )
+              ! vc_u(k, i, n) = vc_u(k, i, n) + wc_gam(k, i, j)*gam_u(j, n
+    !  +          )
               vv_u_diff(k, i, n) = vv_u_diff(k, i, n) + gam_u(j, n)*
      +          wv_gam_diff(k, i, j) + wv_gam(k, i, j)*gam_u_diff(j, n)
               vv_u(k, i, n) = vv_u(k, i, n) + wv_gam(k, i, j)*gam_u(j, n
@@ -194,12 +194,12 @@ SUBROUTINE VELSUM_D()
             ENDDO
           ENDDO
           DO n=1,ncontrol
-            vc_d(k, i, n) = 0.0
+            ! vc_d(k, i, n) = 0.0
             vv_d_diff(k, i, n) = 0.D0
             vv_d(k, i, n) = 0.0
             DO j=1,nvor
-              vc_d(k, i, n) = vc_d(k, i, n) + wc_gam(k, i, j)*gam_d(j, n
-     +          )
+              ! vc_d(k, i, n) = vc_d(k, i, n) + wc_gam(k, i, j)*gam_d(j, n
+    !  +          )
               vv_d_diff(k, i, n) = vv_d_diff(k, i, n) + gam_d(j, n)*
      +          wv_gam_diff(k, i, j) + wv_gam(k, i, j)*gam_d_diff(j, n)
               vv_d(k, i, n) = vv_d(k, i, n) + wv_gam(k, i, j)*gam_d(j, n
@@ -207,11 +207,11 @@ SUBROUTINE VELSUM_D()
             ENDDO
           ENDDO
           DO n=1,ndesign
-            vc_g(k, i, n) = 0.0
+            ! vc_g(k, i, n) = 0.0
             vv_g(k, i, n) = 0.0
             DO j=1,nvor
-              vc_g(k, i, n) = vc_g(k, i, n) + wc_gam(k, i, j)*gam_g(j, n
-     +          )
+              ! vc_g(k, i, n) = vc_g(k, i, n) + wc_gam(k, i, j)*gam_g(j, n
+    !  +          )
               vv_g(k, i, n) = vv_g(k, i, n) + wv_gam(k, i, j)*gam_g(j, n
      +          )
             ENDDO
@@ -228,21 +228,21 @@ SUBROUTINE VELSUM_D()
      +      (2) + wvsrd_u(k, i, 3)*vinf(3) + wvsrd_u(k, i, 4)*wrot(1) + 
      +      wvsrd_u(k, i, 5)*wrot(2) + wvsrd_u(k, i, 6)*wrot(3)
 C--- total velocity at control points and vortex midpoints
-          wc(k, i) = vc(k, i) + wcsrd(k, i)
+          ! wc(k, i) = vc(k, i) + wcsrd(k, i)
           wv_diff(k, i) = vv_diff(k, i) + wvsrd_diff(k, i)
           wv(k, i) = vv(k, i) + wvsrd(k, i)
           DO n=1,numax
-            wc_u(k, i, n) = vc_u(k, i, n) + wcsrd_u(k, i, n)
+            ! wc_u(k, i, n) = vc_u(k, i, n) + wcsrd_u(k, i, n)
             wv_u_diff(k, i, n) = vv_u_diff(k, i, n)
             wv_u(k, i, n) = vv_u(k, i, n) + wvsrd_u(k, i, n)
           ENDDO
           DO n=1,ndmax
-            wc_d(k, i, n) = vc_d(k, i, n)
+            ! wc_d(k, i, n) = vc_d(k, i, n)
             wv_d_diff(k, i, n) = vv_d_diff(k, i, n)
             wv_d(k, i, n) = vv_d(k, i, n)
           ENDDO
           DO n=1,ngmax
-            wc_g(k, i, n) = vc_g(k, i, n)
+            ! wc_g(k, i, n) = vc_g(k, i, n)
             wv_g(k, i, n) = vv_g(k, i, n)
           ENDDO
         ENDDO
diff --git a/src/ad_src/reverse_ad_src/amake_b.f b/src/ad_src/reverse_ad_src/amake_b.f
index abb90cb..96c94a4 100644
--- a/src/ad_src/reverse_ad_src/amake_b.f
+++ b/src/ad_src/reverse_ad_src/amake_b.f
@@ -22,6 +22,7 @@ SUBROUTINE UPDATE_SURFACES_B()
       use avl_heap_inc
       use avl_heap_diff_inc
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       INCLUDE 'AVL_ad_seeds.inc'
       INTEGER ii
       INTEGER isurf
@@ -245,6 +246,7 @@ SUBROUTINE UPDATE_SURFACES_B()
 C
       SUBROUTINE MAKESURF_B(isurf)
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       INCLUDE 'AVL_ad_seeds.inc'
 C
 C
@@ -2237,6 +2239,7 @@ SUBROUTINE MAKESURF_B(isurf)
 C
       SUBROUTINE MAKESURF_MESH_B(isurf)
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       INCLUDE 'AVL_ad_seeds.inc'
 C working variables (AVL original)
       INTEGER isurf
@@ -3949,6 +3952,7 @@ SUBROUTINE MAKESURF_MESH_B(isurf)
 C
       SUBROUTINE SDUPL_B(nn, ypt, msg)
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       INCLUDE 'AVL_ad_seeds.inc'
       CHARACTER*(*) msg
       INTEGER idx_vor
@@ -4463,6 +4467,7 @@ SUBROUTINE SDUPL_B(nn, ypt, msg)
 C and uses the real mesh to compute normals if it is
       SUBROUTINE ENCALC_B()
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       INCLUDE 'AVL_ad_seeds.inc'
 C
       REAL ep(3), eq(3), es(3), eb(3), ec(3), ecxb(3)
diff --git a/src/ainput.f b/src/ainput.f
index 93c42cb..ce0fd98 100644
--- a/src/ainput.f
+++ b/src/ainput.f
@@ -23,6 +23,7 @@ SUBROUTINE INPUT(LUN,FNAME,FERR)
 C     Reads an processes an AVL configuration input file
 C---------------------------------------------------------
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
 C
       CHARACTER*(*) FNAME
       LOGICAL FERR
@@ -604,14 +605,14 @@ SUBROUTINE INPUT(LUN,FNAME,FERR)
         TASEC(2,ISEC, ISURF)  = 0.0
 cc#ifdef USE_CPOML
 C
-        XLASEC(1,ISEC,ISURF) = 0
-        XLASEC(2,ISEC,ISURF) = 0
-        XUASEC(1,ISEC,ISURF) = 0
-        XUASEC(2,ISEC,ISURF) = 0
-        ZLASEC(1,ISEC,ISURF) = 0
-        ZLASEC(2,ISEC,ISURF) = 0
-        ZUASEC(1,ISEC,ISURF) = 0
-        ZUASEC(2,ISEC,ISURF) = 0
+        ! XLASEC(1,ISEC,ISURF) = 0
+        ! XLASEC(2,ISEC,ISURF) = 0
+        ! XUASEC(1,ISEC,ISURF) = 0
+        ! XUASEC(2,ISEC,ISURF) = 0
+        ! ZLASEC(1,ISEC,ISURF) = 0
+        ! ZLASEC(2,ISEC,ISURF) = 0
+        ! ZUASEC(1,ISEC,ISURF) = 0
+        ! ZUASEC(2,ISEC,ISURF) = 0
         CASEC(2,ISEC,ISURF) = 0
         CASEC(2,ISEC,ISURF) = 0
 cc#endif
@@ -709,10 +710,10 @@ SUBROUTINE INPUT(LUN,FNAME,FERR)
             ZF = 0
           ENDIF
           TH = ATAN(SLP)
-          XLASEC(I,ISEC,ISURF) = XF + 0.5*THK*SIN(TH)
-          XUASEC(I,ISEC,ISURF) = XF - 0.5*THK*SIN(TH)
-          ZLASEC(I,ISEC,ISURF) = ZF - 0.5*THK*COS(TH)
-          ZUASEC(I,ISEC,ISURF) = ZF + 0.5*THK*COS(TH)
+        !   XLASEC(I,ISEC,ISURF) = XF + 0.5*THK*SIN(TH)
+        !   XUASEC(I,ISEC,ISURF) = XF - 0.5*THK*SIN(TH)
+        !   ZLASEC(I,ISEC,ISURF) = ZF - 0.5*THK*COS(TH)
+        !   ZUASEC(I,ISEC,ISURF) = ZF + 0.5*THK*COS(TH)
           CASEC(I,ISEC,ISURF) = ZF
 cc#endif
         ENDDO
@@ -724,15 +725,15 @@ SUBROUTINE INPUT(LUN,FNAME,FERR)
           idx_upper = NASEC(ISEC,ISURF) - I +1
           idx_lower = I 
           
-          XSEC(I                  ,isec, isurf) = 
-     &                                   XUASEC(idx_upper, ISEC, ISURF)
-          XSEC(I+NASEC(ISEC,ISURF),isec, isurf) = 
-     &                                   XLASEC(idx_lower, ISEC, ISURF)
+!           XSEC(I                  ,isec, isurf) = 
+!      &                                   XUASEC(idx_upper, ISEC, ISURF)
+!           XSEC(I+NASEC(ISEC,ISURF),isec, isurf) = 
+!      &                                   XLASEC(idx_lower, ISEC, ISURF)
 
-          YSEC(I                  ,isec, isurf) = 
-     &                                   ZUASEC(idx_upper, ISEC, ISURF)
-          YSEC(I+NASEC(ISEC,ISURF),isec, isurf) = 
-     &                                   ZLASEC(idx_lower, ISEC, ISURF)
+!           YSEC(I                  ,isec, isurf) = 
+!      &                                   ZUASEC(idx_upper, ISEC, ISURF)
+!           YSEC(I+NASEC(ISEC,ISURF),isec, isurf) = 
+!      &                                   ZLASEC(idx_lower, ISEC, ISURF)
 
         enddo 
 
@@ -826,10 +827,10 @@ SUBROUTINE INPUT(LUN,FNAME,FERR)
 cc#ifdef USE_CPOML
 C...      lower/upper coordinates (for oml)
           CALL AKIMA(XIN,YIN,NIN,XASEC(I,ISEC,ISURF),ZC,DUMMY)
-          XLASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF)
-          XUASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF)
-          ZLASEC(I,ISEC,ISURF) = ZC - 0.5*TASEC(I,ISEC,ISURF)
-          ZUASEC(I,ISEC,ISURF) = ZC + 0.5*TASEC(I,ISEC,ISURF)
+        !   XLASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF)
+        !   XUASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF)
+        !   ZLASEC(I,ISEC,ISURF) = ZC - 0.5*TASEC(I,ISEC,ISURF)
+        !   ZUASEC(I,ISEC,ISURF) = ZC + 0.5*TASEC(I,ISEC,ISURF)
           CASEC(I,ISEC,ISURF) = ZC
 cc#endif
         END DO
@@ -918,10 +919,10 @@ SUBROUTINE INPUT(LUN,FNAME,FERR)
           SASEC(I,ISEC,ISURF) = 0.0
           TASEC(I,ISEC,ISURF) = 0.0
 cc#ifdef USE_CPOML
-          XLASEC(I,ISEC,ISURF) = 0
-          XUASEC(I,ISEC,ISURF) = 0
-          ZLASEC(I,ISEC,ISURF) = 0
-          ZUASEC(I,ISEC,ISURF) = 0
+        !   XLASEC(I,ISEC,ISURF) = 0
+        !   XUASEC(I,ISEC,ISURF) = 0
+        !   ZLASEC(I,ISEC,ISURF) = 0
+        !   ZUASEC(I,ISEC,ISURF) = 0
           CASEC(I,ISEC,ISURF) = 0
 cc#endif
          ENDDO
@@ -949,10 +950,10 @@ SUBROUTINE INPUT(LUN,FNAME,FERR)
 cc#ifdef USE_CPOML
 C...       lower/upper coordinates (for oml)
            CALL AKIMA(XIN,YIN,NIN,XASEC(I,ISEC,ISURF),ZC,DUMMY)
-           XLASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF)
-           XUASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF)
-           ZLASEC(I,ISEC,ISURF) = ZC - 0.5*TASEC(I,ISEC,ISURF)
-           ZUASEC(I,ISEC,ISURF) = ZC + 0.5*TASEC(I,ISEC,ISURF)
+        !    XLASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF)
+        !    XUASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF)
+        !    ZLASEC(I,ISEC,ISURF) = ZC - 0.5*TASEC(I,ISEC,ISURF)
+        !    ZUASEC(I,ISEC,ISURF) = ZC + 0.5*TASEC(I,ISEC,ISURF)
            CASEC(I,ISEC,ISURF) = ZC 
 cc#endif
          END DO
diff --git a/src/aio.f b/src/aio.f
index 38716ce..0b24d64 100644
--- a/src/aio.f
+++ b/src/aio.f
@@ -1,6 +1,7 @@
       subroutine write_tecplot(file_name, add_time, solution_time)
 C
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_oml.INC'
 C
       LOGICAL LRANGEALL
       INTEGER LU, nChords, nStrips, nPts
diff --git a/src/amake.f b/src/amake.f
index 6571b47..9729faf 100644
--- a/src/amake.f
+++ b/src/amake.f
@@ -24,6 +24,8 @@ SUBROUTINE MAKESURF(ISURF)
 C     using info from configuration input file.
 C--------------------------------------------------------------
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
+      INCLUDE 'AVL_oml.INC'
 C
 C
       REAL XYZLEL(3), XYZLER(3)
@@ -577,35 +579,43 @@ SUBROUTINE MAKESURF(ISURF)
 C...        nodal grid associated with vortex strip (aft-panel nodes)
 C...        NOTE: airfoil in plane of wing, but not rotated perpendicular to dihedral;
 C...        retained in (x,z) plane at this point
-            CALL AKIMA( XLASEC(1,ISEC,ISURF), ZLASEC(1,ISEC,ISURF), NSL,
-     &                  XPT(IVC+1), ZL_L, DSDX )
-            CALL AKIMA( XUASEC(1,ISEC,ISURF), ZUASEC(1,ISEC,ISURF), NSL,
-     &                  XPT(IVC+1), ZU_L, DSDX )
-C
-            CALL AKIMA( XLASEC(1,ISEC+1,ISURF), ZLASEC(1,ISEC+1,ISURF),
-     &                  NSR, XPT(IVC+1), ZL_R, DSDX )
-            CALL AKIMA( XUASEC(1,ISEC+1,ISURF), ZUASEC(1,ISEC+1,ISURF),
-     &                  NSR, XPT(IVC+1), ZU_R, DSDX )
-C
-            XYN1(1,idx_vor) = RLE1(1,idx_strip) + 
-     &                        XPT(IVC+1)*CHORD1(idx_strip)
-            XYN1(2,idx_vor) = RLE1(2,idx_strip)
+
+            ! recompute 
+            ! XLASEC
+            ! XUASEC
+            ! ZLASEC
+            ! ZUASEC
+
+
+!             CALL AKIMA( XLASEC(1,ISEC,ISURF), ZLASEC(1,ISEC,ISURF), NSL,
+!      &                  XPT(IVC+1), ZL_L, DSDX )
+!             CALL AKIMA( XUASEC(1,ISEC,ISURF), ZUASEC(1,ISEC,ISURF), NSL,
+!      &                  XPT(IVC+1), ZU_L, DSDX )
+! C
+!             CALL AKIMA( XLASEC(1,ISEC+1,ISURF), ZLASEC(1,ISEC+1,ISURF),
+!      &                  NSR, XPT(IVC+1), ZL_R, DSDX )
+!             CALL AKIMA( XUASEC(1,ISEC+1,ISURF), ZUASEC(1,ISEC+1,ISURF),
+!      &                  NSR, XPT(IVC+1), ZU_R, DSDX )
+! C
+!             XYN1(1,idx_vor) = RLE1(1,idx_strip) + 
+!      &                        XPT(IVC+1)*CHORD1(idx_strip)
+!             XYN1(2,idx_vor) = RLE1(2,idx_strip)
             
-            ZL =  (1.-F1)*ZL_L + F1 *ZL_R
-            ZU =  (1.-F1)*ZU_L + F1 *ZU_R
+!             ZL =  (1.-F1)*ZL_L + F1 *ZL_R
+!             ZU =  (1.-F1)*ZU_L + F1 *ZU_R
             
-            ZLON1(idx_vor)  = RLE1(3,idx_strip) + ZL*CHORD1(idx_strip)
-            ZUPN1(idx_vor)  = RLE1(3,idx_strip) + ZU*CHORD1(idx_strip)
-C
-            XYN2(1,idx_vor) = RLE2(1,idx_strip) + 
-     &                        XPT(IVC+1)*CHORD2(idx_strip)
-            XYN2(2,idx_vor) = RLE2(2,idx_strip)
+!             ZLON1(idx_vor)  = RLE1(3,idx_strip) + ZL*CHORD1(idx_strip)
+!             ZUPN1(idx_vor)  = RLE1(3,idx_strip) + ZU*CHORD1(idx_strip)
+! C
+!             XYN2(1,idx_vor) = RLE2(1,idx_strip) + 
+!      &                        XPT(IVC+1)*CHORD2(idx_strip)
+!             XYN2(2,idx_vor) = RLE2(2,idx_strip)
             
-            ZL =  (1.-F2)*ZL_L + F2 *ZL_R
-            ZU =  (1.-F2)*ZU_L + F2 *ZU_R
+!             ZL =  (1.-F2)*ZL_L + F2 *ZL_R
+!             ZU =  (1.-F2)*ZU_L + F2 *ZU_R
           
-            ZLON2(idx_vor)  = RLE2(3,idx_strip) + ZL*CHORD2(idx_strip)
-            ZUPN2(idx_vor)  = RLE2(3,idx_strip) + ZU*CHORD2(idx_strip)
+!             ZLON2(idx_vor)  = RLE2(3,idx_strip) + ZL*CHORD2(idx_strip)
+!             ZUPN2(idx_vor)  = RLE2(3,idx_strip) + ZU*CHORD2(idx_strip)
 C
 cc#endif
             idx_vor = idx_vor + 1
@@ -642,6 +652,7 @@ SUBROUTINE MAKESURF(ISURF)
 
       integer function flatidx(idx_x, idx_y, idx_surf)
       include 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       ! store MFRST and  NVC in the common block
       integer idx_x, idx_y, idx_surf
       flatidx = idx_x + (idx_y - 1) * (NVC(idx_surf)+1)
@@ -655,6 +666,8 @@ subroutine makesurf_mesh(isurf)
 C     and the given mesh coordinate array.
 c--------------------------------------------------------------
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
+      INCLUDE 'AVL_oml.INC'
       ! input/output
       integer isurf
 
@@ -1332,51 +1345,51 @@ subroutine makesurf_mesh(isurf)
         !   xptxind2 = (mesh_surf(1,idx_node_yp1+1)
         !  &           - RLE2(1,idx_strip))/CHORD2(idx_strip) 
 
-          ! Interpolate cross section on left side
-          CALL AKIMA( XLASEC(1,iptl,isurf), ZLASEC(1,iptl,isurf),
-     &                NSL,xptxind1, ZL_L, DSDX )
-          CALL AKIMA( XUASEC(1,iptl,isurf), ZUASEC(1,iptl,isurf),
-     &                NSL,xptxind1, ZU_L, DSDX )
+!           ! Interpolate cross section on left side
+!           CALL AKIMA( XLASEC(1,iptl,isurf), ZLASEC(1,iptl,isurf),
+!      &                NSL,xptxind1, ZL_L, DSDX )
+!           CALL AKIMA( XUASEC(1,iptl,isurf), ZUASEC(1,iptl,isurf),
+!      &                NSL,xptxind1, ZU_L, DSDX )
 
-          ! Interpolate cross section on right side
-          CALL AKIMA(XLASEC(1,iptr,isurf),
-     &               ZLASEC(1,iptr,isurf),NSR, xptxind1, ZL_R, DSDX)
+!           ! Interpolate cross section on right side
+!           CALL AKIMA(XLASEC(1,iptr,isurf),
+!      &               ZLASEC(1,iptr,isurf),NSR, xptxind1, ZL_R, DSDX)
                           
-          CALL AKIMA(XUASEC(1,iptr,isurf),
-     &               ZUASEC(1,iptr,isurf),NSR, xptxind1, ZU_R, DSDX)
+!           CALL AKIMA(XUASEC(1,iptr,isurf),
+!      &               ZUASEC(1,iptr,isurf),NSR, xptxind1, ZU_R, DSDX)
 
 
-          ! Compute the left aft node of panel 
-          ! X-point
-          XYN1(1,idx_vor) = RLE1(1,idx_strip) + 
-     &                      xptxind1*CHORD1(idx_strip)
+!           ! Compute the left aft node of panel 
+!           ! X-point
+!           XYN1(1,idx_vor) = RLE1(1,idx_strip) + 
+!      &                      xptxind1*CHORD1(idx_strip)
 
-          ! Y-point
-          XYN1(2,idx_vor) = RLE1(2,idx_strip)
+!           ! Y-point
+!           XYN1(2,idx_vor) = RLE1(2,idx_strip)
 
-          ! Interpolate z from sections to left aft node of panel
-          ZL =  (1.-f1)*ZL_L + f1 *ZL_R
-          ZU =  (1.-f1)*ZU_L + f1 *ZU_R
+!           ! Interpolate z from sections to left aft node of panel
+!           ZL =  (1.-f1)*ZL_L + f1 *ZL_R
+!           ZU =  (1.-f1)*ZU_L + f1 *ZU_R
 
-          ! Store left aft z-point
-          ZLON1(idx_vor)  = RLE1(3,idx_strip) + ZL*CHORD1(idx_strip)
-          ZUPN1(idx_vor)  = RLE1(3,idx_strip) + ZU*CHORD1(idx_strip)
+!           ! Store left aft z-point
+!           ZLON1(idx_vor)  = RLE1(3,idx_strip) + ZL*CHORD1(idx_strip)
+!           ZUPN1(idx_vor)  = RLE1(3,idx_strip) + ZU*CHORD1(idx_strip)
 
-          ! Compute the right aft node of panel 
-          ! X-point
-          XYN2(1,idx_vor) = RLE2(1,idx_strip) + 
-     &                      xptxind1*CHORD2(idx_strip)
+!           ! Compute the right aft node of panel 
+!           ! X-point
+!           XYN2(1,idx_vor) = RLE2(1,idx_strip) + 
+!      &                      xptxind1*CHORD2(idx_strip)
 
-          ! Y-point
-          XYN2(2,idx_vor) = RLE2(2,idx_strip)
+!           ! Y-point
+!           XYN2(2,idx_vor) = RLE2(2,idx_strip)
                 
-          ! Interpolate z from sections to right aft node of panel
-          ZL =  (1.-f2)*ZL_L + f2 *ZL_R
-          ZU =  (1.-f2)*ZU_L + f2 *ZU_R
+!           ! Interpolate z from sections to right aft node of panel
+!           ZL =  (1.-f2)*ZL_L + f2 *ZL_R
+!           ZU =  (1.-f2)*ZU_L + f2 *ZU_R
 
-          ! Store right aft z-point
-          ZLON2(idx_vor)  = RLE2(3,idx_strip) + ZL*CHORD2(idx_strip)
-          ZUPN2(idx_vor)  = RLE2(3,idx_strip) + ZU*CHORD2(idx_strip)
+!           ! Store right aft z-point
+!           ZLON2(idx_vor)  = RLE2(3,idx_strip) + ZL*CHORD2(idx_strip)
+!           ZUPN2(idx_vor)  = RLE2(3,idx_strip) + ZU*CHORD2(idx_strip)
       
 
           idx_vor = idx_vor + 1
@@ -1417,6 +1430,7 @@ subroutine update_surfaces()
 c--------------------------------------------------------------
       use avl_heap_inc
       include 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       integer ii
       
       NSTRIP = 0
@@ -1476,6 +1490,7 @@ SUBROUTINE MAKEBODY(IBODY)
 C     using info from configuration input file.
 C--------------------------------------------------------------
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
 C
 C      REAL XBOD(IBX), YBOD(IBX), TBOD(IBX)
 C
@@ -1574,6 +1589,7 @@ subroutine update_bodies()
 c     Updates all bodies, using the stored data.
 c--------------------------------------------------------------
       include 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       integer IBODY, NBOD, NBLDS
       character*120 upname
 
@@ -1615,6 +1631,7 @@ subroutine  set_section_coordinates(isec,isurf,x,y,n,nin,xfmin,
 c     Sets the airfoil coodinate data for the given section and surface
 c--------------------------------------------------------------
       include 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
 c      input
       integer isec, isurf, n, nin
       real x(n), y(n)
@@ -1658,10 +1675,10 @@ subroutine  set_section_coordinates(isec,isurf,x,y,n,nin,xfmin,
      &               SASEC(i,isec,isurf))
           call AKIMA(xin,tin,nin,XASEC(i,isec,isurf),
      &               TASEC(i,isec,isurf),dummy)
-          XLASEC(i,isec,isurf) = XASEC(i,isec,isurf)
-          XUASEC(i,isec,isurf) = XASEC(i,isec,isurf)
-          ZLASEC(i,isec,isurf) = zc - 0.5*TASEC(i,isec,isurf)
-          ZUASEC(i,isec,isurf) = zc + 0.5*TASEC(i,isec,isurf)
+      !     XLASEC(i,isec,isurf) = XASEC(i,isec,isurf)
+      !     XUASEC(i,isec,isurf) = XASEC(i,isec,isurf)
+      !     ZLASEC(i,isec,isurf) = zc - 0.5*TASEC(i,isec,isurf)
+      !     ZUASEC(i,isec,isurf) = zc + 0.5*TASEC(i,isec,isurf)
           CASEC(i,isec,isurf) = zc
 
         end do
@@ -1674,6 +1691,7 @@ subroutine set_body_coordinates(ibod,xb,yb,nb,nin,storecoords)
 c     Sets the body oml coodinate data for the given section and surface
 c--------------------------------------------------------------
       include 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       integer ibod, nb, nin
       real xb(nb), yb(nb)
       logical storecoords
@@ -1702,6 +1720,8 @@ SUBROUTINE SDUPL(NN, Ypt,MSG)
 C     reflected about y=Ypt.
 C-----------------------------------
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
+      INCLUDE 'AVL_oml.INC'
       CHARACTER*(*) MSG
       integer idx_vor
 C
@@ -1923,6 +1943,7 @@ SUBROUTINE BDUPL(NN,Ypt,MSG)
 C     reflected about y=Ypt.
 C-----------------------------------
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
       CHARACTER*(*) MSG
 C
       NNI = NBODY + 1
@@ -2009,6 +2030,7 @@ SUBROUTINE ENCALC
 C...COMMENTS   
 C
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_surf.INC'
 C
       REAL EP(3), EQ(3), ES(3), EB(3), EC(3), ECXB(3)
       REAL EC_G(3,NDMAX), ECXB_G(3)
diff --git a/src/aoml.f b/src/aoml.f
index 6bfbbd2..8d8dce3 100644
--- a/src/aoml.f
+++ b/src/aoml.f
@@ -23,6 +23,7 @@
       SUBROUTINE CPOML(save_file)
 C
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_oml.INC'
 C
       logical :: save_file
       LOGICAL LRANGEALL
@@ -59,6 +60,7 @@ SUBROUTINE CPTHK
 C...COMMENTS
 C
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_oml.INC'
 C
       PARAMETER (NCMAX=256)
       DIMENSION AICT(NCMAX,NCMAX)       ! AIC for no normal flow BC
@@ -186,6 +188,7 @@ SUBROUTINE CPDUMP(save_file)
 C...COMMENTS    C code reader in read_cpoml.c
 C
       INCLUDE 'AVL.INC'
+      INCLUDE 'AVL_oml.INC'
       logical :: save_file
       
       integer idx_mesh, idx_mesh_surf, idx_LE, idx_LE_lo, idx_LE_up
diff --git a/src/asetup.f b/src/asetup.f
index 90df01c..8ce5400 100644
--- a/src/asetup.f
+++ b/src/asetup.f
@@ -520,34 +520,34 @@ SUBROUTINE VELSUM()
 C
       DO I = 1, NVOR
         DO K = 1, 3
-          VC(K,I) = 0.0
+          ! VC(K,I) = 0.0
           VV(K,I) = 0.0
 C--- h.v. velocity at control points and vortex midpoints
           DO J = 1, NVOR
-            VC(K,I) = VC(K,I) + WC_GAM(K,I,J)*GAM(J)
+            ! VC(K,I) = VC(K,I) + WC_GAM(K,I,J)*GAM(J)
             VV(K,I) = VV(K,I) + WV_GAM(K,I,J)*GAM(J)
           ENDDO
           DO N = 1, NUMAX
-            VC_U(K,I,N) = 0.0
+            ! VC_U(K,I,N) = 0.0
             VV_U(K,I,N) = 0.0
             DO J = 1, NVOR
-              VC_U(K,I,N) = VC_U(K,I,N) + WC_GAM(K,I,J)*GAM_U(J,N)
+              ! VC_U(K,I,N) = VC_U(K,I,N) + WC_GAM(K,I,J)*GAM_U(J,N)
               VV_U(K,I,N) = VV_U(K,I,N) + WV_GAM(K,I,J)*GAM_U(J,N)
             ENDDO
           ENDDO
           DO N = 1, NCONTROL
-            VC_D(K,I,N) = 0.0
+            ! VC_D(K,I,N) = 0.0
             VV_D(K,I,N) = 0.0
             DO J = 1, NVOR
-              VC_D(K,I,N) = VC_D(K,I,N) + WC_GAM(K,I,J)*GAM_D(J,N)
+              ! VC_D(K,I,N) = VC_D(K,I,N) + WC_GAM(K,I,J)*GAM_D(J,N)
               VV_D(K,I,N) = VV_D(K,I,N) + WV_GAM(K,I,J)*GAM_D(J,N)
             ENDDO
           ENDDO
           DO N = 1, NDESIGN
-            VC_G(K,I,N) = 0.0
+            ! VC_G(K,I,N) = 0.0
             VV_G(K,I,N) = 0.0
             DO J = 1, NVOR
-              VC_G(K,I,N) = VC_G(K,I,N) + WC_GAM(K,I,J)*GAM_G(J,N)
+              ! VC_G(K,I,N) = VC_G(K,I,N) + WC_GAM(K,I,J)*GAM_G(J,N)
               VV_G(K,I,N) = VV_G(K,I,N) + WV_GAM(K,I,J)*GAM_G(J,N)
             ENDDO
           ENDDO
@@ -565,18 +565,18 @@ SUBROUTINE VELSUM()
      &               + WVSRD_U(K,I,5)*WROT(2)
      &               + WVSRD_U(K,I,6)*WROT(3)
 C--- total velocity at control points and vortex midpoints
-          WC(K,I) = VC(K,I) + WCSRD(K,I)
+          ! WC(K,I) = VC(K,I) + WCSRD(K,I)
           WV(K,I) = VV(K,I) + WVSRD(K,I)
           DO N = 1, NUMAX
-            WC_U(K,I,N) = VC_U(K,I,N) + WCSRD_U(K,I,N)
+            ! WC_U(K,I,N) = VC_U(K,I,N) + WCSRD_U(K,I,N)
             WV_U(K,I,N) = VV_U(K,I,N) + WVSRD_U(K,I,N)
           ENDDO
           DO N = 1, NDMAX
-            WC_D(K,I,N) = VC_D(K,I,N)
+            ! WC_D(K,I,N) = VC_D(K,I,N)
             WV_D(K,I,N) = VV_D(K,I,N)
           ENDDO
           DO N = 1, NGMAX
-            WC_G(K,I,N) = VC_G(K,I,N)
+            ! WC_G(K,I,N) = VC_G(K,I,N)
             WV_G(K,I,N) = VV_G(K,I,N)
           ENDDO
 C
diff --git a/src/f2py/libavl.pyf b/src/f2py/libavl.pyf
index c5c9157..d847c0f 100644
--- a/src/f2py/libavl.pyf
+++ b/src/f2py/libavl.pyf
@@ -8,6 +8,8 @@ python module libavl ! in
 
         subroutine avl ! in :libavl:avl.f
             include '../includes/AVL.INC'
+            include '../includes/AVL_surf.INC'
+            include '../includes/AVL_oml.INC'
             include '../includes/AVL_ad_seeds.inc'
         
             ! parameters from 
diff --git a/src/includes/ADIMEN.INC b/src/includes/ADIMEN.INC
index cf2a343..3a73bb6 100644
--- a/src/includes/ADIMEN.INC
+++ b/src/includes/ADIMEN.INC
@@ -32,11 +32,35 @@
 !              to make it easier to termine the true size in clear_ad_seeds_fast()
 !              I increased some values to avoid matching
 
+      ! ! ! PARAMETER (NVMAX=4000)
+      ! PARAMETER (NVMAX=5000)
+
+      ! ! PARAMETER (NSMAX=400) ! max number of strips
+      ! ! PARAMETER (NSECMAX=401) ! max number of geometry sections
+      ! PARAMETER (NSMAX=500) ! max number of strips
+      ! PARAMETER (NSECMAX=301) ! max number of geometry sections
+      ! PARAMETER (NFMAX=100) ! max number of surfaces
+      ! PARAMETER (NLMAX=502) ! max number of body nodes
+      ! PARAMETER (NBMAX=20)  ! max number of bodies
+
+      ! PARAMETER (NUMAX=6)  ! number of velocity components
+      ! PARAMETER (NDMAX=30) ! max number of control surfaces
+      ! PARAMETER (NGMAX=21) ! max number of design variables
+
+      ! PARAMETER (NRMAX=25)  ! max number of run cases
+      ! PARAMETER (NTMAX=503) ! max number of time values
+
+      ! PARAMETER (NOBMAX=1) ! max number of off body points
+      ! ! PARAMETER (NOBMAX=1000) ! max number of off body points
+      
+      ! PARAMETER (ICONX=20) ! 
+      ! PARAMETER (IBX=200)  ! max number of airfoil coordinates
+      ! ! PARAMETER (IBX=300)  ! max number of airfoil coordinates
+
+
       ! PARAMETER (NVMAX=4000)
       PARAMETER (NVMAX=5000)
 
-      ! PARAMETER (NSMAX=400) ! max number of strips
-      ! PARAMETER (NSECMAX=401) ! max number of geometry sections
       PARAMETER (NSMAX=500) ! max number of strips
       PARAMETER (NSECMAX=301) ! max number of geometry sections
       PARAMETER (NFMAX=100) ! max number of surfaces
@@ -51,9 +75,9 @@
       PARAMETER (NTMAX=503) ! max number of time values
 
       PARAMETER (NOBMAX=1) ! max number of off body points
-      ! PARAMETER (NOBMAX=1000) ! max number of off body points
+      ! PARAMETER (NOBMAX=1) ! max number of off body points
       
-      PARAMETER (ICONX=20) ! 
+      PARAMETER (ICONX=20) ! max number of control or design variable declaration lines per section 
       PARAMETER (IBX=200)  ! max number of airfoil coordinates
       ! PARAMETER (IBX=300)  ! max number of airfoil coordinates
 
diff --git a/src/includes/AVL.INC b/src/includes/AVL.INC
index fc47c9c..e042fda 100644
--- a/src/includes/AVL.INC
+++ b/src/includes/AVL.INC
@@ -44,10 +44,10 @@ C
      & FEVDEF,         ! default eigenvalue save file
      & TITLE,          ! configuration title
      & STITLE(NFMAX),  ! surface title
-     & AFILES(NSMAX,NFMAX),  ! airfoil file names
+     & AFILES(NSECMAX,NFMAX),  ! airfoil file names
      & BTITLE(NBMAX),  ! body title
      & BFILES(NBMAX),  ! body file names
-     & NACA(NSMAX,NFMAX),  ! naca section names
+     & NACA(NSECMAX,NFMAX),  ! naca section names
      & RTITLE(NRMAX),  ! run case title
      & DNAME(NDMAX),   ! control variable name
      & GNAME(NGMAX),   ! design  variable name
@@ -395,144 +395,6 @@ C
      & CF_LSRF(3,NFMAX), CM_LSRF(3,NFMAX),
      & CL_LSRF(NFMAX), CD_LSRF(NFMAX), CMLE_LSRF(NFMAX)
 
-C      !start added variables for python geometry minipulation
-
-      LOGICAL LDUPL_B
-      COMMON /BODY_GEOM_L/
-     & LDUPL_B(NBMAX)      ! T if body is a duplicated
-c
-      COMMON /BODY_GEOM_I/
-     & NSEC_B(NBMAX),  ! number of sections in body
-     & NVB(NFMAX),     ! number of body elements
-     & NBOD(NBMAX)    ! number of body oml points
-
-      REAL(kind=avl_real) XYZSCAL_B
-      REAL(kind=avl_real) XYZTRAN_B
-      REAL(kind=avl_real) YDUPL_B
-      REAL(kind=avl_real) XYZLES_B 
-      REAL(kind=avl_real) BSPACE
-      REAL(kind=avl_real) XBOD
-      REAL(kind=avl_real) YBOD
-      REAL(kind=avl_real) TBOD
-      REAL(kind=avl_real) XBOD_R
-      REAL(kind=avl_real) YBOD_R
-      COMMON /BODY_GEOM_R/
-     & XYZSCAL_B(3, NBMAX), ! scaling factors for XYZ coordinates
-     & XYZTRAN_B(3, NBMAX), ! translation factors for XYZ coordinates
-     & YDUPL_B(NBMAX), ! y duplicate
-     & XYZLES_B(3, NSMAX, NBMAX), ! leading edge cordinate vector
-     & BSPACE(NBMAX), ! body spacing
-     & XBOD(IBX, NBMAX), ! body oml x-coordinates
-     & YBOD(IBX, NBMAX), ! body oml y-coordinates
-     & TBOD(IBX, NBMAX), ! body oml thickness
-     & XBOD_R(IBX, NBMAX), ! raw input oml x-coordinates
-     & YBOD_R(IBX, NBMAX) ! raw input oml y-coordinates
-c
-      LOGICAL LDUPL
-      LOGICAL LSURFSPACING
-      COMMON /SURF_GEOM_L/
-     & LDUPL(NFMAX),      ! T if surface is a duplicated
-     & LSURFSPACING(NFMAX)   ! surface spacing set under the surface heeading
-
-      COMMON /SURF_GEOM_I/
-     & NSEC(NFMAX),  ! number of sections in surface
-     & NVC(NFMAX), ! number of chordwise elements
-     & NVS(NFMAX), ! number of spanwise elements
-     & NSPANS(NSECMAX, NFMAX), ! number of spanwise elements vector
-     & NASEC(NSECMAX, NFMAX), ! number of points used to represent section geometry
-     & NAPTSSEC(NSECMAX, IBX), ! number of points of input airfoil geometry
-     & ICONTD(ICONX, NSECMAX, NFMAX), ! control variable index
-     & NSCON(NSECMAX, NFMAX), ! number of control variables
-     & IDESTD(ICONX, NSECMAX, NFMAX), ! design variable index
-     & NSDES(NSECMAX, NFMAX) ! number of design variables
-
-      REAL(kind=avl_real) XYZSCAL
-      REAL(kind=avl_real) XYZTRAN
-      REAL(kind=avl_real) YDUPL
-      REAL(kind=avl_real) ADDINC
-      REAL(kind=avl_real) CSPACE
-      REAL(kind=avl_real) SSPACE
-      REAL(kind=avl_real) SSPACES
-      REAL(kind=avl_real) XYZLES
-      REAL(kind=avl_real) XSEC
-      REAL(kind=avl_real) YSEC
-      REAL(kind=avl_real) XFMIN_R
-      REAL(kind=avl_real) XFMAX_R
-      REAL(kind=avl_real) XLASEC
-      REAL(kind=avl_real) ZLASEC
-      REAL(kind=avl_real) XUASEC
-      REAL(kind=avl_real) ZUASEC
-      REAL(kind=avl_real) CHORDS
-      REAL(kind=avl_real) AINCS
-      REAL(kind=avl_real) XASEC
-      REAL(kind=avl_real) CASEC
-      REAL(kind=avl_real) SASEC
-      REAL(kind=avl_real) TASEC
-      REAL(kind=avl_real) CLCDSEC
-      REAL(kind=avl_real) CLCDSRF
-      REAL(kind=avl_real) CLAF
-      REAL(kind=avl_real) XHINGED
-      REAL(kind=avl_real) VHINGED
-      REAL(kind=avl_real) GAIND
-      REAL(kind=avl_real) REFLD
-      REAL(kind=avl_real) GAING  
-      COMMON /SURF_GEOM_R/
-     & XYZSCAL(3, NFMAX), ! scaling factors for XYZ coordinates
-     & XYZTRAN(3, NFMAX), ! translation factors for XYZ coordinates
-     & YDUPL(NFMAX), ! y position of duplicate plane
-     & ADDINC(NFMAX), ! additional incidence angle
-     & CSPACE(NFMAX), ! chordwise spacing
-     & SSPACE(NFMAX), ! spanwise spacing
-     & SSPACES(NSECMAX, NFMAX), ! spanwise spacing vector
-     & XYZLES(3, NSECMAX, NFMAX), ! leading edge cordinate vector
-     & XSEC(IBX, NSECMAX, NFMAX), ! raw section input x coordinates
-     & YSEC(IBX, NSECMAX, NFMAX), ! raw section input y coordinates
-     & XFMIN_R(NSECMAX, NFMAX), ! section input lower x bound 
-     & XFMAX_R(NSECMAX, NFMAX), ! section input upper x bound 
-     & XLASEC(IBX, NSECMAX, NFMAX), ! airfoil lower surface x-coords
-     & ZLASEC(IBX, NSECMAX, NFMAX), ! airfoil lower surface z-coords
-     & XUASEC(IBX, NSECMAX, NFMAX), ! airfoil upper surface x-coords
-     & ZUASEC(IBX, NSECMAX, NFMAX), ! airfoil upper surface z-coords
-     & CHORDS(NSECMAX, NFMAX), ! chord length vectorm
-     & AINCS(NSECMAX, NFMAX), ! incidence angle vector
-     & XASEC(IBX, NSECMAX, NFMAX),  ! the x coordinate aifoil section
-     & CASEC(IBX, NSECMAX, NFMAX),  ! camber line at xasec
-     & SASEC(IBX, NSECMAX, NFMAX),  ! slope of camber line at xasec
-     & TASEC(IBX, NSECMAX, NFMAX),   ! thickness at xasec
-     & CLCDSEC(6, NSECMAX, NFMAX),
-     & CLCDSRF(6, NFMAX),
-     & CLAF(NSECMAX, NFMAX),   ! CL alpha (dCL/da) scaling factor
-     & XHINGED(ICONX, NSECMAX, NFMAX), ! hinge location
-     & VHINGED(3, ICONX, NSECMAX, NFMAX), ! hinge vector
-     & GAIND(ICONX, NSECMAX, NFMAX), ! control surface gain
-     & REFLD(ICONX, NSECMAX, NFMAX), ! control surface reflection
-     & GAING(ICONX, NSECMAX, NFMAX)  ! desgin variable gain
-
-      COMMON /SURF_MESH_I/
-     & MFRST(NFMAX), ! stores the index in the MSHBLK where each surface's mesh begins
-     & IPTSEC(NSMAX,NFMAX) ! stores the iptloc vector for each surface
-
-      REAL(kind=avl_real) MSHBLK
-      REAL(kind=avl_real) RV1MSH
-      REAL(kind=avl_real) RV2MSH
-      REAL(kind=avl_real) RVMSH
-      REAL(kind=avl_real) RCMSH
-      COMMON /SURF_MESH_R/
-     & MSHBLK(3, 4*NVMAX), ! block to store all surface meshes
-     & RV1MSH(3,NVMAX),    !  mesh h.v. vortex left  points
-     & RV2MSH(3,NVMAX),    !  mesh h.v. vortex right  points
-     & RVMSH(3,NVMAX),    ! mesh h.v. vortex center points
-     & RCMSH(3,NVMAX)    ! mesh h.v. control points
-
-      LOGICAL LSURFMSH
-      LOGICAL LMESHFLAT
-      COMMON /SURF_MESH_L/
-     & LSURFMSH(NFMAX), ! T if surface uses a mesh cordinates for its geometry
-     & LMESHFLAT(NFMAX) ! T if the surface's mesh should be flattened when computing vorticies and control points
-
-
-C      !!--- end added variables for python geometry manipulation ---
-
       COMMON /STRP_I/
      & LSSURF(NSMAX),    ! index of surface which contains this strip
      & IJFRST(NSMAX),    ! index of first element in strip
@@ -738,22 +600,6 @@ cc#else
      & RES_U(NVMAX,NUMAX),               ! residual for gamma_u lin sys
      & RES_D(NVMAX,NDMAX)        ! residual for gamma_d lin. sys 
 
-cc#ifdef USE_CPOML
-      REAL(kind=avl_real) XYN1
-      REAL(kind=avl_real) XYN2
-      REAL(kind=avl_real) ZLON1
-      REAL(kind=avl_real) ZLON2
-      REAL(kind=avl_real) ZUPN1
-      REAL(kind=avl_real) ZUPN2
-      REAL(kind=avl_real) XYZSURF
-      REAL(kind=avl_real) CPSURF
-      COMMON /VRTX_S/
-     & XYN1(2,NVMAX), XYN2(2,NVMAX),  ! left/right aft-node of element
-     & ZLON1(NVMAX), ZLON2(NVMAX),    ! left/right lower z-coord of aft-node
-     & ZUPN1(NVMAX), ZUPN2(NVMAX),     ! left/right upper z-coord of aft-node
-     & XYZSURF(3, NVMAX, NFMAX),
-     & CPSURF(NVMAX,NFMAX)
-cc#endif
 C
       COMMON /BODY_I/
      & NL(NBMAX),       ! number of source-line nodes in body
@@ -782,13 +628,13 @@ C
 
       REAL(kind=avl_real) AMACH
       REAL(kind=avl_real) VC
-      REAL(kind=avl_real) VC_U
-      REAL(kind=avl_real) VC_D
-      REAL(kind=avl_real) VC_G
-      REAL(kind=avl_real) WC
-      REAL(kind=avl_real) WC_U
-      REAL(kind=avl_real) WC_D
-      REAL(kind=avl_real) WC_G
+c      REAL(kind=avl_real) VC_U
+c      REAL(kind=avl_real) VC_D
+c      REAL(kind=avl_real) VC_G
+c      REAL(kind=avl_real) WC
+c      REAL(kind=avl_real) WC_U
+c      REAL(kind=avl_real) WC_D
+c      REAL(kind=avl_real) WC_G
       REAL(kind=avl_real) WV
       REAL(kind=avl_real) WV_U
       REAL(kind=avl_real) WV_D
@@ -805,14 +651,14 @@ cc     & AICN(NVMAX,NVMAX),     ! normalwash AIC matrix (and VL system matrix)
 cc     & WC_GAM(3,NVMAX,NVMAX), ! c.p. velocity/Gamma influence matrix
 cc     & WV_GAM(3,NVMAX,NVMAX), ! h.v. velocity/Gamma influence matrix
      & VC(3,NVMAX),           ! h.v. induced velocity at c.p.
-     & VC_U(3,NVMAX,NUMAX),
-     & VC_D(3,NVMAX,NDMAX),
-     & VC_G(3,NVMAX,NGMAX),
+c     & VC_U(3,NVMAX,NUMAX),
+c     & VC_D(3,NVMAX,NDMAX),
+c     & VC_G(3,NVMAX,NGMAX),
 C
-     & WC(3,NVMAX),           ! total induced velocity at c.p. (used in residual)
-     & WC_U(3,NVMAX,NUMAX),
-     & WC_D(3,NVMAX,NDMAX),
-     & WC_G(3,NVMAX,NGMAX),
+c     & WC(3,NVMAX),           ! total induced velocity at c.p. (used in residual)
+C    & WC_U(3,NVMAX,NUMAX),
+C    & WC_D(3,NVMAX,NDMAX),
+C    & WC_G(3,NVMAX,NGMAX),
 C
      & VV(3,NVMAX),           ! h.v. induced velocity at h.v.
      & VV_U(3,NVMAX,NUMAX),
diff --git a/src/includes/AVL_ad_seeds.inc b/src/includes/AVL_ad_seeds.inc
index dae2b9f..6b0d6e5 100644
--- a/src/includes/AVL_ad_seeds.inc
+++ b/src/includes/AVL_ad_seeds.inc
@@ -390,28 +390,28 @@ C
       real(kind=avl_real) SSPACE_DIFF
       real(kind=avl_real) SSPACES_DIFF
       real(kind=avl_real) XYZLES_DIFF
-      real(kind=avl_real) XSEC_DIFF
-      real(kind=avl_real) YSEC_DIFF
+c      real(kind=avl_real) XSEC_DIFF
+c      real(kind=avl_real) YSEC_DIFF
       real(kind=avl_real) XFMIN_R_DIFF
       real(kind=avl_real) XFMAX_R_DIFF
-      real(kind=avl_real) XLASEC_DIFF
-      real(kind=avl_real) ZLASEC_DIFF
-      real(kind=avl_real) XUASEC_DIFF
-      real(kind=avl_real) ZUASEC_DIFF
+c      real(kind=avl_real) XLASEC_DIFF
+c      real(kind=avl_real) ZLASEC_DIFF
+c      real(kind=avl_real) XUASEC_DIFF
+c      real(kind=avl_real) ZUASEC_DIFF
       real(kind=avl_real) CHORDS_DIFF
       real(kind=avl_real) AINCS_DIFF
       real(kind=avl_real) XASEC_DIFF
-      real(kind=avl_real) CASEC_DIFF
+c      real(kind=avl_real) CASEC_DIFF
       real(kind=avl_real) SASEC_DIFF
-      real(kind=avl_real) TASEC_DIFF
-      real(kind=avl_real) CLCDSEC_DIFF
+c      real(kind=avl_real) TASEC_DIFF
+c      real(kind=avl_real) CLCDSEC_DIFF
       real(kind=avl_real) CLCDSRF_DIFF
       real(kind=avl_real) CLAF_DIFF
-      real(kind=avl_real) XHINGED_DIFF
-      real(kind=avl_real) VHINGED_DIFF
-      real(kind=avl_real) GAIND_DIFF
-      real(kind=avl_real) REFLD_DIFF
-      real(kind=avl_real) GAING_DIFF
+c      real(kind=avl_real) XHINGED_DIFF
+c      real(kind=avl_real) VHINGED_DIFF
+c      real(kind=avl_real) GAIND_DIFF
+c      real(kind=avl_real) REFLD_DIFF
+c      real(kind=avl_real) GAING_DIFF
       COMMON /SURF_GEOM_R_DIFF/
      & XYZSCAL_DIFF(3, NFMAX),
      & XYZTRAN_DIFF(3, NFMAX),
@@ -421,28 +421,28 @@ C
      & SSPACE_DIFF(NFMAX),
      & SSPACES_DIFF(NSECMAX, NFMAX),
      & XYZLES_DIFF(3, NSECMAX, NFMAX),
-     & XSEC_DIFF(IBX, NSECMAX, NFMAX),
-     & YSEC_DIFF(IBX, NSECMAX, NFMAX),
+c     & XSEC_DIFF(IBX, NSECMAX, NFMAX),
+c     & YSEC_DIFF(IBX, NSECMAX, NFMAX),
      & XFMIN_R_DIFF(NSECMAX, NFMAX),
      & XFMAX_R_DIFF(NSECMAX, NFMAX),
-     & XLASEC_DIFF(IBX, NSECMAX, NFMAX),
-     & ZLASEC_DIFF(IBX, NSECMAX, NFMAX),
-     & XUASEC_DIFF(IBX, NSECMAX, NFMAX),
-     & ZUASEC_DIFF(IBX, NSECMAX, NFMAX),
+c     & XLASEC_DIFF(IBX, NSECMAX, NFMAX),
+c     & ZLASEC_DIFF(IBX, NSECMAX, NFMAX),
+c     & XUASEC_DIFF(IBX, NSECMAX, NFMAX),
+c     & ZUASEC_DIFF(IBX, NSECMAX, NFMAX),
      & CHORDS_DIFF(NSECMAX, NFMAX),
      & AINCS_DIFF(NSECMAX, NFMAX),
-     & XASEC_DIFF(IBX, NSECMAX, NFMAX),
-     & CASEC_DIFF(IBX, NSECMAX, NFMAX),
-     & SASEC_DIFF(IBX, NSECMAX, NFMAX),
-     & TASEC_DIFF(IBX, NSECMAX, NFMAX),
-     & CLCDSEC_DIFF(6, NSECMAX, NFMAX),
+     & XASEC_DIFF(50, NSECMAX, NFMAX),
+c     & CASEC_DIFF(IBX, NSECMAX, NFMAX),
+     & SASEC_DIFF(50, NSECMAX, NFMAX),
+c     & TASEC_DIFF(IBX, NSECMAX, NFMAX),
+c     & CLCDSEC_DIFF(6, NSECMAX, NFMAX),
      & CLCDSRF_DIFF(6, NFMAX),
-     & CLAF_DIFF(NSECMAX, NFMAX),
-     & XHINGED_DIFF(ICONX, NSECMAX, NFMAX),
-     & VHINGED_DIFF(3, ICONX, NSECMAX, NFMAX),
-     & GAIND_DIFF(ICONX, NSECMAX, NFMAX),
-     & REFLD_DIFF(ICONX, NSECMAX, NFMAX),
-     & GAING_DIFF(ICONX, NSECMAX, NFMAX)
+     & CLAF_DIFF(NSECMAX, NFMAX)
+c     & XHINGED_DIFF(ICONX, NSECMAX, NFMAX),
+c     & VHINGED_DIFF(3, ICONX, NSECMAX, NFMAX),
+c     & GAIND_DIFF(ICONX, NSECMAX, NFMAX),
+c     & REFLD_DIFF(ICONX, NSECMAX, NFMAX),
+c     & GAING_DIFF(ICONX, NSECMAX, NFMAX)
 C
       real(kind=avl_real) MSHBLK_DIFF
       real(kind=avl_real) RV1MSH_DIFF
@@ -594,18 +594,18 @@ C
       real(kind=avl_real) VREFL_DIFF
       real(kind=avl_real) ENC_DIFF
       real(kind=avl_real) ENV_DIFF
-      real(kind=avl_real) ENC_D_DIFF
-      real(kind=avl_real) ENC_G_DIFF
-      real(kind=avl_real) ENV_D_DIFF
-      real(kind=avl_real) ENV_G_DIFF
+c      real(kind=avl_real) ENC_D_DIFF
+c      real(kind=avl_real) ENC_G_DIFF
+c      real(kind=avl_real) ENV_D_DIFF
+c      real(kind=avl_real) ENV_G_DIFF
       real(kind=avl_real) DCP_DIFF
       real(kind=avl_real) DCP_U_DIFF
       real(kind=avl_real) DCP_D_DIFF
       real(kind=avl_real) DCP_G_DIFF
       real(kind=avl_real) GAM_DIFF
       real(kind=avl_real) GAM_U_0_DIFF
-      real(kind=avl_real) GAM_U_D_DIFF
-      real(kind=avl_real) GAM_U_G_DIFF
+c      real(kind=avl_real) GAM_U_D_DIFF
+C      real(kind=avl_real) GAM_U_G_DIFF
       real(kind=avl_real) GAM_U_DIFF
       real(kind=avl_real) GAM_D_DIFF
       real(kind=avl_real) GAM_G_DIFF
@@ -643,17 +643,17 @@ C
      & ENC_DIFF(3,NVMAX),
      & ENV_DIFF(3,NVMAX),
      & ENC_D_DIFF(3,NVMAX,NDMAX),
-     & ENC_G_DIFF(3,NVMAX,NGMAX),
-     & ENV_D_DIFF(3,NVMAX,NDMAX),
-     & ENV_G_DIFF(3,NVMAX,NGMAX),
+c     & ENC_G_DIFF(3,NVMAX,NGMAX),
+c     & ENV_D_DIFF(3,NVMAX,NDMAX),
+c     & ENV_G_DIFF(3,NVMAX,NGMAX),
      & DCP_DIFF(NVMAX),
      & DCP_U_DIFF(NVMAX,NUMAX),
      & DCP_D_DIFF(NVMAX,NDMAX),
      & DCP_G_DIFF(NVMAX,NGMAX),
      & GAM_DIFF(NVMAX),
      & GAM_U_0_DIFF(NVMAX,NUMAX),
-     & GAM_U_D_DIFF(NVMAX,NUMAX,NDMAX),
-     & GAM_U_G_DIFF(NVMAX,NUMAX,NGMAX),
+c     & GAM_U_D_DIFF(NVMAX,NUMAX,NDMAX),
+c      & GAM_U_G_DIFF(NVMAX,NUMAX,NGMAX),
      & GAM_U_DIFF(NVMAX,NUMAX),
      & GAM_D_DIFF(NVMAX,NDMAX),
      & GAM_G_DIFF(NVMAX,NGMAX),
@@ -697,38 +697,38 @@ C
       real(kind=avl_real) AMACH_DIFF
       real(kind=avl_real) VC_DIFF
       real(kind=avl_real) VC_U_DIFF
-      real(kind=avl_real) VC_D_DIFF
-      real(kind=avl_real) VC_G_DIFF
+c      real(kind=avl_real) VC_D_DIFF
+c      real(kind=avl_real) VC_G_DIFF
       real(kind=avl_real) WC_DIFF
       real(kind=avl_real) WC_U_DIFF
-      real(kind=avl_real) WC_D_DIFF
-      real(kind=avl_real) WC_G_DIFF
+c      real(kind=avl_real) WC_D_DIFF
+c      real(kind=avl_real) WC_G_DIFF
       real(kind=avl_real) VV_DIFF
       real(kind=avl_real) VV_U_DIFF
       real(kind=avl_real) VV_D_DIFF
-      real(kind=avl_real) VV_G_DIFF
+c      real(kind=avl_real) VV_G_DIFF
       real(kind=avl_real) WV_DIFF
       real(kind=avl_real) WV_U_DIFF
       real(kind=avl_real) WV_D_DIFF
-      real(kind=avl_real) WV_G_DIFF
+c      real(kind=avl_real) WV_G_DIFF
       COMMON /SOLV_R_DIFF/
      & AMACH_DIFF,
      & VC_DIFF(3,NVMAX),
      & VC_U_DIFF(3,NVMAX,NUMAX),
-     & VC_D_DIFF(3,NVMAX,NDMAX),
-     & VC_G_DIFF(3,NVMAX,NGMAX),
+c     & VC_D_DIFF(3,NVMAX,NDMAX),
+c     & VC_G_DIFF(3,NVMAX,NGMAX),
      & WC_DIFF(3,NVMAX),
      & WC_U_DIFF(3,NVMAX,NUMAX),
-     & WC_D_DIFF(3,NVMAX,NDMAX),
-     & WC_G_DIFF(3,NVMAX,NGMAX),
+c     & WC_D_DIFF(3,NVMAX,NDMAX),
+c     & WC_G_DIFF(3,NVMAX,NGMAX),
      & VV_DIFF(3,NVMAX),
      & VV_U_DIFF(3,NVMAX,NUMAX),
      & VV_D_DIFF(3,NVMAX,NDMAX),
-     & VV_G_DIFF(3,NVMAX,NGMAX),
+c     & VV_G_DIFF(3,NVMAX,NGMAX),
      & WV_DIFF(3,NVMAX),
      & WV_U_DIFF(3,NVMAX,NUMAX),
-     & WV_D_DIFF(3,NVMAX,NDMAX),
-     & WV_G_DIFF(3,NVMAX,NGMAX)
+     & WV_D_DIFF(3,NVMAX,NDMAX)
+c     & WV_G_DIFF(3,NVMAX,NGMAX)
 C
       real(kind=avl_real) ROB_DIFF
       real(kind=avl_real) VOB_DIFF
diff --git a/src/includes/AVL_oml.INC b/src/includes/AVL_oml.INC
new file mode 100644
index 0000000..6cbc72a
--- /dev/null
+++ b/src/includes/AVL_oml.INC
@@ -0,0 +1,16 @@
+cc#ifdef USE_CPOML
+      REAL(kind=avl_real) XYN1
+      REAL(kind=avl_real) XYN2
+      REAL(kind=avl_real) ZLON1
+      REAL(kind=avl_real) ZLON2
+      REAL(kind=avl_real) ZUPN1
+      REAL(kind=avl_real) ZUPN2
+      REAL(kind=avl_real) XYZSURF
+      REAL(kind=avl_real) CPSURF
+      COMMON /VRTX_S/
+     & XYN1(2,NVMAX), XYN2(2,NVMAX),  ! left/right aft-node of element
+     & ZLON1(NVMAX), ZLON2(NVMAX),    ! left/right lower z-coord of aft-node
+     & ZUPN1(NVMAX), ZUPN2(NVMAX),     ! left/right upper z-coord of aft-node
+     & XYZSURF(3, NVMAX, NFMAX),
+     & CPSURF(NVMAX,NFMAX)
+cc#endif
\ No newline at end of file
diff --git a/src/includes/AVL_surf.INC b/src/includes/AVL_surf.INC
new file mode 100644
index 0000000..f7eb152
--- /dev/null
+++ b/src/includes/AVL_surf.INC
@@ -0,0 +1,137 @@
+C      !start added variables for python geometry minipulation
+
+      LOGICAL LDUPL_B
+      COMMON /BODY_GEOM_L/
+     & LDUPL_B(NBMAX)      ! T if body is a duplicated
+c
+      COMMON /BODY_GEOM_I/
+     & NSEC_B(NBMAX),  ! number of sections in body
+     & NVB(NFMAX),     ! number of body elements
+     & NBOD(NBMAX)    ! number of body oml points
+
+      REAL(kind=avl_real) XYZSCAL_B
+      REAL(kind=avl_real) XYZTRAN_B
+      REAL(kind=avl_real) YDUPL_B
+      REAL(kind=avl_real) XYZLES_B 
+      REAL(kind=avl_real) BSPACE
+      REAL(kind=avl_real) XBOD
+      REAL(kind=avl_real) YBOD
+      REAL(kind=avl_real) TBOD
+      REAL(kind=avl_real) XBOD_R
+      REAL(kind=avl_real) YBOD_R
+      COMMON /BODY_GEOM_R/
+     & XYZSCAL_B(3, NBMAX), ! scaling factors for XYZ coordinates
+     & XYZTRAN_B(3, NBMAX), ! translation factors for XYZ coordinates
+     & YDUPL_B(NBMAX), ! y duplicate
+     & XYZLES_B(3, NSMAX, NBMAX), ! leading edge cordinate vector
+     & BSPACE(NBMAX), ! body spacing
+     & XBOD(IBX, NBMAX), ! body oml x-coordinates
+     & YBOD(IBX, NBMAX), ! body oml y-coordinates
+     & TBOD(IBX, NBMAX), ! body oml thickness
+     & XBOD_R(IBX, NBMAX), ! raw input oml x-coordinates
+     & YBOD_R(IBX, NBMAX) ! raw input oml y-coordinates
+c
+      LOGICAL LDUPL
+      LOGICAL LSURFSPACING
+      COMMON /SURF_GEOM_L/
+     & LDUPL(NFMAX),      ! T if surface is a duplicated
+     & LSURFSPACING(NFMAX)   ! surface spacing set under the surface heeading
+
+      COMMON /SURF_GEOM_I/
+     & NSEC(NFMAX),  ! number of sections in surface
+     & NVC(NFMAX), ! number of chordwise elements
+     & NVS(NFMAX), ! number of spanwise elements
+     & NSPANS(NSECMAX, NFMAX), ! number of spanwise elements vector
+     & NASEC(NSECMAX, NFMAX), ! number of points used to represent section geometry
+     & NAPTSSEC(NSECMAX, IBX), ! number of points of input airfoil geometry
+     & ICONTD(ICONX, NSECMAX, NFMAX), ! control variable index
+     & NSCON(NSECMAX, NFMAX), ! number of control variables
+     & IDESTD(ICONX, NSECMAX, NFMAX), ! design variable index
+     & NSDES(NSECMAX, NFMAX) ! number of design variables
+
+      REAL(kind=avl_real) XYZSCAL
+      REAL(kind=avl_real) XYZTRAN
+      REAL(kind=avl_real) YDUPL
+      REAL(kind=avl_real) ADDINC
+      REAL(kind=avl_real) CSPACE
+      REAL(kind=avl_real) SSPACE
+      REAL(kind=avl_real) SSPACES
+      REAL(kind=avl_real) XYZLES
+      REAL(kind=avl_real) XSEC
+      REAL(kind=avl_real) YSEC
+      REAL(kind=avl_real) XFMIN_R
+      REAL(kind=avl_real) XFMAX_R
+c      REAL(kind=avl_real) XLASEC
+c      REAL(kind=avl_real) ZLASEC
+c      REAL(kind=avl_real) XUASEC
+c      REAL(kind=avl_real) ZUASEC
+      REAL(kind=avl_real) CHORDS
+      REAL(kind=avl_real) AINCS
+      REAL(kind=avl_real) XASEC
+      REAL(kind=avl_real) CASEC
+      REAL(kind=avl_real) SASEC
+      REAL(kind=avl_real) TASEC
+      REAL(kind=avl_real) CLCDSEC
+      REAL(kind=avl_real) CLCDSRF
+      REAL(kind=avl_real) CLAF
+      REAL(kind=avl_real) XHINGED
+      REAL(kind=avl_real) VHINGED
+      REAL(kind=avl_real) GAIND
+      REAL(kind=avl_real) REFLD
+      REAL(kind=avl_real) GAING  
+      COMMON /SURF_GEOM_R/
+     & XYZSCAL(3, NFMAX), ! scaling factors for XYZ coordinates
+     & XYZTRAN(3, NFMAX), ! translation factors for XYZ coordinates
+     & YDUPL(NFMAX), ! y position of duplicate plane
+     & ADDINC(NFMAX), ! additional incidence angle
+     & CSPACE(NFMAX), ! chordwise spacing
+     & SSPACE(NFMAX), ! spanwise spacing
+     & SSPACES(NSECMAX, NFMAX), ! spanwise spacing vector
+     & XYZLES(3, NSECMAX, NFMAX), ! leading edge cordinate vector
+     & XSEC(IBX, NSECMAX, NFMAX), ! raw section input x coordinates
+     & YSEC(IBX, NSECMAX, NFMAX), ! raw section input y coordinates
+     & XFMIN_R(NSECMAX, NFMAX), ! section input lower x bound 
+     & XFMAX_R(NSECMAX, NFMAX), ! section input upper x bound 
+c     & XLASEC(IBX, NSECMAX, NFMAX), ! airfoil lower surface x-coords
+c     & ZLASEC(IBX, NSECMAX, NFMAX), ! airfoil lower surface z-coords
+c     & XUASEC(IBX, NSECMAX, NFMAX), ! airfoil upper surface x-coords
+c     & ZUASEC(IBX, NSECMAX, NFMAX), ! airfoil upper surface z-coords
+     & CHORDS(NSECMAX, NFMAX), ! chord length vectorm
+     & AINCS(NSECMAX, NFMAX), ! incidence angle vector
+     & XASEC(50, NSECMAX, NFMAX),  ! the x coordinate aifoil section
+     & CASEC(50, NSECMAX, NFMAX),  ! camber line at xasec
+     & SASEC(50, NSECMAX, NFMAX),  ! slope of camber line at xasec
+     & TASEC(50, NSECMAX, NFMAX),   ! thickness at xasec
+     & CLCDSEC(6, NSECMAX, NFMAX),
+     & CLCDSRF(6, NFMAX),
+     & CLAF(NSECMAX, NFMAX),   ! CL alpha (dCL/da) scaling factor
+     & XHINGED(ICONX, NSECMAX, NFMAX), ! hinge location
+     & VHINGED(3, ICONX, NSECMAX, NFMAX), ! hinge vector
+     & GAIND(ICONX, NSECMAX, NFMAX), ! control surface gain
+     & REFLD(ICONX, NSECMAX, NFMAX), ! control surface reflection
+     & GAING(ICONX, NSECMAX, NFMAX)  ! desgin variable gain
+
+      COMMON /SURF_MESH_I/
+     & MFRST(NFMAX), ! stores the index in the MSHBLK where each surface's mesh begins
+     & IPTSEC(NSMAX,NFMAX) ! stores the iptloc vector for each surface
+
+      REAL(kind=avl_real) MSHBLK
+      REAL(kind=avl_real) RV1MSH
+      REAL(kind=avl_real) RV2MSH
+      REAL(kind=avl_real) RVMSH
+      REAL(kind=avl_real) RCMSH
+      COMMON /SURF_MESH_R/
+     & MSHBLK(3, 4*NVMAX), ! block to store all surface meshes
+     & RV1MSH(3,NVMAX),    !  mesh h.v. vortex left  points
+     & RV2MSH(3,NVMAX),    !  mesh h.v. vortex right  points
+     & RVMSH(3,NVMAX),    ! mesh h.v. vortex center points
+     & RCMSH(3,NVMAX)    ! mesh h.v. control points
+
+      LOGICAL LSURFMSH
+      LOGICAL LMESHFLAT
+      COMMON /SURF_MESH_L/
+     & LSURFMSH(NFMAX), ! T if surface uses a mesh cordinates for its geometry
+     & LMESHFLAT(NFMAX) ! T if the surface's mesh should be flattened when computing vorticies and control points
+
+
+C      !!--- end added variables for python geometry manipulation ---
diff --git a/src/second.f b/src/second.f
index a2daf9f..2ddaf21 100755
--- a/src/second.f
+++ b/src/second.f
@@ -1,12 +1,10 @@
-
       SUBROUTINE SECONDS(TSEC)
       REAL*8 TSEC
 C
-C...SECNDS is a real*4 function that returns seconds.
-C   The value is modified by subtracting the supplied argument.
-C   It acts as in the VMS FORTRAN Manual.
+C...Returns elapsed wall-clock time in seconds
+C   Replacement for non-standard SECNDS intrinsic
 C
-      REAL*4 SECNDS, TIME
-      TIME = 0.0
-      TSEC = SECNDS(TIME)
-      END
+      INTEGER COUNT, COUNT_RATE
+      CALL SYSTEM_CLOCK(COUNT, COUNT_RATE)
+      TSEC = DBLE(COUNT) / DBLE(COUNT_RATE)
+      END
\ No newline at end of file
diff --git a/tests/test_import.py b/tests/test_import.py
index 71c643d..16ea39e 100644
--- a/tests/test_import.py
+++ b/tests/test_import.py
@@ -1,7 +1,79 @@
 # =============================================================================
 # Extension modules
 # =============================================================================
+import struct, sys, sysconfig
+print("Pointer size:", struct.calcsize("P") * 8, "bit")  # 64 for both, less useful
+print("Platform tag:", sysconfig.get_platform())          # THIS is the key one
+print("sys.executable:", sys.executable)
 
+import ctypes
+import ctypes.wintypes
+import os, sys, struct, glob
+
+print("Python arch:", struct.calcsize("P") * 8, "bit")
+print("Executable:", sys.executable)
+
+# Check what VCRUNTIME140.dll Python itself loaded
+import ctypes.util
+kernel32 = ctypes.WinDLL("kernel32", use_last_error=True)
+GetModuleFileNameW = kernel32.GetModuleFileNameW
+GetModuleFileNameW.argtypes = [ctypes.wintypes.HMODULE, ctypes.wintypes.LPWSTR, ctypes.wintypes.DWORD]
+GetModuleFileNameW.restype = ctypes.wintypes.DWORD
+GetModuleHandleW = kernel32.GetModuleHandleW
+GetModuleHandleW.argtypes = [ctypes.wintypes.LPCWSTR]
+GetModuleHandleW.restype = ctypes.wintypes.HMODULE
+
+for dll_name in ["VCRUNTIME140.dll", "msvcp140.dll", "python312.dll"]:
+    h = GetModuleHandleW(dll_name)
+    if h:
+        buf = ctypes.create_unicode_buffer(260)
+        GetModuleFileNameW(h, buf, 260)
+        print(f"{dll_name}: loaded from {buf.value}")
+    else:
+        print(f"{dll_name}: not loaded")
+
+# Search common locations for VCRUNTIME140.dll
+search_dirs = [
+    os.path.dirname(sys.executable),
+    os.path.join(os.environ.get("SystemRoot", r"C:\Windows"), "System32"),
+    os.path.join(os.environ.get("SystemRoot", r"C:\Windows"), "SysWOW64"),
+]
+for d in search_dirs:
+    path = os.path.join(d, "VCRUNTIME140.dll")
+    if os.path.exists(path):
+        # Read PE machine type
+        with open(path, "rb") as f:
+            f.seek(0x3C)  # e_lfanew offset
+            pe_offset = int.from_bytes(f.read(4), "little")
+            f.seek(pe_offset + 4)  # Machine field
+            machine = int.from_bytes(f.read(2), "little")
+        arch = {0x8664: "x86_64", 0xAA64: "ARM64", 0x14C: "x86"}.get(machine, hex(machine))
+        print(f"Found {path} -> {arch}")
+
+# Try a direct load of the .pyd to get better error info
+print("\n--- Direct load test ---")
+import importlib.util
+spec = importlib.util.find_spec("optvl")
+if spec:
+    pyd = glob.glob(os.path.join(spec.submodule_search_locations[0], "libavl*.pyd"))
+    if pyd:
+        print(f"Attempting ctypes load of: {pyd[0]}")
+        try:
+            ctypes.WinDLL(pyd[0])
+            print("SUCCESS: Direct load works")
+        except OSError as e:
+            print(f"FAILED: {e}")
+        
+        # Now test loading from a temp copy (what MExt does)
+        import shutil, tempfile
+        tmpdir = tempfile.mkdtemp()
+        copied = shutil.copy(pyd[0], tmpdir)
+        print(f"Attempting ctypes load of temp copy: {copied}")
+        try:
+            ctypes.WinDLL(copied)
+            print("SUCCESS: Temp copy load works")
+        except OSError as e:
+            print(f"FAILED from temp: {e}")
 # =============================================================================
 # Standard Python Modules
 # =============================================================================
diff --git a/tools/wheels/cibw_test_command.sh b/tools/wheels/cibw_test_command.sh
index 99ebf95..fbb792f 100644
--- a/tools/wheels/cibw_test_command.sh
+++ b/tools/wheels/cibw_test_command.sh
@@ -3,9 +3,13 @@ set -xe
 PROJECT_DIR="$1"
 
 cd $PROJECT_DIR/tests
+python -c "import platform; print(platform.machine())"
+ls ../optvl
+# llvm-objdump -p ../optvl/libavl.cp312-win_arm64.pyd | findstr "DLL Name"
 
 # install tesing dependencies
-pip install "scipy<=1.16" psutil "openmdao!=3.38"
+pip install --only-binary :all: scipy 
+pip install psutil "openmdao!=3.38"
 
 
 #HACK: if the tests are not split up the CI runs out of memory...
diff --git a/tools/wheels/repair_windows.sh b/tools/wheels/repair_windows.sh
index 583e589..58a638d 100644
--- a/tools/wheels/repair_windows.sh
+++ b/tools/wheels/repair_windows.sh
@@ -1,32 +1,45 @@
+# from scipy
+
 set -xe
 
 WHEEL="$1"
 DEST_DIR="$2"
 
-# create a temporary directory in the destination folder and unpack the wheel
-# into there
-pushd $DEST_DIR
-mkdir -p tmp
-pushd tmp
-wheel unpack $WHEEL
-pushd optvl*
 
-# To avoid DLL hell, the file name of libopenblas that's being vendored with
-# the wheel has to be name-mangled. delvewheel is unable to name-mangle PYD
-# containing extra data at the end of the binary, which frequently occurs when
-# building with mingw.
-# We therefore find each PYD in the directory structure and strip them.
+# Skip the strip command based on TARGET_ARCH
+# TARGET_ARCH should be set by the CI environment (e.g., ARM64, AMD64)
+TARGET_ARCH="${TARGET_ARCH:-}" # Default to empty string if not set
+
+if [ "$TARGET_ARCH" = "ARM64" ]; then
+    echo "Skipping stripping for ARM64 target."
+
+else
+    echo "Performing stripping for AMD64 target."
+
+    # create a temporary directory in the destination folder and unpack the wheel
+    # into there
+    pushd $DEST_DIR
+    mkdir -p tmp
+    pushd tmp
+    wheel unpack $WHEEL
+    pushd optvl*
+
+    # To avoid DLL hell, the file name of libopenblas that's being vendored with
+    # the wheel has to be name-mangled. delvewheel is unable to name-mangle PYD
+    # containing extra data at the end of the binary, which frequently occurs when
+    # building with mingw.
+    # We therefore find each PYD in the directory structure and strip them.
+
+    for f in $(find ./optvl* -name '*.pyd'); do strip $f; done
 
-for f in $(find ./optvl* -name '*.pyd'); do strip $f; done
 
+    # now repack the wheel and overwrite the original
+    wheel pack .
+    mv -fv *.whl $WHEEL
 
-# now repack the wheel and overwrite the original
-wheel pack .
-mv -fv *.whl $WHEEL
+    cd $DEST_DIR
+    rm -rf tmp
+fi
 
-cd $DEST_DIR
-rm -rf tmp
+delvewheel repair  -w $DEST_DIR $WHEEL
 
-# the libopenblas.dll is placed into this directory in the cibw_before_build
-# script.
-delvewheel repair --add-path /c/opt/openblas/openblas_dll -w $DEST_DIR $WHEEL