*-----------------------------------------------------------------------
*     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    WMIN, WMAX, LMIN, LMAX
      PARAMETER( I_grd=400, J_grd=200 )
      PARAMETER( WMIN=0., WMAX=40.E+0, LMIN=0., LMAX=20.E+0 )
      character*4 data_number

      REAL    Psi_dcl(I_grd+1,J_grd+1),Omg_dcl(I_grd+1,J_grd+1)
      REAL    Dif_dcl(I_grd+1,J_grd+1)
      PARAMETER ( NBR3=15000*4)

      REAL IBR(NBR3)
      Parameter (I1=95, I2=105, J1=95, J2=105)      !平板の座標

      real   max, min, lelev

*--read data ----
c      do 100 k=0, 2000
            k=100
      call counter(k, data_number)
      open(7,file='omg_'//data_number//'.dat', 
     +          status='old')
            k=2000
      call counter(k, data_number)
      open(8,file='omg_'//data_number//'.dat', 
     +          status='old')

      do 10 j=1,J_grd+1
        do 10 i=1,I_grd+1
              read(7,*) Psi_dcl(i,j)
              read(8,*) Omg_dcl(i,j)
 10   continue
      close(unit=7)
      close(unit=8)

      do 21 j=1,J_grd+1
        do 21 i=1,I_grd+1
       omg_dcl(i,j)=omg_dcl(i,j)-psi_dcl(i,j) 
 21   continue

      max=1.D-10
      do 200 j=1,J_grd+1
       do 200 i=1,I_grd+1
          if (max .lt. abs(omg_dcl(i,j))) then
            max=abs(real(omg_dcl(i,j)))
          end if
 200   continue
   
*--出力ファイル名---
c      CALL SWCSET('FNAME','psi_'//data_number//'.eps)
      CALL SWCSET('FNAME', 'omg_dif_re_46_4020_'
     +//data_number//'')
*--コーナーマーク除去---
      CALL SGLSET('LCORNER', .FALSE.) 
*-- 欠損値処理 ----
      CALL GLLSET( 'LMISS', .TRUE. )
 
      do 20 i=I1, I2
        do 20 j=J1, J2
          psi_dcl(i,j)=999
          omg_dcl(i,j)=999
 20   continue
    
*-- グラフ描画 ----
c      WRITE(*,*) ' WORKSTATION ID (I)  ? ;'
      CALL SGPWSN
c      READ (*,*) IWS
      iws=2

      CALL GROPN( IWS )
c      call sldiv('y',1,2)
      CALL GRFRM

      CALL GRSWND( WMIN, WMAX, LMIN, LMAX )
      CALL GRSVPT(0.05, 0.95, 0.25, 0.7)
c      CALL GRSVPT(0.2, 0.8, 0.45, 0.55)
      CALL USPFIT
      CALL GRSTRF

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

c      CALL UDCNTR( Omg_dcl, I_grd+1, I_grd+1, J_grd+1 )
c      CALL UDCNTZ( Omg_dcl, I_grd+1, I_grd+1, J_grd+1, IBR, NBR3 )

c      CALL UDCNTZ( Psi_dcl, I_grd+1, I_grd+1, J_grd+1, IBR, NBR3 )

c for vorticity
c      write(6,*) max
c      min=real(-1.*max)
c      level=real(max/100.)
c      call udgcla(min/15.,max/15.,level)
c      call uddclv(0.)
      CALL UDCNTZ( omg_dcl, I_grd+1, I_grd+1, J_grd+1, IBR, NBR3 )
c
      CALL UXMTTL( 'T', 'Re=46                        ', 0.0 )
c
      CALL GRCLS

 100  continue

      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 
