program matinv_lu
!http://www.rcs.arch.t.u-tokyo.ac.jp/kusuhara/tips/linux/fortran.html
  implicit none
  integer, parameter :: N = 20, PI = acos(-1.0d0)
  integer :: i, j, ij, ji
  integer, parameter :: ijmax=120, jimax=120
  real(8) :: D(1:N,1:N), workD(1:N,1:N)
  real(8) :: B(1:N,1:N)
  real(8) :: C(1:N,1:N)  
  real(8) :: E(1:N,1:N)
  complex(8) :: Ei(1:N)
  real(8) :: U0, Ujh, Uj, Lamda, z, dz, f, p, l, k, M, beta, dl, dlamda, wimax, wimax2, LamdaMax, LMax
!----------------------------------------------------------------
  integer, parameter :: LDA = N, LDA3 = LDA*3, LWORK_for_eigen = N*4, LDVL=N, LDVR=N
      integer :: LWORK
      real(8) :: A(1:LDA,1:N), Wr(1:LDA), Wi(1:LDA), VL(1:LDA,1:N), VR(1:LDA,1:N), WORK_for_eigen(1:LDA3), Sigma(1:LDA), Sigma1, Sigma2!(1:LDA)!, C(1:LDA,1:N)
!----------------------------------------------------------------
  integer              :: IPIV(1:N), INFO
  real(8)              :: LWORK0
  real(8), allocatable :: WORK(:)

  z = 1.0d+4     ! height [m]
!  Lamda = 1.0d-3!0.     ! vertical shear
  U0 = 0. ! 
  dz = z/dble(N)          ! 
  f = 1.0d-4 ! Coriolis parameter
!  k = 1.0d-6  ! wave number [/m]
!  l = 2.0d0*pi/k    ! wave length [m]
  M = 1.0d-2  ! Brunt-Vaisala frequency [/s]
  beta = 0!1.6d-11 ! beta effect

!  ijmax = 2.0d2
!  dl = 1.25d5
  lmax=2.0d7

 ! jimax = 2.0d2
!  dlamda = 2.0d-5
  lamdamax=5.0d-3

  open(50,file='eigenvalue_lapack_for_unstable.out',status='unknown')
  
  do ji = 1,jimax
     Lamda = (LamdaMax/dble(jimax))*dble(ji)     ! vertical shear
!     Lamda = dlamda*dble(ji)
     
  do ij = 1,ijmax
     l = (LMax/dble(ijmax))*dble(ij)
!     l = dl*dble(ij)    ! wave length [m]
     k = 2.0d0*pi/l  ! wave number [/m]

  do i=1, N
     do j=1,N
        D(j,i)=dble(0)
     end do
  end do
  
  do i = 1,N-1
     D(i,i) = -dble(1) -(M*dz*k/f)**2              ! i = 1 : N
!     U = Lamda*dz*i
  end do
  do i = 1,N-1
     D(i,i+1) = dble(1)                        ! i = 1 : N-1
!     U = Lamda*dz*i
  end do
  do i = 2,N-1
     D(i,1:i-1) =  -(M*dz*k/f)**2        ! i = 2 : N
!     U = Lamda*dz*i
  end do
  do i = 1,N
     D(N,i) = -k**2
  end do


  workd=d

!  close(50)
  
! LAPACK$B$N%k!<%A%s$r8F$S=P$7$F9TNs$r(BLU$BJ,2r$7$^$9!#(B
! K $B$K$OJ,2r$7$?(B L $B$H(B U $B$,$"$o$;$F3JG<$5$l$F$$$^$9!#(B
! $B!J(BL$B$NBP3QMWAG!JA4$F(B1$B!K$O3JG<$5$l$^$;$s!#!K(B
! K $B$r=q$-49$($i$l$?$/$J$$>l9g$OCm0U$,I,MW$G$9!#(B
! IPIV$B$K$O(BLU$BJ,2r$r$KMQ$$$?%T%\%C%H$N>pJs$,3JG<$5$l$^$9!#(B
  call DGETRF( N, N, D, N, IPIV, INFO )
!  write(6,*) INFO
! LU$BJ,2r$5$l$?9TNs$N5U9TNs$r5a$a$k$N$K:GE,$J%o!<%/%9%Z!<%9$NBg$-$5$r5a$a$^$9!#(B
  call DGETRI(N, D, N, IPIV, LWORK0, -1, INFO)

! $B%o!<%/%9%Z!<%9MQ$NG[Ns$r3NJ]$7$^$9!#(B
  LWORK = INT(LWORK0)
  allocate( WORK(1:LWORK) )

! GETRF$B$G(BLU$BJ,2r$7$?7k2L$rMQ$$$F$b$H$N9TNs$N5U9TNs$r5a$a$^$9!#(B
! $B$b$H$N9TNs(B K $B$N5U9TNs$,(B K $B$K3JG<$5$l$^$9!#(B
  call DGETRI(N, D, N, IPIV, WORK, LWORK, INFO)

! $B%o!<%/%9%Z!<%9MQ$NG[Ns$r2rJ|$7$^$9!#(B
  deallocate( WORK )


  do i = 1,N
     do j = 1,N
        B(j,i) = dble(0)
     end do
  end do
  
  do i = 1,N-1
     Ujh = U0+Lamda*dz*dble(i)
     Uj = U0+Lamda*dz*(dble(i)-1./2.)
     B(i,i) = -Ujh -(1./2.)*Lamda*dz -(M*dz*k/f)**2*Uj +(M*dz/f)**2*beta      ! i = 1 : N
  end do
  do i = 1,N-1
     Ujh = U0+Lamda*dz*dble(i)
     B(i,i+1) = Ujh - (1./2.)*Lamda*dz                      ! i = 1 : N-1
  end do
  
  do j = 2,N-1 
   do i = 1,j-1
     Uj = U0+Lamda*dz*(dble(i)-1./2.)
     B(j,i) = -(M*dz*k/f)**2*Uj +(M*dz/f)**2*beta  ! i = 2 : N
   end do
  end do
  
  do i = 1,N
     Uj = U0+Lamda*dz*(dble(i)-1./2.)
     B(N,i) = beta - Uj*k**2
  end do
  
  
  C = matmul(D,B)
      

!c
      call DGEEV('N', 'N', N, C, LDA, Wr, Wi, VL, LDVL, VR, LDVR, WORK_for_eigen, LWORK_for_eigen, INFO)
!c

      wimax=maxval(wi)
      

!      if (wi(1) .eq. 0.) then
!         sigma1=999.
!      else if (wi(2) .eq. 0.) then
!         sigma2=999.
!      else
!         Sigma1=1./(8.64d4*k*Wi(1))
!         Sigma2=1./(8.64d4*k*Wi(2))
!      end if
      do i=1,N
        if (wimax .eq. 0.) then
           Sigma(i)=999.
        else
           Sigma(i)=1./(8.64d4*k*wimax)
        end if
      end do
      
   write(50,'(22E12.4)') l, lamda, (Sigma(i),i=1,N)!, Wi(i),i=1,N)!, Wr(i), i=1,N), k

! 600  format(i15, (4e15.8))

   end do
   end do
   close(50)
!gfortran eigenvalue_lapack_for_unstable_eady.f90 -llapack -lblas -fbounds-check
end program matinv_lu
