subroutine gfrcos(f,n,iop)
****************************
* DFFT at Gauss points ! Must call ginifft.f to initialize !
* iop=0.  physical to spectral, otherwise, spectral to phisical !
****************************
      implicit double precision (a-h,o-z)
      dimension f(0:n)
      common /wsave/wsave(6168)
      if (iop.eq.0) then

      call cosqb(n+1,f,wsave)
      f(0)=f(0)/(4*(n+1.d0))
      do i=1,n
         f(i)=f(i)/(2*(n+1.d0))
      enddo

      else

      tmp=f(0)
      call cosqf(n+1,f,wsave)
      do i=0,n
         f(i)=.5d0*(tmp+f(i))
      enddo
 
      endif
      return
      end