misc, minor code cleanups

src/ATOM/akeato.f

@@ -1,4 +1,4 @@
-      double precision function akeato (i,j,k)
+      double precision function akeato(i,j,k)
 c     angular coefficient by the direct coulomb integral fk for orbitals
 c     i and j
       implicit double precision (a-h,o-z)
@@ -20,7 +20,7 @@ double precision function akeato (i,j,k)
       entry bkeato (i,j,k)
 c     angular coefficient at the exchange coulomb integral gk
 
-      bkeato=0.0d 00
+      bkeato=0.0d0
       if (i .lt. j) then
          bkeato=afgk(j,i,k/2)
       elseif (i.gt.j) then

src/ATOM/fpf0.f

@@ -43,7 +43,8 @@ subroutine fpf0 ( iz, iholep, srho, dr, hx,
 c     always will use first spot to represent initial state
       nosc=1
       np = 251
-
+      xmult1 = 0.0d0
+      xmult2 = 0.0d0
       do 30 iorb =1, norb
         if (xnel(iorb) .gt.0.d0) then
 c         it is core orbital, check if it satisfies dipole selection

src/ATOM/inmuat.f

@@ -85,5 +85,5 @@ subroutine inmuat (ihole, xionin, iunf, xnval, iholep, xmag, iorb)
             if (nre(j).gt.0.or.nre(i).gt.0) ipl=ipl+1
  385     continue
  401  continue
- 999  return
+      return
       end

src/ATOM/muatco.f

@@ -13,10 +13,11 @@ subroutine muatco(xnval)
        external cwig3j
 
       do 511 i=1,30
-      do 511 j=1,30
-      do 511 k=0,3
- 511  afgk(i,j,k)=0.0d 00
- 601  do 701 i=1,norb
+         do 511 j=1,30
+            do 511 k=0,3
+               afgk(i,j,k)=0.0d00
+ 511  continue 
+      do 701 i=1,norb
          li= abs(kap(i))*2-1
          do 701 j=1,i
             lj= abs(kap(j))*2-1

src/ATOM/scfdat.f

@@ -42,9 +42,9 @@ subroutine scfdat ( ipr1, iph, nph, iz, ihole, xion, iunf, vcoul,
       common/comdir/ cl, dz, gg(251), ag(10), gp(251), ap(10), bid(783)
 c  gg,gp are the output from soldir
       common/itescf/ testy, rap(2), teste, nz, norb, norbsc
-      common/mulabk/ afgk
+cXX   common/mulabk/ afgk
       common/inelma/ nem
-      dimension afgk(30, 30, 0:3)
+cXX      dimension afgk(30, 30, 0:3)
       common/messag/ dlabpr, numerr
       character*8 dprlab, dlabpr
       common/ratom1/ xnel(30), en(30), scc(30), scw(30), sce(30),
@@ -72,6 +72,7 @@ subroutine scfdat ( ipr1, iph, nph, iz, ihole, xion, iunf, vcoul,
       endif
 
 c  initialize the data and test parameters
+      open_16 = .false.
       jfail = 0
       ibgp = 10
       numerr = 0

src/ATOM/vlda.f

@@ -14,24 +14,26 @@ subroutine vlda(ia, xnval,srho, srhovl,vtrho, ilast, idfock)
      1nq(30),kap(30),nmax(30)
       common/tabtes/hx,dr(251),test1,test2,ndor,np,nes,method,idim
 
-      do 10 i = 1,251
-        srhovl(i) = 0.0d0
- 10     srho(i) = 0.0d0
-
+      do 10 i = 1, 251
+        srho(i)   = zero
+        srhovl(i) = zero
+ 10   continue 
 c  find total and valence densities. Remove self-interaction if SIC
       do 50 j = 1, norb
-        a = xnel(j)
-        b = xnval(j)
-c      use to test SIC
+         a = xnel(j)
+         b = xnval(j)
+c     use to test SIC
 c       if (j .eq. ia) a=a-1.0d0
 c       if (j .eq. ia) b=b-1.0d0
-      do 50 i = 1,nmax(j)
-        srho(i) = srho(i) + a * (cg(i,j)**2+cp(i,j)**2)
- 50     srhovl(i) = srhovl(i) + b * (cg(i,j)**2+cp(i,j)**2)
+         do 50 i = 1, nmax(j)
+            srho(i)   = srho(i)   + a * (cg(i,j)**2+cp(i,j)**2)
+            srhovl(i) = srhovl(i) + b * (cg(i,j)**2+cp(i,j)**2)
+ 50   continue 
 
 c  constract lda potential. Put your favorite model into vbh.f.
 c  exch=5,6 correspond to 2 ways of core-valence separation of V_xc.
-      do 90 i = 1,251
+      rhoc = zero
+      do 90 i = 1, 251
         rho = srho(i) / (dr(i)**2)
         if (idfock.eq.5) then
 c          for exch=5 valence density*4*pi
@@ -40,18 +42,18 @@ subroutine vlda(ia, xnval,srho, srhovl,vtrho, ilast, idfock)
 c          for exch=6 core density*4*pi
            rhoc = (srho(i)-srhovl(i)) / (dr(i)**2)
         elseif (idfock.eq.1) then
-           rhoc = 0.0d0
+           rhoc = zero
         elseif (idfock.eq.2) then
            rhoc = srho(i) / (dr(i)**2)
         else
             call par_stop(' undefined idfock in subroutine vlda')
         endif
 
-        if (rho .gt. 0.0 ) then
+        if (rho .gt. zero ) then
            rs = (rho/3)**(-third)
-           rsc =101.0
-           if (rhoc .gt.0.0) rsc = (rhoc/3)**(-third)
-           xm = 1.0d0
+           rsc = 101.d0
+           if (rhoc .gt.zero) rsc = (rhoc/3)**(-third)
+           xm = one
 c          vbh and edp in Hartrees
            if (idfock.eq.5 .or. idfock.eq.2) then
 c             for exch=5, 2
@@ -64,7 +66,7 @@ subroutine vlda(ia, xnval,srho, srhovl,vtrho, ilast, idfock)
               vxcvl = vvbh - vdh
            elseif (idfock.eq.1) then
 c          for pure Dirac-Fock
-              vxcvl = 0.0d0
+              vxcvl = zero
            endif
 
 c   contribution to the total energy from V_xc:=\int d^3 r V_xc * rho/2

src/COMMON/fixdsp.f

@@ -13,7 +13,7 @@ subroutine fixdsp (dxorg, dxnew, dgc0, dpc0, dgcx, dpcx, jnew)
 
       dimension xorg(nrptx), xnew(nrptx)
 
-      parameter (xx00 = 8.8)
+      parameter (xx00 = 8.8d0)
 
 c     statement functions to do indexing.  delta is 'dx' for current
 c     grid.  jjj is index of grid point immediately before 'r'
@@ -32,6 +32,7 @@ subroutine fixdsp (dxorg, dxnew, dgc0, dpc0, dgcx, dpcx, jnew)
 
 c     The dgc and dpc arrays are zero beyond a certain point, usually
 c     inside the muffin tin radius.  Find this distance.
+      imax = 0
       do 100  i = 251, 1, -1
          if ( abs(dgc0(i)) .ge. 1.0d-11 .or. 
      1        abs(dpc0(i)) .ge. 1.0d-11 )  then
@@ -43,8 +44,7 @@ subroutine fixdsp (dxorg, dxnew, dgc0, dpc0, dgcx, dpcx, jnew)
    16 continue
 c     jmax is the first point where both dpc and dgc are zero in
 c     the original grid
-      jmax = imax + 1
-      if (jmax.gt.251) jmax = 251
+      jmax = min(imax, 250) + 1
 
       delta = dxorg
       do 10  j = 1, jmax
@@ -68,8 +68,8 @@ subroutine fixdsp (dxorg, dxnew, dgc0, dpc0, dgcx, dpcx, jnew)
 
 c     and zero the arrays past rmax
       do 32  j = jnew+1, nrptx
-         dgcx(j) = 0
-         dpcx(j) = 0
+         dgcx(j) = zero
+         dpcx(j) = zero
    32 continue
 
       return

src/COMMON/itoken.f

@@ -7,6 +7,7 @@ function itoken (word,flname)
       character*20 flname
       integer itoken
 
+      itoken = 0
       w = word(1:4)
       call upper(w)
 

src/COMMON/rdcmt.f

@@ -1,23 +1,22 @@
-      SUBROUTINE rdcmt(iUnt,Cmt)
-      INTEGER iUnt, i1
-      CHARACTER(300) line
-      CHARACTER(4) Cmt
-      CHARACTER TmpCmt(4), ch
-      LOGICAL CmtLin
+      subroutine rdcmt(iunt, cmt)
+      integer iunt, i1
+      character(4) cmt
+      character tmpcmt(4), ch
+      logical cmtlin
 
-      CmtLin = .true.
-      DO i1 = 1, 4
-         TmpCmt(i1) = Cmt(i1:i1)
-      END DO
- 5    CONTINUE
-      READ(iUnt,*,END=10) ch
-      DO i1 = 1, 4
-         IF(ch.eq.TmpCmt(i1)) goto 5
-      END DO
+      cmtlin = .true.
+      do i1 = 1, 4
+         tmpcmt(i1) = cmt(i1:i1)
+      end do
+ 5    continue
+      read(iunt,*,end=10) ch
+      do i1 = 1, 4
+         if(ch.eq.tmpcmt(i1)) goto 5
+      end do
 
-      BACKSPACE(iUnt)
+      backspace(iunt)
 
- 10   CONTINUE
+ 10   continue
 
-      RETURN
-      END
+      return
+      end

src/COMMON/rdpot.f

@@ -71,7 +71,6 @@ subroutine rdpot ( ntitle, title, rnrmav, xmu, vint, rhoint,
       dimension dum(13)
 
   10  format(a)
-   20 format (bn, i15)
 
       open (unit=3, file='pot.bin', status='old')
       read(3,30) ntitle, nph, npadx, nohole, ihole, inters, iafolp,

src/DEBYE/dwpar.h

@@ -3,7 +3,5 @@ c nlegx1 MUST be the same as legtot, the maximum number of scattering
 c       legs in a path
 c nphx1 MUST be the same as nphx, the maximum number of atomic species
 
-      parameter (natxdw = 200)
-      parameter (nlegx1 = 9)
-      parameter (nphx1=7)
+      parameter (natxdw = 200, nlegx1 = 9, nphx1=7)
 

src/DEBYE/sigrem.f

@@ -29,7 +29,7 @@
       subroutine sigem (sig2mx, sig2x, iem, tk, ipath, nleg, rat, sig2)
       implicit double precision (a-h, o-z)
 
-      include 'dwpar.h'
+      parameter (natxdw = 200, nlegx1 = 9, nphx1=7)
       include '../HEADERS/parallel.h'
 
 c feff parameters (from dim.h):
@@ -40,7 +40,7 @@ subroutine sigem (sig2mx, sig2x, iem, tk, ipath, nleg, rat, sig2)
       parameter (natx = natxdw)
 
 c local parameters:
-      parameter (amu0  = 1.660 54)
+      parameter (amu0  = 1.660 54d0)
       parameter (pi = 3.14159 26535 89793 23846 26433d0)
       parameter (nwx = 700)
 
@@ -430,7 +430,7 @@ subroutine dwrdin (rat, iphat, izph, nat,
 
       implicit double precision (a-h, o-z)
 
-      include 'dwpar.h'
+      parameter (natxdw = 200, nlegx1 = 9, nphx1=7)
 
 c feff parameters:
 c     parameter (nphx = 7)
@@ -603,10 +603,10 @@ subroutine dwrdin (rat, iphat, izph, nat,
                call wlog(slog)
                call par_stop('DWRDIN-3')
             endif
-            read(words(2),20,err=900)  izph(iph)
+            read(words(2), 20, err=900)  izph(iph)
 cc          No potential label if user didn't give us one
 cc          Default set above is potlbl=' '
-            if (nwords .ge. 3)  potlbl(iph) = words(3)
+            if (nwords .ge. 3)  potlbl(iph) = words(3)(:6)
          else
             write(slog,'(a,i8)') 
      .        'DWRDIN-4: Mode unrecognized, mode ', mode
@@ -709,7 +709,7 @@ subroutine rdspr(rat1, iz, natom, i0,
 
       implicit double precision (a-h, o-z)
 
-      include 'dwpar.h'
+      parameter (natxdw = 200, nlegx1 = 9, nphx1=7)
       include '../HEADERS/parallel.h'
 
 c feff parameters:
@@ -861,6 +861,8 @@ subroutine rdspr(rat1, iz, natom, i0,
             j=jj+1
             call chekin (jj, natom, line)
             read(words(3),30,err=900) str(i,j)
+            ix = 1
+            jx = 1
             if (str(i,j).lt.strx) then
                strx=str(i,j)
                ix=i
@@ -1196,7 +1198,7 @@ subroutine sang (i, j, k, rat1, si, sj, sk)
 
       implicit double precision (a-h, o-z)
 
-      include 'dwpar.h'
+      parameter (natxdw = 200, nlegx1 = 9, nphx1=7)
 
       parameter (natx = natxdw)
 
@@ -1334,7 +1336,7 @@ subroutine chekin (i, natom, line)
       subroutine sigrm (sig2mx, sig2x,ir1, ir2, tk,ipath,nleg,rat,sig2)
       implicit double precision (a-h, o-z)
 
-      include 'dwpar.h'
+      parameter (natxdw = 200, nlegx1 = 9, nphx1=7)
 
 c feff parameters (from dim.h):
 c     parameter (legtot=9)