c     ************************************************************
c     *                                                          *
c     *     2D contour map (plotted using characters)            *
c     *                                                          *
c     ************************************************************
      implicit real*8 (a-h,o-z)
      parameter(nx=501,ny=501,ms=72)
      common /coor1/ imax,jmax
      common /coor2/ x(nx,ny),y(nx,ny)
      common /aray1/ a(ms,ms)
      common /aray2/ q(nx,ny,4),d(ms,ms,4)
      common /aray3/ m(ms,ms)

      dimension icom(20)
      dimension pcom(20)
      character*1 a
      character*50 file_name


c set arrays
      do j=1,ms
         do i=1,ms
            d(i,j,1) =0.0d0
            d(i,j,2) =0.0d0
            d(i,j,3) =0.0d0
            d(i,j,4) =0.0d0
            m(i,j)   =0
            a(i,j)   =' '
         end do
      end do

c computational mesh
      imax  =101         ! number of grid points in xi-direction
      jmax  =101         ! number of grid points in eta-direction
      xl    =  1.0d0     ! horizontal length of the computational domain
      yl    =  1.0d0     ! vertical length of the computational domain
      dx    =xl/dfloat(imax-1)
      dy    =yl/dfloat(jmax-1)
      do j=1,jmax
         do i=1,imax
            x(i,j) =dx*dfloat(i-1)
            y(i,j) =dy*dfloat(j-1)
         end do
      end do

c size of computational domain
      xmin = 1.0d30
      ymin = 1.0d30
      xmax =-1.0d30
      ymax =-1.0d30
      do j=1,jmax
         do i=1,imax
            xmin =dmin1(xmin,x(i,j))
            ymin =dmin1(ymin,y(i,j))
            xmax =dmax1(xmax,x(i,j))
            ymax =dmax1(ymax,y(i,j))
         end do
      end do


c read field data
      write(6,*) 'input: file name of field data'
      read(5,500) file_name
  500 format(a50)
      open(15,file=file_name,form='formatted')
      read(15,*) (icom(i),i=1,20)
      read(15,*) (pcom(i),i=1,20)
      read(15,*) ((q(i,j,1),q(i,j,2),q(i,j,3),q(i,j,4),i=1,imax),
     &             j=1,jmax)
      close(15)
      loop =icom(3)
      time =pcom(1)
      gamm =pcom(2)


c conversion to plotting area
      size =dmax1(xmax-xmin,ymax-ymin)
      ds   =size/dfloat(ms-1)
      do j=1,jmax-1
         do i=1,imax-1
            xg         =0.25d0*(x(i,j)+x(i+1,j)+x(i,j+1)+x(i+1,j+1))
            yg         =0.25d0*(y(i,j)+y(i+1,j)+y(i,j+1)+y(i+1,j+1))
            mi         =1+dint((xg-xmin)/ds)
            mj         =1+dint((yg-ymin)/ds)
            d(mi,mj,1) =d(mi,mj,1)+q(i,j,1)
            d(mi,mj,2) =d(mi,mj,2)+q(i,j,2)
            d(mi,mj,3) =d(mi,mj,3)+q(i,j,3)
            d(mi,mj,4) =d(mi,mj,4)+q(i,j,4)
            m(mi,mj)   =m(mi,mj)  +1
         end do
      end do
      do j=1,ms
         do i=1,ms
            if(m(i,j).gt.0) then
               d(i,j,1) =d(i,j,1)/dfloat(m(i,j))
               d(i,j,2) =d(i,j,2)/dfloat(m(i,j))
               d(i,j,3) =d(i,j,3)/dfloat(m(i,j))
               d(i,j,4) =d(i,j,4)/dfloat(m(i,j))
               uc       =d(i,j,2)/d(i,j,1)
               vc       =d(i,j,3)/d(i,j,1)
               pc       =(gamm-1.0d0)*(d(i,j,4)-0.5d0*(uc*d(i,j,2)
     &                                                +vc*d(i,j,3)))
               cc2      =gamm*pc/d(i,j,1)
               xm       =dsqrt((uc*uc+vc*vc)/cc2)
               sc       =pc/d(i,j,1)**gamm
               d(i,j,2) =pc
               d(i,j,3) =sc
               d(i,j,4) =xm
            endif
         end do
      end do


c set parameters
 1000 continue
      write(6,*) 'select: ro=1,p=2,s=3,xm=4'
      read(5,*) k
      if(k.le.0) k=1
      if(k.ge.5) k=4
      dmin = 1.0d30
      dmax =-1.0d30
      do j=1,ms
         do i=1,ms
            if(m(i,j).gt.0) then
               dmin =dmin1(dmin,d(i,j,k))
               dmax =dmax1(dmax,d(i,j,k))
            endif
         end do
      end do
      write(6,*) 'dmin=',dmin,' dmax=',dmax
      write(6,*) 'input: tmin and tmax'
      read(5,*) tmin,tmax
      tmin1 =dmin1(tmin,tmax)
      tmax1 =dmax1(tmin,tmax)
      tmin  =tmin1
      tmax  =tmax1


c set characters
      dt  =(tmax-tmin)/9.0d0
      tt0 =tmin
      tt1 =tt0+dt
      tt2 =tt1+dt
      tt3 =tt2+dt
      tt4 =tt3+dt
      tt5 =tt4+dt
      tt6 =tt5+dt
      tt7 =tt6+dt
      tt8 =tt7+dt
      tt9 =tt8+dt

      do j=1,ms
         do i=1,ms
            if(m(i,j).le.0) then
               a(i,j) =' '
            else
               if(d(i,j,k).lt.tt0) then
                  a(i,j) ='0'
               elseif(d(i,j,k).lt.tt1) then
                  a(i,j) ='1'
               elseif(d(i,j,k).lt.tt2) then
                  a(i,j) ='2'
               elseif(d(i,j,k).lt.tt3) then
                  a(i,j) ='3'
               elseif(d(i,j,k).lt.tt4) then
                  a(i,j) ='4'
               elseif(d(i,j,k).lt.tt5) then
                  a(i,j) ='5'
               elseif(d(i,j,k).lt.tt6) then
                  a(i,j) ='6'
               elseif(d(i,j,k).lt.tt7) then
                  a(i,j) ='7'
               elseif(d(i,j,k).lt.tt8) then
                  a(i,j) ='8'
               elseif(d(i,j,k).lt.tt9) then
                  a(i,j) ='9'
               else
                  a(i,j) ='*'
               endif
            endif
         end do
      end do


c display
      do j=ms,1,-1
         write(6,600) (a(i,j),i=1,ms)
  600    format(72a1)
      end do


c continue?
      write(6,*) ' '
      write(6,*) 'select : again(1)   end(-1)'
      read(5,*) ians
      if(ians.ge.0) then
         go to 1000
      endif


c termination
      stop
      end