diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4dc0e47..53a51e0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -29,7 +29,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Removed
 - Removed C-I tests on Microsoft Azure Dev Pipelines
-- Replaced BLAS functions (`WAXPY`, `WCOPY`, `WSCAL`, `WADD`, `WLAMCH`) with pure F90 code from `int/*.f90` integrators (thanks to AI for the help)
+- Replaced BLAS functions (`WAXPY`, `WCOPY`, `WSCAL`, `WADD`, `WLAMCH`, `WDOT`) with pure F90 code from `int/*.f90` integrators (thanks to AI for the help)
 
 ## [3.3.0] - 2025-07-17
 ### Added
diff --git a/util/blas.f90 b/util/blas.f90
index 06c4fa5..c1a9a9b 100644
--- a/util/blas.f90
+++ b/util/blas.f90
@@ -20,70 +20,11 @@
 !%%% (5) WLAMCH_ADD                                        %%%
 !%%% (6) SET2ZERO                                          %%%
 !%%% (7) WADD                                              %%%
+!%%% (8) WDOT                                              %%%
 !%%%                                                       %%%
 !%%% @yantosca, 17 Oct 2025                                %%%
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-!--------------------------------------------------------------
-      KPP_REAL FUNCTION WDOT (N, DX, incX, DY, incY) 
-!--------------------------------------------------------------
-!     dot produce: wdot = x(1:N)*y(1:N) 
-!     only for incX=incY=1
-!     after BLAS
-!     replace this by the function from the optimized BLAS implementation:
-!         CALL SDOT(N,X,1,Y,1) or  CALL DDOT(N,X,1,Y,1)
-!--------------------------------------------------------------
-!      USE messy_mecca_kpp_Precision
-!--------------------------------------------------------------
-      IMPLICIT NONE
-      INTEGER :: N, incX, incY
-      KPP_REAL :: DX(N), DY(N) 
-
-      INTEGER :: i, IX, IY, M, MP1, NS
-                                 
-      WDOT = 0.0D0 
-      IF (N .LE. 0) RETURN 
-      IF (incX .EQ. incY) IF (incX-1) 5,20,60 
-!                                                                       
-!     Code for unequal or nonpositive increments.                       
-!                                                                       
-    5 IX = 1 
-      IY = 1 
-      IF (incX .LT. 0) IX = (-N+1)*incX + 1 
-      IF (incY .LT. 0) IY = (-N+1)*incY + 1 
-      DO i = 1,N 
-        WDOT = WDOT + DX(IX)*DY(IY) 
-        IX = IX + incX 
-        IY = IY + incY 
-      END DO 
-      RETURN 
-!                                                                       
-!     Code for both increments equal to 1.                              
-!                                                                       
-!     Clean-up loop so remaining vector length is a multiple of 5.      
-!                                                                       
-   20 M = MOD(N,5) 
-      IF (M .EQ. 0) GO TO 40 
-      DO i = 1,M 
-         WDOT = WDOT + DX(i)*DY(i) 
-      END DO 
-      IF (N .LT. 5) RETURN 
-   40 MP1 = M + 1 
-      DO i = MP1,N,5 
-          WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) +  &
-                   DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)                   
-      END DO 
-      RETURN 
-!                                                                       
-!     Code for equal, positive, non-unit increments.                    
-!                                                                       
-   60 NS = N*incX 
-      DO i = 1,NS,incX 
-        WDOT = WDOT + DX(i)*DY(i) 
-      END DO 
-
-      END FUNCTION WDOT                                          
-
 !--------------------------------------------------------------
       SUBROUTINE WGEFA(N,A,Ipvt,info)
 !--------------------------------------------------------------
@@ -203,14 +144,14 @@ SUBROUTINE WGESL(Trans,N,A,Ipvt,b)
 
 !        first solve  trans(U)*y = b
          DO k = 1, n
-            t = WDOT(k-1,a(1,k),1,b(1),1)
+            t = DOT_PRODUCT( a(1:k-1, k), b(1:k-1) )
             b(k) = (b(k) - t)/a(k,k)
          END DO
 !        now solve trans(L)*x = y
          IF (n >= 2) THEN
          DO kb = 1, n-1
             k = n - kb
-            b(k) = b(k) + WDOT(n-k,a(k+1,k),1,b(k+1),1)
+            b(k) = b(k) + DOT_PRODUCT( a(k+1:n, k), b(k+1:n) )
             l = Ipvt(k)
             IF (l /= k) THEN
                t = b(l); b(l) = b(k); b(k) = t