program eigenvalu

  use dcl

  integer :: i,j
  integer,parameter :: nx=300, ny=300
  real,parameter :: xmin=0., xmax=2.0d7, ymin=0., ymax=5.0d-4
  real,dimension(0:nx,0:ny) :: Sig
  real(8) :: L, Lam, k, eigen

open(17, file='green.dat', status='unknown')
  do j=1,ny
     do i=1,nx
        L       = ( xmin + (xmax-xmin)*i/dble(nx) )
!         Lam    = 1.0d-3
        Lam     = ( ymin + (ymax-ymin)*j/dble(ny) ) ! [m/s/mb]
        k       = dble(2)*3.14159/L
        
        if (eigen(L,Lam) .eq. 0) then
           Sig(i,j) = 999
           write (17,'(3E12.3)')  L, Lam, Sig(i,j)
        else
           Sig(i,j) = 1/(eigen(L,Lam)*k*9.d4)
           write (17,'(3E12.3)')  L, Lam, Sig(i,j)
        end if
!        write (17,'(3E12.3)')  L, Lam, Sig(i,j)
     end do
  end do
close(17)

  CALL GLRGET( 'RMISS', RMISS )
  CALL GLLSET( 'LMISS', .TRUE. )
  call swcset('FNAME','hoge')
  call DclOpenGraphics()
  call DclNewFrame

  call DclSetWindow( xmin, xmax, ymin, ymax )
  call DclSetViewPort( 0.2, 0.8, 0.2, 0.8 )
  call DclSetTransFunction


  call DclDrawScaledAxis



  call UDSCLV(  1.5, 2, 1,  '1.5', 0.001 )
  call UDSCLV(  2.0, 2, 1,  '2.0', 0.001 )
  call UDSCLV(  3.0, 2, 1,  '3.0', 0.001 )
  call UDSCLV(  4.0, 2, 1,  '4.0', 0.001 )
  call UDSCLV(  5.0, 2, 1,  '5.0', 0.001 )
  call UDSCLV( 13.0, 2, 1, '13.0', 0.001 )
  call UDSCLV( 15.0, 2, 1, '15.0', 0.001 )
  call UDSCLV( 20.0, 2, 1, '20.0', 0.001 )
  call UDSCLV( 25.0, 2, 1, '25.0', 0.001 )
  call UDSCLV( 200.0, 2, 1, '200.0', 0.001 )

  call DclDrawContour( Sig )

  call DclCloseGraphics

end program eigenvalu

function eigen(L,Lam)
  
    implicit none

    real(8) :: L
    real(8) :: Lam  
    real(8) :: eigen

    integer, parameter :: N = 20
    integer :: h,i,j
    real(8) :: D(1:N,1:N), B(1:N,1:N), E(1:N,1:N), DB(1:N,1:N)
    complex :: C(1:N)

    integer              :: IPIV(1:N), LWORK1, INFO
    real(8)              :: LWORK0
    real(8), allocatable :: WORK1(:)

    integer, parameter :: lwork=3*n
    real(8)            :: wi(n), wr(n), work(3*n)
    real(8)            :: vl(n, n), vr(n, n)

    real(8) :: k,S,dp,f,be,P0,a
    real(8) :: U(N)

       
    P0 =1.0d5              ! surface pressure [Pa]
!    L  =1.0d7              ! wave length [m]
!    Lam =3.0d-4            ! vertical shear -dU/dp [m/sec*Pa]
!    Lam =0.   
    k  =dble(2)*3.14159/L   ! wave number [/m]
    S  =2.0d-6             ! measure of the static stability [m^2*/sec^2*Pa^2]
    dp =1.0d5/dble(N)      ! pressure difference between the levels j and j+1 [Pa]
    f  =1.0d-4             ! the Coriolis parameter [/sec]
    be =1.6d-11            ! latitude variation the Coriolis parameter [m/sec]
!    be =0.
    
    a  =S*(dp/f)**2
    
    do i =1,N
       U(i) = Lam*(P0 - dble(i)*dp)
    end do

    D      = 0.
    D(1,1) = -a*k**2 - 1.0d0
    D(1,2) = 1.0d0
    do j=2,N-1
       D(j,1:j-1) = -a*k**2
       D(j,j)     = -a*k**2-1.0d0	
       D(j,j+1)   =  1.0d0
    end do
    D(N,1:N)=k**2
    
    B     =  0.
    B(1,1)= -U(2) + a*(-k**2*U(1) + be)
    B(1,2)=  U(1)
    do j=2,N-1
       do i=1,j-1
          B(j,i)= -a*(k**2*U(i)-be)
       end do
       B(j,j)     = -U(j+1) - a*(k**2*U(j)-be)	
       B(j,j+1)   =  U(j)     
    end do

    do j=1,N
       B(N,j)=k**2*U(j)-be
    end do

!    write(*,'(A)') 'D ='
!    write(*,'(5E11.3)') ( D(i,:), i=1,N )
!    write(*,'(A)') 'B ='
!    write(*,'(5E11.3)') ( B(i,:), i=1,N )

    call DGETRF( N, N, D, N, IPIV, INFO )

    LWORK1 = N
    allocate( WORK1(1:N) )

    call DGETRI(N, D, N, IPIV, WORK1, LWORK1, INFO)

!    write(*,'(A)') 'D^-1 ='
!    write(*,'(10E11.3)') ( D(i,:), i=1,N )

    deallocate( WORK1 )

    DB = matmul(D,B)
    call dgeev('N', 'N', n, DB, n, wr, wi, vl, n, vr, n, work, lwork, info)

    eigen = maxval(wi)
    
end function eigen
