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