PROGRAM loggrd 
************************************************************************
* Lumps and interpolates an x-ordered x-y file to loggrid              *
* From standard input it reads:  
* input data filename
* xmn xmx #points/decade #columns
* The resulting x y file flushed to standard output.                   *
* It interpolates linearly over data-points around the output x-values.*
* Points that fall outside the available data-set are omitted.         *
* If the new grid is actually smaller than the original one,           *
* or if some datapoints are missing, a linear  interpolation is made   *
* between the two neighboring original data-points.                    *
************************************************************************
      parameter(nphys=100000,nc=30)
      real eps,sx,sy(nc),xx(nphys),yy(nphys,nc)
     *    ,xmx,xmn,xlast,xnext,ylast,ynext,wlow,whigh
     *    ,xnw(nphys),tmx,tmn,gmn,gmx
      INTEGER I,J,k,n,mm,nn,nw,ilast,ncol         
      character* 40 flin
      read(*,'(a40)') flin
      read(*,*) gmn,gmx,nw,ncol
      ncol=ncol-1
      gmn = alog10(gmn)
      gmx = alog10(gmx)
      nw = int((gmx-gmn)*nw)
      mm=100000 
      xmx=-1e20
      xmn=1e20
      open (23,file=flin)
      do 10 i=1,mm                                                  
       READ(23,*,END=11) xx(i) , (yy(i,k),k=1,ncol)
       if (xx(i).lt.xmn) xmn=xx(i)
       if (xx(i).gt.xmx) xmx=xx(i)
10    continue   
11    mm=i-1

      do 20 i=1,nw              
       xnw(i)=10**(gmn+(I-1)*(gmx-gmn)/nw)
20    continue   
21    nn=nw

******lump to tsp intervals 
      ilast=1
      eps=1e-8
      do 31 i=1,nn
       if (i.le.(nn-1)) then
        tmx=xnw(i+1)
       else
        tmx=xnw(i)
       endif
       if (i.ge.2) then
        tmn=xnw(i-1)
       else
        tmn=xnw(i)
       endif
       sx=0
       do 17 k=1,ncol
        sy(k)=0
17     continue
       n=0
       do 200 j=ilast,mm
	if (xx(j).lt.tmn) goto 200
        if (xx(j).ge.tmx) goto 210
	do 27 k=1,ncol
         sy(k)=sy(k)+yy(j,k)
27      continue
        sx=sx+xx(j)
        n=n+1
200     continue
210    ilast=j
       if (n.eq.0) then
*****  first omit all points that are outside interval       
        if ((xnw(i).lt.xmn).or.(xnw(i).gt.xmx)) then
	 goto 31
*****   otherwise take appropriate weighted average between
*****   last point and first subsequent point       
        else 
         xlast=xx(ilast-1) 
         xnext=xx(ilast)
	 do 37 k=1,ncol 
          ylast=yy(ilast-1,k)
          ynext=yy(ilast,k)
37       continue
         whigh=(xnw(i)-xlast)/(xnext-xlast)
         wlow=1.-whigh
         y=ylast*wlow+ynext*whigh
	 sx=(xlast+xnext)/2
        endif
       endif
       if (n.ge.1) then
        sx=sx/n
	do 47 k=1,ncol
	 sy(k)=sy(k)/n
47      continue
        write(*,'(30E12.4)') sx,(sy(k),k=1,ncol)
       endif
31    continue                   

      end