	  subroutine vector(flag,NPLANE,NUM_Ds,nz,ny,nx)
C********************************************************************
C  Write out deformation parameters in vector format. Each parameter
C  will be written to its own file.  If the program is being
C  run non-interactively, all parameters of a particular selected
C  type will be written out.  For example, if the User has
C  specified 'gradients of displacement' in the input file, 
C  9 output file will be generated - one for each element of
C  the displacement gradient tensor. If running in an interactive 
C  mode, the User can select which parameters to be written out 
C  (e.g., only the dUx/dz, dUy/dz,dUz/dz components of the
C  displacement gradient tensor or whatever subset is desired).
C
C  See the description of the X,Y,Z,Data formatted files to 
C  find out what coordinate systems the output parameters are
C  calculated for.
C
C  The following describes the output file possibilities for each
C  output type requested.  Note that for in-plane coordinates,
C  x is the strike direction, y is the dip direction, and z is
C  normal to the plane.
C   Output type		Files		     Contents
C 	element		elemds#.???     strike component of displacements 
C  	displacements	elemdd#.???     dip component of displacements 
C     			elemdn#.???     tensile component of displacements
C					#=index number, 1... no. of planes
C	stress	  	Sxx.dat		xx component of the stress tensor
C		  	Sxy.???		xy component of the stress tensor
C			Sxz.???		xz component of the stress tensor
C			Syy.???		yy component of the stress tensor
C			Syz.???		yz component of the stress tensor
C			Szz.???		zz component of the stress tensor
C	strain	  	Exx.???		xx component of the strain tensor
C			Exy.???		xy component of the strain tensor
C			Exz.???		xz component of the strain tensor
C			Eyy.???		yy component of the strain tensor
C			Eyz.???		yz component of the strain tensor
C			Ezz.???		zz component of the strain tensor
C	gradients in
C	 displacements	dUxdx.???	gradient of x-displacement,  x-direction
C			dUxdy.???	gradient of x-displacement,  y-direction
C			dUxdz.???	gradient of x-displacement,  z-direction
C			dUydx.???	gradient of y-displacement,  x-direction
C			dUydy.???	gradient of y-displacement,  y-direction
C			dUydz.???	gradient of y-displacement,  z-direction
C			dUzdx.???	gradient of z-displacement,  x-direction
C			dUzdy.dat	gradient of z-displacement,  y-direction
C			dUzdz.dat	gradient of z-displacement,  z-direction
C   rigid body 
C	 rotations	rbr_deg_x.???	rotation about the X (E) axis
C			rbr_deg_y.???	rotation about the Y (N) axis
C			rbr_deg_z.???	rotation about the Z (vertical) axis
C	orientation of
C    principal axes	prince_max.???	magnitude of the max. principal strain
C			plunge_max.???	plunge wrt horizontal of the max. principal strain 
C			trend_max.???	trend (cw wrt N) of the max. principal strain
C			prince_int.???	magnitude of the intermediate principal strain
C			plunge_int.???	plunge wrt horizontal of the intermediate principal strain
C			trend_int.???	trend (cw wrt N) of the intermediate principal strain
C			prince_min.???	magnitude of the min. principal strain
C			plunge_min.???	plunge wrt horizontal of the min. principal strain
C			trend_min.???	trend (cw wrt N) of the min. principal strain
C   displacements	disp1.???	displacements in the x-direction.
C			disp2.???	displacements in the y-direction.
C			disp3.???	displacements in the z-direction.
C	invariants	delv.???	volumetric strain
C			aftershocks.???	failure stress
C			workdis.???	distortional strain energy density
C			workdil.???	dilational strain energy density
C       failures        strike1.???     strike of 1st failure plane
C                       dip1.???        dip of 1st failure plane
C                       ce1.???         slip vector dir. cos. wrt X(E)
C                       cn1.???         slip vector dir. cos. wrt Y(N)
C                       cz1.???         slip vector dir. cos. wrt Z(up)
C                       strike2.???     strike of 2nd failure plane
C                       dip2.???        dip of 2nd failure plane
C                       ce2.???         slip vector dir. cos. wrt X(E)
C                       cn2.???         slip vector dir. cos. wrt Y(N)
C                       cz2.???         slip vector dir. cos. wrt Z(up)
C********************************************************************
      USE parameters

      logical o_elem,o_stress,o_strain,o_dgrad,o_rbr,
     &  o_orient,o_disp,o_inv,xyz_out,o_fail,o_proj
      character*4 sufx
      common/outputs/sufx,o_elem,o_stress,o_strain,o_dgrad,
     &  o_rbr,o_orient,o_disp,o_inv,xyz_out,o_fail,o_proj

      logical outparam(10)
      character*1 ok,flag

C---- Re-format all output into vector formatted files.	----  
      if(flag.ne.'i'.and.flag.ne.'I') then	 

           write(*,'(/,a)') ' Writing vector output...'
		
	   do 111 i=1,10
111	   outparam(i)=.true.

	   if(o_elem) then
		
	        call elem_disps(num_ds,nplane,outparam)
		write(*,'(/,3a)') 
     1 '  wrote element files elemds#.'//sufx,
     2 ',elemdd#.'//sufx,',elemdn#.'//sufx

           endif

  	   if(o_fail) then
	   
		rewind(IFAIL)
			
                IOUT=IOUTFILE
		open(unit= IOUT,file='strike1.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='dip1.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='ce1.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='cn1.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='cz1.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='strike2.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='dip2.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='ce2.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='cn2.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='cz2.'//sufx,status='unknown')
	
	 	call vect_out(outparam,10,IFAIL,nx,ny,nz)
			
		write(*,'(/,4a)') 
     1 '  wrote failure files strike1.'//sufx,',dip1.'//sufx,
     2 ',ce1.'//sufx,',cn1.'//sufx
		write(*,'(6a)') 
     1 '  cz1.'//sufx,',strike2.'//sufx,',dip2.'//sufx,
     2   ',ce2.'//sufx,',cn2.'//sufx,',cz2.'//sufx
		close(IFAIL)
		IOUT=IOUTFILE
		do 999 i=1,10
		close(IOUT)
 		IOUT=IOUT+1
999	        continue
            endif

  	    if(o_stress) then

		rewind(ISTRESS)
			
                IOUT=IOUTFILE
		open(unit= IOUT,file='Sxx.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='Sxy.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='Sxz.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='Syy.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='Syz.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='Szz.'//sufx,status='unknown')
C
C  code added to output max shear stress, etc. see manual and TOPLANE for details.
C
		iout=iout+1
		open(unit=iout,file='tmax.'//sufx,status='unknown')
		iout=iout+1
		open(unit=iout,file='tmaxa.'//sufx,status='unknown')
		iout=iout+1
		open(unit=iout,file='tmx.'//sufx,status='unknown')
		iout=iout+1
		open(unit=iout,file='tmy.'//sufx,status='unknown')
C
C
	
	 	call vect_out(outparam,10,ISTRESS,nx,ny,nz)
			
		write(*,'(/,3a)') 
     1 '  wrote stress files Sxx.'//sufx,',Sxy.'//sufx,
     2 ',Sxz.'//sufx
		write(*,'(3a)') 
     1 '                     Syy.'//sufx,',Syz.'//sufx,',Szz.'//sufx
		write(*,'(/,3a)') 
     1 '  wrote max shear stress files tmax.'//sufx,',tmx.'//sufx,
     2 ',tmy.'//sufx

		close(ISTRESS)
		IOUT=IOUTFILE
		do 888 i=1,9
		close(IOUT)
 		IOUT=IOUT+1
888	        continue
	    endif		
	    
	    if(o_strain) then

		rewind(ISTRAIN)
	 		
	        IOUT=IOUTFILE
		open(unit= IOUT,file='Exx.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='Exy.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='Exz.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='Eyy.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='Eyz.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='Ezz.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='Evol.'//sufx,status='unknown')

	 	call vect_out(outparam,7,ISTRAIN,nx,ny,nz)

		write(*,'(/,3a)') 
     1 '  wrote strain files Exx.'//sufx,',Exy.'//sufx,',Exz.'//sufx
		write(*,'(3a)') 
     1 '                     Eyy.'//sufx,',Eyz.'//sufx,',Ezz.'//sufx

		close(ISTRAIN)
		IOUT=IOUTFILE
		do 777 i=1,7
		close(IOUT)
 		IOUT=IOUT+1
777		continue

	    endif

            if(o_dgrad)  then

		rewind(IDGRAD)

	        IOUT=IOUTFILE
	  	open(unit= IOUT,file='dUxdx.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='dUxdy.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='dUxdz.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='dUydx.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='dUydy.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='dUydz.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='dUzdx.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='dUzdy.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='dUzdz.'//sufx,status='unknown')
		IOUT=IOUT+1

	 	call vect_out(outparam,9,IDGRAD,nx,ny,nz)

		write(*,'(/,a,3(/,3a))') 
     1'  wrote displacement gradient files ',
     2'             dUxdx.'//sufx,',dUxdy.'//sufx,',dUxdz.'//sufx,
     3'             dUydx.'//sufx,',dUydy.'//sufx,',dUydz.'//sufx,
     4'             dUzdx.'//sufx,',dUzdy.'//sufx,',dUzdz.'//sufx

		close(IDGRAD)
		IOUT=IOUTFILE
		do 666 i=1,9
		close(IOUT)
 		IOUT=IOUT+1
666		continue

	    endif

	    if(o_rbr) then
		rewind(IROTAT)
			
		IOUT=IOUTFILE
		open(unit= IOUT,file='rbr_deg_x.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='rbr_deg_y.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='rbr_deg_z.'//sufx,status='unknown')
		IOUT=IOUT+1

	 	call vect_out(outparam,3,IROTAT,nx,ny,nz)

		write(*,'(/,a,3a)') 
     1 '  wrote rigid body rotation files ',
     2 ' rbr_deg_x.'//sufx,',rbr_deg_y.'//sufx,',rbr_deg_z.'//sufx

		close(IROTAT)
		IOUT=IOUTFILE
		do 555 i=1,3
		close(IOUT)
 		IOUT=IOUT+1
555		continue

	    endif

 	    if(o_orient) then

        	rewind(IORIENT)  

                IOUT=IOUTFILE
	        open(unit= IOUT,file='prince_max.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='plunge_max.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='trend_max.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='prince_int.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='plunge_int.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='trend_int.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='prince_min.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='plunge_min.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='trend_min.'//sufx,status='unknown')
				
	 	call vect_out(outparam,9,IORIENT,nx,ny,nz)
			
		write(*,'(/,a,3(/,3a))') '  wrote principal axes files',
     1'   prince_max.'//sufx,',plunge_max.'//sufx,',trend_max.'//sufx,
     2'   prince_int.'//sufx,',plunge_int.'//sufx,',trend_int.'//sufx,
     3'   prince_min.'//sufx,',plunge_min.'//sufx,',trend_min.'//sufx

		close(IORIENT)
		IOUT=IOUTFILE
		do 444 i=1,9
		close(IOUT)
 		IOUT=IOUT+1
444		continue

	    endif

  	    if(o_disp) then

	        rewind(IDISPL)	  
		IOUT=IOUTFILE
		open(unit= IOUT,file='disp1.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='disp2.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='disp3.'//sufx,status='unknown')
				
	 	call vect_out(outparam,3,IDISPL,nx,ny,nz)

		write(*,'(/,3(a,a3))') 
     1   '  wrote displacement files disp1.'//sufx,
     2   ',disp2.'//sufx,',disp3.'//sufx
			
		close(IDISPL)
		IOUT=IOUTFILE
		do 333 i=1,3
		close(IOUT)
 		IOUT=IOUT+1
333		continue

	    endif
	  
 	    if(o_inv) then
		rewind(INVAR)
			
		IOUT=IOUTFILE 		
		open(unit= IOUT,file='delv.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='aftershocks.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='workdil.'//sufx,status='unknown')
		IOUT=IOUT+1
		open(unit= IOUT,file='workdis.'//sufx,status='unknown')

	 	call vect_out(outparam,4,INVAR,nx,ny,nz)

		write(*,'(/,4a)') 
     1'  wrote invariant files delv.'//sufx,',aftershocks.'//sufx,
     1',workdil.'//sufx,',workdis.'//sufx
			
		close(INVAR)
		IOUT=IOUTFILE
		do 222 i=1,4
		close(IOUT)
 		IOUT=IOUT+1
222		continue

	    endif


C---- Re-format selected subsets of the output into vector formatted files.	
C	  Select these subsets interactively. ----  
	  
       else
	  
	    if(o_elem) then
	  
	        do 9 i=1,3
9               outparam(i)=.false.

                write(*,'(/,2a$)')
     &		' Do you want to output the strike-slip ',
     &           'component of element displacement [n]? '
	        read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') outparam(1)=.true.

                write(*,'(2a$)')
     &          ' Do you want to output the dip-slip component of', 
     &          ' element displacement [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') outparam(2)=.true.

                write(*,'(2a$)')
     &          ' Do you want to output the tensile component of',
     &           ' element displacement [n]? '
	        read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') outparam(3)=.true.

		call elem_disps(num_ds,nplane,outparam)
			
	    endif
		
	    if(o_stress) then
	        rewind(ISTRESS)

		do 8 i=1,6
8		outparam(i)=.false.
			
		IOUT=IOUTFILE
	 		
		write(*,'(/,a$)') 
     &	    ' Do you want to output the Sxx stress component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(1)=.true.
		  open(unit= IOUT,file='Sxx.'//sufx,status='unknown')
	 	  IOUT=IOUT+1
	        endif

	        write(*,'(a$)') 
     &	    ' Do you want to output the Sxy stress component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(2)=.true.
		  open(unit= IOUT,file='Sxy.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output the Sxz stress component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(3)=.true.
		  open(unit= IOUT,file='Sxz.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output the Syy stress component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(4)=.true.
		  open(unit= IOUT,file='Syy.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output the Syz stress component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(5)=.true.
		  open(unit= IOUT,file='Syz.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output the Szz stress component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(6)=.true.
		  open(unit= IOUT,file='Szz.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

	 	call vect_out(outparam,6,ISTRESS,nx,ny,nz)

		close(ISTRESS)
		IOUT=IOUTFILE
		do 88 i=1,6
 		if(outparam(i)) then
		  close(IOUT)
 		  IOUT=IOUT+1
                endif
88		continue

	    endif
			
	    if(o_strain) then
		rewind(ISTRAIN)

		do 7 i=1,7
7		outparam(i)=.false.
			
		IOUT=IOUTFILE
	 		
		write(*,'(/,a$)') 
     &	    ' Do you want to output the Exx strain component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(1)=.true.
		  open(unit= IOUT,file='Exx.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output the Exy strain component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
	          outparam(2)=.true.
		  open(unit= IOUT,file='Exy.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output the Exz strain component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(3)=.true.
		  open(unit= IOUT,file='Exz.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output the Eyy strain component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(4)=.true.
		  open(unit= IOUT,file='Eyy.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output the Eyz strain component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(5)=.true.
		  open(unit= IOUT,file='Eyz.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output the Ezz strain component [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(6)=.true.
		  open(unit= IOUT,file='Ezz.'//sufx,status='unknown')
	          IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output the volumetric strain [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(7)=.true.
		  open(unit= IOUT,file='Evol.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

	        call vect_out(outparam,7,ISTRAIN,nx,ny,nz)

	        close(ISTRAIN)
		IOUT=IOUTFILE
		do 77 i=1,7
 		if(outparam(i)) then
		  close(IOUT)
 		  IOUT=IOUT+1
                endif
77		continue

	    endif

	    if(o_dgrad)  then
	        rewind(IDGRAD)

		do 6 i=1,9
6		outparam(i)=.false.
			
		IOUT=IOUTFILE
	 		
		write(*,'(/,a$)') 
     &	    ' Do you want to output displacement gradient dUx/dx [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(1)=.true.
		  open(unit= IOUT,file='dUxdx.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output displacement gradient dUx/dy [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(2)=.true.
		  open(unit= IOUT,file='dUxdy.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

	        write(*,'(a$)') 
     &	    ' Do you want to output displacement gradient dUx/dz [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(3)=.true.
		  open(unit= IOUT,file='dUxdz.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output displacement gradient dUy/dx [n]? '
	        read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(4)=.true.
		  open(unit= IOUT,file='dUydx.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output displacement gradient dUy/dy [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(5)=.true.
		  open(unit= IOUT,file='dUydy.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output displacement gradient dUy/dz [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(6)=.true.
		  open(unit= IOUT,file='dUydz.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output displacement gradient dUz/dx [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
	          outparam(7)=.true.
		  open(unit= IOUT,file='dUzdx.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output displacement gradient dUz/dx [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(8)=.true.
		  open(unit= IOUT,file='dUzdy.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &	    ' Do you want to output displacement gradient dUz/dx [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(9)=.true.
		  open(unit= IOUT,file='dUzdz.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

	 	call vect_out(outparam,9,IDGRAD,nx,ny,nz)

		close(IDGRAD)
		IOUT=IOUTFILE
		do 66 i=1,9
 		if(outparam(i)) then
		  close(IOUT)
 		  IOUT=IOUT+1
                endif
66	        continue

	    endif

	    if(o_rbr) then
                rewind(IROTAT)
			
		do 515 i=1,9
515		outparam(i)=.false.
			
		IOUT=IOUTFILE
                write(*,'(/,2a$)')
     &' Do you want to output rigid body rotations about', 
     &' the in the X (E) axis [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(1)=.true.
		  open(unit= IOUT,file='rbr_deg_x.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

                write(*,'(2a$)')
     &' Do you want to output rigid body rotations about', 
     &' the in the Y (N) axis [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(1)=.true.
	          open(unit= IOUT,file='rbr_deg_y.'//sufx,status='unknown')
	          IOUT=IOUT+1
                endif

                write(*,'(2a$)')
     &' Do you want to output rigid body rotations about', 
     &' the in the Z (up) axis [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(1)=.true.
		  open(unit= IOUT,file='rbr_deg_z.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

	 	call vect_out(outparam,3,IROTAT,nx,ny,nz)

	        close(IROTAT)
		IOUT=IOUTFILE
		do 55 i=1,3
 		if(outparam(i)) then
		  close(IOUT)
 		  IOUT=IOUT+1
                endif
55	        continue

	    endif

	    if(o_orient) then
		rewind(IORIENT)  

		do 4 i=1,9
4		outparam(i)=.false.
			
		IOUT=IOUTFILE
                write(*,'(/,2a$)')
     &' Do you want to the maximum principal strain,its',
     &' plunge and trend [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(1)=.true.
		  outparam(2)=.true.
		  outparam(3)=.true.
		  open(unit= IOUT,file='prince_max.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='plunge_max.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='trend_max.'//sufx,status='unknown')
		  IOUT=IOUT+1
                endif

                write(*,'(2a$)')
     & ' Do you want to the intermediate principal strain,',
     &' its plunge and trend [n]? '
                read(*,'(a1)') ok
                if(ok.eq.'y'.or.ok.eq.'Y') then
	          outparam(4)=.true.
                  outparam(5)=.true.
                  outparam(6)=.true.
		  open(unit= IOUT,file='prince_int.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='plunge_int.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='trend_int.'//sufx,status='unknown')
		  IOUT=IOUT+1
	        endif

                write(*,'(2a$)')
     & ' Do you want to the minimum principal strain,its',
     &' plunge and trend [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(7)=.true.
		  outparam(8)=.true.
		  outparam(9)=.true.
		  open(unit= IOUT,file='prince_min.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='plunge_min.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='trend_min.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

	 	call vect_out(outparam,9,IORIENT,nx,ny,nz)
			
		close(IORIENT)
		IOUT=IOUTFILE
		do 44 i=1,9
 		if(outparam(i)) then
		  close(IOUT)
 		  IOUT=IOUT+1
                endif
44		continue

	    endif

            if(o_disp) then
	        rewind(IDISPL)	  

		do 3 i=1,3
3		outparam(i)=.false.
			
		IOUT=IOUTFILE
	 		
		write(*,'(/,2a$)') 
     &' Do you want to output displacements in ',
     &'the X (E) or strike direction [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(1)=.true.
		  open(unit= IOUT,file='disp1.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(/,2a$)') 
     &' Do you want to output displacements in ',
     &'the Y (N) or dip direction [n]? '
	        read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(2)=.true.
		  open(unit= IOUT,file='disp2.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif
			

		write(*,'(/,2a$)') 
     &' Do you want to output displacements in ',
     &'the Z (up) or normal direction [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(3)=.true.
		  open(unit= IOUT,file='disp3.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

	 	call vect_out(outparam,3,IDISPL,nx,ny,nz)
			
		close(IDISPL)
		IOUT=IOUTFILE
		do 33 i=1,3
 		if(outparam(i)) then
		  close(IOUT)
 		  IOUT=IOUT+1
                endif
33	        continue

	    endif
	  
            if(o_inv) then
		rewind(INVAR)
			
		do 2 i=1,4
2		outparam(i)=.false.
			
		IOUT=IOUTFILE
	 		
		write(*,'(/,a$)') 
     &          ' Do you want to output volumetric strains [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(1)=.true.
		  open(unit= IOUT,file='delv.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)')
     &          ' Do you want to output failure strains [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(2)=.true.
		  open(unit= IOUT,file='aftershocks.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif
			
		write(*,'(a$)') 
     &          ' Do you want to output octahedral strains [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(3)=.true.
		  open(unit= IOUT,file='oss.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(a$)') 
     &          ' Do you want to output strain energy density [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(4)=.true.
		  open(unit= IOUT,file='work.'//sufx,status='unknown')
		endif

	 	call vect_out(outparam,4,INVAR,nx,ny,nz)
			
		close(INVAR)
		IOUT=IOUTFILE
		do 22 i=1,4
 		if(outparam(i)) then
		  close(IOUT)
 		  IOUT=IOUT+1
                endif
22		continue

	    endif

            if(o_fail) then
		rewind(IFAIL)
			
		do 12 i=1,10
12		outparam(i)=.false.
			
		IOUT=IOUTFILE
	 		
		write(*,'(/,2a$)') 
     &  ' Do you want to output strike,dip,slip vector direction ',
     &  ' cosines of 1st failure plane [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(1)=.true.
		  outparam(2)=.true.
		  outparam(3)=.true.
		  outparam(4)=.true.
		  outparam(5)=.true.
		  open(unit= IOUT,file='strike1.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='dip1.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='ce1.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='cn1.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='cz1.'//sufx,status='unknown')
		  IOUT=IOUT+1
		endif

		write(*,'(/,2a$)') 
     &  ' Do you want to output strike,dip,slip vector direction ',
     &  ' cosines of 2nd failure plane [n]? '
		read(*,'(a1)') ok
		if(ok.eq.'y'.or.ok.eq.'Y') then
		  outparam(6)=.true.
		  outparam(7)=.true.
		  outparam(8)=.true.
		  outparam(9)=.true.
		  outparam(10)=.true.
		  open(unit= IOUT,file='strike2.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='dip2.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='ce2.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='cn2.'//sufx,status='unknown')
		  IOUT=IOUT+1
		  open(unit= IOUT,file='cz2.'//sufx,status='unknown')
		endif

	 	call vect_out(outparam,10,IFAIL,nx,ny,nz)
			
		close(IFAIL)
		IOUT=IOUTFILE
		do 99 i=1,10
 		if(outparam(i)) then
		  close(IOUT)
 		  IOUT=IOUT+1
                endif
99		continue

	    endif
	  
	endif		
	return
	end


	subroutine vect_out(outparam,nparam,INPUT,nx,ny,nz)
C********************************************************************
C  Write out the selected parameters from file INPUT in vector format.
C  Each vector (output file record) corresponds to one row in the 
C  x-direction.  For a gridded plane, each row contains the selected
C  output parameter as they vary along strike at a fixed
C  depth (dip index).  The first line for each plane is the top (closest
C  to the surface) of the plane and lines increment with increasing depth.
C********************************************************************
        USE parameters

        real, dimension(:,:), allocatable :: param
	integer :: alloc_stat,INPUT,IOUT
	integer :: nx,ny,nz,nparam
        real :: x,y,z
	logical outparam(*)

        ALLOCATE(param(nparam,nx),STAT=alloc_stat)
        if(alloc_stat.ne.0) then
          write(*,*) ' Error allocating memory in rtn vect_out!'
          stop
        endif
	  
	do 7 iz=1,nz
	  do 8 iy=1,ny
		  
	    do 9 ix=1,nx
	      read(INPUT,'(3(f9.3,1x),1x,10(g11.4,1x))')x,y,z,
     &        (param(j,ix),j=1,nparam)
 9	    continue		
 
            IOUT=IOUTFILE

	    do 6 j=1,nparam
	    if(outparam(j)) then
              write(IOUT,100) (param(j,ix),ix=1,nx)
	      IOUT=IOUT+1
	    endif
 6	    continue
			
 8	  continue		
 7	continue
 100	format(1000(g11.4,1x))

        DEALLOCATE(param)

	return
	end


       subroutine elem_disps(num_ds,nplane,outparam)
C************************************************************************
C  Subroutine to write displacements on each plane, and its subelements,
C  to matrices that are simply interpreted by MATLAB or other typical 
C  plotting routines. 
C   3 file types are created for each plane -
C     elemds#.???   strike component of displacements 
C     elemdd#.???   dip component of displacements 
C     elemdn#.???   tensile component of displacements 
C    # is the number of the plane as listed in the input file.
C	 Each contains a component of the displacement such that each line
C    of the file corresponds to the elements along strike at a fixed
C    depth (dip index).  The first line for each plane is the top (closest
C    to the surface) of the plane and lines increment with increasing depth.
C    A set of lines are written for each plane, each plane is separated
C    by a blank row.
*************************************************************************
*
      USE parameters

      USE allocatable_arrays
	  
      character*12 filename

      logical outparam(*)
	  
      REAL*8 XO,YO,ZO,C,S,DIP,CDIP,SDIP,BWX1,BWX2
      INTEGER*4 NBX1,NBX2
      COMMON/DEFS/XO(MAX_PLN),YO(MAX_PLN),ZO(MAX_PLN),
     &		   C(MAX_PLN),S(MAX_PLN),DIP(MAX_PLN),
     &             CDIP(MAX_PLN),SDIP(MAX_PLN),
     &             BWX1(MAX_PLN),BWX2(MAX_PLN),
     &		   NBX1(MAX_PLN),NBX2(MAX_PLN)
	  
      logical o_elem,o_stress,o_strain,o_dgrad,o_rbr,
     &  o_orient,o_disp,o_inv,xyz_out,o_fail,o_proj
      character*4 sufx
      common/outputs/sufx,o_elem,o_stress,o_strain,o_dgrad,
     &  o_rbr,o_orient,o_disp,o_inv,xyz_out,o_fail,o_proj

      if(outparam(1)) then	 
C     for the strike displacement
         kk=0				
	 k=0
	  
	 filename='elemds01.'//sufx
      	 do 50 np=1,nplane

	   open(unit= IOUTFILE,file=filename,status='unknown')

	   ialong_strike=3*nbx2(np)
C          start at the top of the plane,increment to the bottom	
	   do 49 jd=nbx2(np),1,-1	  				
C            select the strike component of displacement
	     istart=3*jd-2+k	 	 
	     iend=istart+(nbx1(np)-1)*ialong_strike
			  
C **         increment along strike at a fixed depth (dip) index
	     write(IOUTFILE,'(1000(g15.7,1x))') 
     1       (xmatrix(kk,num_ds),kk=istart,iend,ialong_strike) 					
				
 49	   continue

C          increment to the next plane
	   k=k+nbx1(np)*nbx2(np)*3	
				
           write(filename(7:8),'(i2.2)') np+1

	   close(IOUTFILE)
 50      continue
 
       endif
*

       if(outparam(2)) then	  
C        for the dip displacement
	 kk=0				
	 k=0

	 filename='elemdd01.'//sufx
	   
      	 do 51 np=1,nplane

	   open(unit= IOUTFILE,file=filename,status='unknown')

	   ialong_strike=3*nbx2(np)
	   do 52 jd=nbx2(np),1,-1
C            select the dip component of displacement
	     istart=3*jd-1+k	 	  
	     iend=istart+(nbx1(np)-1)*ialong_strike
			   
C **         increment along strike at a fixed depth (dip) index
	     write(IOUTFILE,'(1000(g15.7,1x))') 
     1       (xmatrix(kk,num_ds),kk=istart,iend,ialong_strike) 					
				
 52	   continue

	   k=k+nbx1(np)*nbx2(np)*3

           write(filename(7:8),'(i2.2)') np+1
	   close(IOUTFILE)
51       continue

       endif
*
       if(outparam(3)) then	  

C        for the normal displacement
	 kk=0				
	 k=0
	 filename='elemdn01.'//sufx
	   
      	 do 53 np=1,nplane

	 open(unit= IOUTFILE,file=filename,status='unknown')

	 ialong_strike=3*nbx2(np)
	   do 54 jd=nbx2(np),1,-1
C          select the normal component of displacement
	   istart=3*jd+k	 	  
	   iend=istart+(nbx1(np)-1)*ialong_strike
			   
C **       increment along strike at a fixed depth (dip) index
	   write(IOUTFILE,'(1000(g15.7,1x))') 
     1     (xmatrix(kk,num_ds),kk=istart,iend,ialong_strike) 					
				
 54	   continue

	 k=k+nbx1(np)*nbx2(np)*3

         write(filename(7:8),'(i2.2)') np+1
	 close(IOUTFILE)
 53      continue
       endif

       return
       end
	  
