Fortran 90 subroutine CYCLOSTATIONARY_WD
This subroutine computes cyclostationary EOFs (CSEOFs), their principal component time
series, and the percent variance explained by each CSEOF. Portions of the code are proto-
typed after Dr.
Kwang-Yul Kim
's cseofx.f for cyclosationary analysis, which was obtained from
ftp://pacific.met.fsu.edu/pub/kkim/eigen/cyclo_eof/ at Florida State University.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! File: SUBR_CYCLOSTATIONARY_WD.f90
! Subprogram: CYCLOSTATIONARY_WD
! Author: David Stepaniak
! Affiliation: National Center for Atmospheric Research
! Climate and Global Dynamics Division
! Climate Analysis Section
! Date Initiated: 19 November 2002
! Last Modified: 25 November 2002
! Dependencies: WEIGHT and INTEGR, which are included in this subprogram -- see the `CONTAINS'
! statement. SSYEV (from LAPACK library).
! Required Lib: LAPACK library (on dataproc use -lcomplib.sgimath; more info
! about dataproc can be found at http://www.scd.ucar.edu/computers/dataproc/)
! Acknowledgments: David Stepaniak acknowledges the gracious assistance of Dr. Kim in answering
! several questions regarding cseofx.f, provided as exchanges of e-mails in
! November of 2002.
!
! This subroutine computes cyclostationary EOFs (CSEOFs), their principal component time
! series, and the percent variance explained by each CSEOF. Portions of the code are proto-
! typed after Dr. Kwang-Yul Kim's cseofx.f for cyclosationary analysis, which was obtained from
! ftp://pacific.met.fsu.edu/pub/kkim/eigen/cyclo_eof/ at Florida State University.
! However, the eigenvectors and eigenvalues of the temporal covariance matrix, a real symmetric
! matrix, are obtained with SSYEV of the public domain LAPACK library, and not with Dr. Kim's
! JACOBI routine (which is from Numerical Recipes). On the other hand, the dependencies WEIGHT
! and INTEGR (see calls to these routines below) are Dr Kim's with the exception that they have
! been updated to the Fortran 90 standard.
!
! In this routine,
! (number of stations) * (nested period) >= (number of temporal samples)
! and thus, as previously mentioned, utilizes cseofx.f as a prototype and the temporal
! covariance matrix contained therein. The temporal covariance matrix has dimensions
! ntime x ntime. The eigenvectors of the covariance matrix, obtained with SSYEV of
! the LAPACK library, are in fact the principal component time series. (The eigenvectors,
! upon return from SSYEV, are stored as columns and must be sorted in reverse order
! so that the leading mode is the first, rather than the last, column.) The principal
! component time series are then scaled by the square root of the number of temporal
! samples, and normalized by their standard deviation (see below).
!
! The basic methodology is as follows. Given a harmonizable time series X(r,t), i.e. a
! time series with an inherent or hypothesized periodicity (nested period) d, where
! r is the spatial dimension and t the temporal dimension we have (from K.-Y. Kim),
!
! "X(r,t) = Sum_k(0,T-1) [ a_k(r,t) exp(tpi ikt/d) ],
! where
! a_k(r,t) = Int [ w(t-s) X(r,s) exp(-tpi iks/d) ds ]
! and
! w(t) = sin(pi t/d) / (pi t),
!
! it then follows that the eigenfunctions are
!
! f_nm(r,t) = exp(2pi*i*n*t/d) * g_m(r,t),
!
! where g_m(r,t) is an eigenfunction of the covariance matrix
!
! C(k,l) = < a_k(r,t) a_l(r,t) >."
!
! This subroutine has been extensively tested and compared with results published by Dr. Kim
! regarding the annual cycle and (quasi-)biennial mode of tropical Pacific SST.
! Annual Cycle:
!
! Kim, K.-Y., and C. Chung, 2001, On the evolution of the annual cycle in the tropical Pacific.
! J. Clim., 14, 991-994.
! Biennial Mode:
! Kim, K.-Y., 2002, Investigation of ENSO variability using cyclostationary EOFs of observational
! data. J. of Met. and Atm. Phys., 81, 149-168.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE CYCLOSTATIONARY_WD(ntime,nstation,nestedPer,nmodes, &
TSER,msgOrig,CSEOF,cseofScale,PCTS,PCVAR)
IMPLICIT NONE
INTEGER, INTENT(IN) :: ntime
! ntime is the number of temporal
! samples.
INTEGER, INTENT(IN) :: nstation
! nstation is the number of spatial
! grid points.
INTEGER, INTENT(IN) :: nestedPer
! nestedPer is the so-called 'nested
! period', representing an inherent
! or hypothesized periodicity of the
! data, say 12 months for the annual
! cycle or 24 months for a biennial
! mode.
INTEGER, INTENT(IN) :: nmodes
! Number of cyclostationary EOFs
! (CSEOFs) to keep. Typically the
! first two or three CSEOFs are
! adequate. (In this routine nmodes
! must be less than or equal to
! ntime.)
REAL, DIMENSION(ntime,nstation), INTENT(IN) :: TSER
! Input data (time series) cast with
! ntime x nstation dimensionality.
! Time series may NOT have missing
! values. Thus, any time series with
! missing values, such as over land
! for an SST data set, should not
! be included.
REAL, INTENT(IN) :: msgOrig
! Missing value, such as _FillValue,
! of data BEFORE data is put into
! ntime x nstation form (see TSER).
! (If there is no missing value,
! supply a value outside the range
! of values represented by TSER.)
REAL, DIMENSION(nmodes,nestedPer,nstation), INTENT(OUT) :: CSEOF
! The CSEOFs. Suppose nestedPer = 24;
! then for each mode (CSEOF), say
! the first, there are 24 spatial
! patterns representing the cyclic
! nature of the given mode. Similarly
! for the second mode, and so on.
REAL, INTENT(IN) :: cseofScale
! Adjustable CSEOF scaling factor;
! (nominally set this to 1.0).
REAL, DIMENSION(nmodes,ntime), INTENT(OUT) :: PCTS
! Principal component time series,
! one for each mode. Each time series
! is normalized such that the
! deviation from the mean is per
! unit standard deviation; that
! is, PCTS(i) = [PCTS(i) - mean(i)]/
! sd(i) where i is the ith mode.
REAL, DIMENSION(nmodes), INTENT(OUT) :: PCVAR
! Percent variance explained by each
! mode.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Miscellaneous local variables:
REAL, DIMENSION(:,:,:), ALLOCATABLE :: TCOF
REAL, DIMENSION(:,:), ALLOCATABLE :: CF
REAL, DIMENSION(:,:), ALLOCATABLE :: SF
REAL, DIMENSION(:,:), ALLOCATABLE :: COV
REAL, DIMENSION(:), ALLOCATABLE :: PC
REAL, DIMENSION(:), ALLOCATABLE :: EOF
REAL, DIMENSION(:), ALLOCATABLE :: WGTS
INTEGER :: i
INTEGER :: im
INTEGER :: j
INTEGER :: k
INTEGER :: kk
INTEGER :: ktime
INTEGER :: npts
INTEGER :: nsub
INTEGER :: nwgt
INTEGER :: sta
REAL :: angle
REAL :: fact
REAL :: freq
REAL :: mn
REAL :: sd
REAL :: pi
REAL :: psum
REAL :: std
REAL :: twopi
REAL :: tvar
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Arguments for SSYEV routine (eigenvectors and eigenvalues of real symmetric matrix) from
! LAPACK library.
CHARACTER(1) :: arg_JOBZ
! Flag for computing only eigenvalues,
! or eigenvalues and eigenvectors.
CHARACTER(1) :: arg_UPLO
! Flag indicating storage of upper or
! lower triangle of arg_A.
INTEGER :: arg_N
! Number of rows (or columns) of arg_A.
REAL, DIMENSION(:,:), ALLOCATABLE :: arg_A
! On entry, real symmetric matrix
! arg_A. On exit, and when computing
! both eigenvalues and eigenvectors,
! contains the orthonormal
! eigenvectors of arg_A with the ith
! column of arg_A holding the
! eigenvector associated with
! arg_W(i).
INTEGER :: arg_LDA
! A parameter (typically N).
REAL, DIMENSION(:), ALLOCATABLE :: arg_W
! Contains the eigenvalues of arg_A in
! ascending order
REAL, DIMENSION(:), ALLOCATABLE :: arg_WORK
! One-dimensional work array.
INTEGER :: arg_LWORK
! A parameter (typically 3*N-1).
INTEGER :: arg_INFO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Diagnose basic restrictions on subprogram arguments:
! Check to make sure that
! (number of stations) * (nested period) >= (number of temporal samples):
IF ( .NOT.( (nstation*nestedPer) >= ntime ) ) THEN
WRITE(*,'(A56)') "(nstation*nestedPer) >= ntime restriction not satisfied"
WRITE(*,'(A33)') "in subroutine CYCLOSTATIONARY_WD."
WRITE(*,'(A61)') "Execution stopped in subroutine CYCLOSTATIONARY_WD: line 214"
STOP
END IF
! Check for missing data:
DO sta = 1, nstation
IF ( ANY( TSER(:,sta) == msgOrig ) ) THEN
WRITE(*,'(A60)') "Missing data found in TSER in subroutine CYCLOSTATIONARY_WD."
WRITE(*,'(A31)') "TSER may not have missing data."
WRITE(*,'(A26,I6)') "First detected at station ", sta
WRITE(*,'(A61)') "Execution stopped in subroutine CYCLOSTATIONARY_WD: line 228"
STOP
END IF
END DO
! Check to make sure that the number of modes to return is less than or equal to
! the number of temporal samples.
IF ( .NOT.( nmodes <= ntime ) ) THEN
WRITE(*,'(A42)') "nmodes <= ntime restriction not satisfied"
WRITE(*,'(A33)') "in subroutine CYCLOSTATIONARY_WD."
WRITE(*,'(A61)') "Execution stopped in subroutine CYCLOSTATIONARY_WD: line 242"
STOP
END IF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Initialize miscellaneous local variables:
npts = nestedPer/2
nsub = 1
nwgt = 2 * nestedPer
pi = 4.0*ATAN(1.0)
twopi = 2.0*pi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Allocate space for each array:
ALLOCATE( TCOF(1:ntime,1:nestedPer,1:nstation) )
ALLOCATE( WGTS(0:nsub*nwgt) )
ALLOCATE( CF(1:ntime,0:nestedPer) )
ALLOCATE( SF(1:ntime,0:nestedPer) )
ALLOCATE( COV(1:ntime,1:ntime) )
ALLOCATE( PC(1:nestedPer) )
ALLOCATE( EOF(1:nestedPer) )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Calculate the weights w(t) (i.e. WGTS), and coefficients a_k(r,t) (i.e. TCOF):
CALL WEIGHT( WGTS, nestedPer, npts, nwgt, nsub)
DO sta = 1, nstation
CALL INTEGR( TSER(:,sta), CF, SF, WGTS, ntime, nestedPer, ntime, npts, nwgt, nsub )
DO ktime = 1, ntime
TCOF(ktime,1,sta) = CF(ktime,0)
END DO
DO k = 2, nestedPer
kk = k/2
IF ( MOD(k,2) == 0 ) THEN
DO ktime = 1, ntime
TCOF(ktime,k,sta) = SQRT(2.0) * CF(ktime,kk)
END DO
ELSE
DO ktime = 1, ntime
TCOF(ktime,k,sta) = (-1.) * SQRT(2.0) * SF(ktime,kk)
END DO
END IF
END DO
IF ( nestedPer == nestedPer ) THEN
DO ktime = 1, ntime
TCOF(ktime,nestedPer,sta) = CF(ktime,npts)
END DO
END IF
END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute the real symmetric temporal covariance matrix:
DO i = 1, ntime
DO j = i, ntime
psum = 0.0
DO sta = 1, nstation
DO k = 1, nestedPer
psum = psum + TCOF(i,k,sta) * TCOF(j,k,sta)
END DO
END DO
COV(i,j) = psum / ( REAL(nstation) * REAL(nestedPer) )
COV(j,i) = COV(i,j)
END DO
END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Total variance (sum of diagonal elements of temporal covariance matrix):
tvar = 0.0
DO ktime = 1, ntime
tvar = tvar + COV(ktime,ktime)
END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Assign arguments for LAPACK SSYEV routine:
arg_JOBZ = "V" ! Compute both eigenvalues and
! eigenvectors.
arg_UPLO = "L" ! Store lower triangle of arg_A.
arg_N = ntime ! arg_N is simply ntime
ALLOCATE( arg_A( 1:ntime, 1:ntime ) )
arg_A = COV ! Assign arg_A as COV.
arg_LDA = ntime ! arg_LDA is simply ntime
ALLOCATE( arg_W( 1:ntime ))
arg_LWORK = 3*ntime - 1
! arg_LWORK is simply 3*ntime - 1.
ALLOCATE( arg_WORK( 1:arg_LWORK ) )
arg_INFO = 0 ! arg_info initialized as zero.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Call SSYEV:
CALL SSYEV( arg_JOBZ, &
arg_UPLO, &
arg_N, &
arg_A, &
arg_LDA, &
arg_W, &
arg_WORK, &
arg_LWORK, &
arg_INFO &
)
IF ( arg_INFO < 0 ) THEN
WRITE (*,'(A22)') "FATAL ERROR CONDITION:"
WRITE (*,'(I4,A20)') -arg_INFO,"th argument of SSYEV"
WRITE (*,'(A21)') "had an illegal value."
WRITE (*,'(A61)') "Execution stopped in subroutine CYCLOSTATIONARY_WD: line 372"
STOP
ELSE IF ( arg_INFO > 0 ) THEN
WRITE (*,'(A22)') "FATAL ERROR CONDITION:"
WRITE (*,'(A41)') "LAPACK SSYEV algorithm failed to converge"
WRITE (*,'(A19,I5,A13)') "after finding only ",arg_INFO-1," eigenvalues."
WRITE (*,'(A61)') "Execution stopped in subroutine CYCLOSTATIONARY_WD: line 378"
STOP
END IF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Sort eigenvectors and eigenvalues in reverse order:
arg_A(:,1:ntime) = arg_A(:,ntime:1:-1)
arg_W(:) = arg_W(ntime:1:-1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Get percent variance contributed by each mode from eigenvalue and total variance:
DO j = 1, nmodes
PCVAR(j) = arg_W(j) * 100.0 / tvar
END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute CSEOFs.
DO im = 1, nmodes
! ---------- PC patterns
DO sta = 1, nstation
DO k = 1, nestedPer
psum = 0.0
DO ktime = 1, ntime
psum = psum + arg_A(ktime,im) * TCOF(ktime,k,sta)
END DO
PC(k) = psum
END DO
! ---------- Bloch functions
DO i = 1, nestedPer
EOF(i) = PC(1)
END DO
DO k = 2, nestedPer
fact = SQRT(2.0)
IF (k == nestedPer) fact = 1.0
kk = k/2
freq = REAL(kk) * twopi / REAL(nestedPer)
IF ( MOD(k,2) == 0 ) THEN
DO i = 1, nestedPer
angle = freq * REAL(i-1)
EOF(i) = EOF(i) + fact * PC(K) * COS(angle)
END DO
ELSE
DO i = 1, nestedPer
angle = freq * REAL(i-1)
EOF(i) = EOF(i) + fact * PC(k) * SIN(angle)
END DO
END IF
END DO
DO i = 1, nestedPer
CSEOF(im,i,sta) = EOF(i)
END DO
END DO
! ---------- Normalization
psum = 0.0
DO sta = 1, nstation
DO i = 1, nestedPer
psum = psum + CSEOF(im,i,sta) * CSEOF(im,i,sta)
END DO
END DO
std = SQRT( psum / ( REAL(nstation) * REAL(nestedPer) ) )
DO i = 1, nestedPer
DO sta = 1, nstation
CSEOF(im,i,sta) = CSEOF(im,i,sta) * cseofScale / std
END DO
END DO
END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Normalize principal component time series.
DO j = 1, nmodes
arg_A(:,j) = arg_A(:,j)*SQRT(REAL(ntime))
mn = SUM( arg_A(:,j), DIM = 1 ) / REAL(ntime)
psum = 0.
DO ktime = 1, ntime
psum = psum + ( arg_A(ktime,j) - mn ) * ( arg_A(ktime,j) - mn )
END DO
sd = SQRT( psum / REAL(ntime-1) )
PCTS(j,:) = ( arg_A(:,j) - mn ) / sd
END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DEALLOCATE( TCOF )
DEALLOCATE( CF )
DEALLOCATE( SF )
DEALLOCATE( COV )
DEALLOCATE( PC )
DEALLOCATE( EOF )
DEALLOCATE( WGTS )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DEALLOCATE( arg_A )
DEALLOCATE( arg_W )
DEALLOCATE( arg_WORK )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Dependencies WEIGHT and INTEGR are contained in this subprogram:
CONTAINS
SUBROUTINE WEIGHT(WGTS, ICYC, NPTS, NWGT, NSUB)
IMPLICIT NONE
REAL, DIMENSION(0:NWGT*NSUB), INTENT(OUT) :: WGTS
INTEGER, INTENT(IN) :: ICYC
INTEGER, INTENT(IN) :: NPTS
INTEGER, INTENT(IN) :: NWGT
INTEGER, INTENT(IN) :: NSUB
REAL :: PI
INTEGER :: I
REAL :: T
REAL :: WGT
REAL :: SWGT
PI = 4.0*ATAN(1.0)
WGTS(0) = 1./REAL(ICYC)
DO I=1, NWGT*NSUB
T = REAL(I)/REAL(NSUB)
WGT = SIN(PI*T/REAL(ICYC))/(PI*T)
WGTS(I) = WGT
END DO
! ------- Normalization
SWGT = REAL(NSUB)
DO I=0, NWGT*NSUB
WGTS(I) = WGTS(I)/SWGT
END DO
END SUBROUTINE WEIGHT
SUBROUTINE INTEGR(TSER, CF, SF, WGTS, NTMAX, ICYC, NTOT, NPTS, NWGT, NSUB)
IMPLICIT NONE
INTEGER, INTENT(IN) :: NTMAX
INTEGER, INTENT(IN) :: ICYC
INTEGER, INTENT(IN) :: NWGT
INTEGER, INTENT(IN) :: NSUB
INTEGER, INTENT(IN) :: NTOT
INTEGER, INTENT(IN) :: NPTS
REAL, DIMENSION(1:NTMAX), INTENT(IN) :: TSER
REAL, DIMENSION(1:NTMAX,0:ICYC), INTENT(OUT) :: CF
REAL, DIMENSION(1:NTMAX,0:ICYC), INTENT(OUT) :: SF
REAL, DIMENSION(0:NWGT*NSUB), INTENT(IN) :: WGTS
REAL :: FRQ
REAL :: T
REAL :: SUM1
REAL :: SUM2
REAL :: TP
REAL :: VAL
REAL :: W1
REAL :: W2
INTEGER :: I
INTEGER :: J
INTEGER :: K
INTEGER :: IJ
INTEGER :: J1
INTEGER :: J2
INTEGER :: JJ
REAL :: PI
REAL :: TPI
INTEGER :: MWGT
PI = 4.0*ATAN(1.0)
TPI = 2.0*PI
MWGT = NWGT*NSUB
! ------- Modification to include seasonal cycle (K.-Y. Kim -- May 16, 2000)
DO K = 0, NPTS
FRQ = TPI*REAL(K)/REAL(ICYC)
DO I = 1, NTOT
T = REAL(I-1)
SUM1 = 0.0
SUM2 = 0.0
DO J = -MWGT, MWGT
TP = T + REAL(J)/REAL(NSUB)
IJ = IABS(J)
J1 = I + J/NSUB
IF (J1 < 1) THEN
JJ = MOD(J1+2*ICYC-1,ICYC)+1
VAL = TSER(JJ)
ELSE IF (J1 > NTOT) THEN
JJ = NTOT-MOD(NTOT-J1+2*ICYC,ICYC)
VAL = TSER(JJ)
ELSE
J2 = J1 + 1
W2 = MOD(J+MWGT,NSUB)/REAL(NSUB)
W1 = 1. - W2
VAL = TSER(J1)*W1 + TSER(J2)*W2
END IF
SUM1 = SUM1 + WGTS(IJ)*VAL*COS(FRQ*TP)
SUM2 = SUM2 - WGTS(IJ)*VAL*SIN(FRQ*TP)
END DO
CF(I,K) = SUM1
SF(I,K) = SUM2
END DO
END DO
END SUBROUTINE INTEGR
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END SUBROUTINE CYCLOSTATIONARY_WD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Back to David Stepaniak's Home Page