sortet tlist output from dtr extract_sqavesort.f90 check the dq criterium, still fixed

program extract_sqavetsort

! very simple extract sqt tables from report tex file and write it to table

! compile: gfortran -o extract_sqavetsort extract_sqavetsort.f90 -lslatec_gfortran

! usage:

!

! extract_sqavetsort < report1234_11.x.dtr  > mylist

!

! Author: Michael Monkenbusch, Forschungszentrum Juelich

!

 implicit none

 integer                :: datr_output = 6

 integer                :: tab_output  = 6

 character(len=128)     :: comment

 character(len=4096)    :: line_buf

 integer                :: ios, i, j, ier

 integer, parameter     :: maxtab =5000

 integer, parameter     :: maxqs  =100

 integer, parameter     :: maxtaus=300

 integer                :: ntab, nq

 double precision       :: tab(7,maxtab)  ! keeping the value enables later sorting if required

 double precision       :: sqt(maxqs, maxtaus,4) = 0d0

 double precision       :: qval0(maxqs)          = 0d0

 double precision       :: qval(maxqs)           = 0d0

 integer                :: nqval(maxqs)          = 1

 integer                :: ntau(maxqs) 

 integer                :: iw_perm(maxtaus)          

 double precision       :: dq = 0.0075d0

 double precision       :: qref, qave, qvar

 ntab = 1

 nq   = 1

 ntau = 1

        read(5,'(a)',iostat=ios) comment

        if(ios .ne. 0) stop "data not readable!"

 r1:  do 

        read(5,'(a)',iostat=ios) line_buf

        if(ios .ne. 0) exit

        read(line_buf,*,iostat=ios) tab(:,ntab)

        if(ios .ne. 0) cycle

!!        write(datr_output,'(4f18.7)') tab(7,ntab), tab(1:3,ntab)

        ntab = ntab+1

      enddo r1

      ntab = ntab-1

      nq = 1

      qval(1)  =  tab(7,1)

      qval0(1) =  tab(7,1)

q1:   do i = 2, ntab

        if(abs(qval0(nq)-tab(7,i)) > dq) then

          nq = nq+1

          qval(nq)  = tab(7,i) 

          qval0(nq) = tab(7,i) 

        else

          qval(nq) = qval(nq)  + tab(7,i) 

          nqval(nq)= nqval(nq) + 1

        endif 

      enddo q1

      qval = qval/nqval

!!      write(*,'(2f12.6)') qval(1:nq), qval0(1:nq)

qt:   do j=1,nq

         ntau(j) = 0  

tt:      do i=1,ntab

             if(abs(qval0(j)-tab(7,i)) > dq) cycle tt

             ntau(j) = ntau(j)+1

             sqt(j,ntau(j),1:3) = tab(1:3,i) 

             sqt(j,ntau(j),  4) = tab(7,i) 

         enddo tt

      enddo qt

!      write(*,'(i8)') ntau(1:nq)

      write(tab_output,'(a)') trim(comment)

      write(tab_output,'(a)') "    tau/ns    s(q,t)/s(q) sqt-error             q-exact"

      do j=1,nq

! q-variation  

        qave = sum(sqt(j,:,4) ) / ntau(j)

        qvar = sqrt(sum( (sqt(j,1:ntau(j),4)-qave)**2) / ntau(j))

        write(tab_output,'("q=",f12.6,"    (",f12.6," +- ",f8.6," variance)")') qval(j), qave, qvar

        iw_perm = [(i,i=1,ntau(j))]

        call  DPSORT (sqt(j,:,1), ntau(j), iw_perm,  1 , ier)

        do i=1,ntau(j)

          write(tab_output,'(f10.4,3x,f10.7,2x,f10.7,10x,f10.5)') sqt(j,iw_perm(i),1:4)

        enddo

      enddo

end program extract_sqavetsort