        SUBROUTINE jpDRIV(IFUNC,RBUF,NBUF,CHR,LCHR)
        integer IFUNC, NBUF,LCHR
        REAL      RBUF(*)
        CHARACTER*50 CHR
        integer ixdim,iydim,lenbuf,ICOL
	include "driv.inc"
        data ixdim,iydim,lenbuf,lun,ICOL/292, 3200, 934400,99,0/
c       data buf/934400*0/
	if(ifunc.eq.28)goto 180	
        GOTO( 10, 20, 30, 40, 50, 60, 70,999, 90,100,
     &       110,120,130,140,150,999,999,180) IFUNC
        GOTO 999
C--- IFUNC= 1, Return device name.
10      CHR='JPS      (jiang, GRAY_pic to ps_file for print)'
	LCHR=LEN(CHR)
        RETURN
C--- IFUNC= 2, Return Physical min and max for plot device.
20      RBUF(1)=0
        RBUF(2)=2303                    !
        RBUF(3)=0
        RBUF(4)=3179                    !
        RBUF(5)=0
        RBUF(6)=1
        NBUF=6
        RETURN
C--- IFUNC= 3, Return device (X and Y) resolution in pixels per inch as
30      RBUF(1)=300.0
        RBUF(2)=300.0
        RBUF(3)=1
        NBUF=3
        LCHR=LEN(CHR)
        RETURN
C--- IFUNC= 4, Return misc device info.
40      CHR='HNNNNNNNNN'
        LCHR=10
        RETURN
C--- IFUNC= 5, Return default file name.
50      CHR='pgplot.ps'
        LCHR=LEN(CHR)
        RETURN
C--- IFUNC= 6, Return default physical size of plot.
60      RBUF(1)=0
        RBUF(2)=2303                    !
        RBUF(3)=0
        RBUF(4)=3179                    !
        RETURN
C--- IFUNC= 7, Return misc defaults.
70      RBUF(1)=1
        NBUF=1
        RETURN
C--- IFUNC= 9, Open workstation.
90      CALL grglun(lun)
	call unlink(CHR(:LCHR))
        OPEN(UNIT=lun,FILE=CHR(:LCHR),STATUS='NEW',IOSTAT=IER)
        IF(IER.NE.0) THEN
            CALL GRWARN('Cannot open graphics device '//CHR(:LCHR))
            RBUF(1)=lun
            RBUF(2)=0
	    close (lun)
		goto 999
        ELSE
            RBUF(1)=0
            RBUF(2)=1
        END IF
        call jpbuf0(934400,PIXMAP)
        RETURN
C--- IFUNC=10, Close workstation.
100	close (lun)
        RETURN
C--- IFUNC=11, Begin Picture.  Note, IXMIN=0 and IXMAX=RBUF(1) so,
C--- IXDIM=IXMAX-IXMIN+1=RBUF(1)+1
110     IxDIM=INT(RBUF(1)/8.)+1  
        IyDIM=RBUF(2)+1           
        LENBUF=IxDIM*IyDIM
        call grjp02(PIXMAP,lenbuf)
        RETURN
C--- IFUNC=12, Draw line.
120     continue
	CALL GRjp01(1,RBUF,ICOL,IxDIM,PIXMAP)
        RETURN
C--- IFUNC=13, Draw dot.
130     CALL GRjp01(0,RBUF,ICOL,IxDIM,PIXMAP)
        RETURN
C--- IFUNC=14, End picture.
140     CALL GRjp04(LUN,PIXMAP,IxDIM,IyDIM)
        RETURN
C--- IFUNC=15, Select pen.  Save pen number (up to 11) for
C--- possible use in pattern interior.
150     ICOL=MAX(0,MIN(NINT(RBUF(1)),11))
        RBUF(1)=MAX(0,MIN(ICOL,1))
        RETURN
c      for gray scale raster
180	call grjp03(rbuf,ixdim,PIXMAP)
	return
C--- Flag function not implemented.
999     NBUF=-1
        RETURN
        END

        subroutine jpbuf0(k,buf)
        byte buf(1)
        do i=1,k
        buf(i)=0
        enddo
        end

        SUBROUTINE GRjp01(LINE,RBUFI,ICOL,IXDIM,QXYBUF)         !
	integer*2 qxybuf(0:*),qmask(0:15),i2
        REAL    RBUFI(4),RBUF(4)  
        DATA   QMASK /32768,16384,8192,4096,
     *              2048,1024,512,256,128,64,32,16,8,4,2,1/

        RBUF(1)=rbufi(1)   
        RBUF(2)=3179-rbufi(2)
        RBUF(3)=rbufi(3)      
        RBUF(4)=3179-rbufi(4)   
	ixd2=ixdim/2
        IF(LINE.GT.0) THEN
          LENGTH=NINT(MAX(ABS(RBUF(3)-RBUF(1)),
     &                    ABS(RBUF(4)-RBUF(2))))
          D = MAX(1,LENGTH)
          XINC = (RBUF(3)-RBUF(1))/D
          YINC = (RBUF(4)-RBUF(2))/D
        ELSE
            LENGTH=0
            XINC=0.
            YINC=0.
        END IF
        XP = RBUF(1)
        YP = RBUF(2)
        IF (ICOL.GT.0) THEN
            DO K=0,LENGTH
                KX = NINT(XP)
                KY = NINT(YP)
		kk=ky*ixd2+kx/16
                QXYBUF(kk) = IOR(QXYBUF(kk),QMASK(MOD(Kx,16)))
                XP = XP + XINC
                YP = YP + YINC
            END DO
        ELSE
            DO K=0,LENGTH
                KX = NINT(XP)
                KY = NINT(YP)
		kk=ky*ixd2+kx/16
                i2=INOT(QMASK(MOD(Kx,16)))
                QXYBUF(kk) = IAND(QXYBUF(kk),i2)
                XP = XP + XINC
                YP = YP + YINC
            END DO
        END IF
        RETURN
        END

	subroutine grjp02(buf,n)
	integer*2 buf(*)
	do 10 i=1,n/2
10	buf(i)=0
	end

	subroutine grjp03(rbuf,ixdim,qxybuf)
	integer*2 qxybuf(0:*),qmask(0:15)
	integer*2 ii(18),jj(18),ix,iy
        REAL    RBUF(*)  
        DATA   QMASK /32768,16384,8192,4096,
     *              2048,1024,512,256,128,64,32,16,8,4,2,1/
	data ii/0,0,1,0,-1,1,1,2,-1,-1,-2,0,0,2,-2,1,-1,0/
	data jj/0,1,0,-1,0,1,-1,0,1,-1,0,2,-2,1,1,2,2,3/
	ixd2=ixdim/2
        ix= rbuf(1)
        iy= 3179-rbuf(2)
	kz= nint(rbuf(3))
	if(kz.lt.1)return
	is=nint(rbuf(4))
	if(kz.gt.is)kz=is
	do 10 i=1,kz
	kx=ix+ii(i)
	ky=iy+jj(i)
	kk=ky*ixd2+kx/16
10      QXYBUF(kk) = IOR(QXYBUF(kk),QMASK(MOD(Kx,16)))
        END

        SUBROUTINE GRjp04(LUN,map,IXDIM,IYDIM)
c map(72,3180), iydim=3180
        integer map(ixdim/4,iydim),c(4)
	character*594 line
	byte bl(594)
	equivalence (line,bl)
	write(lun,1)
1	format('%! Jiang_ps map_dump'/'erasepage initgraphics'/
     *  '[-1 0 0 -1 4700 6500] setmatrix'/'0 0 translate'/
     *  '2 2 scale'/
     *  '/r_data 16 string def'/'/r_zero 16 string def'/
     *  '/r_full 16 string def'/'/r_bkg1 16 string def'/
     *  'currentfile r_full readhexstring ',32('F')/
     *  'currentfile r_bkg1 readhexstring ',32('2')/
     *  'clear'/'/dline {0 exch translate 2304 1 true matrix'/
     *  '{currentfile read pop dup 97 eq'/
     *  '{pop currentfile r_data readhexstring pop}'/
     *  '{dup 98 eq'/'{pop r_zero}'/'{dup 99 eq'/'{pop r_full}'/
     *  '{100 eq'/'{r_bkg1}'/'ifelse}'/'ifelse}'/'ifelse}'/
     *  'ifelse'/'}imagemask} def')
	call decswap(map,ixdim*iydim/4)      ! for DECstation
	ll=0
	do 10 l=1,iydim
	k=0
	m=1
	do 20 j=1,18
	ip=0
	do 30 i=1,4
	k=k+1
	c(i)=map(k,l)
30	if(c(i).ne.0.and.c(i).ne.-1)ip=1
	if(ip.ne.0)then
	write(line(m:),"('a',4z8.8)")c
	do 40 mm=m+1,m+32
40	if(bl(mm).gt.80)bl(mm)=bl(mm)-32
	m=m+33
	else
	if(c(1).eq.0)line(m:m)='b'
	if(c(1).eq.-1)line(m:m)='c'
	m=m+1
	endif
20	continue
	if(line(1:18).eq.'bbbbbbbbbbbbbbbbbb')then
	ll=ll+1
	else
	k=1
	if(ll.gt.9)k=2
	if(ll.gt.99)k=3
	if(ll.gt.999)k=4
c	write(lun,'(i<k>," dline"/a)')ll,line(1:m-1)
	if(k.eq.1)write(lun,'(i1," dline"/a)')ll,line(1:m-1)
	if(k.eq.2)write(lun,'(i2," dline"/a)')ll,line(1:m-1)
	if(k.eq.3)write(lun,'(i3," dline"/a)')ll,line(1:m-1)
	if(k.eq.4)write(lun,'(i4," dline"/a)')ll,line(1:m-1)
	ll=1
	endif
10	continue
	write(lun,'("showpage")')
	end

	subroutine decswap(map,n)
	integer map(n),b
	integer*2 a(2),c
	equivalence (a,b)
	do 10 i=1,n
	b=map(i)
	c=a(1)
	a(1)=a(2)
	a(2)=c
10	map(i)=b
	end
		
