Remove BLAS functions from int/sdirk*.f90 integrators

int/sdirk.f90
int/sdirk4.f90
int/sdirk_adj.f90
- Replaced calls to BLAS functions with pure F90 code
  (as determined by AI)

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

Author: yantosca

int/sdirk.f90

@@ -22,9 +22,7 @@ MODULE KPP_ROOT_Integrator
   USE KPP_ROOT_Global
   USE KPP_ROOT_Parameters
   USE KPP_ROOT_JacobianSP,    ONLY : LU_DIAG
-  USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, &
-                                     WLAMCH,    WCOPY,    WAXPY,    &
-                                     WSCAL,     WADD
+  USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, WLAMCH
 
   IMPLICIT NONE
   PUBLIC
@@ -591,17 +589,18 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr )
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 !~~~>  Starting values for Newton iterations
-       CALL Set2zero(N,Z(1,istage))
-
+!       CALL Set2zero(N,Z(1,istage))
+       G(1:N)        = 0.0_dp
+       Z(1:N,istage) = 0.0_dp
 !~~~>   Prepare the loop-independent part of the right-hand side
-       CALL Set2zero(N,G)
+!      CALL Set2zero(N,G)
        IF (istage > 1) THEN
            DO j = 1, istage-1
                ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj)
-               CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1)
+               G(1:N) = G(1:N) + rkTheta(istage,j) * Z(1:N,j)
                ! Zi(:) = sum_j Alpha(i,j)*Zj(:)
                IF (StartNewton) THEN
-                  CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1)
+                  Z(1:N,istage) = Z(1:N,istage) + rkAlpha(istage,j) * Z(1:N,j)
                END IF
            END DO
        END IF
@@ -613,13 +612,13 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr )
 NewtonLoop:DO NewtonIter = 1, NewtonMaxit
 
 !~~~>   Prepare the loop-dependent part of the right-hand side
-            CALL WADD(N,Y,Z(1,istage),TMP)              ! TMP <- Y + Zi
+            TMP(1:N) = Y(1:N) + Z(1:N,istage)           ! TMP <- Y + Zi
             CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS)      ! RHS <- Fun(Y+Zi)
             ISTATUS(Nfun) = ISTATUS(Nfun) + 1
 !            RHS(1:N) = G(1:N) - Z(1:N,istage) + (H*rkGamma)*RHS(1:N)
-            CALL WSCAL(N, H*rkGamma, RHS, 1)
-            CALL WAXPY (N, -ONE, Z(1,istage), 1, RHS, 1)
-            CALL WAXPY (N, ONE, G,1, RHS,1)
+            RHS(1:N) = RHS(1:N) * (H * rkGamma)
+            RHS(1:N) = RHS(1:N) - Z(1:N,istage)
+            RHS(1:N) = RHS(1:N) + G(1:N)
 
 !~~~>   Solve the linear system
             CALL SDIRK_Solve ( H, N, E, IP, IER, RHS )
@@ -648,7 +647,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr )
             END IF
             NewtonIncrementOld = NewtonIncrement
             ! Update solution: Z(:) <-- Z(:)+RHS(:)
-            CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1)
+            Z(1:N,istage) = Z(1:N,istage) + RHS(1:N)
 
             ! Check error in Newton iterations
             NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
@@ -677,7 +676,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr )
       IF (sdMethod /= BEL) THEN ! All methods but Backward Euler
         CALL Set2zero(N,TMP)
         DO i = 1,rkS
-          IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,TMP,1)
+          IF (rkE(i)/=ZERO) TMP(1:N) = TMP(1:N) + rkE(i) * Z(1:N,i)
         END DO
 
         CALL SDIRK_Solve( H, N, E, IP, IER, TMP )
@@ -704,7 +703,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr )
          T  =  T + H
          ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:)
          DO i = 1,rkS
-            IF (rkD(i)/=ZERO) CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1)
+            IF (rkD(i)/=ZERO) Y(1:N) = Y(1:N) + rkD(i) * Z(1:N,i)
          END DO
 
 !~~~> Update scaling coefficients
@@ -918,10 +917,10 @@ SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS )
       KPP_REAL, INTENT(IN) :: E(LU_NONZERO)
 #endif
       KPP_REAL, INTENT(INOUT) :: RHS(N)
-      KPP_REAL                :: HGammaInv
 
-      HGammaInv = ONE/(H*rkGamma)
-      CALL WSCAL(N,HGammaInv,RHS,1)
+      ! NOTE: This line reproduces the results of the
+      ! previous WAXPY call (@yantosca, 16 Oct 2025)
+      RHS(1:N) = RHS(1:N) * (ONE / (H * rkGamma))
 #ifdef FULL_ALGEBRA
       CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING )
 #else

int/sdirk4.f90

@@ -12,8 +12,7 @@ MODULE KPP_ROOT_Integrator
   USE KPP_ROOT_Global
   USE KPP_ROOT_Parameters
   USE KPP_ROOT_JacobianSP,    ONLY : LU_DIAG
-  USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, &
-                                     WLAMCH,    WAXPY,    WCOPY
+  USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, WLAMCH
 
   IMPLICIT NONE
   PUBLIC
@@ -553,11 +552,11 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, &
       NewtonFactor(istage) = MAX(NewtonFactor(istage),Roundoff)**0.8d0
 
 !~~~>  STARTING VALUES FOR NEWTON ITERATION
-      CALL Set2zero(N,G)
-      CALL Set2zero(N,Z(1,istage))
+      G(1:N)        = 0.0_dp
+      Z(1:N,istage) = 0.0_dp
       IF (istage==1) THEN
           IF (FIRST.OR.NewtonReject) THEN
-              CALL Set2zero(N,Z(1,istage))
+              Z(1:N,istage) = 0.0_dp
           ELSE
               W=ONE+rkGamma*H/Hold
               DO i=1,N
@@ -567,12 +566,11 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, &
       ELSE
           DO j = 1, istage-1
             ! Gj(:) = sum_j Beta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:))
-            CALL WAXPY(N,rkBeta(istage,j),Z(1,j),1,G,1)
-            ! CALL WAXPY(N,H*rkA(istage,j),FV(1,j),1,G,1)
+            G(1:N) = G(1:N) + rkBeta(istage,j) * Z(1:N,j)
             ! Zi(:) = sum_j Alpha(i,j)*Zj(:)
-            CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1)
+            Z(1:N,istage) = Z(1:N,istage) + rkAlpha(istage,j) * Z(1:N,j)
           END DO
-          IF (istage==5) CALL WCOPY(N,Z(1,istage),1,Yhat,1)  ! Yhat(:) <- Z5(:)
+          IF (istage==5) Yhat(1:N) = Z(1:N,istage)  ! Yhat(:) <- Z5(:)
       END IF
 
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -597,7 +595,8 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, &
             TMP(1:N) = Y(1:N) + Z(1:N,istage)
             CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS)
             TMP(1:N) = G(1:N) - Z(1:N,istage)
-            CALL WAXPY(N,HGammaInv,TMP,1,RHS,1) ! RHS(:) <- RHS(:) + HGammaInv*(G(:)-Z(:))
+            ! RHS(:) <- RHS(:) + HGammaInv*(G(:)-Z(:))
+            RHS(1:N) = RHS(1:N) + HGammaInv * TMP(1:N)
 
 !~~~>     SOLVE THE LINEAR SYSTEMS
 #ifdef FULL_ALGEBRA  
@@ -630,8 +629,9 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, &
                 END IF
             END IF
             NewtonErrOld = NewtonErr
-            CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) ! Z(:) <-- Z(:)+RHS(:)
-            
+            ! Z(:) <-- Z(:)+RHS(:)
+            Z(1:N,istage) = Z(1:N,istage) + RHS(1:N)            
+
             END DO Newton
 
 !~~> END OF SIMPLIFIED NEWTON ITERATION
@@ -671,13 +671,13 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, &
          Hold=H
 
 !~~~> COEFFICIENTS FOR CONTINUOUS SOLUTION
-         CALL WAXPY(N,ONE,Z(1,5),1,Y,1) ! Y(:) <-- Y(:)+Z5(:)
-         CALL WCOPY(N,Z(1,5),1,Yhat,1)  ! Yhat <-- Z5
+         Y(1:N)    = Y(1:N) + Z(1:N,5)  ! Y(:) <-- Y(:)+Z5(:)
+         Yhat(1:N) = Z(1:N,5)           ! Yhat <-- Z5
 
          DO i=1,4  ! CONTi <-- Sum_j rkD(i,j)*Zj
-           CALL Set2zero(N,CONT(1,i))
+           CONT(1:N,i) = 0.0_dp
            DO j = 1,5
-             CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1)
+              CONT(1:N,i) = CONT(1:N,i) + rkD(i,j) * Z(1:N,j)
            END DO
          END DO
 
@@ -997,11 +997,7 @@ SUBROUTINE JAC_CHEM( T, Y, JV )
 
 #ifdef FULL_ALGEBRA
       CALL Jac_SP(Y, FIX, RCONST, JS)
-      DO j=1,NVAR
-         DO j=1,NVAR
-          JV(i,j) = 0.0D0
-         END DO
-      END DO
+      JV(1:NVAR,1:NVAR) = 0.0d0
       DO i=1,LU_NONZERO
          JV(LU_IROW(i),LU_ICOL(i)) = JS(i)
       END DO