*-----------------------------------------------------------------------
*     Copyright (C) 2000 GFD Dennou Club. All rights reserved.
*-----------------------------------------------------------------------
c contour plot for Nakamichi's program
c 2015.11.16 Takahiro IWAYAMA
c----------------------------------------------------------------------
      PROGRAM contour_dcl

      integer I_grd, J_grd
      real    x_MIN, x_MAX, y_MIN, y_MAX, pi
      PARAMETER( I_grd=64, J_grd=32 )
      character*4 data_number

      REAL    Psi_dcl(I_grd,0:J_grd+1),Omega_dcl(I_grd,0:J_grd+1)
      REAL    tmp_dcl(I_grd,0:J_grd+1)
      real    x, y
      PARAMETER ( NBR3=15000*4)

      REAL IBR(NBR3)
c      Parameter (I1=45, I2=55, J1=45, J2=55)      !平板の座標
     
      pi=acos(-1.)
      x_MIN=0.
c      x_MAX=dble(2)*pi
      x_MAX=dble(2)
      y_MIN=0.
c      y_MAX=dble(2)*pi
      y_MAX=dble(1)
*--read data ----
      do k=0,40
c      k=20
      call counter(k, data_number)

      open(7,file='psi_'//data_number//'.dat', 
     +          status='old', form='formatted')
      open(8,file='omg_'//data_number//'.dat', 
     +          status='old', form='formatted')
      open(9,file='tmp_'//data_number//'.dat', 
     +          status='old', form='formatted')

      do j=0,J_grd+1
        do i=1,I_grd
              read(7,*) x, y, Psi_dcl(i,j)
              read(8,*) x, y, Omega_dcl(i,j)
              read(9,*) x, y, tmp_dcl(i,j)
        end do
      end do
      close(unit=7)
      close(unit=8)

c      do j=1,J_grd
c        do i=1,I_grd
c              Psi_dcl(i,j)=dsin(dble(i)*dble(2)*3.14/dble(I_grd))
c              read(8,*) x, y, Omega_dcl(i,j)
c        end do
c      end do


*--出力ファイル名---
      CALL SWCSET('FNAME', 'Psi_1000')
*--コーナーマーク除去---
      CALL SGLSET('LCORNER', .FALSE.) 
*-- 欠損値処理 ----
c      CALL GLLSET( 'LMISS', .TRUE. )
 
c      do 20 i=I1, I2
c        do 20 j=J1, J2
c          psi_dcl(i,j)=999
c          omega_dcl(i,j)=999
c 20   continue
    
*-- グラフ描画 ----
c      WRITE(*,*) ' WORKSTATION ID (I)  ? ;'
      CALL SGPWSN
c      READ (*,*) IWS

c      CALL GROPN( IWS )
      CALL GROPN( 1 )
      CALL GRFRM

      CALL GRSWND( x_MIN, x_MAX, y_MIN, y_MAX )
      CALL USPFIT
      CALL GRSTRF

      CALL USSTTL( 'x', '', 'y', '' )
      CALL USDAXS


c      CALL UDCNTZ( Omega_dcl, I_grd, I_grd, J_grd, IBR, NBR3 )
c      call GRFRM
c      CALL UDCNTZ( Psi_dcl, I_grd, I_grd, J_grd+2, IBR, NBR3 )
      CALL UDCNTZ( tmp_dcl, I_grd, I_grd, J_grd+2, IBR, NBR3 )

      CALL GRCLS

      end do

      END
C----------------------------------------------------------------
      subroutine counter(i, data_number)
      integer   i, o_4, o_3, o_2, o_1
      character*10 i_data
      data      i_data/ '0123456789' /

      character*4 data_number

       o_4=i/1000
       o_3=(i-o_4*1000)/100
       o_2=(i-o_4*1000-o_3*100)/10
       o_1= i-o_4*1000-o_3*100-o_2*10
       data_number=i_data(o_4+1:o_4+1)//
     +                 i_data(o_3+1:o_3+1)//
     +                 i_data(o_2+1:o_2+1)//
     +                 i_data(o_1+1:o_1+1) 
      
      return
      end 
