da_eof_decomposition.inc

1 subroutine da_eof_decomposition( kz, bx, e, l )
2    
3    !---------------------------------------------------------------------------
4    ! Purpose: Compute eigenvectors E and eigenvalues L of vertical covariance 
5    !          matrix
6    !          B_{x} defined by equation:  E^{T} B_{x} E = L, given input kz x kz 
7    !          BE field.
8    !---------------------------------------------------------------------------
9    
10    implicit none
11 
12    integer, intent(in)      :: kz               ! Dimension of error matrix. 
13    real, intent(in)         :: bx(1:kz,1:kz)    ! Vert. background error.
14    real, intent(out)        :: e(1:kz,1:kz)     ! Eigenvectors of Bx.
15    real, intent(out)        :: l(1:kz)          ! Eigenvalues of Bx.
16 
17    integer                  :: work             ! Size of work array.
18    integer                  :: m                ! Loop counters
19    integer                  :: info             ! Info code.
20 
21    real                     :: work_array(1:3*kz-1)
22    real                     :: ecopy(1:kz,1:kz)
23    real                     :: lcopy(1:kz)
24 
25    !-------------------------------------------------------------------------
26    ! [5.0]: Perform global eigenvalue decomposition using LAPACK software:
27    !-------------------------------------------------------------------------
28    
29    work = 3 * kz - 1   
30    ecopy(1:kz,1:kz) = bx(1:kz,1:kz)
31    lcopy(1:kz) = 0.0
32 
33    call dsyev( 'V', 'U', kz, ecopy, kz, lcopy, work_array, work, info )
34    
35    if ( info /= 0 ) then
36       write(unit=message(1),fmt='(A,I4)') &
37          "Error in decomposition, info = ", info
38       call da_error(__FILE__,__LINE__,message(1:1))
39    end if
40    
41    ! Swap order of eigenvalues, vectors so 1st is one with most variance:
42    
43    do m = 1, kz
44       l(m) = lcopy(kz+1-m)
45       e(1:kz,m) = ecopy(1:kz,kz+1-m)
46    end do
47    
48 end subroutine da_eof_decomposition
49 
50