subroutine gr3cos(x,wk,n,iflag)
**************************************************
*** 3-D Chebyshev transform !! MUST call GINIFFT(N) to initialize !!
***  iflag=1: Spectral  to Physical
***  iflag=0: Physical  to Spectral  
***       wk: work array of dimension wk(0:n)
************* n  <=  2048  *****************

      implicit double precision (a-h,o-z)
      dimension x(0:n,0:n,0:n),wk(0:n)
      common /wsave/wsave(6168)

      do 10 j=0,n
 10   call gr2cos(x(0,0,j),wk,n,iflag)
      do 20 k=0,n
         do 20 i=0,n
         do 30 j=0,n
 30      wk(j)=x(k,i,j)
         call gfrcos(wk,n,iflag)
         do 40 j=0,n
 40      x(k,i,j)=wk(j)
 20   continue
      return
      end