subroutine deofcovcor( data, nx, nt, icovcor, covariance, &
                               correlation, spacked, doswitched, iminnxnt )

!c	-------------------------------------------------
!c	Form the covariance or correlation matrix, put it 
!c	into 's'.  Note that 's' is always symmetric --
!c	the correlation between X and Y is the same as 
!c	between Y and X -- so use packed storage for 's'.
!c	The packed storage scheme we use is the same as
!c	LAPACK uses so we can pass 's' directly to the
!c	solver routine: the matrix is lower triangular,
!c	and s(i,j) = spacked(i+(j-1)*(1*n-j)/2).
!c
!c	Inputs:
!c		data(nx,nt): The basic data array.  THIS MUST
!c			BE ANOMALIES.
!c
!c		nx, nt: size of data array
!c			[INTEGER]
!c
!c		icovcor: if .eq. covariance, then calculate
!c			the covariance array;
!c			 if .eq. correlation, then calculate
!c			the correlation array.
!c			[INTEGER]
!c
!c		covariance, correlation: integer values to
!c			indicate each of these options.
!c			[INTEGER]
!c
!c		doswitched: if .TRUE., then calculate the
!c			'switched' array (which is of size
!c			(nt,nt)); if .FALSE., then calculate
!c			the normal array of size (nx,nx).
!c			[LOGICAL]
!c
!c		iminnxnt: min(nt,nx).  Used to dimension
!c			'spacked'.
!c
!c	Outputs:
!c
!c		spacked(iminnxnt*(iminnxnt+1)/2): the covariance 
!c			or correlation array.  This is in packed 
!c			form corresponding to LAPACK's lower
!c			triangular form.  
!c
!c	-------------------------------------------------

       implicit none

!c	-----------------
!c	Passed parameters
!c	-----------------
       integer :: nx, nt, icovcor, covariance, correlation, iminnxnt
       double precision :: data(nx,nt), spacked(iminnxnt*(iminnxnt+1)/2)
       logical :: doswitched

!c	---------------
!c	Local variables
!c	---------------
       integer :: i, j, k
       double precision :: sum, sum2, sum3, fact

       if( nx .le. 1 ) then
           write(0,*) 'deofcovcor.F: error: nx too small!! nx=', nx
           call exit(-1)
       endif
       if( doswitched ) then
         do j=1, nt
            do i=j, nt
               sum  = 0.0
               sum2 = 0.0
               sum3 = 0.0
               do k=1, nx
                  sum  = sum  + data(k,i)*data(k,j)
                  sum2 = sum2 + data(k,i)*data(k,i)
                  sum3 = sum3 + data(k,j)*data(k,j)
               enddo
               if( icovcor .eq. covariance ) then
                   fact = 1.0/float(nt-1)           ! Chuan theo nt (Wills)
          !T        fact = 1.0/float(nx-1)
               else
                  fact = 1.0/(sqrt(sum2)*sqrt(sum3))
               endif
               spacked(i+(j-1)*(2*nt-j)/2) = sum*fact
            enddo
         enddo
       else
         do j=1, nx
            do i=j, nx
               sum  = 0.0
               sum2 = 0.0
               sum3 = 0.0
               do k=1, nt
                  sum  = sum  + data(j,k)*data(i,k)
                  sum2 = sum2 + data(i,k)*data(i,k)
                  sum3 = sum3 + data(j,k)*data(j,k)
               enddo
               if( icovcor .eq. covariance ) then
                  fact = 1.0/float(nt-1)
               else
                  fact = 1.0/(sqrt(sum2)*sqrt(sum3))
               endif
               spacked(i+(j-1)*(2*nx-j)/2) = sum*fact
            enddo
         enddo
       endif

       return
       end