#include <params.h>
subroutine nmmatrix(zb ,zcr1 ,bm1 ,bmi ) 1,10
!-----------------------------------------------------------------------
!
! Purpose:
! Determine eigenvalues and left/right eigenvectors of the hydrostatic
! matrix. Results used to solve for two-time-level SLD, semi-implicit
! Coriolis term
!
! Author: J. Olson
!
!-----------------------------------------------------------------------
!
! $Id: nmmatrix.F90,v 1.2.42.2 2004/01/27 17:10:51 rosinski Exp $
! $Author: rosinski $
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid
use pspect
use abortutils, only: endrun
implicit none
!------------------------------Parameters-------------------------------
!
integer, parameter :: lwork = 100*plev ! work array dimension
!
!------------------------------Arguments--------------------------------
!
real(r8), intent(in) :: zb (plev,plev) ! semi-implicit matrix in d equation
real(r8), intent(out) :: zcr1(plev) ! eigenvalues (real)
real(r8), intent(out) :: bm1 (plev,plev) ! transpose of right eigenvector (normalized)
real(r8), intent(out) :: bmi (plev,plev) ! transpose of inverse of bm1
! ! ( = left eigenvectors -- normalized)
!
!---------------------------Local variables-----------------------------
!
integer k ! index
integer kk ! index
integer kkk ! index
integer istat ! status of SGEEV call
real(r8) eps ! used for epsilon tests
real(r8) sum ! tmp
real(r8) fmax ! tmp
real(r8) x (plev,plev) ! local copy of "zb"
real(r8) bmlft (plev,plev) ! left eigenvector
real(r8) bmrgt (plev,plev) ! right eigenvector
real(r8) zci (plev) ! eigenvalues (complex)
real(r8) work (lwork) ! work array
real(r8) diag ! tmp
logical lcomp ! error flag
logical lneg ! error flag
logical lsame ! error flag
logical lnonzero ! error flag
logical flag(plev) ! flag used in sorting eigenvalues/vectors
!
!-----------------------------------------------------------------------
!
! Matrix "zb" is stored 1st index = column; 2nd index = row of matrix.
! Transpose "zb" and store in work array before calling library routine
!
do k = 1,plev
do kk = 1,plev
x(k,kk) = zb(kk,k)
end do
end do
!
! Determine eigenvalues and left/right eigenvectors of "zb"
!
call sgeev('v' ,'v' ,plev ,x ,plev , &
zcr1 ,zci ,bmlft ,plev ,bmrgt , &
plev ,work ,lwork ,istat )
!
if (istat .ne. 0) then
write(6,*) 'ERROR in NMMATRIX: SGEEV returned "istat" = ',istat
call endrun ()
endif
if (int(work(1)) .gt. lwork) then
write(6,*) 'ERROR in NMMATRIX: dimension of work array must be at least', int(work(1))
call endrun ()
endif
!
! Find "eps" based upon the real eigenvalue with the largest magnitude
!
fmax = -1.e36
do kk = 1,plev
if(zcr1(kk) .gt. fmax) fmax = zcr1(kk)
end do
eps = fmax*epsilon(fmax)*10.
!
! Test that all eigenvalues are positive, real, and distinct
!
lcomp = .false.
lneg = .false.
lsame = .false.
do k = 1,plev
if( zcr1(k) .lt. 0.) lneg = .true.
if(abs(zci(k )) .gt. eps) lcomp = .true.
do kk = k+1,plev
if(abs( (zcr1(k) - zcr1(kk)) ) .lt. eps) lsame = .true.
end do
end do
if (lneg) then
call endrun ('ERROR in NMMATRIX: negative eigenvalue(s)')
endif
if (lcomp) then
call endrun ('ERROR in NMMATRIX: complex eigenvalue(s)')
endif
if (lsame) then
call endrun ('ERROR in NMMATRIX: non-distinct eigenvalue(s)')
endif
!
! Re-order the eigenvalues (and associated eigenvectors) in descending
! order. Use "zci", "bmi", and "bm1" as temporary arrays.
!
do k = 1,plev
flag(k) = .true.
end do
do k = 1,plev
fmax = -1.e35
do kk = 1,plev
if(flag(kk) .and. zcr1(kk) .gt. fmax) then
kkk = kk
fmax = zcr1(kk)
endif
end do
flag(kkk) = .false.
zci(k) = zcr1(kkk)
do kk = 1,plev
bmi(kk,k) = bmlft(kk,kkk)
bm1(kk,k) = bmrgt(kk,kkk)
end do
end do
!
! Copy arrays from temporaries back to original arrays
!
do k = 1,plev
zcr1(k) = zci(k)
do kk = 1,plev
bmlft(kk,k) = bmi(kk,k)
bmrgt(kk,k) = bm1(kk,k)
end do
end do
!
! Make sure the element of largest magnitude within each eigenvector
! is positive (in this case, we will choose the LEFT eigenvector).
! Adjust the signs of all other elements of the left AND associated
! RIGHT eigenvector, accordingly.
!
! NOTE: This is done only to ensure that the signs of all eigenvectors
! will be consistent across all computing platforms. This step is
! not necessary except that it makes debugging the "vertical normal
! mode" code MUCH easier.
!
do k = 1,plev
fmax = -1.e35
do kk = 1,plev
if(abs(bmlft(kk,k)) .gt. fmax) then
kkk = kk
fmax = abs(bmlft(kk,k))
endif
end do
if(bmlft(kkk,k) .lt. 0.) then
do kk = 1,plev
bmlft(kk,k) = -bmlft(kk,k)
bmrgt(kk,k) = -bmrgt(kk,k)
end do
endif
end do
!
! Normalize "bm1" and determine its inverse (and confirm).
!
! NOTE: both "bm1" and "bmi" are stored in transpose form (1st index =
! column; 2nd index = row of matrix) for later use in coding vector
! inner products using these matrices.
!
lnonzero = .false.
do k = 1,plev
sum = 0.
do kk = 1,plev
sum = sum + bmlft(kk,k)*bmrgt(kk,k)
end do
!
if(sum .lt. 0.) then
sum = sqrt(-sum)
do kk = 1,plev
bmi(kk,k) = -bmlft(kk,k)/sum
bm1(k,kk) = bmrgt(kk,k)/sum
end do
else
sum = sqrt(sum)
do kk = 1,plev
bmi(kk,k) = bmlft(kk,k)/sum
bm1(k,kk) = bmrgt(kk,k)/sum
end do
endif
end do
do k = 1,plev
do kk = 1,plev
diag = 0.
do kkk = 1,plev
diag = diag + bmi(kkk,kk)*bm1(k,kkk)
end do
if(kk .ne. k) then
if(diag .gt. eps) lnonzero = .true.
endif
end do
end do
if (lnonzero) then
call endrun ('NMMATRIX: Y(T)X has non-zero off-diagonal elements')
endif
!
return
end subroutine nmmatrix