program matinv_lu
!http://www.rcs.arch.t.u-tokyo.ac.jp/kusuhara/tips/linux/fortran.html
  implicit none
  integer, parameter :: N = 5, PI = acos(-1.0d0)
  integer :: i, j, ij, ijmax, ji, jimax
  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, f0, y
!----------------------------------------------------------------
  integer, parameter :: LDA = N, LDA3 = LDA*4, 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)!, C(1:LDA,1:N)
!----------------------------------------------------------------
  integer              :: IPIV(1:N), INFO
  real(8)              :: LWORK0
  real(8), allocatable :: WORK(:)

  z = 1.0d+4     ! height [m]
  Lamda = 0.!4.0d-3     ! vertical shear
  U0 = 10. ! 
  dz = z/dble(N)          ! 
  f0 = 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 = 1.6d-11 ! beta effect

  f=f0+beta*dl*ijmax
  
  ijmax = 1.0d2
  dl = 2.5d5

  jimax = 1.0d2
  dlamda = 5.0d-2
  
!  do ji = 1,jimax
!     Lamda = dlamda*dble(ji)     ! vertical shear

     
  do ij = 1,ijmax
     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
!  do i = 1,N-2
!     D(i,i+2:N) = dble(0)                      ! i = 1 : N-2
!     U = Lamda*dz*i
!  end do


  
!!!  write(*,'(A)') 'D ='
!!!  write(*,'(5E12.4)') ( D(i,:), i=1,N )

  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 )

  
!!!  write(*,'(A)') 'D-1 ='
!!!  write(*,'(5E12.4)') ( D(i,:), i=1,N )

  
!!!  write(*,'(A)') 'D*D-1 ='
!!!  write(*,'(5E12.4)') matmul(D,workD)
  
!  real(8) :: B(1:N,1:N)

  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
!     write(6,*) U
!     stop
  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-1
!     U = Lamda*dz*(dble(i)+1./2.)
!     B(i,i+1:N) = dble(0)                                   ! i = 1 : N-2
!  end do
  
  do i = 1,N
     Uj = U0+Lamda*dz*(dble(i)-1./2.)
     B(N,i) = beta - Uj*k**2
  end do
  

  
!!!  write(*,'(A)') 'B ='
!!!  write(*,'(5E12.4)') ( B(i,:), i=1,N )
  
  C = matmul(D,B)

!!!  write(*,'(A)') 'D-1*B ='  
!!!  write(*,'(5E12.4)') ( C(i,:), i=1,N )
      

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


   open(50,file='eigenvalue_lapack_for_neutral.out',status='unknown')
    
   !      write(50,*) 'ij ='
!      wimax=maxval(Wi)   
      write(50,'(22E12.4)') l, lamda, (Wr(i), i=1,N) !, (Wr(i), i=1,N), (8.64d4/k/Wi(i), i=1,N)
      
!      write(50,*) 'i =, Wr =, Wi ='
!      do i = 1, N
!         write(50,600) i, Wr(i), Wi(i)
!      end do
!c
! 600  format(i15, (4e15.8))
!c
   end do
!   end do
   close(50)
!gfortran eigenvalue_lapack_for_neutral_10.f90 -llapack -lblas -fbounds-check
end program matinv_lu
