!     ==================================================================
      program main
!     ==--------------------------------------------------------------==
!     == read the input data from unit nfinp                          ==
!     == and charge density and generate file(s) in XSF format        ==
!     ==--------------------------------------------------------------==
      implicit none
      integer                      :: katm,natm,natm2,ktyp,ntyp,&
                                      ncord,ninv,brav,&
                                      iline,nline,idummy
      integer                      :: ii,jj,kk,iatm,it
      integer, allocatable         :: iwei(:),imdtyp(:),ityp(:)
      real(kind=8)                 :: pi,bohr,eps0,etoc
      real(kind=8)                 :: aaa,bbb,ccc,alpha,beta,gamma,&
                                      dummy,tmp
      real(kind=8)                 :: univol,sbc
      real(kind=8)                 :: altv(3,3),aux(3)
      real(kind=8)                 :: ca(3),cb(3),cc(3)
      real(kind=8)                 :: dipole_layer
      real(kind=8), allocatable    :: cps(:,:),pos(:,:),atomn(:)
      character(len=2),allocatable :: type(:)
      integer                      :: nfinp,nfout,ios
      logical                      :: fexist
      character(len=128)           :: filen,string
      character(len=2)             :: atmnm
!
      integer                      :: iu1,iu2,iu3,iout,iogp
      character(len=128)           :: infil1,infil2,infil3,&
                                      rhofil,drhfil,zdrfil
      character(len=128)           :: ctrl
      integer                      :: stdin=05,stdout=06,stderr=00
!
      integer :: i1,i2,i3,ia,i_gap,skip=0
      integer :: nl,nm,nn,il,im,in,itmp
      integer :: zval,typ, islab=0
      real(kind=8) :: aa1,aa2,dz,znow,ztot,mutot,mutot1,mutot2,mutot3,&
                      drhotot,a1,a2,cpst,&
                      muion,mui1,mui2,mui3
      real(kind=8) :: dpion,dpele,dptot,&
                      zmin,zmax, rhomin,rhomax,dmumin,dmumax,&
                      xmin,xmax,ymin,ymax,y2min,y2max,&
                      val,&
                      eang2debye ! e*Angstrom to Debye
      real(kind=8), allocatable :: rho1(:,:,:),rho2(:,:,:),rho3(:,:,:),&
                                   drho(:,:,:),drhot(:,:,:)
      real(kind=8), allocatable :: rho1t(:,:,:),rho2t(:,:,:),&
                                   rho3t(:,:,:)
      real(kind=8), allocatable :: drhoz(:),dmuz(:),drhoint(:)
      real(kind=8), allocatable :: rhoz(:)
      integer, allocatable :: na(:)
      integer :: idt,ndt,ispin,nspin,ilevel,nlevel
      integer :: nerr
      integer, allocatable :: errarr(:)
      character(len=128) :: rhonam
      real(kind=8) :: rhotmp
      integer, allocatable :: spin(:),level(:)
      character(len=8)     :: updw,occ,suffix
      integer            :: iarg,iargc,narg
      character(len=128) :: option
      character(len=128) :: codnam
!     ==--------------------------------------------------------------==
! ... code name
      codnam='chg2xsf'
!     ==--------------------------------------------------------------==
! ... default options
      skip=1  ! reduce grid data
      islab=1 ! for slab geometry
!     ==--------------------------------------------------------------==
      narg=0
      narg=iargc()
      if(narg.ge.1)then
        do iarg=1,narg
          call getarg(iarg,option)
          if(trim(option).eq.'-fine')then
            skip=0
          elseif(trim(option).eq.'-crystal'.or. &
     &           trim(option).eq.'-crys')then
            islab=0 
          elseif(trim(option).eq.'-h')then
            write(*,'(a,a,a)') &
     &      'Usage: ',trim(codnam),' [-fine][-crystal]'
            stop
          else
            write(*,'(a,a,a)') &
     &      trim(codnam),': illegal option ',trim(option)
            write(*,'(a,a,a)') &
     &      'Usage: ',trim(codnam),' [-fine][-crystal]'
            stop
          endif
        enddo
      endif
!     ==--------------------------------------------------------------==
! ... i/o
      nfinp=10
      nfout=20 
      iu1=12
      iu2=14
      iu3=16
      iout=30
!     ==--------------------------------------------------------------==
! ... constants
      pi=4.d0*atan(1.d0)
      bohr=0.529177d0
      eps0=8.85418781700d-12
      etoc=1.60217733000d-19
      eang2debye=1.d0/0.208194d0
!     ==--------------------------------------------------------------==
!     == read file names                                              ==
!     ==--------------------------------------------------------------==
      write(stderr,'(a,$)')'Enter the name of the STATE input file> '
      read(*,'(a)')filen
      write(stderr,'(a,$)')'Enter the name of the charge density file> '
      read(*,'(a)')infil1
      write(stderr,'(a,$)')'Enter the prefix for output file(s)> '
      read(*,'(a)')rhonam
!     ==--------------------------------------------------------------==
!     == show options                                                 ==
!     ==--------------------------------------------------------------==
      if(skip.eq.0)then
        write(*,'(/,a)')' *** FINE MESH WILL BE USED FOR DATAGRID ***'
      endif
      if(islab.eq.0)then
        write(*,'(/,a)')' *** CRYSTAL GEOMETRY ***'
      endif
!     ==--------------------------------------------------------------==
!     == read input file                                              ==
!     ==--------------------------------------------------------------==
      inquire(file=filen,exist=fexist)
      if(.not.fexist)then
        write(*,'(a,a,a)')'*** ERROR: ',trim(filen),' does not exist.'
        stop
      endif
      write(stderr,'(/,a,a,a,$)')'Reading ',trim(filen),' ... '
      open(nfinp,file=filen,iostat=ios,status='old')
      read(nfinp,*) idummy
      read(nfinp,*) dummy, dummy, ntyp, natm,natm2
      ktyp=ntyp
      katm=natm2
      allocate(cps(3,katm),pos(3,katm))
      allocate(iwei(katm),imdtyp(katm),ityp(katm))
      allocate(atomn(ktyp))
      allocate(type(ktyp))
      read(nfinp,*) idummy, brav
      read(nfinp,'(a)')string
      string=adjustl(string)
      if(string(1:9).eq.'Cartesian'.or.&
         string(1:9).eq.'cartesian'.or.&
         string(1:9).eq.'CARTESIAN'     )then
        do jj=1,3
          read(nfinp,*)(altv(ii,jj),ii=1,3)
        enddo        
      else
        backspace(nfinp)
        read(nfinp,*)aaa,bbb,ccc,alpha,beta,gamma
        altv(:,:)=0.d0
        alpha=alpha/180.d0*pi
        beta =beta /180.d0*pi
        gamma=gamma/180.d0*pi
        altv(1,1)=aaa
        altv(1,2)=bbb*cos(gamma)
        altv(2,2)=bbb*sin(gamma)
        tmp=(cos(alpha)-cos(beta)*cos(gamma))/sin(gamma)
        altv(1,3)=ccc*cos(beta)
        altv(2,3)=ccc*tmp
        altv(3,3)=sqrt(ccc**2-altv(1,3)**2+altv(2,3)**2)
      endif
!     read(nfinp,*) idummy, idummy, idummy
!     read(nfinp,*) ncord,ninv
      read(nfinp,'(a)') string
      string=adjustl(string)
      if(string(1:6).eq.'Dipole' .or.&
     &   string(1:6).eq.'dipole' .or.&
     &   string(1:6).eq.'DIPOLE'     )then
        read(nfinp,*) dipole_layer,dummy
        read(nfinp,*) idummy,idummy,idummy 
      else
        read(string,*)idummy,idummy,idummy
      endif
      read(nfinp,*)ncord,ninv
      if(ncord.eq.1)then
        do iatm=1,katm
          read(nfinp,*)(cps(ii,iatm),ii=1,3),&
                       iwei(iatm),imdtyp(iatm),ityp(iatm)
        enddo
      elseif(ncord.eq.0)then
        do iatm=1,katm
          read(nfinp,*)(pos(ii,iatm),ii=1,3),&
                       iwei(iatm),imdtyp(iatm),ityp(iatm)
        enddo
        do iatm=1,katm
          cps(1,iatm)= altv(1,1)*pos(1,iatm)&
                      +altv(1,2)*pos(2,iatm)&
                      +altv(1,3)*pos(3,iatm)
          cps(2,iatm)= altv(2,1)*pos(1,iatm)&
                      +altv(2,2)*pos(2,iatm)&
                      +altv(2,3)*pos(3,iatm)
          cps(3,iatm)= altv(3,1)*pos(1,iatm)&
                      +altv(3,2)*pos(2,iatm)&
                      +altv(3,3)*pos(3,iatm)
        enddo
      endif 
      do it=1,ktyp
        read(nfinp,*)atomn(it)
      enddo
!     do it=1,ktyp
!       type(it)=atmnm(int(atomn(it)))
!     enddo
      close(nfinp)
      allocate(na(ktyp))
      na(:)=0
      do iatm=1,natm
        na(ityp(iatm))=na(ityp(iatm))+1 
      enddo
!
      ca(:)=altv(:,1)
      cb(:)=altv(:,2)
      cc(:)=altv(:,3)
      if(brav.eq.1)then
! ... body centered lattice ...
        altv(:,1)=0.5d0*(-ca(:)+cb(:)+cc(:))
        altv(:,2)=0.5d0*(+ca(:)-cb(:)+cc(:))
        altv(:,3)=0.5d0*(+ca(:)+cb(:)-cc(:))
      elseif(brav.eq.2)then
! ... face centered lattice ...
        altv(:,1)=0.5d0*(+cb(:)+cc(:))
        altv(:,2)=0.5d0*(+ca(:)+cc(:))
        altv(:,3)=0.5d0*(+ca(:)+cb(:))
      elseif(brav.eq.3)then
! ... A-face centered ...
        altv(:,1)=ca(:)
        altv(:,2)=cb(:)
        altv(:,3)=0.5d0*(cb(:)+cc(:))
      elseif(brav.eq.4)then
! ... B-face centered ...
        altv(:,1)=ca(:)
        altv(:,2)=cb(:)
        altv(:,3)=0.5d0*(ca(:)+cc(:))
      elseif(brav.eq.5)then
! ... C-face centered ...
        altv(:,1)=ca(:)
        altv(:,2)=0.5d0*(ca(:)+cb(:))
        altv(:,3)=cc(:)
      elseif(brav.eq.6)then
! ... Rhombohedral center ...
        altv(:,1)=(2.0d0/3.0d0)*ca(:)+&
                  (1.0d0/3.0d0)*cb(:)+& 
                  (1.0d0/3.0d0)*cc(:)
        altv(:,2)=(1.0d0/3.0d0)*ca(:)+&
                  (1.0d0/3.0d0)*cb(:)+& 
                  (1.0d0/3.0d0)*cc(:)
        altv(:,3)=(1.0d0/3.0d0)*ca(:)-&
                  (2.0d0/3.0d0)*cb(:)+& 
                  (1.0d0/3.0d0)*cc(:)
      endif
!
! ... calculate the unit cell volume
      aux(:)=0.d0
      aux(1)=altv(2,2)*altv(3,3)-altv(3,2)*altv(2,3)
      aux(2)=altv(3,2)*altv(1,3)-altv(1,2)*altv(3,3)
      aux(3)=altv(1,2)*altv(2,3)-altv(2,2)*altv(1,3)
      univol=altv(1,1)*aux(1)+altv(2,1)*aux(2)+altv(3,1)*aux(3)
      univol=abs(univol)
      write(stderr,'(a)')'done.'
!     ==--------------------------------------------------------------==
!     == read charge densities                                        ==
!     ==--------------------------------------------------------------==
      inquire(file=infil1,exist=fexist)
      if(.not.fexist)then
        write(stderr,'(a,a)')trim(infil1),' not found.'
        stop
      endif
      write(stderr,'(/,a,a,a,$)')'Reading ',trim(infil1),' .'
      open(iu1,file=infil1,status='old')
! ... first check the number of data
      nspin=1;nlevel=1
      ndt=0
      do while(.true.)
        read(iu1,*,iostat=ios),nl,nm,nn,ispin,ilevel
        if(ios.gt.0)then
          nspin=1
          nlevel=1
          ndt=1
          exit
        elseif(ios.lt.0)then
          exit
        endif
        read(iu1,*,iostat=ios)(((rhotmp,il=1,nl),im=1,nm),in=1,nn) 
        if(ios.ne.0)then
          write(stderr,'(/,a)')&
          &'*** Error occurs while reading the charge density ***'
          stop
        endif
        ndt=ndt+1 
        if(ispin.eq.2)  nspin=2
        if(ilevel.eq.-1)nlevel=2
      enddo
      if(ndt.eq.0)then
        write(stderr,'(a)')'*** ERROR: DATA NOT FOUND.'
        stop
      endif
      write(stderr,'(a,$)')'.'
      rewind(iu1)
      allocate(spin(ndt),level(ndt))
      allocate(errarr(ndt))
      errarr(:)=0
      do idt=1,ndt
        read(iu1,*,iostat=ios),nl,nm,nn,spin(idt),level(idt)
        if(ios.gt.0)then
          exit
        endif
        read(iu1,*,iostat=ios)(((rhotmp,il=1,nl),im=1,nm),in=1,nn) 
        if(ios.ne.0)then
          write(stderr,'(/,a,a,i1)')&
          &'*** Error occurs while reading the charge density',&
     &     ' in dataset #',idt
          errarr(idt)=1
        endif
      enddo
!     ==--------------------------------------------------------------==
      write(stderr,'(a)')'. done.'
      write(*,'(/,a,a,a,i1)')'# of data blocks in ',trim(infil1),' : ',&
     &ndt
      do idt=1,ndt
        write(*,'(a,i1,1x,a,i2,1x,a,i2)')&
     &  'data #',idt,' : ispin :',spin(idt),' : ilevel:',level(idt)
      enddo
!     ==--------------------------------------------------------------==
      rewind(iu1)
      DATASET: do idt=1,ndt
        read(iu1,*,iostat=ios)nl,nm,nn
        if(ios.ne.0)then
          write(stderr,'(/,a,a,i1)')&
          &'*** Error occurs while reading the charge density',&
     &     ' in dataset #',idt
          errarr(idt)=1
          cycle DATASET
        endif
        if(.not.allocated(rho1))  allocate(rho1(nl,nm,nn))
        if(.not.allocated(rho1t)) allocate(rho1t(nl,nm,nn))
        write(stderr,'(/,a,i1,a,a,$)')&
     &  'Re-reading data #',idt,' in ',trim(infil1)
        read(iu1,*,iostat=ios)(((rho1(il,im,in),il=1,nl),im=1,nm),in=1,nn)
        if(ios.ne.0)then
          write(stderr,'(/,a,a,i1)')&
          &'*** Error occurs while reading the charge density',&
     &     ' in dataset #',idt
          errarr(idt)=1
          cycle DATASET
        endif
        rho1t(:,:,:)=0.d0
        do il=1,nl
          do im=1,nm
            do in=1,nn
              itmp=in+(nn/2*islab)
              if(itmp.le.0)itmp=itmp+nn
              if(itmp.gt.nn)itmp=itmp-nn
              rho1t(il,im,itmp)=rho1(il,im,in)
            enddo
          enddo
        enddo
!     ==--------------------------------------------------------------==
!     == generate a grid data in XSF format                           ==
!     ==--------------------------------------------------------------==
        suffix='.xsf'
! ... generate a file name
        if(nspin.eq.1)then
          if(nlevel.eq.2)then
            if(level(idt).eq.1)then
              rhofil=trim(rhonam)//'_occ'//trim(suffix)
            else
              rhofil=trim(rhonam)//'_unocc'//trim(suffix)
            endif
          else
            rhofil=trim(rhonam)//trim(suffix)
          endif
        elseif(nspin.eq.2)then
          if(nlevel.eq.2)then
            if(level(idt).eq.1)then
              if(spin(idt).eq.1)then
                rhofil=trim(rhonam)//'_up'//'_occ'//trim(suffix)
              else
                rhofil=trim(rhonam)//'_dw'//'_occ'//trim(suffix)
              endif
            else
              if(spin(idt).eq.1)then
                rhofil=trim(rhonam)//'_up'//'_unocc'//trim(suffix)
              else
                rhofil=trim(rhonam)//'_dw'//'_unocc'//trim(suffix)
              endif
            endif
          else
            if(spin(idt).eq.1)then
              rhofil=trim(rhonam)//'_up'//trim(suffix)
            else
              rhofil=trim(rhonam)//'_dw'//trim(suffix)
            endif
          endif
        endif
        write(*,'(a,a)')' ... Creating ',trim(rhofil)
        open(iout,file=rhofil,status='unknown')
        if(islab.eq.1)then
          write(iout,'(a)')' SLAB'
        else
          write(iout,'(a)')' CRYSTAL'
        endif
        write(iout,'(a)')' PRIMVEC'
        write(iout,'(3(f15.9,2x))')(bohr*altv(ii,1),ii=1,3)
        write(iout,'(3(f15.9,2x))')(bohr*altv(ii,2),ii=1,3)
        write(iout,'(3(f15.9,2x))')(bohr*altv(ii,3),ii=1,3)
        write(iout,'(a)')' PRIMCOORD'
        write(iout,'(i3,i2)')natm2,1
        do ia=1,natm2
          write(iout,'(i5,3(f15.9,2x))')&
          int(atomn(ityp(ia))),(bohr*cps(ii,ia),ii=1,3) 
        enddo
        write(iout,'(a)')' BEGIN_BLOCK_DATAGRID_3D'
        write(iout,'(a)')' data_grid'
        write(iout,'(a,i1)')' DATAGRID_3D_DENSITY#',1
        write(iout,'(3i4)')nl/(skip+1)+1,nm/(skip+1)+1,nn/(skip+1)+1
        write(iout,'(2x,3(1x,e17.10e2))')&
        0.0, 0.0, -0.5d0*altv(3,3)*bohr*islab
        write(iout,'(3(f15.9,2x))')(bohr*altv(ii,1),ii=1,3)
        write(iout,'(3(f15.9,2x))')(bohr*altv(ii,2),ii=1,3)
        write(iout,'(3(f15.9,2x))')(bohr*altv(ii,3),ii=1,3)
        do in=0,nn,skip+1
          do im=0,nm,skip+1
            do il=0,nl,skip+1
              i1=mod(il,nl)+1
              i2=mod(im,nm)+1
              i3=mod(in,nn)+1
!             write(iout,'(e17.10e2)')rho1(i1,i2,i3)
              write(iout,'(e17.10e2)')rho1t(i1,i2,i3)
            enddo
            write(iout,*)
          enddo
        enddo
        write(iout,'(a,i1)')' END_DATAGRID_3D_DENSITY#',1
        write(iout,'(a)')' END_BLOCK_DATAGRID_3D'
        close(iout)
!     ==--------------------------------------------------------------==
      enddo DATASET
      nerr=0
      do idt=1,ndt
        nerr=nerr+errarr(idt)
      enddo
!     ==--------------------------------------------------------------==
!     == deallocate                                                   ==
!     ==--------------------------------------------------------------==
      deallocate(cps,pos,atomn)
      deallocate(iwei,imdtyp,ityp)
!     ==--------------------------------------------------------------==
      deallocate(rho1)
      deallocate(spin,level,errarr)
!     ==--------------------------------------------------------------==
      if(nerr.eq.0)then
        write(*,'(/,a)')'Program successfully ended'
!       write(*,'(a,a,a)')&
!       'Generation of ',trim(rhofil),' successfully done.'
      else
        write(*,'(/,a)')'Program failed due to a reading error.'
      endif
!     ==--------------------------------------------------------------==
      stop
      end
!     ==================================================================