      PROGRAM SlabVThick
c     *************************************************************************
c     *** Usage: slabvthick {-v vacuum_thickness}|{-s vacuum_thickness} \
c     ***                   struct-file 
c     ***
c     ***         -f
c     ***            vacuum_thickness is in fractional units 
c     ***            regarding the vacuum_thickness in old 
c     ***            struct file
c     ***
c     ***         -s
c     ***            vacuum thickness is in fractional units
c     ***            accoring to the thickness of slab
c     ***            (without vacuum)
c     *************************************************************************
      implicit real*8 (a-h,o-z)
      parameter (nato=1000)
      character*80  title
      character*4   typ, mode
      character*10  aname(nato)
      character*2   flag
      character*256 file, char_vacT_new
      dimension x(3,nato),klz(nato), rmt(nato), z(nato), r0(nato),
     $     rot_loc(nato,3,3)
      integer   ai(nato), npt(nato), mult(nato), isplit(nato)
      logical linv, lf

      if (iargc().ne.3) then
         print *,
     $'Usage: slabvthick {-v vacuum_thickness}|{-s vacuum_thickness} ',
     $'struct-file'
         stop
      endif

      Call GetArg(1,flag)
      lf=.false.
      if ( flag(1:2) .eq. '-f' ) then
         lf=.true.
      elseif ( flag(1:2) .eq. '-s' ) then
         lf=.false.
      else
         print *,
     $'Usage: slabvthick {-v vacuum_thickness}|{-s vacuum_thickness} ',
     $'struct-file'
         stop
      endif
      
      Call GetArg(2,char_vacT_new)
      vacT_new = c2d(char_vacT_new)
      print *,'VAC_T: ', vacT_new

      Call GetArg(3,file)
      len = index(file, ' ') - 1
      Open(unit=1, file=file(1:len), status='old')

      read(1,99900) title
      read(1,99901) typ,nat
      read(1,99902) mode
      read(1,99903) a,b,c,alpha,beta,gamma

      na=0
      do 10 i=1,nat
         na=na+1
         if(na.gt.nato) goto 9990
         read(1,99990) ai(i), x(1,na),x(2,na),x(3,na)
         read(1,99991) mult(i), isplit(i)
         do 11 im=2,mult(i)
            na=na+1
            if(na.gt.nato) goto 9990
            read(1,99990) ai(i), x(1,na),x(2,na),x(3,na)        
 11      continue
         read(1,99992) aname(i), npt(i), r0(i), rmt(i), z(i)
         read(1,99993) (rot_loc(i,1,iy), iy=1,3)
         read(1,99993) (rot_loc(i,2,iy), iy=1,3)
         read(1,99993) (rot_loc(i,3,iy), iy=1,3)
 10   continue

      goto 9999
 9990 write(*,*)'Too many atoms. Increase parameter nato in source ',
     *     'code.'
      stop
 9999 continue

c     ************************************************************
c     We must checked if slab is located on the middle of the slab
c     or not !!!
c     IF on the middle than it should remain in the middle
c
c     IF it is located on the middle the Zaverage must be 0.5 !!!!
c     ************************************************************

c     find an atom with lower Z and an average Z 
      Zsum=x(3,1)
      zmin=x(3,1)
      zmax=x(3,1)
      linv=.false.
      do i=2,na
         Zsum=Zsum + x(3,i)
         if (x(3,i).lt.zmin) zmin=x(3,i)
         if (x(3,i).gt.zmax) zmax=x(3,i)         
      enddo
      if ( dabs(Zsum/dble(na) - 0.5d0) .lt. 1.0d-5 ) linv=.true.
      slabT = zmax - zmin
      vacT  = 1.0d0 - slabT
      if (lf) then
         vacT_new = vacT * vacT_new
      else
         vacT_new = slabT * vacT_new
      endif

c     change c and x(3,*):
      c = c * vacT_new / vacT
      if (linv) then
c     *** we have center of inversion
         do i=1,na
            if ( x(3,i) .gt. 0.5d0 )
     $           x(3,i) = 0.5d0 + (x(3,1)-0.5d0) * vacT_new / vacT 
            if ( x(3,i) .lt. 0.5d0 )
     $           x(3,i) = 0.5d0 - (0.5d0-x(3,i)) * vacT_new / vacT
         enddo            
      else
c     *** no center of inversion
         do i=1,na
            x(3,i) = X(3,i) * vacT_new / vacT
         enddo
      endif

99900 format(a79)      
99901 format(a4,24x,i2)      
99902 format(13x,a4)
99903 format(6F10.7)
99904 format(2(3F15.9/),3f15.9)

99990 format(5x,i3,4x,f10.7,3x,f10.7,3x,f10.7)
99991 format(15x,i2,17x,i2)
99992 format(a10,5x,i5,5x,f10.8,5x,f10.5,5x,f5.2)
99993 format(20x,3f10.8)
               
c     ****************************
c     now write WIEN97 struct file
c     ****************************
      write(*,'(a)')
     $     'WIEN97 struct file generated by XCrySDen program'
      write(*,'(a4,a24,i2)') typ, 'LATTICE.NONEQUIV. ATOMS:',nat
      write(*,'(a13,a4)') 'MODE OF CALC=', mode
      write(*,'(6f10.6)') a, b, c, alpha, beta, gamma
      
      na=0
      do i=1,nat
         na=na+1
         write(*,'(a5,i3,a4,f10.8,a3,f10.8,a3,f10.8)')
     $        'ATOM=', ai(i),
     $        ': X=', x(1,na), ' Y=', x(2,na), ' Z=', x(3,na)
         write(*,'(a15,i2,a17,i2)') 'MULT=', mult(i),
     $        'ISPLIT=', isplit(i)
         do im=2,mult(i)
            na=na+1
            write(*,'(a5,i3,a4,f10.8,a3,f10.8,a3,f10.8)')
     $           'ATOM=', ai(i),
     $           ': X=', x(1,na), ' Y=', x(2,na), ' Z=', x(3,na)
         enddo
         read(1,99992) aname(i), npt(i), r0(i), rmt(i), z(i)
         write(*,'(a10,a5,i5,a5,f10.8,a5,f10.4,a5,f5.1)')
     $        aname(i),'NPT=',npt(i),'R0=',r0(i),
     $        'RMT=',rmt(i),'Z:', z(i)
         
         write(*,'(a20,3f10.7)') 'LOCAL ROT MATRIX:   ',
     $        (rot_loc(i,1,iy), iy=1,3)
         write(*,99993) (rot_loc(i,2,iy), iy=1,3)
         write(*,99993) (rot_loc(i,3,iy), iy=1,3)
      enddo
      write(*,'(i4,A)') 1,'      NUMBER OF SYMMETRY OPERATIONS'
      write(*,'(3i2,f10.7)') 1, 0, 0, 0.0
      write(*,'(3i2,f10.7)') 0, 1, 0, 0.0
      write(*,'(3i2,f10.7)') 0, 0, 1, 0.0
      write(*,'(i8)') 1
      END

c     =========================
      REAL*8 FUNCTION C2D(ARG)
      implicit real*8 (a-h,o-z)
      character*(*) arg

      print *,'c2d> arg:',arg
c     INTEGER PART
      i   = index(arg,'.')-1
      print *,'c2d> i(.):',i
      c2d = 0.0d0
      if(i .eq. 0)then
         i = len(arg)
         do m=1,i
            c2d = c2d + dble((ichar(arg(m:m)) - 48)*10**(i-m))
         enddo
         return
      endif

      do m=1,i
         c2d = c2d + dble((ichar(arg(m:m)) - 48)*10**(i-m))
      enddo
      print *,'c2d> c2d(INT):',c2d
c     REAL PART
      j   = index(arg,' ')-1
      print *,'c2d> j( )',j
      if(j .eq. 0) j=len(arg)
      do m=i+2,j
         c2d = c2d + dble(ichar(arg(m:m)) - 48)/dble(10**(m-(i+1)))
      enddo
      print *,'c2d> c2d(ALL):',c2d
      return 
      end
