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(nx,ny,4),e(ms,ms,4)
      common /aray3/ m(ms,ms)

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

c parameters
      imax  =301         ! number of grid points in xi-direction
      jmax  =201         ! number of grid points in eta-direction
      xl    =  1.5d0     ! horizontal length of the computational domain
      yl    =  1.0d0     ! vertical length of the computational domain

      jb =(1+jmax-jetwidth)/2
      je =(1+jmax+jetwidth)/2

      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
      xmin0 = 1.0d30
      ymin0 = 1.0d30
      xmax0 =-1.0d30
      ymax0 =-1.0d30
      do j=1,jmax
         do i=1,imax
            xmin0 =dmin1(xmin0,x(i,j))
            ymin0 =dmin1(ymin0,y(i,j))
            xmax0 =dmax1(xmax0,x(i,j))
            ymax0 =dmax1(ymax0,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 variables
      do j=1,jmax-1
         do i=1,imax-1
            rc       =q(i,j,1)
            uc       =q(i,j,2)/rc
            vc       =q(i,j,3)/rc
            pc       =(gamm-1.0d0)*(q(i,j,4)-0.5d0*(uc*q(i,j,2)
     &                                             +vc*q(i,j,3)))
            cc2      =gamm*pc/rc
            xm       =dsqrt((uc*uc+vc*vc)/cc2)
            sc       =pc/rc**gamm
            d(i,j,1) =rc
            d(i,j,2) =pc
            d(i,j,3) =sc
            d(i,j,4) =xm
         end do
      end do

c conversion to plotting area
      size0 =dmax1(xmax0-xmin0,ymax0-ymin0)

      do j=1,ms
         do i=1,ms
            e(i,j,1) =0.0d0
            e(i,j,2) =0.0d0
            e(i,j,3) =0.0d0
            e(i,j,4) =0.0d0
            m(i,j)   =0
            a(i,j)   =' '
         end do
      end do

      xor  =0.5d0*(xmax0+xmin0)
      yor  =0.5d0*(ymax0+ymin0)
      ds   =size0/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 =ms/2+dint((xg-xor)/ds)
            mj =ms/2+dint((yg-yor)/ds)
            if(mi.ge.1.and.mi.le.ms) then
               if(mj.ge.1.and.mj.le.ms) then
                  e(mi,mj,1) =e(mi,mj,1)+d(i,j,1)
                  e(mi,mj,2) =e(mi,mj,2)+d(i,j,2)
                  e(mi,mj,3) =e(mi,mj,3)+d(i,j,3)
                  e(mi,mj,4) =e(mi,mj,4)+d(i,j,4)
                  m(mi,mj)   =m(mi,mj)  +1
               endif
            endif
         end do
      end do
      do j=1,ms
         do i=1,ms
            if(m(i,j).gt.0) then
               e(i,j,1) =e(i,j,1)/dfloat(m(i,j))
               e(i,j,2) =e(i,j,2)/dfloat(m(i,j))
               e(i,j,3) =e(i,j,3)/dfloat(m(i,j))
               e(i,j,4) =e(i,j,4)/dfloat(m(i,j))
            endif
         end do
      end do

c set parameters
      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
      emin = 1.0d30
      emax =-1.0d30
      do j=1,ms
         do i=1,ms
            if(m(i,j).gt.0) then
               emin =dmin1(emin,e(i,j,k))
               emax =dmax1(emax,e(i,j,k))
            endif
         end do
      end do
      write(6,*) 'emin=',emin,' emax=',emax
      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(e(i,j,k).lt.tt0) then
                  a(i,j) ='0'
               elseif(e(i,j,k).lt.tt1) then
                  a(i,j) ='1'
               elseif(e(i,j,k).lt.tt2) then
                  a(i,j) ='2'
               elseif(e(i,j,k).lt.tt3) then
                  a(i,j) ='3'
               elseif(e(i,j,k).lt.tt4) then
                  a(i,j) ='4'
               elseif(e(i,j,k).lt.tt5) then
                  a(i,j) ='5'
               elseif(e(i,j,k).lt.tt6) then
                  a(i,j) ='6'
               elseif(e(i,j,k).lt.tt7) then
                  a(i,j) ='7'
               elseif(e(i,j,k).lt.tt8) then
                  a(i,j) ='8'
               elseif(e(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 termination
      stop
      end