subroutine histogram(npoints,datalabel,datavalue,ngroup,
     &grouplabel,nhistmark,histmark,nfreq,ncumumark,cumumark,
     &ncumu)
      implicit none
!npoints: the total number of points of the whole series
!datalabel: the lable of each datum
!datavalue: the value of each datum
!ngroup: the number of different data labeles, each label represents a group
!grouplabel: the label of each group
!nhistmark: the number of histogram band marks in each group
!histmark: the value of each histogram mark
!nfreq: the number of points in each histogram bands
!ncumumark: the number of marks for cumulative distribution
!cumumark: the marks for the cumulative distribution
!ncumu: the cumulative distribution in each group
      integer npoints,ngroup,nfreq(ngroup,npoints),nhistmark(ngroup),
     &ncumumark(ngroup),ncumu(ngroup,npoints),iorder(npoints),i,j,n
      character*50 datalabel(npoints),grouplabel(ngroup)
      double precision datavalue(npoints),histmark(ngroup,npoints),
     &cumumark(ngroup,npoints),vector(npoints)
      do i=1,ngroup
        do j=1,npoints
          nfreq(i,j)=0
          ncumu(i,j)=0
        enddo
        ncumumark(i)=0
      enddo
      do k=1,npoints
        do i=1,ngroup
          if(trim(datalabel(k)).eq.trim(grouplabel(i))then
!initially we set the ncumumark to the total number of points in each group. Later we will
!merge points with the same value
            ncumumark(i)=ncumumark(i)+1
            cumumark(i,ncumumark(i))=datavalue(k)
            goto 10
          endif
        enddo
10      continue
      enddo
      do i=1,ngroup
        call sort_shell(ncumumark(i),cumumark(i:i,1:ncumumark(i)),
     &iorder)
        do k=1,ncumumark(i)
          do j=1,nhistmark(i)-1
            if(j.eq.(nhistmark(i)-1))then
              if(cumumark(i,k).ge.histmark(i,j).and.cumumark(i,k).le.
     &histmark(i,j+1))then
              nfreq(i,j)=nfreq(i,j)+1
              goto 20
            else
              if(cumumark(i,k).ge.histmark(i,j).and.cumumark(i,k).lt.
     &histmark(i,j+1))then
                nfreq(i,j)=nfreq(i,j)+1
                goto 20
              endif
            endif
          enddo
20        continue
        enddo
!cumulative distribution. we have to merge points with equal values
        do j=1,ncumumark(i)
          vector(j)=cumumark(i,j)
        enddo
        k=ncumumark(i)
        ncumumark(i)=1
        cumumark(i,1)=vector(1)
        do j=2,k
          if(vector(j).ne.cumumark(i,ncumumark(i)))then
            ncumumark(i)=ncumumark(i)+1
            cumumark(i,ncumumark(i))=vector(j)
          endif
        enddo
        do j=1,ncumumark(i)
          ncumu(i,j)=0
          do n=1,k
            if(ivector(n).le.cumumark(i,j))then
              ncumu(i,j)=ncumu(i,j)+1
            endif
          enddo
        enddo
      enddo
      return
      end