w3ft12.f

Go to the documentation of this file.

C> @file
C> @brief Fast fourier for 2.5 degree grid.
C> @author Joe Sela @date 1980-11-21
 
C> Fast fourier to compute 145 grid values at desired
C> latitude from 31 complex fourier coefficients. This subroutine
C> is special purpose for converting coefficients to a 2.5 degree
C> lat,lon grid.
C>
C> ### Program History Log:
C> Date | Programmer | Comment
C> -----|------------|--------
C> 1980-11-21 | Joe Sela | Initial.
C> 1984-06-21 | Ralph Jones | Change to ibm vs fortran.
C> 1993-04-12 | Ralph Jones | Change to cray cft77 fortran.
C>
C> @param[in] COEF 31 complex fourier coefficients.
C> @param[in] TRIGS 216 trig functions assumed precomputed by w3fa13() before
C> first call to w3ft12().
C> @param[in] WORK 144 real work space
C> @param[out] GRID 145 grid values, grid(1)=grid(145)
C>
C> @author Joe Sela @date 1980-11-21
       SUBROUTINE w3ft12(COEF,WORK,GRID,TRIGS)
       REAL            COEF( 62 )
       REAL            GRID(145)
       REAL            TRIGS(216)
       REAL            WORK(144)
C
       SAVE
C
       DATA  sin60/0.866025403784437/
C
      grid(1) = coef(1)
      grid(2) = coef(1)
      k = 147
      j = 143
      DO 100 i=3, 61 ,2
      temp = coef(i)*trigs(k+1) - coef(i+1)*trigs(k)
      grid(i) = coef(i) - temp
      grid(j) = coef(i) + temp
      temp = coef(i)*trigs(k) + coef(i+1)*trigs(k+1)
      grid(i+1) = temp - coef(i+1)
      grid(j+1) = temp + coef(i+1)
      k = k + 2
      j = j - 2
100   CONTINUE
      DO 110 i= 63 , 84
      grid(i) = 0.0
110   CONTINUE
C
      a0 = grid(1) + grid(73)
      a2 = grid(1) - grid(73)
      b0 = grid(2) + grid(74)
      b2 = grid(2) - grid(74)
      a1 = grid(37) + grid(109)
      a3 = grid(37) - grid(109)
      b1 = grid(38) + grid(110)
      b3 = grid(38) - grid(110)
      work(1) = a0 + a1
      work(5) = a0 - a1
      work(2) = b0 + b1
      work(6) = b0 - b1
      work(3) = a2 - b3
      work(7) = a2 + b3
      work(4) = b2 + a3
      work(8) = b2 - a3
      kb = 3
      kc = 5
      kd = 7
      j = 75
      k = 39
      l = 111
      m = 9
      DO 300 i=3,35,2
      a0 = grid(i) + grid(j)
      a2 = grid(i) - grid(j)
      b0 = grid(i+1) + grid(j+1)
      b2 = grid(i+1) - grid(j+1)
      a1 = grid(k) + grid(l)
      a3 = grid(k) - grid(l)
      b1 = grid(k+1) + grid(l+1)
      b3 = grid(k+1) - grid(l+1)
      work(m  ) = a0 + a1
      work(m+4) = a0 - a1
      work(m+1) = b0 + b1
      work(m+5) = b0 - b1
      work(m+2) = a2 - b3
      work(m+6) = a2 + b3
      work(m+3) = b2 + a3
      work(m+7) = b2 - a3
      temp = work(m+2)*trigs(kb) - work(m+3)*trigs(kb+1)
      work(m+3) = work(m+2)*trigs(kb+1) + work(m+3)*trigs(kb)
      work(m+2) = temp
      temp = work(m+4)*trigs(kc) - work(m+5)*trigs(kc+1)
      work(m+5) = work(m+4)*trigs(kc+1) + work(m+5)*trigs(kc)
      work(m+4) = temp
      temp = work(m+6)*trigs(kd) - work(m+7)*trigs(kd+1)
      work(m+7) = work(m+6)*trigs(kd+1) + work(m+7)*trigs(kd)
      work(m+6) = temp
      j = j + 2
      k = k + 2
      l = l + 2
      kb = kb + 2
      kc = kc + 4
      kd = kd + 6
      m = m + 8
300   CONTINUE
C
      i = 1
      j = 1
      k = 73
      DO 440 l=1,4
      grid(i) = work(j) + work(k)
      grid(i+8) = work(j) - work(k)
      grid(i+1) = work(j+1) + work(k+1)
      grid(i+9) = work(j+1) - work(k+1)
      i = i + 2
      j = j + 2
      k = k + 2
440   CONTINUE
      DO 500 kb=9,65,8
      i = i + 8
      DO 460 l=1,4
      grid(i) = work(j) + work(k)
      grid(i+8) = work(j) - work(k)
      grid(i+1) = work(j+1) + work(k+1)
      grid(i+9) = work(j+1) - work(k+1)
      temp = grid(i+8)*trigs(kb) - grid(i+9)*trigs(kb+1)
      grid(i+9) = grid(i+8)*trigs(kb+1) + grid(i+9)*trigs(kb)
      grid(i+8) = temp
      i = i + 2
      j = j + 2
      k = k + 2
460   CONTINUE
500   CONTINUE
C
      i = 1
      l = 1
      kc = 1
      j = 49
      k = 97
      m = 17
      n = 33
      DO 660 ll=1,8
      a1 = grid(j) + grid(k)
      a3 = sin60*(grid(j)-grid(k))
      b1 = grid(j+1) + grid(k+1)
      b3 = sin60*(grid(j+1)-grid(k+1))
      work(l) = grid(i) + a1
      a2 = grid(i) - 0.5*a1
      work(l+1) = grid(i+1) + b1
      b2 = grid(i+1) - 0.5*b1
      work(n) = a2 + b3
      work(m) = a2 - b3
      work(m+1) = b2 + a3
      work(n+1) = b2 - a3
      i = i + 2
      j = j + 2
      k = k + 2
      l = l + 2
      m = m + 2
      n = n + 2
660   CONTINUE
      DO 700 kb=17,33,16
      l = l + 32
      m = m + 32
      n = n + 32
      kc = kc + 32
      DO 680 ll=1,8
      a1 = grid(j) + grid(k)
      a3 = sin60*(grid(j)-grid(k))
      b1 = grid(j+1) + grid(k+1)
      b3 = sin60*(grid(j+1)-grid(k+1))
      work(l) = grid(i) + a1
      a2 = grid(i) - 0.5*a1
      work(l+1) = grid(i+1) + b1
      b2 = grid(i+1) - 0.5*b1
      work(n) = a2 + b3
      work(m) = a2 - b3
      work(m+1) = b2 + a3
      work(n+1) = b2 - a3
      temp = work(m)*trigs(kb) - work(m+1)*trigs(kb+1)
      work(m+1) = work(m)*trigs(kb+1) + work(m+1)*trigs(kb)
      work(m) = temp
      temp = work(n)*trigs(kc) - work(n+1)*trigs(kc+1)
      work(n+1) = work(n)*trigs(kc+1) + work(n+1)*trigs(kc)
      work(n) = temp
      i = i + 2
      j = j + 2
      k = k + 2
      l = l + 2
      m = m + 2
      n = n + 2
680   CONTINUE
700   CONTINUE
C
      j = 49
      k = 97
      l = 144
      m = 96
      n = 48
      DO 900 i=1,47,2
      a1 = work(j) + work(k)
      a3 = sin60 * (work(j)-work(k))
      b3 = sin60 * (work(j+1)-work(k+1))
      b1 = work(j+1) + work(k+1)
      grid(l+1) = work(i) + a1
      a2 = work(i) - 0.5*a1
      b2 = work(i+1) - 0.5*b1
      grid(l) = work(i+1) + b1
      grid(n+1) = a2 + b3
      grid(m+1) = a2 - b3
      grid(m) = b2 + a3
      grid(n) = b2 - a3
      j = j + 2
      k = k + 2
      l = l - 2
      m = m - 2
      n = n - 2
900   CONTINUE
      grid(1) = grid(145)
C
      RETURN
      END