Index: src/Mat_util.f90
===================================================================
--- src/Mat_util.f90	(revision 1)
+++ src/Mat_util.f90	(revision 2)
@@ -80,7 +80,6 @@
       IMPLICIT NONE
 
       INTEGER(KINT), INTENT(IN) :: N,NREAL
-      INTEGER(KINT), PARAMETER :: Max_it=500
       REAL(KREAL) :: AP(N*(N+1)/2)
       REAL(KREAL), ALLOCATABLE ::  A(:,:)
       REAL(KREAL) :: V(Nreal,Nreal),D(Nreal)
@@ -265,7 +264,7 @@
   Use VarTypes
   IMPLICIT NONE
 
-  integer(KINT) :: i,j,k,nb_niv,max_niv, nreal
+  integer(KINT) :: i, j, k, nb_niv, nreal
   real(KREAL) :: ene(nreal),psi(nreal,nreal)
   real(KREAL) :: a
 
@@ -320,7 +319,7 @@
   INTEGER(KINT), INTENT(IN) :: N
   REAL(KREAL), INTENT(OUT) :: EigVal(N), Eigvec(N,N)
 
-  integer(KINT) :: i,j,k
+  integer(KINT) :: i
 
   INTEGER(KINT), ALLOCATABLE :: Rank(:) !N
   REAL(KREAL), ALLOCATABLE :: EigValT(:) !n
Index: src/upcase.f90
===================================================================
--- src/upcase.f90	(revision 1)
+++ src/upcase.f90	(revision 2)
@@ -3,8 +3,9 @@
 ! This subroutine converts a string into its Upcase
 ! version
 
+        IMPLICIT NONE
+
         integer, parameter :: KINT = kind(1)
-        integer, parameter :: KREAL = kind(1.0d0)
 
         CHARACTER(*), INTENT(INOUT) :: String
         INTEGER(KINT) :: i,ic, length
Index: src/PathCreate.f90
===================================================================
--- src/PathCreate.f90	(revision 1)
+++ src/PathCreate.f90	(revision 2)
@@ -24,7 +24,6 @@
   REAL(KREAL), PARAMETER :: Inf=1e32
   INTEGER(KINT), INTENT(IN) :: IOpt
 
-  INTEGER(KINT) ::  at,long
   REAL(KREAL), ALLOCATABLE :: val_zmat(:,:), val_zmatTmp(:,:) ! (3,Nat)
   REAL(KREAL), ALLOCATABLE :: X0(:), Y0(:), Z0(:) ! Nat
   REAL(KREAL), ALLOCATABLE ::  XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)
@@ -35,14 +34,13 @@
   ! NPrim=Number of primitive coordinates.
   REAL(KREAL), ALLOCATABLE :: XGeom(:) ! NGeomI
   REAL(KREAL), ALLOCATABLE :: XGeomF(:) ! NGeomF
-  REAL(KREAL), ALLOCATABLE :: Hess_int(:,:) ! Hess_int (3*Nat,3*Nat), used for Baker case
   INTEGER(KINT) :: NGeomS
 
-  REAL(KREAL) :: Rmsd,MRot(3,3), Delta, normDeltaX_int
-  REAL(KREAL), PARAMETER ::  LimAngle=145.d0,TwoPi=360.d0
-  REAL(KREAL)  ::s,ds,dist,a_val,d
-  INTEGER(KINT) :: I,J,K,NbCycle,igeom,nb,ig,iat
-  INTEGER(KINT) :: Itmp, Jat, L
+  REAL(KREAL) :: Rmsd, MRot(3, 3), Delta
+  REAL(KREAL), PARAMETER ::  TwoPi=360.d0
+  REAL(KREAL)  ::s, dist, a_val, d
+  INTEGER(KINT) :: I, J, igeom, iat
+  INTEGER(KINT) :: Itmp, Jat
   INTEGER(KINT) :: Idx,IBeg
   CHARACTER(LCHARS) :: Title
   REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)
Index: src/egrad_test.f90
===================================================================
--- src/egrad_test.f90	(revision 1)
+++ src/egrad_test.f90	(revision 2)
@@ -12,22 +12,20 @@
        REAL(KREAL), INTENT(IN) :: Geom(Nat,3)
 
 !     Bohr --> Angstr
-      real(KREAL), parameter :: BOHR   = 0.52917726D+00
+!      real(KREAL), parameter :: BOHR   = 0.52917726D+00
 !
+
 ! Parameters to define the surface
       INTEGER(KINT), DIMENSION(6), PARAMETER :: IECOEF = (/-1,9,-45,45,-9,1/)
       INTEGER(KINT), DIMENSION(6), PARAMETER :: ISCOEF = (/-3,-2,-1,1,2,3/)
       REAL(KREAL), PARAMETER :: hh=0.001d0
       
 ! Variables
-       INTEGER(KINT) :: i,j,iat,jat
+       INTEGER(KINT) :: i, iat, jat
        REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:)
-       REAL(KREAL) :: xp,yp
 
-       CHARACTER(132) :: Line
        LOGICAL :: Debug
   
-      real(KREAL), parameter :: zero = 0.0_KREAL,   one = 1.0_KREAL
       REAL(KREAL), external :: Exyz 
 
 
@@ -107,16 +105,14 @@
       INTEGER(KINT) :: i
       REAL(KREAL) :: dCH,dNH,dCN
 
-      real(KREAL) :: autoA,autoeV
-      parameter (autoA=0.52917715d0)
+!      real(KREAL) :: autoA,autoeV
+      real(KREAL) :: autoeV
+!      parameter (autoA=0.52917715d0)
       parameter (autoeV=27.21183d0)
 
 
-      real(KREAL) :: Q(3)
-      real(KREAL) :: TT,gR,pR
 
-      real(KREAL) :: mB,mC,mH
-      real(KREAL) :: T,RN,RC,r1,r2,r3
+      real(KREAL) :: r1, r2, r3
 
       real(KREAL) :: Z1,Z12,V1
       real(KREAL) :: Z2,Z22,V2
@@ -228,26 +224,17 @@
       REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3)
       REAL(KREAL) :: Exyz
 
-     INTEGER(KINT) :: i,j,iat,jat
+     INTEGER(KINT) :: iat, jat
      REAL(KREAL) :: Exyz2
 
-      real(KREAL) :: autoA,autoeV
-      parameter (autoA=0.52917715d0)
-      parameter (autoeV=27.21183d0)
+!      real(KREAL) :: autoA,autoeV
+!      parameter (autoA=0.52917715d0)
+!      parameter (autoeV=27.21183d0)
 
 
-      real(KREAL) :: Q(3)
-      real(KREAL) :: TT,gR,pR
 
-      real(KREAL) :: mB,mC,mH
-      real(KREAL) :: T,RN,RC,r1,r2,r3
 
-      real(KREAL) :: Z1,Z12,V1
-      real(KREAL) :: Z2,Z22,V2
-      real(KREAL) :: Z3,Z32,V3
 
-      real(KREAL) :: S1,S2,S3,S12,S22,S32,POLY
-      real(KREAL) :: E1,E2,E3,HY1,HY2,HY3,SWITCH,V123
 
       Exyz=0.
       do jat=1,natoms
Index: src/Write_vasp.f90
===================================================================
--- src/Write_vasp.f90	(revision 1)
+++ src/Write_vasp.f90	(revision 2)
@@ -7,16 +7,14 @@
 
   CHARACTER(SCHARS), INTENT(IN) :: poscar
 
-  INTEGER(KINT) :: IGeom, Idx, Iat, I,J,NFr
+  INTEGER(KINT) :: IGeom, Idx
   REAL(KREAL), ALLOCATABLE :: GeomTmpC(:), GeomTmp(:)
-  REAL(KREAL) :: MRot(3,3), Rmsd
   CHARACTER(SCHARS) :: Line,Line2
 
   ! For VMD printing
   CHARACTER(SCHARS), PARAMETER :: FirstWord="mol "
   CHARACTER(LCHARS), PARAMETER :: LastWords=" type POSCAR first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all"
 
-  REAL(KREAL) :: X,Y,Z
   CHARACTER(LCHARS) :: TITLE
   LOGICAL :: Debug
 
Index: src/Calc_zmat_constr_frag.f90
===================================================================
--- src/Calc_zmat_constr_frag.f90	(revision 1)
+++ src/Calc_zmat_constr_frag.f90	(revision 2)
@@ -45,7 +45,7 @@
   real(KREAL) ::  vx4,vy4,vz4,norm4
   real(KREAL) ::  vx5,vy5,vz5,norm5
   real(KREAL) ::  val,val_d, d12, d13,d23,d
-  Logical PasFini,Debug
+  Logical Debug
   LOGICAL, ALLOCATABLE :: DejaFait(:), FCaf(:) !(na)
   Logical, ALLOCATABLE ::  FrozAt(:)    !T if this atom is frozen
   LOGICAL F1213, F1223,F1323
@@ -492,7 +492,7 @@
         WRITE(*,*) "Dans le caca car tous les atoms de ce block sont alignés!"
      END IF
 
-     CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
+     CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
           vx5,vy5,vz5,norm5)
 
 
@@ -534,7 +534,7 @@
         ! les autres fragments ne seraient additonnés que si l'on ne trouve pas notre bonheur dans le premier bloc
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-        CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
+        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
              vx4,vy4,vz4,norm4)
         val_d=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
              vx2,vy2,vz2,norm2)
@@ -558,7 +558,7 @@
            if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
               n1=JAt
               CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-              CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
+              CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
                    vx4,vy4,vz4,norm4)
               val_d=angle_d(vx4,vy4,vz4,norm4, &
                    vx5,vy5,vz5,norm5, &
@@ -578,8 +578,7 @@
               if (.NOT.DejaFait(JAt)) THEN
                  n1=JAt
                  CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-                 CALL produit_vect(vx1,vy1,vz1,norm1, &
-                      vx2,vy2,vz2,norm2, &
+                 CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
                       vx4,vy4,vz4,norm4)
                  val_d=angle_d(vx4,vy4,vz4,norm4, &
                       vx5,vy5,vz5,norm5, &
@@ -697,16 +696,14 @@
      n4=LiaisonsBis(Iat,2)
      CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
      CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
-     CALL produit_vect(vx3,vy3,vz3,norm3, &
-          vx2,vy2,vz2,norm2, &
+     CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
           vx5,vy5,vz5,norm5)
 
      DO While ((ITmp.LE.FrozBlock(Iat,0)).AND.(sDihe.LE.0.09d0))
         ITmp=ITmp+1
         n1=FrozBlock(Iat,Itmp)
         CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-        CALL produit_vect(vx1,vy1,vz1,norm1, &
-             vx2,vy2,vz2,norm2, &
+        CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, &
              vx4,vy4,vz4,norm4)
         val_d=angle_d(vx4,vy4,vz4,norm4,  &
              vx5,vy5,vz5,norm5,           &
@@ -725,8 +722,7 @@
            if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
               n1=JAt
               CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-              CALL produit_vect(vx1,vy1,vz1,norm1,  &
-                   vx2,vy2,vz2,norm2,               &
+              CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2,   &
                    vx4,vy4,vz4,norm4)
               val_d=angle_d(vx4,vy4,vz4,norm4,    &
                    vx5,vy5,vz5,norm5,              &
@@ -746,8 +742,7 @@
               if (.NOT.DejaFait(JAt)) THEN
                  n1=JAt
                  CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-                 CALL produit_vect(vx1,vy1,vz1,norm1,   &
-                      vx2,vy2,vz2,norm2,                 &
+                 CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2,   &
                       vx4,vy4,vz4,norm4)
                  val_d=angle_d(vx4,vy4,vz4,norm4,       &
                       vx5,vy5,vz5,norm5,                 &
@@ -889,8 +884,7 @@
      n4=I1
      CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
      CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
-     CALL produit_vect(vx3,vy3,vz3,norm3,  &
-          vx2,vy2,vz2,norm2,    &
+     CALL produit_vect(vx3,vy3,vz3, vx2,vy2,vz2,  &
           vx5,vy5,vz5,norm5)
 
      !     We will look at the other frozen atoms
@@ -901,8 +895,7 @@
         if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
            n1=JAt
            CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-           CALL produit_vect(vx1,vy1,vz1,norm1, &
-                vx2,vy2,vz2,norm2, &
+           CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, &
                 vx4,vy4,vz4,norm4)
            val_d=angle_d(vx4,vy4,vz4,norm4,  &
                 vx5,vy5,vz5,norm5,   &
@@ -922,8 +915,7 @@
            if (.NOT.DejaFait(JAt)) THEN
               n1=JAt
               CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-              CALL produit_vect(vx1,vy1,vz1,norm1,  &
-                   vx2,vy2,vz2,norm2, &
+              CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, &
                    vx4,vy4,vz4,norm4)
               val_d=angle_d(vx4,vy4,vz4,norm4,   &
                    vx5,vy5,vz5,norm5,           &
@@ -1965,9 +1957,9 @@
         val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
 
         CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
-        CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
+        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
              vx4,vy4,vz4,norm4)
-        CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
+        CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
              vx5,vy5,vz5,norm5)
 
         val_d=angle_d(vx4,vy4,vz4,norm4, vx5,vy5,vz5,norm5, &
Index: src/Calc_mixed_frag.f90
===================================================================
--- src/Calc_mixed_frag.f90	(revision 1)
+++ src/Calc_mixed_frag.f90	(revision 2)
@@ -46,18 +46,18 @@
   real(KREAL) ::  vx3,vy3,vz3,norm3
   real(KREAL) ::  vx4,vy4,vz4,norm4
   real(KREAL) ::  vx5,vy5,vz5,norm5
-  real(KREAL) ::  val,val_d, d12, d13,d23,d
-  Logical PasFini,Debug
+  real(KREAL) ::  val, val_d, d
+  Logical Debug
   LOGICAL, ALLOCATABLE :: DejaFait(:),FCaf(:) !(na)
   Logical, ALLOCATABLE ::  FrozAt(:)    !T if this atom is frozen
-  LOGICAL F1213, F1223,F1323
 
 
-  INTEGER(KINT) :: I,J,n0,n1,n2,n3,n4,IAt,IL,JL,IFrag,ITmp, K, KMax
-  INTEGER(KINT) :: IMax, I3,I1,Ip, IFragFroz,IBeg
+
+  INTEGER(KINT) :: I, J, n1, n2, n3, n4, IAt, IL, JL, IFrag, ITmp
+  INTEGER(KINT) :: I3, I1, Ip, IFragFroz, IBeg
   INTEGER(KINT) :: I0, Izm, JAt,Izm0,Idx
 
-  REAL(KREAL) :: sDihe, Pi, Sang
+  REAL(KREAL) :: Pi, Sang
 
   INTERFACE
      function valid(string) result (isValid)
@@ -938,9 +938,9 @@
      val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
 
      CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
-     CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
+     CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
           vx4,vy4,vz4,norm4)
-     CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
+     CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
           vx5,vy5,vz5,norm5)
 
      val_d=angle_d(vx4,vy4,vz4,norm4, vx5,vy5,vz5,norm5, &
Index: src/TestCart.f90
===================================================================
--- src/TestCart.f90	(revision 1)
+++ src/TestCart.f90	(revision 2)
@@ -22,8 +22,7 @@
 
   LOGICAL, INTENT(IN) :: SetCart
 
-  INTEGER(KINT) :: IGeom, I,J,K,Iat
-  CHARACTER(SCHARS) :: Line,Line2
+  INTEGER(KINT) :: I, J, K, Iat
 
   REAL(KREAL) :: d
   LOGICAL, ALLOCATABLE :: TobeCart(:) ! Nat
Index: src/gaussj.f90
===================================================================
--- src/gaussj.f90	(revision 1)
+++ src/gaussj.f90	(revision 2)
@@ -27,7 +27,8 @@
                   icol=k
                 endif
               else if (ipiv(k).gt.1) then
-                pause 'singular matrix in gaussj'
+                Write(*,*) 'singular matrix in gaussj'
+                STOP
               endif
 12          continue
           endif
@@ -47,7 +48,10 @@
         endif
         indxr(i)=irow
         indxc(i)=icol
-        if (a(icol,icol).eq.0.d0) pause 'singular matrix in gaussj'
+        if (a(icol,icol).eq.0.d0) THEN
+           WRITE(*,*)  'singular matrix in gaussj'
+           STOP
+        END IF
         pivinv=1.d0/a(icol,icol)
         a(icol,icol)=1.d0
         do 16 l=1,n
Index: src/CalcBMat_mixed.f90
===================================================================
--- src/CalcBMat_mixed.f90	(revision 1)
+++ src/CalcBMat_mixed.f90	(revision 2)
@@ -1,4 +1,4 @@
-subroutine CalcBMat_mixed (natoms, xyzatm, indzmat, dzdc, atmass,atome)
+subroutine CalcBMat_mixed (natoms, xyzatm, indzmat, dzdc)
 
   use Io_module
   use Path_module, only : MW, NCart
@@ -6,8 +6,9 @@
   IMPLICIT NONE
 
 
-  integer(KINT) :: natoms, indzmat(natoms,5),atome(natoms)
-  real(KREAL)   :: xyzatm(3,natoms), dzdc(3*natoms,3*natoms),atmass(natoms)
+  integer(KINT), INTENT(IN) :: natoms, indzmat(natoms,5)
+  real(KREAL), INTENT(IN)   :: xyzatm(3,natoms)
+  real(KREAL), INTENT(OUT)   :: dzdc(3*natoms,3*natoms)
 
 
 !!!!!!!!!!!!!!!!!!!!!
Index: src/Calc_zmat.f90
===================================================================
--- src/Calc_zmat.f90	(revision 1)
+++ src/Calc_zmat.f90	(revision 2)
@@ -20,26 +20,19 @@
       integer(KINT), INTENT(OUT) :: ind_zmat(na,5)
       real(KREAL),INTENT(OUT) ::  val_zmat(na,3)
 
-      character(2) :: ATOM
-      integer(KINT) :: at,long
 
-      real(KREAL) ::  vx,vy,vz,dist
       INTEGER(KINT), ALLOCATABLE :: LIAISONS(:,:) !(na,0:NMaxL
       INTEGER(KINT), ALLOCATABLE :: LiaisonsBis(:,:) ! na,0:NMaxL
       INTEGER(KINT), ALLOCATABLE :: LieS_CP(:,:),LieS_CF(:,:) ! (na,0:NMaxL)
       INTEGER(KINT), ALLOCATABLE :: NbAt3(:,:) !(Na,2)
       INTEGER(KINT), ALLOCATABLE :: CycleAt(:) !(Na)
-      INTEGER(KINT) :: Nbli,Nblj
       INTEGER(KINT), ALLOCATABLE :: CAFaire(:) ! Na
 
       Integer(KINT) :: AtTypCycl(Max_Z)
       Integer(KINT) :: NbCycle
       real(KREAL) ::  vx1,vy1,vz1,norm1
       real(KREAL) ::  vx2,vy2,vz2,norm2
-      real(KREAL) ::  vx3,vy3,vz3,norm3
-      real(KREAL) ::  vx4,vy4,vz4,norm4
-      real(KREAL) ::  vx5,vy5,vz5,norm5
-      real(KREAL) ::  val,val_d
+      real(KREAL) ::  val
       Logical PasFini,AtTerm,Debug,Bicycle,FLie3,FirstCycle
       LOGICAL, ALLOCATABLE ::  Former3(:), DejaFait(:) ! Na
       Logical FTmp
@@ -182,8 +175,6 @@
       END IF
 
 
- 1036 FORMAT(I2)
-
 !     Premiere etape : calcul des connectivites
       DO I=1,na
          DejaFait(I)=.FALSE.
@@ -958,7 +949,7 @@
             END IF
             IndAFaire=IndAFaire+1
          END DO
-12345    CONTINUE
+!12345    CONTINUE
       END DO
 
       if (debug) THEN
Index: src/Step_GEDIIS.f90
===================================================================
--- src/Step_GEDIIS.f90	(revision 1)
+++ src/Step_GEDIIS.f90	(revision 2)
@@ -22,10 +22,10 @@
       LOGICAL DEBUG, PRINT, ci_lt_zero
       INTEGER(KINT), SAVE :: MSET ! mth Iteration
     REAL(KREAL) :: ci(MRESET), ci_tmp(MRESET) ! MRESET = maximum number of iterations.
-      INTEGER(KINT) :: NGEDIIS, MPLUS, INV, ITERA, MM, cis_zero, IXX, JXX, MSET_minus_cis_zero
-      INTEGER(KINT) :: I,J,K, JJ, KJ, JNV, II, IONE, IJ, INK,ITmp, IX, JX, KX, MSET_C_cis_zero
-    INTEGER(KINT) :: current_size_B_mat, MyPointer, Iat
-      REAL(KREAL) :: XMax, XNorm, S, DET, THRES, tmp, ER_star, ER_star_tmp
+      INTEGER(KINT) :: NGEDIIS, MPLUS, INV, ITERA, MM, cis_zero
+      INTEGER(KINT) :: I, J, K, JJ, JNV, II, IONE, IJ, IX, JX, KX
+    INTEGER(KINT) :: current_size_B_mat, MyPointer
+      REAL(KREAL) :: XMax, XNorm, DET, THRES, tmp, ER_star, ER_star_tmp
 
       DEBUG=.TRUE.
       PRINT=.FALSE.
Index: src/EgradPath.f90
===================================================================
--- src/EgradPath.f90	(revision 1)
+++ src/EgradPath.f90	(revision 2)
@@ -1,4 +1,4 @@
-SUBROUTINE EGradPath(Iopt,Flag_Opt_Geom)
+SUBROUTINE EGradPath(Iopt)
 
   ! This subroutine calculates the energy+gradient for all points of
   ! the path
@@ -13,7 +13,6 @@
   IMPLICIT NONE
 
   INTEGER(KINT), INTENT(IN) :: Iopt
-  LOGICAL, INTENT(IN) :: Flag_Opt_Geom
 
   INTEGER(KINT) :: IGeom,IGeom0,IGeomF
   CHARACTER(LCHARS) :: Line,FileInp,RunCommand
Index: src/Hinvup_DFP.f90
===================================================================
--- src/Hinvup_DFP.f90	(revision 1)
+++ src/Hinvup_DFP.f90	(revision 2)
@@ -46,7 +46,7 @@
   REAL(KREAL), ALLOCATABLE  :: Hup(:,:)
   REAL(KREAL), ALLOCATABLE  :: Num(:,:), Num2(:,:)
   REAL(KREAL) :: Den, Den2,NormDS
-  INTEGER(KINT) :: i,j,k, Idx
+  INTEGER(KINT) :: i, j, Idx
 
   LOGICAL :: debug
 
Index: src/addligne.f90
===================================================================
--- src/addligne.f90	(revision 1)
+++ src/addligne.f90	(revision 2)
@@ -2,8 +2,8 @@
 
       use Path_module, only : Nat, KINT, KREAL
 
-      integer(KINT) :: na,atome(Nat),at,ind_zmat(Nat,5),long
-      real(KREAL) ::  x(Nat),y(Nat),z(Nat),fact,val_zmat(Nat,3)
+      integer(KINT) :: ind_zmat(Nat, 5)
+      real(KREAL) ::  x(Nat), y(Nat), z(Nat), val_zmat(Nat, 3)
 
       real(KREAL) :: vx1,vy1,vz1,norm1
       real(KREAL) ::  vx2,vy2,vz2,norm2
@@ -21,11 +21,11 @@
            vx2,vy2,vz2,norm2)
 
       CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
-      CALL produit_vect(vx1,vy1,vz1,norm1, &
-           vx2,vy2,vz2,norm2, &
+      CALL produit_vect(vx1,vy1,vz1, &
+           vx2,vy2,vz2, &
            vx4,vy4,vz4,norm4)
-      CALL produit_vect(vx3,vy3,vz3,norm3, &
-           vx2,vy2,vz2,norm2, &
+      CALL produit_vect(vx3,vy3,vz3, &
+           vx2,vy2,vz2, &
            vx5,vy5,vz5,norm5)
 
       val_d=angle_d(vx4,vy4,vz4,norm4, &
@@ -33,7 +33,7 @@
            vx2,vy2,vz2,norm2)
 
 !     write(*,11) n1,n2,norm1,n3,val,n4,val_d
- 11   format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))
+! 11   format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))
 
       ind_zmat(i,1)=n1
       ind_zmat(i,2)=n2
Index: src/ConvXyzMixed.f90
===================================================================
--- src/ConvXyzMixed.f90	(revision 1)
+++ src/ConvXyzMixed.f90	(revision 2)
@@ -12,7 +12,7 @@
   REAL(KREAL), INTENT(INOUT) :: val_zmat(nb,3)
 
   integer(KINT) :: n1,n2,n3,n4
-  integer(KINT) :: i,i1,i2,ibeg
+  integer(KINT) :: i, ibeg
 
   real(KREAL) ::  vx1,vy1,vz1,norm1
   real(KREAL) ::  vx2,vy2,vz2,norm2
@@ -119,12 +119,12 @@
 
      !        write(*,*) 'Debut procedure Xyz_Zmat43'
      CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
-     CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2,vx4,vy4,vz4,norm4)
-     CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2,vx5,vy5,vz5,norm5)
+     CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
+     CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,vx5,vy5,vz5,norm5)
      val_d=angle_d(vx4,vy4,vz4,norm4, vx5,vy5,vz5,norm5,vx2,vy2,vz2,norm2)
 
      !               write(*,11) n1,n2,norm1,n3,val,n4,val_d
-11   format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))
+!11   format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))
 
      val_zmat(i,1)=norm1
      val_zmat(i,2)=val
Index: src/AlignPartial.f90
===================================================================
--- src/AlignPartial.f90	(revision 1)
+++ src/AlignPartial.f90	(revision 2)
@@ -35,11 +35,11 @@
   Real(KREAL) ::  xc1,yc1,zc1, xc2,yc2,zc2
 
 
-  integer(KINT) :: i, j, nd, ia
+  integer(KINT) :: i, j, ia
   Real(KREAL) :: x_norm, y_norm, lambda
   Real(KREAL) :: Rmatrix(3,3)
-  Real(KREAL) :: S(4,4), dS(4,4), q(4)
-  Real(KREAL) :: tmp(3),tmp1(4),tmp2(4),EigVec(4,4), EigVal(4)
+  Real(KREAL) :: S(4, 4)
+  Real(KREAL) :: EigVec(4, 4), EigVal(4)
 
   INTERFACE
      function valid(string) result (isValid)
Index: src/Extrapol_cart.f90
===================================================================
--- src/Extrapol_cart.f90	(revision 1)
+++ src/Extrapol_cart.f90	(revision 2)
@@ -21,13 +21,12 @@
   REAL(KREAL), INTENT(IN) :: Xgeom(NGeomI),Coef(NGeomI,NCoord)
 ! Iopt: Number of the cycles for the optimization
   INTEGER(KINT), INTENT(IN) :: Iopt
-  INTEGER(KINT) :: IdxGeom, I, J, K, Idx,Iat, IGeom
+  INTEGER(KINT) :: IdxGeom, I, J, K, Iat
 ! NSpline is the number of points along the interpolating path
   INTEGER(KINT) :: NSpline
 ! FileSpline: Filename to save the interpolating path coordinates
   CHARACTER(LCHARS) :: FileSpline,TmpChar
-  REAL(KREAL) :: Rmsd,MRot(3,3), ds, u, v
-  REAL(KREAL) :: a_val, d
+  REAL(KREAL) :: Rmsd, MRot(3, 3), ds, u
 
   REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)
   REAL(KREAL), ALLOCATABLE :: DerXyz(:,:) ! Nat,3
Index: src/egrad_chamfre.f90
===================================================================
--- src/egrad_chamfre.f90	(revision 1)
+++ src/egrad_chamfre.f90	(revision 2)
@@ -22,7 +22,7 @@
        REAL(KREAL), INTENT(IN) :: Geom(Nat,3)
 
 !     Bohr --> Angstr
-      real(KREAL), parameter :: BOHR   = 0.52917726D+00
+!      real(KREAL), parameter :: BOHR   = 0.52917726D+00
 
 !! PFL 2010 Nov. 28 ->
 ! for now I keep this old F77 scheme. This has to be replaced by ALLOCATABLE
@@ -30,13 +30,12 @@
 ! <- PFL 2010 Nov 28.
   REAL(KREAL) :: xa(100),ya(100),za(100),dvdr(100,3)
   REAL(KREAL) :: rshort,rlong,zshort,zlong
-  REAL(KREAL), SAVE :: rnc,rnn,rhh,rch,rnh,ann(100),anc(100),anh(100),ach(100),ahh(100),a(100)
+  REAL(KREAL), SAVE :: rnc, rnn, rhh, rch, rnh, ann(100), anc(100), anh(100), ach(100), ahh(100)
   REAL(KREAL) :: dv(100,10,3),df(100,10,3)
   REAL(KREAL) :: px(9000),py(9000),pz(9000),tmp,temp,rtemp,yf
   INTEGER(KINT) :: ka(9000),kx(9000),ky(9000),itx,ity,nf(100),mf(100,500)
-  REAL(KREAL) :: xi(3),xj(3),xk(3),rij,rik,rjk,rji,rki,er,eb
-  REAL(KREAL) :: fcut,tf,bodij,bodji,bod,cost,gtheta,dgdt,xt,t1,t2,t3
-  REAL(KREAL) :: zi,fcut2
+  REAL(KREAL) :: xi(3), xj(3), rij, er, eb
+  REAL(KREAL) :: fcut, tf
   INTEGER(KINT) :: np,nq,nh,nc
   INTEGER(KINT) :: i, j, nn, m, n, mid, n1, im, jm, n2,n3,n4, iat
   LOGICAL, SAVE :: First=.TRUE.
@@ -584,7 +583,7 @@
 subroutine fctNiNi(r2,r,fcut)
 
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
 
  REAL(KREAL) :: r2,r,fcut,r1,pi
@@ -605,7 +604,7 @@
 !!************
 subroutine dfcutNiNi(r2,r,df)
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
   REAL(KREAL) :: r, r2,df,pi,r1
 pi=4.0*atan(1.0)
@@ -627,7 +626,7 @@
 !!************
 subroutine fctNiH(r2,r,fcut)
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
 
  REAL(KREAL) ::r2,r,fcut
@@ -655,7 +654,7 @@
 
 
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
 
 
@@ -683,7 +682,7 @@
 subroutine fctNiC(r2,r,fcut)
 
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
 
 
@@ -711,7 +710,7 @@
 subroutine dfcutNiC(r2,r,df)
 
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
 
      
@@ -740,7 +739,7 @@
 
 
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
 
 
@@ -769,7 +768,7 @@
 
 
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
 
      
@@ -796,7 +795,7 @@
 subroutine fctHH(r2,r,fcut)
 
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
 
 
@@ -819,7 +818,7 @@
 
 
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
 
 
@@ -844,7 +843,7 @@
 
 
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
 
 
@@ -871,7 +870,7 @@
 subroutine fct2(zshort,zlong,r,fcut)
 
    IMPLICIT NONE
-     integer, parameter :: KINT = kind(1)
+!     integer, parameter :: KINT = kind(1)
      integer, parameter :: KREAL = kind(1.0d0)
 
 
Index: src/Step_RFO_all.f90
===================================================================
--- src/Step_RFO_all.f90	(revision 1)
+++ src/Step_RFO_all.f90	(revision 2)
@@ -69,17 +69,17 @@
 
   REAL(KREAL), PARAMETER :: tiny=1e-20, zero=0., dele3=1e-6, eps=1e-12
   REAL(KREAL), PARAMETER :: crit=1e-8
-  INTEGER(KINT) :: n, i, j, isignx, im, idx,k, isch
-  REAL(KREAL) :: grd, evl, eval,e1,e2, norm
-  REAL(KREAL) :: dg, dx
-  CHARACTER(120) :: fmt, fmt2, Line
+  INTEGER(KINT) :: i, j, idx, k, isch
+  REAL(KREAL) :: grd, eval, norm
+  REAL(KREAL) :: dx
+  CHARACTER(120) :: fmt2
 
   LOGICAL :: Debug
 
   REAL(KREAL), ALLOCATABLE :: HFree(:,:) ! NFree,NFree
   REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree
 
-  LOGICAL, SAVE :: LFirst=.TRUE.
+!  LOGICAL, SAVE :: LFirst=.TRUE.
 
   INTERFACE
      function valid(string) result (isValid)
@@ -168,6 +168,7 @@
         STOP
      END IF
      
+     if (debug) WRITE(*,*) "Deallocating Tanf"
      DEALLOCATE(Tanf)
 
      NFree=Idx
@@ -252,7 +253,20 @@
   END DO
 ! <- PFL 04 Apr 2008
      
-  DEALLOCATE(eigval,eigvec,eigvli,hfree,htmp,grad_f,step_f) 
+  if (debug) WRITE(*,*) "Deallocataing Eigval"
+  DEALLOCATE(eigval) 
+  if (debug) WRITE(*,*) "Deallocataing Eigvec"
+  DEALLOCATE(eigvec) 
+  if (debug) WRITE(*,*) "Deallocataing Eigvli"
+  DEALLOCATE(eigvli) 
+  if (debug) WRITE(*,*) "Deallocataing hfree"
+  DEALLOCATE(hfree) 
+  if (debug) WRITE(*,*) "Deallocataing htmp"
+  DEALLOCATE(htmp) 
+  if (debug) WRITE(*,*) "Deallocataing grad_f"
+  DEALLOCATE(grad_f) 
+  if (debug) WRITE(*,*) "Deallocataing step_f"
+  DEALLOCATE(step_f) 
 
   if (debug) THEN
      WRITE(*,*) 'DBG Step_RFO_All Step'
Index: src/Opt_Geom.f90
===================================================================
--- src/Opt_Geom.f90	(revision 1)
+++ src/Opt_Geom.f90	(revision 2)
@@ -289,7 +289,7 @@
                      ! GradTmp now becomes "step" below:
                   GradTmp = Geom - GeomOld
               CASE ('GDIIST')
-          Call Step_GDIIS_Simple_Err(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)
+          Call Step_GDIIS_Simple_Err(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,NCoord,FRST)
                ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
                      Geom=0.d0
             DO I=1, NCoord
@@ -321,7 +321,7 @@
         if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep
         GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all.
 
-        Call Check_step(IGeom,Coord,Nat,NCoord,GeomOld,GradTmp)
+        Call Check_step(Coord,Nat,NCoord,GeomOld,GradTmp)
        
         Geom=GeomOld+GradTmp ! GradTmp is step size here, from Step_RFO_all.
 
Index: src/IntCoord_der.f90
===================================================================
--- src/IntCoord_der.f90	(revision 1)
+++ src/IntCoord_der.f90	(revision 2)
@@ -30,10 +30,8 @@
       REAL(KREAL)   ,INTENT(IN) :: R0(3,NAT)
       REAL(KREAL)   ,INTENT(OUT):: B(3,NAT)
 
-      INTEGER(KINT)            :: I,J
-      REAL(KREAL)               :: DIS,DI,DJ
-      REAL(KREAL)               :: DELTAK
-      REAL(KREAL)               :: SVAR
+      INTEGER(KINT)            :: I
+      REAL(KREAL)               :: DIS, DI
       real(KREAL)               :: dr(3)
 !     ******************************************************************
       DR(:)=R0(:,IAT2)-R0(:,IAT1)
@@ -94,11 +92,11 @@
       INTEGER(KINT),INTENT(IN) :: IAT3
       REAL(KREAL)   ,INTENT(OUT):: B(3,NAT)
       REAL(KREAL)   ,INTENT(IN) :: R0(3,NAT)
-      REAL(KREAL)               :: X(9),DI(9)
+      REAL(KREAL)               :: X(9)
       real(KREAL)               :: phi
       REAL(KREAL)               :: X1X2,Y1Y2,Z1Z2,X2X3,Y2Y3,Z2Z3,D12,D22,  &
                                    D1,D2,D13,D23
-      integer(KINT)            :: i,j
+      integer(KINT)            :: i
 !     ******************************************************************
 !
 ! --- TRANSFER R0 TO X ARRAY TO MAKE HANDLING AND PASSING EASIER
@@ -155,9 +153,9 @@
       INTEGER(KINT),INTENT(IN) :: IAT3
       REAL(KREAL)   ,INTENT(OUT):: B(3,NAT)
       REAL(KREAL)   ,INTENT(IN) :: R0(3,NAT)
-      REAL(KREAL)               :: X(9),DI(9),DIDJ(9,9),B1(9)
+      REAL(KREAL)               :: X(9), DI(9)
       real(KREAL)               :: phi
-      integer(KINT)            :: i,j,ichk, iat
+      integer(KINT)            :: i
       REAL(KREAL)               :: UU    ! U*U
       REAL(KREAL)               :: VV    ! V*V
       REAL(KREAL)               :: UV    ! U*V
@@ -167,11 +165,9 @@
       REAL(KREAL)               :: DCDX(9)
       REAL(KREAL)               :: DDDX(9)
       REAL(KREAL)               :: DJDX(9)
-      REAL(KREAL)               :: DVDX(9)
       REAL(KREAL)               :: DGDX(9)    ! D(A*B)/DX
       REAL(KREAL)               :: AB,AB2N,V,PREFAC
       REAL(KREAL)               :: COSPHI,SINPHI
-      REAL(KREAL)               :: DIVDX,AI
 !
 ! Variable used in case of a linear molecule with three atoms call A,B and C
 ! Phi is the ABC angle
@@ -382,17 +378,14 @@
       INTEGER(KINT),INTENT(IN) :: IAT4
       REAL(KREAL)   ,INTENT(IN) :: R0(3,NAT)
       REAL(KREAL)   ,INTENT(OUT):: B(3,NAT)
-      REAL(KREAL)               :: X(12),DI(12),DIDJ(12,12),B1(12)
+      REAL(KREAL)               :: X(12), DI(12)
       REAL(KREAL)               :: P(3),Q(3),T(3),U(3),R(3),W(3)
-      real(KREAL)               :: phi
-      integer(KINT)            :: i,j,ichk,k
+      integer(KINT)            :: i, j
 
       REAL(KREAL)   :: TORS, PP,QQ, PQ
       REAL(KREAL) :: DADX(12),DBDX(12),DCDX(12)
-      REAL(KREAL) :: DDDX(12),DJDX(12),DVDX(12),DGDX(12),DPDX(3,12),DQDX(3,12)
-      REAL(KREAL) :: D2QDX2(3,12,12),D2PDX2(3,12,12)
-      REAL(KREAL) :: D2BDX2(12,12),D2ADX2(12,12),D2CDX2(12,12),D2DDX2(12,12)
-      real(KREAL) :: d,coverd,abyb2n,abyb,ai,divdx,prefac,v
+      REAL(KREAL) :: DDDX(12), DJDX(12), DGDX(12), DPDX(3, 12), DQDX(3, 12)
+      real(KREAL) :: d, coverd, abyb2n, abyb, prefac, v
 
 !     ******************************************************************
 !
@@ -572,20 +565,15 @@
       INTEGER(KINT),INTENT(IN) :: IAT4
       REAL(KREAL)   ,INTENT(IN) :: R0(3,NAT)
       REAL(KREAL)   ,INTENT(OUT):: B(3,NAT)
-      REAL(KREAL)               :: X(12),DI(12),DIDJ(12,12),B1(12)
-      REAL(KREAL)               :: P(3),Q(3),T(3),U(3),R(3),W(3)
+      REAL(KREAL)               :: X(12), DI(12)
+      REAL(KREAL)               :: T(3)
       REAL(KREAL)  :: Rij(3), rkj(3), rkl(3),rmj(3),rnk(3)
-      REAL(KREAL)  :: Nrij, Nrkj, Nrkl, Nrmj, Nrnk
+      REAL(KREAL)  :: Nrkj, Nrmj, Nrnk
       REAL(KREAL) :: SignPhi
-      real(KREAL)               :: phi
-      integer(KINT)            :: i,j,ichk,k
+      integer(KINT)            :: i
 
-      REAL(KREAL)   :: TORS, PP,QQ, PQ
-      REAL(KREAL) :: DADX(12),DBDX(12),DCDX(12)
-      REAL(KREAL) :: DDDX(12),DJDX(12),DVDX(12),DGDX(12),DPDX(3,12),DQDX(3,12)
-      REAL(KREAL) :: D2QDX2(3,12,12),D2PDX2(3,12,12)
-      REAL(KREAL) :: D2BDX2(12,12),D2ADX2(12,12),D2CDX2(12,12),D2DDX2(12,12)
-      real(KREAL) :: d,coverd,abyb2n,abyb,ai,divdx,prefac,v
+      REAL(KREAL)   :: TORS
+      real(KREAL) :: coverd
 
       LOGICAL :: Debug=.TRUE.
 
@@ -691,11 +679,8 @@
       INTEGER(KINT),INTENT(IN) :: IAT4
       REAL(KREAL)   ,INTENT(IN) :: R0(3,NAT)
       REAL(KREAL)   ,INTENT(OUT):: DZDC(3,NAT)
-      REAL(KREAL)               :: DI(12)
       REAL(KREAL)               :: A(3),B(3),C(3),Vec(3)
       REAL(KREAL)               :: U(3),V(3)
-      real(KREAL)               :: phi
-      integer(KINT)            :: i,j,ichk,k
 
       REAL(KREAL)   :: norm1,norm2,norm3
       real(KREAL) :: s,bleng
@@ -727,7 +712,7 @@
 ! --- We check that PQR and QRS are not aligned.
 ! If they are, we return 0. as derivatives instead of just crashing.
 ! ----------------------------------------------------------------
-    call produit_vect (b(1),b(2),b(3),Norm1,c(1),c(2),c(3),Norm2, vec(1),vec(2),vec(3),Norm3)
+    call produit_vect (b(1),b(2),b(3),c(1),c(2),c(3),vec(1),vec(2),vec(3),Norm3)
      norm1=sqrt(dot_product(b,b))
      norm2=sqrt(dot_product(c,c))
     Failed=.FALSE.
@@ -737,7 +722,7 @@
            write (*,*) 'QRS aligned',iat2,iat3,iat4
            Failed=.TRUE.
      end if
-     call produit_vect (b(1),b(2),b(3),Norm1,a(1),a(2),a(3),Norm2, vec(1),vec(2),vec(3),Norm3)
+     call produit_vect (b(1),b(2),b(3),a(1),a(2),a(3),vec(1),vec(2),vec(3),Norm3)
      if (norm3 <= eps) then
            WRITE(*,*) '*** WARNING: ILL-DEFINED DIHEDRAL ANGLE(S)'
            write (*,*) 'PQR aligned',iat1,iat2,iat3
@@ -749,28 +734,27 @@
 ! ----------------------------------------------------------------
 
      bleng = sqrt (b(1)**2 + b(2)**2 + b(3)**2)
-     call produit_vect (a(1),a(2),a(3),Norm1,b(1),b(2),b(3),Norm2, vec(1),vec(2),vec(3),Norm3)
+     call produit_vect (a(1),a(2),a(3),b(1),b(2),b(3),vec(1),vec(2),vec(3),Norm3)
      s = bleng / (vec(1)**2 + vec(2)**2 + vec(3)**2)
 !     dzdc(:,ip,3,ip) = s * vec
       dzdc(:,iat1) = s * vec
 
-     call produit_vect (b(1),b(2),b(3),Norm1,c(1),c(2),c(3),Norm2, vec(1),vec(2),vec(3),Norm3)
+     call produit_vect (b(1),b(2),b(3),c(1),c(2),c(3),vec(1),vec(2),vec(3),Norm3)
      s = bleng / (vec(1)**2 + vec(2)**2 + vec(3)**2)
 !    dzdc(:,is,3,ip) = s * vec
      dzdc(:,iat4) = s * vec
 
 
-     call produit_vect (c(1),c(2),c(3),Norm1,    &
-          dzdc(1,iat4),  dzdc(2,iat4),  dzdc(3,iat4), Norm2, &
+     call produit_vect (c(1),c(2),c(3),   &
+          dzdc(1,iat4),  dzdc(2,iat4),  dzdc(3,iat4), &
           u(1),u(2),u(3),Norm3)
      vec = a - b
-     call produit_vect (vec(1),vec(2),vec(3),Norm1,    &
-          dzdc(1,iat1),  dzdc(2,iat1),  dzdc(3,iat1), Norm2, &
+     call produit_vect (vec(1),vec(2),vec(3),    &
+          dzdc(1,iat1),  dzdc(2,iat1),  dzdc(3,iat1),  &
           v(1),v(2),v(3),Norm3)
      vec = u + v
-     call produit_vect (vec(1),vec(2),vec(3),Norm1,    &
-          b(1), b(2), b(3), Norm2, &
-          u(1),u(2),u(3),Norm3)
+     call produit_vect (vec(1),vec(2),vec(3),    &
+          b(1), b(2), b(3), u(1),u(2),u(3),Norm3)
 
 !     dzdc(:,iq,3,ip) = u / bleng**2
 !     dzdc(:,ir,3,ip) = - (dzdc(:,ip,3,ip) + dzdc(:,iq,3,ip) + dzdc(:,is,3,ip))
Index: src/Step_GDIIS_Simple_Err.f90
===================================================================
--- src/Step_GDIIS_Simple_Err.f90	(revision 1)
+++ src/Step_GDIIS_Simple_Err.f90	(revision 2)
@@ -1,8 +1,8 @@
 !C  HEAT is never used, not even in call of Space(...)
 !C  Geom = input parameter vector (Geometry).
 !C  Grad = input gradient vector.
-      SUBROUTINE Step_GDIIS_Simple_Err(NewGeom,Geom,NewGrad,GRAD,HP,HEAT,Hess,NCoord,FRST)
-!      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SUBROUTINE Step_GDIIS_Simple_Err(NewGeom,Geom,NewGrad,GRAD,HP,HEAT,NCoord,FRST)
+
       IMPLICIT NONE
       integer, parameter :: KINT = kind(1)
       integer, parameter :: KREAL = kind(1.0d0)
@@ -10,7 +10,7 @@
 !      INCLUDE 'SIZES'
 
       INTEGER(KINT) :: NCoord
-      REAL(KREAL) :: NewGeom(NCoord), Geom(NCoord), NewGrad(NCoord), GRAD(NCoord), Hess(NCoord*NCoord)
+      REAL(KREAL) :: NewGeom(NCoord), Geom(NCoord), NewGrad(NCoord), GRAD(NCoord)
       REAL(KREAL) :: HEAT, HP
       LOGICAL :: FRST
     
@@ -57,12 +57,12 @@
       REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:), GradSet(:), ERR(:) ! MRESET*NCoord
       REAL(KREAL) :: ESET(MRESET)
       REAL(KREAL), ALLOCATABLE, SAVE :: DX(:), GSAVE(:) !NCoord
-      REAL(KREAL) :: B(M2), BS(M2), BST(M2)
+      REAL(KREAL) :: B(M2), BS(M2)
       LOGICAL DEBUG, PRINT
       INTEGER(KINT), SAVE :: MSET
       INTEGER(KINT) :: NDIIS, MPLUS, INV, ITERA, MM, I, J, K
-      INTEGER(KINT) :: JJ, KJ, JNV, II, IONE, IJ, INK,ITmp
-      REAL(KREAL) :: XMax, XNorm, S, DET, THRES
+      INTEGER(KINT) :: JJ, JNV, II, IONE, IJ, INK, ITmp
+      REAL(KREAL) :: XMax, XNorm, DET, THRES
 
       DEBUG=.TRUE.
       PRINT=.TRUE.
Index: src/Calc_baker_allGeomF.f90
===================================================================
--- src/Calc_baker_allGeomF.f90	(revision 1)
+++ src/Calc_baker_allGeomF.f90	(revision 2)
@@ -27,21 +27,11 @@
   REAL(KREAL), ALLOCATABLE :: XPrimRef(:) ! NPrim
   INTEGER(KINT) :: IGeom
   
-  real(KREAL) ::  vx,vy,vz,dist, Norm
-  real(KREAL) ::  vx1,vy1,vz1,norm1
-  real(KREAL) ::  vx2,vy2,vz2,norm2
-  real(KREAL) ::  vx3,vy3,vz3,norm3
-  real(KREAL) ::  vx4,vy4,vz4,norm4
-  real(KREAL) ::  vx5,vy5,vz5,norm5
-  real(KREAL) ::  val,val_d, Q, T
 
-  INTEGER(KINT) :: I,J, n1,n2,n3,n4,IAt,IL,JL,IFrag,ITmp, K, KMax
-  INTEGER(KINT) :: I0, IOld, IAtTmp, Izm, JAt, Kat, Lat, L, NOUT
+  INTEGER(KINT) :: I, J, K
 
-  REAL(KREAL) :: sAngleIatIKat, sAngleIIatLat
-  REAL(KREAL) :: DiheTmp
 
-  LOGICAL :: debug, bond, AddPrimitiveCoord, FAIL
+  LOGICAL :: debug
   LOGICAL :: DebugPFL
 
   INTERFACE
@@ -120,42 +110,6 @@
      z(1:Nat) = XyzGeomF(IGeom,3,1:Nat)
      XPrimREf=XPrimitiveF(IGeom,:)
      Call Calc_XPrim(nat,x,y,z,Coordinate,NPrim,XPrimitiveF(IGeom,:),XPrimRef)
-!    ScanCoord=>Coordinate
-!    I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
-!    DO WHILE (Associated(ScanCoord%next))
-!      I=I+1
-!      SELECT CASE (ScanCoord%Type)
-!       CASE ('BOND')
-!       Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)
-!            XprimitiveF(IGeom,I) = Norm2  
-!       CASE ('ANGLE')
-!     Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx1,vy1,vz1,Norm1)
-!     Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)
-!          XprimitiveF(IGeom,I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
-!     CASE ('DIHEDRAL')
-!     Call vecteur(ScanCoord%At3,ScanCoord%At2,x,y,z,vx2,vy2,vz2,Norm2)
-!     Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1)
-!     Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3)
-!     Call produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
-!                         vx5,vy5,vz5,norm5)
-!     Call produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
-!                         vx4,vy4,vz4,norm4)
-!                 DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
-!                         vx2,vy2,vz2,norm2)
-!                 XprimitiveF(IGeom,I) = DiheTmp*Pi/180.
-!! We treat large dihedral angles differently as +180=-180 mathematically and physically
-!! but this causes lots of troubles when converting baker to cart.
-!! So we ensure that large dihedral angles always have the same sign
-!                    if (abs(ScanCoord%SignDihedral).NE.1) THEN
-!                        ScanCoord%SignDihedral=Int(Sign(1.d0,DiheTmp))
-!                    ELSE
-!                         If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN
-!                          XprimitiveF(IGeom,I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi
-!                        END IF
-!                    END IF
-!      END SELECT
-!      ScanCoord => ScanCoord%next
-!    END DO ! matches DO WHILE (Associated(ScanCoord%next))
   END DO ! matches DO IGeom=1, NGeomF
    
   ALLOCATE(BprimT(3*Nat,NPrim))
Index: src/Calc_Xprim.f90
===================================================================
--- src/Calc_Xprim.f90	(revision 1)
+++ src/Calc_Xprim.f90	(revision 2)
@@ -34,17 +34,14 @@
 
   Type (ListCoord), POINTER :: ScanCoord
 
-  real(KREAL) :: vx,vy,vz,dist, Norm
   real(KREAL) :: vx1,vy1,vz1,norm1
   real(KREAL) :: vx2,vy2,vz2,norm2
   real(KREAL) :: vx3,vy3,vz3,norm3
   real(KREAL) :: vx4,vy4,vz4,norm4
   real(KREAL) :: vx5,vy5,vz5,norm5
-  real(KREAL) :: val,val_d, Q, T
 
-  INTEGER(KINT) :: I,J, n1,n2,n3,n4
+  INTEGER(KINT) :: I
 
-  REAL(KREAL) :: sAngleIatIKat, sAngleIIatLat
   REAL(KREAL) :: DiheTmp
 
   LOGICAL :: debug, debugPFL
@@ -119,9 +116,9 @@
            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1)
            Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2)
            Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3)
-           Call produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
+           Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
                 vx4,vy4,vz4,norm4)
-           Call produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
+           Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
                 vx5,vy5,vz5,norm5)
 
                  DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
Index: src/Hupdate_Bofill.f90
===================================================================
--- src/Hupdate_Bofill.f90	(revision 1)
+++ src/Hupdate_Bofill.f90	(revision 2)
@@ -27,7 +27,6 @@
   REAL(KREAL), ALLOCATABLE  :: HupInt(:)
   REAL(KREAL), ALLOCATABLE  :: HS(:), Num(:,:), Num2(:,:)
   REAL(KREAL) :: Phi, Num3, Den, Den2,NormDS
-  CHARACTER(120) :: fmt,fmt2,Line
   INTEGER(KINT) :: i,j, idx
 
   LOGICAL debug,valid
Index: src/CalcTangent.f90
===================================================================
--- src/CalcTangent.f90	(revision 1)
+++ src/CalcTangent.f90	(revision 2)
@@ -18,14 +18,14 @@
   REAL(KREAL), ALLOCATABLE :: PathTan(:,:) ! NGeomF,NCoord
   REAL(KREAL), ALLOCATABLE :: DbgTgt(:),GeomTmp(:) ! NCoord
 
-  INTEGER(KINT) :: i, j, n3at, Idx, iat,igeom,k
+  INTEGER(KINT) :: i, j, n3at, k
   INTEGER(KINT), SAVE :: icall=-1
   INTEGER(KINT) :: NMaxPt
-  REAL(KREAL), PARAMETER :: Inf=1e32, h=0.001d0
+!  REAL(KREAL), PARAMETER :: Inf=1e32, h=0.001d0
   CHARACTER(120) :: Line,File
   CHARACTER(LCHARS) :: TITLE
   LOGICAL :: Debug=.False.
-  REAL(KREAL) :: utmp,der
+  REAL(KREAL) :: utmp
 
   LOGICAL, EXTERNAL :: valid
 
Index: src/Space_all.f90
===================================================================
--- src/Space_all.f90	(revision 1)
+++ src/Space_all.f90	(revision 2)
@@ -23,7 +23,7 @@
       REAL(KREAL) :: GeomSet(NGeomF,MRESET*NCoord),GradSet(NGeomF,MRESET*NCoord), ESET(MRESET)
         ! Geom and GRAD for all iterations.
 
-      INTEGER(KINT) :: I,J,K,NMK,MI,NI
+      INTEGER(KINT) :: I, K, NMK, MI, NI
       INTEGER(KINT), SAVE :: NRESET
 !C
 !C     UPDATE PARAMETER AND GRADIENT SUBSPACE
Index: src/Extrapol_mixed.f90
===================================================================
--- src/Extrapol_mixed.f90	(revision 1)
+++ src/Extrapol_mixed.f90	(revision 2)
@@ -25,8 +25,7 @@
   REAL(KREAL), INTENT(IN) :: Xgeom(NGeomI),Coef(NGeomI,NCoord)
 
   INTEGER(KINT) :: IdxGeom, I, J, K, Idx
-  REAL(KREAL) :: Rmsd,MRot(3,3), ds, u, v
-  REAL(KREAL) :: a_val, d
+  REAL(KREAL) :: Rmsd, MRot(3, 3), ds, u
 
   REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)
   REAL(KREAL), ALLOCATABLE :: IntCoordTmp(:),DerInt(:) ! (NCoord)
Index: src/CalcHess.f90
===================================================================
--- src/CalcHess.f90	(revision 1)
+++ src/CalcHess.f90	(revision 2)
@@ -50,9 +50,9 @@
       logical           :: debug
 
       REAL(KREAL) :: Ener
-      INTEGER(KINT) :: I, J, IAT, JAT, IJ0, IJ1, IConst
+      INTEGER(KINT) :: I, J, IAT, IConst
 
-       REAL(KREAL), ALLOCATABLE :: Temp(:),Temp2(:,:) 
+       REAL(KREAL), ALLOCATABLE :: Temp2(:, :) 
 
 
        REAL(KREAL), ALLOCATABLE :: Value(:), ValueOld(:), GradTMP(:)
Index: src/Int2cart.f90
===================================================================
--- src/Int2cart.f90	(revision 1)
+++ src/Int2cart.f90	(revision 2)
@@ -12,7 +12,7 @@
   REAL(KREAL), INTENT(IN) :: IntCoord(3*Nat-6)
   REAL(KREAL), INTENT(OUT) :: XyzCoord(Nat,3)
 
-  INTEGER(KINT) :: I,J,K,L,Idx
+  INTEGER(KINT) :: I
   REAl(KREAL), ALLOCATABLE :: ValZmat(:,:),XyzTmp(:,:)
   REAL(KREAL) :: d, a_val
   LOGICAL :: debug
Index: src/Space.f90
===================================================================
--- src/Space.f90	(revision 1)
+++ src/Space.f90	(revision 2)
@@ -23,7 +23,7 @@
       REAL(KREAL) :: GeomSet(MRESET*NVAR), GradSet(MRESET*NVAR), ESET(MRESET)
         ! Geom and GRAD for all iterations.
 
-      INTEGER(KINT) :: I,J,K,NMK,MI,NI
+      INTEGER(KINT) :: I, K, NMK, MI, NI
       INTEGER(KINT), SAVE :: NRESET
 !C
 !C     UPDATE PARAMETER AND GRADIENT SUBSPACE
Index: src/ConvertBakerInternal_cart.f90
===================================================================
--- src/ConvertBakerInternal_cart.f90	(revision 1)
+++ src/ConvertBakerInternal_cart.f90	(revision 2)
@@ -43,7 +43,7 @@
 
 
   ! IdxGeom: geometry index.
-  Integer(KINT) :: n1,n2,n3,I,J,K,Counter
+  Integer(KINT) :: I, J, Counter
   REAL(KREAL) :: normDeltaX_int
 
   REAL(KREAL), ALLOCATABLE :: DeltaX_int(:)  !DeltaX_int(NPrim,3*Nat-6), 3*Nat-6??
@@ -53,17 +53,8 @@
   REAL(KREAL), ALLOCATABLE :: GMat(:,:) !(NPrim,NPrim)
   REAL(KREAL), ALLOCATABLE :: EigVec(:,:), EigVal(:) ! EigVec(NPrim,NPrim)
   REAL(KREAL), ALLOCATABLE :: UMat_tmp(:,:)
-  REAL(KREAL), ALLOCATABLE :: V(:,:) ! (NCoord,NCoord)
 
-  REAL(KREAL) ::  vx,vy,vz,dist,Norm
-  REAL(KREAL) :: vx1,vy1,vz1,norm1
-  REAL(KREAL) :: vx2,vy2,vz2,norm2
-  REAL(KREAL) :: vx3,vy3,vz3,norm3
-  REAL(KREAL) :: vx4,vy4,vz4,norm4
-  REAL(KREAL) :: vx5,vy5,vz5,norm5
-  REAL(KREAL) :: val,val_d, Q, T
-  REAL(KREAL) :: Angle, angle_d, tol,Delta,DiheTmp
-  REAL(KREAL) :: MRot(3,3), Rmsd
+  REAL(KREAL) :: Angle, angle_d
   REAL(KREAL), PARAMETER :: Crit=1e-08
 
   EXTERNAL Angle, angle_d
@@ -193,8 +184,8 @@
 !            Call vecteur(ScanCoord%At2,ScanCoord%At1,x_k,y_k,z_k,vx1,vy1,vz1,Norm1)
 !            Call vecteur(ScanCoord%At2,ScanCoord%At3,x_k,y_k,z_k,vx2,vy2,vz2,Norm2)
 !            Call vecteur(ScanCoord%At3,ScanCoord%At4,x_k,y_k,z_k,vx3,vy3,vz3,Norm3)
-!            Call produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2,vx5,vy5,vz5,norm5)
-!            Call produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2,vx4,vy4,vz4,norm4)
+!            Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,vx5,vy5,vz5,norm5)
+!            Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
 ! !!! PFL 28 Feb 2008 Test to be sure that angle does not change by more than Pi
 !            !           Xprimitive_t(I)=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5,vx2,vy2,vz2,norm2)*Pi/180.
 !            DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
@@ -605,8 +596,8 @@
 !            Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2)
 !            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1)
 !            Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3)
-!            Call produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2,vx5,vy5,vz5,norm5)
-!            Call produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2,vx4,vy4,vz4,norm4)
+!            Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,vx5,vy5,vz5,norm5)
+!            Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
 ! !!! PFL 28 Feb 2008 Test to be sure that angle does not change by more than Pi
 !            !           Xprimitive_t(I)=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5,vx2,vy2,vz2,norm2)*Pi/180.
 !            DiheTmp=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5,vx2,vy2,vz2,norm2)
Index: src/Calc_zmat_frag.f90
===================================================================
--- src/Calc_zmat_frag.f90	(revision 1)
+++ src/Calc_zmat_frag.f90	(revision 2)
@@ -36,7 +36,6 @@
 
   IMPLICIT NONE
 
-  CHARACTER(5) ::  AtName
   integer(KINT), INTENT(IN) ::  na,atome(na)
   real(KREAL), INTENT(IN)   ::  x(Na),y(Na),z(Na),fact
   real(KREAL), INTENT(IN)   ::  r_cov(0:Max_Z)
@@ -48,7 +47,7 @@
 ! Frozen contains the indices of frozen atoms
 ! INTEGER(KINT) Frozen(*),NFroz
 ! Number of fragment, Index of the current fragment for loops
-  INTEGER(KINT) :: NbFrag,IdxFrag
+  INTEGER(KINT) :: NbFrag
 ! Fragment(I) contains the fragment index to which I belongs
 ! NbAtFrag(j) contains the number of atoms of fragment j
   INTEGER(KINT), ALLOCATABLE :: Fragment(:),NbAtFrag(:) !na
@@ -75,14 +74,14 @@
   real(KREAL) ::  vx4,vy4,vz4,norm4
   real(KREAL) ::  vx5,vy5,vz5,norm5
   real(KREAL) ::  val,val_d, d12, d13,d23,d
-  Logical PasFini,Debug, FirstAt, DebugGaussian
+  Logical Debug, FirstAt, DebugGaussian
   LOGICAL, ALLOCATABLE :: DejaFait(:), FCaf(:) !(na)
 !  Logical, ALLOCATABLE ::  FrozAt(:)    !T if this atom is frozen
   LOGICAL F1213, F1223,F1323
 
   INTEGER(KINT) :: I,J,n0,n1,n2,n3,n4,IAt,IL,JL,IFrag,ITmp, K, KMax
-  INTEGER(KINT) :: IMax, I3,I1,Ip, IFragFroz
-  INTEGER(KINT) :: I0, Izm, JAt,Izm0
+  INTEGER(KINT) :: I3, I1, Ip
+  INTEGER(KINT) :: I0, Izm, JAt
   INTEGER(KINT) :: OrderZmat
 
   REAL(KREAL) :: sDihe
@@ -392,8 +391,7 @@
     STOP
      END IF
 
-     CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
-          vx5,vy5,vz5,norm5)
+     CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,  vx5,vy5,vz5,norm5)
 
 
 ! We search for the fourth atom while checking that it is not aligned with 1st and 2nd atoms.
@@ -439,8 +437,7 @@
         ! les autres fragments ne seraient additonnés que si l'on ne trouve pas notre bonheur dans le premier bloc
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-        CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
-             vx4,vy4,vz4,norm4)
+        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)
         val_d=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
              vx2,vy2,vz2,norm2)
         sDihe=abs(sin(val_d*pi/180.d0))
@@ -463,8 +460,7 @@
 !           if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
 !              n1=JAt
 !              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-!              CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
-!                   vx4,vy4,vz4,norm4)
+!              CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
 !              val_d=angle_d(vx4,vy4,vz4,norm4, &
 !                   vx5,vy5,vz5,norm5, &
 !                   vx2,vy2,vz2,norm2)
@@ -484,9 +480,7 @@
 !              if (.NOT.DejaFait(JAt)) THEN
 !                 n1=JAt
 !                 CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-!                 CALL produit_vect(vx1,vy1,vz1,norm1, &
-!                      vx2,vy2,vz2,norm2, &
-!                      vx4,vy4,vz4,norm4)
+!                 CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
 !                 val_d=angle_d(vx4,vy4,vz4,norm4, &
 !                      vx5,vy5,vz5,norm5, &
 !                      vx2,vy2,vz2,norm2)
@@ -627,17 +621,13 @@
 !     n4=Liaisons(Iat,2)
 !     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
 !     CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
-!     CALL produit_vect(vx3,vy3,vz3,norm3, &
-!          vx2,vy2,vz2,norm2, &
-!          vx5,vy5,vz5,norm5)
+!     CALL produit_vect(vx3,vy3,vz3, vx2,vy2,vz2, vx5,vy5,vz5,norm5)
 !
 !     DO While ((ITmp.LE.FrozBlock(Iat,0)).AND.(sDihe.LE.0.09d0))
 !        ITmp=ITmp+1
 !        n1=FrozBlock(Iat,Itmp)
 !        CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-!        CALL produit_vect(vx1,vy1,vz1,norm1, &
-!             vx2,vy2,vz2,norm2, &
-!             vx4,vy4,vz4,norm4)
+!        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)
 !        val_d=angle_d(vx4,vy4,vz4,norm4,  &
 !             vx5,vy5,vz5,norm5,           &
 !             vx2,vy2,vz2,norm2)
@@ -655,9 +645,7 @@
 !           if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
 !              n1=JAt
 !              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-!              CALL produit_vect(vx1,vy1,vz1,norm1,  &
-!                   vx2,vy2,vz2,norm2,               &
-!                   vx4,vy4,vz4,norm4)
+!              CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)
 !              val_d=angle_d(vx4,vy4,vz4,norm4,    &
 !                   vx5,vy5,vz5,norm5,              &
 !                   vx2,vy2,vz2,norm2)
@@ -676,9 +664,7 @@
 !              if (.NOT.DejaFait(JAt)) THEN
 !                 n1=JAt
 !                 CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-!                 CALL produit_vect(vx1,vy1,vz1,norm1,   &
-!                      vx2,vy2,vz2,norm2,                 &
-!                      vx4,vy4,vz4,norm4)
+!                 CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)
 !                 val_d=angle_d(vx4,vy4,vz4,norm4,       &
 !                      vx5,vy5,vz5,norm5,                 &
 !                      vx2,vy2,vz2,norm2)
@@ -830,9 +816,7 @@
 !     n4=I1
 !     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
 !     CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
-!     CALL produit_vect(vx3,vy3,vz3,norm3,  &
-!          vx2,vy2,vz2,norm2,    &
-!          vx5,vy5,vz5,norm5)
+!     CALL produit_vect(vx3,vy3,vz3, vx2,vy2,vz2,vx5,vy5,vz5,norm5)
 !
 !     !     We will look at the other frozen atoms
 !     !     we might improve the search so as to take the atom closest to IAt
@@ -842,9 +826,7 @@
 !        if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
 !           n1=JAt
 !           CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-!           CALL produit_vect(vx1,vy1,vz1,norm1, &
-!                vx2,vy2,vz2,norm2, &
-!                vx4,vy4,vz4,norm4)
+!           CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)
 !           val_d=angle_d(vx4,vy4,vz4,norm4,  &
 !                vx5,vy5,vz5,norm5,   &
 !                vx2,vy2,vz2,norm2)
@@ -863,9 +845,7 @@
 !           if (.NOT.DejaFait(JAt)) THEN
 !              n1=JAt
 !              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
-!              CALL produit_vect(vx1,vy1,vz1,norm1,  &
-!                   vx2,vy2,vz2,norm2, &
-!                   vx4,vy4,vz4,norm4)
+!              CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)
 !              val_d=angle_d(vx4,vy4,vz4,norm4,   &
 !                   vx5,vy5,vz5,norm5,           &
 !                   vx2,vy2,vz2,norm2)
@@ -1175,10 +1155,8 @@
         val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
 
         CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
-        CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
-             vx4,vy4,vz4,norm4)
-        CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
-             vx5,vy5,vz5,norm5)
+        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
+        CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,vx5,vy5,vz5,norm5)
 
         val_d=angle_d(vx4,vy4,vz4,norm4, vx5,vy5,vz5,norm5, &
              vx2,vy2,vz2,norm2)
@@ -1238,7 +1216,7 @@
   INTEGER(KINT), INTENT(IN) :: ind_zmat(Na,5)
   real(KREAL), INTENT(IN)   :: val_zmat(Na,3)
 
-  character(6) :: at1,at2,at3,at4,at5,d,a,dh
+  character(6) :: at1, at2, at3, at4, d, a, dh
   character(SCHARS), ALLOCATABLE :: tab(:) ! 3*na
   character(LCHARS) :: ligne
   
Index: src/PrintGeomVasp.f90
===================================================================
--- src/PrintGeomVasp.f90	(revision 1)
+++ src/PrintGeomVasp.f90	(revision 2)
@@ -12,7 +12,7 @@
   LOGICAL, INTENT(IN) :: FrozAt(Nat)
   INTEGER(KINT), INTENT(IN) :: IoUnit
 
-  INTEGER(KINT) :: IGeom, Idx, Iat, I,J
+  INTEGER(KINT) :: Iat, I, J
   REAL(KREAL), ALLOCATABLE :: GeomTmpC2(:,:),GeomTmpC(:,:)
   REAL(KREAL) :: MRot(3,3)
 
Index: src/valid.f90
===================================================================
--- src/valid.f90	(revision 1)
+++ src/valid.f90	(revision 2)
@@ -26,8 +26,9 @@
       CHARACTER(132) :: LINE, ValTMP, LString
       LOGICAL       :: ValueTMP, Found
       LOGICAL,SAVE       :: Tall=.False., FValidFile, ValAll
-      INTEGER(4) :: ios, i,j,k
+      INTEGER(4) :: ios, i, j
       INTEGER(4), SAVE ::  NEntries
+      LOGICAL, SAVE :: Debug=.TRUE.
 
       IF (.not.TInit) THEN
          NEntries=0
@@ -39,12 +40,12 @@
 ! The file exists so we will read it !
             DO While (ios==0)
                READ(IOTMP,'(A)',ERR=999,END=999) LINE
-!               WRITE(*,*) 'LINE READ',trim(Line)
+               if (Debug) WRITE(*,*) 'LINE READ',trim(Line)
                if (ios==0) THEN
                   Line=AdjustL(Line)
-!               WRITE(*,*) 'LINE ADJ',trim(Line)
+                  if (debug) WRITE(*,*) 'LINE ADJ',trim(Line)
                   call Upcase(Line)
-!               WRITE(*,*) 'LINE UP',trim(Line)
+                  if (debug) WRITE(*,*) 'LINE UP',trim(Line)
                   if (Line(1:1)/='#') THEN
 ! This is not a comment, we go on
                      i=index(Line,' ')
@@ -52,8 +53,8 @@
                      ValTMP=Line(1:i-1)
                      j=len_trim(line)
                      Line=adjustL((Line(i+1:j)))
-!                     WRITE(*,*) 'VALID INIT ',NEntries,trim(ValTMP),
-!     &   ' - ',trim(Line)
+                     if (debug) WRITE(*,*) 'VALID INIT ',NEntries,trim(ValTMP), &
+        ' - ',trim(Line)
                      if (ValTMP(1:1)=='.') ValTmp=ValTMP(2:)
                      if (ValTMP(1:1)=='F') ValueTMP=.False.
                      if (ValTMP(1:1)=='T') ValueTMP=.True.
@@ -72,11 +73,13 @@
             END DO
  999     CONTINUE
          CLOSE(IOTMP)
-!         WRITE(*,*) "DBG Valid, " // Trim(DebugFile) // "read, with "
-!    &    ,NEntries
-!         DO I=1,NEntries
-!            WRITE(*,*) List(I),Value(I)
-!         END DO
+         if (debug) THEN
+            WRITE(*,*) "DBG Valid, " // Trim(DebugFile) // "read, with " &
+                 ,NEntries
+            DO I=1,NEntries
+               WRITE(*,*) List(I),Value(I)
+            END DO
+         END IF
        END IF
       END IF
 
@@ -111,6 +114,10 @@
             isValid=.FALSE.
          END IF
          NEntries=Nentries+1
+         If (NEntries>NMax) THEN
+            WRITE(*,*) 'Error: NEntries > NMax in Valid ! recompile'
+            STOP
+         END IF
          List(Nentries)=Trim(LString)
          Value(Nentries)=Valall
       END IF
Index: src/egrad_mopac.f90
===================================================================
--- src/egrad_mopac.f90	(revision 1)
+++ src/egrad_mopac.f90	(revision 2)
@@ -22,14 +22,10 @@
   logical           :: debug
   LOGICAL, SAVE :: first=.true.
 
-  REAL(KREAL), ALLOCATABLE :: GradTmp(:)
+  REAL(KREAL) :: Pi
 
-  REAL(KREAL) :: d, a_val, Pi
-
-  REAL(KREAL) :: coef,x
-  INTEGER(KINT) :: iat, jat, kat, i, j, k,n3at
+  INTEGER(KINT) :: iat, i, j, k, n3at
   INTEGER(KINT) :: ITmp1,ITmp2, Idx
-  INTEGER(KINT), PARAMETER :: IOLOG=12, IOCOM=11
 
   CHARACTER(132) :: FileIn, FileOut, FileChk
 
Index: src/Calc_baker.f90
===================================================================
--- src/Calc_baker.f90	(revision 1)
+++ src/Calc_baker.f90	(revision 2)
@@ -40,21 +40,19 @@
   REAL(KREAL), ALLOCATABLE :: UMat_tmp(:,:)
   INTEGER(KINT) :: NbBonds, NbAngles, NbDihedrals, IGeom
 
-  real(KREAL) :: vx,vy,vz,dist, Norm
+  real(KREAL) :: vx, vy, vz, Norm
   real(KREAL) :: vx1,vy1,vz1,norm1
   real(KREAL) :: vx2,vy2,vz2,norm2
   real(KREAL) :: vx3,vy3,vz3,norm3
   real(KREAL) :: vx4,vy4,vz4,norm4
   real(KREAL) :: vx5,vy5,vz5,norm5
-  real(KREAL) :: val,val_d, Q, T
 
-  INTEGER(KINT) :: I,J, n1,n2,n3,n4,IAt,IL,JL,IFrag,ITmp, K, KMax
-  INTEGER(KINT) :: I0, IOld, IAtTmp, Izm, JAt, Kat, Lat, L, NOUT
+  INTEGER(KINT) :: I, J, IAt, K
+  INTEGER(KINT) :: Kat, Lat, L
 
   REAL(KREAL) :: sAngleIatIKat, sAngleIIatLat
-  REAL(KREAL) :: DiheTmp
 
-  LOGICAL :: debug, bond, AddPrimitiveCoord, FAIL
+  LOGICAL :: debug, bond, AddPrimitiveCoord
 
   INTERFACE
      function valid(string) result (isValid)
@@ -124,10 +122,7 @@
 !!!!!!!! PFL Debug
 !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-  REAL(KREAL), ALLOCATABLE :: V(:,:) ! (NCoord,NCoord)
-  REAL(KREAL) :: tol
   INTEGER(KINT) :: I1,I2,I3,I4
-  LOGICAL ::  FALOBBT,FALOBPrimt,FAloBBTInv
   LOGICAL :: debugPFL
   REAL(KREAL), ALLOCATABLE :: val_zmat(:,:)
   INTEGER(KINT), ALLOCATABLE :: indzmat0(:,:)
@@ -283,9 +278,9 @@
               Call vecteur(I2,I1,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1)
               Call vecteur(I2,I3,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2)
               Call vecteur(I3,I4,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3)
-                       CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
+                       CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
                             vx4,vy4,vz4,norm4)
-                       CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
+                       CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
                             vx5,vy5,vz5,norm5)
 !                       CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
 !                            vx2,vy2,vz2,norm2)*Pi/180.
@@ -519,9 +514,9 @@
 
                        ! IF (AddPrimitiveCoord) THEN
                        !   NbDihedrals=NbDihedrals+1
-                       !  CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
+                       !  CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
                             !      vx5,vy5,vz5,norm5)
-                       !CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
+                       !CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
                             !    vx4,vy4,vz4,norm4)
                        !         CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
                        !             vx2,vy2,vz2,norm2)*Pi/180.
@@ -598,9 +593,9 @@
                     IF (AddPrimitiveCoord) THEN
                        NbDihedrals=NbDihedrals+1
                        Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3)
-                       CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
+                       CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
                             vx5,vy5,vz5,norm5)
-                       CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
+                       CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
                             vx4,vy4,vz4,norm4)
                        CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
                             vx2,vy2,vz2,norm2)*Pi/180.
@@ -702,9 +697,9 @@
 !            Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1)
 !            Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2)
 !            Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3)
-!            Call produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
+!            Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
 !                 vx4,vy4,vz4,norm4)
-!            Call produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
+!            Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
 !                 vx5,vy5,vz5,norm5)
 
 !                  DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
Index: src/Space_GEDIIS_all.f90
===================================================================
--- src/Space_GEDIIS_all.f90	(revision 1)
+++ src/Space_GEDIIS_all.f90	(revision 2)
@@ -17,7 +17,7 @@
       REAL(KREAL) :: GeomSet(NGeomF,MRESET*NCoord), GradSet(NGeomF,MRESET*NCoord), ESET(NGeomF,MRESET)
       ! Geom and GRAD for all iterations.
 
-      INTEGER(KINT) :: I,J,K,NMK,MI,NI
+      INTEGER(KINT) :: I, K, NMK, MI, NI
       INTEGER(KINT), SAVE :: NRESET
 
       ! UPDATE PARAMETER AND GRADIENT SUBSPACE:
Index: src/Check_step.f90
===================================================================
--- src/Check_step.f90	(revision 1)
+++ src/Check_step.f90	(revision 2)
@@ -2,7 +2,7 @@
 ! v1.0
 ! Only valence angles are checked to remain in to 0:Pi range.
 
-SUBROUTINE Check_step(IGeom,Coord,Nat,NCoord,GeomOld,Step)
+SUBROUTINE Check_step(Coord,Nat,NCoord,GeomOld,Step)
 
   Use Path_module, only : pi,ncart
 
@@ -12,15 +12,15 @@
   INTEGER, PARAMETER :: KINT=KIND(1)
   INTEGER, PARAMETER :: KREAL=KIND(1.D0)
 
-  INTEGER(KINT), INTENT(IN) :: Nat,NCoord,IGeom
+  INTEGER(KINT), INTENT(IN) :: Nat,NCoord
   CHARACTER(32), INTENT(IN) :: Coord
   REAL(KREAL), INTENT(INOUT) :: Step(NCoord)
   REAL(KREAL), INTENT(IN) :: GeomOld(NCoord)
 
   ! As this subroutine is here for debugging, debug=T !
   LOGICAL, PARAMETER :: debug=.TRUE.
-  LOGICAL :: Fini
 
+
   ! Variables
   INTEGER(KINT) :: IBeg, Idx,Iat
   REAL(KREAL) ::  Fact
Index: src/ConvGrad_Cart2Int.f90
===================================================================
--- src/ConvGrad_Cart2Int.f90	(revision 1)
+++ src/ConvGrad_Cart2Int.f90	(revision 2)
@@ -36,7 +36,6 @@
       logical    debug, init
       REAL(KREAL), parameter :: zero = 0d0, tol = 1d-8,  eps = 1d-12
  
-      REAL(KREAL) , ALLOCATABLE :: AEigVec(:,:), AEigVal(:),ATA(:)
       INTEGER(KINT) :: n3, iat, k
       INTEGER(KINT) :: i
 !
@@ -45,7 +44,7 @@
 
 !     .. Local Arrays  for LAPACK
       INTEGER(KINT) :: Info, LWork, Rank
-      REAL(KREAL), ALLOCATABLE :: S(:), Work(:), GradCTmp(:,:)
+      REAL(KREAL), ALLOCATABLE :: Work(:), GradCTmp(:, :)
 
 !     IWORK dimension should be at least 3*MIN(M,N)*NLVL + 11*MIN(M,N),
 !     where NLVL = MAX( 0, INT( LOG_2( MIN(M,N)/(SMLSIZ+1) ) )+1 )
Index: src/minv.f90
===================================================================
--- src/minv.f90	(revision 1)
+++ src/minv.f90	(revision 2)
@@ -39,12 +39,12 @@
             IZ=N*(J-1)
             DO 20 I=K,N
                IJ=IZ+I
-   10          IF (ABS(BIGA).LT.ABS(A(IJ)))THEN
+               IF (ABS(BIGA).LT.ABS(A(IJ)))THEN
                   BIGA=A(IJ)
                   L(K)=I
                   M(K)=J
                ENDIF
-   20    CONTINUE
+20          CONTINUE
 !
 !     INTERCHANGE ROWS
 !
Index: src/egrad_vasp.f90
===================================================================
--- src/egrad_vasp.f90	(revision 1)
+++ src/egrad_vasp.f90	(revision 2)
@@ -22,22 +22,18 @@
 
   character(LCHARS) :: LINE
 
-  REAL(KREAL) :: d, a_val, Pi
+  REAL(KREAL) :: Pi
 
-  INTEGER(KINT) :: iat, jat, kat, i, j, n3at, Nfr
+  INTEGER(KINT) :: iat, i, n3at
   REAL(KREAL) :: RTmp1, RTmp2, RTmp3
 
   !
   CHARACTER(132), PARAMETER :: FileIn='POSCAR', FileOut='OUTCAR'
 
   CHARACTER(VLCHARS), SAVE :: RunCommand
-  INTEGER(kint) :: NbLists,LastIt,Nat4,NStep, Firstit
 
-  INTEGER(KINT) :: IGeom, Idx
 !  REAL(KREAL), ALLOCATABLE :: GeomTmpC(:), GeomTmp(:)
-  REAL(KREAL) :: MRot(3,3), Rmsd
 
-  REAL(KREAL) :: X,Y,Z
   LOGICAL :: Debug
 
   INTERFACE
Index: src/ConvXyzZmat.f90
===================================================================
--- src/ConvXyzZmat.f90	(revision 1)
+++ src/ConvXyzZmat.f90	(revision 2)
@@ -4,7 +4,7 @@
 
   IMPLICIT NONE
 
-  integer(KINT) :: i,i1,i2
+  integer(KINT) :: i
 
   integer(KINT) :: n1,n2,n3,n4,nb
   real(KREAL) ::  x(Nb),y(Nb),z(Nb)
@@ -106,18 +106,14 @@
 
               !        write(*,*) 'Debut procedure Xyz_Zmat43'
               CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
-              CALL produit_vect(vx1,vy1,vz1,norm1, &
-                   vx2,vy2,vz2,norm2, &
-                   vx4,vy4,vz4,norm4)
-              CALL produit_vect(vx3,vy3,vz3,norm3, &
-                   vx2,vy2,vz2,norm2, &
-                   vx5,vy5,vz5,norm5)
+              CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
+              CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, vx5,vy5,vz5,norm5)
               val_d=angle_d(vx4,vy4,vz4,norm4, &
                    vx5,vy5,vz5,norm5, &
                    vx2,vy2,vz2,norm2)
 
               !               write(*,11) n1,n2,norm1,n3,val,n4,val_d
-11            format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))
+!11            format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))
 
               val_zmat(i,1)=norm1
               val_zmat(i,2)=val
Index: src/Energy_GEDIIS.f90
===================================================================
--- src/Energy_GEDIIS.f90	(revision 1)
+++ src/Energy_GEDIIS.f90	(revision 2)
@@ -8,7 +8,7 @@
     REAL(KREAL) :: ci(MRESET), GeomSet(MRESET*NCoord), GradSet(MRESET*NCoord), ESET(MRESET)
     REAL(KREAL) :: ER_star, tmp
     
-    INTEGER(KINT) :: I, J, IX, JX, INV, JNV, KX
+    INTEGER(KINT) :: I, IX, JX, INV, JNV, KX
     
     LOGICAL PRINT
     
Index: src/Space_GEDIIS.f90
===================================================================
--- src/Space_GEDIIS.f90	(revision 1)
+++ src/Space_GEDIIS.f90	(revision 2)
@@ -18,7 +18,7 @@
       REAL(KREAL) :: GeomSet(MRESET*NCoord), GradSet(MRESET*NCoord), ESET(MRESET)
       ! Geom and GRAD for all iterations.
 
-      INTEGER(KINT) :: I,J,K,NMK,MI,NI
+      INTEGER(KINT) :: I, K, NMK, MI, NI
       INTEGER(KINT), SAVE :: NRESET
 
       ! UPDATE PARAMETER AND GRADIENT SUBSPACE:
Index: src/CalcBMat_int.f90
===================================================================
--- src/CalcBMat_int.f90	(revision 1)
+++ src/CalcBMat_int.f90	(revision 2)
@@ -1,13 +1,14 @@
-subroutine CalcBmat_int(natoms, xyzatm, indzmat, dzdc, atmass,atome)
+subroutine CalcBmat_int(natoms, xyzatm, indzmat, dzdc)
 
   use Io_module
-  use Path_module, only : MW,Pi
+  use Path_module, only : MW,Pi, atome
 
   IMPLICIT NONE
 
 
-  integer(KINT) :: natoms, indzmat(natoms,5),atome(natoms)
-  real(KREAL)   :: xyzatm(3,natoms), dzdc(3*natoms,3*natoms),atmass(natoms)
+  integer(KINT), INTENT(IN) :: natoms, indzmat(natoms,5)
+  real(KREAL), INTENT(IN)   :: xyzatm(3,natoms)
+  real(KREAL), INTENT(OUT)   :: dzdc(3*natoms,3*natoms)
 
 
 !!!!!!!!!!!!!!!!!!!!!
Index: src/Egrad.f90
===================================================================
--- src/Egrad.f90	(revision 1)
+++ src/Egrad.f90	(revision 2)
@@ -45,10 +45,7 @@
   REAL(KREAL) ::  Pi
 
   INTEGER(KINT) :: iat, i, j, IBeg,Idx
-  INTEGER(KINT), PARAMETER :: IOLOG=12, IOCOM=11
 
-  INTEGER(KINT), PARAMETER :: NbExtName=4
-
   REAL(KREAL), ALLOCATABLE :: XPrimRef(:), XPrim(:) ! NPrim
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Index: src/CheckPeriodicBound.f90
===================================================================
--- src/CheckPeriodicBound.f90	(revision 1)
+++ src/CheckPeriodicBound.f90	(revision 2)
@@ -13,7 +13,6 @@
 
   INTEGER(KINT) :: IGeom, I,J,K
   REAL(KREAL), ALLOCATABLE :: GeomRef(:,:), GeomTmp(:,:) ! (3,Nat)
-  CHARACTER(SCHARS) :: Line,Line2
 
   REAL(KREAL) :: Latrloc(3,3),Bloc(3),Xt,Yt,Zt
   LOGICAL :: Debug
Index: src/Hupdate_MS.f90
===================================================================
--- src/Hupdate_MS.f90	(revision 1)
+++ src/Hupdate_MS.f90	(revision 2)
@@ -40,12 +40,11 @@
   REAL(KREAL), INTENT(IN) :: dGrad(Nfree)
 
 
-  REAL(KREAL), ALLOCATABLE  :: HessOld(:,:) , HupMS(:,:)
+  REAL(KREAL), ALLOCATABLE  :: HessOld(:, :) 
   REAL(KREAL), ALLOCATABLE  :: Hup(:,:)
   REAL(KREAL), ALLOCATABLE  :: HupInt(:)
   REAL(KREAL), ALLOCATABLE  :: HS(:), Num(:,:)
-  REAL(KREAL) :: Phi, Num3, Den, Den2,NormDS
-  CHARACTER(120) :: fmt,fmt2,Line
+  REAL(KREAL) :: Den, NormDS
   INTEGER(KINT) :: i,j, idx
 
   LOGICAL :: debug
Index: src/Extrapol_baker.f90
===================================================================
--- src/Extrapol_baker.f90	(revision 1)
+++ src/Extrapol_baker.f90	(revision 2)
@@ -29,9 +29,8 @@
   ! XGeomF(NGeomF): Final geometries.
   REAL(KREAL), INTENT(OUT) :: XGeomF(NGeomF)
 
-  INTEGER(KINT) :: IdxGeom, I, J, K, Idx, IdxAtom
+  INTEGER(KINT) :: IdxGeom, I, J, K
   REAL(KREAL) :: Rmsd,MRot(3,3), ds, u, v
-  REAL(KREAL) :: a_val, d
 
   REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:), DerInt(:) ! (Nat,3)
   REAL(KREAL), ALLOCATABLE :: Xyz_k(:,:) ! (Nat,3)
@@ -43,7 +42,6 @@
   LOGICAL, EXTERNAL :: valid
 
   INTEGER(KINT) :: NSpline
-  CHARACTER(LCHARS) :: FileSpline,TmpChar
 
 
   ! We will calculate the length of the path, in MW coordinates...
Index: src/Decomp_frag.f90
===================================================================
--- src/Decomp_frag.f90	(revision 1)
+++ src/Decomp_frag.f90	(revision 2)
@@ -11,7 +11,7 @@
   LOGICAL,       INTENT(IN) :: ListAt(na)
   INTEGER(KINT), INTENT(OUT) :: NbFrag,Fragment(na),NbAtFrag(na),FragAt(na,na)
 
-  INTEGER(KINT) :: i,j,k,iat,jat
+  INTEGER(KINT) :: i, j, iat, jat
   INTEGER(KINT) :: IdxFrag, IdxtoDo, IFrag,IaFaire
   INTEGER(KINT), ALLOCATABLE :: ToDo(:)
 
Index: src/Write_path.f90
===================================================================
--- src/Write_path.f90	(revision 1)
+++ src/Write_path.f90	(revision 2)
@@ -7,7 +7,7 @@
 
   INTEGER(KINT), INTENT(IN) :: Iopt
 
-  INTEGER(KINT) :: IGeom, Idx, Iat, I,J
+  INTEGER(KINT) :: IGeom, Iat, I
   REAL(KREAL), ALLOCATABLE :: GeomTmpC(:), GeomTmp(:),GeomTmpC2(:,:)
   CHARACTER(SCHARS) :: Line
   CHARACTER(LCHARS) ::Title
Index: src/ConvertZmat_cart_3.f90
===================================================================
--- src/ConvertZmat_cart_3.f90	(revision 1)
+++ src/ConvertZmat_cart_3.f90	(revision 2)
@@ -10,19 +10,15 @@
 
         IMPLICIT NONE
 
-        integer(KINT) :: iat,n1,n2,n3
+        integer(KINT) :: iat, n1, n2
         real(KREAL)    :: x(nat),y(nat),z(nat)
         real(KREAL) ::  d,a_val
         integer(KINT) :: ind_zmat(Nat,5)
 
         real(KREAL) ::  vx1,vy1,vz1,norm1
         real(KREAL) ::  vvx1,vvy1,vvz1,normv1
-        real(KREAL) ::  vx2,vy2,vz2,norm2
-        real(KREAL) ::  vvx2,vvy2,vvz2,normv2
-        real(KREAL) ::  vx3,vy3,vz3,norm3
-        real(KREAL) ::  vx4,vy4,vz4,norm4
+        real(KREAL) ::  vx4, vy4, vz4
         real(KREAL) ::  a11_z1,a12_z1
-        real(KREAL) ::  a11_z2,a12_z2
         real(KREAL) ::  a11_y,a12_y
 
 
@@ -90,3 +86,4 @@
 
         END
 
+
Index: src/Splin1D.f90
===================================================================
--- src/Splin1D.f90	(revision 1)
+++ src/Splin1D.f90	(revision 2)
@@ -96,7 +96,6 @@
          REAL(KREAL) :: x,y
 
 
-       INTEGER(KINT) :: ind
        LOGICAL debug
 
 ! Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a
@@ -147,7 +146,10 @@
 ! klo and khi now bracket the input value of x.
           h=xa(khi)-xa(klo)
 ! The xa's must be distinct.
-          if (h.eq.0.) pause 'bad xa input in splint'
+          if (h.eq.0.) THEN
+             WRITE(*,*) 'bad xa input in splint'
+             STOP
+          END IF
 ! Cubic spline polynomial is now evaluated.
            a=(xa(khi)-x)/h
            b=(x-xa(klo))/h
@@ -175,7 +177,6 @@
          REAL(KREAL) :: x,y
 
 
-       INTEGER(KINT) :: ind
        LOGICAL debug
 
 ! Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a
@@ -220,7 +221,10 @@
 ! klo and khi now bracket the input value of x.
           h=xa(khi)-xa(klo)
 ! The xa's must be distinct.
-          if (h.eq.0.) pause 'bad xa input in splint'
+          if (h.eq.0.) THEN
+             WRITE(*,*) 'bad xa input in splinder'
+             STOP
+          END IF
 ! Cubic spline polynomial is now evaluated.
            a=(xa(khi)-x)/h
            b=(x-xa(klo))/h
@@ -255,7 +259,6 @@
          REAL(KREAL) :: x,y, der
 
 
-       INTEGER(KINT) :: ind
        LOGICAL debug
 
        INTEGER(KINT) :: k,khi,klo
@@ -301,7 +304,11 @@
 ! klo and khi now bracket the input value of x.
           h=xa(khi)-xa(klo)
 ! The xa's must be distinct.
-          if (h.eq.0.) pause 'bad xa input in splint'
+          if (h.eq.0.) THEN
+             WRITE(*,*) 'bad xa input in splintDer'
+             STOP
+          END IF
+
 ! Cubic spline polynomial is now evaluated.
            a=(xa(khi)-x)/h
            b=(x-xa(klo))/h
@@ -318,9 +325,8 @@
       SUBROUTINE LinearInt(x,y,der,N,xa,ya)
 
 ! Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a
-! function (with the xai's in order), and given the array y2a(1:n),
-! which is the output from spline above, and given a value of x,
-! this routine returns a cubic-spline interpolated value y.
+! function (with the xai's in order),  and given a value of x,
+! this routine returns a linear interpolated value y and its derivative der
 
          Use Vartypes
 
@@ -334,7 +340,6 @@
          REAL(KREAL) :: x,y, der
 
 
-       INTEGER(KINT) :: ind
        LOGICAL debug
 
        INTEGER(KINT) :: k,khi,klo
@@ -380,7 +385,11 @@
 ! klo and khi now bracket the input value of x.
           h=xa(khi)-xa(klo)
 ! The xa's must be distinct.
-          if (h.eq.0.) pause 'bad xa input in splint'
+          if (h.eq.0.) THEN
+             WRITE(*,*) 'bad xa input in LinearIntDer'
+             STOP
+          END IF
+
 ! Linear int now evaluated.
            a=(xa(khi)-x)/h
            b=(x-xa(klo))/h
Index: src/Hinvup_BFGS.f90
===================================================================
--- src/Hinvup_BFGS.f90	(revision 1)
+++ src/Hinvup_BFGS.f90	(revision 2)
@@ -46,8 +46,7 @@
   REAL(KREAL), ALLOCATABLE  :: Hup(:,:),ssT(:,:)
   REAL(KREAL), ALLOCATABLE  :: HupInt(:)
   REAL(KREAL), ALLOCATABLE  :: Num(:,:), Num2(:,:)
-  REAL(KREAL) :: Phi, Num3, Den, Den2,NormDS
-  CHARACTER(120) :: fmt,fmt2,Line
+  REAL(KREAL) :: Den, NormDS
   INTEGER(KINT) :: i,j,k, Idx
 
   LOGICAL debug,valid
Index: src/Step_DIIS_all.f90
===================================================================
--- src/Step_DIIS_all.f90	(revision 1)
+++ src/Step_DIIS_all.f90	(revision 2)
@@ -18,7 +18,7 @@
     INTEGER(KINT) :: NGeomF,IGeom
     INTEGER(KINT), INTENT(IN) :: NCoord
     
-    REAL(KREAL) :: Geom_new(NCoord),Geom(NCoord),Grad(NCoord)
+    REAL(KREAL) :: Geom(NCoord), Grad(NCoord)
     REAL(KREAL) :: Hess(NCoord*NCoord),Step(NCoord)
     REAL(KREAL) :: HEAT,HP
     LOGICAL :: allocation_flag
@@ -67,7 +67,7 @@
       REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:,:),GradSet(:,:),ERR(:,:) ! MRESET*NCoord
       REAL(KREAL), SAVE :: ESET(MRESET)
     REAL(KREAL), ALLOCATABLE, SAVE :: DXTMP(:,:),GSAVE(:,:) !NCoord, why DXTMP has SAVE attribute??
-      REAL(KREAL) :: B(M2),BS(M2),BST(M2)
+      REAL(KREAL) :: B(M2), BS(M2)
       LOGICAL DEBUG, PRINT
       INTEGER(KINT), ALLOCATABLE, SAVE :: MSET(:)
     LOGICAL, ALLOCATABLE, SAVE :: FRST(:)
Index: src/Makefile
===================================================================
--- src/Makefile	(revision 1)
+++ src/Makefile	(revision 2)
@@ -5,7 +5,7 @@
 # You might also have to edit the location of
 # some libraries (like mkl for ifort)
 
-Machine=gfortran
+Machine=arqg
 
 ###########################################
 #                                         #
@@ -135,8 +135,8 @@
 ifeq ($(Machine),arqP)
 # Flags for arq Ifort
 COMP=ifort
-F90=${COMP}   -g -check bounds  -check format -check uninit  -traceback
-#F90=${COMP}   -g -check all  -traceback
+#F90=${COMP}   -g -check bounds  -check format -check uninit  -traceback
+F90=${COMP}   -g -check all  -traceback
 F77=${F90}
 #F90=${COMP}  
 LINK=${COMP} -lguide -lpthread   -L/usr/lib/ \
@@ -148,7 +148,7 @@
 
 ifeq ($(Machine),arqg)
 # Flags for arq GFortran
-COMP=gfortran -fbacktrace -fbounds-check -Wall
+COMP=gfortran -g -fbacktrace -fbounds-check -Wall -Wextra
 F90=${COMP}  
 F77=${F90}
 #F90=${COMP}  
Index: src/bib_oper.f90
===================================================================
--- src/bib_oper.f90	(revision 1)
+++ src/bib_oper.f90	(revision 2)
@@ -100,14 +100,14 @@
 !       produit vectoriel
 !================================================================
 
-       SUBROUTINE produit_vect(v1x,v1y,v1z,norm1, &
-                               v2x,v2y,v2z,norm2, &
+       SUBROUTINE produit_vect(v1x,v1y,v1z, &
+                               v2x,v2y,v2z, &
                                v3x,v3y,v3z,norm3)
 
         use Path_module, only :  Pi,KINT, KREAL
 
-        real(KREAL) ::  v1x,v1y,v1z,norm1 ! what do you do with norm1, norm2???
-        real(KREAL) ::  v2x,v2y,v2z,norm2
+        real(KREAL) ::  v1x,v1y,v1z ! what do you do with norm1, norm2???
+        real(KREAL) ::  v2x,v2y,v2z
         real(KREAL) ::  v3x,v3y,v3z,norm3
 
         v3x= v1y*v2z-v1z*v2y
Index: src/Step_DIIS.f90
===================================================================
--- src/Step_DIIS.f90	(revision 1)
+++ src/Step_DIIS.f90	(revision 2)
@@ -59,7 +59,7 @@
       REAL(KREAL), ALLOCATABLE, SAVE :: XSET(:),GSET(:),ERR(:) ! MRESET*NVAR
       REAL(KREAL) :: ESET(MRESET)
       REAL(KREAL), ALLOCATABLE, SAVE :: DX(:),GSAVE(:) !NVAR
-      REAL(KREAL) :: B(M2),BS(M2),BST(M2)
+      REAL(KREAL) :: B(M2), BS(M2)
       LOGICAL DEBUG, PRINT
       INTEGER(KINT), SAVE :: MSET
       INTEGER(KINT) :: NDIIS, MPLUS, INV, ITERA, MM
Index: src/cmshft.f90
===================================================================
--- src/cmshft.f90	(revision 1)
+++ src/cmshft.f90	(revision 2)
@@ -22,7 +22,7 @@
   INTEGER(KINT) :: IShift
 
   REAL(KREAL), parameter :: zero=0.d0,one=1.d0
-  INTEGER(KINT) :: i,j
+  INTEGER(KINT) :: j
   REAL(KREAL) :: TmassI, tmass
 
   if(natom.le.0) then
Index: src/ConvertZmat_cart.f90
===================================================================
--- src/ConvertZmat_cart.f90	(revision 1)
+++ src/ConvertZmat_cart.f90	(revision 2)
@@ -17,9 +17,8 @@
         real(KREAL) ::  vvx1,vvy1,vvz1,normv1
         real(KREAL) ::  vx2,vy2,vz2,norm2
         real(KREAL) ::  vvx2,vvy2,vvz2,normv2
-        real(KREAL) ::  vx3,vy3,vz3,norm3
-        real(KREAL) ::  vx4,vy4,vz4,norm4
-        real(KREAL) ::  d,a_val,a_dih,c1,c2
+        real(KREAL) ::  vx4, vy4, vz4
+        real(KREAL) ::  d, a_val, a_dih
         real(KREAL) ::  a11_z1,a12_z1
         real(KREAL) ::  a11_z2,a12_z2
         real(KREAL) ::  a11_y,a12_y
Index: src/Step_GEDIIS_All.f90
===================================================================
--- src/Step_GEDIIS_All.f90	(revision 1)
+++ src/Step_GEDIIS_All.f90	(revision 2)
@@ -7,7 +7,7 @@
 
       INTEGER(KINT) :: NGeomF,IGeom
       INTEGER(KINT), INTENT(IN) :: NCoord
-      REAL(KREAL) :: Geom(NCoord), Grad(NCoord), Hess(NCoord*NCoord), Step(NCoord), Geom_new(NCoord)
+      REAL(KREAL) :: Geom(NCoord), Grad(NCoord), Hess(NCoord*NCoord), Step(NCoord)
     REAL(KREAL) :: HEAT ! HEAT= Energy
     LOGICAL :: allocation_flag
     REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord)
@@ -22,10 +22,10 @@
       INTEGER(KINT), ALLOCATABLE, SAVE :: MSET(:) ! mth Iteration
     LOGICAL, ALLOCATABLE, SAVE :: FRST(:)
     REAL(KREAL) :: ci(MRESET), ci_tmp(MRESET) ! MRESET = maximum number of iterations.
-      INTEGER(KINT) :: NGEDIIS, MPLUS, INV, ITERA, MM, cis_zero, IXX, JXX, MSET_minus_cis_zero
-      INTEGER(KINT) :: I,J,K, JJ, KJ, JNV, II, IONE, IJ, INK,ITmp, IX, JX, KX, MSET_C_cis_zero
-    INTEGER(KINT) :: current_size_B_mat, MyPointer, Iat, Isch, NFree, Idx
-      REAL(KREAL) :: XMax, XNorm, S, DET, THRES, tmp, ER_star, ER_star_tmp, Norm
+      INTEGER(KINT) :: NGEDIIS, MPLUS, INV, ITERA, MM, cis_zero
+      INTEGER(KINT) :: I, J, K, JJ, JNV, II, IONE, IJ, IX, JX, KX
+    INTEGER(KINT) :: current_size_B_mat, MyPointer, Isch, NFree, Idx
+      REAL(KREAL) :: XMax, XNorm, DET, THRES, tmp, ER_star, ER_star_tmp, Norm
     REAL(KREAL), PARAMETER :: eps=1e-12
       REAL(KREAL), PARAMETER :: crit=1e-8
     REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord
Index: src/Header.f90
===================================================================
--- src/Header.f90	(revision 1)
+++ src/Header.f90	(revision 2)
@@ -5,7 +5,6 @@
 
          CHARACTER(*) :: String
          CHARACTER(70) :: Head1,Sep
-         INTEGER, PARAMETER :: I=70
          INTEGER :: LenS,Len1,Pos1
 
          Head1="====================================================================="
Index: src/CalcRmsd.f90
===================================================================
--- src/CalcRmsd.f90	(revision 1)
+++ src/CalcRmsd.f90	(revision 2)
@@ -60,11 +60,11 @@
       real(KREAL) ::  x0c1,y0c1,z0c1, xc2,yc2,zc2
       
 
-      INTEGER(KINT) :: i, j, nd, ia
+      INTEGER(KINT) :: i, j, ia
       real(KREAL) :: x_norm, y_norm, lambda
       real(KREAL) :: Rmatrix(3,3)
-      real(KREAL) :: S(4,4), dS(4,4), q(4)
-      real(KREAL) :: tmp(3),tmp1(4),tmp2(4),EigVec(4,4), EigVal(4)
+      real(KREAL) :: S(4, 4)
+      real(KREAL) :: EigVec(4, 4), EigVal(4)
 
       INTERFACE
        function valid(string) result (isValid)
Index: src/ZmatBuild.f90
===================================================================
--- src/ZmatBuild.f90	(revision 1)
+++ src/ZmatBuild.f90	(revision 2)
@@ -14,37 +14,53 @@
 ! after the atom was added to ind_zmat for the three first atoms,
 ! but before for all others. Now, it starts at zero and it
 ! is increased before the atom is added to ind_zmat
-!     IMPLICIT NONE
 
-      use Path_module, only : Nat, NMaxL, max_Z, KINT, KREAL
+      use Path_module, only : Nat, NMaxL, max_Z, KINT, KREAL,Pi
+      use Io_module
 
+      IMPLICIT NONE
 
-      INTEGER(KINT) :: Izm
-      integer(KINT) :: na,atome(NAt),at,ind_zmat(NAt,5),long
-      real(KREAL) ::  x(NAt),y(NAt),z(NAt),fact
+
+      INTEGER(KINT) :: Izm,I,IAf,IAt,IAt0,J,Jat
+      INTEGER(KINT) :: n1, NMaxCF, NMinCP, VCF
+      INTEGER(KINT) :: IAta,IAtd,IAtDh,IAti,IAtv,ICat
+      INTEGER(KINT) :: ICP0, IdAt0, IdAtl
+      INTEGER(KINT) :: Idx, IIat, IIatDh, IndAFaire, IndFin
+      INTEGER(KINT) :: IdxAt0
+      INTEGER(KINT) :: ITmp, IZm0, IZm1, IZm2, IZm3, IZm4, IzmTmp
+      integer(KINT) :: na, atome(NAt), ind_zmat(NAt, 5)
+      real(KREAL) ::  x(NAt), y(NAt), z(NAt)
       real(KREAL) ::  val_zmat(NAt,3)
 !     ListAt contains the indices of frozen atoms
-      LOGICAL ListAT(NAt),FIzm1,FFirst
+      LOGICAL ListAT(NAt), FFirst
       INTEGER(KINT) :: Natreat
 
-      real(KREAL) ::  vx,vy,vz,dist
-      INTEGER(KINT) :: LIAISONS(NAt,0:NMaxL),Nbli,Nblj
-      INTEGER(KINT) :: LiaisonsBis(NAt,0:NMaxL)
+      real(KREAL) ::  dist
+      INTEGER(KINT) :: LIAISONS(NAt, 0:NMaxL)
       INTEGER(KINT) :: CAFaire(NAt)
       INTEGER(KINT) :: LieS_CP(NAt,0:NMaxL),LieS_CF(NAt,0:NMaxL)
       INTEGER(KINT) :: AtCP0(NAt),NbAt0,NbAt0r
-      Integer(KINT) :: AtTypCycl(Max_Z), NbAt3(NAt,2),CyCleAt(NAt)
       Integer(KINT) :: NbCycle
       real(KREAL) ::  vx1,vy1,vz1,norm1
       real(KREAL) ::  vx2,vy2,vz2,norm2
-      real(KREAL) ::  vx3,vy3,vz3,norm3
-      real(KREAL) ::  vx4,vy4,vz4,norm4
-      real(KREAL) ::  vx5,vy5,vz5,norm5
-      real(KREAL) ::  val,val_d
-      Logical AtTerm,Debug,Done
+      real(KREAL) ::  val
+      Logical Debug, Done
       Logical DejaFait(NAt)
-      LOGICAL F1213, F1223,F1323,FTmp
 
+
+      INTERFACE 
+         FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
+           use Path_module, only :  Pi,KINT, KREAL
+           real(KREAL) ::  v1x,v1y,v1z,norm1
+           real(KREAL) ::  v2x,v2y,v2z,norm2
+           real(KREAL) ::  angle
+         END FUNCTION ANGLE
+      END INTERFACE
+
+
+
+
+
 !     As we may have to treat only partial molecules, it may happen
 !     that there are no central atoms... so we first check this by looking
 !     for the least number of CP atoms
@@ -809,7 +825,7 @@
             END IF
             IndAFaire=IndAFaire+1
          END DO
-12345    CONTINUE
+!12345    CONTINUE
       END DO
 
       WRITE(*,*) 'TOTOTOTOTOTOTOTOT'
Index: src/Path.f90
===================================================================
--- src/Path.f90	(revision 1)
+++ src/Path.f90	(revision 2)
@@ -1557,7 +1557,7 @@
   END IF
 
   IOpt=0
-  Call EgradPath(IOpt,Flag_Opt_Geom)
+  Call EgradPath(IOpt)
 
   if (debug)  WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt
   Call Write_path(IOpt)
@@ -2252,7 +2252,7 @@
 
         ! We have the new path, we calculate its energy and gradient
 
-        Call EgradPath(IOpt,Flag_Opt_Geom)
+        Call EgradPath(IOpt)
         !IF(IOPT .EQ. 10) Then
         !         Print *, 'Stopping at my checkpoint.'
         !           STOP !This is my temporary checkpoint.
@@ -2428,10 +2428,13 @@
      END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc)
   END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
 
-  IF (PROG=="GAUSSIAN") THEN
-     DEALLOCATE(Gauss_paste)
-     DEALLOCATE(Gauss_root,Gauss_comment,Gauss_end,current)
-  END IF
+!   IF (PROG=="GAUSSIAN") THEN
+!      DEALLOCATE(Gauss_paste)
+!      DEALLOCATE(Gauss_root)
+!      DEALLOCATE(Gauss_comment)
+!      DEALLOCATE(Gauss_end)
+!   END IF
+
   DEALLOCATE(XyzGeomF, IntCoordF)
   DEALLOCATE(XyzGeomI, IntCoordI)
   DEALLOCATE(XyzTangent,IntTangent)
@@ -2465,7 +2468,7 @@
         previous => current  ! make current point to new head
      END DO
 
-     DEALLOCATE(current)
+
   END IF
 
   DEALLOCATE(GeomTmp, Hess, GradTmp)
Index: src/egrad_ext.f90
===================================================================
--- src/egrad_ext.f90	(revision 1)
+++ src/egrad_ext.f90	(revision 2)
@@ -20,42 +20,21 @@
 
   ! ======================================================================
 
-  character(LCHARS) :: LINE
 
   logical           :: debug
-  LOGICAL :: fexist, FSim
-  LOGICAL, SAVE :: first=.true.
-  LOGICAL, ALLOCATABLE :: Done(:)
 
-  REAL(KREAL), SAVE :: Eold=1.e6
-  REAL(KREAL) :: d, a_val, Pi
+  REAL(KREAL) :: Pi
 
-  REAL(KREAL) :: coef,x
-  INTEGER(KINT) :: iat, jat, kat, i, j, n3at, absidg, absidg2
-  INTEGER(KINT) :: kTmp, Istart,ITmp1,ITmp2, Idx
-  INTEGER(KINT), PARAMETER :: IOLOG=12, IOCOM=11
+  INTEGER(KINT) :: iat, i, n3at
 
   !
   CHARACTER(132) :: FileIn,FileOut
 
-  CHARACTER(LCHARS) :: ListName, TitleTmp, CH32SVAR1
-  CHARACTER(VLCHARS), SAVE :: RstrtCopy, RunCommand
-  LOGICAL, SAVE :: FCopyRstrt=.False., FOrderChecked=.False.
-  LOGICAL   :: FRdyn,FStopNucl, FTmp, Tchk,Tchk1
-  INTEGER(kint) :: NbLists,LastIt,Nat4,NStep, Firstit
-  INTEGER(KINT) :: NStepAd,NStepTmp
-  INTEGER(KINT) :: IShowTmp
-  REAL(KREAL)   :: FricPsiMin,FricPsiT, FricNuclTmp
+  CHARACTER(VLCHARS), SAVE :: RunCommand
+!  LOGICAL, SAVE :: FCopyRstrt=.False., FOrderChecked=.False.
 
 
-  INTEGER(KINT), PARAMETER :: NbExtName=4
 
-  INTEGER(KINT) :: ICouc
-
-  INTEGER(KINT) :: ILine
-  INTEGER(KINT) :: NPulses, PulseLen, NStepWarm,NStepAv,NStepEq,NStepTot
-  REAL(KREAL) :: TempWarm,Mean,Slope,Dev,RTmp1,RTmp2
-
   ! ======================================================================
 
   LOGICAL, EXTERNAL :: valid
@@ -112,8 +91,8 @@
 
   RETURN
 
-999 CONTINUE
-  if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!'
-  STOP
+!999 CONTINUE
+!  if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!'
+!  STOP
   ! ======================================================================
 end subroutine egrad_ext
Index: src/Mixed2cart.f90
===================================================================
--- src/Mixed2cart.f90	(revision 1)
+++ src/Mixed2cart.f90	(revision 2)
@@ -12,7 +12,7 @@
   REAL(KREAL), INTENT(IN) :: IntCoord(*)
   REAL(KREAL), INTENT(OUT) :: XyzCoord(Nat,3)
 
-  INTEGER(KINT) :: I,J,K,L,Idx,IBeg,Iat
+  INTEGER(KINT) :: I, J, Idx, IBeg, Iat
   REAl(KREAL), ALLOCATABLE :: ValZmat(:,:),XyzTmp(:,:)
   REAL(KREAL) :: d, a_val
 
Index: src/egrad_turbomole.f90
===================================================================
--- src/egrad_turbomole.f90	(revision 1)
+++ src/egrad_turbomole.f90	(revision 2)
@@ -21,18 +21,15 @@
 
   logical           :: debug
 
-  REAL(KREAL), ALLOCATABLE :: GradTmp(:)
 
-  REAL(KREAL) :: d, a_val, Pi
+  REAL(KREAL) :: Pi
 
-  REAL(KREAL) :: coef,x
-  INTEGER(KINT) :: iat, jat, kat, i, j, k,n3at
-  INTEGER(KINT) :: ITmp1,ITmp2, Idx
-  INTEGER(KINT), PARAMETER :: IOLOG=12, IOCOM=11
+  INTEGER(KINT) :: iat, i, n3at
+  INTEGER(KINT) :: Idx
 
-  CHARACTER(132) :: FileIn, FileOut, FileChk
+  CHARACTER(132) :: FileIn, FileOut
 
-  CHARACTER(3) :: SepIn=' ', SepOut=' '
+!  CHARACTER(3) :: SepIn=' ', SepOut=' '
   CHARACTER(VLCHARS), SAVE ::  RunCommand
 
   ! ======================================================================
@@ -135,8 +132,8 @@
 
   RETURN
 
-999 CONTINUE
-  WRITE(*,*) 'We should not be here !!!!'
-  STOP
+!999 CONTINUE
+!  WRITE(*,*) 'We should not be here !!!!'
+!  STOP
   ! ======================================================================
 end subroutine egrad_TurboMole
Index: src/Tensor.f90
===================================================================
--- src/Tensor.f90	(revision 1)
+++ src/Tensor.f90	(revision 2)
@@ -15,13 +15,13 @@
 
   IMPLICIT NONE
 
-  REAL(KREAL), parameter :: zero=0.d0,one=1.d0
+  REAL(KREAL), parameter :: zero=0.d0
 
   INTEGER(KINT), INTENT(IN) :: Natom
   real(KREAL) :: rx(Natom),ry(Natom),rz(Natom),rmass(Natom)
   real(KREAL) :: a,b,c,Tinert(3,3),CoordG(3)
 
-  INTEGER(KINT) :: jover, I, J
+  INTEGER(KINT) :: jover, J
   REAL(KREAL) :: Txx, Tyy, Tzz, Txy, Txz, Tyz, Txxi
   REAL(KREAL) :: Ryb, Rzc, Rxa, Rmj2, ryb2, rzc2, rxa2, Rmj
 
@@ -48,8 +48,7 @@
      CoordG(1)=a
      CoordG(2)=b
      CoordG(3)=c
-1000 FORMAT(3F15.3)
-     !c       WRITE(IOOUT,1000) (CoordG(i),i=1,3)
+
      jover=mod(natom,2)
      if(jover.eq.0) then
         txx=zero
Index: src/CalcCnct.f90
===================================================================
--- src/CalcCnct.f90	(revision 1)
+++ src/CalcCnct.f90	(revision 2)
@@ -32,8 +32,7 @@
 
   use Path_module, only : NMaxL, max_Z,Nom,Prog, KINT, KREAL
 
-  character(2) :: ATOM
-  integer(KINT) :: na,atome(Na),at,ind_zmat(Na,5),long
+  integer(KINT) :: na, atome(Na)
   real(KREAL) ::  x(Na),y(Na),z(Na),fact
   real(KREAL) :: vx,vy,vz,dist
   REAL(KREAL) :: r_cov(0:Max_Z)
Index: src/egrad_gaussian.f90
===================================================================
--- src/egrad_gaussian.f90	(revision 1)
+++ src/egrad_gaussian.f90	(revision 2)
@@ -20,41 +20,20 @@
   character(LCHARS) :: LINE
 
   logical           :: debug, G03Ok
-  LOGICAL :: fexist, FSim
   LOGICAL, SAVE :: first=.true.
-  LOGICAL, ALLOCATABLE :: Done(:)
 
-  REAL(KREAL), SAVE :: Eold=1.e6
-  REAL(KREAL), ALLOCATABLE :: GradTmp(:)
+  REAL(KREAL) :: Pi
 
-  REAL(KREAL) :: d, a_val, Pi
-
-  REAL(KREAL) :: coef,x
   INTEGER(KINT) :: ItG03
-  INTEGER(KINT) :: iat, jat, kat, i, j, n3at, absidg, absidg2
-  INTEGER(KINT) :: kTmp, Istart,ITmp1,ITmp2, Idx
-  INTEGER(KINT), PARAMETER :: IOLOG=12, IOCOM=11
+  INTEGER(KINT) :: iat, i, n3at
+  INTEGER(KINT) :: ITmp1, ITmp2
 
   CHARACTER(132) :: FileIn, FileOut, FileChk
 
   CHARACTER(3) :: SepIn=' < ', SepOut=' > '
-  CHARACTER(LCHARS) :: ListName, TitleTmp, CH32SVAR1
-  CHARACTER(VLCHARS), SAVE :: RstrtCopy, RunCommand
-  LOGICAL, SAVE :: FCopyRstrt=.False., FOrderChecked=.False.
-  LOGICAL   :: FRdyn,FStopNucl, FTmp, Tchk,Tchk1
-  INTEGER(kint) :: NbLists,LastIt,Nat4,NStep, Firstit
-  INTEGER(KINT) :: NStepAd,NStepTmp
-  INTEGER(KINT) :: IShowTmp
-  REAL(KREAL)   :: FricPsiMin,FricPsiT, FricNuclTmp
+  CHARACTER(VLCHARS), SAVE :: RunCommand
+!  LOGICAL, SAVE :: FCopyRstrt=.False., FOrderChecked=.False.
 
-  INTEGER(KINT), PARAMETER :: NbExtName=4
-
-  INTEGER(KINT) :: ICouc
-
-  INTEGER(KINT) :: ILine
-  INTEGER(KINT) :: NPulses, PulseLen, NStepWarm,NStepAv,NStepEq,NStepTot
-  REAL(KREAL) :: TempWarm,Mean,Slope,Dev,RTmp1,RTmp2
-
   ! ======================================================================
 
   LOGICAL, EXTERNAL :: valid
@@ -224,8 +203,8 @@
 
   RETURN
 
-999 CONTINUE
-  if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!'
-  STOP
+!999 CONTINUE
+!  if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!'
+!  STOP
   ! ======================================================================
 end subroutine egrad_gaussian
Index: utils/AnaPath
===================================================================
--- utils/AnaPath	(revision 1)
+++ utils/AnaPath	(revision 2)
@@ -6,9 +6,9 @@
 
 Fout=$1
 
-ItMax=`grep MAXCYC $Fout | tail -1 | awk '{print $3}'`
-Nom=`grep PATHNAME $Fout | tail -1 | awk '{print $3}'`
-NGeomF=`grep NGEOMF $Fout | tail -1 | awk '{print $3}'`
+ItMax=`grep MAXCYC $Fout | tail -1 | awk '{print $NF*1}'`
+Nom=`grep PATHNAME  $Fout | tail -1 | awk -F "=" '{print $2}' | sed 's/"//g; s/,//g; s/ //g'`
+NGeomF=`grep NGEOMF $Fout | tail -1 | awk '{print $NF*1}'`
 
 xcol=1
 if [ $# -eq 2 ]; then
Index: utils/AnaPathrel
===================================================================
--- utils/AnaPathrel	(revision 1)
+++ utils/AnaPathrel	(revision 2)
@@ -5,11 +5,12 @@
 fi
 
 Fout=$1
-ItMax=`grep MAXCYC $Fout | tail -1 | awk '{print $3}'`
-Nom=`grep PATHNAME $Fout | tail -1 | awk '{print $3}'`
-NGeomF=`grep NGEOMF $Fout | tail -1 | awk '{print $3}'`
 
+ItMax=`grep MAXCYC $Fout | tail -1 | awk '{print $NF*1}'`
+Nom=`grep PATHNAME  $Fout | tail -1 | awk -F "=" '{print $2}' | sed 's/"//g; s/,//g; s/ //g'`
+NGeomF=`grep NGEOMF $Fout | tail -1 | awk '{print $NF*1}'`
 
+
 export LANG=C
 
 echo "#ItMax=$ItMax"
