parameter (maxmodel=50) 
      parameter (maxres=1000)
      parameter (maxclust=1000)
      character*32 model(maxmodel)
      integer imodel(maxmodel),iindex(maxmodel)
      integer ifragpair(maxres,maxmodel,maxmodel),
     &    ifragpair_new(maxres,maxmodel,maxmodel),
     & iflag(maxres),isumfmax(maxmodel),kmax(maxmodel),
     & isumf(maxmodel),isumfmaxmin(maxmodel),kmaxmin(maxmodel),
     & jmaxmin(maxmodel),isump(maxmodel),ifragpair_temp(maxres),
     & ifragpair_temp1(maxres),ninclust(maxclust),nlenclust(maxclust),
     & icluster(maxclust,maxclust),ifrag(maxres,maxclust),
     & imodelclust(maxclust)
      logical overlap
      character*1 charm
      character*8 ctemp
      character*128 outfile
 
      narg = iargc()
      if (narg.lt.2) then
        print '(2a)', "Usage: klapaucjusz nres nmodel",
     &   " [lencut] [minclust_size] [outfile] < infile > logfile"
        stop
      endif
      call getarg(1,ctemp)
      read(ctemp,*) nres
      call getarg(2,ctemp)
      read(ctemp,*) nmodel
      if (narg.gt.2) then
        call getarg(3,ctemp)
        read(ctemp,*) lencut
      else
        lencut=20
      endif
      if (narg.gt.3) then
        call getarg(4,ctemp)
        read(ctemp,*) minclust_size
      else
        minclust_size = 5
      endif
      if (narg.gt.4) then
        call getarg(5,outfile)
      else
        outfile = "klapaucjusz.out"
      endif
      open(2,file=outfile,status="unknown")
      print *,"nres",nres," nmodel",nmodel
      print *,"lencut",lencut," minclust_size",minclust_size
      print *," outfile ",outfile
      ifragpair_new=0 
      nclust = 0

      do i=1,nmodel
        do j=1,nmodel
          if (i.eq.j) cycle
          print *,"i",i," j",j
          read(*,'(a32,32x,1000i2)')model(i),(ifragpair(k,i,j),k=1,nres)
        enddo
      enddo
      do i=1,nmodel
        do j=i+1,nmodel
          write(*,'(2a32,1000i2)')model(i),model(j),
     &    (ifragpair(k,i,j),k=1,nres)
        enddo
      enddo
      isumcut=100
      iclust=0
      DO WHILE (isumcut.ge.lencut)
c search the longest fragment (number 1)
      iflag=0
      isummax=0
      nfrag=0
      do i=1,nmodel
        do j=1,i-1
c          write(*,'(2a32,1000i1)')model(i),model(j),(ifragpair(k,i,j)
c     &       -ifragpair(k,j,i),k=1,nres)
          isump=0
          do k=1,nres
            kk = ifragpair(k,i,j)
            if (kk.gt.0) isump(kk)=isump(kk)+1
            if (kk.gt.nfrag) nfrag=kk
          enddo
          isum=0
          do k=1,nfrag
            if (isump(k).gt.isum) then
              isum=isump(k)
              kk=k
            endif
          enddo 
c          write(*,'(2a32,i10)')model(i),model(j),isum
          if (isum.gt.isummax) then
            isummax=isum
            imodel(1)=i
            imodel(2)=j
            kfmax=kk
          endif
        enddo
      enddo
      write (*,*) "nfrag",nfrag," kfmax",kfmax
c      icut = max0(isummax*2/3,lencut)
      icut = lencut
      nflag=2
      iflag(imodel(1))=i
      iflag(imodel(2))=j
      write(*,*)"Maximum"
      write(*,'(2a32,i10)')model(imodel(1)),model(imodel(2)),isummax
      i1 = imodel(1)
      i2 = imodel(2)
      isumcut=isummax
      if (isummax.ge.lencut) then
          kkmod=kfmax
          kmaxminmax = kfmax
          ifragpair_temp=0
          ifragpair_temp1=0
          do k=1,nres
            if (ifragpair(k,i1,i2).eq.kfmax) then
              ifragpair_temp(k)=ifragpair(k,i1,i2)
              ifragpair_temp1(k)=ifragpair(k,i1,i2)
            endif
          enddo
          overlap=.true.
      else
        overlap=.false.
      endif
      write (*,*) ">>>>>> dowhile"
      do while (overlap)
        do i=1,nmodel
          if (iflag(i).gt.0) cycle
          do j=1,nflag
            jj = imodel(j)
            isumf=0
            do k=1,nres
              kk = ifragpair(k,i,jj)
              if (ifragpair_temp1(k).ne.kfmax .or. kk.eq.0) cycle
              isumf(kk)=isumf(kk)+1
            enddo 
            isumfmax(j)=0
            do k=1,nfrag
              if (isumf(k).gt.isumfmax(j)) then
                isumfmax(j)=isumf(k)
                kmax(j)=k
              endif
            enddo
          enddo ! j
          isumfmaxmin(i)=1000
          kmaxmin(i)=0
          jmaxmin(i)=0
          do j=1,nflag
            if (isumfmax(j).lt.isumfmaxmin(i)) then
              isumfmaxmin(i)=isumfmax(j)
              kmaxmin(i)=kmax(j)
              jmaxmin(i)=imodel(j)
            endif
          enddo
        enddo ! i
        isumfmaxminmax=0
        kmaxminmax=0
        imaxminmax=0
        jmaxminmax=0
        do i=1,nmodel
          if (iflag(i).gt.0) cycle
          if (isumfmaxmin(i).gt.isumfmaxminmax) then
            isumfmaxminmax=isumfmaxmin(i) 
            imaxminmax=i
            kmaxminmax=kmaxmin(i)
            jmaxminmax=jmaxmin(i)
          endif
        enddo
        print *,"isumfmaxminmax",isumfmaxminmax," icut",icut
        if (isumfmaxminmax.lt.icut) then
          overlap=.false.
        else
          isumcut = isumfmaxminmax
          iimod=imaxminmax
          jjmod=jmaxminmax
          kkmod=kmaxminmax
          nflag=nflag+1
          iflag(imaxminmax)=imaxminmax
          imodel(nflag)=imaxminmax 
          write(*,'(2a32,i10)')model(imodel(1)),model(imodel(2)),isummax
          write(*,'(1000i1)')(ifragpair(k,imodel(1),imodel(2)),k=1,nres)
          write(*,*) "isumfmaxminmax",isumfmaxminmax
          print *,"kmaxminmax",kmaxminmax
          write(*,'(a32)')(model(imodel(i)),i=1,nflag)
          write(*,'(1000i1)')
     &     (ifragpair(k,imaxminmax,jmaxminmax),k=1,nres)
          ifragpair_temp=0
          do k=1,nres
            if (ifragpair(k,imaxminmax,jmaxminmax).eq.kmaxminmax
     &       .and. ifragpair_temp1(k).gt.0)
     &      ifragpair_temp(k)=ifragpair(k,imaxminmax,jmaxminmax)
          enddo
          do k=1,nres
            ifragpair_temp1(k)=ifragpair_temp(k)
          enddo
         print *,"ifragpair_temp,ifragpair_temp1"
         print '(1000i1)',(ifragpair_temp(k),k=1,nres)
         print '(1000i1)',(ifragpair_temp1(k),k=1,nres)
         do i=1,nmodel
           do j=i+1,nmodel
             write(*,'(2a32,1000i1)')model(i),model(j),
     &       (ifragpair(k,i,j),k=1,nres)
           enddo
         enddo
        endif
      enddo ! while
      print *,"iimod",iimod," jjmod",jjmod," kkmod",kkmod,
     &   " isumcut",isumcut
      if (isumcut.ge.icut) then
          print *,"ifragpair_temp1 before zeroing"
          print '(64x,1000i1)',(ifragpair_temp1(k),k=1,nres)
          do i=1,nflag
            do j=1,i-1
              ii = imodel(i)
              jj = imodel(j)
              do k=1,nres
c                if (ifragpair_temp(k).eq.kmaxminmax .and. 
c     &           ifragpair_temp1(k).eq.kfmax) then
                if (ifragpair_temp1(k).gt.0) then 
                  ifragpair(k,ii,jj)=0
                  ifragpair(k,jj,ii)=0
                endif
              enddo
            enddo
          enddo
         print *,">>>>>>>>>>>>>>> IFRAGPAIR AFTER ZEROING <<<<<<<<<<<"
         do i=1,nmodel
           do j=i+1,nmodel
             write(*,'(2a32,1000i1)')model(i),model(j),
     &       (ifragpair(k,i,j),k=1,nres)
           enddo
         enddo

        if (nflag.ge.minclust_size) then

        nclust=nclust+1
        print *,"cluster",nclust
        ninclust(nclust)=nflag
        nlenclust(nclust)=isumcut
        do i=1,nflag
          icluster(i,nclust)=imodel(i)
        enddo
        ii=0
        do k=1,nres
c          if (ifragpair_temp(k).eq.kkmod .and. 
c     &     ifragpair_temp1(k).eq.kfmax) then
          if (ifragpair_temp1(k).gt.0) then
            ii=ii+1
c            print *,k
            ifrag(ii,nclust)=k
          endif
        enddo
        if (ii.ne.nlenclust(nclust)) write(*,*) "CHUJ NASTAPIL!!!",
     &    ii,nlenclust(nclust)
        do i=1,nflag
          do j=1,i-1
            do k=1,nres
              if (ifragpair_temp1(k).gt.0) then
                ifragpair_new(k,imodel(i),imodel(j))=nclust
                ifragpair_new(k,imodel(j),imodel(i))=nclust
              endif
            enddo
          enddo
        enddo

        endif ! nflag > 2

      endif
      ENDDO ! while
      do i=1,nmodel
        do j=i+1,nmodel
          write(*,'(2a32,1000a1)')model(i),model(j),
     &      (charm(ifragpair_new(k,i,j)),k=1,nres)
        enddo
      enddo
      write (*,*) "nclust",nclust
      do i=1,nclust
        write(*,*) "Cluster",i," ninclust",ninclust(i),
     &   " length",nlenclust(i)
        write(*,'(a32)') (model(icluster(j,i)),j=1,ninclust(i))
        write(*,'(i5)') (ifrag(j,i),j=1,nlenclust(i))
      enddo
      imodelclust=0
      do i=1,nclust
        do j=1,ninclust(i)
          imodelclust(icluster(j,i))=1
        enddo
      enddo
      print *,"Models in clusters"
      nmodel_out = 0
      do i=1,nmodel
        if (imodelclust(i).eq.1) then
          print '(a32)',model(i)
          nmodel_out = nmodel_out + 1
        endif
      enddo
      print *,"Models excluded"
      do i=1,nmodel
        if (imodelclust(i).eq.0) print '(a32)',model(i)
      enddo
c Write model information
      write(2,'(i5)') nmodel_out,nclust
      ii=0
      iindex=0
      do i=1,nmodel
        if (imodelclust(i).eq.1) then
          ii=ii+1
          write(2,'(a32)') model(i)
          iindex(i)=ii 
        endif
      enddo
c Write cluster information
      do i=1,nclust
        write(2,'(2i5)') ninclust(i),nlenclust(i)
        write(2,'(16i5)') (iindex(icluster(k,i)),k=1,ninclust(i))
        write(2,'(16i5)') (ifrag(k,i),k=1,nlenclust(i))
      enddo
      end
c---------------------------------------------------------------------------------
      character function charm(i)
      integer i
      if (i.eq.0) then
        charm = "."
      else
        charm = char(ichar("A")+i-1)
      endif
      return
      end