	  SUBROUTINE GRID(NUM_Ds,NPLANE)
C******************************************************************************
C   This code uses 3-d dislocation solutions for the stress tensor, displacement
C   gradients and the displacement vectors as a function of position. 
C
C	The input file is assumed to be still open.  If the program is being run
C	in a non-interactive mode, the user specifies all the options describing
C	the grid parameters, etc. to the bottom of the input file.  If it is run 
C	interactively, then the user must enter grid parameters interactively 
C   (there are lots of questions to answer!).
C
C   All routines that control output to X,Y,Z,Data formatted files are 
C   contained in this file.  These output files are read by routines in
C	the file OUTPUT_VECTOR.F and re-formatted in a vector format (e.g.
C   as required by plotting programs such as MATLAB); this re-formatting
C   is of course optional.
C******************************************************************************
          USE parameters

	  INTEGER*4 NUM_Ds,NPLANE
 
          character*1 answer
	  character*1 flag
	  	  
	  read(INFILE,'(a1)') flag
	  if(flag.eq.'i'.or.flag.eq.'I')then
C         Interactive option chosen
            write(*,'(/,a)') 
     &' You can examine a volume,plane or user supplied points.'
            write(*,'(a)') 
     &'  Examine a volume [v],plane [p],or user points [u] '
            write(*,'(a$)') 
     &      '  (default is a plane - just hit return)? '
            read(*,'(a1)') answer

            if(answer.eq.'v'.or.answer.eq.'V')then
	        call volume(flag,NPLANE,NUM_Ds)
            elseif(answer.eq.'u'.or.answer.eq.'U')then
	        call userco(flag,NPLANE,NUM_Ds)
            else
	      call plane(flag,NPLANE,NUM_Ds)
	    endif
		
	else
	    read(INFILE,'(a1)') answer
	    if(answer.eq.'v'.or.answer.eq.'V')then
	      call volume(flag,NPLANE,NUM_Ds)
            elseif(answer.eq.'u'.or.answer.eq.'U')then
	        call userco(flag,NPLANE,NUM_Ds)
	    else
	      call plane(flag,NPLANE,NUM_Ds)
	    endif
	endif
 
	return
	end

	  subroutine volume(flag,NPLANE,NUM_Ds) 
C******************************************************************************
C	Subroutine to enter the coordinates of a volume with horizontal and
C	vertical sides and the increment size of grid units. x,y, and z are 
C	transient coordinates.  nx,ny, and nz
C	are the number of increments along the sides of the volume.
C
C   Note that all output parameters, except element relative displacements,
C   will be in global coordinates.
C******************************************************************************
          USE parameters
	  
	  character*1 flag,dummy

	  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
	  
	  real*8 x1,y1,z1,x2,y2,z2,x,y,z,xinc,yinc,zinc
	  real*8 xlong,ywide,zdeep
	  
C         initialize output flags - assume no output
	  call init_opts		

C---- Read input from the input file. -------	  
C         Read parameters from INPUT file
	  if(flag.ne.'i'.and.flag.ne.'I')then	

C               Determine what information to write to output files
		call get_fopts				
		
C           only write out element displacements
	    if(o_elem.and..not.o_stress.and..not.o_strain.and.
     &	.not.o_dgrad.and..not.o_rbr.and..not.o_orient.and.
     &	.not.o_disp.and..not.o_inv.and..not.o_fail) go to 33		

C               Upper NW and lower SE corners of the volume
		read(infile,*) x1,y1,z1,x2,y2,z2	
C               Descriptor line
	  	read(infile,'(a1)') dummy			
C               Increments in X (E), Y (N),Z (up) directions
		read(infile,*) xinc,yinc,zinc		
 
 9		if(z1.gt.0.0)then
		  write(*,*) ' z must be 0.0 or less. Try again.'
		  write(*,*)' Enter the z value, or to quit enter 999'
		  read(*,*)z1
			if(z1.eq.999.) stop
		  go to 9	
		endif
 7	        if(z2.gt.z1) then
 
	          write(*,*)
     &          ' The second z value must be less than the first.'
	  	  write(*,*) ' Enter a new z or z=999 to quit'
		  read(*,*) z2
		  if(z2.eq.999.) return
		  go to 7
	        endif

C---- Read input from the terminal. -------	    

	  else	
	  
C               ask the user what kind of output files to write
	  	call get_iopts		
		
C           only write out element displacements
	    if(o_elem.and..not.o_stress.and..not.o_strain.and.
     &	.not.o_dgrad.and..not.o_rbr.and..not.o_orient.and.
     &	.not.o_disp.and..not.o_inv.and..not.o_fail) go to 33		

	    write(*,'(/,a)') ' You will need to enter two sets of coordinates.'
        write(*,'(2a)') 
     >  ' Let the positive coordinate axes, x,y,and z, be east, north',
     >  ' and up, respecively.'
	    write(*,'(2a)') ' The first set of coordinates should be ',
     >  ' at the northwestern corner of the top surface'
	    write(*,'(2a)') ' the second at ',
     >  'the southeastern corner of the lower surface.'
	    write(*,'(2a)') ' The volume ',
     >  'must have horizontal and vertical sides.'

            write(*,*)
	    write(*,*)'                 1                    z'
	    write(*,*)'               *   *                  z'
	    write(*,*)'            *         *               z            (N)'
	    write(*,*)'         *               *            z          y'
	    write(*,*)'      *                     *         z       y'
	    write(*,*)'      *  *                     *      z    y'
	    write(*,*)'      *     *                     *   z y'
	    write(*,*)'      *        *                      *'
	    write(*,*)'      *           *                 * *   x'
	    write(*,*)'         *           *           *    *       x'
	    write(*,*)'            *           *     *       *           x'
	    write(*,*)'               *           *          *             (E)'
	    write(*,*)'                  *        *        *'
	    write(*,*)'                     *     *     *'
	    write(*,*)'                        *  *  *'
	    write(*,*)'                           2'
	    write(*,*)
	    write(*,*)
	    write(*,'(a$)') ' Enter the upper NW corner <x,y,z>: '
	    read(*,*) x1,y1,z1
 
 4	    if(z1.gt.0.0) then
	      write(*,*)' z must be 0.0 or less. Try again.'
	      write(*,*)' Enter the z value, or to quit enter 999'
	      read(*,*) z1
	      if(z1.eq.999.) return
	      go to 4
	    endif
 
	    write(*,'(a$)') ' Enter the lower SE corner <x,y,z>: '
	    read(*,*) x2,y2,z2
 5	    if(z2.gt.z1) then
	      write(*,*) ' The second z value must be less than the first.'
	      write(*,*) ' Enter a new z or z=999 to quit'
	      read(*,*) z2
	      if(z2.eq.999.) return
	      go to 5
	    endif
 
 3	    write(*,'(a$)')' Enter the x (E/W) increment length: '
      	    read(*,*) xinc
	    if(xinc.le.0.0.or.xinc.gt.xlong) then
	      write(*,*)'Invalid increment - try again'
	      go to 3
	    endif
 6	    write(*,'(a$)')' Enter the y (N/S) increment length: '
      	    read(*,*) yinc
	    if(yinc.le.0.0.or.yinc.gt.ywide) then
	      write(*,*)'Invalid increment - try again'
	      go to 6
	    endif
 8	    write(*,'(a$)')' Enter the vertical increment length: '
      	    read(*,*) zinc
	    if(zinc.le.0.0.or.zinc.gt.zdeep) then
	      write(*,*)'Invalid increment - try again'
	      go to 8
	    endif
				
	  endif
        
C---  Be sure that no projection is requested
          o_proj=.false.
	  

C---- Reset the number of increments and initialize origin. -------	    
	  xlong=x2-x1
	  ywide=y1-y2
	  zdeep=z1-z2

	  nz=int(zdeep/zinc)+1			
	  if((nz*zinc).lt.zdeep) nz=nz+1
          z=z1						
 										
	  nx=int(xlong/xinc)+1
	  if((nx*xinc).lt.xlong) nx=nx+1
	  ny=int(ywide/yinc)+1
	  if((ny*yinc).lt.ywide) ny=ny+1
	  
	  x=x1+.0001				
	  y=y1		  
 
	  x0=x
	  y0=y
	  z0=z

C---- Open appropriate files. -------	      
	  call open_opts


C---- Do calculations on the grid and write the results. -------
      ntotal=nz*ny*nx
	  do 44 iz=1,nz
	  y=y0
		do 55 iy=1,ny
		x=x0
			do 66 ix=1,nx
  	  		call XYZ_Results(x,y,z,NPLANE,NUM_Ds)
			
	  		x=x+xinc

 66		  continue
		  y=y-yinc
			
 55		continue
		z=z-zinc

		percent=float(iz*ny*nx*100/ntotal)
		write(*,*)'working.....',percent,' percent done'
		
 44	  continue


C---- Write the element displacements and Vector output if requested  	      	  
C     and close all files. -------	      	  
 33	  call final_write(flag,NPLANE,NUM_Ds,nz,ny,nx)
 
 	  return
	  end

	  subroutine plane(flag,NPLANE,NUM_Ds)
C****************************************************************************
C  Subroutine to read parameters for output on a gridded plane.
C
C  Note that output parameters that are referenced to the global 
C  coordinate system include:
C		rigid body rotations
C		orientation of principal strain axes
C               failure plane descriptions
C  output parameters that are in the coordinate system of the gridded plane 
C  (X=strike,Y=dip,Z=normal directions) include:
C  		stress tensor elements
C  		strain tensor elements
C  		displacement gradients tensor elements
C		displacements
C  Element relative displacements correspond to the strike,dip,normal directions 
C  of each element's plane.
C***************************************************************************

          USE parameters
	  character*1 flag,dummy
	  
	  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

	  real*8 x,y,z,xo,yo,zo,stk,dip,xlength,width

	  REAL*8 SSTK,CSTK,SDIP,CDIP
	  COMMON/TOPLN/SSTK,CSTK,SDIP,CDIP

          real*8 stress,ewe,dgrten,pr,ym,friction,tmax,tmaxo,tmaxa
	  common/mike/stress(6),ewe(3),dgrten(3,3),pr,ym,friction,
     &                tmax,tmaxo(2),tmaxa
	  
	  real*8 e
	  common/temps/e(3,3)

C         scaling from degrees to radians
	  DATA TORAD/.017453293/		

C         initialize output flags - assume no output
	  call init_opts				


C---- Read input from the input file. -------	  
	
C         Read parameters from INPUT
	  if(flag.ne.'i'.and.flag.ne.'I')then		

C               Determine what information to write to output files
		call get_fopts			

C           only write out element displacements
	    if(o_elem.and..not.o_stress.and..not.o_strain.and.
     &	.not.o_dgrad.and..not.o_rbr.and..not.o_orient.and.
     &	.not.o_disp.and..not.o_inv.and..not.o_fail) go to 33		

		read(infile,*) xo,yo,zo,stk,dip,xlength,width
	        read(infile,'(a1)') dummy		
		read(infile,*) nl,nw
		xo=xo+.0001

 9		if(zo.gt.0.0)then
		   write(*,*) ' z must be 0.0 or less. Try again.'
		   write(*,'(a$)') ' Enter the z value, or to quit enter 999: '
		   read(*,*) zo
		   if(zo.eq.999.) return
		   go to 9
		endif

		

C---- Read input from the terminal. -------	    
	  else				

C               ask the user what kind of output files to write
		call get_iopts			

	    if(o_elem.and..not.o_stress.and..not.o_strain.and.
     &	.not.o_dgrad.and..not.o_rbr.and..not.o_orient.and.
     &	.not.o_disp.and..not.o_inv.and..not.o_fail) go to 33		

	 write(*,'(/,a)') ' Enter reference corner <Xo,Yo,Zo>, '
      	write(*,'(a,/,a)')  '       strike (cw from N) & dip (0-90), ',
     >  '       length <L> along strike, width along dip <W>.'
            write(*,*) ' '
	    write(*,*)'                          +z'
	    write(*,*)'                           z'
	    write(*,*)'                  x        z        (N)'
	    write(*,*)'                     x     z      y' 
	    write(*,*)'                        x  z    y'
	    write(*,*)'                           x y'            
	    write(*,*)'                          y   x'
	    write(*,*)'                       y   z |   x'           
	    write(*,*)'                    y      z d     (E)'    
	    write(*,*)'                           z e '           
	    write(*,*)'                           z p '                              
	    write(*,*)'  Xo,Yo,Zo                 z t'
	    write(*,*)'      * * * * * * * */     z h '
	    write(*,*)'     *             */.     z | '
	    write(*,*)'    *             */.      z | '
	    write(*,*)'   *             * width      '
	    write(*,*)'  *             */.            '
	    write(*,*)' * * * * * * * */ '
	    write(*,*)' |---length----|'
	    write(*,*)

      	    write(*,*)
     &      ' Enter <Xo,Yo,Zo> (all z''s must be 0 or negative):'
      	    read(*,*) xo,yo,zo
	    xo=xo+.0001
 
 4	    if(zo.gt.0.0)then
	        write(*,*) ' z must be 0.0 or less. Try again.'
		write(*,'(a$)') ' Enter the z value, or to quit enter 999'
		read(*,*) zo
		if(zo.eq.999.) return
		go to 4
	    endif
  
	    write(*,'(a$)') ' Enter strike,dip,L,W: '
	    read(*,*) stk,dip,xlength,width
 
	    write(*,'(a$)') ' Enter number of increments along L, along W: '
	    read(*,*) nl,nw
  		
      endif


C---- Reset the number of increments and define grid spacing. -------	    
	  
	  xinc=xlength/float(nl)
	  yinc=width/float(nw)
	  nl=nl+1
	  nw=nw+1
 
	  rdip=dip*torad
	  rstk=stk*torad
	  sstk=sin(rstk)
	  cstk=cos(rstk)
	  sdip=sin(rdip)
	  cdip=cos(rdip)

C	zinc and yinc are redefined to be the increment lengths projected 
C	to the horizontal surface.
	  zinc=yinc*sdip
	  yinc=yinc*cdip


C---- Open appropriate files and do calculations on the grid. -------	      
 	  call open_opts

C---- Do calculations on the grid and write the results. -------	      	  
	  do 2 j=1,nw
		x = xo+float(j-1)*cstk*yinc
		y = yo-float(j-1)*sstk*yinc
		z = zo-float(j-1)*zinc
                percent=float(j)*100./float(nw)
		write(*,*)'working.....',percent,' percent done'
	   	do 2 i=1,nl

  	  		call XYZ_Results(x,y,z,NPLANE,NUM_Ds)
 
			x = x+sstk*xinc
			y = y+cstk*xinc
 
 2	  continue
 
 
C---- Write the element displacements and Vector output if requested  	      	  
C     and close all files. -------	      	  
33    call final_write(flag,NPLANE,NUM_Ds,1,nw,nl)

      return
      end
	  
	  
	  subroutine final_write(flag,NPLANE,NUM_Ds,nz,ny,nx)
C****************************************************************************
C  Write out element displacements in X,Y,Z,Data format if desired.  Read
C  X,Y,Z,Data formatted files are write new ones in vector format if
C  desired.
C
C  Close all files.
C****************************************************************************
          USE parameters
	  
	  character*1 flag
	  
	  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(xyz_out) then
C------- XYZ output	-------	  
	    if(o_elem) then
	      open(unit= IELEM,file='elements.'//sufx,status='unknown')
	      call XYZ_elems(NPLANE,NUM_Ds)
	      close(IELEM)
	    endif
	    call close_opts
		
	  else
C------- Vector (MATLAB) output	-------
	    call vector(flag,NPLANE,NUM_Ds,nz,ny,nx)
	  
	  endif

	  return
	  end


        subroutine get_fopts
C****************************************************************************
C This routine reads the input file options and sets output file options
C accordingly.  The input file should look like
C       * local coordinate output (global otherwise) - only for planar grids
C       global
C  	* Output information (each line results in a separate output file)
C	stress
C	strain
C	gradients of displacement
C	rigid body rotations
C	orientation of principal strain axes
C	displacements
C       failure plane 
C	invarients of the strain field
C	element relative displacements
C       * Output file suffix
C       dat1
C Remember that lines beginning with * are only descriptions and aren't used
C by the program (they can say anything).  The list of files containing various
C characteristics of the calculated deformation field can be in any order. 
C The word on each line sets an appropriate flag so that an output file 
C containing the parameters described by the word will be generated.  
C The program only interprets the first letter (the first 4 letters for
C 'stress' and 'strain') so that their is some flexibility in what words
C are used.  For example, you must write 'gradients of displacement' or
C simply 'gradients' to generate an output file containing the components
C of the displacement gradient tensor; use of 'displacement gradients'
C would be incorrect (the key word must start with a 'g') and result in
C an output file containing displacements.
C****************************************************************************
        USE parameters
	
	character*4 o_temp
	  
	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
		
C Get output file type [XYZ or Vector]
        read(infile,'(a)') o_temp
	if(o_temp(1:1).ne.'x'.and.o_temp(1:1).ne.'X') xyz_out=.false.
	if(xyz_out) then
	   write(*,'(/,a,/)') ' Writing only xyz output files.'
        else
	   write(*,'(/,a,/)') ' Writing xyz and vector output files.'
        endif
 
C       descriptor line
	read(infile,'(a)') o_temp		
C       projection type if planar output
	read(infile,'(a)') o_temp		
        if(o_temp(1:1).ne.'g'.and.o_temp(1:1).ne.'G') o_proj=.true.

C       descriptor line
	read(infile,'(a)') o_temp		
		
C Get output file contents 
12	read(infile,'(a)') o_temp
	if(o_temp(1:1).ne.'*') then
	  if(o_temp(1:4).eq.'stre'.or.o_temp(1:4).eq.'STRE') then
	    o_stress=.true.
	    write(*,*) ' writing stress tensor'
	  elseif(o_temp(1:4).eq.'stra'.or.o_temp(1:4).eq.'STRA') then
	    o_strain=.true.
	    write(*,*) ' writing strain tensor'
          elseif(o_temp(1:1).eq.'g'.or.o_temp(1:1).eq.'G') then
	    o_dgrad=.true.
	    write(*,*)
     &         ' writing displacement gradient tensor'
	  elseif(o_temp(1:1).eq.'r'.or.o_temp(1:1).eq.'R') then
	    o_rbr=.true.
	    write(*,*) ' writing rigid body rotations'
	  elseif(o_temp(1:1).eq.'o'.or.o_temp(1:1).eq.'O') then
	    o_orient=.true.
	    write(*,*) ' writing principal axes'
	  elseif(o_temp(1:1).eq.'d'.or.o_temp(1:1).eq.'D') then
	    o_disp=.true.
	    write(*,*) ' writing displacements'
	  elseif(o_temp(1:1).eq.'i'.or.o_temp(1:1).eq.'I') then
	    o_inv=.true.
	    write(*,*) ' writing stress/strain invariants'
	  elseif(o_temp(1:1).eq.'f'.or.o_temp(1:1).eq.'F') then
	    o_fail=.true.
	    write(*,*) ' writing failure planes'
	  elseif(o_temp(1:1).eq.'e'.or.o_temp(1:1).eq.'E') then
	    o_elem=.true.
	    write(*,*) ' writing element displacements'
	  endif
	  goto 12
	endif		

       read(infile,'(a4)') sufx
       read(infile,'(a)') o_temp
	   
        write(*,'(a)') ' '
        return
	end
	  
	subroutine get_iopts
C****************************************************************************
C This polls the user and sets output file options accordingly.  The user
C may just hit return to answer the questions if the default answer (in [])
C is okay. 
C****************************************************************************
		
	character*1 o_temp
	  
	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
		
*       Get output file type 
        write(*,'(/,2a)') 
     &   ' Do you want an X,Y,Z,DATA or Vector output ',
     &   ' format [xyz or vector]'
	write(*,'(a$)') 
     &   ' (xyz is the default - just hit return)? '
 	read(*,'(a1)') o_temp
	if(o_temp(1:1).eq.'v'.or.o_temp(1:1).eq.'V') xyz_out=.false.

*       Get output file suffix 
        write(*,'(/,a$)') 
     &   ' Enter output file suffix: '
        read(*,'(a4)') sufx

*	The next few lines of the input file indicate the type of info
*	wanted and written to files.
	
	write(*,'(/,a)') 
     &  ' Here are all your options -  Answer yes or no.'
	write(*,'(a)') 
     &  ' (y or Y for yes.  Default is always no - just hit return.)'
	
        write(*,'(/,a)')
     & ' Project output to local coord system of your gridded plane '
        write(*,'(a$)')
     & '  (this is meaningless for output on a gridded volume) [n]? '
        read(*,'(a1)') o_temp
	if(o_temp(1:1).eq.'y'.or.o_temp(1:1).eq.'Y')  o_proj=.true.

        write(*,'(/,a$)')
     & ' Do you want to record element displacements (relative) [n]? '
        read(*,'(a1)') o_temp
	if(o_temp(1:1).eq.'y'.or.o_temp(1:1).eq.'Y')  o_elem=.true.

	write(*,'(a$)')' Do you want to record stress components [n]? '
	read(*,'(a1)') o_temp
	if(o_temp(1:1).eq.'y'.or.o_temp(1:1).eq.'Y')  o_stress=.true.

	write(*,'(a$)')' Do you want to record strain components [n] ?'
	read(*,'(a1)') o_temp
	if(o_temp(1:1).eq.'y'.or.o_temp(1:1).eq.'Y')  o_strain=.true.

	write(*,'(a$)')
     &  ' Do you want to record displacement components [n]? '
        read(*,'(a1)') o_temp
	if(o_temp(1:1).eq.'y'.or.o_temp(1:1).eq.'Y')  o_disp=.true.

	write(*,'(a$)')
     & ' Do you want to record displacement gradient components [n]? '
        read(*,'(a1)') o_temp
	if(o_temp(1:1).eq.'y'.or.o_temp(1:1).eq.'Y')  o_dgrad=.true.

	write(*,'(a$)')
     &  ' Do you want to record rigid body rotation components [n]? '
        read(*,'(a1)') o_temp
	if(o_temp(1:1).eq.'y'.or.o_temp(1:1).eq.'Y')  o_rbr=.true.

	write(*,'(a$)')
     &  ' Do you want orientations of principal strains [n]? '
	read(*,'(a1)') o_temp
	if(o_temp(1:1).eq.'y'.or.o_temp(1:1).eq.'Y') o_orient=.true.

	write(*,'(a)')' Do you want stress & strain invariants [n]'
	write(*,'(a,a$)')
     &  '  (volumetric strain,failure & octahedral stresses,',
     1   'strain energy density)? '
        read(*,'(a1)') o_temp	
	if(o_temp(1:1).eq.'y'.or.o_temp(1:1).eq.'Y')  o_inv=.true.
	  
	write(*,'(a)')' Do you want failure plane descriptions '
	write(*,'(a,a$)')
     &  ' (strike,dip,slip vector direction cosines',
     1   ' for 2 planes) [n]? '
        read(*,'(a1)') o_temp	
	if(o_temp(1:1).eq.'y'.or.o_temp(1:1).eq.'Y') o_fail=.true.
	return
	end

	  	  
	  subroutine init_opts
C****************************************************************************
C  Set all output flags to .false. except that specifying generation of
C  X,Y,Z,Data formatted files and output of element displacements.
C****************************************************************************
  
	  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

          o_proj=.false.
	  o_stress=.false.
	  o_strain=.false.
	  o_dgrad=.false.
	  o_fail=.false.
	  o_rbr=.false.
	  o_orient=.false.
	  o_disp=.false.
	  o_inv=.false.
	  o_elem=.false.
	  xyz_out=.true.
	  
	  return
	  end

	  subroutine open_opts
C****************************************************************************
C  Open all flagged X,Y,Z,Data formatted output files.
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
	  
	  if(o_stress) open(unit= ISTRESS,file='stressten.'//sufx,
     &      status='unknown')
	  if(o_fail) open(unit= IFAIL,file='failure.'//sufx,status='unknown')
	  if(o_strain) open(unit= ISTRAIN,file='strainten.'//sufx,
     &      status='unknown')
	  if(o_dgrad)  open(unit= IDGRAD,file='dgten.'//sufx,
     &      status='unknown')
	  if(o_rbr) open(unit= IROTAT,file='rbr.'//sufx,
     &      status='unknown')
	  if(o_orient) 
     &    open(unit= IORIENT,file='strainorient.'//sufx,  
     &      status='unknown')
	  if(o_disp) 
     &    open(unit= IDISPL,file='displacements.'//sufx,	  
     &      status='unknown')
	  if(o_inv) open(unit= INVAR,file='invariants.'//sufx,	  
     &      status='unknown')

	  return
	  end

          subroutine close_opts
C****************************************************************************
C  Close all flagged X,Y,Z,Data formatted output files.
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
		
	  	if(o_stress) close(ISTRESS)
	  	if(o_fail) close(IFAIL)
	  	if(o_strain) close(ISTRAIN)
	  	if(o_dgrad)  close(IDGRAD)
	  	if(o_rbr) close(IROTAT)
	  	if(o_orient) close(IORIENT)  
	  	if(o_disp) close(IDISPL)	  
	  	if(o_inv) close(INVAR)	  

	  return
	  end
	  
	  subroutine XYZ_Results(x,y,z,NPLANE,NUM_Ds)
C******************************************************************************
C  Subroutine to calculate various parameters describing the deformation field.
C  These are output to ascii files. 
C
C  In this routine seven files are written to; at each point X,Y,Z in global 
C  coordinates the following is output (the coordinates of each point are also
C  written to the file):
C    File 		         Contents
C    stressten.??? - Six independent elements of the stress tensor, 
C		     Sxx,Sxy,Sxz,Syy,Syz,Szz. If output on a plane is 
C                    requested these correspond to the in-plane coordinates 
C		     (e.g. Szz is the stress normal to the fault plane,
C                    Sxy is the traction in the strike direction,
C                        is the traction in the up-dip direction, etc.)
C    strainten.??? - Six independent elements of the strain tensor, Exx,
C		     Exx,Exy,Exz,Eyy,Eyz,Ezz, and the volumetric strain 
C		     (Exx+Eyy+Ezz,independent of the coordinate syst
C		     If output on a plane is requested these correspond to  
C                    the in-plane coordinates as in stressten.out.
C   displacements.??? - Displacements in the x,y,z directions. If output 
C	             on a plane is requested these correspond to thespond to the 
C		     in-plane coordinates; Ux is in the strike direction,directio
C                    Uy in the up-dip direction, Uz is the normal direction.l to the plane.
C    dgten.??? -     Displacement gradient tensor; dUx/dx,dUx/dy,dUx/dz,
C		     dUy/dx,dUy/dy,dUy/dz, dUz/dx,dUz/dy,dUz/dz.  As in,
C                    the above, these correspond to the in-plane coordinates
C                    if the user has selected the plane option.
C    strainorient.??? -	 Principal strains,their plunge and trend with respect
C		     to the global axes.  They are listed from max. to min.from
C    invariants.??? -Strain & stress invarients: the volume change, failure 
C	             stresses,octahedral shear,strain energy density.
C    rbr.??? -	     Rotations about the global x,y,and z axes
C    failure.??? -   Failure plane descriptions: strike,dip,slip vector
C                    direction cosines for each plane.
C  Remember that if you change the format of these files, you must also
C  change the routines in OUTPUT_VECTOR.f since they must read these files.
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

	  real*8  x,y,z,delV
	  REAL*8 zd
	  
	  real*8 rbrots(3)
	  
          REAL*8 E
          COMMON/TEMPS/e(3,3)

	  REAL*8 STRESS,EWE,dgrten,pr,ym,friction,tmax,tmaxo,tmaxa
	  COMMON/MIKE/STRESS(6),EWE(3),dgrten(3,3),
     &                pr,ym,friction,tmax,tmaxo(2),tmaxa

	  REAL*8 SSTRIKE,CSTRIKE,SDYP,CDYP
	  COMMON/TOPLN/SSTRIKE,CSTRIKE,SDYP,CDYP

	  real*8 ex,ey,ez,px,tx,py,ty,pz,tz
	  real*8  volchg,critic,workdil,workdis
	  real*8  sp1,d1,rake1,sp2,d2,rake2
	  
	  DATA TODEG/57.295779/		

	  zd=z	  
	  
C         Gets stress and displacements
	  call CALCPO(x,y,zd,NPLANE,NUM_Ds,ierr)	
C         if(ierr.ne.0) return
 
C         Gets strain tensor and vol. change
	  call strain(e,delV)		

C Gets principal strains and orientations and write to output file
	  if(o_orient) then
		call third_results(ex,ey,ez,px,tx,py,ty,pz,tz) 
		write(IORIENT,'(3(f9.3,1x),1x,9(g11.4,1x))') x,y,z,
     &      ex,px,tx,ey,py,ty,ez,pz,tz
	  endif
 	
C rbrots is a column vector (3,1) that contains the rotations
C	about the x-axis (1),y-axis (2), and the z-axis (3)
C	(See Ramsay, 1967, p.170.)    
	  if(o_rbr)then
	 	 rbrots(1)=(dgrten(3,2)-dgrten(2,3))*.5
	  	 rbrots(1)=atan(rbrots(1))*TODEG
	  	 rbrots(2)=(dgrten(1,3)-dgrten(3,1))*.5
	  	 rbrots(2)=atan(rbrots(2))*TODEG
	  	 rbrots(3)=(dgrten(2,1)-dgrten(1,2))*.5
	  	 rbrots(3)=atan(rbrots(3))*TODEG
 	  	 write(IROTAT,'(3(f9.3,1x),1x,3(g11.4,1x))') x,y,z,
     &       rbrots(1),rbrots(2),rbrots(3)
	  endif

C     Rotate to in-plane coordinates 
          if(o_proj) CALL TOPLANE(SSTRIKE,CSTRIKE,SDYP,CDYP)
C
C  Code added 2/15/94 to get max shear stress and its orientation on the inspection 
C  plane.  See subroutine TOPLANE in 3dmain.f and manual for more info.
C
	tmaxa=tmaxa*TODEG
	  
C----- Now write to external files if so deemed -----
	  if(o_stress) 
     & write(ISTRESS,'(3(f9.3,1x),1x,10(g11.4,1x))') x,y,z,
     &   STRESS(1),STRESS(2),STRESS(3),
     &   STRESS(4),STRESS(5),STRESS(6),
     &   tmax,tmaxa,tmaxo(1),tmaxo(2)

	  if(o_strain) 
     & write(ISTRAIN,'(3(f9.3,1x),1x,7(g11.4,1x))') x,y,z,
     &   e(1,1),e(1,2),e(1,3),e(2,2),e(2,3),e(3,3),delV

	  if(o_disp) 
     & write(IDISPL,'(3(f9.3,1x),1x,3(g11.4,1x))') x,y,z,
     &   EWE(1),EWE(2),EWE(3)
	 
	  if(o_dgrad) 
     & write(IDGRAD,'(3(f9.3,1x),1x,9(g11.4,1x))') x,y,z,
     &   dgrten(1,1),dgrten(1,2),dgrten(1,3),
     &   dgrten(2,1),dgrten(2,2),dgrten(2,3),
     &   dgrten(3,1),dgrten(3,2),dgrten(3,3)

	  if(o_inv) then
	    call invariants(volchg,critic,workdil,workdis)
 	    write(INVAR,'(3(f9.3,1x),1x,4(g11.4,1x))') x,y,z,
     &      volchg,critic,workdil,workdis
	  endif
 	
	  if(o_fail) then
	    call failure_planes(sp1,d1,rake1,sp2,d2,rake2)
 	    write(IFAIL,'(3(f9.3,1x),1x,6(g11.4,1x))') x,y,z,
     &      sp1,d1,rake1,sp2,d2,rake2
	  endif

          return
	  end

	  subroutine XYZ_elems(NPLANE,NUM_Ds)
C********************************************************************
C Subroutine to write out element displacements (relative) in
C X,Y,Z,Data format.  Each line of the output file will contain the
C coordinates of the center of the element in global and in-plane
C coordinate systems,and the relative
C displacements in the strike, dip, and tensile (normal)
C directions.
C********************************************************************

      USE parameters

      USE allocatable_arrays

      REAL*8 UG2P,SG2P
      COMMON/TRANS/UG2P(3,3,MAX_PLN),SG2P(3,6,MAX_PLN)

      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)
			
      IBC=0
      DO 45 NP=1,NPLANE		  
C----- sub-element half-widths -----
      BX1P=BWX1(NP)*.5		
      BX2P=BWX2(NP)*.5

      DO 45 JP=1,NBX1(NP)
      DO 45 IP=1,NBX2(NP)

C----- Find center of this sub-element in global coordinates -----    
C      - center of sub-element IP,JP in local NP coords	 
C                center in local strike direction
	 	 XNP=BX1P*(2*JP-1)				
	 	 ZNP=BX2P*(2*IP-1)*SDIP(NP)
C                center in local dip direction
	 	 YNP=BX2P*(2*IP-1)*CDIP(NP)		
      	 
C      - center of sub-element IP,JP in global coords (Zo shifted)
	 	 Z=ZNP
	 	 X=S(NP)*XNP-C(NP)*YNP
	 	 Y=C(NP)*XNP+S(NP)*YNP

C      - center of sub-element IP,JP in in-plane coords
                 XNP=UG2P(1,1,NP)*X+UG2P(1,2,NP)*Y+UG2P(1,3,NP)*Z
                 YNP=UG2P(2,1,NP)*X+UG2P(2,2,NP)*Y+UG2P(2,3,NP)*Z
                 ZNP=UG2P(3,1,NP)*X+UG2P(3,2,NP)*Y+UG2P(3,3,NP)*Z

	 	 X=X+XO(NP)
	 	 Y=Y+YO(NP)
	 	 Z=-ZO(NP)+Z
		 
C----- Write out relative displacements -----    
		write(IELEM,'(2(3(f9.3,1x),1x),3(g15.4,1x))') x,y,z,
     &      XNP,YNP,ZNP,(XMATRIX(IBC+J,NUM_Ds),J=1,3)

 45	  IBC=IBC+3
	  
	  return
	  end

	  subroutine userco(flag,NPLANE,NUM_Ds)
C****************************************************************************
C  Subroutine to read parameters for output at user coordinates
c        found on the input file
C
C  Note that output parameters that are referenced to the global 
C  coordinate system include:
C		rigid body rotations
C		orientation of principle strain axes
C               failure plane descriptions
C  output parameters that are in the coordinate system of the gridded plane 
C  (X=strike,Y=dip,Z=normal directions) include:
C  		stress tensor elements
C  		strain tensor elements
C  		displacement gradients tensor elements
C		displacements
C  Element relative displacements correspond to the strike,dip,normal directions 
C  of each element's plane.
c c.meertens 24feb93
C***************************************************************************
          USE parameters

	  character*1 flag
	  
	  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

	  real*8 x,y,z

          real*8 stress,ewe,dgrten,pr,ym,friction,tmax,tmaxo,tmaxa
	  common/mike/stress(6),ewe(3),dgrten(3,3),pr,ym,friction,
     &                tmax,tmaxo(2),tmaxa
	  
	  real*8 e
	  common/temps/e(3,3)
c user coordinates
        real*8, DIMENSION(:),ALLOCATABLE :: xu,yu,zu
       
        integer :: iuco,nuco,alloc_stat

C         scaling from degrees to radians
	  DATA TORAD/.017453293/		

C         initialize output flags - assume no output
	  call init_opts				

C---- Read input from the input file. -------

C         Read parameters from INPUT
          if(flag.ne.'i'.and.flag.ne.'I')then

C               Determine what information to write to output files
                call get_fopts

C           only write out element displacements
            if(o_elem.and..not.o_stress.and..not.o_strain.and.
     &  .not.o_dgrad.and..not.o_rbr.and..not.o_orient.and.
     &  .not.o_disp.and..not.o_inv.and..not.o_fail) go to 33
C---- Read input from the terminal. -------
          else
C               ask the user what kind of output files to write
                call get_iopts

            if(o_elem.and..not.o_stress.and..not.o_strain.and.
     &  .not.o_dgrad.and..not.o_rbr.and..not.o_orient.and.
     &  .not.o_disp.and..not.o_inv.and..not.o_fail) go to 33
	endif

C---- Ensure that no coordinate projections are attempted
          o_proj=.false.

C---- Open appropriate files and do calculations on the grid. -------
          call open_opts

C read user coordinates
          READ(INFILE,*) nuco
          ALLOCATE(xu(1:nuco),STAT=alloc_stat)
          if(alloc_stat.ne.0) then
            write(*,*) ' Error allocating memory in rtn userco!'
            stop
          endif
          ALLOCATE(yu(1:nuco),STAT=alloc_stat)
          if(alloc_stat.ne.0) then
            write(*,*) ' Error allocating memory in rtn userco!'
            stop
          endif
          ALLOCATE(zu(1:nuco),STAT=alloc_stat)
          if(alloc_stat.ne.0) then
            write(*,*) ' Error allocating memory in rtn userco!'
            stop
          endif
          do 198 iu=1,nuco
          read(INFILE,*,end=199) iuco,xu(iu),yu(iu),zu(iu)
          itot=iu
198       continue
199       continue
c check if correct number of user coordinates
          if(itot.gt.nuco.or.itot.lt.1) then
                write(6,*) 'INCORRECT NUMBER OF USER COORDINATES'
                write(6,*) 'EXPECTED = ',nuco,' ACTUAL = ',itot
          endif

C---- Do calculations on the grid and write the results. -------	      	  
	  nw=nuco
	  do 2 j=1,nw
		x = xu(j)
		y = yu(j)
		z = zu(j)
                percent=float(j)*100./float(nw)
		write(*,*)'working.....',percent,' percent done'
  		call XYZ_Results(x,y,z,NPLANE,NUM_Ds)
 
 2	  continue
 
 
C---- Write the element displacements and Vector output if requested  	      	  
C     and close all files. -------	      	  
33    call final_write(flag,NPLANE,NUM_Ds,nw,1,1)
      deallocate(xu)
      deallocate(yu)
      deallocate(zu)

      return
      end
