c***************************************************************
      subroutine initgraphics
c
      include 'plots.cmn'
c
      p_dographics=.true.
      p_dofile=.false.
      if (p_dographics) then
	call opngks
	call reset
      end if
      if (p_dofile) then
	open(unit=p_fileno,file='g.dat',status='NEW',form='FORMATTED')
	write(p_fileno,*) 'Graphics data output for p2d!'
	write(p_fileno,*)
	close(p_fileno)
      end if
c
      p_log_min=1.e-5
      p_graph_flags(0)=1
      p_graph_flags(1)=1
      p_graph_flags(2)=1
      p_graph_flags(3)=1
      p_graph_flags(4)=1
      p_surf_scale=.5
      p_surf_eye1=2.
      p_surf_eye2=-2.
      p_surf_eye=2.
      p_contour_flag=-15.
      p_title=' '
      p_xtitle=' '
      p_ytitle=' '
c
      return
      end
c***************************************************************
      subroutine endgraphics
      include 'plots.cmn'
c
      if (p_dographics) then
	call clsgks
      end if
cc      if (p_dofile) close(p_fileno)
c
      return
      end
c************************************************************
      subroutine graph(x,y,n)
      implicit none
c
      include 'plots.cmn'
c
      integer n,i
      real*4 x(n),y(n)
      integer*2 dshl
c
      if (p_dographics) then
	call anotat(p_xtitle,p_ytitle,0,0,-1,dshl)
	call ezxy (x,y,n,p_title)
      end if
      if (p_dofile) then
	open(unit=p_fileno,file='g.dat',status='OLD')
	write(p_fileno,*)'c==========================================='
	write(p_fileno,*) 'type: graph'
	write(p_fileno,*) p_title
	write(p_fileno,*) p_xtitle
	write(p_fileno,*) p_ytitle
	write(p_fileno,*) 'start_data'
	write(p_fileno,*) n
	do i=1,n
	  write(p_fileno,'(i2,1x,1pe9.2,1x,1pe9.2)') i,x(i),y(i)
	end do
	write(p_fileno,*) 'end_data'
	close(p_fileno)
      end if
      p_title=' '
      p_xtitle=' '
      p_ytitle=' '
c
      return
      end
c*************************************************************
      subroutine graphs(nvars,npts,di,x1,y1,y2,y3,y4)
      implicit none
c
c     nvars specify number of curves... 
c     di specify the slice from a larger matrix. 
c
      include 'plots.cmn'
      integer nvars,npts,di
      real*4 x1(npts),y1(npts*di),y2(npts*di),y3(npts*di),y4(npts*di)
      real*4 norm(p_max,4),mx(4)
      integer*2 dshl(10)
      integer i
c
      dshl(1)=2**4-1
      dshl(2)=2**6-1
      dshl(3)=2**8-1
      dshl(4)=2**10-1
c
      do i=1,4
	mx(i)=0.
      end do
      call normalize(y1,npts,di,p_graph_flags(1),norm(1,1),mx(1))
      if (nvars.ge.2) then
	if (p_graph_flags(2).eq.2) then
	  mx(2)=mx(1)
	  p_graph_flags(2)=p_graph_flags(1)
	end if
	call normalize(y2,npts,di,p_graph_flags(2),norm(1,2),mx(2))
      end if
      if (nvars.ge.3) then
	if (p_graph_flags(3).eq.2) then
	  mx(3)=mx(2)
	  p_graph_flags(3)=p_graph_flags(2)
	end if
	call normalize(y3,npts,di,p_graph_flags(3),norm(1,3),mx(3))
      end if
      if (nvars.ge.4) then
	if (p_graph_flags(4).eq.2) then
	  mx(4)=mx(3)
	  p_graph_flags(4)=p_graph_flags(3)
	end if
	call normalize(y4,npts,di,p_graph_flags(4),norm(1,4),mx(4))
      end if
c
      write(p_title,10) mx(1),mx(2),mx(3),mx(4)
   10 format(e9.3,' ',e9.3,' ',e9.3,' ',e9.3)
c
      if (p_dographics) then
	call anotat(p_xtitle,p_ytitle,0,0,-1,dshl)
	call ezmxy(x1,norm(1,1),p_max,nvars,npts,p_title)
      end if
c
      if (p_dofile) then
	open(unit=p_fileno,file='g.dat',status='OLD')
	write(p_fileno,*)'c==========================================='
	write(p_fileno,*) 'type: graphs'
	write(p_fileno,*) p_title
	write(p_fileno,*) p_xtitle
	write(p_fileno,*) p_ytitle
	write(p_fileno,*) 'start_data'
	write(p_fileno,*) nvars,npts
	do i=1,npts
	  write(p_fileno,'(i2,5(1x,1pe9.2))')
     1	   i,x1(i),y1(i),y2(i),y3(i),y4(i)
	end do
	write(p_fileno,*) 'end_data'
	close(p_fileno)
      end if
c
      p_title=' '
      p_xtitle=' '
      p_ytitle=' '
      return
      end
c*********************************************************
      subroutine surface(y,x1,n1,x2,n2,dim1,flag)
c
c     Same as countour for flag
c
      implicit none
      include 'plots.cmn'
      integer n1,n2,dim1,i,j,flag
      real*4 y(dim1,n1),x1(n1),x2(n2)
      real*4 y_tmp(p2_max,p2_max),yscale,s(p2_max,p2_max)
      real*4 scratch(2*p2_max*p2_max)
      real*4 stereo,sight(6),mxval,minval,avg
c
      do i=1,n1
	do j=1,n2
	  if (flag.eq.0) then
	    if (y(i,j).gt.p_log_min) then
	      y_tmp(i,j)=log10(y(i,j))
	    else
	      y_tmp(i,j)=log10(p_log_min)
	    end if
	  else
	    y_tmp(i,j)=y(i,j)
	  end if
	end do
      end do
      call info2d(y_tmp,n1,n2,p2_max,mxval,minval,avg)
      yscale=p_surf_scale*sqrt(1.*x1(n1)*x2(n2))
      do i=1,n1
	do j=1,n2
	  s(i,j)=yscale*(y_tmp(i,j)-minval)/(mxval-minval)
	end do
      end do
c eye position
      sight(1)=p_surf_eye1*x1(n1)
      sight(2)=p_surf_eye2*x2(n2)
      sight(3)=p_surf_eye*yscale
c where to look at
      sight(4)=.5*x1(n1)
      sight(5)=.5*x2(n2)
      sight(6)=yscale*(avg-minval)/(mxval-minval)
c
      write(p_xtitle,'(e9.3,2x,e9.3,2x,e9.3)') mxval,minval,avg
c
      stereo=0.0
c
      if (p_dographics) then
	call wtstr(.3955,.9765,p_title,2,0,0)
	call wtstr(.3955,.9300,p_xtitle,2,0,0)
	call srface(x1,x2,s,scratch,p2_max,n1,n2,sight,stereo)
      end if
      if (p_dofile) then
	open(unit=p_fileno,file='g.dat',status='OLD')
	write(p_fileno,*)'c==========================================='
	write(p_fileno,*) 'type: surface'
	write(p_fileno,*) p_title
	write(p_fileno,*) p_xtitle
	write(p_fileno,*) p_ytitle
	call blob_o_data(y,n1,n2,dim1)
	close(p_fileno)
      end if
      return
      end
c************************************************************
      subroutine blob_o_data(var,n1,n2,dim1)
      include 'plots.cmn'
c
      integer n1,n2,dim1,i,j
      real var(dim1,n2)
c
      write(p_fileno,*) 'start_data'
      write(p_fileno,*) n1,n2
      if (n2.eq.1) then
	do i=0,(n1-1)/5
	  write (p_fileno,'(5(1pe9.2,2x))') var(5*i+1,1),
     1	   var(5*i+2,1),var(5*i+3,1),var(5*i+4,1),var(5*i+5,1)
	end do
      else
	do j=1,n1
	  do i=0,(n2-1)/5
	    write(p_fileno,'(5(1pe9.2,2x))') var(j,5*i+1),
     1	     var(j,5*i+2),var(j,5*i+3),var(j,5*i+4),var(j,5*i+5)
	  end do
	end do
      end if
      write(p_fileno,*) 'end_data'
      return
      end
c************************************************************
      subroutine contour(f,n1,n2,dim1,flag,ipage)
c
c     flag=0 -> log plot       flag=1 -> linear
c
      implicit none
      include 'plots.cmn'
      integer n1,n2,dim1,ipage
      real*4 f(dim1,n2)
      real*4 flo,hi,finc,logf(p2_max,p2_max)
      integer nset,nhi,ndot,flag,i,j
c
      flo=0.
      hi=0.
      finc=-abs(p_contour_flag)
      nset=0
      nhi=0
      if ((p_contour_flag.lt.0.).and.(flag.ne.0)) then
	ndot=-2**7-2**6-2**5
      else
	ndot=0
      end if
c
      if (p_dographics) then
        call setzer
        call plchlq(.5,.95,p_title,20.,0.,0.)
c	call wtstr(.3955,.9765,p_title,2,0,0)
	if (flag.eq.0) then
	  do i=1,n1
	    do j=1,n2
	      if (f(i,j).gt.p_log_min) then
		logf(i,j)=log10(f(i,j))
	      else
		logf(i,j)=log10(p_log_min)
	      end if
	    end do
	  end do
	  call conrec(logf,p2_max,n1,n2,flo,hi,finc,nset,nhi,ndot)
	else
	  call conrec(f,dim1,n1,n2,flo,hi,finc,nset,nhi,ndot)
	end if
        call finish(ipage)
c	call frame
      end if
c
      if (p_dofile) then
	open(unit=p_fileno,file='g.dat',status='OLD')
	write(p_fileno,*)'c==========================================='
	write(p_fileno,*) 'type: contour'
	write(p_fileno,*) p_title
	write(p_fileno,*) p_xtitle
	write(p_fileno,*) p_ytitle
	call blob_o_data(f,n1,n2,dim1)
	close(p_fileno)
      end if
c
      p_title=' '
      p_xtitle=' '
      p_ytitle=' '
      return
      end
c************************************************************
      subroutine vector(v1,v2,n1,n2,dim1)
      implicit none
      include 'plots.cmn'
      integer n1,n2,dim1
      real*4 v1(dim1,n2),v2(dim1,n2)
      real*4 flo,hi,spv(2)
      integer nset,length,ispv
c
      flo=0.
      hi=0.
      nset=0
      length=0
      ispv=0
c
      if (p_dographics) then
	call wtstr(.3955,.9765,p_title,2,0,0)
	call velvct(v1,dim1,v2,dim1,n1,n2,flo,hi,nset,length,ispv,spv)
	call frame
      end if
c
      if (p_dofile) then
	open(unit=p_fileno,file='g.dat',status='OLD')
	write(p_fileno,*)'c==========================================='
	write(p_fileno,*) 'type: vector'
	write(p_fileno,*) p_title
	write(p_fileno,*) p_xtitle
	write(p_fileno,*) p_ytitle
	call blob_o_data(v1,n1,n2,dim1)
	call blob_o_data(v2,n1,n2,dim1)
	close(p_fileno)
      end if
      p_title=' '
      return
      end
c***********************************************************
      subroutine normalize(y,n,di,flag,norm,mx)
      implicit none
      include 'plots.cmn'
      integer n,di,j,k,flag
      real*4 y(n*di),norm(n),mx
      logical*1 set_max
c
      set_max=.true.
      if (mx.gt.0.) set_max=.false.
      do j=1,n
	k=1+(j-1)*di
	if (flag.eq.1) then
	  if (set_max) mx=max(mx,abs(y(k)))
	  norm(j)=y(k)
	else
	  if (y(k).gt.p_log_min) then
	    if (set_max) mx=max(mx,abs(log10(y(k))))
	    norm(j)=log10(y(k))
	  else
	    norm(j)=log10(p_log_min)
	    if (set_max) mx=abs(norm(j))
	  end if
	end if
      end do
      if (mx.eq.0.) mx=1.
      do j=1,n
	norm(j)=max(-1.,min(1.,norm(j)/mx))
      end do
      return
      end
c************************************************************
      subroutine info2d(y,n1,n2,d1,big,small,avg)
      integer n1,n2,d1,i,j
      real*4 y(d1,n2),big,small,avg
c
      big=y(1,1)
      small=y(1,1)
      avg=0.0
      do i=1,n1
	do j=1,n2
	  avg=avg+y(i,j)
	  big=max(big,y(i,j))
	  small=min(small,y(i,j))
	end do
      end do
      avg=avg/n1/n2
      if (big.eq.small) then
	if (big.eq.0) then
	  big=1.
	else
	  big=big+.1*abs(big)
	end if
      end if
      return
      end
c************************************************************
