      subroutine fourt (data,nn,ndim,isign,iform,work)

      dimension       data(1)    ,nn(1)      ,ifact(32)  ,work(1)
      data np0/0/,nprev/0/
      data twopi/6.2831853071796/,rthlf/0.70710678118655/
      if (ndim-1) 232,101,101
  101	ntot = 2
      do 103 idim=1,ndim
         if (nn(idim)) 232,232,102
  102    ntot = ntot*nn(idim)
  103 continue
c
c     maipn loop for each dimension
c
      np1 = 2
      do 231 idim=1,ndim
         n = nn(idim)
         np2 = np1*n
         if (n-1) 232,230,104
c
c     is n a power of two and if not, what are its factors
c
  104  m = n
       ntwo = np1
       if = 1
       idiv = 2
  105  iquot = m/idiv
       irem = m-idiv*iquot
       if (iquot-idiv) 113,106,106
  106  if (irem) 108,107,108
  107  ntwo = ntwo+ntwo
       ifact(if) = idiv
       if = if+1
       m = iquot
       go to 105
  108  idiv = 3
       inon2 = if
  109  iquot = m/idiv
       irem = m-idiv*iquot
       if (iquot-idiv) 115,110,110
  110  if (irem) 112,111,112
  111  ifact(if) = idiv
       if = if+1
       m = iquot
       go to 109
  112  idiv = idiv+2
       go to 109
  113  inon2 = if
       if (irem) 115,114,115
  114  ntwo = ntwo+ntwo
       go to 116
  115  ifact(if) = m
c
c     separate four cases--
c        1. complex transform or real transform for the 4th, 9th,etc.
c           dimensions.
c        2. real transform for the 2nd or 3rd dimension.  method--
c           transform half the data, supplying the other half by con-
c           jugate symmetry.
c        3. real transform for the 1st dimension, n odd.  method--
c           set the imaginary parts to zero.
c        4. real transform for the 1st dimension, n even.  method--
c           transform a complex array of length n/2 whose real parts
c           are the even numbered real values and whose imaginary parts
c           are the odd numbered real values.  separate and supply
c           the second half by conjugate symmetry.
c
  116  icase = 1
       ifmin = 1
       i1rng = np1
       if (idim-4) 117,122,122
  117  if (iform) 118,118,122
  118  icase = 2
       i1rng = np0*(1+nprev/2)
       if (idim-1) 119,119,122
  119    icase = 3
         i1rng = np1
         if (ntwo-np1) 122,122,120
  120    icase = 4
         ifmin = 2
         ntwo = ntwo/2
         n = n/2
         np2 = np2/2
         ntot = ntot/2
         i = 1
         do 121 j=1,ntot
            data(j) = data(i)
            i = i+2
  121    continue
c
c     shuffle data by bit reversal, since n=2**k.  as the shuffling
c     can be done by simple interchange, no working array is needed
c
  122    if (ntwo-np2) 132,123,123
  123    np2hf = np2/2
         j = 1
         do 131 i2=1,np2,np1
            if (j-i2) 124,127,127
  124       i1max = i2+np1-2
            do 126 i1=i2,i1max,2
               do 125 i3=i1,ntot,np2
                  j3 = j+i3-i2
                  tempr = data(i3)
                  tempi = data(i3+1)
                  data(i3) = data(j3)
                  data(i3+1) = data(j3+1)
                  data(j3) = tempr
                  data(j3+1) = tempi
  125          continue
  126       continue
  127       m = np2hf
  128       if (j-m) 130,130,129
  129       j = j-m
            m = m/2
            if (m-np1) 130,128,128
  130       j = j+m
  131    continue
         go to 142
c
c     shuffle data by digit reversal for general n
c
  132    nwork = 2*n
         do 141 i1=1,np1,2
            do 140 i3=i1,ntot,np2
               j = i3
               do 138 i=1,nwork,2
                  if (icase-3) 133,134,133
  133             work(i) = data(j)
                  work(i+1) = data(j+1)
                  go to 135
  134             work(i) = data(j)
                  work(i+1) = 0.
  135             ifp2 = np2
                  if = ifmin
  136             ifp1 = ifp2/ifact(if)
                  j = j+ifp1
                  if (j-i3-ifp2) 138,137,137
  137             j = j-ifp2
                  ifp2 = ifp1
                  if = if+1
                  if (ifp2-np1) 138,138,136
  138          continue
               i2max = i3+np2-np1
               i = 1
               do 139 i2=i3,i2max,np1
                  data(i2) = work(i)
                  data(i2+1) = work(i+1)
                  i = i+2
  139          continue
  140       continue
  141    continue
c
c     main loop for factors of two.  perform fourier transforms of
c     length four, with one of length two if needed.  the twiddle factor
c     w=exp(isign*2*pi*sqrt(-1)*m/(4*mmax)).  check for w=isign*sqrt(-1)
c     and repeat for w=w*(1+isign*sqrt(-1))/sqrt(2).
c
  142    if (ntwo-np1) 174,174,143
  143    np1tw = np1+np1
         ipar = ntwo/np1
  144    if (ipar-2) 149,146,145
  145    ipar = ipar/4
         go to 144
 146	 do 148 i1=1,i1rng,2
            do 147 k1=i1,ntot,np1tw
               k2 = k1+np1
               tempr = data(k2)
               tempi = data(k2+1)
               data(k2) = data(k1)-tempr
               data(k2+1) = data(k1+1)-tempi
               data(k1) = data(k1)+tempr
               data(k1+1) = data(k1+1)+tempi
  147       continue
  148    continue
  149    mmax = np1
  150    if (mmax-ntwo/2) 151,174,174
  151    lmax = max0(np1tw,mmax/2)
         do 173 l=np1,lmax,np1tw
            m = l
            if (mmax-np1) 156,156,152
  152       theta = -twopi*float(l)/float(4*mmax)
            if (isign) 154,153,153
  153       theta = -theta
  154       wr = cos(theta)
            wi = sin(theta)
  155       w2r = wr*wr-wi*wi
            w2i = 2.*wr*wi
            w3r = w2r*wr-w2i*wi
            w3i = w2r*wi+w2i*wr
  156       do 169 i1=1,i1rng,2
               kmin = i1+ipar*m
               if (mmax-np1) 157,157,158
  157          kmin = i1
  158          kdif = ipar*mmax
  159          kstep = 4*kdif
               if (kstep-ntwo) 160,160,169
  160          do 168 k1=kmin,ntot,kstep
                  k2 = k1+kdif
                  k3 = k2+kdif
                  k4 = k3+kdif
                  if (mmax-np1) 161,161,164
  161             u1r = data(k1)+data(k2)
                  u1i = data(k1+1)+data(k2+1)
                  u2r = data(k3)+data(k4)
                  u2i = data(k3+1)+data(k4+1)
                  u3r = data(k1)-data(k2)
                  u3i = data(k1+1)-data(k2+1)
                  if (isign) 162,163,163
  162             u4r = data(k3+1)-data(k4+1)
                  u4i = data(k4)-data(k3)
                  go to 167
  163             u4r = data(k4+1)-data(k3+1)
                  u4i = data(k3)-data(k4)
                  go to 167
  164             t2r = w2r*data(k2)-w2i*data(k2+1)
                  t2i = w2r*data(k2+1)+w2i*data(k2)
                  t3r = wr*data(k3)-wi*data(k3+1)
                  t3i = wr*data(k3+1)+wi*data(k3)
                  t4r = w3r*data(k4)-w3i*data(k4+1)
                  t4i = w3r*data(k4+1)+w3i*data(k4)
                  u1r = data(k1)+t2r
                  u1i = data(k1+1)+t2i
                  u2r = t3r+t4r
                  u2i = t3i+t4i
                  u3r = data(k1)-t2r
                  u3i = data(k1+1)-t2i
                  if (isign) 165,166,166
  165             u4r = t3i-t4i
                  u4i = t4r-t3r
                  go to 167
  166             u4r = t4i-t3i
                  u4i = t3r-t4r
  167             data(k1) = u1r+u2r
                  data(k1+1) = u1i+u2i
                  data(k2) = u3r+u4r
                  data(k2+1) = u3i+u4i
                  data(k3) = u1r-u2r
                  data(k3+1) = u1i-u2i
                  data(k4) = u3r-u4r
                  data(k4+1) = u3i-u4i
  168          continue
               kdif = kstep
               kmin = 4*(kmin-i1)+i1
               go to 159
  169       continue
            m = m+lmax
            if (m-mmax) 170,170,173
  170       if (isign) 171,172,172
  171       tempr = wr
            wr = (wr+wi)*rthlf
            wi = (wi-tempr)*rthlf
            go to 155
  172       tempr = wr
            wr = (wr-wi)*rthlf
            wi = (tempr+wi)*rthlf
            go to 155
  173    continue
         ipar = 3-ipar
         mmax = mmax+mmax
         go to 150
c
c     main loop for factors not equal to two.  apply the twiddle factor
c     w=exp(isign*2*pi*sqrt(-1)*(j1-1)*(j2-j1)/(ifp1+ifp2)), then
c     perform a fourier transform of length ifact(if), making use of
c     conjugate symmetries.
c
  174    if (ntwo-np2) 175,201,201
  175    ifp1 = ntwo
         if = inon2
         np1hf = np1/2
  176    ifp2 = ifact(if)*ifp1
         j1min = np1+1
         if (j1min-ifp1) 177,177,184
  177    do 183 j1=j1min,ifp1,np1
            theta = -twopi*float(j1-1)/float(ifp2)
            if (isign) 179,178,178
  178       theta = -theta
  179       wstpr = cos(theta)
            wstpi = sin(theta)
            wr = wstpr
            wi = wstpi
            j2min = j1+ifp1
            j2max = j1+ifp2-ifp1
            do 182 j2=j2min,j2max,ifp1
               i1max = j2+i1rng-2
               do 181 i1=j2,i1max,2
                  do 180 j3=i1,ntot,ifp2
                     tempr = data(j3)
                     data(j3) = data(j3)*wr-data(j3+1)*wi
                     data(j3+1) = tempr*wi+data(j3+1)*wr
  180             continue
  181          continue
               tempr = wr
               wr = wr*wstpr-wi*wstpi
               wi = tempr*wstpi+wi*wstpr
  182       continue
  183    continue
  184    theta = -twopi/float(ifact(if))
         if (isign) 186,185,185
  185    theta = -theta
  186    wstpr = cos(theta)
         wstpi = sin(theta)
         j2rng = ifp1*(1+ifact(if)/2)
         do 200 i1=1,i1rng,2
            do 199 i3=i1,ntot,np2
               j2max = i3+j2rng-ifp1
               do 197 j2=i3,j2max,ifp1
                  j1max = j2+ifp1-np1
                  do 193 j1=j2,j1max,np1
                     j3max = j1+np2-ifp2
                     do 192 j3=j1,j3max,ifp2
                        jmin = j3-j2+i3
                        jmax = jmin+ifp2-ifp1
                        i = 1+(j3-i3)/np1hf
                        if (j2-i3) 187,187,189
  187                   sumr = 0.
                        sumi = 0.
                        do 188 j=jmin,jmax,ifp1
                           sumr = sumr+data(j)
                           sumi = sumi+data(j+1)
  188                   continue
                        work(i) = sumr
                        work(i+1) = sumi
                        go to 192
  189                   iconj = 1+(ifp2-2*j2+i3+j3)/np1hf
                        j = jmax
                        sumr = data(j)
                        sumi = data(j+1)
                        oldsr = 0.
                        oldsi = 0.
                        j = j-ifp1
  190                   tempr = sumr
                        tempi = sumi
                        sumr = twowr*sumr-oldsr+data(j)
                        sumi = twowr*sumi-oldsi+data(j+1)
                        oldsr = tempr
                        oldsi = tempi
                        j = j-ifp1
                        if (j-jmin) 191,191,190
  191                   tempr = wr*sumr-oldsr+data(j)
                        tempi = wi*sumi
                        work(i) = tempr-tempi
                        work(iconj) = tempr+tempi
                        tempr = wr*sumi-oldsi+data(j+1)
                        tempi = wi*sumr
                        work(i+1) = tempr+tempi
                        work(iconj+1) = tempr-tempi
  192                continue
  193             continue
                  if (j2-i3) 194,194,195
  194             wr = wstpr
                  wi = wstpi
                  go to 196
  195             tempr = wr
                  wr = wr*wstpr-wi*wstpi
                  wi = tempr*wstpi+wi*wstpr
  196             twowr = wr+wr
  197          continue
               i = 1
               i2max = i3+np2-np1
               do 198 i2=i3,i2max,np1
                  data(i2) = work(i)
                  data(i2+1) = work(i+1)
                  i = i+2
  198          continue
  199       continue
  200    continue
         if = if+1
         ifp1 = ifp2
         if (ifp1-np2) 176,201,201
c
c     complete a real transform in the 1st dimension, n even, by con-
c     jugate symmetries.
c
  201    go to (230,220,230,202),icase
  202    nhalf = n
         n = n+n
         theta = -twopi/float(n)
         if (isign) 204,203,203
  203    theta = -theta
  204    wstpr = cos(theta)
         wstpi = sin(theta)
         wr = wstpr
         wi = wstpi
         imin = 3
         jmin = 2*nhalf-1
         go to 207
  205    j = jmin
         do 206 i=imin,ntot,np2
            sumr = (data(i)+data(j))/2.
            sumi = (data(i+1)+data(j+1))/2.
            difr = (data(i)-data(j))/2.
            difi = (data(i+1)-data(j+1))/2.
            tempr = wr*sumi+wi*difr
            tempi = wi*sumi-wr*difr
            data(i) = sumr+tempr
            data(i+1) = difi+tempi
            data(j) = sumr-tempr
            data(j+1) = -difi+tempi
            j = j+np2
  206    continue
         imin = imin+2
         jmin = jmin-2
         tempr = wr
         wr = wr*wstpr-wi*wstpi
         wi = tempr*wstpi+wi*wstpr
  207    if (imin-jmin) 205,208,211
  208    if (isign) 209,211,211
  209    do 210 i=imin,ntot,np2
            data(i+1) = -data(i+1)
  210    continue
  211    np2 = np2+np2
         ntot = ntot+ntot
         j = ntot+1
         imax = ntot/2+1
  212    imin = imax-2*nhalf
         i = imin
         go to 214
  213    data(j) = data(i)
         data(j+1) = -data(i+1)
  214    i = i+2
         j = j-2
         if (i-imax) 213,215,215
  215    data(j) = data(imin)-data(imin+1)
         data(j+1) = 0.
         if (i-j) 217,219,219
  216    data(j) = data(i)
         data(j+1) = data(i+1)
  217    i = i-2
         j = j-2
         if (i-imin) 218,218,216
  218    data(j) = data(imin)+data(imin+1)
         data(j+1) = 0.
         imax = imin
         go to 212
  219    data(1) = data(1)+data(2)
         data(2) = 0.
         go to 230
c
c     complete a real transform for the 2nd or 3rd dimension by
c     conjugate symmetries.
c
  220    if (i1rng-np1) 221,230,230
  221    do 229 i3=1,ntot,np2
            i2max = i3+np2-np1
            do 228 i2=i3,i2max,np1
               imin = i2+i1rng
               imax = i2+np1-2
               jmax = 2*i3+np1-imin
               if (i2-i3) 223,223,222
  222          jmax = jmax+np2
  223          if (idim-2) 226,226,224
  224          j = jmax+np0
               do 225 i=imin,imax,2
                  data(i) = data(j)
                  data(i+1) = -data(j+1)
                  j = j-2
  225          continue
  226          j = jmax
               do 227 i=imin,imax,np0
                  data(i) = data(j)
                  data(i+1) = -data(j+1)
                  j = j-np0
  227          continue
  228       continue
  229    continue
c
c     end of loop on each dimension
c
  230    np0 = np1
         np1 = np2
         nprev = n
  231 continue
  232 return
c
c revision history---
c
c january 1978     deleted references to the  *cosy  cards and
c                  added revision history
c-----------------------------------------------------------------------
      end
