C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
!---------------------------------------------------------------------------------------------------
!  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*2, 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

!---------------------------------------------------------------------------------------------------
!  BOP
!
!  MODULE: cblock
!
      module cblock
!
!  USES:
C	 use cblock
!
!  PUBLIC TYPES
C	implicit none
!
!  PUBLIC DATA MEMBERS
      real*8, allocatable, dimension(:) :: sxc,szc,sdx,sdz,     ! these 's' variables are the key-sorted  arrays assigned 
     & sv, scab, sca, stauvp, stauss                            ! to the individual processes. the variable purpose, for example
                                                                ! stauss, is given by the base  variable definition, in this example
      integer, allocatable, dimension(:) :: skey1               ! by tauss. these are sub arrays, that is they are local to each process
                                                                ! and their dimension is the number of observations assigned to that processor

!  DESCRIPTION:
!	the multipole library requires that all the observation points (grid points) be sorted by key, and in the parallel implementation 
!       the observation points
!       of process n must have keys that are larger than those of process n-1 and smaller than process n+1. these 's' variables are the 
!       key-sorted  data associated with the observations
!  EOP
!-----------------------------------------------------------------------------------------------------
      end module cblock

!---------------------------------------------------------------------------------------------------
!  BOP
!
!  MODULE: masters
!
       module masters
!
!  USES:
C	use masters
!
!
!  PUBLIC DATA MEMBERS

      real*8, allocatable, dimension(:) :: rxc,rzc,rdx,rdz,   ! these are all of the key-sorted observation data. these arrays are local only
     & rv, rcab, rca, rtauvp, rtauss                          ! to the master process and have dimension ne (total number of grid points or elements)
                                                              ! the 'r' variables, as with the 's' variables, have purpose given by their base
                                                              ! counterpart. for example rxc is the array xc in key sorted order.   
      integer, allocatable, dimension(:) :: kindex, rkey1     ! kindex holds the index of each point in the base reference frame.
                                                              !
!  DESCRIPTION:
!       in the original program structure, the data associated with a grid point was held in a series of arrays, e.g., xc, zc, tauvp etc. because of the
!       requirements of the multipole library, the data in all of the arrays must be in keysorted order. the 'r' variables are the key-sorted versions
!       of the base variables
!  EOP
!-----------------------------------------------------------------------------------------------------
      end module

!---------------------------------------------------------------------------------------------------
!  BOP
!
!  MODULE: cderivs
!
         module cderivs
!  USES:
C	use cderivs
!
!
!  PUBLIC DATA MEMBERS
      real*8, allocatable, dimension(:) :: tauvp, ca, cab, 
     & dx, dz, tauss

      real*4, allocatable, dimension(:) :: sa, sab 

!  DESCRIPTION:
!	 these are variables that were in the common block derivs in the serial implementation. 
!         to save space they are declared as allocatable but only used and dimensioned by the 
!        master
!  EOP
!-----------------------------------------------------------------------------------------------------
      end module

!---------------------------------------------------------------------------------------------------
!  BOP
!
!  MODULE: cpat
!
        module cpat
!  USES:
C	use cpat
!
!
!  PUBLIC DATA MEMBERS
!
      real*8, allocatable, dimension(:) :: xc, zc  ! all of these variables are defined in the header to the main program
   
      integer, allocatable, dimension(:) :: iec, iecinv

!  DESCRIPTION:
!	 these are variables that were in the common block pat in the serial implementation. 
!         to save space they are declared as allocatable but only used and dimensioned by the 
!        master
!
!  EOP
!-----------------------------------------------------------------------------------------------------
      end module

!-----------------------------------------------------------------------------------------------------

C***********************************************************************
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C***********************************************************************
!---------------------------------------------------------------------------------------------------
!  BOP
!
!  MAIN PROGRAM:  PARK
!
!  USES:
      USE Green_sshr		! Make data in module "Green_sshr" available
      use cblock                ! use the variables and data in cblock module
      use masters               ! use the variables and data in masters module
      use cderivs               ! use variables in cderivs
      use cpat
!
        

c	implicit none   ! Not yet converted to using "implicit none" and explicitly typing all variables
      implicit real*8 (a-h,o-z)

      include 'mpif.h'


      INTEGER :: a_and_amb      ! Flag to say if a and amb are output for plotting check

      real*4 seps,t4
      parameter (NF_MAX=160003, NE2_MAX=320006)

      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

      real*8, allocatable, dimension (:) :: w, d, celar, u, uold, du,
     & tv, dtvdt, tvscal, tvold

      real*4, allocatable, dimension (:) ::  svln

      integer, allocatable, dimension (:) :: ibctmp
      integer, allocatable, dimension (:,:) :: nabor
 
      dimension out(100000,4,7)                        ! ++ increase if over 100,000 time steps
      dimension ibca(4),ibcb(4)

      dimension ncrp(10),nstep(10)

      dimension g1(10),g2(10) ! Increase if use > 10 random surf. pts

      integer numprocs, myid, ierr, master,  astat
      integer  status(MPI_STATUS_SIZE), istart, iend, anstype, sender
      integer finish, ib(4)

      dimension  b(6)

      real*8, allocatable, dimension (:) :: stv,sdtvdt,rtv
      integer, allocatable, dimension (:) ::  disps, rcounts

      common /deriv/ e10,xldis,vpl,xmu,ne,ne2,MMP

      common /pat/ vlim,npat,nb,ibc(20)

      common /newderiv/ p1(4), nnn, myid, umm(4)

      external derivs
!
!  DESCRIPTION:
c     Program PARK (main)
c
c     Earthquake simulation program
c
c     calls:
c      READ_FILE_NAME
c      cprogtime
c      READ_CELL_INFO
c      GREENS_FUNCTION
c      KOMPRS
c      PATCH5 (or PATCH4)
c      CREEP
c      DERIVS
c      RKQC
c      WRSTEP 
c
c     Link it with Ctm.o, upcase.o, numrec.o and appropriate derivs subroutine 
C		  file, e.g. derivs_slip.o or derivs_slowness.o or derivs_composite.o,
c         and with tree.o, print.o
c     The numrec.o file needs to be created from a numrec.f file that contains
c         two subroutines from Vol. 1 of the second edition of Numerical 
c         Recipes (NR) by Press et al. (Press, W. H., S. A. Teukolsky, 
c         W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran 77, 
c         The Art of Scientific Computing, Second Edition, Volume 1 of Fortran 
c         Numerical Recipes, Cambridge University Press, Cambridge, New York and 
c         Melbourne, 933 pages, 2001 - originally published 1992.) These two 
c         subroutines are rkqs and rkck, and they are found on pages 712-714. 
c         These routines are not included in this package of earthquake programs 
c         distributed freely on the web site because the NR progras are copyrighted. 
c         However, they can be obtained at minimal cost from the publisher. Two ways
c         to do this are 1) Go to the Numerical Recipies On-Line Software Store at 
c         http://www.nr.com, or 2) Get from Cambridge University Press by telephone
c         at (800) 872-7423 in North America, or by email to orders@cup.org for
c         North America or tread@cup.cam.ac.us for the rest of the world, or on the
c         web at http://www.cup.org for North America or http://w2ww.cup.cam.ac.uk 
c         for  the rest of the world.
c     We used the Fortran77 versions of the two Runge Kutta subroutines, rkqs and rkck,
c         namely the verison in Volume 1 of 2nd Ed. of NR rather than the version in 
c         Volume 2. Nevertheless, we are using it with other Fortran90 code and in fact
c         we make one modification to it that requeires it to be compiled as F90. Several
c         minor modifications have to be made to both of these subroutines, either due
c         to the need to make them run with MPI in parallel, or to be sure the
c         variable types are done consistently with the rest of the earthquake code.
c     These changes follow. Note that the lines below that are not indented (i.e. those 
c         with the comments that have the code starting in column 8), are to be inserted 
c         as code and uncommented. Those with a c! at the start are intended to be 
c         comment lines. 
c
c The changes that need to be made in rkqs are:
c         
c          Add the following two lines right after the "SUBROUTINE rkqs(y, dydx, . . )" line:
c      implicit real*8 (a-h,o-z)
c      include 'mpif.h'
c          In the 2 REAL type declaration lines, replace "REAL" with "REAL*8"
c          Remove the "PARAMETER (NMAX=50)" line.
c          In the second REAL type declaration line, replace "NMAX" with "n"
c		   (These two changes would be invalid in F77, but are allowed in F90 - this 
c           is the "automatic array" feature of F90)
c          Add a third REAL type declaration line:
c          	real*8 errscratch      ! scratch variable added for mpi call
c          Add the following lines right after the line "errmax=errmax/eps":
c!     parallel: check other processors to get the 
c!     maximum values of errmax
c!
c      call MPI_ALLREDUCE(errmax, errscratch,1,
c     + MPI_DOUBLE_PRECISION,MPI_MAX, MPI_COMM_WORLD,ierr) 
c
c      errmax=errscratch
c
c End of changes in rkqs.
c 
c The changes that need to be made in rkck are:
c
c          Add the following line right after the "SUBROUTINE rkck(y,dydy, . . )" line:
c      implicit real*8 (a-h,o-z)
c          In the 2 REAL type declaration lines, replace "REAL" with "REAL*8"
c          Remove the "PARAMETER (NMAX=50)" line.
c          In the second REAL type declaration line, replace "NMAX" with "n" (in 6 places)
c		   (These two changes would be invalid in F77, but are allowed in F90 - this 
c           is the "automatic array" feature of F90)
c
c End of changes in rkck.
c
c     Use Makefile provided to compile and link.
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 program uses 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
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 cash karp Runge Kutta programs from
c			the second edition of Numerical Recipies from Press et al. (subroutines rkqs and rkck).
c			This can be improved upon - the Runge Kutta approach is probably
c			the best one, but more efficient routines may 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 developed
c			by him and Michael Warren of Los Alamos National Laboratory. The
c           multipole library contained in the directory t17-7 is copyrighted
c           by John Salmon and Michael Warren.
c        5. The programs resulting from this collection of workers has
c           been cleaned up and additionally commented by Terry Tullis.
c        6. the program was converted to run in parallel using MPI by nick beeler of the 
c           United States Geological Survey and Brown University, 5/03 - 8/03 
c
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 outside 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 is the dimension of many arrays: Set by parameter NF_MAX
C        2nf is the dimension of a few arrays: Set by parameter NE2_MAX
C             This second parameter can probalby be only 2*ne long, not 2*nf long
C              but to be safe it is set to 2*nf at present - very little difference.
C             This explains the confusionin its name and size. ++ TET fix
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( , ).
!
!  EOP
!---------------------------------------------------------------------------------------------------


c
c     start mpi and mp library throgh c function
c
        call initializestuff()     ! this initializes MPI and the mplib
c
c     get processor ids and total number
c
      call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
      call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr ) 

      master=0

c
c     do all the io that sets up the problem on the master
c

9998  if (myid.eq.master)then

      print *, 'numprocs is : ', numprocs

c
c     allocate variables for input output and set up
c
      allocate (tauvp(NF_MAX), ca(NF_MAX), cab(NF_MAX), dx(NF_MAX),
     &  dz(NF_MAX), tauss(NF_MAX), sa(NF_MAX), sab(NF_MAX), 
     & xc(NF_MAX), zc(NF_MAX), iec(NF_MAX), iecinv(NF_MAX), STAT=astat) 

      allocate (w(NF_MAX), d(NF_MAX), celar(NF_MAX), u(NF_MAX), 
     & uold(NF_MAX), du(NF_MAX), svln(NF_MAX), ibctmp(NF_MAX), 
     & STAT=astat)

      allocate (tv(NE2_MAX), dtvdt(NE2_MAX), tvscal(NE2_MAX), 
     & tvold(NE2_MAX), STAT=astat)

      allocate (nabor(NF_MAX,25), STAT=astat)

      CALL READ_FILE_NAME(inputfile)

      print *,' number of steps between output writes for clocktime !! '
	  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=9
      read (1,'(a)') afile         ! 8 unit=20
      read (1,'(a)') ambfile       ! 9 unit=21
      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

      write (26,1310)  timedate
      write (27,1310)  timedate


      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 for the present MP version (MMP=1), 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.


          ALLOCATE (sshr(1,1), STAT=isshr_status)      ! degenerate array for MP

      if (isshr_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

                 dtv = dtv + shear_stress     ! shear_stress is a scalar, reused in MP version

         end do
         dtv = vpl*dtv
         tauvp(i) = tauvp(i) + dtv
   54 continue

      nbb=nb

C     These next 4 lines are to check if a and a-b are done correctly - a and amb written later - designed
c      to be in same format of binary file as the *.flt file so can use same ploting program to look at

      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,factinit)  
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
          t = 0.0
	  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

   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    allocate the sorted variables- these are the masters variables 
c

      allocate (rxc(ne), rzc(ne), rdx(ne), rdz(ne), rv(ne),
     & rcab(ne), rca(ne), rtauvp(ne), rtauss(ne), 
     & rtv(ne+ne), STAT=astat)

      allocate (rkey1(ne), kindex(ne), STAT=astat)
c
c   get keys and index for sorting from c routine
c
      do 42  i=1,ne
      rxc(i)=xc(i)
      rzc(i)=zc(i)
      rdx(i)=dx(i)
      rdz(i)=dz(i)
      rv(i)=tv(ne+i)
 42   continue

          CALL getkeys(ne,rxc(1),rzc(1),rdx(1),rdz(1), rv(1), 
     +             rkey1(1), kindex(1), umm(1)) ! Multipole version

c
c     figure out how many element per processor:
c       master gets the first few elements (iremain)
c       the slaves get idn each
c 
          if(numprocs.gt.1)then

      idn=int(floor( real((ne-1)/(numprocs-1))))
      itest1= ne-idn*(numprocs-1)
      idn=int(ceiling( real( (ne-2.*(numprocs-1))/numprocs+2)))
      itest2= ne-idn*(numprocs-1)
      iremain=max(itest1,itest2)
      idn=(ne-iremain)/(numprocs-1)

      else                ! special case 1 processor
         iremain=ne
         idn=0
      endif

      print *, 'number of elements = ', ne, ', number of slaves = ',
     + numprocs-1, ', slaves get ', idn,
     +               ', master gets = ', iremain

      write (6,41)     ! this writes line to the std output file

c
c     allocate the buffer variable and the sorted variables needed by master
c

      allocate (sxc(iremain), szc(iremain), sdx(iremain), 
     & sdz(iremain), sv(iremain),
     & scab(iremain), sca(iremain), stauvp(iremain), stauss(iremain), 
     & stv(iremain+iremain), sdtvdt(iremain+iremain), STAT=astat)

      allocate (skey1(iremain), STAT=astat)
c
c     populate the stresses and velocities of the key sorted
c     dependent variables 
c

      do 43 i=1,ne
         k=i+ne
         rtv(i)=tv(1+kindex(i))
         rtv(k)=tv(1+kindex(i)+ne)
         rcab(i)=cab(1+kindex(i))
         rca(i)=ca(1+kindex(i))
         rtauvp(i)=tauvp(1+kindex(i))
 43   continue

c
c     set values for master
c
      stv(1:iremain)=rtv(1:iremain)
      stv(iremain+1:iremain+iremain)=rtv(ne+1:ne+iremain)
      skey1(1:iremain)=rkey1(1:iremain)
      sxc(1:iremain)=rxc(1:iremain)
      szc(1:iremain)=rzc(1:iremain)
      sdx(1:iremain)=rdx(1:iremain)
      sdz(1:iremain)=rdz(1:iremain)
      sv(1:iremain)=rv(1:iremain)
      stauvp(1:iremain)=rtauvp(1:iremain)
      sca(1:iremain)=rca(1:iremain)
      scab(1:iremain)=rcab(1:iremain)
c
c     send the necessary arrays to the slaves
c
      if(numprocs.gt.1) then

      call MPI_BCAST(idn,1 , MPI_INTEGER, master,
     &                MPI_COMM_WORLD, ierr)           ! broadcast idn to all slaves

      call MPI_BCAST(iremain,1 , MPI_INTEGER, master,
     &                MPI_COMM_WORLD, ierr)

      call MPI_BCAST(umm,4 , MPI_DOUBLE_PRECISION, master,
     &                MPI_COMM_WORLD, ierr)           ! broadcast universal mins and maxs to all slaves


      do 47 i = 1,numprocs-1

        istart=iremain+1+(i-1)*idn
        iend=iremain+i*idn

        call MPI_SEND(rtv(istart:iend), idn, MPI_DOUBLE_PRECISION, i,
     &                 iend, MPI_COMM_WORLD, ierr)

        call MPI_SEND(rtv(istart+ne:iend+ne), idn, 
     &   MPI_DOUBLE_PRECISION, i, iend, MPI_COMM_WORLD, ierr)

        call MPI_SEND(rkey1(istart:iend), idn, MPI_INTEGER, i,
     &                 iend, MPI_COMM_WORLD, ierr)

        call MPI_SEND(rxc(istart:iend), idn, MPI_DOUBLE_PRECISION, i,
     &                 iend, MPI_COMM_WORLD, ierr)

        call MPI_SEND(rzc(istart:iend), idn, MPI_DOUBLE_PRECISION, i,
     &                 iend, MPI_COMM_WORLD, ierr)

        call MPI_SEND(rdx(istart:iend), idn, MPI_DOUBLE_PRECISION, i,
     &                 iend, MPI_COMM_WORLD, ierr)

        call MPI_SEND(rdz(istart:iend), idn, MPI_DOUBLE_PRECISION, i,
     &                 iend, MPI_COMM_WORLD, ierr)

        call MPI_SEND(rv(istart:iend), idn, MPI_DOUBLE_PRECISION, i,
     &                 iend, MPI_COMM_WORLD, ierr)

        call MPI_SEND(rtauvp(istart:iend), idn, MPI_DOUBLE_PRECISION, i,
     &                 iend, MPI_COMM_WORLD, ierr)

        call MPI_SEND(rca(istart:iend), idn, MPI_DOUBLE_PRECISION, i,
     &                 iend, MPI_COMM_WORLD, ierr)

        call MPI_SEND(rcab(istart:iend), idn, MPI_DOUBLE_PRECISION, i,
     &                 iend, MPI_COMM_WORLD, ierr)


 47   continue

c
c   send array b containing constants to slaves
c
      b(1)=e10
      b(2)=xldis
      b(3)=vpl
      b(4)=xmu
      b(5)=dt
      b(6)=eps

      call MPI_BCAST(b,6 , MPI_DOUBLE_PRECISION, master,
     &                MPI_COMM_WORLD, ierr)
c
c     send array ib containing time step numbers
c
      ib(1)=nt
      ib(2)=ninit
      ib(3)=ndt

      call MPI_BCAST(ib,3 , MPI_INTEGER, master,
     &                MPI_COMM_WORLD, ierr)

      endif

         nnn=iremain

c   this is the end of master's set up work that started at line label 9998

      else          ! this is the slaves set-up work

c
c    allocate the sorted variables- these are the slaves variables 
c     first get length (idn) from master

      call MPI_BCAST(idn,1 , MPI_INTEGER, master,
     &                MPI_COMM_WORLD, ierr)

      call MPI_BCAST(iremain,1 , MPI_INTEGER, master,
     &                MPI_COMM_WORLD, ierr)


      call MPI_BCAST(umm,4 , MPI_DOUBLE_PRECISION, master,
     &                MPI_COMM_WORLD, ierr)          

      allocate (sxc(idn), szc(idn), sdx(idn), sdz(idn), sv(idn),
     & scab(idn), sca(idn), stauvp(idn), stauss(idn), 
     & stv(idn+idn), sdtvdt(idn+idn), tvscal(idn+idn),  STAT=astat)

      allocate (skey1(idn), STAT=astat)
c
c     allocate the degenerate arrays unused by slaves
c
      allocate (tauvp(1), ca(1), cab(1), dx(1),
     &  dz(1), tauss(1), sa(1), sab(1), 
     & xc(1), zc(1), iec(1), iecinv(1), STAT=astat) 

      allocate (w(1), d(1), celar(1), u(1), 
     & uold(1), du(1), svln(1), ibctmp(1), 
     & STAT=astat)

      allocate (tv(1), dtvdt(1),
     & tvold(1), STAT=astat)

      allocate (nabor(1,1), STAT=astat)

c
c     slaves receive arrays from the master
c
         call MPI_RECV(stv(1:idn), idn, MPI_DOUBLE_PRECISION, master,
     &                  MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

         call MPI_RECV(stv(idn+1:idn+idn), idn,MPI_DOUBLE_PRECISION, 
     &     master, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

         call MPI_RECV(skey1, idn, MPI_INTEGER, master,
     &                  MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

c         print 49,  (stv(i),stv(i+idn), skey1(i),i=1,2)

         call MPI_RECV(sxc, idn, MPI_DOUBLE_PRECISION, master,
     &                  MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

         call MPI_RECV(szc, idn, MPI_DOUBLE_PRECISION, master,
     &                  MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

         call MPI_RECV(sdx, idn, MPI_DOUBLE_PRECISION, master,
     &                  MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

         call MPI_RECV(sdz, idn, MPI_DOUBLE_PRECISION, master,
     &                  MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

         call MPI_RECV(sv, idn, MPI_DOUBLE_PRECISION, master,
     &                  MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

         call MPI_RECV(stauvp, idn, MPI_DOUBLE_PRECISION, master,
     &                  MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

         call MPI_RECV(sca, idn, MPI_DOUBLE_PRECISION, master,
     &                  MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

         call MPI_RECV(scab, idn, MPI_DOUBLE_PRECISION, master,
     &                  MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

            finish = status(MPI_TAG)     ! finish is the master element number of the last point in the slave's array 
c
c        slaves receive array b from master, b contains values of 
c        e10, xldis, vpl, and xmu that are used in derivs 
c
         call MPI_BCAST(b, 6, MPI_DOUBLE_PRECISION, master,
     &                   MPI_COMM_WORLD, ierr)

      e10=b(1)
      xldis=b(2)
      vpl=b(3)
      xmu=b(4)
      dt=b(5)
      eps=b(6)
    
         call MPI_BCAST(ib, 3, MPI_INTEGER, master,
     &                   MPI_COMM_WORLD, ierr)


       nt=ib(1)
       ninit=ib(2)
       ndt=ib(3)
         
c     here's the end of slaves setup work

            nnn=idn

      endif           ! this is the end of all set up work 

      p1(1)=e10
      p1(2)=xldis
      p1(3)=vpl
      p1(4)=xmu


c
c      populate the arrays for the gather operation
c
      allocate (disps(numprocs), rcounts(numprocs), STAT=astat)

         rcounts(1)=iremain+iremain
         disps(1)=0

      do 177 i=2, numprocs

         rcounts(i)=idn+idn
         disps(i)=rcounts(i-1)+disps(i-1)
 177    continue

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
c     master does first set up work
c
         if(myid.eq.master) then

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                   ! not sure we use tm any more ++ TET fix
      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)

       endif                     ! end of masters first set up work

c       call MPI_BARRIER(MPI_COMM_WORLD, ierr)       ! stop the slaves from running away during debugging

      CALL derivs (t,stv,sdtvdt)   ! all processors get the starting derivatives for this time step

c
c      master does a second bit of set up work
c
      if(myid.eq.master) then

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

	   write (7,35) n,xc(n),zc(n),dx(n),dz(n)
           

C           The next lines are determining the cell mumber where A-B is nearly the 
C             most pos within the small group of cells to the North of the same depth 
C             where it is nearly the most negative. The cell where it is the most negative 
C             already has been 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

      endif                ! this is the end of the masters second bit of set up work


      do 87 n=1,nnn+nnn
   87     tvscal(n) = abs(stv(n)) + abs(dt*sdtvdt(n)) + tiny  ! ++ improve 

      CALL rkqs(stv,sdtvdt,nnn+nnn,t,dt,eps,tvscal,
     &  dtdid,dtnxt,derivs)  ! Integrate to the next time
      
       call MPI_Gatherv(stv,nnn+nnn, MPI_DOUBLE_PRECISION, rtv, rcounts,
     &   disps, MPI_DOUBLE_PRECISION, master,  MPI_COMM_WORLD, ierr)

c
c     if then else group where slaves send stv array back to the master
c      master recieves and puts all results into analogous un-key-sorted array tv     
c
      if(myid.eq.master)then      ! if master then

c     receive the results from the slaves
c
c     write them back into unsorted dependent variables
c     and their derivatives
c
       
         ii=1
         jj=1
       do 55 i=1,numprocs

          if(i.eq.1)then
             nn=iremain
            else
                nn=idn
             endif

         do 56 j=1,nn
         k=ii+nn
         tv(1+kindex(jj))=rtv(ii)
         tv(1+kindex(jj)+ne)=rtv(k)
         ii=ii+1
         jj=jj+1
         
c         dtvdt(1+kindex(i))=sdtvdt(i)
c         dtvdt(1+kindex(i)+ne)=sdtvdt(k)
c         tauss(1+kindex(i))=stauss(i)

 56      continue
         ii=k+1
55     continue

      else                   ! make sure slave knows the next time step size

        dt=dtnxt             ! set the next time step length to recommended value

      endif                  ! this is the end of if then else
c
c     rest of the work in the time loop is the masters
c
      if(myid.eq.master) then

      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)

      endif        ! end of masters work

  100 continue
c
C     END OF MAIN LOOP OVER TIME
C------------------------------------------------------------------------
      if(myid.eq.master)then   ! rest is masters work

      call cprogtime(timedate)
      write (6,1310)  timedate
      write (7,1310)  timedate
      write (8,1310)  timedate
      write (26,1310)  timedate
      write (27,1310)  timedate

! writing for a tecplot file
      write (50,*)'TITLE="History at Pos A-B and Most Neg A-b elements"'
      write (50,*)'VARIABLES="nt" "Time" "Element No." "Stress" "Velocit
     &y" "Displacement" "SS stress"'
      do 101 j=1,2
         if (j.eq.1) write (50,*)'ZONE T="Pos A-B"'       
         if (j.eq.2) write (50,*)'ZONE T="Most Neg A-B"'       
         do 101 i=1,nt-1
 101         write (50,*) (out(i,j,k),k=1,7)  ! writing for a tecplot file
      CLOSE (UNIT = 50)

      write (6,102) nwf1cnt
      write (8,102) nwf1cnt

      endif         ! end of masters work

  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,4f8.4,3f8.3)

 48   format (8(e10.4,1x), i12)
 49   format (2(e10.4,1x), i12)
c 49   format (e10.4)

   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      dx    ' 
     . '  dz     tauvp     sshr')
   77 format (/' non-bc cells:'//'       cell    xc      zc      dx    ' 
     . '  dz')
    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')


 9999  continue  ! this is used to stop the program during debugging

c      close(18)

       call MPI_BARRIER(MPI_COMM_WORLD, ierr)       ! stop the slaves from running away before master is done
                                                    ! writing output
       call MPI_ABORT(MPI_COMM_WORLD,errcode,ierr)  ! program will hang if mpi_finalize is used to end- fix this

      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)    ! returned - 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
      integer*4 sstrlen
      integer*4 lenfilenm,lendefalt
      character*30 filenm,defalt
!
!  NECESSARY PRE-CONDITIONS:
!   none
!
!  SIDE EFFECTS:
!	none
!
!  POST_CONDITIONS UPON EXIT
!	value of filenm is either changed to value read in or is left as the defalt
!
!  DESCRIPTION:
!   Called by main
!   Calls function sstrlen
!   Reads in an alphanumeric string from standard input, with length up to as long 
!    as 30 characters. Requests at standard output that this string be input and
!    tells one the name of the default value that will be used
!   Determines the length of the input alphanumeric string by calling function sstrlen,
!    a separate external function immediately following this subroutine.
!   If no string is entered from standard input, the default value prk.input is used.
!
!  EOP
!---------------------------------------------------------------------------------------------------
      lendefalt=14
      defalt='prk.input'
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***********************************************************************
!---------------------------------------------------------------------------------------------------
!  BOP
!
!  ROUTINE:  sstrlen
!
!  INTERFACE:
      function sstrlen(
     &   string,    ! supplied - alphanumeric string
     &   ilen)      ! returned - number of characters in that string

!
!  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
      integer sstrlen,imax,ilen
      character string*80

!
!  NECESSARY PRE-CONDITIONS:
!   none
!
!  SIDE EFFECTS:
!	none
!
!  POST_CONDITIONS UPON EXIT
!	value of filenm is either changed to value read in or is left as the defalt
!
!  DESCRIPTION:
!   Called by READ_FILE_NAME
!   Calls nothing
!   Determines the number of characers in the input alphanumeric string, max of 80 allowed
!   Returns that number of characters
!
!  EOP
!---------------------------------------------------------------------------------------------------
	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***********************************************************************
!---------------------------------------------------------------------------------------------------
!  BOP
!
!  ROUTINE:  READ_CELL_INFO
!
!  INTERFACE:
c  --- input variables this code --------

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

      SUBROUTINE READ_CELL_INFO(
     .   ibca,      ! returned - no. of first cell in each of ngrp's of bc cells
     .   ibcb,      ! returned - no. of last cell in each of ngrp's of bc cells
     .   w,         ! returned - width, km of each cell 
     .   d,         ! returned - depth, km of each cell
     .   dx,        ! returned - cell width, km for all cells in this block
     .   dz,        ! returned - cell height, km for all cells in this block   
     .   title,     ! returned - title of model
     .   xmu,       ! returned - rigidity, poisson solid (Megabars)
     .   nblk,      ! returned - no. of blocks of identical cells 
     .   nf,        ! returned - total no. of fault cells, tot. of bc and non bc; nf = ne + nb
     .   ndt,       ! returned - maximum no. of time steps
     .   nbc,       ! returned - =1 some cells are b.c. cells w/ vpl
                    !            =0 entire fault outside patch slips at vpl
     .   nwp1,      ! returned - step at which to start writing print file prk.prt (7); =0, no print
     .   wpdt,      ! returned - write prk.prt if ge this interval in yrs has passed; wpdt=0, write every t step
     .   nwc1,      ! returned - write surface creep rate if > 0; =0, no write
     .   wcdt,      ! returned - write prk.crp (9) at intervals ge wcdt yrs
     .   nwf1,      ! returned - step at which to start writing binary file prk.flt (4); =0, no write
     .   wfdt,      ! returned - write .flt at intervals ge wfdt yrs
     .   dt1,       ! returned - trial dt for first t step, yrs
     .   eps,       ! returned - max err, see tvscal, subr rkqc  ca. 1.e-4 or 1.e-5
     .   vpl,       ! returned - plate velocity used for boundary conditions
     .   factinit,  ! returned - factor times vlim that is used for initial condition of
                    !            velocity in velocity weakening area - typically use 0.01
     .   ngrp)      ! returned - no. of cell range pairs for cells having vpl bndy. cond.

!
!  USES:
	use cpat    
!
!  RETURN VALUE:
c	implicit none   ! Not yet converted to using "implicit none" and explicitly typing all variables
c	<type> :: <name>			! Description
CC        nblk is the dimension of several arrays: currently is 40
C
      implicit real*8 (a-h,o-z)
      parameter (NF_MAX=160003)

      common /pat/ vlim,npat,nb,ibc(20)

      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)

!
!  NECESSARY PRE-CONDITIONS:
!   none
!
!  SIDE EFFECTS:
!	none
!
!  POST_CONDITIONS UPON EXIT
!	value of filenm is either changed to value read in or is left as the defalt
!
!  DESCRIPTION:
!   Called by main
!   Calls nothing
!	This subroutine reads in some of the data from the .dat file.
!	 This includes the cell information, the "Simulation Parameters", 
!	 and the list of which cells are the boundary-condition cells. 
!	 The rest of the information from the .dat file, the friction law 
!	 parameters, how they are distributed on the fault and the 
!	 creepmeter names and corresponding cell numbers, are read 
!	 in by the main program.
!
!  EOP
!---------------------------------------------------------------------------------------------------

      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
C         read(2,*) nx(n),nz(n),dxb(n),dzb(n),xul(n),zul(n)  ! Old way of inputting xul and zul
         read(2,*) nx(n),nz(n),dxb(n),dzb(n),xulnum,xuldenom,  ! New way of inputting xul and zul
     .                                       zulnum,zuldenom   ! allows inputting fractions w/ precision
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
         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)
	       xul(n) = xulnum/xuldenom
	       zul(n) = zulnum/zuldenom
               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***********************************************************************
!---------------------------------------------------------------------------------------------------
!  BOP
!
!  ROUTINE:  GREENS_FUNCTION
!
!  INTERFACE:
 	  SUBROUTINE GREENS_FUNCTION(
     &   m,         ! supplied - index number of source element
     &   n)         ! supplied - index number of observation element
!
!  USES:
      USE Green_sshr	! Make data in module "Green_sshr" available
      use cderivs
      use cpat
!
!  RETURN VALUE:
c	implicit none   ! Not yet converted to using "implicit none" and explicitly typing all variables
c	<type> :: <name>			! Description

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

      common /deriv/ e10,xldis,vpl,xmu,ne,ne2,MMP

      common /pat/vlim,npat,nb,ibc(20)

!
!  NECESSARY PRE-CONDITIONS:
!   none
!
!  SIDE EFFECTS:
!	none
!
!  POST_CONDITIONS UPON EXIT
!	The value of sshr(m,n) in the Green_sshr module is changed
!
!  DESCRIPTION:
!   Called by main and, if MMP = 0, by DERIVS subroutine
!   Calls STRN
!   For unit left lateral slip in each cell, find the shear stress at all cell centers
!
!  EOP
!---------------------------------------------------------------------------------------------------

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***********************************************************************
!---------------------------------------------------------------------------------------------------
!  BOP
!
!  ROUTINE:  STRN
!
!  INTERFACE:
      SUBROUTINE STRN (
     &   xw2,       ! supplied - half width of source cell (i.e. horizontally along fault)
     &   xd2,       ! supplied - half height of source cell (i.e. in depth direction)
     &   x3,        ! supplied - center of source cell coordinate, downward
     &   y1,        ! supplied - observation point coordinate, along fault
     &   y3,        ! supplied - observation point coordinate, downward
     &   e12)       ! returned - tensor strain e12 on fault plane at obs. pt.
!
!  RETURN VALUE:
c	implicit none   ! Not yet converted to using "implicit none" and explicitly typing all variables
c	<type> :: <name>			! Description

      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./
!
!  NECESSARY PRE-CONDITIONS:
!   none
!
!  SIDE EFFECTS:
!	none
!
!  POST_CONDITIONS UPON EXIT
!   none
!
!  DESCRIPTION:
!   Called by GREENS_FUNCTION
!   Calls nothing
!   Determines shear strain at observation point due to unit slip on source cell
!       Equations from Chinnery (1963)   
!       analytic soln pers. com. 1983 WDS; xlam=xmu
!       Coordinates: 1 parallel to flt,  2 normal to flt,  3 down
!       x (source cell center),  y (obervation pt)  
!       1 and 2 coordinates of source point = 0 as is 2 coordinate of obs. pt.
!        (Namely, only relative distances along fault matter, but absolute depths matter,
!        so for source point only x3 need be given, and y1-x1 can become y1 with x1=0.)
!
!  EOP
!---------------------------------------------------------------------------------------------------

      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***********************************************************************
!---------------------------------------------------------------------------------------------------
!  BOP
!
!  ROUTINE:  KOMPRS
!
!  INTERFACE:
      SUBROUTINE KOMPRS (
     &   nb,        ! supplied - # of bc cells
     &   ibctmp,    ! supplied and returned - ibctmp(i) is the index of the ith one of the bc cells
     &   nf)        ! supplied - # of bc and non-bc cells

!
!  USES:
      USE Green_sshr	! Make data in module "Green_sshr" available
      use cderivs
!
!  RETURN VALUE:
c	implicit none   ! Not yet converted to using "implicit none" and explicitly typing all variables
c	<type> :: <name>			! Description

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

      common /deriv/ e10,xldis,vpl,xmu,ne,ne2,MMP

      dimension ibctmp(1)
!
!  NECESSARY PRE-CONDITIONS:
!   none
!
!  SIDE EFFECTS:
!	none
!
!  POST_CONDITIONS UPON EXIT
!	The sshr(m,n) array in the Green_sshr module is changed
!   ibctmp(i) is updated to reflect the changed indices
!
!  DESCRIPTION:
!   Called by main
!   Calls nothing
!   Compress sshr matrix by removing booundary condition rows and columns
!    (However, it is convenient to number the bc cells last, so then this 
!     is not really needed. Doing that initial numbering saves confusion 
!     between cell's pre- and post-compression indices.)
!
!  EOP
!---------------------------------------------------------------------------------------------------

      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***********************************************************************
!---------------------------------------------------------------------------------------------------
!  BOP
!
!  ROUTINE:  PATCH4
!
!  INTERFACE:
      SUBROUTINE PATCH4 (
     &   tv,        ! returned - array of stresses and velocities
     &   factinit)  ! supplied - factor times vlim that is used for initial condition of
                    !            velocity in velocity weakening area - typically use 0.01
       use cderivs
       use cpat
!
!  RETURN VALUE:
c	implicit none   ! Not yet converted to using "implicit none" and explicitly typing all variables
c	<type> :: <name>			! Description

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

      common /deriv/  e10,xldis,vpl,xmu,ne,ne2,MMP

      common /pat/ vlim,npat,nb,ibc(20)

      dimension tv(1), pxo(5),pzo(5),ax(5),az(5),abmn(5),abmx(5)
      parameter (pi=3.1415926)
!
!  NECESSARY PRE-CONDITIONS:
!   none
!
!  SIDE EFFECTS:
!	none
!
!  POST_CONDITIONS UPON EXIT
!	The tv arry is set up for the initial conditions
!
!  DESCRIPTION:
!   Called by main
!   Calls INIT
!
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
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( , ). This is no longer used I believe,
C        now that the radiation damping approach is used.
!
!  EOP
!---------------------------------------------------------------------------------------------------

      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***********************************************************************
!---------------------------------------------------------------------------------------------------
!  BOP
!
!  ROUTINE:  PATCH5
!
!  INTERFACE:
      SUBROUTINE PATCH5 (
     &   tv,        ! returned - array of stresses and velocities
     &   iout6,     ! supplied - a flag - 1 if output is made to unit 6, 0 if not
     &   iout7,     ! supplied - a flag - 1 if output is made to unit 7, 0 if not
     &   irestart,  ! supplied - a flag - 1 if this run is to start from some time step with data output
     &              !            from an previous program, 0 if just starting from beginning
     &   factinit,  ! supplied - factor times vlim that is used for initial condition of
                    !            velocity in velocity weakening area - typically use 0.01
     &   a_and_amb, ! supplied - flag to say if a and amb are output for plotting check: 1 yes, 0 no
     &   nminc)     ! returned - cell no. of min A-B

      use cderivs
      use cpat
!
!  RETURN VALUE:
c	implicit none   ! Not yet converted to using "implicit none" and explicitly typing all variables
c	<type> :: <name>			! Description

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

      INTEGER :: a_and_amb      ! Flag to say if a and amb are output for plotting check
      parameter (NF_MAX=160003)

      common /deriv/  e10,xldis,vpl,xmu,ne,ne2,MMP

      common /pat/ vlim,npat,nb,ibc(20)

      dimension tv(1),Temp(NF_MAX)
!
!  NECESSARY PRE-CONDITIONS:
!   none
!
!  SIDE EFFECTS:
!	none
!
!  POST_CONDITIONS UPON EXIT
!	The tv arry is set up for the initial conditions
!   The values of A and A-B for each cell are optionally put in a binary file for plotting
!   THe cell number of the cell with minimum A-B id returned
!
!  DESCRIPTION:
!   Called by main
!   Calls INIT
!
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 or Dc of Dieterich, etc.
c       depthfact - scale factor applied to L&S z(T); commonly = .73
!
!  EOP
!---------------------------------------------------------------------------------------------------

      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
       
       use cderivs
       use cpat

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

      common /deriv/ e10,xldis,vpl,xmu,ne,ne2,MMP

      common /pat/ vlim,npat,nb,ibc(20)

      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 ( NE2_MAX=320006)
      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 velocity 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***********************************************************************
C  What follows below is the NASA format for Prologues for different types of 
C    program units. These have been implementd, at least to a first 
C    approximation, in the preceeding module, main program, and subrountines.
C  The lines beginning with C! are comments in the prologues. The lines
C    beginning with ! would be functioning lines in the prologues.
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