!>============================================================================< ! FILE: /fs/cgd/home0/davestep/f90_routines/SUB_eof_scalar.f90 ! ! AUTHOR: David Stepaniak, NCAR/CGD/CAS ! DATE INITIATED: Mon Jan 26 13:15:21 MST 1998 ! LAST MODIFIED: Fri Mar 13 08:31:36 MST 1998 ! ! SYNOPSIS: Computes EOFs, eigenvalues, and coefficients of EOFs for ! a scalar data field. ! ! DESCRIPTION: To fix ideas, suppose that one has constructed a data ! matrix consisting of M rows, each row corresponding to ! a spatial grid point, and N columns, each column ! corresponding to a temporal observation. Each row is ! simply a time series of observations for a particular ! grid point, and each column is simply a set of grid ! point observations at a particular time. We may call ! this data matrix F. ! ! The time series of F may consist of raw data, or anomalies ! computed using a climatology, or some other user-defined ! form. ! ! In what follows it will be assumed that, as a final step, ! the overall mean of each time series of F has been subtracted ! from its corresponding time series so as to form deviations ! from the mean. (This allows efficient matrix based ! calculation of the variance-covariance matrix and its ! relatives.) The matrix of deviations will be denoted by Fd. ! ! The M Empirical Orthogonal Functions (EOFs) of this data ! Fd may be obtained by computing the eigenvectors of the MxM ! variance-covariance matix R = Fd * Fd_T / N, where Fd_T is ! the transpose of Fd. Since it is most often the case that ! M >> N with M numbering in the thousands, the task of ! computing the eigenvectors (and eigenvalues) of R is ! computationally prohibitive. ! ! An alternative and efficient method put forth by von ! Storch and Hannoschock (1984) computes the N eigenvalues ! and eigenvectors of the NxN matrix formed by Fd_T * Fd / N, ! a much less computationally intensive exercise. From this ! NxN eigensystem it is straight forward to obtain the first ! N of M eigenvalues and eigenvectors of the original MxM ! system where the corresponding pairs of eigenvalues and ! eigenvectors are arranged in descending order according to ! eigenvalue. This subroutine implements the von Storch and ! Hannoschock algorithm. ! ! INPUT: ucM. Number of spatial grid points. ! ucN. Number of temporal observations. ! ucFd. Data matrix in which rows are time series of ! deviations from their corresponding overall time ! mean. ! ! OUTPUT: ucE. EOFs arranged column by column in order of their ! corresponding set of eigenvalues. ! ucL_N1. Eigenvalues arranged in descending order. ! ucC. Coefficent matrix of EOFs. ! ! NOTES: Input and output arguments are more fully described in the ! argument list below. ! ! Time series containing missing values are not allowed. It ! is left to the user to make provisions for time series ! containing missing values. This usually involves employing ! a grid point mask to assist in the packing and bookkeeping ! of the input data array and the unpacking and bookkeeping ! of the output EOF array. ! ! It is assumed that M > N. ! ! Utilizes SSYEV routine from LAPACK to compute eigenvalues ! and eigenvectors of a real symmetric matrix. On winterpark ! compile with -r4 and -lcomplib.sgimath flags where r4 ! indicates 32 bit (4 byte) storage of real variables and ! and complib.sgimath is the Silicon Graphics Scientific ! Mathematical Library which contains LAPACK routines among ! others. ! ! REFERENCE: von Storch, H., and G. Hannoschock (1984) Comments on ! "Emipirical Orthogonal Function Analysis of Wind Vectors ! Over the Tropical Pacific Ocean". Bulletin of the ! Meteorological Society of America, 65, 162. ! ! [Appeared as a letter to the editor concerning: Legler, ! D. M. (1983) Empirical Orthogonal Function Analysis of Wind ! Vectors over the Tropical Pacific Region. Bulletin of the ! Meteorological Society of America, 64, 234-241.] ! !>============================================================================< SUBROUTINE eof_scalar( ucM, ucN, ucFd, ucE, ucL_N1, ucC ) IMPLICIT NONE !>----------------------------------------------------------------------------< ! Argument list in order of appearance: INTEGER, INTENT(IN) :: ucM ! Number of spatial grid points to be ! represented in MxN data matrix. (ucM ! refers to upper case M.) INTEGER, INTENT(IN) :: ucN ! Number of temporal observations of a ! time series. (ucN refers to uppercase ! N.) It is assumed that M > N. REAL, DIMENSION( 1:ucM, 1:ucN ), INTENT(IN) :: ucFd ! MxN data array in which there are M ! rows (M grid points) each of which ! contains N columns (N temporal ! observations). Thus each row is ! simply a time series of observations ! at a particular grid point. !---------------- NOTE! ---------------+ ! Irregardless of how the data (in | ! particular the time series) are | ! prepared, AS A FINAL STEP the overall| ! mean of each time series should be | ! subtracted from its corresponding | ! time series so as to form deviations | ! from the mean. This allows for | ! efficient matrix-based calculation of| ! the variance-covariance matrix and | ! its relatives. | !--------------------------------------+ REAL, DIMENSION( 1:ucM, 1:ucN ), INTENT(OUT) :: ucE ! MxN array containing N EOFs ! (eigenvectors.) ! Each EOF has M components which ! correspond to the M spatial grid ! points. The EOFs are stored in ! columns in the same order as their ! corresponding eigenvalues arranged ! from greatest to least. REAL, DIMENSION( 1:ucN ), INTENT(OUT) :: ucL_N1 ! Nx1 one-dimensional array containing ! eigenvalues arranged from greatest ! to least. REAL, DIMENSION( 1:ucN, 1:ucN ), INTENT(OUT) :: ucC ! NxN array containing raw coefficients ! of EOFs. The raw coefficients ! form time series of the temporal ! behavoir of each EOF, and are stored ! row by row with the first row ! corresponding to the first EOF and so ! on. !>----------------------------------------------------------------------------< ! Declare useful transposes: REAL, DIMENSION(:,:), ALLOCATABLE :: ucFd_T ! NxM transpose array of ucFd. REAL, DIMENSION(:,:), ALLOCATABLE :: ucE_T ! NxM transpose array of ucE. !>----------------------------------------------------------------------------< ! Declare arrays and scalars to implement von Storch and Hannoschock (1984) ! algorithm: REAL, DIMENSION(:,:), ALLOCATABLE :: ucQ ! NxN array formed by the multiplication ! ucFd_T * ucFd / ucN => ucQ. Thus ! [NxM] * [MxN] => [NxN]. ! Note that ucQ is a real symmetric ! matrix. REAL, DIMENSION(:,:), ALLOCATABLE :: ucZ ! NxN array containing N eigenvectors ! of ucQ. The eigenvectors each have N ! components, and are stored in ! columns in the same order as their ! corresponding eigenvalues arranged ! from greatest to least. INTEGER :: col ! Simply a DO loop control variable for ! incrementing over columns. REAL, DIMENSION(:), ALLOCATABLE :: ucFdz ! Mx1 column vector (array) formed by the ! multiplication of ucFd and a single ! column (eigenvector) of ucZ. Thus ! ucFd * ucZ(:,col) => ucFdz, ! [MxN]* [Nx1] => [Mx1] REAL :: euc_norm_ucFdz ! Scalar. Euclidean norm of ucFdz. That ! is, the square root of the sum of the ! element by element multiplication of ! ucFdz and ucFdz. !>----------------------------------------------------------------------------< ! Arguments for LAPACK SSYEV routine. SSYEV computes the eigenvectors and ! eigenvalues of a real symmetric matrix. 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 ! Error trapping parameter diagnosed in ! code as below. ! +------------------------+ !>------------------------| BEGIN EXECUTABLE CODE: |--------------------------< ! +------------------------+ ! ! Form the NxN array Q = ucFd_T * ucFd / ucN: ALLOCATE( ucFd_T( 1:ucN, 1:ucM ) ) ucFd_T = TRANSPOSE( ucFd ) ALLOCATE( ucQ( 1:ucN, 1:ucN ) ) ucQ = MATMUL( ucFd_T, ucFd ) / REAL( ucN ) DEALLOCATE( ucFd_T ) !>----------------------------------------------------------------------------< ! 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 = ucN ! arg_N is simply ucN. ALLOCATE( arg_A( 1:ucN, 1:ucN ) ) arg_A = ucQ ! Assign arg_A as ucQ. DEALLOCATE( ucQ ) arg_LDA = ucN ! arg_LDA is simply ucN. ALLOCATE( arg_W( 1:ucN )) arg_LWORK = 3*ucN - 1 ! arg_LWORK is simply 3*ucN - 1. ALLOCATE( arg_WORK( 1:arg_LWORK ) ) arg_INFO = 0 ! arg_info initialized as zero. !>----------------------------------------------------------------------------< !> Compute eigenvalues and eigenvectors of ucQ: 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 (*,* ) "FATAL ERROR CONDITION: " WRITE (*,100) -arg_INFO,"th argument of SSYEV " WRITE (*,*) "had an illegal value. " WRITE (*,*) "Execution halted in subroutine eof_scalar." STOP ELSE IF ( arg_INFO > 0 ) THEN WRITE (*,*) "FATAL ERROR CONDITION: " WRITE (*,*) "LAPACK SSYEV algorithm failed to converge " WRITE (*,200) "after finding only ",arg_INFO-1, & " eigenvalues. " WRITE (*,*) "Execution halted in subroutine eof_scalar." STOP END IF DEALLOCATE( arg_WORK ) !>----------------------------------------------------------------------------< ! Extract ucE (EOFs) from ucZ: ALLOCATE( ucZ( 1:ucN, 1:ucN ) ) ucZ(:,1:ucN) = arg_A(:,ucN:1:-1) ! Arrange eigenvectors of ucZ so as to ! have same order as corresponding set ! of descending eigenvalues. DEALLOCATE( arg_A ) ALLOCATE( ucFdz( 1:ucM ) ) DO col = 1, ucN ! Loop over all eigenvectors of ucZ. ucFdz(:) = MATMUL( ucFd(:,:), ucZ(:,col) ) ! Compute ucFd * ucZ for a single ! eigenvector of ucZ. euc_norm_ucFdz = SQRT( & DOT_PRODUCT( ucFdz(:),ucFdz(:) ) & ) ! Compute Euclidean norm of ucFdz. ucE(:,col) = ucFdz(:) / euc_norm_ucFdz ! Normalize ucFdz by its Euclidean ! norm to get corresponding Mx1 ! eigenvector (EOF) of ucE. END DO DEALLOCATE( ucZ ) DEALLOCATE( ucFdz ) ALLOCATE( ucE_T( 1:ucN,1:ucM ) ) ucE_T = TRANSPOSE( ucE ) ucC = MATMUL( ucE_T, ucFd ) ! Compute coefficient matrix. DEALLOCATE( ucE_T ) ucL_N1(:) = arg_W(ucN:1:-1) ! Assign eigenvalues in descending order. DEALLOCATE( arg_W ) !>----------------------------------------------------------------------------< ! Format statements: 100 FORMAT( I1,A23 ) 200 FORMAT( A19,I4,A37 ) !>----------------------------------------------------------------------------< END SUBROUTINE eof_scalar !>----------------------------------------------------------------------------<