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

************* n  <=  2048  *****************

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

      if (iflag.ne.1) then
         val=1.d0/(n*n*n)
         do 55 k=1,n-1
            do 55 j=1,n-1
               do 55 i=1,n-1
 55      x(i,j,k)=val*x(i,j,k)
         val=.5d0*val
         do 65 j=1,n-1
            do 65 i=1,n-1
            x(i,j,0)=val*x(i,j,0)
            x(i,j,n)=val*x(i,j,n)
            x(i,0,j)=val*x(i,0,j)
            x(i,n,j)=val*x(i,n,j)
            x(0,i,j)=val*x(0,i,j)
 65      x(n,i,j)=val*x(n,i,j)
         val=.5d0*val
         do 70 i=1,n-1
            x(i,0,0)=val*x(i,0,0)
            x(i,n,0)=val*x(i,n,0)
            x(i,0,n)=val*x(i,0,n)
            x(i,n,n)=val*x(i,n,n)
            x(0,i,0)=val*x(0,i,0)
            x(n,i,0)=val*x(n,i,0)
            x(0,i,n)=val*x(0,i,n)
            x(n,i,n)=val*x(n,i,n)
            x(0,0,i)=val*x(0,0,i)
            x(n,0,i)=val*x(n,0,i)
            x(0,n,i)=val*x(0,n,i)
            x(n,n,i)=val*x(n,n,i)
 70      continue
         val=.5d0*val
         x(0,0,0)=val*x(0,0,0)
         x(0,0,n)=val*x(0,0,n)
         x(0,n,0)=val*x(0,n,0)
         x(0,n,n)=val*x(0,n,n)
         x(n,0,0)=val*x(n,0,0)
         x(n,0,n)=val*x(n,0,n)
         x(n,n,0)=val*x(n,n,0)
         x(n,n,n)=val*x(n,n,n)
      endif

      return
      end