da_1d_eigendecomposition.inc

1 subroutine da_1d_eigendecomposition( bx, e, l )
2 
3    !------------------------------------------------------------------------------
4    !  Purpose: Compute eigenvectors E and eigenvalues L of vertical covariance matrix
5    !           B_{x} defined by equation:  E^{T} B_{x} E = L, given input 3D field of
6    !           errors (sum over all horizontal locations).
7    !------------------------------------------------------------------------------
8 
9    implicit none
10    
11    real, intent(in)         :: bx(:,:)          ! Global vert. background error.
12    
13    real, intent(out)        :: e(:,:)           ! Eigenvectors of Bx.
14    real, intent(out)        :: l(:)             ! Global eigenvalues of Bx.
15    
16    integer                  :: kz               ! Size of 3rd dimension.
17    integer                  :: m                ! Loop counters
18    integer                  :: work             ! Size of work array.
19    integer                  :: info             ! Info code.
20    
21    real, allocatable        :: ecopy(:,:)
22    real, allocatable        :: lcopy(:)
23    real, allocatable        :: work_array(:)
24 
25    !-------------------------------------------------------------------------
26    ! [1.0]: Initialise:
27    !-------------------------------------------------------------------------
28    
29    kz = size(bx, dim=1)
30    
31    !-------------------------------------------------------------------------
32    ! [5.0]: Perform global eigenvalue decomposition using LAPACK software:
33    !-------------------------------------------------------------------------
34    
35    work = 3 * kz - 1
36    allocate( work_array(1:work) )
37    
38    allocate( ecopy(1:kz,1:kz) )
39    allocate( lcopy(1:kz) )
40    
41    ecopy(:,:) = bx(:,:)
42 
43    lcopy = 0.0
44 
45    call dsyev( 'V', 'U', kz, ecopy, kz, lcopy, &
46               work_array, work, info )
47    
48    if ( info /= 0 ) then
49       write(unit=message(1),fmt='(A,I4,A)') &
50          ' da_1d_eigendecomposition: info = ', &
51          info,' - error in decomposition.'
52       call da_error(__FILE__,__LINE__,message(1:1))
53    end if
54    
55    !--Swap order of eigenvalues, vectors so 1st is one with most
56    !  variance, etc:
57    
58    do m=1,kz
59       l(m) = lcopy(kz+1-m)
60       e(1:kz,m) = ecopy(1:kz,kz+1-m)
61    end do
62    
63    deallocate( work_array )
64    deallocate( ecopy )
65    deallocate( lcopy )
66    
67 end subroutine da_1d_eigendecomposition
68 
69