PROGRAM uscale
      IMPLICIT none
      integer*4 NAT,N,i,j,ia,ix,ja,jx,ii,jj,ntype,nsum,no,nc,nh,nn,ico
      parameter (ntype=12)
      real*8,allocatable::f(:,:),s(:)
      integer*4,allocatable::q(:),bt(:,:),nt(:)
      character*20 OK
      real*8 a,sf(ntype)
      integer*4 nu(ntype)
      logical fscale
C
      write(6,609)
609   format(/,'          Harmonic force field universal scaling',/,/,
     1         '                   Input: FILE.X',/,
     1         '                          OLD.FC',/,
     1         '                          USCALE.PAR',/,
     1         '                  Output: FILE.FC',/,/)
      
      OK='f'
      if(iargc().gt.0)then
       call getarg(1,OK)
      endif
      if(OK.eq.'f')then
       write(6,*)'full scaling'
       fscale=.true.
      else
       write(6,*)'diagonal scaling only'
       fscale=.false.
      endif

      do 102 i=1,ntype
      nu(i)=0
102   sf(i)=1.0d0
      open(9,file='USCALE.PAR',status='old')
101   read(9,'(a)',end=901,err=901)OK
      if(OK(1:1).eq.'#')then
       read(9,*)i,a
       if(i.lt.1.or.i.gt.ntype)call report('SF out of range')
       sf(i)=a
      endif
      goto 101
901   close(9)
      write(6,*)'Parameters from USCALE.PAR read'
      write(6,*)

      open(9,file='FILE.X',status='old')
      read(9,*)
      read(9,*)nat
      allocate(q(nat),bt(nat,4),nt(nat),s(nat))
      do 2 i=1,nat
      read(9,*)q(i),a,a,a,(bt(i,j),j=1,4)
      nt(i)=0
      do 2 j=1,4
2     if(bt(i,j).ne.0)nt(i)=nt(i)+1
      close(9)
      write(6,*)'FILE.X ',nat,' atoms'

      do 3 ia=1,nat
      if(q(ia).eq.1)then
       s(ia)=sf(1)
       nu(1)=nu(1)+1       
      else
       if(q(ia).eq.6)then
        if(nt(ia).eq.2)then
         s(ia)=sf(2)
         nu(2)=nu(2)+1       
        else
         if(nt(ia).eq.3)then
          ico=0
          no=0
          nc=0
          nh=0
          nn=0
          do 31 j=1,nt(ia)
          if(q(bt(ia,j)).eq.1)nh=nh+1
          if(q(bt(ia,j)).eq.6)nc=nc+1
          if(q(bt(ia,j)).eq.7)nn=nn+1
          if(q(bt(ia,j)).eq.8)then
           no=no+1
           if(nt(bt(ia,j)).eq.1)ico=1
          endif
 31       continue
          if(ico.eq.1)then
           s(ia)=sf(12)
           nu(12)=nu(12)+1
          else
           s(ia)=sf(3)
           nu(3)=nu(3)+1       
          endif
         else
          if(nt(ia).eq.4)then
           s(ia)=sf(4)
           nu(4)=nu(4)+1       
          endif
         endif
        endif
       else
        if(q(ia).eq.7)then
         if(nt(ia).eq.2)then
          s(ia)=sf(5)
          nu(5)=nu(5)+1       
         else
          if(nt(ia).eq.3)then
           s(ia)=sf(6)
           nu(6)=nu(6)+1       
          else
           if(nt(ia).eq.4)then
            s(ia)=sf(7)
            nu(7)=nu(7)+1       
           endif
          endif
         endif
        else
         if(q(ia).eq.8)then
          if(nt(ia).eq.1)then
           s(ia)=sf(8)
           nu(8)=nu(8)+1       
          else
           if(nt(ia).eq.2)then
            s(ia)=sf(9)
            nu(9)=nu(9)+1       
           endif
          endif
         else
c         less frequent cases:
          if(q(ia).eq.16.and.nt(ia).eq.1)then
           s(ia)=sf(10)
           nu(10)=nu(10)+1       
          endif
          if(q(ia).eq.16.and.nt(ia).eq.2)then
           s(ia)=sf(11)
           nu(11)=nu(11)+1       
          endif
         endif
        endif
       endif
      endif
3     continue

      write(6,*)
      write(6,903)nu( 1),'  H- hydrogens'
      write(6,903)nu( 2),' -C- carbons  '
      write(6,903)nu( 3),' >C- carbons  '
      write(6,903)nu( 4),' >C< carbons  '
      write(6,903)nu( 5),' =N- nitrogens'
      write(6,903)nu( 6),' >N- nitrogens'
      write(6,903)nu( 7),' >N< nitrogens'
      write(6,903)nu( 8),' =O  oxygens  '
      write(6,903)nu( 9),' -O- oxygens  '
      write(6,903)nu(10),' =S  sulfurs  '
      write(6,903)nu(11),' -S- sulfurs  '
      write(6,903)nu(12),' >C= in CO    '
903   format(i6,a14)
      nsum=0
      do 103 i=1,ntype
103   nsum=nsum+nu(i)
      write(6,*)
      write(6,*)nat-nsum,' unasigned atoms'
      write(6,*)
      
      N=3*NAT
      allocate(f(N,N))
      OPEN(20,FILE='OLD.FC',STATUS='OLD')
      CALL READFF(N,f)
      CLOSE(20)
      write(6,*)'OLD.FC with force fields read'

      if(fscale)then
       do 104 ia=1,nat
       do 104 ja=1,nat
       a=dsqrt(s(ia)*s(ja))
       do 104 ix=1,3
       ii=3*(ia-1)+ix
       do 104 jx=1,3
       jj=3*(ja-1)+jx
104    f(ii,jj)=f(ii,jj)*a
      else
       do 105 ia=1,nat
       a=s(ia)
       do 105 ix=1,3
       ii=3*(ia-1)+ix
105    f(ii,ii)=f(ii,ii)*a
      endif
      write(6,*)'FF scaled'

      OPEN(20,FILE='FILE.FC')
      CALL WRITEFF(N,f)
      CLOSE(20)
      write(6,*)'FILE.FC written'

      end

      SUBROUTINE READFF(N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(4X,5D14.6)
      DO 31 I=1,N
      DO 31 J=I+1,N
   31 FCAR(I,J)=FCAR(J,I)
      RETURN
      END

      SUBROUTINE WRITEFF(N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(I4,5D14.6)
      RETURN
      END

      SUBROUTINE report(s)
      character*(*) s
      write(6,*)s
      stop
      end