subroutine gr2cos(x,wk,n,iflag)
**************************************************
*** 2-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)
***
*** Not working yet. May need rewrite the case for iflag=0 !
************* n  <=  2048  *****************

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

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