subroutine vr2cos(x,wk,n,iflag)
**************************************************
*** 2-D Chebyshev transform !! MUST call VINIFFT(N) to initialize !!
***
*** x(i,j)  <-----> sum_{k,l}=0^n x_kl T_k(x_i) T_l(x_j)
*** 
*** PS. x_i=cos(i*pi/N) ! So x_0=1 and x_N=-1 !!
*** 
*** Inputs: x: 2-D array to be transformed
***   iflag=1: Spectral  to Physical
***   iflag=0: Physical  to Spectral  
***        wk: work array of dimension wk(0:n,0:n)
*** Output:
***         x: 2-D array after the transformation
************* n  <=  2048  *****************

      implicit double precision (a-h,o-z)
      dimension x(0:n,0:n),wk(0:n,0:n)
      common /wsave/wsave(6165)
      if (iflag.eq.1) then
*        val=.5d0*real(n)
      val=.25d0
         do 10 j=1,n-1
            do 10 i=1,n-1
 10      x(i,j)=val*x(i,j)
         val=2.d0*val
         do 15 i=1,n-1
            x(i,0)=val*x(i,0)
            x(i,n)=val*x(i,n)
            x(0,i)=val*x(0,i)
 15      x(n,i)=val*x(n,i)
         val=2.d0*val
         x(0,0)=val*x(0,0)
         x(n,0)=val*x(n,0)
         x(0,n)=val*x(0,n)
         x(n,n)=val*x(n,n)
         endif
         do 31 i=0,n
 31         call cost(n+1,x(0,i),wsave)
         do 30 i=0,n
            do 30 j=0,n
 30      wk(j,i)=x(i,j)
            do 32 i=0,n
         call cost(n+1,wk(0,i),wsave)
 32      continue
         do 40 i=0,n
            do 40 j=0,n
 40      x(j,i)=wk(i,j)
         if (iflag.ne.1) then
*         val=2.d0/real(n)
         val=1.d0/(1.d0*n*n)
         do 20 j=1,n-1
            do 20 i=1,n-1
 20      x(i,j)=val*x(i,j)
         val=.5d0*val
         do 25 i=1,n-1
            x(i,0)=val*x(i,0)
            x(i,n)=val*x(i,n)
            x(0,i)=val*x(0,i)
 25      x(n,i)=val*x(n,i)
         val=0.5d0*val
         x(0,0)=val*x(0,0)
         x(n,0)=val*x(n,0)
         x(0,n)=val*x(0,n)
         x(n,n)=val*x(n,n)
         endif
      return
      end