MODULE FFT72 6,1
!@sum FFT72 calculates the Fast Fourier Transform
!@auth Gary Russell
!@ver 1.0 (for KM=72)
USE CONSTANT
, only : twopi,rt2,rt3
IMPLICIT NONE
SAVE
!@param KM length of input array; to change KM => rewrite module !!!
INTEGER, PARAMETER :: KM=72
REAL*8, PARAMETER :: BYKM=1d0/KM !@param BYKM 1/KM
REAL*8, PARAMETER :: BYKMH=2d0/KM !@param BYKMH 1/(KM/2)
REAL*8, PARAMETER :: BYKM2=1d0/(2*KM)!@param BYKM2 1/(2*KM)
!@var C,S cos/sin evaluated on grid points
REAL*8 :: C(0:KM),S(0:KM)
!@var CH,SH cos/sin evaluated on half points
REAL*8 :: CH(KM/2-1),SH(KM/2-1)
!@var C240,C241,S241,C8,S8 intermediate sums for FFT
!@var C41,C42,C43,C44,S41,S42,S43,S44 intermediate sums for FFT
!@var C21,C22,S21,S22 intermediate sums for FFT
REAL*8 :: C240(24), C241(24), S241(24)
REAL*8 :: C8(8,0:4),S8(8,4)
REAL*8 :: C41(0:9),C42(0:9),C43(0:9),C44(0:9),
* S41( 9),S42( 9),S43( 9),S44( 9)
REAL*8 :: C21(0:18),C22(0:18),S21(0:18),S22(0:18)
COMMON /FFTCOM/ C240,C241,S241,C8,S8,C41,C42,C43,C44,
* S41,S42,S43,S44,C21,C22,S21,S22
!$OMP THREADPRIVATE(/FFTCOM/)
END MODULE FFT72
C****
SUBROUTINE FFT0 (IM) 1,2
!@sum FFT0 initializes sines and cosines used by FFT routines.
!@auth Gary Russell
!@ver 1.0
USE FFT72
IMPLICIT NONE
INTEGER IQ,N !@var IQ,N loop variables
INTEGER, INTENT(IN) :: IM !@var IM size of arrays (must=KM)
C
IF(IM.NE.KM) GO TO 100
DO N=0,KM/4
C(N) = COS(TWOPI*N/REAL(KM,KIND=8))
S(KM/4 -N) = C(N)
C(KM/2 -N) = -C(N)
S(KM/4 +N) = C(N)
C(KM/2 +N) = -C(N)
S(3*KM/4 -N) = -C(N)
C(KM -N) = C(N)
S(3*KM/4 +N) = -C(N)
END DO
DO N=1,KM/2-1
CH(N) = COS(TWOPI*N/REAL(2*KM,KIND=8))
SH(N) = SIN(TWOPI*N/REAL(2*KM,KIND=8))
END DO
RETURN
100 WRITE (6,*) ' This version of FFT is for ',KM,'. IM =',IM
call stop_model('stopped in FFT72.f',255)
END SUBROUTINE FFT0
C****
SUBROUTINE FFT (F,A,B) 10,2
!@sum FFT calculates fast fourier transform of an input array F
!@auth Gary Russell
!@ver 1.0
!@calls DOCALC
C****
C**** FFT calculates a fast fourier transform of the input array F,
C**** producing the cosine and sine coefficients in the output arrays
C**** A and B. F is dimensioned by KM; A and B are dimensioned
C**** 0:KM/2. Upon entering FFT, the total energy is:
C**** .5*sum(F(K)**2)
C**** with the sum being taken over all K from 1 to KM. The
C**** Fourier coefficients are defined by:
C**** A(N)+i*B(N) = sum(F(K)*exp(-2*PI*i*N*K/KM))/KMH
C**** with the sum being taken over all K from 1 to KM. KMH = KM
C**** when N = 0 or KM/2, and KMH = KM/2 otherwise. In the
C**** program's notation, CPQN means:
C**** CPQN = sum(F(K)*cos(2*PI*N*K/KM))
C**** with the sum being taken over all K = Q mod(P). SPQN has a
C**** similar definition, but cos is replaced by sin. The notation
C**** A=10, B=11, etc. is used for P, Q and N. The same total
C**** energy can be calculated from the spectral coefficients as:
C**** .5*sum(A(N)**2+B(N)**2)*KMH
C**** with the sum being taken over all wave numbers N from 0 to KM/2.
C****
USE FFT72
IMPLICIT NONE
REAL*8, INTENT(IN) :: F(KM) !@var F input gridpoint array
REAL*8, INTENT(OUT) :: A(0:KM/2) !@var A fourier coeffs. (cos)
REAL*8, INTENT(OUT) :: B(0:KM/2) !@var B fourier coeffs. (sin)
INTEGER IQ,N !@var IQ,N loop variables
CALL DOCALC(F)
C**** Calculate final coefficients of fourier expansion
A(0) = (C22(0)+C21(0))*BYKM
B(0) = 0.
DO N=1,KM/4
A(N) = (C22(N)+C21(N))*BYKMH
B(N) = (S22(N)+S21(N))*BYKMH
END DO
DO N=1,KM/4-1
A(KM/2-N) = (C22(N)-C21(N))*BYKMH
B(KM/2-N) = (S21(N)-S22(N))*BYKMH
END DO
A(KM/2) = (C22(0)-C21(0))*BYKM
B(KM/2) = 0.
C****
RETURN
END SUBROUTINE FFT
C****
SUBROUTINE FFTI (A,B,F) 1,1
!@sum FFTI performs an inverse fast fourier transform
!@auth Gary Russell
!@ver 1.0 (KM=72)
USE FFT72
IMPLICIT NONE
REAL*8, INTENT(OUT) :: F(KM) !@var F output gridpoint array
REAL*8, INTENT(IN) :: A(0:KM/2) !@var A fourier coeffs. (cos)
REAL*8, INTENT(IN) :: B(0:KM/2) !@var B fourier coeffs. (sin)
INTEGER IQ,N !@var IQ,N loop variables
C**** These have been divided by 18
C22(0) = ( C21(0) = ( DO N=1,18
C22(N) = A(N)+A(36-N)
C21(N) = A(N)-A(36-N)
S22(N) = B(N)-B(36-N)
S21(N) = B(N)+B(36-N)
END DO
C**** Now multiply by 2, so FACTOR = 2/18 = 9
DO N=0,9
C42(N)=C22(N)-C22(18-N)
C44(N)=C22(N)+C22(18-N)
END DO
DO N=0,8
C41(N)=C21(N)+S21(18-N)
C43(N)=C21(N)-S21(18-N)
END DO
DO N=1,9
S42(N)=S22(N)+S22(18-N)
S44(N)=S22(N)-S22(18-N)
S41(N)=S21(N)+C21(18-N)
S43(N)=S21(N)-C21(18-N)
END DO
C**** Multiply by 2 again, so FACTOR = 2*2/18 = 2/9
DO N=0,4
C8(8,N)=C44(N)+C44(9-N)
C8(4,N)=C44(N)-C44(9-N)
C8(2,N)=C42(N)+S42(9-N)
C8(6,N)=C42(N)-S42(9-N)
END DO
DO N=1,4
S8(2,N)=S42(N)+C42(9-N)
S8(6,N)=S42(N)-C42(9-N)
S8(8,N)=S44(N)-S44(9-N)
S8(4,N)=S44(N)+S44(9-N)
S8(5,N)=S41(N)+(S41(9-N)-C41(9-N))*S(9)
S8(1,N)=S41(N)-(S41(9-N)-C41(9-N))*S(9)
S8(3,N)=S43(N)+(S43(9-N)+C43(9-N))*S(9)
S8(7,N)=S43(N)-(S43(9-N)+C43(9-N))*S(9)
C8(5,N)=C41(N)-(S41(9-N)+C41(9-N))*S(9)
C8(1,N)=C41(N)+(S41(9-N)+C41(9-N))*S(9)
C8(3,N)=C43(N)+(S43(9-N)-C43(9-N))*S(9)
C8(7,N)=C43(N)-(S43(9-N)-C43(9-N))*S(9)
END DO
C8(5,0)=C41(0)-S41(9)*RT2
C8(1,0)=C41(0)+S41(9)*RT2
C8(3,0)=C43(0)+S43(9)*RT2
C8(7,0)=C43(0)-S43(9)*RT2
C**** Multiply by 3, so FACTOR = 3*2*2/18 = 3*2/9 = 2/3
C240(24)= C8(8,0)+C8(8,3)*2.
C240(16)= C8(8,0)-C8(8,3)-S8(8,3)*RT3
C240(8) = C8(8,0)-C8(8,3)+S8(8,3)*RT3
C240(18)= C8(2,0)-S8(2,3)*2.
C240(2) = C8(2,0)+S8(2,3)+C8(2,3)*RT3
C240(10)= C8(2,0)+S8(2,3)-C8(2,3)*RT3
C240(12)= C8(4,0)-C8(4,3)*2.
C240(4) = C8(4,0)+C8(4,3)+S8(4,3)*RT3
C240(20)= C8(4,0)+C8(4,3)-S8(4,3)*RT3
C240(6) = C8(6,0)+S8(6,3)*2.
C240(14)= C8(6,0)-S8(6,3)-C8(6,3)*RT3
C240(22)= C8(6,0)-S8(6,3)+C8(6,3)*RT3
C240(9) = C8(1,0)-(C8(1,3)-S8(1,3))*RT2
C240(1) = C8(1,0)+C8(1,3)*2.*S(15)+S8(1,3)*2.*S(3)
C240(17)= C8(1,0)-C8(1,3)*2.*S(3)-S8(1,3)*2.*S(15)
C240(3) = C8(3,0)+(C8(3,3)+S8(3,3))*RT2
C240(11)= C8(3,0)-C8(3,3)*2.*S(15)+S8(3,3)*2.*S(3)
C240(19)= C8(3,0)+C8(3,3)*2.*S(3)-S8(3,3)*2.*S(15)
C240(21)= C8(5,0)+(C8(5,3)-S8(5,3))*RT2
C240(13)= C8(5,0)-C8(5,3)*2.*S(15)-S8(5,3)*2.*S(3)
C240(5) = C8(5,0)+C8(5,3)*2.*S(3)+S8(5,3)*2.*S(15)
C240(15)= C8(7,0)-(C8(7,3)+S8(7,3))*RT2
C240(23)= C8(7,0)+C8(7,3)*2.*S(15)-S8(7,3)*2.*S(3)
C240(7) = C8(7,0)-C8(7,3)*2.*S(3)+S8(7,3)*2.*S(15)
C241(24)= C8(8,2)+C8(8,4)+C8(8,1)
S241(24)= -S8(8,2)+S8(8,4)+S8(8,1)
C241(8) =-(C8(8,2)+C8(8,4))*S(6)+(S8(8,2)+S8(8,4))*S(12)+C8(8,1)
C241(16)=-(C8(8,2)+C8(8,4))*S(6)-(S8(8,2)+S8(8,4))*S(12)+C8(8,1)
S241(8) = (C8(8,2)-C8(8,4))*S(12)+(S8(8,2)-S8(8,4))*S(6)+S8(8,1)
S241(16)=-(C8(8,2)-C8(8,4))*S(12)+(S8(8,2)-S8(8,4))*S(6)+S8(8,1)
C241(18)= -S8(2,2)-S8(2,4)+C8(2,1)
S241(18)= -C8(2,2)+C8(2,4)+S8(2,1)
C241(2) = (C8(2,2)+C8(2,4))*S(12)+(S8(2,2)+S8(2,4))*S(6)+C8(2,1)
C241(10)=-(C8(2,2)+C8(2,4))*S(12)+(S8(2,2)+S8(2,4))*S(6)+C8(2,1)
S241(2) = (C8(2,2)-C8(2,4))*S(6)-(S8(2,2)-S8(2,4))*S(12)+S8(2,1)
S241(10)= (C8(2,2)-C8(2,4))*S(6)+(S8(2,2)-S8(2,4))*S(12)+S8(2,1)
S241(12)= S8(4,2)-S8(4,4)+S8(4,1)
C241(12)= -C8(4,2)-C8(4,4)+C8(4,1)
S241(4) = (C8(4,2)-C8(4,4))*S(12)-(S8(4,2)-S8(4,4))*S(6)+S8(4,1)
S241(20)=-(C8(4,2)-C8(4,4))*S(12)-(S8(4,2)-S8(4,4))*S(6)+S8(4,1)
C241(4) = (C8(4,2)+C8(4,4))*S(6)+(S8(4,2)+S8(4,4))*S(12)+C8(4,1)
C241(20)= (C8(4,2)+C8(4,4))*S(6)-(S8(4,2)+S8(4,4))*S(12)+C8(4,1)
C241(6) = S8(6,2)+S8(6,4)+C8(6,1)
S241(6) = C8(6,2)-C8(6,4)+S8(6,1)
C241(22)= (C8(6,2)+C8(6,4))*S(12)-(S8(6,2)+S8(6,4))*S(6)+C8(6,1)
C241(14)=-(C8(6,2)+C8(6,4))*S(12)-(S8(6,2)+S8(6,4))*S(6)+C8(6,1)
S241(22)=-(C8(6,2)-C8(6,4))*S(6)-(S8(6,2)-S8(6,4))*S(12)+S8(6,1)
S241(14)=-(C8(6,2)-C8(6,4))*S(6)+(S8(6,2)-S8(6,4))*S(12)+S8(6,1)
C241(9) =(-C8(1,2)-C8(1,4)+S8(1,2)+S8(1,4))*S(9)+C8(1,1)
S241(9) = (C8(1,2)-C8(1,4)+S8(1,2)-S8(1,4))*S(9)+S8(1,1)
C241(1) = (C8(1,2)+C8(1,4))*S(15)+(S8(1,2)+S8(1,4))*S(3)+C8(1,1)
C241(17)=-(C8(1,2)+C8(1,4))*S(3)-(S8(1,2)+S8(1,4))*S(15)+C8(1,1)
S241(1) = (C8(1,2)-C8(1,4))*S(3)-(S8(1,2)-S8(1,4))*S(15)+S8(1,1)
S241(17)=-(C8(1,2)-C8(1,4))*S(15)+(S8(1,2)-S8(1,4))*S(3)+S8(1,1)
C241(3) = (C8(3,2)+C8(3,4)+S8(3,2)+S8(3,4))*S(9)+C8(3,1)
S241(3) = (C8(3,2)-C8(3,4)-S8(3,2)+S8(3,4))*S(9)+S8(3,1)
C241(11)=-(C8(3,2)+C8(3,4))*S(15)+(S8(3,2)+S8(3,4))*S(3)+C8(3,1)
S241(19)=-(C8(3,2)-C8(3,4))*S(15)-(S8(3,2)-S8(3,4))*S(3)+S8(3,1)
C241(19)= (C8(3,2)+C8(3,4))*S(3)-(S8(3,2)+S8(3,4))*S(15)+C8(3,1)
S241(11)= (C8(3,2)-C8(3,4))*S(3)+(S8(3,2)-S8(3,4))*S(15)+S8(3,1)
C241(21)= (C8(5,2)+C8(5,4)-S8(5,2)-S8(5,4))*S(9)+C8(5,1)
S241(21)=(-C8(5,2)+C8(5,4)-S8(5,2)+S8(5,4))*S(9)+S8(5,1)
C241(13)=-(C8(5,2)+C8(5,4))*S(15)-(S8(5,2)+S8(5,4))*S(3)+C8(5,1)
S241(13)=-(C8(5,2)-C8(5,4))*S(3)+(S8(5,2)-S8(5,4))*S(15)+S8(5,1)
C241(5) = (C8(5,2)+C8(5,4))*S(3)+(S8(5,2)+S8(5,4))*S(15)+C8(5,1)
S241(5) = (C8(5,2)-C8(5,4))*S(15)-(S8(5,2)-S8(5,4))*S(3)+S8(5,1)
C241(15)=-(C8(7,2)+C8(7,4)+S8(7,2)+S8(7,4))*S(9)+C8(7,1)
S241(15)=(-C8(7,2)+C8(7,4)+S8(7,2)-S8(7,4))*S(9)+S8(7,1)
C241(23)= (C8(7,2)+C8(7,4))*S(15)-(S8(7,2)+S8(7,4))*S(3)+C8(7,1)
C241(7) =-(C8(7,2)+C8(7,4))*S(3)+(S8(7,2)+S8(7,4))*S(15)+C8(7,1)
S241(23)=-(C8(7,2)-C8(7,4))*S(3)-(S8(7,2)-S8(7,4))*S(15)+S8(7,1)
S241(7) = (C8(7,2)-C8(7,4))*S(15)+(S8(7,2)-S8(7,4))*S(3)+S8(7,1)
C**** Multiply by FACTOR = 3/2
DO IQ=1,24
F(IQ+24) = (S241(IQ)*S(IQ+24)+C241(IQ)*C(IQ+24))+C240(IQ)*.5
F(IQ) = (S241(IQ)*C(IQ+48)-C241(IQ)*S(IQ+48))*RT3+F(IQ+24)
F(IQ+48) =-(S241(IQ)*C(IQ) -C241(IQ)*S(IQ)) *RT3+F(IQ+24)
END DO
C****
RETURN
END SUBROUTINE FFTI
C****
SUBROUTINE FFTE (F,E) 4,2
!@sum FFTE calcs. the spectral energy E from input gridpoint values F
!@auth Gary Russell
!@ver 1.0
!@calls DOCALC
USE FFT72
IMPLICIT NONE
REAL*8, INTENT(IN) :: F(KM) !@var F input gridpoint array
REAL*8, INTENT(OUT) :: E(0:KM/2) !@var E spectral energy
INTEGER IQ,N !@var IQ,N loop variables
CALL DOCALC(F)
C****
E(0) = (C22(0)+C21(0))*(C22(0)+C21(0))*BYKM2
DO N=1,KM/4
E(N) = ((C22(N)+C21(N))*(C22(N)+C21(N)) +
* (S22(N)+S21(N))*(S22(N)+S21(N)))*BYKM
END DO
DO N=1,KM/4-1
E(KM/2-N) = ((C22(N)-C21(N))*(C22(N)-C21(N)) +
* (S21(N)-S22(N))*(S21(N)-S22(N)))*BYKM
END DO
E(KM/2) = (C22(0)-C21(0))*(C22(0)-C21(0))*BYKM2
C****
RETURN
END SUBROUTINE FFTE
C****
SUBROUTINE FFT2 (F,A,B),2
!@sum FFT2 calculates Fourier coefficients on secondary grid
!@auth Gary Russell
!@ver 1.0
!@calls FFT
C****
C**** Values in F(K) are located a half space right of those in FFT.
C****
USE FFT72
IMPLICIT NONE
REAL*8, INTENT(IN) :: F(KM) !@var input array
REAL*8, INTENT(OUT) :: A(0:KM/2) !@var fourier coefficients (cos)
REAL*8, INTENT(OUT) :: B(0:KM/2) !@var fourier coefficients (sin)
REAL*8 ANEW !@var ANEW dummy variable
INTEGER N !@var N loop variable
CALL FFT (F,A,B)
DO N=1,KM/2-1
ANEW = A(N)*CH(N)-B(N)*SH(N)
B(N) = A(N)*SH(N)+B(N)*CH(N)
A(N) = ANEW
END DO
B(KM/2) = A(KM/2)
A(KM/2) = 0.
RETURN
END SUBROUTINE FFT2
C****
SUBROUTINE DOCALC(F) 2,1
!@sum DOCALC calculate intermediate expressions for FFT
!@auth Gary Russell
!@ver 1.0 (KM=72)
USE FFT72
IMPLICIT NONE
REAL*8, INTENT(IN) :: F(KM) !@var F input grid point array
INTEGER IQ,N !@var IQ,N loop variables
C**** Calculate expressions summed by increments of 24
DO IQ=1,24
C240(IQ) = F(IQ) + F(IQ+24) + F(IQ+48)
C241(IQ) = F(IQ)*C(IQ) + F(IQ+24)*C(IQ+24) + F(IQ+48)*C(IQ+48)
S241(IQ) = F(IQ)*S(IQ) + F(IQ+24)*S(IQ+24) + F(IQ+48)*S(IQ+48)
END DO
C**** Calculate expressions summed by increments of 8
DO IQ = 1,8
C8(IQ,0) = C240(IQ)+C240(IQ+8)+C240(IQ+16)
C8(IQ,1) = C241(IQ)+C241(IQ+8)+C241(IQ+16)
S8(IQ,1) = S241(IQ)+S241(IQ+8)+S241(IQ+16)
END DO
C8(1,2)= (C241(1)-S241(17))*S(15)+(S241(9)-C241(9))*S(9)+(S241(1)
* -C241(17))*S(3)
C8(2,2)= (C241(2)-C241(10))*S(12)+(S241(2)+S241(10))*S(6)-S241(18)
C8(3,2)= (C241(3)+S241(3))*S(9)-(C241(11)+S241(19))*S(15)
* +(C241(19)+S241(11))*S(3)
C8(4,2)= (C241(4)+C241(20))*S(6)+(S241(4)-S241(20))*S(12)-C241(12)
C8(5,2)= (C241(5)-S241(13))*S(3)+(C241(21)-S241(21))*S(9)
* -(C241(13)-S241(5))*S(15)
C8(6,2)= (C241(22)-C241(14))*S(12)-(S241(22)+S241(14))*S(6)
* +S241(6)
C8(7,2)= (C241(23)+S241(7))*S(15)-(S241(15)+C241(15))*S(9)
* -(S241(23)+C241(7))*S(3)
C8(8,2)= C241(24)-(C241(8)+C241(16))*S(6)+(S241(8)-S241(16))*S(12)
C8(1,3)= C240(1)*S(15)-C240(9)*S(9)-C240(17)*S(3)
C8(2,3)= (C240(2)-C240(10))*S(12)
C8(3,3)= C240(3)*S(9)-C240(11)*S(15)+C240(19)*S(3)
C8(4,3)= (C240(4)+C240(20))*S(6)-C240(12)
C8(5,3)= C240(5)*S(3)-C240(13)*S(15)+C240(21)*S(9)
C8(6,3)= (C240(22)-C240(14))*S(12)
C8(7,3)= C240(23)*S(15)-C240(15)*S(9)-C240(7)*S(3)
C8(8,3)= C240(24)-(C240(8)+C240(16))*S(6)
C8(1,4)= (C241(1)+S241(17))*S(15)-(S241(9)+C241(9))*S(9)-(S241(1)
* +C241(17))*S(3)
C8(2,4)= (C241(2)-C241(10))*S(12)-(S241(2)+S241(10))*S(6)+S241(18)
C8(3,4)= (C241(3)-S241(3))*S(9)-(C241(11)-S241(19))*S(15)
* +(C241(19)-S241(11))*S(3)
C8(4,4)= (C241(4)+C241(20))*S(6)-(S241(4)-S241(20))*S(12)-C241(12)
C8(5,4)= (C241(5)+S241(13))*S(3)+(C241(21)+S241(21))*S(9)
* -(C241(13)+S241(5))*S(15)
C8(6,4)= (C241(22)-C241(14))*S(12)+(S241(22)+S241(14))*S(6)
* -S241(6)
C8(7,4)= (C241(23)-S241(7))*S(15)+(S241(15)-C241(15))*S(9)
* +(S241(23)-C241(7))*S(3)
C8(8,4)= C241(24)-(C241(8)+C241(16))*S(6)-(S241(8)-S241(16))*S(12)
S8(1,2)= (C241(1)+S241(17))*S(3)+(C241(9)+S241(9))*S(9)-(C241(17)
* +S241(1))*S(15)
S8(2,2)= (C241(10)+C241(2))*S(6)+(S241(10)-S241(2))*S(12)-C241(18)
S8(3,2)= (C241(3)-S241(3))*S(9)+(S241(11)-C241(19))*S(15)
* +(C241(11)-S241(19))*S(3)
S8(4,2)= (C241(4)-C241(20))*S(12)-(S241(4)+S241(20))*S(6)+S241(12)
S8(5,2)= (C241(5)+S241(13))*S(15)-(S241(21)+C241(21))*S(9)
* -(C241(13)+S241(5))*S(3)
S8(6,2)= C241(6)-(C241(14)+C241(22))*S(6)+(S241(14)-S241(22))
* *S(12)
S8(7,2)= (C241(7)-S241(23))*S(15)+(S241(15)-C241(15))*S(9)
* -(C241(23)-S241(7))*S(3)
S8(8,2)= (C241(8)-C241(16))*S(12)+(S241(8)+S241(16))*S(6)-S241(24)
S8(1,3)= C240(1)*S(3)+C240(9)*S(9)-C240(17)*S(15)
S8(2,3)= (C240(10)+C240(2))*S(6)-C240(18)
S8(3,3)= C240(3)*S(9)-C240(19)*S(15)+C240(11)*S(3)
S8(4,3)= (C240(4)-C240(20))*S(12)
S8(5,3)= C240(5)*S(15)-C240(21)*S(9)-C240(13)*S(3)
S8(6,3)= C240(6)-(C240(14)+C240(22))*S(6)
S8(7,3)= C240(7)*S(15)-C240(15)*S(9)-C240(23)*S(3)
S8(8,3)= (C240(8)-C240(16))*S(12)
S8(1,4)= (C241(1)-S241(17))*S(3)+(C241(9)-S241(9))*S(9)-(C241(17)
* -S241(1))*S(15)
S8(2,4)= (C241(10)+C241(2))*S(6)-(S241(10)-S241(2))*S(12)-C241(18)
S8(3,4)= (C241(3)+S241(3))*S(9)-(S241(11)+C241(19))*S(15)
* +(C241(11)+S241(19))*S(3)
S8(4,4)= (C241(4)-C241(20))*S(12)+(S241(4)+S241(20))*S(6)-S241(12)
S8(5,4)= (C241(5)-S241(13))*S(15)+(S241(21)-C241(21))*S(9)
* -(C241(13)-S241(5))*S(3)
S8(6,4)= C241(6)-(C241(14)+C241(22))*S(6)-(S241(14)-S241(22))
* *S(12)
S8(7,4)= (C241(7)+S241(23))*S(15)-(S241(15)+C241(15))*S(9)
* -(C241(23)+S241(7))*S(3)
S8(8,4)= (C241(8)-C241(16))*S(12)-(S241(8)+S241(16))*S(6)+S241(24)
C**** Calculate expressions summed by increments of 4
DO N=0,4
C41(N)=C8(1,N)+C8(5,N)
C42(N)=C8(2,N)+C8(6,N)
C43(N)=C8(3,N)+C8(7,N)
C44(N)=C8(8,N)+C8(4,N)
END DO
DO N=1,4
S41(N)=S8(1,N)+S8(5,N)
S42(N)=S8(2,N)+S8(6,N)
S43(N)=S8(3,N)+S8(7,N)
S44(N)=S8(8,N)+S8(4,N)
C41(9-N) =(C8(1,N)-C8(5,N)+S8(1,N)-S8(5,N))*C(9)
C42(9-N) = S8(2,N)-S8(6,N)
C43(9-N) =(C8(7,N)-C8(3,N)+S8(3,N)-S8(7,N))*C(9)
C44(9-N) = C8(8,N)-C8(4,N)
S41(9-N) =(C8(1,N)-C8(5,N)+S8(5,N)-S8(1,N))*C(9)
S42(9-N) = C8(2,N)-C8(6,N)
S43(9-N) =(C8(3,N)-C8(7,N)+S8(3,N)-S8(7,N))*C(9)
S44(9-N) = S8(4,N)-S8(8,N)
END DO
C41(9) =(C8(1,0)-C8(5,0))*C(9)
C42(9) = 0.
C43(9) =(C8(7,0)-C8(3,0))*C(9)
C44(9) = C8(8,0)-C8(4,0)
S41(9) =(C8(1,0)-C8(5,0))*C(9)
S42(9) = C8(2,0)-C8(6,0)
S43(9) =(C8(3,0)-C8(7,0))*C(9)
S44(9) = 0.
C**** Calculate expressions summed by increments of 2
DO N=0,9
C21(N) = C41(N)+C43(N)
C22(N) = C44(N)+C42(N)
END DO
DO N=1,8
S21(N) = S41(N)+S43(N)
S22(N) = S44(N)+S42(N)
END DO
DO N=1,9
C21(18-N) = S41(N)-S43(N)
C22(18-N) = C44(N)-C42(N)
S21(18-N) = C41(N)-C43(N)
S22(18-N) = S42(N)-S44(N)
END DO
C21(18) = 0.
C22(18) = C44(0)-C42(0)
S21(18) = C41(0)-C43(0)
S22(18) = 0.
C****
RETURN
END SUBROUTINE DOCALC