!---------------------------------------------------------------------------------------------------
!  BOP
!
!  MODULE:  Green_sshr
!
      MODULE Green_sshr
!
!  PUBLIC TYPES
!	implicit none			! uncomment this later
      INTEGER, PARAMETER :: double = SELECTED_REAL_KIND(p=12)
      SAVE

!  PUBLIC DATA MEMBERS
      REAL(KIND=double), ALLOCATABLE, DIMENSION (:,:) :: sshr  ! Greens function that relates
!                unit slip on cell of first subscript to stress on cell of second subscript
!      REAL(KIND=double), ALLOCATABLE, DIMENSION (:) :: tauvp2  ! Not needed? ++ TET fix
      INTEGER, ALLOCATABLE, DIMENSION (:,:) :: igf_test        ! Flag to say if given Greens 
!                function has been calculated yet ot not
      REAL(KIND=double) :: shear_stress  ! Same as above, but a scalar for MP verion (MMP=1) 
!                and used only to calculate the effect of BC slip on all other elements only 
!                one time in main program (this might change if use arrary storage for sshr 
!                in future version of MP version)
!	<type> ::  <variable>			!  Variable description

!  DESCRIPTION:
!	To share data among the main program and several subroutines
!
!  EOP
!-----------------------------------------------------------------------------------------------------
      END MODULE Green_sshr
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
c     Program PARK
c
c     Earthquake simulation program
c
c     Link it with Ctm.o, upcase.o, and appropriate DERIVS subroutine file
c
c         This program simulates earthquakes using rate and state 
c           friction, the friction being incorporated in the DERIVS
c           subroutine file. Other friction could be used with appropriate
c           changes in that subroutine, althought the program is designed
c			to deal with a friction constitutive law in which slip velocity
c			and shear stress are related, so a simple static/dynamic friction
c			or a slip-weakening with no velocity dependence presumably would 
c			require changes in the main program too. The parts of the program
c			that assign the distribution of rate and state parameters on the
c			fault plane would also have to be modified if altrenate friction
c			laws were used.
c		  The programuses the boundary element method in which slip on each
c			cell (element) contributes to stress on every other cell, so with
c			n cells there are n-squared relations, Greens functions, that
c			relate unit slip on each cell to stress on each cell. Stress 
c			builds up due to slip at an assigned rate on certain boundary 
c			condition cells, and is relieved by slip, either slow or rapid, 
c			on every other cell.
c		  Initial conditions are also required. Presently they are set so that
c			the cells that are velocity weakening (negative A-B in the rate and
c			state constitutive law) are set to be slipping at some fraction of
c			1 m/s and the velocity strengthening cells are slipping at the
c			rate of plate motion, namely the slip rate assigned to the boundary
c			condition cells, which all have the same slip rate in the model.
c		  The program is designed to use a Fast Multipole Method to deal with
c			the fact that with n cells there are n squared Greens functions to
c			treat the problem "exactly," and with large n the compute time to do
c			this is prohibitive. In the Fast Multipole Method, instead of using 
c			all n souce cells to determine the stress at each observation cell,
c			the source cells far away are grouped together and the area-weighted
c			slip on the group is used as the source contribution. The farther 
c			away the cells, the more grouping occurs, in a heirarchical manner. The
c			term "Multipole" comes from the idea of using Greens functions for
c			higher order poles (e.g dipoles, quadrploes in electromagnetics) to
c			approximate the non-uniform distribution of slip on the grouped cells.
c			However, experience shows that it is typically more effricient to simply
c			use more and smaller groupings than to use more and higher order
c			Greens's functions (and for our problem higher order Greens functions
c			may not have been derived anway). There are a variety of things that
c			makes the Multipole Method used here "Fast." These include the option of 
c			using an analytical analytical approximation to the Greens functions 
c			(that fall off as 1/r^3) to decide when the next higher order of 
c			grouping should be done, and numbering the grouped cells so 
c			that cells that are close in space are close in memory.
c		  There is an input parameter "MMP" that allows one to choose to use the
c			multipole version of the program or to skip that and to simply do the 
c			n-squared case. This parameter is used in "if" statements to choose 
c			alternate parts of the program.
c         The program uses radiation damping as a quasi-dynamic
c           approximation to dealing with elastodynamics. This accounts
c           for inertial forces retarding the slip velocity at each
c           point on the fault, but does not include wave effects.
c         This program is presently specificlly focused on modeling
c           earthquakes at Parkfield, California, hence its name.
c			However, there is nothing about it that requires it be used
c           only for that. 
c		  The program presently is set up to deal with a vertical
c           planar strike slip fault in which only the strike slip component
c           of slip and shear stress is involved. Generalizing it for
c           inclined, non-planar, fault with more than those components
c           of slip and stress would require added Greens functions and
c           more complicated sets of differential equations.
c		  The Greens functions that relate the slip on a source cell to the
c			stress at the observation cell comes from the dislocation equations
c			from Mike Chinnery (1963) and 1983 personal communication to W. D. Stuart.
c		  The forward intergating of the coupled ordinary differential equations
c			is presently done using the 5th order Runge Kutta programs from
c			Numerical Recipies from Press et al. (subroutines RKQC and RK4).
c			This can be improved upon - the Runge Kutta approach is probably
c			the best one, but more efficient routines exist.
c
c     The geneology of the program is as follows:
c        1. Began with a set of programs written by William D. Stuart of 
c           the U.S. Geological Survey in about 1988, using the approach of
c			Tse and Rice (1986).
c        2. Modifications were made by Terry E. Tullis of Brown Univeristy
c           during the 1990's, primarily to change the way in which rate
c           and state friction parameters are assigned to different parts
c           of the fault surface, trying to remain as consistant with
c           laboratory values of friction paramenters as possible.
c        3. Radiation damping ws added by Naouyki Kato of the Geologcial 
c           survey of Japan and the Earthquake Institute of the University 
c           of Tokyo while he was a visiting scholar at Brown Univeristy,
C			following the approach of Rice (1993).
c        4. Incorporation of Fast Multipole methods was done with the help
c           of John Salmon of Sandpiper Networks, using a library deveolped
c			by him and Michael Warren of Los Alamos National Laboratory.
c        5. The programs resulting from this collection of workers has
c           been cleaned up and additionally commented by Terry Tullis.
c     References relevant to stages 1 and 2 in the above geneology are:
c         Blanpied, M.L., D.A. Lockner, and J.D. Byerlee, Fault stability 
c			inferred from granite sliding experiments at hydrothermal 
c			conditions, Geophys. Res. Let., v. 18, no. 4, pg. 609-612, 1991.
c		  Chinnery, M.A., The stress changes that accompany strike-slip 
c			faulting, Bulletin of the Seismological Society of America, 
c			v. 53, pg. 921-932, 1963.
c		  Tse, S. T. and J. R. Rice, Crustal earthquake instability in 
c			relation to the depth variation of frictional slip properties,
c			J. Geophys. Res., v. 91, pg. 9452-9472, 1986.
c         Stuart, W. D. and T. E. Tullis, Fault model for preseismic
c           deformation at Parkfield, California, J. Geophys. Res., v. 100,
c           pg. 24079-24099, 1995.
c         Tullis, T. E., Rock friction and its implications for earthquake
c           prediction examined via models of Parkfield earthquakes, in
c           "Earthquake Prediction: The Scientific Challenge" ed. L. Knopoff,
c           K. Aki, J. R. Rice, and L. R. Sykes, Prod. Nat. Acad. Sci. USA, 
c           v. 93, pg. 3803-3810, 1996.
c     A reference relevant to stage 3 above is:
c         Rice, J. R., Spatio-temporal complexity of slip on a fault, 
c           J. Geophys. Res., v. 98, pg. 9885-9907, 1993.
c
c	  Requirements for setting up the geometry of the cells:
c       Let no cell centers be along extensions of the rectangular cell edges 
c		(four dislocation lines define the rectangles) because singularities
c		exist along those lines, and the stresses are calculated at the cell
c		centers. This is most easily ensured by changing cell sizes by odd 
c		factors, e.g., use dx(i+1)=3*dx(i) (or 5 instead of 3 - any odd number).
c		The boundary condition cells are excepted from this since stresses are
c		not calculated at the center of bc cells.
c
c     Run err "log 0 or neg" in subr derivs means time step is too big
c       eg dt1 too big or stable faulting w/ big dt; reduce eps or dt1,
c       or improve definition of tvscal; or make ac bigger
c
c     t step is just a counter
c
c     units:    displ,mm   dist,km   vel, mm/yr   time, years
c               stress, bars   strain, microstrain 
c               rigidity, mbars   constit law const, bars
c
c     files:    2-prk.dat      
c				4-prk.flt
c               5,6-in,out term  
c				7-prk.prt  
c				8-prk.sum       
c				9-prk.crp
c
c -- input data in .dat file --
c
c     title

c     < >       <block geometry info for prktr1>
c               <first cell block encompasses patch, unstable slip>

c     ndt       max no. of t steps
c     nbc       =1, some cells are b.c. vpl, distant fault locked
c               =0, entire fault ourside patch to infinity slips at vpl
c                   ie. savage procedure of removing rigid slip
c     npat      no. fault of patches, only npat=1 works 

c     nwp1      (7) write print file prk.prt starting  w/ 
c                   t step nwp1; ; nwp1=0 for no print
c     wpdt      write .prt at intervals ge wpdt yrs
c               if wpdt=0, write every t step

c     nwc1      (9) write surface creep rate prk.crp
c                   >0 write; nwc1=0 no write
c     wcdt      write .crp at intervals ge wcdt yrs

c     nwf1      (4) write binary file prk.flt for plotting starting w/ 
c                   t step nwf1; ; nwf1=0 for no write
c     wfdt      write .flt at intervals ge wfdt yrs
c               if wfdt=0, write every t step

c     dt1       trial dt for first t step, yr
c     vlim      instab defined if any v>vlim, mm/yr
c                 1 m/sec = 3.156 e10 mm/yr
c     eps       max err, see tvscal, subr rkqc  ca. 1.e-4 or 1.e-5
c     vpl       plate velocity, >0, =vstar 
c	  factinit  factor times vlim that is used for initial condition of
c				velocity in velocity weakening area - typically use 0.01

c     ngrp      no. of cell range pairs for cells having vpl bndy. cond.

c     ibca,ibcb indices of first, last cells having vpl, 4 pairs max

c               patch info read by subr patch:

c     patnam    name of patch subr: patch4-wds  patch5-tet lab
c     pxio,pzo  position of patch center, km
c     xmn,xmx   x range, km; used for multiple patches
c     zmn,zmx   depth range, km
c     ax,az     patch size, half width, half height, km

c     ac        val of A at z=0, ca. 0.4, >0; A inst strength, bar
c     arat      linear rate of increase of A w/ z, bar/km
c     abmn,abmx min, max value of A-B, bar;  <0 brittle,   >0 visc
c     xldis     slip weakening displ, mm

c     ncr       no. of trace or other cells vor rate prnt on .crp
c     ncrp      list of the cell indices

c ------------------ other variables -----------------------------

c     nf        no. fault cells, bc + non-bc
c     ne        no. non-bc cells, = no. pde's
c     nb        no. of bc cells, nf=ne+nb
c     t         elapsed time, yr
c     dt        time step size, years
c     x,y,z     (x1,x2,x3) x1 +rt,par flt  x2 nor flt  x3 dwn
c     tv(2ne)   first ne are tau, 2nd ne are vel
c     dtvdt     derivatives of tv wrt/ t
c     u(nf)     fault slip u=u(+)-u(-), left lat pos
c     du        slip during instab, same time; also ea. t step
c     sshr      (nf,nf) - for no MP version - green functions fault shear stress
c     mapx,mapz =no. chars x,y dir for coseis plot  subr seis plot
c
c   NEXT VARIABLES ARE FROM THE prkf.input.run FILE:
c     nsteps - number of time steps at which formatted I/O is wanted
c     irestart - 1 if this run is to start from some time step with data output
c                form an previosu program, 0 if just starting from beginning
c     ireend - 1 if this run is to output data when it stops at the end that 
c              will allow restarting a subsequent run from that endpoint, 0 if 
c              doesn't set up for a later run to continue
c     ninit - step number to begin from if restartign at the end of a previous 
c             run (used if irestart=1)
c     nend - step number to write out data near end of run if setting up for a 
c            laterrestart (used if ireend=1)

C     Dimensioning changes as change model:     ! TET Mod
C        nf (or ne?) is the dimension of many arrays: currently is 15003
C                                                   2427 was too big for 96Megs
C        2nf (or 2ne?) is the dimension of a few arrays: currently is 30006
C
C        nblk is the dimension of several arrays: currently is 40
C
C        The number of neighbor cells (+1) that a cell can have was set 
C        to 10 by WDS, but I changed it to 25 to accomodate a 9x reduction
C        at a boundary on up to two sides of a cell, plus 3x on the other
C        two sides and one for the number of entries in the array. This is
C        the second dimension of nabor( , ).

      USE Green_sshr		! Make data in module "Green_sshr" available
      implicit real*8 (a-h,o-z)
      INTEGER :: a_and_amb      ! Flag to say if a and amb are output for plotting check
      real*4 sa,sab
      real*4 seps,svln,t4
      parameter (NF_MAX=15003, NE2_MAX=30006)
      character*5  crpnam(10)
      character*6  patnam
      character*30 timedate
      character*80 title1
      character*80 title
      character*30 inputfile,datafile,prtfile,sumfile,fltfile,
     &          crpfile,afile,ambfile,alleqtimes,restarbegin,restarend
      character*30 outputfile,clocktimefile,tecplotfile
      dimension out(100000,4,7)                        ! ++ increase if over 100,000 time steps
      dimension w(NF_MAX),d(NF_MAX),ibca(4),ibcb(4),ibctmp(NF_MAX),
     &                 celar(NF_MAX)
      dimension u(NF_MAX),uold(NF_MAX),du(NF_MAX)
      dimension tv(NE2_MAX),dtvdt(NE2_MAX),tvscal(NE2_MAX),
     &          tvold(NE2_MAX)
      dimension ncrp(10),nstep(10)
      dimension nabor(NF_MAX,25)
      dimension svln(NF_MAX)
      dimension g1(10),g2(10) ! Increase if use > 10 random surf. pts
      common /deriv/ tauvp(NF_MAX),ca(NF_MAX),cab(NF_MAX),dx(NF_MAX),
     . dz(NF_MAX),tauss(NF_MAX),sa(NF_MAX),sab(NF_MAX),
     . e10,xldis,vpl,xmu,ne,ne2,MMP
      common /pat/ xc(NF_MAX),zc(NF_MAX),vlim,npat,nb,ibc(20),
     & iec(NF_MAX),iecinv(NF_MAX)

c     zc,xc, are in common /pat/ dx,dz, are in READ_CELL_INFO, sshr is in MODULE Green_sshr

      CALL READ_FILE_NAME(inputfile)

      print *,' number of steps between output writes for clocktime !! ' !  NK
	  read (*,*) ioutno
c
      open(unit=1,status='old',form='formatted',file=inputfile)
      read (1,2) title1
      read (1,*) numlines
      read (1,'(a)') datafile      ! 1 unit=2
      read (1,'(a)') outputfile    ! 2 unit=26
      read (1,'(a)') clocktimefile ! 3 unit=27
      read (1,'(a)') prtfile       ! 4 unit=7   
      read (1,'(a)') sumfile       ! 5 unit=8
      read (1,'(a)') fltfile       ! 6 unit=4
      read (1,'(a)') crpfile       ! 7 unit=
      read (1,'(a)') afile         ! 8 unit=
      read (1,'(a)') ambfile       ! 9 unit=
      read (1,'(a)') tecplotfile   ! 10 unit=50
c     ---- skip lines used only by later surface strain program that also reads this file ---
      do k=1,numlines-12           
         read (1,*)
      end do
      read (1,'(a)') restarbegin   ! file to read that contains the info. needed for this execution to restart in the middle of a simulation
      read (1,'(a)') restarend     ! file to write that contains the info. needed for a subsequent execution to restart where this one stops
      read (1,*) nsteps, ne_dumm, vpl_dumm, a_and_amb, MMP   ! number of steps for formatted I/O; 2nd and 3rd used in later strain program
      n=-4
      do i=1,nsteps,5              ! read in steps for formatted I/O, 5 per row
	      n=n+5
	      m=n+4
	      read (1,*) (nstep(j),j=n,m)  ! read in specific step #s for writing formated I/O
      end do                       ! end lines used only by next program
      read (1,*)                   ! skipping line giving strain component
      read (1,*)  xog,yog,dxx,dy,nx,ny,nbm   ! obs pts on grid for later surface strain program
      if (nbm .ne. 0) then                   ! obs pts irregular locations for surface strain by later program
         read (1,1009) (g1(i+nxy),i=1,nbm)
         read (1,1009) (g2(i+nxy),i=1,nbm)
      end if
      read (1,*)irestart,ireend,ninit,nend   ! determine whether restart file is read or written and step numbers at which this is done
      if(irestart.eq.0) ninit=1
c
      write (6,2) title1
      open (unit=2,status='old',form='formatted',file=datafile)
      open (unit=7,status='new',form='formatted',file=prtfile)
      open (unit=8,status='new',form='formatted',file=sumfile)
      open (unit=26,status='new',form='formatted',file=outputfile)
      open (unit=27,status='new',form='formatted',file=clocktimefile)
      open (unit=50,status='new',form='formatted',file=tecplotfile)

      if(irestart.eq.1)open (unit=12,status='old',form='unformatted',
     &   file=restarbegin)
      if(ireend.eq.1) open (unit=13,status='new',form='unformatted',
     &   file=restarend)

      timedate='                              '  ! Initialize w/ 30 blanks
      call cprogtime(timedate)
      nwf1cnt = 0                                ! count t steps written
c
      Call READ_CELL_INFO      (ibca,ibcb,w,d,dx,dz,title,xmu,
     .nblk,nf,ndt,nbc,nwp1,wpdt,nwc1,wcdt,nwf1,wfdt,dt1,eps,
     .vpl,factinit,ngrp) 

c     -----
      seps=sngl(eps)
      dt=dt1                         ! time step size, for tvscal
      t=0.d0                         ! initial value of time
      tiny=1.e-30                    ! for rkqc

      write (6,*)
      write (6,*) 'a_and_amb =',a_and_amb,'  MMP =',MMP
      write (6,*)
      write (6,1310)  timedate
 1310 format (1x,a30)
      write (6,5)
      write (6,14)  ndt,nbc,npat,nwp1,wpdt,nwc1,wcdt,nwf1,wfdt,nf,
     . mapx,mapz
      write (6,11)
      write (6,32)  vlim,vpl,xmu,dt1,seps,factinit

      write (7,*)
      write (7,*) 'a_and_amb =',a_and_amb,'  MMP =',MMP
      write (7,*)
      write (7,1310)  timedate
      write (7,5)
      write (7,14)  ndt,nbc,npat,nwp1,wpdt,nwc1,wcdt,nwf1,wfdt,nf,
     . mapx,mapz
      write (7,11)
      write (7,32)  vlim,vpl,xmu,dt1,seps,factinit

      write (8,17)
      write (8,2)   title
      write (8,*)
      write (8,1310)  timedate
      write (8,*)
      write (8,*) 'a_and_amb =',a_and_amb,'  MMP =',MMP
      write (8,5)
      write (8,14)  ndt,nbc,npat,nwp1,wpdt,nwc1,wcdt,nwf1,wfdt,nf,
     . mapx,mapz
      write (8,11)
      write (8,32)  vlim,vpl,xmu,dt1,seps,factinit


      if (nwf1 .eq. 0) go to 22
      open (unit=4,status='new',form='unformatted',file=fltfile)
      write (4)  title,timedate,ndt,nf,seps

c     ----- boundary conditions, cell slip vel -----

   22 write (6,9)
      write (7,9)
c
c	Question - what is the difference between cell numbers and cell indices?
c	Answer - They are the same thing. The cell index or cell number is the 
c	         sequential numbers given to the cell when it was initially defined 
c            in the sequence of defining the locations of ALL the cells
c          - The other sort of list that exists is the sequence number of the
c   	     non-bc cells (ne of them), or the sequence number of the bc cells
c            (nb) of them. In a typical definiton of bc cells where they are all
c            put at the end of the list of ALL cells, the cell number or index of
c            the non-bc cells would be the same as the sequence number of those
c            cells in the ne-long list of non-bc cells. In this case, where also
c            typically there are only a few bc cells (say, e.g. 5) the sequence
c            number of the be cells would range from 1-5 but the cell numbers for
c	         those bc cells would be larger than the sequence number by ne, the 
c            number of non-bc cells, all with lower number.
c
c     -- store indices of bc cells

      nb=0

      write (6,16) ngrp,(ibca(j),ibcb(j),j=1,ngrp)

      do 28 j=1,ngrp
         l=ibcb(j)-ibca(j)+1      ! l is number of bc cells in this group
         do 28 k=1,l
            nb=nb+1
            ibc(nb)=ibca(j)+k-1   ! ibc is the index of the ith one of the nb bc cells
            ibctmp(nb)=ibc(nb)	  ! ibctmp is used in the KOMPRS routine
   28 continue

c     -- store indices of non-bc cells

   26 ne=0
      do 30 n=1,nf
         do 29 j=1,nb
            if (n .eq. ibc(j)) go to 30
   29    continue
         ne=ne+1
         iec(ne)=n      ! iec is the index of the ith one of the ne non-bc cells
         iecinv(n)=ne   ! iecinv(i)is the sequence number in the list of the ne non-bc cells
						!        that goes with cell number i	
   30 continue

      fltarea=0.                             ! area of non bc part of fault
      do j=1,ne                              ! Note: WDS had sum over bc cells too
          n=iec(j)
          celar(n) = dx(n)*dz(n)             ! cell area
          fltarea  = fltarea+celar(n)
      end do

      write (6,*)
      write (6,*)'fault area (sq km) =', fltarea
      write (7,*)
      write (7,*)'fault area (sq km) =', fltarea

      write (6,*)'from main: npat=',npat

      write (7,31) nb
      write (7,1)  (ibc(i),i=1,nb)    ! list of bc cell numbers or indices, 
                                      !   referenced by sequence number in the nb 
                                      !   long list of bc cells
      write (7,33) ne
      write (7,1)  (iec(i),i=1,ne)    ! list of non-bc cell numbers or indices
                                      !   referenced by sequence number in the ne 
                                      !   long list of non-bc cells
      write (7,*) ' iecinv'
      write (7,1)  (iecinv(i),i=1,ne) ! list of non-bc sequence numbers (iec indices)
                                      !   referenced by cell number in the ne 
                                      !   long list of non-bc cells 
      if (nwf1 .eq. 0) go to 34

   34 ne2=2*ne
      nev1=ne+1
      nev2=ne2

c	  The first subscript in the sshr array is for the source point, 
c		the second for the observation point.
c		Both of these need to be nf for noMP version. They are both set
c	        just below to 1 in the present MP version, since sshr not used, and
c		the array is huge if it is n by n and n is large. 
c		However, in a future version, hopefully can set the second to n and the first to <<n if 
c		can figure out how big the first needs to be from the MP library and then store the needed 
c		Greens functions in this smaller array. If the first can be small enough, then may be 
c		able to keep sshr all in memory and won't have to compute it every single call to DERIVS
c		and so execution time may go way down. John Salmon says some people at Yale
c		are working on something that is like this, so need to investigate that.

      if (MMP .eq. 1) then
          ALLOCATE (sshr(1,1), STAT=isshr_status)      ! degenerate array for MP
      else
          ALLOCATE (sshr(nf,nf), STAT=isshr_status)    ! nf by nf array for noMP
      end if

      if (sshr_status .ne. 0) then
          write(6,*) 'stop because isshr_status = ',isshr_status,
     &                                                '; insuf memory?'
          STOP
      end if


c     -- calc stress at non-bc cells due to vp over only bc cells
      do 54 i=1,ne				 ! loop over all non-bc cells as obs pts
         tauvp(i)=0.
         if (nbc .eq. 0) go to 54       ! rigid displ, inf fault
         n=iec(i)                ! iec(i) is cell number of ith non-bc cell
         dtv=0.
         do j=1,nb
             m=ibc(j)			 ! ibc(j) is cell number of jth bc cell
	     call GREENS_FUNCTION(m,n)        ! m is the source, n is the obs pt
             if (MMP .eq. 1) then             ! Don't use array if MP version
                 dtv = dtv + shear_stress     ! shear_stress is a scalar, reused in MP version
             else
                 dtv = dtv + sshr(m,n)        ! m is the source, n is the obs pt
             end if
         end do
         dtv = vpl*dtv
         tauvp(i) = tauvp(i) + dtv
   54 continue

      if (MMP .ne. 1) then  ! Init flag to determine if a particular Green's Function has been calculated
          ALLOCATE (igf_test(nf,nf), STAT=igf_test_status) ! nf by nf array for noMP
          if (sshr_status .ne. 0) then
              write(6,*) 'stop because isshr_status = ',isshr_status,
     &                                                '; insuf memory?'
              STOP
          end if

          Do n=1,ne                                                      
	      Do i=1,ne                                                    
                  igf_test(i,n)=1
	      End do
          End do
      end if

      nbb=nb
	  if (MMP .ne. 1) then  ! don't need to Compress or use sshr again if is MP version, i.e. if MMP=1
c	      ++ TET fix - Note that if later use stored array of sshr in MP version, may want to Compress
          if (nbc .gt. 0) CALL KOMPRS (nbb,ibctmp,nf)  ! This removes bc rows and columns from sshr
      end if

C     These next 4 lines are to check if a and a-b are done correctly - a and amb written later

      if(a_and_amb .eq. 1) then      
          open (unit=20,status='new',form='unformatted',file=afile)
          open (unit=21,status='new',form='unformatted',file=ambfile)
          write (20)  title,timedate,ndt,nf,seps
          write (21)  title,timedate,ndt,nf,seps
      end if

c     These not used in current version, but useful in future for different way of assigning friction
c      if (patnam .eq. 'patch4') CALL PATCH4 (tv)  
c      write (6,*)'from main: npat=',npat
c      if (patnam .eq. 'patch5') CALL PATCH5 (tv,1,1,irestart,factinit,a_and_amb) ! line will be too long

      CALL PATCH5 (tv,1,1,irestart,factinit,a_and_amb,nminc)

      read (2,*)   ncr                     ! for .crp file if calculate creep at surface
      read (2,23) (crpnam(n),n=1,ncr)      ! have already read in rest of data in subroutines
      read (2,4)    (ncrp(n),n=1,ncr)
      close (unit=2)

      write (6,*)

      tplast = 0.
      tclast = 0.
      tflast = 0.

      if (nwc1 .eq. 0) go to 44
      open (unit=9,status='new',form='formatted',  file=crpfile)
      CALL CREEP (nt,t,tv,xc,ne,1,ncr,ncrp,crpnam)

c     This section is used for restart possibility
 44   if(irestart.eq.0)then
	      do i=1,ne
	          u(i)=0.
          end do
      else !if restarting from earlier run, read all but tv (tv read earlier in init)
	      read (12) nninit
	      if(nninit.ne.ninit)then
	          write (6,*)'Stopping because ninit not equal in',
     &                   ' prkf.input.run and restartbegin files'
	          STOP
	      end if
	      read (12) t
	      read (12) dt
	      read (12) (u(i),i=1,ne)
	      CLOSE (UNIT = 12)
      end if

      write (6,41)
   41 format ('  t step      t         dt     v ave (mm/yr)')
      mode  = 0 ! mode is left over from a previous version; written in .flt file so keep for now
c
c	THIS IS THE START OF THE MAIN TIME STEP LOOP - subr rkqc integrates to next time
c
      do 100 nt=ninit,ndt                ! allows restart at step # ninit 
c           the following is to allow later restarting at step # nend
	  if(ireend.eq.1.and.nt.eq.nend) then  ! set up for restart if desired
	      write (13) (tv(i),i=1,ne2)
	      write (13) nt
	      write (13) t
	      write (13) dt
	      write (13) (u(i),i=1,ne)
	      CLOSE (UNIT = 13)
	  end if
c
      ntm   = nt-1                ! save values of last t step ++ fix
      tm    = t
      dtold = dtdid
      smomold = smom
      uavold  = uav

      do 52 i=1,ne
          k=i+ne
          du(i)    = 0.
          tvold(i) = tv(i)
          tvold(k) = tv(k)
   52     uold(i)  = u(i)

      CALL DERIVS (tv,dtvdt)

      if (nt.eq.ninit) then
         write (7,7)              ! writes heading for table to follow
         iflag=0
         iflagpos=0
         iflagneg=0
         do 46 i=1,ne
            n=iec(i)

	    if (MMP .ne. 1) then
                write (7,35) n,xc(n),zc(n),tauvp(i),sshr(i,i)  ! Will give good values, since sshr just calculated
	    end if

C           The next lines are determining the cell mumber where A-B is nearly the 
C             most pos within the to the North of the same depth where it is nearly the most 
C             negative. The cell where it is the most negative already has been 
C             determined in PATCH5 and it is cell nminc. 
C           The purpose of this is to compare behavior between various models with 
C             different numbers of cells, different values of the analytic error 
C             function for the Multipole method, Multipole, vs. noMultipole, etc.
c           The time history at this pos A-B cell and the most neg A-B cell are
C             printed out near the end of the large time loop as an ascii file at
C             every time step for plotting by Tecplot.
            if(iflag.eq.0 .and.xc(n).ge.-5.5 .and.zc(n).ge. zc(nminc))    ! is within smaller cell groups
     &                                                              then
               nmaxc=n
               iflag=1
            end if
C           The next lines are setting up to print out a time history at the point where
C             A-B is nearly the most neg and to the North at the same depth where it is 
C             nearly the most positive within the group of smallest cells. This will be
C             removed soon, and is replaced by the above if statement, but is here still 
C             to compare with previous plots. compare behavior between various models with ndfferent numbers of cells,
C             different values of the analytic error function for the Multipole method,
C             Multipole, vs. noMultipole, etc.
            if(iflagpos. eq.0.and.xc(n).ge.-5.5 .and.zc(n).ge. 7.5) then
               icellpos=n
               iflagpos=1
            end if
            if(iflagneg. eq.0.and.xc(n).ge.13.5 .and.zc(n).ge. 7.5) then
               icellneg=n
               iflagneg=1
            end if
   46    continue
      end if

      do 87 n=1,ne2
   87     tvscal(n) = abs(tv(n)) + abs(dt*dtvdt(n)) + tiny  ! ++ improve 

      CALL RKQC   (tv,dtvdt,ne2,dt,eps,tvscal,dtdid,dtnxt)  ! Integrate to the next time

      t  = t+dtdid
      t4 = sngl(t)

   62 do 64 i=1,ne
          k=ne+i
          du(i) = .5*dtdid* (tvold(k)+tv(k)) ! ++ ? replace w/ higher order
   64 u(i)  = u(i) + du(i)
      
c     ----- calc total moment, non-bc cells, dyn-cm -----

   66 vav = 0.                           ! average fault velocity
      uav = 0.                           ! (upl = t4*vpl)

      do 68 i=1,ne                       ! non bc cells
         n   = iec(i)  ! ++ TET fix - not sure of effect of KOMPRES on this in noMP version
         vav = vav + (tv(i+ne)-vpl)*celar(n) 
   68 uav = uav + u(i)*celar(n)

      vav  = vav/fltarea
      uav  = uav/fltarea
      smom = 1.e21*xmu*uav

c     ----- write t, fields, this t step -----

c     -- write .prt every wpdt yrs, just after instab, nf1, and ndt
c        wfdt=0, write every t step
c        nyes=1 returned if time since last write, tlast, ge dt years

c     -- .flt, 4

      if (nwf1 .eq. 0 .or. nwf1 .gt. nt) go to 76
      if (  nt .eq. 1 .or. nt .eq.  ndt) go to 81
      if (wfdt .eq. 0.) go to 81
      CALL WRSTEP (t,tflast,wfdt,nyes)
      if (nyes .ne. 1) go to 76

   81 nwf1cnt = nwf1cnt + 1
      write (4) nt,t,mode     ! ++ TET fix - mode is left over from a previous version; keep for now 
      do 1375 i=1,ne
	      svln(i)=sngl( log(tv(i+ne)/vpl) ) ! sngl. prec. output nat. log normalized vel
 1375 continue
      write (4) (svln(i),i=1,ne)

c     -- .log, 6

   76 if( mod( nt, ioutno ) .eq. 0 .or. nt .eq.1 ) then
          call cprogtime(timedate)
          write (27,*) nt,t,dtdid,timedate ! ++ TET fix - need to think abt. what is written out in 26 & 27
      end if

      write (6,65) nt,t,dtdid,vav

      if( mod(nt,100).eq.0 .or. nt.eq.1 .or. nt.eq.ndt) then      ! NK mod
          call cprogtime(timedate)
          write (26,*) nt,t,dtdid,timedate
          year1sec = 3.1536e7
          v100 = tv( ne + 100 ) / year1sec
          v200 = tv( ne + 200 ) / year1sec
          v500 = tv( ne + 500 ) / year1sec
c	  find cell with maximum velocity
          vmax00 = tv( ne + 1 )
	      iimax = 1
	      do iii = 2, ne
	          if( tv( ne + iii ) .gt. vmax00 ) then
		          vmax00 = tv( ne + iii )
		          iimax = iii
	          end if 
	      end do
	      vmax00 = vmax00 / year1sec 
c	       print *,' v100=',v100,' v200=',v200,' v500=',v500
c	       print *,'          vmax=', vmax00,' at i =',iimax  ! NK
	      write (26,*)'          vmax=', vmax00,' in cell ',iimax  ! NK
C	      write (26,*) (tv(iii), iii=1,ne)            ! TET C'd out
C	      write (26,*) (tv(ne+iii)/year1sec, iii=1,ne) ! TET C'd out
      end if

c     -- .prt, 7

      if (nwp1 .eq. 0 .or. nwp1 .gt. nt)   go to 72
      if (nt .le. 2 .or. nt .eq. ndt)      go to 73
      if (wpdt .eq. 0.) go to 73        ! This means write each step
      CALL WRSTEP (t,tplast,wpdt,nyes)  
      if (nyes .ne. 1) go to 72         ! skip printing if not time yet

   73 write (7,67) nt,t,dtdid
      write (7,69)
      do 70 i=1,ne
          write (7,71) iec(i),tv(i),tv(ne+i),u(i),tauss(i)
   70 continue

c     -- .crp, 9

   72 if (nwc1 .eq. 0)  go to 90
      CALL WRSTEP (t,tclast,wcdt,nyes)  
      if (nyes .ne. 1) go to 90
      CALL CREEP  (nt,t,tv,xc,ne,0,ncr,ncrp,crpnam)

   90 dt=dtnxt

      out(nt,1,1)=float(nt)
      out(nt,1,2)=t
      out(nt,1,3)=float(nmaxc)
      out(nt,1,4)=tv(nmaxc)
      out(nt,1,5)=tv(ne+nmaxc)
      out(nt,1,6)=u(nmaxc)
      out(nt,1,7)=tauss(nmaxc)

      out(nt,2,1)=float(nt)
      out(nt,2,2)=t
      out(nt,2,3)=float(nminc)
      out(nt,2,4)=tv(nminc)
      out(nt,2,5)=tv(ne+nminc)
      out(nt,2,6)=u(nminc)
      out(nt,2,7)=tauss(nminc)

      out(nt,3,1)=float(nt)
      out(nt,3,2)=t
      out(nt,3,3)=float(iec(icellpos))
      out(nt,3,4)=tv(icellpos)
      out(nt,3,5)=tv(ne+icellpos)
      out(nt,3,6)=u(icellpos)
      out(nt,3,7)=tauss(icellpos)

      out(nt,4,1)=float(nt)
      out(nt,4,2)=t
      out(nt,4,3)=float(iec(icellneg))
      out(nt,4,4)=tv(icellneg)
      out(nt,4,5)=tv(ne+icellneg)
      out(nt,4,6)=u(icellneg)
      out(nt,4,7)=tauss(icellneg)

  100 continue

C     END OF MAIN LOOP OVER TIME

      do 101 j=1,2
         write (50,*)'ZONE'                    ! writing for a tecplot file
         do 101 i=1,nt-1
 101         write (50,*) (out(i,j,k),k=1,7)  ! writing for a tecplot file
      CLOSE (UNIT = 50)

      do 103 j=3,4               ! writing to fort51 on purpose - will remove soon
         write (51,*)'ZONE'                  ! writing for a tecplot file
         do 103 i=1,nt-1
 103         write (51,*) (out(i,j,k),k=1,7) ! writing for a tecplot file

      write (6,102) nwf1cnt
      write (8,102) nwf1cnt
  102 format (/i6,' total t step writes to prk.flt')

c     .......... end time steps ..........

    1 format (1x,8i10)
    2 format (a80)
    4 format (15i5)
    6 format (1x,15i5)
    8 format (3f10.2,2f5.1)
 1009 format (10f5.1)
   10 format (1x,2a10,5x,2a10)
   13 format (4(i5,f5.1))
   14 format (1x,3i5,3(i5,f7.3),3i5)
   15 format (8i5)
   16 format (i5,8i10)
   18 format (2i5,2f5.0)
   19 format (a6)
   21 format (1x,10i8)
   23 format (15a5)
   32 format (1x,1p10e10.2)
   35 format (i10,1x,2f8.4,3f8.3)
   65 format (1x,i7,f11.4,2ES10.2)
   71 format (i10, f10.3,1pe10.2, 0pf8.1,f8.2) 
   17 format (' prk.sum'/)
    5 format (/'   ndt   nbc npat nwp1  wpdt  nwc1  wcdt  nwf1  wfdt',
     . ' cells mapx mapz')
    7 format (/' non-bc cells:'//'       cell    xc      zc     tauvp',
     .         '   sshr')
    9 format (/'  ngrp  <grp1>   <grp2>   <grp3>   <grp4>')
   11 format (/'     vlim      vpl       xmu       dt1       eps       
     .factinit')
   31 format (/1x,i4,' bc cells:')
   33 format (/1x,i4,' non-bc cells:')
   57 format (/' last fields before unstab slip')
   67 format (/' -- t step       t             dt --',
     .        / i8,5x, f17.12, 1p2e9.2)
   69 format (/'       cell     tau       vel      u     tauss')
      stop
      end			    			! This is the end of the main program
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
!---------------------------------------------------------------------------------------------------
!  BOP
!
!  ROUTINE:  READ_FILE_NAME
!
!  INTERFACE:
      SUBROUTINE READ_FILE_NAME (
     &   filenm)         ! filename of file that is being read in from standard input


!
!  USES:
c	use <module>    Not yet converted to Fortran 90/95 use of modules
!
!  RETURN VALUE:
c	implicit none   ! Not yet converted to using "implicit none" and explicitly typing all variables
c	<type> :: <name>			! Description
!
!  PARAMETERS:
c	<type, intent> :: <parameter>		! Description
!
!  NECESSARY PRE-CONDITIONS:
!	e.g., specific values being set in Module global variables; a specific object
!	having been created; etc.
!
!  SIDE EFFECTS:
!	e.g., modifying global variables; creating objects that persist beyond 
!	scope of the routine/function.
!
!  POST_CONDITIONS UPON EXIT
!	e.g., any files being closed, objects being deleted, etc.
!
!  DESCRIPTION:  Describe the function of the routine and algorithm(s) used
!			in the routine.  Include any applicable external references.
!
!  EOP
!---------------------------------------------------------------------------------------------------
      integer*4 sstrlen
      integer*4 lenfilenm,lendefalt
      character*30 filenm,defalt
      character hold*4,string*500
c
      lendefalt=14
      defalt='prkf.input.run'
c                                       prompt user for filename
      write (6,710) defalt(1:lendefalt)
 710  format (' Give name of data file. Default is [',a,']')
      read (5,'(a)') filenm
c
      lenfilenm=sstrlen(filenm,30)
c                                       equate to default if no input
      if (lenfilenm.le.0) filenm(1:30)=defalt(1:30)
      RETURN
      END
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
	function sstrlen(string,ilen)
	integer sstrlen,imax,ilen
	character string*80
	if(ilen.gt.80) then
	   imax=80
	else
	   imax=ilen
	endif
c don't return a negative number:
	if(imax.le.0) then
	   sstrlen=0
	   return
	endif
c count backwards to find the last ascii character in the string:
	i=imax+1
5       i=i-1
	if((i.eq.0).or.((ichar(string(i:i)).gt.32).and.
     .        (ichar(string(i:i)).le.127))) goto 20
	goto 5
20      sstrlen=i
	return
	end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
      SUBROUTINE READ_CELL_INFO(ibca,ibcb,w,d,dx,dz,title,xmu,
     .nblk,nf,ndt,nbc,nwp1,wpdt,nwc1,wcdt,nwf1,wfdt,dt1,eps,
     .vpl,factinit,ngrp) 

c	This subroutine reads in some of the data from the .dat file.
c	 This includes the cell information, the "Simulation Parameters", 
c	 and the list of which cells are the boundary-condition cells. 
c	 The rest of the information from the .dat file, the friction law 
c	 parameters, how they are distributed on the fault and the 
c	 creepmeter names and corresponding cell numbers, are read 
c	 in by the main program.

c  --- input variables this code --------
c     title    title
c     nblk     no. of blocks of identical cells 
c     nx,nz    number of identical cells x,z dir this block
c     dxb,dzb  cell width, height, km for all cells in this block
c              if negative, cell dimensions = 1/|dxb|, 1/|dzb|
c     xul,zul  coords upper left corner this cell block, km

c ---- other variables
c     nf       no. fault cells, total of bc and non bc
c     u        fault slip  u=u(+)-u(-)  lft lat +
c     x,y,z    (x1,x2,x3) x1 +rt,par flt  x2 nor flt  x3 dwn
c     xc,zc    coords cell centers, km
c     dx,dz    full cell width, ht, km
c     xmu      rigidity, poisson solid

C     Dimensioning changes as change model:
C        nf is the dimension of many arrays: currently is 15003
C
C        2ne is the dimension of a few arrays: currently is 30006 
C
C        nblk is the dimension of several arrays: currently is 40
C
      implicit real*8 (a-h,o-z)
      parameter (NF_MAX=15003, NE2_MAX=30006)
      common /pat/ xc(NF_MAX),zc(NF_MAX),vlim,npat,nb,ibc(20),
     & iec(NF_MAX),iecinv(NF_MAX)
      character*80 title
      character*30 inputfile,datafile,gffile,prt1file
      dimension nx(40),nz(40),dxb(40),dzb(40),xul(40),zul(40)
      dimension dx(NF_MAX),dz(NF_MAX),ibca(4),ibcb(4)
      xmu=0.3
      read (2,6)  title
      write (7,17)
      write (7,6)  title
      write (7,*)
      write (6,6)  title
c ---   read cell block info, calc coords cell centers   -------
      read (2,*) nblk
      nf=0
      nc1=1
      write (7,16)
      do 50 n=1,nblk
         read  (2,*) nx(n),nz(n),dxb(n),dzb(n),xul(n),zul(n)
         if (dxb(n) .lt. 0.) dxb(n) = -1./dxb(n)
         if (dzb(n) .lt. 0.) dzb(n) = -1./dzb(n)
         do 42 k=1,nz(n)
            do 40 i=1,nx(n)
               nf=nf+1             ! nf is determined by this process
               dx(nf) = dxb(n)
               dz(nf) = dzb(n)
               xc(nf) = xul(n) + (i-0.5)*dxb(n)
               zc(nf) = zul(n) + (k-0.5)*dzb(n)
   40       continue
   42    continue

         nc2=nf       ! what is the purpose of nc2? May no longer be needed? ++ TET fix
         write (7,13) n,nc1,nc2,nx(n),nz(n),dxb(n),dzb(n),xul(n),zul(n)
         nc1=nc2+1    ! what is the purpose of nc1? May no longer be needed? ++ TET fix
   50 continue
      read (2,*)
      read (2,*) ndt,nbc,npat
      read (2,*) nwp1,wpdt,nwc1,wcdt,nwf1,wfdt
      read (2,*) dt1,vlim,eps,vpl,factinit
      read (2,*)
      read (2,*) ngrp                         ! bc cells
      read (2,*) (ibca(n),ibcb(n),n=1,ngrp)   ! cell number of first (ibca) and last (ibcb)
      write (7,4) ngrp
      write (7,4) (ibca(n),ibcb(n),n=1,ngrp)
      write (6,2) nf
    2 format (/i5,' cells total')
    3 format (2i5,6f5.2)
    4 format (10i7)
    6 format (a)
    7 format (1x,a30)
   13 format (i4,4i7,8f10.2)
   15 format (i4,4f8.3,f10.4)
    1 format (/' cell     x       z      dx      dz    self strs')
   14 format (' u=1 on cell',i5)
   16 format (/' cell block(s):'/' blk     cells         nx     nz ',
     & '     dxb       dzb       xul       zul ')
   17 format (' prk.prt'/)
      return
      end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
 	  SUBROUTINE GREENS_FUNCTION(m,n)
c
c       for unit lft lat slip each cell, find shear stress
c       at all cell centers
c
      USE Green_sshr	! Make data in module "Green_sshr" available
      implicit real*8 (a-h,o-z)   
      real*4 sa,sab
      parameter (NF_MAX=15003, NE2_MAX=30006)
      common /deriv/ tauvp(NF_MAX),ca(NF_MAX),cab(NF_MAX),dx(NF_MAX),
     . dz(NF_MAX),tauss(NF_MAX),sa(NF_MAX),sab(NF_MAX),
     . e10,xldis,vpl,xmu,ne,ne2,MMP
      common /pat/ xc(NF_MAX),zc(NF_MAX),vlim,npat,nb,ibc(20),
     & iec(NF_MAX),iecinv(NF_MAX)

c   Source fault is indexed by m
       x3  = zc(m)
       dx2 = dx(m)/2.
       dz2 = dz(m)/2.
c   Obs point is indexed by n
       y1 = xc(n)-xc(m)
       y3 = zc(n)
       call STRN (dx2,dz2,x3,y1,y3,e12)
	   if (MMP .eq. 1) then
               shear_stress = 2.*xmu*e12  ! shear_stress is a scalar that is reused in MP version
	   else
	       sshr(m,n) = 2.*xmu*e12     ! sshr is an array that is referred to repeatedly in noMP version
	   end if
C	Unit 17 no longer defined. Probably do not want to write sshr, but may want to in some cases ++ TET fix
C       write (17,73) n,m,sshr(n,m)  
C   73  format (1x,2i6,10f10.4)
	return
	end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
      SUBROUTINE STRN (xw2,xd2,x3,y1,y3,e12)
c       Equations from Chinnery (1963)   y1 par flt  y2 nor flt  y3 dwn
c       analytic soln pers. com. 1983 WDS; xlam=xmu  x1=x2=y2=0
c       tensor strain e12 on fault plane
c       x (src cell center)  y (obs pt)  in disloc coords x

      implicit real*8 (a-h,o-z)
      dimension sgn(4),fx1(4),fx3(4)
      parameter (cpi=1./(4.*3.1415926), c23=2./3.)
      data sgn /1.,-1.,-1.,1./

      fx1(1)=  +xw2              ! xw2,xd2 are half width, half height
      fx1(2)=  +xw2
      fx1(3)=  -xw2
      fx1(4)=  -xw2
      fx3(1)=x3+xd2
      fx3(2)=x3-xd2
      fx3(3)=x3+xd2
      fx3(4)=x3-xd2
      e12=0.

      do 20 i=1,4
      t=fx1(i)-y1
      q=fx3(i)-y3
      p=fx3(i)+y3
      tt=t*t
      qq=q*q
      pp=p*p
      s1=sqrt(tt+qq)
      s2=sqrt(tt+pp)
      s1q=s1+q
      s2p=s2+p

      f1 = 1./(s1*s1q) + 1./(s2*s2p)
      f2 = (s2/4.+q)/(s2*s2p*s2p) - (pp-qq)*(2.*s2+p)/(2.*s2**3*s2p*s2p)
      f4 = q/(s1*(s1+t)) + p/(s2*(s2+t))

      etmp = c23*t*(f1 + f2) + f4/2.
   20 e12  = e12 + sgn(i)*cpi*etmp
      return
      end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
      SUBROUTINE KOMPRS (nb,ibctmp,nf)

c       compress k matrices, remove bc rows, cols

      USE Green_sshr	! Make data in module "Green_sshr" available
      implicit real*8 (a-h,o-z)
      real*4 sa,sab
      parameter (NF_MAX=15003, NE2_MAX=30006)
      common /deriv/ tauvp(NF_MAX),ca(NF_MAX),cab(NF_MAX),dx(NF_MAX),
     . dz(NF_MAX),tauss(NF_MAX),sa(NF_MAX),sab(NF_MAX),
     . e10,xldis,vpl,xmu,ne,ne2,MMP
      dimension ibctmp(1)
c											nf = # of bc and non-bc cells
c											ne = # of non-bc cells
c											nb = # of bc cells
      do 45 n=1,nb
          i=ibctmp(n)
          do j=i,nf-1                      ! move rows up
              do k=1,nf
                  sshr(j,k)=sshr(j+1,k)
		      end do
          end do
c
          do k=1,nf                        ! fill in last row with zeros
              sshr(nf,k)=0.
		  end do
c
          do k=i,nf-1                      ! move columns left
              do j=1,nf
                  sshr(j,k)=sshr(j,k+1)
		      end do
		  end do
c
          do j=1,nf                        ! fill in last column with zeros
              sshr(j,nf)=0.
		  end do
c
          do j=n,nb                        ! use indices of new matrix
              ibctmp(j)=ibctmp(j)-1
          end do
   45 continue
      return
      end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
      SUBROUTINE PATCH4 (tv,factinit)

c       called from main, but not called in current way of assigning
c			constitutive parameters - may want to use in the future, 
c			perhps in addition to PATCH5 that is currently used, in 
c			order to assign local variations in paraamenter to create
c			microearthquakes.

c       14 august 92 long cylindrical patch - WDS
c       dble gaussian lft, rt ends; z gauss. center

c       A-B inverted gaussians; min=abmn, max=abmx to inf

c       tau, vel, u, vp always pos
c       initial tv set in subr init
c       A linear w/ depth only

c       ca    A
c       cab   A-B, d tau ss/d vel, <0 for instab,  >0 visc
c             ca(i),cab(i), ... corresponds to tv(i),tv(ne+i)
C     Dimensioning changes as change model:
C        nf (or ne?) is the dimension of many arrays: currently is 15003
C
C        2nf (or 2ne?) is the dimension of a few arrays: currently is 30006
C          but none of these seem to exist in prkf2s.f, only in prkf2.f
C
C        nblk is the dimension of several arrays: currently is 40
C
C        The number of neighbor cells (+1) that a cell can have was set 
C        to 10 by WDS, but I changed it to 25 to accomodate a 9x reduction
C        at a boundary on up to two sides of a cell, plus 3x on the other
C        two sides and one for the number of entries in the array. This is
C        the second dimension of nabor( , ).

      implicit real*8 (a-h,o-z)
      real*4 sa,sab
      parameter (NF_MAX=15003, NE2_MAX=30006)
      common /deriv/ tauvp(NF_MAX),ca(NF_MAX),cab(NF_MAX),dx(NF_MAX),
     . dz(NF_MAX),tauss(NF_MAX),sa(NF_MAX),sab(NF_MAX),
     . e10,xldis,vpl,xmu,ne,ne2,MMP
      common /pat/ xc(NF_MAX),zc(NF_MAX),vlim,npat,nb,ibc(20),
     & iec(NF_MAX),iecinv(NF_MAX)
      dimension tv(1), pxo(5),pzo(5),ax(5),az(5),abmn(5),abmx(5)
      parameter (pi=3.1415926)

      do 8 n=1,ne
          ca(n)  = 0.
    8     cab(n) = 0.
      do 10 n=1,ne2     ! tau, vel
   10     tv(n) = 0.
      read (2,*)
      do 14 k=1,2
          read  (2,4)  pxo(k),pzo(k),xmn,xmx,zmn,zmx,ax(k),az(k)
          read  (2,4)  ac,arat,abmn(k),abmx(k),xldis

          write (6,3)
          write (6,7)  pxo(k),pzo(k),xmn,xmx,zmn,zmx,ax(k),az(k)
          write (6,5)
          write (6,6)  ac,arat,abmn(k),abmx(k),xldis

          write (7,3)
          write (7,7)  pxo(k),pzo(k),xmn,xmx,zmn,zmx,ax(k),az(k)
          write (7,5)
          write (7,6)  ac,arat,abmn(k),abmx(k),xldis

          write (8,3)
          write (8,7)  pxo(k),pzo(k),xmn,xmx,zmn,zmx,ax(k),az(k)
          write (8,5)
          write (8,6)  ac,arat,abmn(k),abmx(k),xldis
   14 continue

      ampl   = abmx(1)-abmn(1)
      do 20 n=1,ne
          i      = iec(n)                ! i cell no., n tv indx
          ca(n)  = ac     + arat*zc(i)
          cab(n) = abmx(1)
          rz     = (zc(i)-pzo(1))/az(1)
          erz    = exp(-rz*rz)
          rx     = (xc(i)-pxo(1))/ax(1)
          if (xc(i) .le. pxo(1)) go to 16
          if (xc(i) .lt. pxo(2)) go to 18

c         lft or rt dbl gauss. patch end
          rx     = (xc(i)-pxo(2))/ax(1)
   16     cab(n) = cab(n) - ampl*erz*exp(-rx*rx)
          go to 20

c         sngl gauus. z dir patch center
   18     cab(n) = cab(n) - ampl*erz
   20 continue

      CALL INIT (tv,irestart,factinit)

      write (7,42)           ! ++ ?
      caneg=1.
      do 40 n=1,ne
          i=iec(n)
          if (ca(n) .le. 0.) caneg=0.
          write (7,44) n,i,ca(n),cab(n),tv(n),tv(ne+n)
   40 continue

      if (caneg .ne. 0.) return
      write (6,41)
      write (7,41)
   41 format (/' err subr patch, bad input data: some fault law a=0')
      stop

    4 format (15f5.4)
    6 format (1x,10f8.2)
    7 format (1x,6f8.2,2f8.4)
   44 format (1x,2i5, 2f8.2, f8.2,1pe9.1)
    3 format (/'      pxo     pzo    xmn     xmx     zmn     zmx',
     &          '     ax      az')
    5 format (/'      ac     arat    abmn    abmx   xldis')
   42 format (/'     n  cell     A     A-B    tau(0)  vel(0)')
      end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
      SUBROUTINE PATCH5 (tv,iout6,iout7,irestart,factinit,a_and_amb,
     &                                                           nminc)

c       called from main

c       THIS IS THE VERSION OF SUBR PATCH THAT USES THE Tse-Rice (1986)
c       METHOD OF ASSIGNING A AND A-B.  From TET 1/93

c       Stresses in bars for Parkfield Version here
c       fault normal stress is positive this subr.

c       Uses the A-B vs. depth relation from Blanpied et al. (1991)
c         instead of the patch model of Stuart or the Tse & Rice depth data.
c       Used Geotherm from Sass, Galanis, and Lachenbruch (Eos, 1992) 
c		  poster, closer to "Coast Range Average" than "Parkfield" from 
c         talking to Sass and Lachenbruch. I picked
c         a geotherm that gave T of 10 deg. at surface and 760 deg. at
c         30.0375 km and has about the right amount of curvature.

c       In this Parkfield version a-b in friction units is
c       converted to A-B (bars) at all Temps.

c       T       - temperature, function of depth zc, + downward
c       rsc     - dist along trace, lft to rt, same as xc() here
c       iec(i)  - actual cell number for ca(i) position
c       a,b,c   - coefficients in temperature vs depth eqn
c       ca      - A
c       cab     - A-B, d tau ss/d vel, <0 for instab,  >0 visc
c       ca(i),  - corresponds to tvrn(i),tvrn(ne+i),...,tvrn(ne3+i)

c       ----- input vars -----

c       rsco      - right strike coord of band center, =xc()
c       xl,xr     - distance beyond lft,rt side of band where A-B=0
c       wx        - half width patch band
c       frict0    - not used

c       ac        - linear outside band: a(zc=0)
c       amx       -                      da/dz
c       abc       -                      (a-b)(zc=0)
c       abmx      -                      d(a-b)/dz
c       xldis     - L of t&r
c       depthfact - scale factor applied to L&S z(T); = .73

      implicit real*8 (a-h,o-z)
      real*4 sa,sab
      INTEGER :: a_and_amb      ! Flag to say if a and amb are output for plotting check
      parameter (NF_MAX=15003, NE2_MAX=30006)
      common /deriv/ tauvp(NF_MAX),ca(NF_MAX),cab(NF_MAX),dx(NF_MAX),
     . dz(NF_MAX),tauss(NF_MAX),sa(NF_MAX),sab(NF_MAX),
     . e10,xldis,vpl,xmu,ne,ne2,MMP
      common /pat/ xc(NF_MAX),zc(NF_MAX),vlim,npat,nb,ibc(20),
     & iec(NF_MAX),iecinv(NF_MAX)
      dimension tv(1),Temp(NF_MAX)

      do 18 n=1,ne
          ca(n)  = 0.
   18     cab(n) = 0.
      do 16 n=1,ne2     ! tau, vel
   16     tv(n)  = 0.
c  ------  Allow only 1 patch for now  ---------  ++ fix if needed
      write (6,*)'from PATCH5: npat=',npat
      if (npat .eq. 1) go to 10
      write (6,12) npat
   12 format (/' err: subr patch5 - only npat=1 allowed this subr',i5)
      stop

   10 read (2,*)
      read (2,*) rsco,xl,xr,wx,frict0
      read (2,*) ac,amx,abc,abmx,xldis,depthfact

      if (iout6 .eq. 1) then
          write (6,3)
          write (6,7) rsco,xr,xl,wx,frict0
          write (6,5)
          write (6,6) ac,amx,abc,abmx,xldis,depthfact
      end if

      if (iout7 .ne. 0) then
          write (7,3)
          write (7,7) rsco,xr,xl,wx,frict0
          write (7,5)
          write (7,6) ac,amx,abc,abmx,xldis,depthfact
      end if

      write (8,3)
      write (8,7) rsco,xr,xl,wx,frict0
      write (8,5)
      write (8,6) ac,amx,abc,abmx,xldis,depthfact

c       z(T):  z = a * T**2 + b * T  + c   (+z down)
c       Multiplying by depthfact is to make the depth of everything more
c       appropriate for observed eq depths. T's normally" reached at
c       a depth of 10 km will be reached at depth of 10*depthfact km

      a =  5.0e-6
      b =  3.62e-2
      c = -0.3625

      a = depthfact * a
      b = depthfact * b
      c = depthfact * c

      z150 = a*150.*150. + b*150. + c         ! z150 is depth where T=150
      z300 = a*300.*300. + b*300. + c
      z350 = a*350.*350. + b*350. + c
      zbig = 30.                              ! zbig > z(T=350)

      T000   = 1./(2.*a) * (-b + sqrt( b*b-4*a*c        ) ) ! for check
      T150   = 1./(2.*a) * (-b + sqrt( b*b-4*a*(c-z150) ) )
      T300   = 1./(2.*a) * (-b + sqrt( b*b-4*a*(c-z300) ) )
      T350   = 1./(2.*a) * (-b + sqrt( b*b-4*a*(c-z350) ) )
      Tbig   = 1./(2.*a) * (-b + sqrt( b*b-4*a*(c-zbig) ) )

c ++  Tbig (zc=30) = 1005 but should = 760 according to tet ++ TET fix - check to see if there is a problem

      ca000  = 0.0
      cab000 = 0.0    ! ++ TET fix? - unlike Blanpied et al. but perhaps should be >0 as they have it

      ca150  = 0.012*18.*z150*10.
      cab150 = 0.004 - T150 * 5.33333333e-5
      cab150 = cab150*18.*z150*10.

      ca300  = 0.012*18.*z300*10.
      cab300 = -0.004
      cab300 = cab300*18.*z300*10.

      ca350  = 0.012*18.*z350*10.
      cab350 = -0.100 + T350 * 3.2e-4
      cab350 = cab350*18.*z350*10.

      cabig  = 18.*z350* (0.012 + (Tbig - 350.) * 0.000089) ! ++ ?
      cabig  = cabig * 10.
      cabbig = cabig

      z000 = 0.           ! ++  fix - z000 is unused

      if (iout6 .ne. 1) go to 30
      write (6,28)
      write (6,14) T000,z000,ca000,cab000
      write (6,14) T150,z150,ca150,cab150
      write (6,14) T300,z300,ca300,cab300
      write (6,14) T350,z350,ca350,cab350
      write (6,14) Tbig,zbig,cabig,cabbig

   30 if (iout7 .ne. 1) go to 32
      write (7,28)
      write (7,14) T000,z000,ca000,cab000
      write (7,14) T150,z150,ca150,cab150
      write (7,14) T300,z300,ca300,cab300
      write (7,14) T350,z350,ca350,cab350
      write (7,14) Tbig,zbig,cabig,cabbig

      write (8,28)
      write (8,14) T000,z000,ca000,cab000
      write (8,14) T150,z150,ca150,cab150
      write (8,14) T300,z300,ca300,cab300
      write (8,14) T350,z350,ca350,cab350
      write (8,14) Tbig,zbig,cabig,cabbig

   14 format (1x,4f9.2)
   28 format (/' Piecewise linear a, a-b vs. temp from lab'/
     .         ' T, A, and A-B nonlinear w/ depth'/
     . '       T        z        A       A-B')

c     ----- Cell loop: Assign A-B, A vs. depth, centerline of patch -----

c       assign A-B, A all cells; overwrite cells outside band later

c       Use temperature to determine the constitutive parameters
c       a and a-b (in friction units) at and below the temperature
c       at which b=0. For T's greater than that use A=A-B,
c       in stress units (MPa) - (Bars used in this Parkfield version).

c       The temperature at which b=0 is set to 350. degrees C.

c       T(z) from binomial theorem, = inverse of quadratic fit

c     ----- in do 20 first calc A (2 segs), then A-B (4 segs) -----

   32 cabmin = 1.e6
      do 20 n=1,ne          ! n is sequence # within iec(), pde list tv
          i=iec(n)          ! i is actual cell number, nf total
          Temp(n) = 1./(2.*a) * (-b + sqrt( b*b-4*a*(c-zc(i)) ) )
          if (Temp(n) .le. 350.) then
              ca(n) = 0.012                   ! this is a (=A/signorm)
              ca(n) = 0.012*18.*zc(n)*10.     ! this is A (a*signorm), bars
          end if

c         A (=a*signorm) to match a*18 at z350 and then incr. w/ T

          if (Temp(n) .gt. 350.) then
              ca(n) = 18.*z350* (0.012 + (Temp(n) - 350.) * 0.000089)
              ca(n) = ca(n) * 10.             ! MPa to bars, no dep. on signorm
          end if

c         Piecewise linear A-B law:          T .le. 150  deg C
c                                   150 .lt. T .le. 300
c                                   300 .lt. T .le. 350
c                                   350 .lt. T           :B = 0, A-B = A

c       Next lines are fit to Blanpied et al. (1991, Fig 2)
c ++    Mike Blanpied says inconsistency, should be a=0.014 ++ TET fix -check what I use and what it should be
c ++    TET uses A-B=0 at z=0, unlike Blanpied et al. (1991) A-B>0 at z=0, based on TET data
c ++     However, making A-B>0 should help have creep at Earth's surface, so >0  might be right!

          if (Temp(n) .le. 150.) then
              cab(n) = 0.004 - Temp(n) * 5.33333333e-5  ! this is a-b (mu units)
              cab(n) = cab(n)*18.*zc(n)*10.             ! this is A-B (bars)
          end if
          if (Temp(n) .gt. 150. .and. Temp(n) .le. 300.) then
              cab(n) = -0.004
              cab(n) = cab(n)*18.*zc(n)*10.
          end if

          if (Temp(n) .gt. 300.  .and. Temp(n) .le. 350.) then
              cab(n) = -0.100 + Temp(n) * 3.2e-4
              cab(n) = cab(n)*18.*zc(n)*10.
          end if

          if (Temp(n) .gt. 350.) cab(n)=ca(n)

          if (cab(n) .ge. cabmin) go to 20
          cabmin = cab(n)                     ! min val of A-B
          nminc  = iec(n)                     ! cell no. of min A-B
   20 continue

      if (iout6 .eq. 1) write (6,22) cabmin,zc(nminc),nminc
      if (iout7 .ne. 0) write (7,22) cabmin,zc(nminc),nminc
                        write (8,22) cabmin,zc(nminc),nminc ! ++ TET fix - am I writing too much into unit 8? - it is a summary
   22 format (/' (A-B)min = ',f6.2, ' at depth of', f6.2,
     &  ' (= depth of cell ',i4,')')

c     ----- end center patch profile A, A-B -----

c     ----- start gaussian horizontal variation outside band -----

c     The scheme is that there is a Gaussian change in the strike direction
c     from the depth-determined value inside band to the maximum value of
c     A-B if the x point is outside the band where only the depth
c     determined value is used. ++ ?

c     Inside a band with half width = wx, only the depth determined value
c     is used.

c     Outside that band the maximum value of A-B that the Gaussian decays 
c     to is a linear function of depth.

c     Gaussian is A * exp{-(delx/w)**2} , where del x is distance from peak,
c     or in our case the beginning of the decline, and w is the distance
c     it takes to fall to 1/e of its value, A being the peak value.

c     In the band at the depth where A-B is the minimum value, the Gaussian
c     is set up so that it will fall to zero in a distance of xr on right
c     and  xl on the left.  It keeps falling to the depth determined value
c     of A-B.  At all other depths, the Gaussian decay has the same
c     e-folding length scale as at the depth of the most negative A-B.

c     This means that the zero crossing and the rate of decay in
c     absolute A-B units will differ somewhat, since the starting point
c     of the decay will be closer to zero (or positive) and the total
c     amplitude of the decay will depend on the difference between
c     the depth determined value in the veloctiy weakening patch and the
c     depth determined value (abmax) outside of that.

c     The scheme for A is similar.  The same e-folding length is
c     used for A as for A-B. The values in the band and the maximum
c     values outside the band can be determined in somewhat different
c     ways.

c     xxr   - e-decay distance to make A-B=0 at xl, xr
c     abmax - A-B outside band at zc of A-B min inside

c       Calc e-folding decay distance of Gaussian
c       Determine max  value of A-B at depth of min value of A-B ++ ?

      abmax = abc + abmx*18.*zc(nminc)
      abmax = abmax * 10.                   ! convert to bars
      abmm  = abmax - cabmin                ! amplitude of variation

      xxl = xl / sqrt (log(abmm/abmax))    ! same for left side
      xxr = xr / sqrt (log(abmm/abmax))    ! right side decay dist

      write (8,24) abmax,abmm,xxl,xxr
   24 format (/'   abmax     abmm    xxl      xxr'/ 4f9.2)

c       Now use these special xxr, xxl for general cells outside band

      do 25 n=1,ne
          i = iec(n)
          delrsc = rsco-xc(i)                  ! signed dist from band center
          if (abs(delrsc) .le. wx) go to 25    ! skip inside band
c         ---- Set up variation of a-b from depth det. value to outside ----
c         abmx = 0.025 Serpentine
          abmax  = abc + abmx*18.*zc(i)        ! linear var. w/ depth outside
          abmax  = abmax*10.                   ! convert to bars
          abmm   = abmax - cab(n)              ! cab = val inside band
c         ---- Set up variation of a same way ----
          amax   = ac + amx*18.*zc(i)          ! linear var. w/ depth
          amax   = amax*10.                    ! bars
          amm    = amax - ca(n)
          if (delrsc .gt. wx) then             ! rt side Gaussian
              xrt    = (delrsc-wx)/xxr
              cab(n) = abmax - abmm * exp(-xrt*xrt)
              ca(n)  = amax  - amm  * exp(-xrt*xrt) ! A decay same as A-B
          else                                  ! lft side Gaussian
              xlf    = -(wx+delrsc)/xxl             ! sum here = - diff of abs values
              cab(n) = abmax - abmm * exp(-xlf*xlf)
              ca(n)  = amax  - amm  * exp(-xlf*xlf)
          end if
   25 continue
c     ----- end horiz. var. loop -----

c     This section is writing out a and a-b into a binary file
c       so it can be plotted in color to check that it looks OK
      do i=1,ne
	      sa(i)=0.0                     ! single precision output
          sab(i)=0.0
      end do
      do i=1,ne
          sa(i)=sngl( ca(i) )
          sab(i)=sngl( cab(i) )
      end do
      ndum1=1
      dum2=0.0                          ! double precision output
      ndum3=0
	  if(a_and_amb .eq. 1) then      
          write (20) ndum1,dum2,ndum3
          write (20) (sa(i),i=1,ne)         ! single precision output
          write (21) ndum1,dum2,ndum3
          write (21) (sab(i),i=1,ne)        ! single precision output
          CLOSE (UNIT = 20)
          CLOSE (UNIT = 21)
	  end if
c     end of writing out a and a-b section

      CALL INIT (tv,irestart,factinit)

      if (iout7 .ne. 0) write (7,42)
      caneg = 1.
      do 40 n=1,ne
          i = iec(n)
          if (ca(n) .le. 0.) caneg=0.
          write (7,44) n,i,ca(n),cab(n),tv(n),tv(ne+n)
   44     format (1x,2i10, 2f8.2, f8.2,1pe9.1)
   40 continue

      if (caneg .ne. 0.) return
      write (6,41)
      if (iout7 .ne. 0) write (7,41)
   41 format (/' err subr patch, bad input data: some fault law a=0')
      stop

    4 format (15f5.4)
    6 format (1x,10f10.4)
    7 format (1x,6f10.4)
    3 format (/' TET vals A and A-B from Blanpied et al.',
     & ' (1991) a-b values and Sass et al (1992) geotherm',//
     & '      rsco       xr       xl        wx       frict0')
    5 format (/'        ac        amx       abc      abmx    xldis
     &   depthfact')
   42 format (/'          n       cell     A     A-B    tau(0)  vel(0)')
      end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
      SUBROUTINE INIT (tv,irestart,factinit)

c ++    add smoothing to ca, cab
c       initial conditions on tau, sig, vel
c       coseis zone approx immed post instab v=vlim
c       deep viscous zone approx linear v to bc vpl  ++fix
c       coseis-deep boundary is approx where A-B=0, get from subr patch
c       tau from these values of v

      implicit real*8 (a-h,o-z)                 ! TET Mod
      real*4 sa,sab
      parameter (NF_MAX=15003, NE2_MAX=30006)
      common /deriv/ tauvp(NF_MAX),ca(NF_MAX),cab(NF_MAX),dx(NF_MAX),
     . dz(NF_MAX),tauss(NF_MAX),sa(NF_MAX),sab(NF_MAX),
     . e10,xldis,vpl,xmu,ne,ne2,MMP
      common /pat/ xc(NF_MAX),zc(NF_MAX),vlim,npat,nb,ibc(20),
     & iec(NF_MAX),iecinv(NF_MAX)
      dimension tv(1)

      e10=exp(-10.)
      if(irestart.eq.0)then
         do 10 n=1,ne
             i=iec(n)
             nv=ne+n
    6        tv(nv) = vpl                         ! viscous area outside patch
             if (cab(n) .gt. 0.) go to 8
             tv(nv) = factinit * vlim             ! brittle area inside patch
    8        tv(n)  = -cab(n)*log(vpl/tv(nv)+e10) ! tauss(vel init)
   10    continue
      else       ! if is restarting from end of earlier run, read in tv
	     read (12) (tv(i),i=1,ne2)
      end if
c ++  CALL SMOOTH (tv,ne)         ! smooth initial tv ++ fix - not called now
      return
      end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
      SUBROUTINE SMOOTH (tv,ne)   ! This not called at present - needs finishing

c       smooth initial tau, vel field on cell block 1 only  9/92
c ++    finish, check
      implicit real*8 (a-h,o-z)
      parameter (NF_MAX=15003, NE2_MAX=30006)
      dimension tv(1),t(NE2_MAX)

      mapx=39                         ! ++ pass arg
      mapz=15
      ne2=2*ne
      do 8 n=1,ne2                    ! need for map edges of copy back
    8     t(n)=tv(n)

      do 10 i=2,mapz-1                ! row
          do 10 j=2,mapx-1            ! col
              do 10 k=1,2             ! first tau, then vel
                  m = (i-1)*mapx+j + (k-1)+ne
                  mup = m-39
                  mdn = m+39
                  t(m) = ( tv(m-1) + tv(m+1) + tv(mup) + tv(mdn) )/4.  ! tau
c                 t(m+ne) =                          ! vel - ++ this not completed
   10 continue
      do 20 n=1,ne2                   ! copy back smoothed tau, vel
   20     tv(n)=t(n)
      return
      end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
      SUBROUTINE CREEP (nt,t,tv,xc,ne,m,ncr,ncrp,crpnam)

c       print cell vel for assigned trace (or other) cells

      implicit real*8 (a-h,o-z)
      character*5  crpnam(1)
      dimension    crp(10),tv(1),xc(1),ncrp(1)

      if (m .eq. 0) go to 10
      write (9,4) (ncrp(i),i=1,ncr)
      write (9,6) (xc(ncrp(i)),i=1,ncr)
      write (9,5) (crpnam(i),i=1,ncr)
      write (9,9)
    4 format (/6x,'   cell',15i5)
    6 format ( 6x,'   x   ',15f5.1)
    5 format (/6x,'   name',15a5)
    9 format (/' t step  year  ----------- creep rates, mm/yr --------')
      return
   10 do 12 i=1,ncr
          k = ncrp(i)
          crp(i) = tv(k+ne)
          if (crp(i) .gt. 999.) crp(i)=999. ! avoid format overrun
   12 continue
      write (9,8) nt,t,(crp(i),i=1,ncr)
    8 format (1x,i5,f7.2,15f5.1)
      return
      end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
      SUBROUTINE WRSTEP (t,tlast,dt,nyes)

c       return nyes=1 if time since last write tlast ge dt years

      implicit real*8 (a-h,o-z)

      nyes=0
      if (t-tlast .lt. dt) return
      nyes=1
      tlast=t
      return
	  end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
      SUBROUTINE RKQC (y,dydt,ny,dtold,eps,yscal,dtdid,dtnxt)

c       adaptive step size control for rk4, from press et al
c       fifth order runge kutta  dimension all by no. of pde's
c++     errcon=(4/safety)**(1/pgrow)

      implicit real*8 (a-h,o-z)
      parameter (NF_MAX=15003, NE2_MAX=30006)
      parameter (pgrow=-0.20,pshrnk=-0.25,safety=0.9,errcon=6.e-4)
      dimension y(1),dydt(1),yscal(1)
      dimension ytmp(NE2_MAX),ysav(NE2_MAX),dysav(NE2_MAX)

      do 11 i=1,ny                              ! for initial step
          ysav(i)=y(i)
   11     dysav(i)=dydt(i)
      dt=dtold
    1 dtt=0.5*dt
      CALL RK4 (ysav,dysav,ny,dtt,ytmp)         ! first half step
      CALL DERIVS (ytmp,dydt)
      CALL RK4 (ytmp,dydt,ny,dtt,y)             ! second half step
      CALL RK4 (ysav,dysav,ny,dt,ytmp)          ! the large step
      errmax=0.
      do 12 i=1,ny
          ytmp(i)=y(i)-ytmp(i)                      ! err est
   12     errmax=max(errmax,abs(ytmp(i)/yscal(i)))

      errmax=errmax/eps                         ! scale to required tol
      if (errmax .gt. 1.) then                  ! trunc err too large
          dt=safety*dt*(errmax**pshrnk)         ! reduce step size
          go to 1                               ! try smaller dt
      else
          dtdid=dt                              ! this dt ok
          if (errmax .gt. errcon) then
            dtnxt=safety*dt*(errmax**pgrow)
          else
            dtnxt=4.*dt
          endif
      endif
      do 13 i=1,ny
   13     y(i)=y(i)+ytmp(i)/15.                 ! get fifth order est
      return
      end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
      SUBROUTINE RK4 (y,dydt,n,h,yout)

c       fourth order runge-kutta from press et al.
c       dimension all by 2*nf

      implicit real*8 (a-h,o-z)
      parameter (NF_MAX=15003, NE2_MAX=30006)
      dimension y(1),dydt(1),yout(1)
      dimension yt(NE2_MAX),dyt(NE2_MAX),dym(NE2_MAX)

      hh=h*0.5
      h6=h/6.
      do 11 i=1,n                               ! first step
   11     yt(i)=y(i)+hh*dydt(i)

      CALL DERIVS (yt,dyt)                      ! second step
      do 12 i=1,n
   12     yt(i)=y(i)+hh*dyt(i)

      CALL DERIVS (yt,dym)                      ! third step
      do 13 i=1,n
          yt(i)=y(i)+h*dym(i)
   13     dym(i)=dyt(i)+dym(i)

      CALL DERIVS (yt,dyt)                      ! fourth step
      do 14 i=1,n                               ! accumulate incrmnts
   14     yout(i)=y(i)+h6*(dydt(i)+dyt(i)+2.*dym(i))     ! w/ weights
      return
      end
C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************


C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************

CPrologue for Functions or Subroutines

C!---------------------------------------------------------------------------------------------------
C!  BOP
C!
C!  ROUTINE:  <Function name>
C!
C!  INTERFACE:
C	function <name> (
C		<arg1>,			! <description>
C		<arg2>,			! <description>
C		   .
C		   .
C		<argn>)			! <description>
C
Cc  ++ TET fix - The commented lines below are the ones that will eventually be converted to Fortran 90/95, but not yet done.
C!
C!  USES:
Cc	use <module>    ++ TET fix - Not yet converted to Fortran 90/95 use of modules
C!
C!  RETURN VALUE:
Cc	implicit none   ! ++ TET fix - Not yet converted to using "implicit none" and explicitly typing all variables
Cc	<type> :: <name>			! Description
C!
C!  PARAMETERS:
Cc	<type, intent> :: <parameter>		! Description
C!
C!  NECESSARY PRE-CONDITIONS:
C!	e.g., specific values being set in Module global variables; a specific object
C!	having been created; etc.
C!
C!  SIDE EFFECTS:
C!	e.g., modifying global variables; creating objects that persist beyond 
C!	scope of the routine/function.
C!
C!  POST_CONDITIONS UPON EXIT
C!	e.g., any files being closed, objects being deleted, etc.
C!
C!  DESCRIPTION:  Describe the function of the routine and algorithm(s) used
C!			in the routine.  Include any applicable external references.
C!
C!  EOP
C!---------------------------------------------------------------------------------------------------


CPrologue for a Module

C!---------------------------------------------------------------------------------------------------
C!  BOP
C!
C!  MODULE:  <Module name>
C!
C!  USES:
C	use <module>
C!
C!  PUBLIC TYPES
C	implicit none
C	<type declaration>

C!  PUBLIC MEMBER FUNCTIONS
C!	<function>				!  Description
C!
C!  PUBLIC DATA MEMBERS
C	<type> ::  <variable>			!  Variable description

C!  DESCRIPTION:
C!	<Describe the function of the module>
C!
C!  EOP
C!-----------------------------------------------------------------------------------------------------


CPrologue for a Header File

C!---------------------------------------------------------------------------------------------------
C!  BOP
C!
C!  !INCLUDE:	  <Header file name>
C!
C!  !DEFINED PARAMETERS
C	<type>  ::  <parameter>		! parameter description
C!
C!  DESCRIPTION:
C!	<Describe the contents of the Header file>
C!
C!  EOP
C!----------------------------------------------------------------------------------------------------

C++ TET fix - add what routines are called by each routine and what routines call each routine