      program matinv_lu
!http://www.rcs.arch.t.u-tokyo.ac.jp/kusuhara/tips/linux/fortran.html
      implicit none
      integer N, PI, nx, ny, IWS
      parameter (N = 20, PI = acos(-1.0d0), nx=120, ny=120)
      integer i, j, ij, ji, ni, nj, k
      real*8 Lamda, l
      real*8 imax, imin, jmax, jmin
!----------------------------------------------------------------
      real*8 Wr(1:N), Wi(1:N), sigma(1:N,1:Nx,1:Ny),
     &     sigma1(1:Nx,1:Ny), sigma2(1:Nx,1:Ny)
      real  work1(1:nx,1:ny), work2(1:nx,1:ny),work3(1:nx,1:ny),
     &     work5(nx,ny), work4(nx,ny), work6(1:nx,1:ny)
      real  xmax, xmin, ymax, ymin, vxmax, vxmin, vymax, vymin
      real*8 rmiss1, rmiss2, rmiss
!---missing value------------------------------------------------
      rmiss=999. 
!---set window---------------------------------------------------
      xmin=0.
      xmax=2.0d7!2.5d7
      ymin=0.
      ymax=5.0d-3!4.0d-3
!---set viewport-------------------------------------------------
      vxmax=0.8d0 
      vxmin=0.2d0  
      vymax=0.8d0  
      vymin=0.2d0  
!----------------------------------------------------------------      
      open(50, file='eigenvalue_lapack_for_unstable.out')
      do j=1,ny
         do i=1,nx
            read(50,100) l, lamda,(Sigma(k,i,j),k=1,N)!l, lamda, 
             work1(i,j)=sigma(1,i,j)

             if ( work1(i,j) .ge. 500. ) then
                work1(i,j) = RMISS             
!             else if ( work1(i+1,j+1) .eq. 500. ) then
              !  work1(i,j) = RMISS
             end if             

          end do
      end do

!      open(20, file='eigenvalue_lapack_for_unstable_sigma1.out')
!      do j=1, ny
!         do i=1,nx
      !   work(i,j)=(sigma(1,i,j),sigma(2,i,j),sigma(3,i,j),sigma(4,i,j),sigma(5,i,j))
!          write(20,100) sigma(1,i,j)
          !  work2(i,j)=sigma(2,i,j)
          !  work3(i,j)=sigma(3,i,j)
!             work(i,j)=(dble(i)+dble(j))/2.!10.*cos(2.*pi*real(i)/real(Nx))
!         end do
!      end do
!      close(20)
      close(50)

      
      CALL GLRGET( 'RMISS', RMISS )
      CALL GLLSET( 'LMISS', .TRUE. )

!      CALL GLRGET( 'RMISS', RMISS2 )
!      CALL GLLSET( 'LMISS', .TRUE. )      
      
      WRITE(*,*) ' WORKSTATION ID (I)  ? ;'
      CALL SGPWSN
      READ (*,*) IWS

      CALL GROPN( IWS )
      CALL GRFRM

      CALL GRSWND( XMIN, XMAX, YMIN, YMAX )
      CALL GRSVPT( 0.2, 0.8, 0.2, 0.8 )       
      CALL GRSTRN( 1 )
!      CALL USPFIT      
      CALL GRSTRF
      CALL USSTTL('Wave Length','m','Vertical Shear','/s')

      CALL USDAXS
     
      CALL UDSFMT( '(F4.1)' )
!      CALL UDGCLA( 1.5, 4.5, 1.5 )
!      CALL UDGCLA( 15, 45, 15 )
!--- clearly -----------------------------------------------------
!      CALL UDSCLV( 0., 1, 1, '0', 0.028 )
!      CALL UDSCLV( 1., 1, 1, '1', 0.028 )
      CALL UDSCLV( 1.5, 1, 1, '1.5', 0.028 )
      CALL UDSCLV( 2., 1, 1, '2', 0.028 )      
      CALL UDSCLV( 3., 1, 1, '3', 0.028 )
      CALL UDSCLV( 4., 1, 1, '4', 0.028 )
      CALL UDSCLV( 5., 1, 1, '5', 0.028 )
!--- dusty by conterlines ----------------------------------------
      CALL UDSCLV( 13., 1, 1, '13', 0.028 )
      CALL UDSCLV( 15., 1, 1, '15', 0.028 )
      CALL UDSCLV( 20., 1, 1, '20', 0.028 )
      CALL UDSCLV( 25., 1, 1, '25', 0.028 )
      CALL UDSCLV( 30., 1, 1, '30', 0.028 ) ! about 10-80 (~~~)
      CALL UDSCLV( 40., 1, 1, '40', 0.028 )
      CALL UDSCLV( 45., 1, 1, '45', 0.028 )
      CALL UDSCLV( 50., 1, 1, '50', 0.028 )
      CALL UDSCLV( 55., 1, 1, '55', 0.028 )
      CALL UDSCLV( 60., 1, 1, '60', 0.028 )
      CALL UDSCLV( 65., 1, 1, '65', 0.028 )
      CALL UDSCLV( 70., 1, 1, '70', 0.028 )
      CALL UDSCLV( 75., 1, 1, '75', 0.028 )
      CALL UDSCLV( 80., 1, 1, '80', 0.028 )
      CALL UDSCLV( 85., 1, 1, '85', 0.028 )
!--- dusty upper-right -------------------------------------------
      CALL UDSCLV( 150., 1, 1, '150', 0.028 )
      CALL UDSCLV( 120., 1, 1, '120', 0.028 )
      CALL UDSCLV( 100., 1, 1, '100', 0.028 )
      CALL UDSCLV( 180., 1, 1, '180', 0.028 )
      CALL UDSCLV( 200., 1, 1, '200', 0.028 )
      CALL UDSCLV( 300., 1, 1, '300', 0.028 )
      CALL UDSCLV( 400., 1, 1, '400', 0.028 )
      CALL UDSCLV( 500., 1, 1, '500', 0.028 )
      CALL UDSCLV( 600., 1, 1, '600', 0.028 )
      CALL UDSCLV( 700., 1, 1, '700', 0.028 )
!     CALL UDDCLV( 1.5 )
!      CALL UDDCLV( 0. )
!      CALL UDDCLV( 50. )
!      CALL UDDCLV( 4. )
!      CALL UDDCLV( 5. )
!      CALL UDDCLV( 999. )
      CALL UDCNTR( work1, NX, NX, NY )       
            
      CALL GRCLS
      
 100  format(22E12.4)
      end program matinv_lu



      
