@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
***********************************************************************
*
*  Menu driven program to operate volume rendering routine in file
*    'VRENDER FORTRAN'
*
*  Requires the following additional subroutines be available
*    vrinit - initializes data routines
*    vrfunc - computes value of function at given set of points
*    vrgrad - computes the gradient of function at a given point
*
*    startit- performs any necessary initial graphics commands and
*              returns the device number
*    vrndpl - plots the picture
*    getout - performs any necessary closing graphics commands
*
*    vrrays - volume renders isovalue surfaces of a 3-D function
*
*                      Written by Daniel Kartch
*
*      Version 6.4                               September 26, 1988
*
*            Modified for parallel execution by Hugh Caffey
*      Version ?.?                                     13 July 1989
***********************************************************************
      PROGRAM vrend
*
      INCLUDE (newvrcm2)
*
      REAL*8 btime(4), etime(4)
      INTEGER*4 devtyp,opt,ierr,i
      LOGICAL qdidit,qexist
      CHARACTER fname*8,ftype*8,fd*1,cmdstr*80
*
      dparm = 0.
      anno1 = ' '
      qdidit = .false.
      INCLUDE (timing) 1
      call dinit
*
*  main menu
   10 continue
      rewind(5)
      INCLUDE (cmsclr) 1
      print*,'  0> Exit'
      print*,'  1> Data/device initialization'
      if (qdidit) then
         if (filmon.eq.0) then
           print*,'  2> Plot picture'
         else
           print*,'  2> Film from ',filmfr,' to ',filmto
         endif
         print*,'  3> Viewing distance: ',viewd
         print*,'  4> Look at: ',(lookat(i),i=1,3)
         print*,'  5> View angle: ',vtheta,vphi
         print*,'  6> Up angle: ',upang
         print*,'  7> Screen height: ',scrht
         print*,'  8> Near and far clipping distances: ',near,far
         print*,'  9> Near and far fade distances: ',nrfade,frfade
         print*,' 10> Ray segment size: ', delta
         print*,' 11> Film setup'
         print*,' 12> Save settings'
         print*,' 13> Load settings'
         if (numshl .eq. 0) then
            print*,'*14> SET SHELL LEVELS AND COLORS!'
         else
            print*,' 14> Set shell levels and colors'
         endif
         print*,' 15> Ambient and diffuse reflection coeff:',
     &           (reflct(i),i=1,2)
         print*,' 16> Specular reflection coeff and exp:',
     &           (reflct(i),i=3,4)
         print*,' 17> Light angle:  ',ltheta,lphi
         print*,' 18> Annotation: ',anno1
      endif
*
      read(5,*,err=10,end=10)option
*
*  quit
      if (option.eq.0) then
        call getout
*
*  initialize function
      else if (option.eq.1) then
 100    continue
        rewind(5)
        INCLUDE (cmsclr) 1
        print*,' 0> Return to main menu'
        print*,' 1> Read data file'
        print*,' 2> Select output device'
        read(5,*,err=100,end=100)opt
        if (opt .eq. 0) then
           if (qdidit .and. dparm .ne. 0.) goto 10
           if (.not. qdidit) print*,'No data!'
           if (dparm .eq. 0.) print*,'No device!'
           print *,'Please select again'
        else if (opt .eq. 1) then
           call vrinit(datmin,datmax)
           qdidit = .true.
        else if (opt .eq. 2) then
           call startit
        else
           print*,'Invalid selection'
        endif
        goto 100
*
*  plot picture
      else if (option.eq.2) then
        if (filmon.eq.0) then
          if (dparm .ge. 2. .and. dparm .le. 4.) then
             ftype = 'PIXELS'
             ftype(7:7) = char(int(dparm)+ichar('0')-1)
200          print*,'File name (up to 8 characters)?'
             read(5,'(a8)')fname
             INCLUDE (cmsio1) 1
             INCLUDE (aixio1) 2
             INCLUDE (mvsio1) 3
          endif
*
          print*,'calculating picture'
          call vrrays
*
          print*,'plotting picture'
          call vrndpl(screen, npixel, asprat)
          if (dparm .ge. 2. .and. dparm .le. 4.) close(2)
        else
          do 210 i = filmfr, filmto
            INCLUDE (cmsclr) 1
            print*,'frame #',i
*
            call frame(i)
*
            print*,'calculating picture'
            call vrrays
*
            print*,'plotting picture'
            call vrndpl(screen, npixel, asprat)
  210     continue
        endif
        goto 10
*
*  set viewing distance
      else if (option.eq.3) then
300     print*,'enter new viewing distance:'
        read(5,*,err=300,end=10)viewd
        goto 10
*
*  set look at point
      else if (option.eq.4) then
400     print*,'enter new look at point (x,y,z):'
        read(5,*,err=400,end=10)(lookat(i),i=1,3)
        goto 10
*
*  set viewing angle
      else if (option.eq.5) then
500     print*,'enter new viewing angle (theta,phi):'
        read(5,*,err=500,end=10)vtheta,vphi
        goto 10
*
*  set up angle
      else if (option.eq.6) then
600     print*,'enter new up angle (radians from theta dir):'
        read(5,*,err=600,end=10)upang
        goto 10
*
*  set screen height
      else if (option.eq.7) then
700     print*,'enter new screen height:'
        read(5,*,err=700,end=10)scrht
        goto 10
*
*  set clipping distances
      else if (option.eq.8) then
800     print*,'enter new clipping distances (near,far):'
        read(5,*,err=800,end=10)near,far
        goto 10
*
*  set fade distances
      else if (option.eq.9) then
900     print*,'enter new fade distances (near,far):'
        read(5,*,err=900,end=10)nrfade,frfade
        goto 10
*
*  set ray segment size
      else if (option.eq.10) then
1000    print*,'enter new ray segment size:'
        read(5,*,err=1000,end=10)delta
        goto 10
*
*  film setup
      else if (option.eq.11) then
        call filmse
        goto 10
*
*  save settings
      else if (option.eq.12) then
        call varsve
        goto 10
*
*  save settings
      else if (option.eq.13) then
        call varld
        goto 10
*
*  set shells
      else if (option.eq.14) then
        call setshl
        goto 10
*
*  set ambient and diffuse reflection
      else if (option.eq.15) then
1500    print*,'enter new ambient reflection coefficient:'
        read(5,*,err=1500,end=10)reflct(1)
1510    print*,'enter new diffuse reflection coefficient:'
        read(5,*,err=1510,end=10)reflct(2)
        goto 10
*
*  set specular reflection
      else if (option.eq.16) then
1600    print*,'enter new specular reflection coefficient:'
        read(5,*,err=1600,end=10)reflct(3)
1610    print*,'enter new specular reflection exponent:'
        read(5,*,err=1610,end=10)reflct(4)
        goto 10
*
*  set light angle
      else if (option.eq.17) then
1700    print*,'enter new light angle (theta,phi):'
        read(5,*,err=1700,end=10)ltheta,lphi
        goto 10
*
*  set annotation text
      else if (option.eq.18) then
        print*,'enter annotation text:'
        read(5,'(a60)',end=10)anno1
        anno = anno1
        call cvtstr(anno,60)
        goto 10
*
*
      else
        goto 10
      endif
*
      INCLUDE (timing2) 1
*
      STOP
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
@PROCESS DC(datblk)
      SUBROUTINE dinit
*
      IMPLICIT NONE
      INCLUDE (newvrcm2)
      INCLUDE (datacom2)
      REAL*4 xtemp,ytemp,ztemp
*
      xtemp = float(xsize)
      ytemp = float(ysize)
      ztemp = float(zsize)
c
      lookat(1) = xtemp * 0.5
      lookat(2) = ytemp * 0.5
      lookat(3) = ztemp * 0.5
      viewd  = anint(2.*(xtemp*xtemp+ytemp*ytemp+ztemp*ztemp)**0.5)
      vtheta =   0.E0
      vphi   =   0.E0
      upang  =   3.14
      ltheta =   0.7854
      lphi   =   0.E0
      reflct(1) = 0.35
      reflct(2) = 0.9
      reflct(3) = 0.4
      reflct(4) = 50.E0
      scrht     = amax1(xtemp,ytemp,ztemp)
      near    = -viewd
      far     =  viewd
      delta   =   1.E0
      nrfade  = -viewd
      frfade  =  viewd
      numshl  =   0
      thresh  =   7.8E-3
      filmst  =   0
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
      SUBROUTINE filmse
* set up for filming
*
      IMPLICIT NONE
      INCLUDE (newvrcm2)
      INTEGER*4 option
*
 100  continue
      INCLUDE (cmsclr) 1
*
      print*,'0> Return to main menu'
      if (filmst.eq.0) then
        print*,'1> Create new film file'
        print*,'2> Load film file'
      else
        print*,'1> Close film file'
        if (filmon.eq.0) then
          print*,'2> Film mode: off'
        else
          print*,'2> Film mode: on'
        endif
        print*,'3> Film from: ',filmfr
        print*,'4> Film to:   ',filmto
      endif
*
 200  continue
      rewind(5)
      read(5,*,err=200,end=200) option
*
*  quit
      if (option.eq.0) then
        goto 2000
*
*  create file
      else if ((option.eq.1).and.(filmst.eq.0)) then
        call filmcr
        goto 100
*
*  open file
      else if ((option.eq.2).and.(filmst.eq.0)) then
        call filmld
        filmst = 1
        goto 100
*
*  close file
      else if ((option.eq.1).and.(filmst.eq.1)) then
        close(42)
        filmst = 0
        filmon = 0
        goto 100
*
*  toggle film on/off
      else if ((option.eq.2).and.(filmst.eq.1)) then
        filmon = 1 - filmon
        goto 100
*
*  set first frame
      else if ((option.eq.3).and.(filmst.eq.1)) then
 300    continue
        rewind(5)
        print*,'start at what frame?'
        read(5,*,err=300,end=300)filmfr
        if (filmfr.gt.filmnd) then
          print*,' Out of range.  Enter value from 1 to',filmnd
          goto 300
        endif
        goto 100
*
*  set last frame
      else if ((option.eq.4).and.(filmst.eq.1)) then
 400    continue
        rewind(5)
        print*,'end after what frame?'
        read(5,*,err=400,end=400)filmto
        if (filmto.gt.filmnd) then
          print*,' Out of range.  Enter value from 1 to',filmnd
          goto 400
        endif
        goto 100
*
      else
        print*,'There is no option #',option
        goto 200
      endif
*
 2000 continue
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
      SUBROUTINE filmcr
*  Create new film file
*
      IMPLICIT NONE
      INCLUDE (newvrcm2)
      INTEGER*4 i,j,len,ierr
      CHARACTER*8 fname, varnam(nfilmv)
      CHARACTER*1 fdisk
*
      DATA varnam(1)/'DATAPARM'/
      DATA (varnam(i),i=2,3)/'VTHETA  ','VPHI    '/
      DATA varnam(4)/'NUMSHL  '/
      DATA ((varnam(i+j),i=5,nshlv*5,5),j=0,4)/nshlv*'LEVEL   ',
     &     nshlv*'OPACITY ',nshlv*'HUE     ',nshlv*'LIGHT   ',
     &     nshlv*'SAT     '/
*
 10   FORMAT(I5)
 20   FORMAT(I5,2x,A8,2I10,2F15.5)
*
      print*,'enter length of new film (0 to cancel):'
      read(5,*)filmnd
      if (filmnd.eq.0) goto 2000
 30   print*,'enter name of new file:'
      read(5,'(A8)')fname
      INCLUDE (cmsio2) 1
      INCLUDE (aixio2) 2
      INCLUDE (mvsio2) 3
      if (ierr .ne. 0) then
        print*,'Error opening film file!'
        goto 30
      endif
*
*
      do 200 i = 1,2
        filmv(i,1) = 0.0
        filmv(i,2) = vtheta
        filmv(i,3) = vphi
        filmv(i,4) = numshl
        do 190 j = 1,nshlv
          filmv(i,5*j)   = shells(j)
          filmv(i,5*j+1) = opac(j)
          filmv(i,5*j+2) = hlscol(1,j)
          filmv(i,5*j+3) = hlscol(2,j)
          filmv(i,5*j+4) = hlscol(3,j)
  190   continue
  200 continue
*
  220 write(42,10)filmnd
      do 300 i = 1,nfilmv
        write(42,20)i, varnam(i), 1, filmnd, filmv(1,i), filmv(2,i)
  300 continue
      close(42)
*
 2000 continue
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
      SUBROUTINE filmld
*  load film file
*
      IMPLICIT NONE
      INCLUDE (newvrcm2)
      REAL*4      v1, v2
      INTEGER*4   vnum, f1, f2, j, i
      LOGICAL qexist
      CHARACTER*8 fname, vname
      CHARACTER*1 fdisk
*
 10   FORMAT(I5)
 20   FORMAT(I5,2x,A8,2I10,2F15.5)
*
 30   print*,'enter name of file to load:'
      read(5,'(A8)')fname
      INCLUDE (cmsio3) 1
      INCLUDE (aixio3) 2
      INCLUDE (mvsio3) 3
*
      do 100 j = 1, nfilmv
        filmrn(1,j) = 0
        filmrn(1,j) = 0
 100  continue
*
      read(42,10)filmnd
 200  read(42,20,end=500)vnum,vname,f1,f2,v1,v2
      if (f1.eq.1) then
        filmrn(1,vnum) = f1
        filmrn(2,vnum) = f2
        filmv(1,vnum) = v1
        filmv(2,vnum) = v2
      endif
      goto 200
*
 500  continue
*
      filmfr = 1
      filmto = filmnd
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
      SUBROUTINE frame(framno)
*  set variables for current frame
*
      IMPLICIT NONE
      INCLUDE (newvrcm2)
      REAL*4 rdummy,datprm
      INTEGER*4 framno, idummy, i
*
      datprm = 0.
      call setvar(framno, datprm, idummy, 1)
      call setvar(framno, vtheta, idummy, 2)
      call setvar(framno, vphi, idummy, 3)
      call setvar(framno, rdummy, numshl, 4)
      do 10 i = 1, nshlv
        call setvar(framno, shells(i), idummy, 5*i)
        call setvar(framno, opac(i), idummy, 5*i+1)
        call setvar(framno, rdummy, hlscol(1,i), 5*i+2)
        call setvar(framno, rdummy, hlscol(2,i), 5*i+3)
        call setvar(framno, rdummy, hlscol(3,i), 5*i+4)
        if(i.le.numshl)
     &    call hlsrgb(float(hlscol(1,i)), float(hlscol(2,i))/32767.0,
     &                  float(hlscol(3,i))/32767.0,
     &                  rgbcol(1,i),rgbcol(2,i),rgbcol(3,i))
   10 continue
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
      SUBROUTINE setvar(fr, var, ivar, vnum)
*  set specific variable
*
      IMPLICIT NONE
      INCLUDE (newvrcm2)
      REAL*4 var, dum5, dum6
      INTEGER*4 fr, vnum, dum1, dum3, dum4, ivar
      CHARACTER*8 dum2
*
 10   FORMAT(I5)
 20   FORMAT(I5,2x,A8,2I10,2F15.5)
*
      if (filmrn(1,vnum).ne.0) then
        if ((fr.lt.filmrn(1,vnum)).or.(fr.gt.filmrn(2,vnum))) then
          rewind(42)
          read(42,10)dum1
 100      read(42,20)dum1,dum2,dum3,dum4,dum5,dum6
          if ((dum1.eq.vnum).and.(dum3.le.fr).and.(dum4.ge.fr)) then
            filmrn(1,vnum) = dum3
            filmrn(2,vnum) = dum4
            filmv(1,vnum) = dum5
            filmv(2,vnum) = dum6
          else
            goto 100
          endif
        endif
*
        if (filmrn(1,vnum).ne.filmrn(2,vnum)) then
          var = filmv(1,vnum) + (filmv(2,vnum) - filmv(1,vnum))
     &        * float(fr - filmrn(1,vnum)) / float(filmrn(2,vnum)
     &        - filmrn(1,vnum))
        else
          var = filmv(1,vnum)
        endif
        ivar = nint(var)
      endif
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
      SUBROUTINE varsve
* save current settings
*
      IMPLICIT NONE
      INCLUDE (newvrcm2)
      INTEGER*4 i,j,ierr,len
      CHARACTER*8 vfile
      CHARACTER*1 vdisk
*
 30   print*,'enter name of new save file:'
      read(5,'(A8)')vfile
      INCLUDE (cmsio4) 1
      INCLUDE (aixio4) 2
      INCLUDE (mvsio4) 3
* error on open
      if (ierr .ne. 0) then
         print*,'Error opening save file!'
         goto 30
      endif
*
 1010 format(' Look at: ',3(2x,F9.5))
 1020 format(' Viewing angle: ',2(2x,F8.5))
 1030 format(' Viewing distance: ',2x,F10.5)
 1040 format(' Up angle: ',2x,F8.5)
 1050 format(' Screen height: ',2x,F10.5)
 1060 format(' Near and far clipping distances: ',2(5x,F10.5))
 1070 format(' Near and far fade distances:     ',2(5x,F10.5))
 1080 format(' Ray segment size: ',2x,F8.5)
 1100 format(' Number of shells: ',2x,I2)
 1110 format(' Shell #',I2)
 1120 format(' Level: ',2x,E10.4,15x,'Opacity: ',2x,F7.5)
 1130 format(' Hue: ',2x,I3,6x,'Lightness: ',2x,I5,6x,'Saturation: ',
     &       2x,I5)
 1140 format(' Light angle: ',2(2x,F8.5))
 1150 format(' Ambient & diffuse reflection coefficients: ',2(2x,F7.5))
 1160 format(' Specular reflection coefficient & exponent: ',2x,F7.5,
     &       2x,F9.5)
 1170 format(' Opacity cutoff: ',2x,F7.5)
 1180 format(' Annotation: ',2x,a60,a1)
*
      write(43,1010) (lookat(i),i=1,3)
      write(43,1020) vtheta,vphi
      write(43,1030) viewd
      write(43,1040) upang
      write(43,1050) scrht
      write(43,1060) near,far
      write(43,1070) nrfade,frfade
      write(43,1080) delta
      write(43,1100) numshl
*
      do 100 i = 1,numshl
        write(43,1110) i
        write(43,1120) shells(i),opac(i)
        write(43,1130) (hlscol(j,i),j=1,3)
 100  continue
*
      write(43,1140) ltheta,lphi
      write(43,1150) reflct(1),reflct(2)
      write(43,1160) reflct(3),reflct(4)
      write(43,1170) thresh
      write(43,1180) anno1,'.'
*
      close(unit=43)
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
      SUBROUTINE varld
* load current settings
*
      IMPLICIT NONE
      INCLUDE (newvrcm2)
      INTEGER*4 i,j,dummy
      LOGICAL qexist
      CHARACTER*8 vfile
      CHARACTER*1 vdisk
*
 30   print*,'enter name of file to load:'
      read(5,'(A8)')vfile
      INCLUDE (cmsio5) 1
      INCLUDE (aixio5) 2
      INCLUDE (mvsio5) 3
*
 1010 format(' Look at: ',3(2x,F9.5))
 1020 format(' Viewing angle: ',2(2x,F8.5))
 1030 format(' Viewing distance: ',2x,F10.5)
 1040 format(' Up angle: ',2x,F8.5)
 1050 format(' Screen height: ',2x,F10.5)
 1060 format(' Near and far clipping distances: ',2(5x,F10.5))
 1070 format(' Near and far fade distances:     ',2(5x,F10.5))
 1080 format(' Ray segment size: ',2x,F8.5)
 1100 format(' Number of shells: ',2x,I2)
 1110 format(' Shell #',I2)
 1120 format(' Level: ',2x,E10.4,15x,'Opacity: ',2x,F7.5)
 1130 format(' Hue: ',2x,I3,6x,'Lightness: ',2x,I5,6x,'Saturation: ',
     &       2x,I5)
 1140 format(' Light angle: ',2(2x,F8.5))
 1150 format(' Ambient & diffuse reflection coefficients: ',2(2x,F7.5))
 1160 format(' Specular reflection coefficient & exponent: ',2x,F7.5,
     &       2x,F9.5)
 1170 format(' Opacity cutoff: ',2x,F7.5)
 1180 format(' Annotation: ',2x,a60)
*
      anno1 = ' '
      read(43,1010) (lookat(i),i=1,3)
      read(43,1020) vtheta,vphi
      read(43,1030) viewd
      read(43,1040) upang
      read(43,1050) scrht
      read(43,1060) near,far
      read(43,1070) nrfade,frfade
      read(43,1080) delta
      read(43,1100) numshl
*
      do 100 i = 1,numshl
        read(43,1110) dummy
        read(43,1120) shells(i),opac(i)
        read(43,1130) (hlscol(j,i),j=1,3)
        call hlsrgb(float(hlscol(1,i)),float(hlscol(2,i))/32767.0,
     &                float(hlscol(3,i))/32767.0,
     &                rgbcol(1,i),rgbcol(2,i),rgbcol(3,i))
 100  continue
*
      read(43,1140) ltheta,lphi
      read(43,1150) reflct(1),reflct(2)
      read(43,1160) reflct(3),reflct(4)
      read(43,1170) thresh
      read(43,1180,end=200) anno1
 200  anno = anno1
      call cvtstr(anno,60)
*
      close(43)
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
c	debug unit(6),subchk
c	end debug

      SUBROUTINE setshl
*  set shell levels and colors
*
      IMPLICIT NONE
      INCLUDE (newvrcm2)
      REAL*4    temp(5),drange
      INTEGER*4 option, i, j, k, l, m
*
      DATA l/0/
	integer*4 bogus(100)
*
	continue
  42  continue
      INCLUDE (cmsclr) 1
 100  format(1x,'Data minimum value: ',e10.4,3x,'maximum value: ',e10.4,
     &    //,1x,' #',3x,'   value  ',5x,'opac',8x,'hue',3x,'light',
     &       3x,' sat ')
 105  format(1x,2('-'),3x,10('-'),4x,7('-'),6x,3('-'),3x,5('-'),3x,
     &       5('-'))
 110  format(1x,I2,3x,E10.4,5x,F5.3,7x,I3,3x,I5,3x,I5)
      write(*,100)datmin,datmax
      write(*,105)
      do 150 i = l+1, min(10+l,numshl)
        write(*,110)i,shells(i),opac(i),hlscol(1,i),hlscol(2,i),
     &              hlscol(3,i)
 150  continue
      write(*,105)
*
      print*,'0> Return to main menu'
      print*,'1> Linear auto shell calculation'
      if (numshl.lt.nshlv) print*,'2> Add shell'
      if (numshl.gt.0) then
        print*,'3> Remove shell'
        print*,'4> Change shell'
        print*,'5> Change opac'
        print*,'6> Change color'
        if (l.gt.0) print*,'7> Previous page'
        if (numshl.gt.l+10) print*,'8> Next page'
      endif
*
 1100 continue
      rewind(5)
      read(5,*,err=1100,end=1100) option
*
*  quit
      if (option.eq.0) then
        goto 3000
*
*  auto shell calculation at 10%, 30%, 50%, 70%,and 90% of data range
      else if (option.eq.1) then
          numshl = 5
          drange = datmax - datmin
          do 1150 i=1,numshl
             shells(i) = datmin + drange * (float(i-1)*0.2 + 0.1)
             opac(i) = float(i)*0.1
             if (i .eq. 1) then
                hlscol(1,i) = 0
             else
                hlscol(1,i) = 360 - (i-1)*60
             endif
             hlscol(2,i) = 16000
             hlscol(3,i) = 32000
             call hlsrgb(float(hlscol(1,i)),float(hlscol(2,i))/32767.0,
     +                   float(hlscol(3,i))/32767.0,
     +                   rgbcol(1,i),rgbcol(2,i),rgbcol(3,i))
 1150     continue
          goto 42
*
*  new shell
      else if (option.eq.2) then
 1200   continue
        rewind(5)
        print*,'enter new shell value:'
        read(5,*,err=1200,end=1200)temp(1)
*
        if (numshl .eq. 0) then
           shells(1) = temp(1)
           numshl = 1
        else
           do 1230 i=1,numshl
              if (temp(1) .eq. shells(i)) then
                 print*,'shell value already exists.  enter new value:'
                 goto 1200
              else if (temp(1) .lt. shells(i)) then
                 do 1220 j=numshl+1,i+1,-1
                    shells(j) = shells(j-1)
                    opac(j) = opac(j-1)
                    do 1210 k=1,3
                       hlscol(k,j) = hlscol(k,j-1)
                       rgbcol(k,j) = rgbcol(k,j-1)
 1210               continue
 1220            continue
                 shells(i) = temp(1)
                 numshl = numshl + 1
	print *,'got thru newshell'
                 goto 42
              endif
 1230      continue
*
           numshl = numshl + 1
           shells(numshl) = temp(1)
        endif
        goto 42
*
*  remove shell
      else if ((option.eq.3).and.(numshl.ge.1)) then
 1300   continue
        rewind(5)
        print*,'remove which shell # (0 to cancel):'
        read(5,*,err=1300,end=1300)m
        if (m.le.0) then
          goto 42
        else if (m.gt.numshl) then
          print*,'Shell # out of range, try again'
          goto 1300
        else if (m .eq. numshl) then
          shells(m) = 0.
          opac(m) = 0.
          do 1310 k=1,3
            hlscol(k,m) = 0.
            rgbcol(k,m) = 0.
 1310     continue
        else
           do 1330 j = m, numshl-1
             shells(j) = shells(j+1)
             opac(j) = opac(j+1)
             do 1320 k = 1,3
               hlscol(k,j) = hlscol(k,j+1)
               rgbcol(k,j) = rgbcol(k,j+1)
 1320        continue
 1330      continue
        endif
        numshl = numshl - 1
        if ((l.eq.numshl).and.(l.ne.0)) l = l - 10
        goto 42
*
*  change shell value
      else if ((option.eq.4).and.(numshl.ge.1)) then
 1400   continue
        rewind(5)
        print*,'change value of which level#? (0 to quit)'
        read(5,*,err=1400,end=1400)m
        if (m .gt. numshl) then
           print*,'shell value out of range.'
           goto 1400
        else if (m .gt. 0) then
 1410      continue
           rewind(5)
           print*,'enter new shell value'
           read(5,*,err=1410,end=1410)temp(1)
           do 1420 i=1,numshl
              if (i .eq. m) goto 1420
              if (temp(1) .eq. shells(i)) then
                 print*,'shell value already exists.'
                 goto 1400
              endif
 1420      continue
           shells(m) = temp(1)
        endif
        goto 42
*
*  change opac
      else if (option.eq.5) then
 1500   continue
        rewind(5)
        print*,'change opac of which level#?'
        read(5,*,err=1500,end=1500)m
        if (m*(m-numshl).le.0) then
 1510     continue
          rewind(5)
          print*,'enter new opac:'
          read(5,*,err=1510,end=1510)opac(m)
        endif
        goto 42
*
*  change color
      else if (option.eq.6) then
 1600   continue
        rewind(5)
        print*,'change color of which level #?'
        read(5,*,err=1600,end=1600)m
*
        if (m*(m-numshl).le.0) then
 1610     continue
          rewind(5)
          print*,'enter new color (hue,lightness,saturation)'
          read(5,*,err=1610,end=1610)hlscol(1,m),hlscol(2,m),hlscol(3,m)
*
          call hlsrgb(float(hlscol(1,m)), float(hlscol(2,m))/32767.0,
     &                  float(hlscol(3,m))/32767.0,
     &                  rgbcol(1,m),rgbcol(2,m),rgbcol(3,m))
        endif
        goto 42
*
*  previous page
      else if ((option.eq.7).and.(l.gt.0)) then
        l = l - 10
        goto 42
*
*  next page
      else if ((option.eq.8).and.(numshl.gt.l+10)) then
        l = l + 10
        goto 42
*
      else
        goto 42
      endif
*
*  exit
 3000 continue
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
      SUBROUTINE hlsrgb(hue, light, sat,
*  convert hls to rgb
     &                    red, green, blue)
*
      IMPLICIT NONE
      REAL*4 hue, light, sat
      REAL*4 red, green, blue
      REAL*4 m1,m2,rgbval
*
      if (light.le.0.5) then
        m2 = light * (1.0 + sat)
      else
        m2 = light + sat - light*sat
      endif
      m1 = 2.0 * light - m2
*
      if (sat.eq.0.0) then
        red   = 255.0 * light
        green = red
        blue  = green
      else
        red   = 255.0 * rgbval(m1,m2,hue)
        green = 255.0 * rgbval(m1,m2,hue-120.0)
        blue  = 255.0 * rgbval(m1,m2,hue+120.0)
      endif
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
      REAL FUNCTION rgbval(n1,n2,hue)
      IMPLICIT NONE
      REAL*4 n1,n2,hue,h
*
      if (hue.gt.360.0) then
        h = hue - 360.0
      else if (hue.lt.0.0) then
        h = hue + 360.0
      else
        h = hue
      endif
*
      if (h.lt.60.0) then
chmc    rgbval = n1 + (n2-n1)*h/60.0
        rgbval = n1 + (n2-n1)*h*0.01666667
      else if (h.lt.180) then
        rgbval = n2
      else if (h.lt.240.0) then
chmc    rgbval = n1 + (n2-n1)*(240.0-h)/60
        rgbval = n1 + (n2-n1)*(240.0-h) * 0.01666667
      else
        rgbval = n1
      endif
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
@PROCESS DC(raydat)
************************************************************************
*  A vectorized subroutine to volume render isovalue surfaces of a
*   three-dimensional function
*
*            Written by Daniel Kartch
*
*
*        Version 6.5              September 6, 1988
***********************************************************************
      SUBROUTINE vrrays
*
      IMPLICIT NONE
      INCLUDE (newvrcm2)
*
*  vcosth, vsinth, vcosph, vsinph are the sines and cosines of the
*   viewing angles
*
*  ucos and usin are the sine and cosine of the up angle
*
*  view contains the view vectors for each screen pixel
*
*  oldval contains the previous function values for each pixel
*  value contains the current function values for each pixel
*
      REAL*4     vcosth, vsinth, vcosph, vsinph, ucos, usin
      REAL*4     view(3, npixel*npixel), vmag, rvmag, rhmag
      REAL*4     oldval(npixel*npixel), value(npixel*npixel)
      REAL*4     rayop(npixel*npixel), raycol(3, npixel*npixel)
      REAL*4     grad1, grad2, grad3
      REAL*4     pixht, datprm
      REAL*4     eps
chmc  REAL*4     op, adj1, adj2, gmag, gdotl, gdoth, gdotv, hmag
      REAL*4     op, adj1, adj2, gmag, rgmag, gdotl,gdoth,gdotv,hmag
      REAL*4     d, u1, u2
      PARAMETER(eps = 1.0E-50)
*
      INTEGER*4  raycrs(npixel*npixel)
      INTEGER*4  hnpix, irc
      INTEGER*4  i, j, l, m, n, shl1, shl2, inc, ofs
*
      COMMON /RAYDAT/ VIEW,OLDVAL,VALUE,RAYOP,RAYCOL,RAYCRS
*
*  turn off underflow interrupt-handling
      call xuflow(0)
      datprm = 0.0
*
*  calculate sines and cosines of viewing and up angles
      vcosth = cos(vtheta)
      vsinth = sin(vtheta)
      vcosph = cos(vphi)
      vsinph = sin(vphi)
      ucos   = cos(upang)
      usin   = sin(upang)
*
*  determine eye point
      eyepnt(1) = lookat(1) + viewd * vcosph * vsinth
      eyepnt(2) = lookat(2) + viewd * vsinph * vsinth
      eyepnt(3) = lookat(3) + viewd * vcosth
*
*  determine up vector
      upvect(1) = vcosph * ucos * vcosth - usin * vsinph
      upvect(2) = vsinph * ucos * vcosth + usin * vcosph
      upvect(3) = -ucos * vsinth
*
*  determine horizontal vector
      hvect(1) =  ucos * vsinph + usin * vcosph * vcosth
      hvect(2) = -ucos * vcosph + usin * vsinph * vcosth
      hvect(3) = -usin * vsinth
*
*  determine light vector
      light(1) = -cos(lphi) * sin(ltheta)
      light(2) = -sin(lphi) * sin(ltheta)
      light(3) = -cos(ltheta)
*
*
      hnpix = npixel/2
*
*  initialize opacities and colors of the rays
CHMC REPLACE THIS ZEROING LOOP WITH CALLS TO ZEROST?
      do 110 n = 1, npixel*npixel
        rayop(n) = 1.0
        raycol(1,n) = 0.0
        raycol(2,n) = 0.0
        raycol(3,n) = 0.0
 110  continue
chmc  irc = 0
chmc  call zerost(raycol,3*npixel*npixel,8,irc)
chmc  if (irc .ne. 0) stop 'Problem zeroing RAYCOL'
*
*  calculate view and light vectors and find initial ray points for all
*   rays
      pixht = scrht/float(npixel)
chmc this should be a bit faster because it avoids divisions
chmc within the loops
chmc  rnpixel = FLOAT(1.D0 / npixel)
      do 300 n = 1, npixel*npixel
        m = (n-1)/npixel + 1
chmc    m = ((n-1) * rnpixel) + 1
        l = n - (m-1)*npixel
*
*  calculate view vectors
        vmag  = 0.0
        do 230 i = 1, 3
          view(i,n) = lookat(i) + pixht * (float(m-hnpix)
chmc substitute rasprt for asprat?
     &              * upvect(i) + hvect(i) * float(l-hnpix)/asprat)
     &              - eyepnt(i)
          vmag = vmag + view(i,n)*view(i,n)
 230    continue
        vmag  = sqrt(vmag)
        rvmag = 1.E0 / vmag
        do 240 i = 1, 3
chmc      view(i,n) = view(i,n)/vmag
          view(i,n) = view(i,n) * rvmag
  240   continue
  300 continue
*
*  determine initial function value for all rays
      call vrfunc(view, viewd+near, value, npixel, datprm)
*
*
*  simultaneously fire rays for the entire screen
chmc for the present code, viewd = 100
chmc                       near  = -50
chmc                       far   =  50
chmc                       delta =   1
chmc so this loop:  51, 150,1
      do 2000 d = viewd + near + delta, viewd + far, delta
*
*  save values for previous set of ray points
        do 1110 n = 1, npixel*npixel
          oldval(n) = value(n)
 1110   continue
*
*  determine function value at current ray positions
        call vrfunc(view, d, value, npixel, datprm)
*
*  see if any rays have crossed a shell and if so, adjust their colors
*   accordingly
*
        do 1620 n = 1, npixel*npixel
          raycrs(n) = 0
          do 1610 j = 1, numshl
            if((shells(j)-oldval(n))*(shells(j)-value(n)).le.0.0) then
              if (raycrs(n).eq.0) then
                raycrs(n) = 101*j
              else
                raycrs(n) = raycrs(n) + 1
              endif
            endif
 1610     continue
          if (oldval(n).gt.value(n)) raycrs(n) = -raycrs(n)
 1620   continue
*
        do 1750 n = 1, npixel*npixel
          inc = sign(1, raycrs(n))
          raycrs(n) = abs(raycrs(n))
          if (raycrs(n).eq.0) then
            shl1 = numshl
            shl2 = 0
          else if (inc.gt.0) then
chmc        shl1 = raycrs(n)/100
            shl1 = raycrs(n) * 0.01
            shl2 = raycrs(n) - shl1*100
          else
chmc        shl2 = raycrs(n)/100
            shl2 = raycrs(n) * 0.01
            shl1 = raycrs(n) - shl2*100
          endif
*
          do 1740 j = shl1, shl2, inc
*
          if (rayop(n).gt.thresh) then
            u2 = (shells(j) - oldval(n))/(value(n)-oldval(n))
            u1 = 1.0 - u2
            call vrgrad(view(1,n), d - u1*delta, datprm,
     &                  grad1, grad2, grad3)
chmc        gmag = sqrt(grad1*grad1 + grad2*grad2 + grad3*grad3) + eps
            rgmag = 1.E0/sqrt(grad1*grad1+grad2*grad2+grad3*grad3)+eps
chmc        grad1 = grad1/gmag
chmc        grad2 = grad2/gmag
chmc        grad3 = grad3/gmag
            grad1 = grad1 * rgmag
            grad2 = grad2 * rgmag
            grad3 = grad3 * rgmag
*
            gdotl = abs(grad1*light(1)+grad2*light(2)+grad3*light(3))
*
chmc        hmag = sqrt((view(1,n)+light(1))**2+(view(2,n)+light(2))**2
chmc &                  + (view(3,n)+light(3))**2)
*
            rhmag = 1.E0/sqrt((view(1,n)+light(1))*(view(1,n)+light(1))
     &                      + (view(2,n)+light(2))*(view(2,n)+light(2))
     &                      + (view(3,n)+light(3))*(view(3,n)+light(3)))
*
            gdoth = abs (grad1*(view(1,n)+light(1))
     &                 + grad2*(view(2,n)+light(2))
     &                 + grad3*(view(3,n)+light(3)) )  *  rhmag
*
            gdotv = abs(grad1*view(1,n) + grad2*view(2,n)
     &                + grad3*view(3,n))
*
            op = 1.0 - (1.0 - opac(j))**(1.0/(gdotv + eps))
            if ((d-viewd-u1*delta).lt.nrfade) then
              op = op * ((d-viewd-u1*delta)-near)/(nrfade - near)
            else if ((d-viewd-u1*delta).gt.frfade) then
              op = op * ((d-viewd-u1*delta)-far)/(frfade - far)
            endif
            adj1 = rayop(n)*op * (reflct(1) + reflct(2) * gdotl)
            adj2 = rayop(n)*op * 255.0*reflct(3) * gdoth**reflct(4)
*
            do 1730 i = 1, 3
              raycol(i,n) = raycol(i,n) + adj1 * rgbcol(i,j) + adj2
 1730       continue
            rayop(n) = rayop(n) * (1.0 - op)
          endif
*
 1740     continue
 1750   continue
*
 2000 continue
*
      do 3010 n = 1, npixel*npixel
        do 3000 i = 1, 3
          if (raycol(i,n).lt.255.0) then
            screen(i,n) = int(raycol(i,n))
          else
            screen(i,n) = 255
          endif
 3000   continue
 3010 continue
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
      SUBROUTINE startit
*
      IMPLICIT NONE
      INCLUDE (newvrcm2)
      REAL*4 rasprt
      INTEGER*2 ix1,ix2
*
*  turn off raster tech if it has already been initialized
*
      INCLUDE (cmsdev6) 1
*
10    continue
      rewind(5)
      print*
      print*,'Select an output device:'
      print*,'   1> Raster Tech (CMS only)'
      print*,'   2> Run-length encoded bitmap (640x480)'
      print*,'   3> Run-length encoded bitmap (1280x1024)'
      print*,'   4> Run-length encoded bitmap (1024x1024)'
      read(5,*,err=10,end=10)dparm
      INCLUDE (cmsdev1) 1
      INCLUDE (aixdev1) 2
      INCLUDE (mvsdev1) 3
      if (dparm .eq. 2. .or. dparm .eq. 3.) then
         asprat    = 0.8
         rasprt  = 1.25
      endif
      if (dparm .eq. 4.) then
         asprat    = 1.E0
         rasprt  = 1.E0
      endif
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,pixs)
      SUBROUTINE vrndpl(rgb, npixel, asprat)
*
      IMPLICIT NONE
      REAL*4 dparm, asprat, view(3), lookat(3), upvect(3), scrht
      REAL*4    uh1a,uh1b,uh2a,uh2b,uv1,uv2,rh,rv
      INTEGER*4 npixel,nsize,nrows,ncols,fixit
      INTEGER*4 i,j,k1,k2,k,la,lb,m,pict,irow
      INTEGER*4 rda,rdb,bla,blb,gra,grb
      PARAMETER (nsize=655360)
      INTEGER*2 ix1,ix2,iy1,iy2,izero,rlen,glen,blen
      INTEGER*2 rgb(3,npixel,npixel),nr,nc
      INTEGER*2 red(nsize),grn(nsize),blu(nsize)
      INTEGER*2 encred(1280),encgrn(1280),encblu(1280)
      CHARACTER*60 anno,anno1
      CHARACTER*20 wfmt
      COMMON /parm1/ dparm,anno,anno1
      COMMON /pixs/ red,grn,blu
*
      INCLUDE (cmsdev2) 1
      IF (dparm .eq. 2.) THEN
         nrows = 480
         ncols = 640
         irow = nrows
         INCLUDE (aixdev3) 2
      ENDIF
      IF (dparm .eq. 3.) THEN
         nrows = 1024
         ncols = 1280
         irow = nrows
         INCLUDE (aixdev4) 2
      ENDIF
      IF (dparm .eq. 4.) THEN
         nrows = 1024
         ncols = 1024
         irow = nrows
         INCLUDE (aixdev5) 2
      ENDIF
*
      rh = float(npixel)/float(ncols)
      rv = float(npixel)/float(nrows)
      do 150 j=1,irow
         m = npixel - (j-1)*rv
         uv2 = float(npixel) - (float(j)-1.0)*rv - float(m)
         uv1 = 1.0 - uv2
         do 140 i = 1,ncols,2
            la = (i-1)*rh + 1
            lb = i*rh + 1
            uh2a = (float(i)-1.0)*rh +1.0 - float(la)
            uh1a = 1.0 - uh2a
            uh2b = float(i)*rh +1.0 - float(lb)
            uh1b = 1.0 - uh2b
            k = ishft(((j-1)*ncols + i + 1),-1)
c
            rda=int(uh1a*uv1*rgb(1,la,m)+uh2a*uv1*rgb(1,la+1,m)
     &          +uh1a*uv2*rgb(1,la,m+1)+uh2a*uv2*rgb(1,la+1,m+1))
            rdb=int(uh1b*uv1*rgb(1,lb,m)+uh2b*uv1*rgb(1,lb+1,m)
     &          +uh1b*uv2*rgb(1,lb,m+1)+uh2b*uv2*rgb(1,lb+1,m+1))
            red(k) = ishft(fixit(rda),8) + fixit(rdb)
c
            gra=int(uh1a*uv1*rgb(2,la,m)+uh2a*uv1*rgb(2,la+1,m)
     &          +uh1a*uv2*rgb(2,la,m+1)+uh2a*uv2*rgb(2,la+1,m+1))
            grb=int(uh1b*uv1*rgb(2,lb,m)+uh2b*uv1*rgb(2,lb+1,m)
     &          +uh1b*uv2*rgb(2,lb,m+1)+uh2b*uv2*rgb(2,lb+1,m+1))
            grn(k) = ishft(fixit(gra),8) + fixit(grb)
c
            bla=int(uh1a*uv1*rgb(3,la,m)+uh2a*uv1*rgb(3,la+1,m)
     &          +uh1a*uv2*rgb(3,la,m+1)+uh2a*uv2*rgb(3,la+1,m+1))
            blb=int(uh1b*uv1*rgb(3,lb,m)+uh2b*uv1*rgb(3,lb+1,m)
     &          +uh1b*uv2*rgb(3,lb,m+1)+uh2b*uv2*rgb(3,lb+2,m+1))
            blu(k) = ishft(fixit(bla),8) + fixit(blb)
c
  140    continue
         IF (dparm .ge. 2. .and. dparm .le. 4.) THEN
            INCLUDE (cmsdev3) 1
            INCLUDE (aixdev2) 2
            INCLUDE (mvsdev2) 3
         ENDIF
  150 continue
      INCLUDE (cmsdev4) 1
*
      RETURN
      END
@PROCESS OPT(3)
      integer function fixit(ivar)
      integer*4 ivar
*
      if (ivar .lt. 0) then
         fixit = 0
      else if (ivar .gt. 255) then
         fixit = 255
      else
         fixit = ivar
      endif
      return
      end
c***********************************************************************
      subroutine encode(rawval,encval,enclen,length)
      implicit none
      integer*4 icount,curpos,icase,length,i
      logical qval,qcount,qloop
      integer*2 enclen
      character*1 rawval(length),encval(length*2),oldval
c
c     enclen returns the number of bytes in the encoded string
c     length is the length of each scan line
c
      curpos = 0
      icount = 0
      oldval = rawval(1)
c
      do 100 i=1,length
         qval = rawval(i) .eq. oldval
         qcount = icount .eq. 256
         qloop = i .eq. length
*
         if      (      qval .and. .not. qcount .and. .not. qloop) then
            icount = icount + 1
*
         else if (      qval .and. .not. qcount .and.       qloop) then
            curpos = curpos + 1
            encval(curpos) = char(icount)
            curpos = curpos + 1
            encval(curpos) = oldval
*
         else if (.not. qloop) then
            curpos = curpos + 1
            if (icount .eq. 0) icount = 1
            encval(curpos) = char(icount-1)
            curpos = curpos + 1
            encval(curpos) = oldval
            icount = 1
            oldval = rawval(i)
*
         else if (qloop) then
            curpos = curpos + 1
            if (icount .eq. 0) icount = 1
            encval(curpos) = char(icount-1)
            curpos = curpos + 1
            encval(curpos) = oldval
*
* encode that last bugger
*
            curpos = curpos + 1
            encval(curpos) = char(0)
            curpos = curpos + 1
            encval(curpos) = oldval
         endif
100   continue
      enclen = curpos
      return
      end
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
      SUBROUTINE getout
      IMPLICIT NONE
      INCLUDE (newvrcm2)
      INTEGER*2 izero
*
      INCLUDE (cmsdev5) 1
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(datblk)
      SUBROUTINE vrinit(datmin,datmax)
*  load data array
*
      IMPLICIT NONE
      INCLUDE (datacom2)
*
      REAL*4 datmin,datmax
      INTEGER*4 i,j,k,n,xdum,ydum,zdum,iform
      CHARACTER fname*80,frmstr*80
      LOGICAL qexist
*
*******************************************************************
****  modify this section for particular method of array storage
***
**
*
*
      INCLUDE (cmsio6) 1
      INCLUDE (aixio6) 2
      INCLUDE (mvsio6) 3
      if (.not. qexist) then
         print*,'Data file not found!'
         goto 10
      endif
c
15    print *,' 1> Binary data file'
      print *,' 2> Formatted data'
      read(5,*,err=15,end=15)iform
c
      if (iform .eq. 1) then
         INCLUDE (cmsio7) 1
         INCLUDE (aixio7) 2
         INCLUDE (mvsio7) 3
         read(2,err=20)
     +         (((datary(i,j,k),i=0,xsize-1),j=0,ysize-1),k=0,zsize-1)
         close(2)
      else if (iform .eq. 2) then
         print *,' Enter format string (including parentheses):'
         read(5,'(a)')frmstr
         INCLUDE (cmsio8) 1
         INCLUDE (aixio8) 2
         INCLUDE (mvsio8) 3
         read(2,frmstr,err=20)
     +         (((datary(i,j,k),i=0,xsize-1),j=0,ysize-1),k=0,zsize-1)
         close(2)
      else
         goto 15
      endif
      goto 30
*
20    continue
      print*,'Error reading data file'
      goto 10
*
**
***
****
*********************************************************************
*
30    continue
*
*     find data min and max (will be used on shell value submenu)
*
      datmin = datary(0,0,0)
      datmax = datary(0,0,0)
      do 120 k=0,zsize-1
         do 110 j=0,ysize-1
            do 100 i=0,xsize-1
               datmin = amin1(datary(i,j,k),datmin)
               datmax = amax1(datary(i,j,k),datmax)
  100       continue
  110    continue
  120 continue
*
      do 1100 k=0,zsize-1
        do 1000 j=0,ysize-1
          datary(-1,j,k) = datary(0,j,k)
          datary(xsize,j,k) = datary(xsize-1,j,k)
 1000   continue
 1100 continue
*
      do 2100 k = 0,zsize-1
        do 2000 i = -1,xsize
          datary(i,-1,k) = datary(i,0,k)
          datary(i,ysize,k) = datary(i,ysize-1,k)
 2000   continue
 2100 continue
*
      do 3100 j = -1,ysize
        do 3000 i = -1,xsize
          datary(i,j,-1) = datary(i,j,0)
          datary(i,j,zsize) = datary(i,j,zsize-1)
 3000   continue
 3100 continue
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
@PROCESS DC(datblk, edblk)
      SUBROUTINE vrfunc(view, d, values, npts, reset)
*  calculate value at a given set of n*n points using triinear interp.
*
      IMPLICIT NONE
*
      INCLUDE (newvrcm2)
      INCLUDE (datacom2)
*
      INTEGER*4 npts
      INTEGER*4 n,a,b,c
*
      REAL*4 view(3,npts*npts),d
      REAL*4 reset, values(npts*npts)
      REAL*4 nredge(npixel*npixel), fredge(npixel*npixel)
      REAL*4 t,x,y,z,u1,u2,v1,v2,w1,w2
*
      COMMON/edblk/nredge,fredge
*
      if (reset.eq.0.0) then
        call edset(view, npts*npts)
        reset=1.0
      endif
*
      do 100 n = 1, npts*npts
        t = d
        if (t.lt.nredge(n)) t = nredge(n)
        if (t.gt.fredge(n)) t = fredge(n)
*
        x = eyepnt(1) + t*view(1,n)
        y = eyepnt(2) + t*view(2,n)
        z = eyepnt(3) + t*view(3,n)
*
        if ((x*(x-xsize+1).gt.0.0).or.(y*(y-ysize+1).gt.0.0).or.
     &      (z*(z-zsize+1).gt.0.0)) then
          x = 0.0
          y = 0.0
          z = 0.0
        endif
*
        a  = int(x)
        u2 = x - a
        u1 = 1.0 - u2
        b  = int(y)
        v2 = y - b
        v1 = 1.0 - v2
        c  = int(z)
        w2 = z - c
        w1 = 1.0 - w2
*
        values(n) = w1 * (v1 * (u1 * datary(a,b,c)
     &            +             u2 * datary(a+1,b,c))
     &            +       v2 * (u1 * datary(a,b+1,c)
     &            +             u2 * datary(a+1,b+1,c)))
     &            + w2 * (v1 * (u1 * datary(a,b,c+1)
     &            +             u2 * datary(a+1,b,c+1))
     &            +       v2 * (u1 * datary(a,b+1,c+1)
     &            +             u2 * datary(a+1,b+1,c+1)))
*
  100 continue
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
@PROCESS DC(datblk, edblk)
      SUBROUTINE edset(view, npts)
*  set near and far ends of rays
*
      IMPLICIT NONE
*
      INCLUDE (newvrcm2)
      INCLUDE (datacom2)
*
      INTEGER*4 npts
      INTEGER*4  i,j,n
*
      REAL*4 view(3,npts)
      REAL*4     nredge(npixel*npixel), fredge(npixel*npixel)
      REAL*4     t,edge(6),x,y,z
      REAL*4     plndst
*
      COMMON/edblk/nredge,fredge
*
      do 500 n = 1, npts
        edge(1) = plndst(eyepnt,view(1,n),1,0)
        edge(4) = plndst(eyepnt,view(1,n),1,xsize-1)
        edge(2) = plndst(eyepnt,view(1,n),2,0)
        edge(5) = plndst(eyepnt,view(1,n),2,ysize-1)
        edge(3) = plndst(eyepnt,view(1,n),3,0)
        edge(6) = plndst(eyepnt,view(1,n),3,zsize-1)
*
        do 400 i = 1, 5
          do 300 j = 6, i+1, -1
            if (edge(i).gt.edge(j)) then
              t       = edge(j)
              edge(j) = edge(i)
              edge(i) = t
            endif
 300      continue
 400    continue
*
        nredge(n) = edge(3) * 1.0001
        fredge(n) = edge(4) * 0.9999
*
 500  continue
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
      REAL FUNCTION plndst(point, vector, axis, bound)
*
      REAL*4     point(3), vector(3)
      INTEGER*4  axis, bound
*
      if (vector(axis).ne.0.0) then
        plndst = (float(bound) - point(axis))/vector(axis)
      else if (bound.eq.0) then
        plndst = 100000000.0
      else
        plndst = -100000000.0
      endif
*
      RETURN
      END
@PROCESS OPT(3) VEC(LEV(2) REP(LIST))
@PROCESS DC(parm1,view1,view2,light1,scrn1,ray1,shell1,film1,film2)
@PROCESS DC(datblk)
      SUBROUTINE vrgrad(view, d, dummy, gradx, grady, gradz)
*  calculate the gradient at a point
*
      IMPLICIT NONE
*
      INCLUDE (newvrcm2)
      INCLUDE (datacom2)
*
      REAL*4 view(3), d, dummy
      REAL*4 gradx, grady, gradz
      REAL*4 x,y,z,u1,u2,v1,v2,w1,w2,hu1,hu2,hv1,hv2,hw1,hw2
*
      INTEGER*4 a,b,c,ha,hb,hc
*
      x = eyepnt(1) + d * view(1)
      y = eyepnt(2) + d * view(2)
      z = eyepnt(3) + d * view(3)
*
      ha = int(x - .5)
      hu2 = (x - .5) - ha
      hu1 = 1.0 - hu2
      hb = int(y - .5)
      hv2 = (y - .5) - hb
      hv1 = 1.0 - hv2
      hc  = int(z - .5)
      hw2 = (z - .5) - hc
      hw1 = 1.0 - hw2
*
      a  = int(x)
      u2 = x - a
      u1 = 1.0 - u2
      b  = int(y)
      v2 = y - b
      v1 = 1.0 - v2
      c  = int(z)
      w2 = z - c
      w1 = 1.0 - w2
*
      gradx = w1 * (v1 * (hu1 * (-datary(ha,b,c)
     &      +                     datary(ha+1,b,c))
     &      +             hu2 * (-datary(ha+1,b,c)
     &      +                     datary(ha+2,b,c)))
     &      +       v2 * (hu1 * (-datary(ha,b+1,c)
     &      +                     datary(ha+1,b+1,c))
     &      +             hu2 * (-datary(ha+1,b+1,c)
     &      +                     datary(ha+2,b+1,c))))
     &      + w2 * (v1 * (hu1 * (-datary(ha,b,c+1)
     &      +                     datary(ha+1,b,c+1))
     &      +             hu2 * (-datary(ha+1,b,c+1)
     &      +                     datary(ha+2,b,c+1)))
     &      +       v2 * (hu1 * (-datary(ha,b+1,c+1)
     &      +                     datary(ha+1,b+1,c+1))
     &      +             hu2 * (-datary(ha+1,b+1,c+1)
     &      +                     datary(ha+2,b+1,c+1))))
*
      grady = w1 * (u1 * (hv1 * (-datary(a,hb,c)
     &      +                     datary(a,hb+1,c))
     &      +             hv2 * (-datary(a,hb+1,c)
     &      +                     datary(a,hb+2,c)))
     &      +       u2 * (hv1 * (-datary(a+1,hb,c)
     &      +                     datary(a+1,hb+1,c))
     &      +             hv2 * (-datary(a+1,hb+1,c)
     &      +                     datary(a+1,hb+2,c))))
     &      + w2 * (u1 * (hv1 * (-datary(a,hb,c+1)
     &      +                     datary(a,hb+1,c+1))
     &      +             hv2 * (-datary(a,hb+1,c+1)
     &      +                     datary(a,hb+2,c+1)))
     &      +       u2 * (hv1 * (-datary(a+1,hb,c+1)
     &      +                     datary(a+1,hb+1,c+1))
     &      +             hv2 * (-datary(a+1,hb+1,c+1)
     &      +                     datary(a+1,hb+2,c+1))))
*
      gradz = v1 * (u1 * (hw1 * (-datary(a,b,hc)
     &      +                     datary(a,b,hc+1))
     &      +             hw2 * (-datary(a,b,hc+1)
     &      +                     datary(a,b,hc+2)))
     &      +       u2 * (hw1 * (-datary(a+1,b,hc)
     &      +                     datary(a+1,b,hc+1))
     &      +             hw2 * (-datary(a+1,b,hc+1)
     &      +                     datary(a+1,b,hc+2))))
     &      + v2 * (u1 * (hw1 * (-datary(a,b+1,hc)
     &      +                     datary(a,b+1,hc+1))
     &      +             hw2 * (-datary(a,b+1,hc+1)
     &      +                     datary(a,b+1,hc+2)))
     &      +       u2 * (hw1 * (-datary(a+1,b+1,hc)
     &      +                     datary(a+1,b+1,hc+1))
     &      +             hw2 * (-datary(a+1,b+1,hc+1)
     &      +                     datary(a+1,b+1,hc+2))))
*
      RETURN
      END
C****************************************************************
@PROCESS OPT(3)
      SUBROUTINE WRITIM(LUNITO,TIME1,TIME2)
      IMPLICIT NONE
C This routine simply writes timing output to logical unit, LUNITO.
      REAL*8      TIME1(4), TIME2(4)
      INTEGER*4   LUNITO
C
      WRITE(LUNITO,*) ' '
C
      WRITE(LUNITO,10) TIME2(1)-TIME1(1)
   10 FORMAT(1X,'Wallclock time: ',F13.3)
C
      WRITE(LUNITO,20) TIME2(3)-TIME1(3)
   20 FORMAT(1X,'Total CPU time: ',F13.3)
C
      WRITE(LUNITO,30) TIME2(2)-TIME1(2)
   30 FORMAT(1X,'Virt. CPU time: ',F13.3)
C
      WRITE(LUNITO,40) TIME2(4)-TIME1(4)
   40 FORMAT(1X,'Vect. CPU time: ',F13.0)
C
      WRITE(LUNITO,50) 100.0*(TIME2(4)-TIME1(4))/(TIME2(2)-TIME1(2))
   50 FORMAT(1X,'Percent Vector: ',F13.3)
C
      WRITE(LUNITO,*) ' '
      WRITE(LUNITO,*) '****************************************'
C
      RETURN
      END
      subroutine cvtstr(strng,ilen)
      implicit none
      integer ilen,eb2as,i
      character*(*) strng
      character*1 str1,str2
      do 10 i=1,ilen,2
        if (i+1 .le. ilen) then
            str1 = char(eb2as(ichar(strng(i+1:i+1))))
        else
            str1 = ' '
        endif
        str2 = char(eb2as(ichar(strng(i:i))))
        strng(i:i+1) = str1//str2
  10  continue
      return
      end
      INTEGER FUNCTION EB2AS (EBCDIC)
      IMPLICIT NONE
C
C     EB2AS RETURNS THE ASCII REPRESENTATION OF THE EBCDIC
C     CHARACTER CONTAINED IN EBCDIC
C
      INTEGER    EBCDIC
      INTEGER    ICTOAS(255)
      DATA ICTOAS /
     +      1,  2,  3, 22,  9, 22,127, 22, 22, 22, 11, 12, 13, 14, 15,
     + 16, 17, 18, 22, 22, 10,  8, 22, 24, 25, 22, 22, 28, 29, 30, 31,
     + 22, 22, 22, 22, 22, 10, 23, 27, 22, 22, 22, 22, 22,  5,  6,  7,
     + 22, 22, 22, 22, 22, 30, 22,  4, 22, 22, 22, 22, 20, 21, 22, 26,
     + 32, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 46, 60, 40, 43, 94,
     + 38, 22, 22, 22, 22, 22, 22, 22, 22, 22, 33, 36, 42, 41, 59, 94,
     + 45, 47, 22, 22, 22, 22, 22, 22, 22, 22,124, 44, 37, 95, 62, 63,
     + 22, 22, 22, 22, 22, 22, 22, 22, 22, 96, 58, 35, 64, 39, 61, 34,
     + 22, 97, 98, 99,100,101,102,103,104,105, 22,123, 22, 22, 22, 22,
     + 22,106,107,108,109,110,111,112,113,114, 22,125, 22, 22, 22, 22,
     + 22,126,115,116,117,118,119,120,121,122, 22, 22, 22, 91, 22, 22,
     + 94, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 93, 22, 22,
     + 22, 65, 66, 67, 68, 69, 70, 71, 72, 73, 22, 22, 22, 22, 22, 22,
     + 22, 74, 75, 76, 77, 78, 79, 80, 81, 82, 22, 22, 22, 22, 22, 22,
     + 92, 22, 83, 84, 85, 86, 87, 88, 89, 90, 22, 22, 22, 22, 22, 22,
     + 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 22, 22, 22, 22, 22, 22/
C
      EB2AS = ICTOAS(EBCDIC)
      RETURN
      end
