var/da/da_vtox_transforms/da_check_eof_decomposition.inc

   1 subroutine da_check_eof_decomposition(be_eigenval, be_eigenvec, name)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Check eigenvectors E of vertical covariance matrix B_{x} which
   5    ! have been read in from NMC-statistics file.
   6    !
   7    ! Method:  E and L (eigenvalues) are computed using LAPACK/BLAS software in
   8    ! the NMC code using the definition E^{T} B_{x} E = L. This routine checks
   9    ! for eigenvector orthogonality and completeness as defined below.
  10    !---------------------------------------------------------------------------
  11
  12    implicit none
  13
  14    real*8, intent(in)            :: be_eigenval(:)   ! Back. error eigenvalues.
  15    real*8, intent(in)            :: be_eigenvec(:,:) ! Back. error eigenvector
  16    character(len=*), intent(in)  :: name             ! Variable name.
  17
  18    integer                       :: kz               ! Size of 3rd dimension.
  19    integer                       :: k, k1, k2, m     ! Loop counters
  20    real                          :: tot_variance     ! Total variance.
  21    real                          :: cum_variance     ! Cumulative variance.
  22    real                          :: max_off_diag     ! Maximum off-diagonal.
  23
  24    real, allocatable             :: matrix2(:,:)     ! 2D Work matrix.
  25    logical, allocatable          :: array_mask(:)    ! Array mask for MAXVAL.
  26
  27    if (trace_use) call da_trace_entry("da_check_eof_decomposition")
  28
  29    !----------------------------------------------------------------------
  30    ! [1.0]: Initialise:
  31    !----------------------------------------------------------------------
  32
  33    kz = size(be_eigenval)
  34
  35    !----------------------------------------------------------------------
  36    ! [2.0]: Print out global eigenvalues (used for truncation):
  37    !----------------------------------------------------------------------
  38
  39    cum_variance = 0.0
  40    tot_variance = sum(be_eigenval(1:kz))
  41
  42    write(unit=stdout,fmt='(A)') 'Name Mode  Eigenvalue Normalised Variance'
  43    do k = 1, kz
  44       cum_variance =  be_eigenval(k) + cum_variance
  45       write(unit=stdout,fmt='(A,I4,e14.4,f10.4)') &
  46          trim(name), k, be_eigenval(k), cum_variance / tot_variance
  47    end do
  48    write(unit=stdout,fmt=*) ' '
  49    write(unit=stdout,fmt='(A,e13.5)') 'Total variance = Tr(Bv) = ', tot_variance
  50    write(unit=stdout,fmt=*) ' '
  51
  52    !--------------------------------------------------------------------------
  53    ! [2.0]: Test global eigenvectors:
  54    !--------------------------------------------------------------------------
  55
  56    ! [2.1] Print out global eigenvectors:
  57
  58    write(unit=stdout,fmt='(2A)') 'Domain eigenvectors for ', trim(name)
  59
  60    write(unit=stdout,fmt='(50i13)')(m, m=1,kz)
  61    do k = 1, kz
  62       write(unit=stdout,fmt='(I3,50e13.5)')k, (be_eigenvec(k,m), m=1,kz)
  63    end do
  64    write(unit=stdout,fmt='(A)') " "
  65
  66    ! [2.2]: Test eigenvector orthogonality: sum_k (e_m(k) e_n(k)) = delta_mn:
  67
  68    allocate(array_mask(1:kz))
  69    allocate(matrix2(1:kz,1:kz))
  70
  71    write(unit=stdout,fmt='(2A)') &
  72       'Eigenvector orthogonality check for ', trim(name)
  73    write(unit=stdout,fmt='(A)')' Mode     Diagonal         Maximum off-diagonal'
  74    do k1 = 1, kz
  75       do k2 = 1, kz
  76          matrix2(k1,k2) = sum(be_eigenvec(:,k1) * be_eigenvec(:,k2))
  77       end do
  78
  79       array_mask(1:kz) =.true.
  80       array_mask(k1) = .false.
  81       max_off_diag = MAXVAL(ABS(matrix2(k1,:)),mask=array_mask(:))
  82       write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, matrix2(k1,k1), &
  83          max_off_diag
  84    end do
  85    write(unit=stdout,fmt=*) ' '
  86
  87    ! [2.3] Test eigenvectors completeness - sum_m (e_m(k1) e_m(k2)) = delta_k1k2
  88
  89    write(unit=stdout,fmt='(2A)') &
  90       'Eigenvector completeness check for ', trim(name)
  91    write(unit=stdout,fmt='(A)')' Level    Diagonal         Maximum off-diagonal'
  92    do k1 = 1, kz
  93       do k2 = 1, kz
  94          matrix2(k1,k2) = sum(be_eigenvec(k1,:) * be_eigenvec(k2,:))
  95       end do
  96
  97       array_mask(1:kz) =.true.
  98       array_mask(k1) = .false.
  99       max_off_diag = MAXVAL(ABS(matrix2(k1,:)),mask=array_mask(:))
 100       write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, matrix2(k1,k1), &
 101          max_off_diag
 102    end do
 103    write(unit=stdout,fmt=*) ' '
 104
 105    !-------------------------------------------------------------------------
 106    ! [3.0]: Tidy up:
 107    !-------------------------------------------------------------------------
 108
 109    deallocate(matrix2)
 110    deallocate(array_mask)
 111
 112    if (trace_use) call da_trace_exit("da_check_eof_decomposition")
 113
 114 end subroutine da_check_eof_decomposition
 115
 116