From 68257ae8a7ffdc9c33b069448599f22d8f7d30ed Mon Sep 17 00:00:00 2001
From: Huub Van Dam <hvandam@bnl.gov>
Date: Wed, 2 Apr 2014 06:26:22 +0000
Subject: [PATCH] HvD: While checking the dispersion corrections in the DFT I
 noticed that the bond lengths in H2O in different orientations agreed only in
 8 digits. This seemed like there was a REAL*8 -> REAL -> REAL*8 conversion
 going on somewhere some place to do with geometries. In response I have:

- eliminated all implicit typing
- replaced all hardwired definitions of Pi with computing Pi to machine
  precision when needed
- ensured that all definitions of Bohr are the same and are the one that
  lives in util_params.fh

Unfortunately this did not fix the discrepancy I was seeing but it cleans
things up a bit anyway.
---
 src/geom/geom.F            |  33 ++-
 src/geom/geom_hnd.F        | 405 +++++++++++++++++++++++++------------
 src/geom/geom_input.F      | 249 +++++++++++++++--------
 src/geom/geom_print_ecce.F |   1 +
 4 files changed, 476 insertions(+), 212 deletions(-)

diff --git a/src/geom/geom.F b/src/geom/geom.F
index 31c3cf902d..6ef2e8a3ef 100644
--- a/src/geom/geom.F
+++ b/src/geom/geom.F
@@ -2885,7 +2885,7 @@ c
 c
       if (.not. geom_tag_to_element(tag, symbol, element, atn)) then
 c
-c     Is not an atom.  Try removing Bq or X.
+c        Is not an atom.  Try removing Bq or X.
 c
          if (inp_compare(.false., tag(1:1), 'x')) then
             ttag = tag(2:)
@@ -2894,8 +2894,25 @@ c
          else
             return              ! Nothing recognizable
          endif
-         if (.not. geom_tag_to_element(ttag, symbol, element, atn))
-     $        return
+         if (.not. geom_tag_to_element(ttag, symbol, element, atn)) then
+c
+c           We found a "Bq" or "X" but it is not labeled with an element
+c           (e.g. "XH" or "BqN") so we cannot associate any atomic
+c           properties. 
+c
+            geom_tag_to_covalent_radius = .false.
+            return
+         else
+c
+c           We found a "Bq" or "X" with an atomic type indication and 
+c           we now have the atomic number indicated, so we will go 
+c           with that.
+c
+         endif
+      else
+c
+c        We found an atom and now know its atomic number
+c
       endif
 c
 c     atn should be set to something sensible
@@ -3598,6 +3615,7 @@ c
 C> \return Return .true. if successfull, and .false. otherwise.
 c
       logical function geom_nuc_charge(geom, total_charge)
+      implicit none
 #include "nwc_const.fh"
 #include "geomP.fh"
       integer geom                  !< [Input] the geometry handle
@@ -3729,6 +3747,7 @@ C> \return Return .true. if the system type was found successfully, and
 C> .false. otherwise.
 c
       logical function geom_systype_get(geom, itype)
+      implicit none
 #include "nwc_const.fh"
 #include "geomP.fh"
       integer geom              !< [Input] the geometry handle
@@ -3865,7 +3884,7 @@ c
 c     
       integer geom,i,j               ! [input]
       double precision lattice(6)    ! [output]
-      real*8 rad,scale
+      double precision rad,scale
       logical  geom_check_handle,geom_get_user_scale
       external geom_check_handle,geom_get_user_scale
 
@@ -4140,11 +4159,11 @@ c
       end
       subroutine xlattice_abc_abg(a,b,c,alpha,beta,gamma,lattice_unita)
       implicit none
-      real*8 a,b,c
-      real*8 alpha,beta,gamma,lattice_unita(3,3)
+      double precision a,b,c
+      double precision alpha,beta,gamma,lattice_unita(3,3)
 
 *     *** local variables ****
-      real*8 d2,pi
+      double precision d2,pi
 
 *     **** determine a,b,c,alpha,beta,gmma ***
       pi = 4.0d0*datan(1.0d0)
diff --git a/src/geom/geom_hnd.F b/src/geom/geom_hnd.F
index 3e9d283200..2f06057cd6 100644
--- a/src/geom/geom_hnd.F
+++ b/src/geom/geom_hnd.F
@@ -101,31 +101,37 @@ c
 c     
       end
       SUBROUTINE HND_BMAT(GEOM,NZVAR,NCART,B,ZMAT,IZMAT,NZMAT,C)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "errquit.fh"
 #include "global.fh"
 #include "mafdecls.fh"
 #include "stdio.fh"
 #include "geom.fh"
+#include "util_params.fh"
       integer geom
 C     
 C     ----- CONSTRUCT THE B MATRIX -----
 C     
+      integer na,nzvar,ncart,nzmat
       PARAMETER (    NA=10 )
       LOGICAL     DBUG
       LOGICAL     OUT
       LOGICAL     GEOM_ZMT_GET_IZMAT
       LOGICAL     GEOM_ZMT_GET_NIZMAT
       CHARACTER*8 ERRMSG
-      DIMENSION B(NCART,*)
-      DIMENSION ZMAT(*),IZMAT(*)
-      DIMENSION A(3,NA)
-      DIMENSION C(3,*)
+      double precision B(NCART,*)
+      double precision ZMAT(*)
+      integer          IZMAT(*)
+      double precision A(3,NA)
+      double precision C(3,*)
       DIMENSION ERRMSG(3)
+      double precision eqval,pi2,degree,bohr,zero,four,pi,pio2,rtod
+      integer iw,iadd,i,itype,izvar,j,ndum
       DATA ERRMSG /'PROGRAM ','STOP IN ','- BMAT -'/
-      DATA PI2    /6.28318530717958D+00/
+c     DATA PI2    /6.28318530717958D+00/
       DATA DEGREE /360.0D+00/
-      DATA BOHR   /0.52917715D+00/
+      DATA BOHR   /cau2ang/
       DATA ZERO   /0.0D+00/
       DATA FOUR   /4.0D+00/
 C     
@@ -134,6 +140,8 @@ C
       OUT   =.FALSE. 
       OUT   =OUT.OR.DBUG
 C     
+      pi    =acos(-1.0d0)
+      pi2   =2.0d0*pi
       PIO2  =PI2/FOUR
       RTOD  =PI2/DEGREE
 C     
@@ -472,13 +480,19 @@ c
 c
       END
       SUBROUTINE HND_BSTR(EQVAL,NOINT,I,J,C,B,NCART,BOHR)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 C
 C     -----THIS ROUTINE COMPUTES THE B MATRIX ELEMENTS FOR A
 C          BOND STRETCH AS DEFINED BY WILSON (SEE WDC P.55) -----
 C
-      DIMENSION C(3,1),B(NCART,1)
-      DIMENSION RIJ(3)
+      double precision eqval,bohr
+      integer          noint,i,j,ncart
+      double precision C(3,1),B(NCART,1)
+      integer          nocol1,nocol2,m
+      double precision dijsq
+      double precision RIJ(3)
+      double precision ZERO
       DATA ZERO /0.0D+00/
       DIJSQ = ZERO
       DO 100 M = 1,3
@@ -493,7 +507,8 @@ C
       RETURN
       END
       SUBROUTINE HND_BEND(EQVAL,NOINT,I,J,K,C,B,NCART,RTOD,BOHR)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "errquit.fh"
 C
 C     -----THIS ROUTINE COMPUTES THE B MATRIX ELEMENTS OF A
@@ -503,8 +518,15 @@ C
 C     -----I AND K ARE THE NUMBERS OF THE END ATOMS.  J IS THE
 C          NUMBER OF THE CENTRAL ATOM -----
 C
-      DIMENSION C(3,1),B(NCART,1)
-      DIMENSION RJI(3),RJK(3),EJI(3),EJK(3)
+      double precision eqval,rtod,bohr
+      integer          noint,i,j,k,ncart
+      integer          m,nocol1,nocol2,nocol3
+      double precision C(3,1),B(NCART,1)
+      double precision RJI(3),RJK(3),EJI(3),EJK(3)
+      double precision zero,one,tol,pideg
+      double precision djisq,djksq
+      double precision dotj,dji,djk,dot
+      double precision sinj
       DATA ZERO,ONE  /0.0D+00,1.0D+00/
       DATA TOL,PIDEG /5.0D-05,180.00D+00/
 C
@@ -547,7 +569,8 @@ C
       RETURN
       END
       SUBROUTINE HND_TORS(EQVAL,NOINT,I,J,K,L,C,B,NCART,RTOD,BOHR)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 C
 C     ----- THIS ROUTINE COMPUTES THE B MATRIX ELEMENTS FOR THE
 C           TORSION AS DEFINED BY WILSON.  SEE WDC P60.
@@ -560,12 +583,19 @@ C           ALONG THE VECTOR J-->K WITH J NEARER THE OBSERVER, THEN
 C           THE CLOCKWISE ROTATION OF J-->I WHICH SUPERPOSES J--> I
 C           WITH K-->L IS A NEGATIVE TORSION ANGLE.
 C
+      double precision eqval,rtod,bohr
+      integer          noint,i,j,k,l,ncart
       CHARACTER*8 ERRMSG
-      DIMENSION C(3,1),B(NCART,1)
-      DIMENSION RIJ(3),RJK(3),RKL(3)
-      DIMENSION EIJ(3),EJK(3),EKL(3)
-      DIMENSION CR1(3),CR2(3),D(3)
+      double precision C(3,1),B(NCART,1)
+      double precision RIJ(3),RJK(3),RKL(3)
+      double precision EIJ(3),EJK(3),EKL(3)
+      double precision CR1(3),CR2(3),D(3)
+      double precision DKLSQ,DKL,DIJ,DIJSQ,DJKSQ,DJK,DOTPK,DOT
+      double precision DOTPJ,F2,F1,SML
+      double precision SINPJ,SMI,SINPK,SMJ
+      integer          iw,nocol1,nocol2,nocol3,nocol4,m
       DIMENSION ERRMSG(3)
+      double precision ZERO,ONE,TOLRD,TOL,PIDEG
       DATA ERRMSG    /'PROGRAM ','STOP IN ','-TORS  -'/
       DATA ZERO,ONE  /0.0D+00,1.0D+00/
       DATA TOLRD     /1.0001D+00/
@@ -694,7 +724,8 @@ cedo  CALL HND_HNDERR(3,ERRMSG)
      1       I3,' - ',I3,' - ',I3,' - ',I3,'. RUN STOPPED.')
       END
       SUBROUTINE HND_OPLA(EQVAL,NOINT,I,J,K,L,C,B,NCART,RTOD,PIO2,BOHR)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "errquit.fh"
 C
 C     -----THIS ROUTINE COMPUTES THE B MATRIX ELEMENTS FOR AN
@@ -703,10 +734,16 @@ C          SEE WDC P58. -----
 C     -----I IS THE END ATOM, J IS THE APEX ATOM, AND
 C          K AND L ARE THE ANCHOR ATOMS -----
 C
-      DIMENSION B(NCART,1),C(3,1)
-      DIMENSION RJI(3),RJK(3),RJL(3)
-      DIMENSION EJI(3),EJK(3),EJL(3)
-      DIMENSION C1(3),C2(3),C3(3)
+      double precision eqval,rtod,pio2,bohr
+      integer          noint,i,j,k,l,ncart
+      double precision B(NCART,1),C(3,1)
+      double precision RJI(3),RJK(3),RJL(3)
+      double precision EJI(3),EJK(3),EJL(3)
+      double precision C1(3),C2(3),C3(3)
+      double precision zero,one,tol,pideg
+      double precision DJISQ,DJKSQ,DJLSQ,DJI,DJK,DJL,DET,COST,DOTI,DOT
+      double precision SMK,SINI,SMI,SINT,TANT,SML,THETA
+      integer          m,nocol1,nocol2,nocol3,nocol4
       DATA ZERO,ONE  /0.0D+00,1.0D+00/
       DATA TOL,PIDEG /5.0D-05,180.00D+00/
 C
@@ -780,7 +817,8 @@ C
       RETURN
       END
       SUBROUTINE HND_LIBE(EQVAL,NOINT,I,J,K,L,C,B,NCART,RTOD,AA,NA,BOHR)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 C
 C     ----- THIS ROUTINE COMPUTES THE B MATRIX ELEMENTS FOR A
 C           PAIR OF PERPENDICULAR LINEAR BENDING COORDINATES.  SEE
@@ -797,17 +835,23 @@ C           PLANE AND THE SECOND IN A PLANE PERPENDICULAR TO THE FIRST
 C           THROUGH POINTS I,J, AND K -----
 C
 #include "stdio.fh"
+      double precision eqval,rtod,bohr
+      integer noint,i,j,k,l,ncart,na
       CHARACTER*8 ERRMSG
-      DIMENSION B(NCART,1),C(3,1)
-      DIMENSION AA(3,1),A(3)
-      DIMENSION RJI(3),RJK(3)
-      DIMENSION UNIT(3),UP(3),UN(3)
-      DIMENSION EJI(3),EJK(3)
+      double precision B(NCART,1),C(3,1)
+      double precision AA(3,1),A(3)
+      double precision RJI(3),RJK(3)
+      double precision UNIT(3),UP(3),UN(3)
+      double precision EJI(3),EJK(3)
       DIMENSION ERRMSG(3)
       DATA ERRMSG    /'PROGRAM ','STOP IN ','-LIBE  -'/
+      double precision zero,one,tol,pideg,pt0001
       DATA ZERO,ONE  /0.0D+00,1.0D+00/
       DATA TOL,PIDEG /5.0D-05,180.00D+00/
       DATA PT0001    /0.0001D+00/
+      integer iw,m,nocol1,nocol2,nocol3,no2
+      double precision DJISQ,DJKSQ,DAJSQ,DJI,DJK,DAJ,DOTJ,DOTP,TEST
+      double precision DUM
 c
       IW = LuOut
 C
@@ -899,15 +943,23 @@ C
       END
       SUBROUTINE HND_DIHPLA(DIHANG,NOINT,IZMAT,CARTC,BMAT,NCART,
      1                      DTORAD)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "stdio.fh"
+      double precision dihang,dtorad
+      integer noint,ncart
       CHARACTER*8 ERRMSG
-      DIMENSION IZMAT(5), CARTC(3,1), BMAT(NCART,1),
+      integer          IZMAT(5)
+      double precision CARTC(3,1), BMAT(NCART,1),
      1          A(3), B(3), C(3), D(3), E1(3), E2(3), E3(3)
       DIMENSION ERRMSG(3)
       DATA ERRMSG   /'PROGRAM ','STOP IN ','-DIHPLA-'/
+      double precision zero,one,tol
       DATA ZERO,ONE /0.0D+00,1.0D+00/
       DATA TOL      /1.0D-06/
+      integer          i,j,k,l,m,n,iw,iatom,jatom,katom,latom,matom,ixyz
+      double precision b1,b2,b3,b4,b5,ADOTE2,F1,F2,F4,F5,E1DE2,CDOTE1
+      double precision E2MAG,E1MAGI,E1MAG,E2MAGI,SINDI
 C
 C
 C     COMPUTE THE B MATRIX AND DIHEDRAL ANGLE BETWEEN 5 ATOMS
@@ -1058,7 +1110,8 @@ C
      1          'INTERNAL COORDINATE',I4,' USES ATOMS',5I4)
       END
       logical function geom_zmtmak(rtdb, geom, oprint)
-      implicit double precision (a-h,o-z)
+c     implicit double precision (a-h,o-z)
+      implicit none
 #include "errquit.fh"
 #include "stdio.fh"
 #include "global.fh"
@@ -1067,8 +1120,11 @@ C
 #include "util.fh"
 #include "nwc_const.fh"
 #include "rtdb.fh"
+#include "util_params.fh"
+      integer  rtdb
       integer  geom
       logical  oprint
+      integer  max_zcoord
       parameter (max_zcoord=160)
 c     
 c     Make a Z-matrix for the current geometry.  If successful store
@@ -1077,7 +1133,10 @@ c
 c     This version does not use the old HONDO common blocks so
 c     it may be called without preamble or fear of side-effects.
 c     
-      parameter (toangs=0.52917724924d+00)
+      double precision toangs
+      integer          mxatom,mxcart,mxzmat,mxcoor,mxizmt,mxbond,mxbnds
+      integer          mxangs,mxtors,mxoopa,mxlinb,mxlnba,mxseg
+      parameter (toangs=cau2ang)
       parameter (mxatom=nw_max_atom)
       parameter (mxcart=3*mxatom)
       parameter (mxzmat=nw_max_zmat)    
@@ -1154,17 +1213,25 @@ c
       dimension endatm(mxatom)
       dimension endmod(mxatom)
       dimension  idone(mxatom)
-      dimension iatseg(mxatom,mxseg)
-      dimension namseg(       mxseg)
-      dimension lenseg(       mxseg)
+      integer   iatseg(mxatom,mxseg)
+      integer   namseg(       mxseg)
+      integer   lenseg(       mxseg)
 c
       character*8 zvarname(mxcoor)
       double precision zvarsign(mxcoor)
 c
       double precision zero
       parameter (zero = 0.0d0)
+c
+      integer iat,jat,kat,nseg,iadd,i,j,k,l,lcon,ianchr,ii,iatom,jatom
+      integer ipass,itype,izvar,icon,jcon,kcon
+      integer m,mlnba,mbond,mangs,max_tor_per_bond,mbnds,mlinb,ncvrpass
+      integer moopa,mseg,mtors,mxconn
+      integer nbnds,npass,nlnba,numang,nuconn,numods,numbnd,numlnb
+      integer numijxtra,numtor,numoop,nuser
+      double precision dist,cvr_scaling,cvfac,rij,rcv,radius
 c     
-      dist(iat,jat)=sqrt((c(1,iat)-c(1,jat))**2+
+      dist(iat,jat)=dsqrt((c(1,iat)-c(1,jat))**2+
      1     (c(2,iat)-c(2,jat))**2+
      2     (c(3,iat)-c(3,jat))**2) * toangs
 c     
@@ -1949,9 +2016,15 @@ c
      3                      ijklop,numoop,
      4                      ijklnb,numlnb,xyzlnb,
      $     max_tor_per_bond, oprint)
-      implicit double precision (a-h,o-z)
+c     implicit double precision (a-h,o-z)
+      implicit none
 #include "stdio.fh"
 #include "nwc_const.fh"
+      integer nzvar,nizmat,mxizmt
+      integer nat,nbnds,numbnd,numang,numtor,numoop,numlnb
+      integer max_tor_per_bond
+      integer i
+      integer mxatom,mxbnds,mxangs,mxtors,mxoopa,mxlinb
       parameter (mxatom=nw_max_atom)        
       parameter (mxbnds= 8*mxatom)
       parameter (mxangs= 8*mxatom)
@@ -1963,13 +2036,14 @@ c
       logical     zdone, oprint
       logical     dolinb
       integer iw
-      dimension c(3,*)
-      dimension izmat(*)   
-      dimension ijbnds(2,*)
-      dimension iibnds(  *)
-      dimension nibnds(  *)
-      dimension ijkang(3,*),ijklto(4,*),ijklop(4,*)
-      dimension ijklnb(4,*),xyzlnb(3,*)
+      double precision c(3,*)
+      integer   izmat(*)   
+      integer   ijbnds(2,*)
+      integer   iibnds(  *)
+      integer   nibnds(  *)
+      integer   ijkang(3,*),ijklto(4,*),ijklop(4,*)
+      integer   ijklnb(4,*)
+      double precision xyzlnb(3,*)
 c
       iw=luout
       dbug=.false.
@@ -2121,22 +2195,25 @@ c
       end
       subroutine hnd_zmtyp1(nzvar,nizmat,izmat,mxizmt,
      1                      nat,nbnds,ijbnds,iibnds,nibnds,numbnd)
-      implicit double precision (a-h,o-z)
+c     implicit double precision (a-h,o-z)
+      implicit none
 #include "errquit.fh"
 #include "stdio.fh"
 #include "nwc_const.fh"
+      integer nzvar,nizmat,mxizmt,nat,nbnds,numbnd
+      integer ione,mxatom,mxbnds
       parameter (ione=1)
       parameter (mxatom=nw_max_atom)        
       parameter (mxbnds= 8*mxatom)
       logical     dbug
-      integer iw
-      dimension izmat(*)   
-      dimension ijbond(2,mxbnds)
-      dimension iibond(  mxatom)
-      dimension nibond(  mxatom)
-      dimension ijbnds(2,*)
-      dimension iibnds(  *)
-      dimension nibnds(  *)
+      integer iw,ij,iat,i,ii,ji,jj,ibnd,jbnd,ijbnd,nbnd,nbond
+      integer izmat(*)   
+      integer ijbond(2,mxbnds)
+      integer iibond(  mxatom)
+      integer nibond(  mxatom)
+      integer ijbnds(2,*)
+      integer iibnds(  *)
+      integer nibnds(  *)
 c
       iw=luout
       dbug=.false.
@@ -2204,10 +2281,13 @@ c
       subroutine hnd_zmtyp2(nzvar,nizmat,izmat,mxizmt,c,
      1     nat,nbnds,ijbnds,iibnds,nibnds,ijkang,numang,
      2     dolinb, zdone, oprint)
-      implicit double precision (a-h,o-z)
+c     implicit double precision (a-h,o-z)
+      implicit none
 #include "errquit.fh"
 #include "stdio.fh"
 #include "nwc_const.fh"
+      integer nzvar,nizmat,nat,nbnds,numang,mxizmt
+      integer itwo,mxatom,mxbnds,mxangs
       parameter (itwo=2)
       parameter (mxatom=nw_max_atom)        
       parameter (mxbnds= 8*mxatom)
@@ -2218,18 +2298,25 @@ c
       logical     linbnd
       logical     oprint
       integer iw
-      dimension c(3,*)
-      dimension izmat(*)   
-      dimension ijkang(3,*)     
-      dimension ijbond(2,mxbnds)
-      dimension iibond(  mxatom)
-      dimension nibond(  mxatom)
-      dimension ijbnds(2,*)
-      dimension iibnds(  *)
-      dimension nibnds(  *)
+      double precision c(3,*)
+      integer          izmat(*)   
+      integer          ijkang(3,*)     
+      integer          ijbond(2,mxbnds)
+      integer          iibond(  mxatom)
+      integer          nibond(  mxatom)
+      integer          ijbnds(2,*)
+      integer          iibnds(  *)
+      integer          nibnds(  *)
+      double precision tol,two
       data tol    /1.0d-01/     ! 5.7 degrees
 *     data tol    /1.0d-08/
       data two    /2.0d+00/
+      integer iat,jat,kat,nbond,ii,jk,i,icon,iang,ij
+      integer ijbnd,ik,jcon,jang,jj,ji,klbnd,kcon,jkbnd,lmbnd,lat,mat
+      integer mang,nang,nfound,nuser_specified
+      double precision cosb,rijsq,rij,rik,riksq,rjksq,rilsq,rjk,rklsq
+      double precision rkl
+      double precision distsq,dotprd
 c     
       distsq(iat,jat)=(c(1,iat)-c(1,jat))**2+
      1     (c(2,iat)-c(2,jat))**2+
@@ -2482,10 +2569,13 @@ c
       subroutine hnd_zmtyp3(nzvar,nizmat,izmat,mxizmt,c,
      1     numatom,nbnds,ijbnds,iibnds,nibnds, ! note rename nat -> numatom
      2     ijklto,numtor, max_tor_per_bond, zdone, oprint)
-      implicit double precision (a-h,o-z)
+c     implicit double precision (a-h,o-z)
+      implicit none
 #include "errquit.fh"
 #include "stdio.fh"
 #include "nwc_const.fh"
+      integer nzvar,nizmat,numatom,nbnds,numtor, max_tor_per_bond
+      integer ithree,mxatom,mxbnds,mxtors,mxizmt
       parameter (ithree=3)
       parameter (mxatom=nw_max_atom)        
       parameter (mxbnds= 8*mxatom)
@@ -2496,20 +2586,29 @@ c
       logical zdone, oprint, oopb, oi3d, oj3d
       logical geom_add_tor, bonds_span3d
       integer iw
-      dimension c(3,*)
-      dimension izmat(*)   
-      dimension ijklto(4,*)      
-      dimension ijbond(2,mxbnds)
-      dimension iibond(  mxatom)
-      dimension nibond(  mxatom)
-      dimension ijbnds(2,*)
-      dimension iibnds(  *)
-      dimension nibnds(  *)
+      double precision c(3,*)
+      integer          izmat(*)   
+      integer          ijklto(4,*)      
+      integer          ijbond(2,mxbnds)
+      integer          iibond(  mxatom)
+      integer          nibond(  mxatom)
+      integer          ijbnds(2,*)
+      integer          iibnds(  *)
+      integer          nibnds(  *)
+      double precision tol,two
       data tol    /1.0d-01/     ! sin(5.7 degrees) = 0.1
 *     data tol    /1.0d-08/
 c$$$  data zero   /0.0d+00/
 c$$$  data one    /1.0d+00/
       data two    /2.0d+00/
+      integer iat,jat,kat
+      double precision distsq,dotprd
+      double precision cosbk,cosbj,rjksq,rik,rijsq,rij,rjk,riksq
+      double precision rjl,rjlsq,rklsq,rkl
+      integer i,ii,ij,ijbnd,ik,il,jtor,jj,itor,ji,jk,jkbnd,jlbnd,jl
+      integer lat,klbnd,lmbnd,mat,mnbnd,mtor,nat,nbond,ntornonlin
+      integer ntor_for_ij,ntor
+
 c     
       distsq(iat,jat)=(c(1,iat)-c(1,jat))**2+
      1     (c(2,iat)-c(2,jat))**2+
@@ -3154,7 +3253,9 @@ c$$$ 9992 FORMAT(' NO LINEAR-BENDS FOUND ')
 c$$$ 9991 FORMAT(' LINEAR-BENDS NOT SET UP CURRENTLY . STOP . ')
 c$$$      END
       subroutine hnd_dparsc(a,la,c,lc)
-      implicit double precision (a-h,o-z)
+c     implicit double precision (a-h,o-z)
+      implicit none
+      integer la,lc,ic
       character*(*) a
       character*(*) c
       character*1   blk
@@ -3170,7 +3271,9 @@ c$$$      END
       return
       end
       subroutine hnd_dparsi(a,la,n)
-      implicit double precision (a-h,o-z)
+c     implicit double precision (a-h,o-z)
+      implicit none
+      integer la,i1,i2,isign,n,ia,ib,i
       character*(*) a
       character*1   char(12)
       data char     /'0','1','2','3','4','5','6','7','8','9',
@@ -3191,7 +3294,7 @@ c
       else
          isign= 1
       endif
-      na=i2-i1+1
+c     na=i2-i1+1
 c
       n=0
       do ia=i1,i2
@@ -3207,7 +3310,10 @@ c
       return
       end
       subroutine hnd_dparsr(a,la,x)
-      implicit double precision (a-h,o-z)
+c     implicit double precision (a-h,o-z)
+      implicit none
+      integer la,i1,i2,isign,ie1,ie2,ie,iexp,itmp,i,j
+      double precision dum,zero,ten,x
       logical rep
       character*(*) a
       character*1   char(17)
@@ -3288,9 +3394,11 @@ c
       return
       end
       SUBROUTINE HND_HNDERR(LERR,ERRMSG)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "errquit.fh"
 #include "stdio.fh"
+      integer lerr,iw
       CHARACTER*8  ERRMSG
       DIMENSION    ERRMSG(LERR)
 c
@@ -3306,8 +3414,10 @@ C
  9997 FORMAT(/,1X,72(1H-))
       END
       SUBROUTINE HND_GEOCLS(NFT)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "stdio.fh"
+      integer nft,iw
       CHARACTER*8  ERRMSG
       DIMENSION    ERRMSG(3)
       DATA ERRMSG  /'PROGRAM ','STOP IN ','-GEOCLS-'/
@@ -3320,8 +3430,10 @@ C
  9999 FORMAT(/,' ----- ERROR CLOSING UNIT ',I3,' IN -GEOCLS- . STOP .')
       END
       SUBROUTINE HND_GEOOPN(NFT,GEOFIL)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "stdio.fh"
+      integer nft,iw
       CHARACTER*80 GEOFIL
       CHARACTER*8  ERRMSG
       DIMENSION    ERRMSG(3)
@@ -3341,7 +3453,8 @@ C
      2                NUMVAR,PRSVAR,VARCHR,FRZVAR,FRZVAL,LST,
      3                IZMAT,IZ,IZFRZ,MXIZMT,NZMOD,DBUG,
      $     zvarname, zvarsign)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 C
 C     THIS ROUTINE TAKES THE CHARACTER STRING VALUES FOR THE
 C     -Z- MATRIX INPUT FOR THE INDICES, AND THE VALUES AND TRANSFORMS
@@ -3378,47 +3491,46 @@ C     S. CHIN: 11/08/90 - IBM KINGSTON, NY
 C
 #include "stdio.fh"
 #include "nwc_const.fh"
+      integer       IAT,NAT,NVAR,NZMOD,IZ
+      integer       MXATOM,MAXGEO,MAXZMT,MAXVAR,MXIZMT
       PARAMETER     (MXATOM=nw_max_atom)
       PARAMETER     (MAXGEO=MXATOM+1,MAXZMT=40,MAXVAR=nw_max_zmat)
       LOGICAL       DBUG
       LOGICAL       LST
       LOGICAL       READY
-      INTEGER     FLGZMT
       CHARACTER*80  PRSZMT
-      INTEGER     ZMTCHR
       CHARACTER*80  PRSVAR
-      INTEGER     VARCHR
-      INTEGER     ZMT
-      LOGICAL       FRZVAL
-      LOGICAL       FRZVAR
       LOGICAL       CART
-      LOGICAL       IZFRZ
       CHARACTER*8   ATNAME
       CHARACTER*1   PLUS
       CHARACTER*1   MINUS
       CHARACTER*8   ERRMSG
+      double precision xx,yy,zz,atnum
       COMMON/HND_XYZGEO/XX(MAXGEO),YY(MAXGEO),ZZ(MAXGEO),
      1                  ATNAME(MAXGEO),ATNUM(MAXGEO),CART(MAXGEO)
-      DIMENSION     ERRMSG(3)
-      DIMENSION     NUMZMT(       MAXGEO)
-      DIMENSION     PRSZMT(MAXZMT,MAXGEO)
-      DIMENSION     FLGZMT(MAXZMT,MAXGEO)
-      DIMENSION     ZMTCHR(MAXZMT,MAXGEO)
-      DIMENSION     PRSVAR(MAXZMT,MAXVAR)
-      DIMENSION     VARCHR(MAXZMT,MAXVAR)
-      DIMENSION     FRZVAR(       MAXVAR)
-      DIMENSION     NUMVAR(       MAXVAR)
-      DIMENSION     FRZVAL(     3,MAXGEO)
-      DIMENSION       ZVAL(     3,MAXGEO)
-      DIMENSION        ZMT(     5,MAXGEO)
-      DIMENSION       ZLST(     3,MAXGEO)
-      DIMENSION      IZMAT(MXIZMT)
-      DIMENSION      IZFRZ(MXIZMT)
+      DIMENSION         ERRMSG(3)
+      integer           NUMZMT(       MAXGEO)
+      DIMENSION         PRSZMT(MAXZMT,MAXGEO)
+      integer           FLGZMT(MAXZMT,MAXGEO)
+      integer           ZMTCHR(MAXZMT,MAXGEO)
+      DIMENSION         PRSVAR(MAXZMT,MAXVAR)
+      integer           VARCHR(MAXZMT,MAXVAR)
+      logical           FRZVAR(       MAXVAR)
+      integer           NUMVAR(       MAXVAR)
+      logical           FRZVAL(     3,MAXGEO)
+      double precision    ZVAL(     3,MAXGEO)
+      integer              ZMT(     5,MAXGEO)
+      double precision    ZLST(     3,MAXGEO)
+      integer            IZMAT(MXIZMT)
+      logical            IZFRZ(MXIZMT)
 
+      integer          i,j,k,JVAR,JAT
       character*(*) zvarname(*)
       double precision zvarsign(*)
 
       DATA ERRMSG /'PROGRAM ','STOP IN ','- ZDAT -'/
+      double precision ZERO,TWO,PIDEG
+      integer          iw
       DATA ZERO   /0.0D+00/
       DATA TWO    /2.0D+00/
       DATA PIDEG  /180.0D+00/
@@ -3917,7 +4029,8 @@ C
  9983 FORMAT(' SAME ATOM REFERRED TO MORE THAN ONCE ON THIS LINE')
       END
       SUBROUTINE HND_ZXYZ(NAT,NZMT,ZVAL)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 C
 C     ROUTINE COMPUTES :
 C
@@ -3939,26 +4052,35 @@ C          =3     ATOMS B-C-D ARE COLLINEAR
 C
 #include "stdio.fh"
 #include "nwc_const.fh"
+      integer       nat
+      integer       MXATOM,MAXGEO,MAXWRD,MAXVAR
       PARAMETER     (MXATOM=nw_max_atom)
       PARAMETER     (MAXGEO=MXATOM+1,MAXWRD=40,MAXVAR=nw_max_zmat)
       LOGICAL       DBUG
       LOGICAL       CART
       CHARACTER*8   ATNAME
       CHARACTER*8   ERRMSG
+      double precision xx,yy,zz,atnum
       COMMON/HND_XYZGEO/XX(MAXGEO),YY(MAXGEO),ZZ(MAXGEO),
      1                  ATNAME(MAXGEO),ATNUM(MAXGEO),CART(MAXGEO)
-      DIMENSION   NZMT(5,MAXGEO)
-      DIMENSION   ZVAL(3,MAXGEO)
+      integer            NZMT(5,MAXGEO)
+      double precision   ZVAL(3,MAXGEO)
       DIMENSION   ERRMSG(3)
-      double precision numd
+      double precision numd,zero,one,two,three,rcd,ccosp,ssinp
+      double precision costd,sintd,ccos,ssin,ccost,ssint
+      integer iw,na,nb,nc,nd,iat
       DATA ERRMSG    /'PROGRAM ','STOP IN ','- ZXYZ -'/
       DATA ZERO,ONE  /0.0D+00,1.0D+00/
       DATA TWO,THREE /2.0D+00,3.0D+00/
+      double precision t11,t12,t13,t21,t22,t23,t31,t32,t33
+      double precision phi,bet,alp,gam,dum,dot,rca,pifac,rab,rcb
+      double precision theta,thbcd
+      double precision xa,xb,xxd,yyd,zzd,ya,yb,za,zb
 C
       DBUG=.FALSE.
       IW = luout
 C
-      PIFAC=3.1415926536D+00/180.0D+00
+      PIFAC=acos(-1.0d0)/180.0D+00
       COSTD=-ONE/THREE
       SINTD= TWO/THREE* SQRT(TWO)
       NA   =0
@@ -4174,14 +4296,19 @@ C
  9997 FORMAT(' IAT = ',I4,' IS ALREADY SPECIFIED IN CARTESIAN SPACE.')
       END
       SUBROUTINE HND_PRSQ(V,M,N,NDIM)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "stdio.fh"
 C
 C     ----- PRINT OUT A SQUARE MATRIX -----
 C
 ***      COMMON/HND_LISTNG/LIST
-      DIMENSION V(NDIM,1)
-      DIMENSION IC(5),C(5)
+      integer m,n,ndim
+      double precision V(NDIM,1)
+      integer IC(5)
+      double precision C(5)
+      double precision vtol
+      integer icmax,iw,list,imax,imin,i,j,ii,idum,max
       DATA VTOL  /1.5D-01/
       DATA ICMAX /5/
 c
@@ -4242,16 +4369,22 @@ C
  9348 FORMAT(5(I5,F11.5))
       END
       SUBROUTINE HND_PREV(V,E,M,N,NDIM)
-      IMPLICIT REAL*8 (A-H,O-Z)
+c     IMPLICIT REAL*8 (A-H,O-Z)
+      implicit none
 #include "stdio.fh"
 C
 C     ----- PRINT OUT E AND V-MATRICES
 C
 ***      COMMON/HND_LISTNG/LIST
-      DIMENSION V(NDIM,1),E(1)
-      DIMENSION IC(5),C(5)
+      integer m,n,ndim
+      double precision V(NDIM,1),E(1)
+      integer IC(5)
+      double precision C(5)
+      double precision vtol
+      integer icmax
       DATA VTOL  /1.5D-01/
       DATA ICMAX /5/
+      integer iw,list,imax,imin,max,i,j,ii,idum
 C
       IW = LuOut
       LIST=1
@@ -4316,13 +4449,15 @@ C
  9348 FORMAT(5(I5,F11.5))
       END
       SUBROUTINE HND_PRTR(D,N)
-      IMPLICIT REAL*8 (A-H,O-Z)
+c     IMPLICIT REAL*8 (A-H,O-Z)
+      implicit none
 #include "stdio.fh"
 C
 C     ----- PRINT OUT A TRIANGULAR MATRIX -----
 C
 ***      COMMON/HND_LISTNG/LIST
-      DIMENSION D(1),DD(10)
+      integer n,iw,list,max,imax,imin,i,j,k,ii,jj,ij
+      double precision D(1),DD(10)
 C
       IW = LuOut
       LIST=1
@@ -4362,11 +4497,18 @@ C
  9248 FORMAT(I5,1X,7E15.8)
       END
       SUBROUTINE HND_TFTR(H,F,Q,T,IA,M,N,NDIM)
-      IMPLICIT REAL*8 (A-H,O-Z)
+c     IMPLICIT REAL*8 (A-H,O-Z)
+      implicit none
 C
 C     ----- H(M,M) = Q(DAGGER)(N,M) * F(N,N) * Q(N,M) -----
 C
-      DIMENSION H(1),F(1),Q(NDIM,1),T(1),IA(1)
+      integer m,n,ndim
+      double precision H(1),F(1),Q(NDIM,1),T(1)
+      integer IA(1)
+      integer ij,j,ik,max,i,k
+      double precision small,zero,dum,qij,hij
+      double precision ddot
+      external ddot
       DATA SMALL /1.0D-11/
       DATA ZERO  /0.0D+00/
       IJ = 0
@@ -4395,7 +4537,8 @@ C
       RETURN
       END
       SUBROUTINE HND_DIAGIV(A,VVEC,EIG,IA,NVEC,N,NDIM)
-      IMPLICIT REAL*8 (A-H,O-Z)
+c     IMPLICIT REAL*8 (A-H,O-Z)
+      implicit none
 #include "errquit.fh"
 #include "stdio.fh"
 C
@@ -4405,7 +4548,10 @@ C
 #include "mafdecls.fh"
 #include "global.fh"
       logical status
-      DIMENSION A(1),VVEC(NDIM,1),EIG(1),IA(1)
+      integer nvec,n,ndim
+      double precision A(1),VVEC(NDIM,1),EIG(1)
+      integer IA(1)
+      integer i_idbl,idbl,i_iint,k_iint,iw
 c
 c     ----- memory -----
 c
@@ -4438,13 +4584,22 @@ c
       END
       SUBROUTINE HND_GIVDIA(A,VEC,EIG,IA,N,NDIM,
      1           w,gamma,beta,betasq,p,q,iposv,ivpos,iord)
-      IMPLICIT REAL*8 (A-H,O-Z)
+c     IMPLICIT REAL*8 (A-H,O-Z)
+      implicit none
 C        -----  A GIVENS HOUSHOLDER MATRIX DIAGONALIZATION  -----
 C        -----  ROUTINE SAME AS EIGEN BUT WORKS WITH A      -----
 C        -----  LINEAR ARRAY.                               -----
-      DIMENSION A(1),VEC(NDIM,1),EIG(1),IA(1)
-      DIMENSION W(*),GAMMA(*),BETA(*),BETASQ(*)
-      DIMENSION P(*),Q(*),IPOSV(*),IVPOS(*),IORD(*)
+      integer n,ndim,n1,n2,nr,ik
+      integer IA(1)
+      integer IPOSV(*),IVPOS(*),IORD(*)
+      double precision A(1),VEC(NDIM,1),EIG(1)
+      double precision W(*),GAMMA(*),BETA(*),BETASQ(*)
+      double precision P(*),Q(*)
+      double precision zero,pt5,one,two,rhosq
+      double precision b,s,sum,sgn,temp,d,sqrts,g,a2,cosa,COSAP,DIF
+      double precision DIA,QJ,PPBR,PP,PPBS,SINA2,SINA,R1,R2,R12,SHIFT
+      double precision U,WTAW,WJ
+      integer i1,i,il,ij,ii,j,itemp,jj,k,lv,l,m,nrr,np,npas,nv,nt
       DATA ZERO,PT5,ONE,TWO /0.0D+00,0.5D+00,1.0D+00,2.0D+00/
       DATA RHOSQ /1.0D-22/
 C
diff --git a/src/geom/geom_input.F b/src/geom/geom_input.F
index 5e91dc9dc2..878b0a3afb 100644
--- a/src/geom/geom_input.F
+++ b/src/geom/geom_input.F
@@ -1833,15 +1833,17 @@ c
 c     
       end
       subroutine geom_autoz_input(geom,oprint)
-      implicit double precision (a-h,o-z)
+c     implicit double precision (a-h,o-z)
+      implicit none
 #include "errquit.fh"
 #include "nwc_const.fh"
 #include "geomP.fh"
+#include "util_params.fh"
       integer geom
       logical oprint
       logical geom_get_user_scale
       double precision scale, bohr
-      parameter (bohr=0.52917715d+00) !!!!!!!!!!!!!!!!!!!!!!!!!
+      parameter (bohr=cau2ang)
 c
       if (.not. geom_get_user_scale(geom, scale))
      $     call errquit('geom_autoz_input: scale?',0, GEOM_ERR)
@@ -2258,6 +2260,7 @@ c
 #include "nwc_const.fh"
 #include "geomP.fh"
 #include "stdio.fh"
+#include "util_params.fh"
       integer geom
       double precision zmat(*)  ! [output] nzvar of these
 c     
@@ -2266,10 +2269,10 @@ c     them in zmat.
 c     
       double precision b(3*max_cent) ! scratch dummy zmatrix
       double precision a(3*10)
-      double precision pi2, degree, bohr, zero, four, pio2, rtod
-      parameter (pi2=6.28318530717958d+00)
+      double precision pi2, degree, bohr, zero, four, pio2, rtod,pi
+c     parameter (pi2=6.28318530717958d+00)
       parameter (degree=360.0d+00)
-      parameter (bohr=0.52917715d+00) ! !!!!!!!!!!!!!!!!!!!!!!!!
+      parameter (bohr=cau2ang)
       parameter (zero=0.0d+00)
       parameter (four=4.0d+00)
 c     
@@ -2282,6 +2285,8 @@ c
       geom_compute_zmatrix = geom_check_handle(geom,'compute zmat')
       if (.not. geom_compute_zmatrix) return
 c     
+      pi    = dacos(-1.0d0)
+      pi2   = 2.0d0*pi
       pio2  =pi2/four
       rtod  =pi2/degree
 c     
@@ -2899,7 +2904,8 @@ c
      1                        IZ,IZMAT,NZMOD,ICFRZ,NCFRZ,
      2                        UNITS,zvarname,zvarsign,
      3                        found_cart)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "errquit.fh"
 C
 C     -----  PARAMETERS DEFINING MAXIMUM VALUES -----
@@ -2924,12 +2930,15 @@ C      ZVAL        =  DISTANCE, ANGLE AND TORSION ANGLE VALUE
 C       ZMT        =  I, J, K, L INDICES
 C
 #include "nwc_const.fh"
+      integer MXATOM,MXCOOR,MAXGEO,MAXWRD,MAXVAR,MAXPRM,MXIZMT,MAXLST
       PARAMETER     (MXATOM=nw_max_atom)
       PARAMETER     (MXCOOR=nw_max_coor)
       PARAMETER     (MAXGEO=MXATOM+1,MAXWRD=40,MAXVAR=nw_max_zmat)
       PARAMETER     (MAXPRM=100)
       PARAMETER     (MXIZMT=nw_max_izmat)
       PARAMETER     (MAXLST=10+1)
+      integer       ncenter,iz,izmat,nzmod,icfrz,ncfrz
+      double precision units
       CHARACTER*16  TAGS
       LOGICAL       GEOM_LST_PUT_COORD
       EXTERNAL      GEOM_LST_PUT_COORD
@@ -2990,32 +2999,36 @@ C
       CHARACTER*10  WR2VAR
       CHARACTER*10  WR1CON
       CHARACTER*10  WR2CON
+      integer ir,iw
+      double precision atnum
+      integer numchr,numwrd
+      double precision xx,yy,zz
       COMMON/HND_IOFILE/IR,IW
       COMMON/HND_FREERD/PRSWRD(40),NUMCHR(40),FLGWRD(40),NUMWRD
       COMMON/HND_XYZGEO/XX(MAXGEO),YY(MAXGEO),ZZ(MAXGEO),
      1                  ATNAME(MAXGEO),ATNUM(MAXGEO),CART(MAXGEO)
-      DIMENSION COORDS(3,*)
-      DIMENSION CHARGE(  *)
+      double precision COORDS(3,*)
+      double precision CHARGE(  *)
       DIMENSION   TAGS(  *)
       DIMENSION  IZMAT(  *)   
       DIMENSION  ICFRZ(  *)    
       DIMENSION  IZFRZ(MXCOOR)
       DIMENSION  CART0(MAXGEO)
-      DIMENSION  XXLST(MAXGEO,     2)
-      DIMENSION  YYLST(MAXGEO,     2)
-      DIMENSION  ZZLST(MAXGEO,     2)
+      double precision  XXLST(MAXGEO,     2)
+      double precision  YYLST(MAXGEO,     2)
+      double precision  ZZLST(MAXGEO,     2)
       DIMENSION  GHOST(       MAXGEO)
-      DIMENSION NUMZMT(       MAXGEO)
+      integer   NUMZMT(       MAXGEO)
       DIMENSION PRSZMT(MAXWRD,MAXGEO)
       DIMENSION FLGZMT(MAXWRD,MAXGEO)
       DIMENSION ZMTCHR(MAXWRD,MAXGEO)
       DIMENSION FRZVAL(     3,MAXGEO)
-      DIMENSION   ZVAL(     3,MAXGEO)
+      double precision   ZVAL(     3,MAXGEO)
       DIMENSION    ZMT(     5,MAXGEO)
-      DIMENSION   ZLST(     3,MAXGEO)
-      DIMENSION   ZSTP(     3,MAXGEO)
+      double precision   ZLST(     3,MAXGEO)
+      double precision   ZSTP(     3,MAXGEO)
       DIMENSION FRZVAR(       MAXVAR)
-      DIMENSION NUMVAR(       MAXVAR)
+      integer   NUMVAR(       MAXVAR)
       DIMENSION PRSVAR(MAXWRD,MAXVAR)
       DIMENSION FLGVAR(MAXWRD,MAXVAR)
       DIMENSION VARCHR(MAXWRD,MAXVAR)
@@ -3025,8 +3038,12 @@ C
       DIMENSION BLK(80)
       DIMENSION GH(4)
       DIMENSION ERRMSG(3)
+      integer i,iat,icenter,ilst,ilen,iflg,idigit,ierr,isymbl,irsav,ivar
+      integer iwrd,j,izmod,nvar,nat,nchr
       external inp_compare
       logical inp_compare 
+      double precision zero
+      double precision xxiat,yyiat,zziat
       EQUIVALENCE (BLNK80,BLK(1))
       EQUIVALENCE (CHREND,WRDEND)
       DATA ERRMSG /'PROGRAM ','STOP IN ','- ZGEO -'/
@@ -4046,29 +4063,37 @@ c
 c
       end
       SUBROUTINE HND_AUTSYM(ODONE,RTDB,THREQUIV,groupname,geom)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "nwc_const.fh"
       INTEGER     RTDB
       integer geom
+      double precision THREQUIV
       LOGICAL     ODONE
+      integer MXATOM,MXSYM
       PARAMETER   (MXATOM=nw_max_atom)
       PARAMETER   (MXSYM =120)
       CHARACTER*8 groupname
       LOGICAL     SOME
       LOGICAL     OUT
       LOGICAL     DBUG
+      integer ir,iw
+      double precision c,zan,v
+      integer nat
+      integer invt,nt,ntmax,ntwd,nosym
       COMMON/HND_IOFILE/IR,IW
       COMMON/HND_MOLXYZ/C(3,MXATOM),ZAN(MXATOM),V(3,MXATOM),NAT
       COMMON/HND_SYMTRY/INVT(MXSYM),NT,NTMAX,NTWD,NOSYM
-      DIMENSION   CI(3),AI(3,3)
-      DIMENSION   C1(3,MXATOM)
-      DIMENSION   C2(3,MXATOM)
-      DIMENSION   C3(3,MXATOM)
-      DIMENSION   V2(3,MXATOM)
-      DIMENSION   V3(3,MXATOM)
-      DIMENSION   AXS(3,3),EIG(3)
-      DIMENSION   RT(3,3)
-      DIMENSION   TR(3)
+      double precision   CI(3),AI(3,3)
+      double precision   C1(3,MXATOM)
+      double precision   C2(3,MXATOM)
+      double precision   C3(3,MXATOM)
+      double precision   V2(3,MXATOM)
+      double precision   V3(3,MXATOM)
+      double precision   AXS(3,3),EIG(3)
+      double precision   RT(3,3)
+      double precision   TR(3)
+      integer iat,i
 C
       DBUG=.FALSE.
       OUT =.FALSE. 
@@ -4149,7 +4174,8 @@ C
       END
       SUBROUTINE HND_MOLOPS(C,C1,C2,C3,NAT,AXS,EIG,RT,TR,ODONE,
      $     THREQUIV,groupname,geom,V,V2,V3)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "errquit.fh"
 C
 C     ----- ROUTINE DETECTS MOLECULAR SYMMETRY OPERATIONS -----
@@ -4161,11 +4187,13 @@ C
 #include "util.fh"
 #include "nwc_const.fh"
 #include "stdio.fh"
+      integer     nat
       integer     geom
       LOGICAL     ODONE
+      integer     mxordr, mxatom
       PARAMETER   (MXORDR=24)
       PARAMETER   (MXATOM=nw_max_atom)
-      CHARACTER*8 ERRMSG,groupname
+      CHARACTER*8 groupname
       LOGICAL     MUCH
       LOGICAL     DBUG
       LOGICAL     OUT
@@ -4174,9 +4202,7 @@ C
       LOGICAL     GOTX,GOTY,GOTZ
       LOGICAL     GOTONE
       LOGICAL     NOSYM
-      LOGICAL     PROPER
       LOGICAL     PRPAXS
-      LOGICAL     IMPROP
       LOGICAL     IMPAXS
       LOGICAL     SYMINV
       LOGICAL     INVERS
@@ -4192,11 +4218,8 @@ C
       LOGICAL     DEGNR3
       LOGICAL     CUBIC
       LOGICAL     C2ROT
-      LOGICAL     C2AXS
       LOGICAL     C4ROT
-      LOGICAL     C4AXS
       LOGICAL     S4ROT
-      LOGICAL     S4AXS
       LOGICAL     GRPOH
       LOGICAL     GRPTH
       LOGICAL     GRPTD
@@ -4205,44 +4228,60 @@ C
       logical     stopordrx
       logical     stopordry
       logical     stopordrz
-      INTEGER     AXORDR 
-      COMPLEX*16  QX,QY,QZ
-      COMPLEX*16  QZERO
-      COMPLEX*16  QDUMX,QDUMY,QDUMZ
-      COMPLEX*16  QDUMI,QDUMJ
+      double COMPLEX  QZERO
+      double COMPLEX  QDUMX,QDUMY,QDUMZ
+      double COMPLEX  QDUMI,QDUMJ
       character*16 atmlab
+      integer ir,iw
+      integer nuc
       COMMON/HND_IOFILE/IR,IW
       COMMON/HND_MOLNUC/NUC(MXATOM)
       common/hnd_mollab/atmlab(mxatom)
-      DIMENSION IEQU(MXATOM)
-      DIMENSION   QX(MXORDR),  QY(MXORDR),  QZ(MXORDR)
-      DIMENSION XMAG(MXORDR),YMAG(MXORDR),ZMAG(MXORDR)
-      DIMENSION PROPER(3)
-      DIMENSION IMPROP(3)
-      DIMENSION AXORDR(3)
-      DIMENSION NEWAXS(3)
-      DIMENSION NUORDR(3)
-      DIMENSION AXS(3,3),EIG(3)
-      DIMENSION  RT(3,3), TR(3)
-      DIMENSION PRM(3,3)
-      DIMENSION  C(3,*)
-      DIMENSION C1(3,*)
-      DIMENSION C2(3,*)
-      DIMENSION C3(3,*)
-      DIMENSION  V(3,*)
-      DIMENSION V2(3,*)
-      DIMENSION V3(3,*)
-      DIMENSION CM(3),AM(3,3)
-      DIMENSION AXM(3,3)
-      DIMENSION C2AXS(3)
-      DIMENSION C4AXS(3)
-      DIMENSION S4AXS(3)
-      DIMENSION ERRMSG(3)
-      double precision big,deter3
-      external deter3
+      integer IEQU(MXATOM)
+      double complex   QX(MXORDR),  QY(MXORDR),  QZ(MXORDR)
+      double precision XMAG(MXORDR),YMAG(MXORDR),ZMAG(MXORDR)
+      double precision fulltheta
+      logical PROPER(3)
+      logical IMPROP(3)
+      integer AXORDR(3)
+      integer NEWAXS(3)
+      integer NUORDR(3)
+      integer kat,k,JJC4,JJS4,jjc2,j,i
+      integer kaxis,jaxis,isteps,istep,iordr,iaxis
+      integer iiat,jjat,iaxm
+      double precision AXS(3,3),EIG(3)
+      double precision  RT(3,3), TR(3)
+      double precision PRM(3,3)
+      double precision  C(3,*)
+      double precision C1(3,*)
+      double precision C2(3,*)
+      double precision C3(3,*)
+      double precision  V(3,*)
+      double precision V2(3,*)
+      double precision V3(3,*)
+      double precision CM(3),AM(3,3)
+      double precision AXM(3,3)
+      logical C2AXS(3)
+      logical C4AXS(3)
+      logical S4AXS(3)
+      character*8 ERRMSG
+      dimension ERRMSG(3)
+      integer naxm,nsteps
+      double precision tmp,twopi,theta,tenm05,tenm04,three,four
+      double precision temp,two,step,rr,one
+      double precision pi,pi2,pi4
+      double precision dx,dy,dz,dirx,diry,dirz,dd
+      double precision d2,degree,d1,d3
+      double precision ethresh,dum,eps,disti,distj
+      double precision big,deter3,dlamch
+      integer nzer
+      external deter3,dlamch
       logical sameatm,ngotax(3)
-      complex*16 geom_powcmpl
+      double complex geom_powcmpl
       external geom_powcmpl
+      double precision xx,yy,zz,zero
+      double precision THREQUIV
+      integer iat,jat
       EQUIVALENCE (DEGNR3,CUBIC)
       DATA ERRMSG /'PROGRAM ','STOP IN ','-MOLOPS-'/
       DATA ZERO   /0.0D+00/
@@ -5548,7 +5587,8 @@ c
      2                     PROPRX,PROPRY,PROPRZ,
      3                     CUBIC,GRPOH,GRPTH,GRPTD,GRPT,GRPO,
      4                     AXORDR,TR,RT,ODONE,groupname)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 #include "util.fh"
       LOGICAL  ODONE
       LOGICAL  DBUG
@@ -5591,10 +5631,16 @@ c
       LOGICAL  MIRRYZ
       LOGICAL  MIRRZX
       LOGICAL  MIRRXY
+      integer  ir,iw
+      integer  igroup,naxis,linear
+      double precision complx,cmxsav
+      double precision dlamch
+      integer  igrsav,naxsav,linsav
+      integer  i,j
       COMMON/HND_IOFILE/IR,IW
       COMMON/HND_SYMMOL/COMPLX,IGROUP,NAXIS,LINEAR
       COMMON/HND_SYMNAM/GROUP
-      DIMENSION TR(3),RT(3,3)
+      double precision TR(3),RT(3,3)
       DATA BLANK  /'        '/
       DATA B      /' '/
       DATA S      /'S'/
@@ -5866,16 +5912,22 @@ C
      4       '                            ',3F10.5)
       END
       SUBROUTINE HND_MOLAXS(A,VEC,EIG,NVEC,N,NDIM)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
 C
 C     ----- ROUTINE TO SUBSTITUTE DIAGIV FOR DIAGONALIZATION -----
 C           OF SYMMETRIC 3X3 MATRIX A  IN TRIANGULAR FORM
 C
       CHARACTER*8 ERRMSG
+      integer ir,iw,n,ndim,nvec
       COMMON/HND_IOFILE/IR,IW
-      DIMENSION A(3,3),VEC(3,3),EIG(3)
-      DIMENSION AA(6),IA(3)
+      double precision A(3,3),VEC(3,3),EIG(3)
+      double precision AA(6)
+      integer IA(3)
       DIMENSION ERRMSG(3)
+      double precision zero,one,conv
+      double precision test,t
+      integer i,j
       DATA ERRMSG /'PROGRAM ','STOP IN ','-MOLAXS-'/
       DATA ZERO   /0.0D+00/
       DATA ONE    /1.0D+00/
@@ -5942,7 +5994,13 @@ C
       END
       SUBROUTINE HND_MOLGRP(XPT0,YPT0,ZPT0,XPT1,YPT1,ZPT1,         
      1                      XPT2,YPT2,ZPT2,XPT3,YPT3,ZPT3)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
+      double precision XPT0,YPT0,ZPT0,XPT1,YPT1,ZPT1
+      double precision XPT2,YPT2,ZPT2,XPT3,YPT3,ZPT3
+      integer       mxsym
+      integer       mxsymt
+      integer       mxioda
       PARAMETER     (MXSYM =120)
       PARAMETER     (MXSYMT=120*9)
       PARAMETER     (MXIODA=255)
@@ -5955,6 +6013,17 @@ C
       CHARACTER*8   CHRGRP
       CHARACTER*8   CHRDIR
       CHARACTER*8   CHRBLK
+      double precision   GRP
+      double precision   DRC,direct,blank,conv
+      integer       ir, iw
+      integer       idaf, nav, ioda
+      integer       invt, nt, ntmax, ntwd, nosym
+      double precision XSMAL,YSMAL,ZSMAL,XNEW,YNEW,ZNEW,XP,YP,ZP
+      double precision U1,U2,U3,V1,V2,V3,W1,W2,W3,X0,Y0,Z0
+      double precision complx
+      integer IGROUP,NAXIS,LINEAR
+      double precision group
+      double precision t
       COMMON/HND_IOFILE/IR,IW
       COMMON/HND_DAFILE/IDAF,NAV,IODA(MXIODA)
       COMMON/HND_SYMTRY/INVT(MXSYM),NT,NTMAX,NTWD,NOSYM
@@ -5969,11 +6038,19 @@ C
       DIMENSION     GRP(19)
       DIMENSION     DRC(2)
       DIMENSION     ERRMSG(3)
+      double precision pi2,zero,pt5,one,three,pi
+      integer i,j
+      double precision rho
+      double precision z1,z2,z3,z02,y1,y2,y3,y02,x1,x2,x3,x02,ww
+      double precision uu,u,test,sina,sinb,sign,dum,cosa,cosb
+      double precision beta,alpha,alph
+      integer n1,n2,n,itr,it,nn,ii
       EQUIVALENCE   (CHRGRP,GROUP )
       EQUIVALENCE   (CHRDIR,DIRECT)
       EQUIVALENCE   (CHRBLK,BLANK)
       EQUIVALENCE   (GRP(1),GRPCHR(1))
       EQUIVALENCE   (DRC(1),DRCCHR(1)) 
+      double precision tol
       DATA ERRMSG   /'PROGRAM ','STOP IN ','- PTGRP-'/
       DATA GRPCHR   /'C1      ','CS      ','CI      ',
      1               'CN      ','S2N     ','CNH     ',
@@ -5985,7 +6062,7 @@ C
       DATA CHRBLK   /'        '/
       DATA DRCCHR   /'NORMAL  ','PARALLEL'/
       DATA TOL      /1.0D-10/
-      DATA PI2      /6.28318530717958D+00/
+c     DATA PI2      /6.28318530717958D+00/
       DATA ZERO     /0.0D+00/
       DATA PT5      /0.5D+00/
       DATA ONE      /1.0D+00/
@@ -5997,6 +6074,8 @@ C
 C
 C     ----- GROUP INFO -----
 C
+      pi = acos(-1.0d0)
+      pi2 = 2.0d0*pi
       LINEAR=0
 C
       IGROUP=20
@@ -6571,15 +6650,23 @@ C
  9978 FORMAT(' THE MOLECULE IS LINEAR ')
       END
       SUBROUTINE HND_PRTSYM(N1,N2)
-      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+c     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      implicit none
+      integer n1,n2
+      integer MXSYM, MXSYMT, MXIODA
       PARAMETER (MXSYM =120)
       PARAMETER (MXSYMT=120*9)
       PARAMETER (MXIODA=255)
+      integer IR,IW
+      integer IDAF,NAV,IODA
+      integer INVT,NT,NTMAX,NTWD,NOSYM
+      double precision t
       COMMON/HND_IOFILE/IR,IW
       COMMON/HND_DAFILE/IDAF,NAV,IODA(MXIODA)
       COMMON/HND_SYMTRY/INVT(MXSYM),NT,NTMAX,NTWD,NOSYM
       COMMON/HND_SYMMAT/T(MXSYMT)
-      DIMENSION NN(MXSYM)
+      integer NN(MXSYM)
+      integer imax,imin,j,i,ni,nj
       IMAX=N1-1
   100 IMIN=IMAX+1
       IMAX=IMAX+4
@@ -6607,6 +6694,7 @@ C
 #include "global.fh"
 #include "nwc_const.fh"
 #include "stdio.fh"
+#include "util_params.fh"
       integer geom
       double precision alpha, ds(*)
       double precision err      ! [output] Returns the error
@@ -6652,8 +6740,8 @@ c
 c
 c     Form the target internals in bohr & degrees
 c
-      bohr = 0.52917715d0
-      deg  = 0.52917715d0*180d0/(4d0*atan(1d0))
+      bohr = cau2ang
+      deg  = bohr*180d0/(4d0*atan(1d0))
       if (.not. geom_compute_zmatrix(geom, q))
      $     call errquit('driver_u_c_f_i: geom?',0, GEOM_ERR)
       call geom_zmat_ico_scale(geom, dq, bohr, deg)
@@ -6785,6 +6873,7 @@ c
 #include "global.fh"
 #include "nwc_const.fh"
 #include "stdio.fh"
+#include "util_params.fh"
       integer geom
       external impose_constraints
 c     
@@ -6822,8 +6911,8 @@ c
      $     call errquit('driver_u_c_f_i: geom?',0, GEOM_ERR)
       ncart = nat * 3
 c     
-      bohr = 0.52917715d0
-      deg  = 0.52917715d0*180d0/(4d0*atan(1d0))
+      bohr = cau2ang
+      deg  = bohr*180d0/(4d0*atan(1d0))
 c     
       if (.not. ma_push_get(mt_dbl,ncart*nzvar,'mem bi',l_bi,k_bi))
      $     call errquit('opt_int_to_cart: ma', ncart*nzvar, MA_ERR)
@@ -7176,7 +7265,7 @@ c
       double precision theta, c2(3,nat)
       double precision pi, rr, dd, dz
       double precision threquiv
-      COMPLEX*16  QDUMI, QDUMJ,ctheta
+      double COMPLEX  QDUMI, QDUMJ,ctheta
 c
 #ifdef DEBUG
       some = .false.
@@ -7260,7 +7349,7 @@ c
       enddo
 c
       end
-      complex*16 function geom_powcmpl(c,s,n,big)
+      double complex function geom_powcmpl(c,s,n,big)
       implicit none
       double precision c,s,big
       integer n
diff --git a/src/geom/geom_print_ecce.F b/src/geom/geom_print_ecce.F
index aa8caa5f1a..04fc0ec0aa 100644
--- a/src/geom/geom_print_ecce.F
+++ b/src/geom/geom_print_ecce.F
@@ -1,4 +1,5 @@
       subroutine geom_print_rtdb_ecce(rtdb)
+      implicit none
 *
 * $Id$
 *