c--- program hazSUBXv3 last changed 3/25/03
c--- reads sd7 for Sadigh in second line
c----added distance taper of weighting
c---- Sadigh et al. for M9; set it to M8.5
c--- fixed M9 rupture to not float
c--- modified to specify fault length for subduction zone EQ
c--- also increased various array dimensions
c---  hazSUBX.f
c--- calculates mean annual frequency of exceedance for dipping faulta
c--- written by Art Frankel and Rob Wesson 
c--- this version handles multiple periods and gnd motion relations
c--- to compile: cc -c iosubs.c
c                 f77 -o hazSUBX hazSUBX.f iosubs.o
c--- to run: hazFX inputfile
c-----  your choice of attenuation relations: -BJF, Toro, Sadigh, Campbell
c-----  Atkinson and Boore(95) look up table, Crouse, Geomatrix
c--- output of this program is input into hazallX to get probabilistic
c    ground motions
c--- ground motion levels should be in units of g
c--- period=0 indicates PGA
      dimension x(500,2),y(500,2)
      dimension dmin2(50,1001,9)
      dimension dmin(9),dminr(9),nmag0(500)
      type header_type
         character*30 :: name(6)
         real*4 :: period
         integer*4 :: nlev
         real*4 :: xlev(20)
         real*4 :: extra(10)
      end type header_type
      type(header_type) headr
      real magmin(500),magmax(500)
      integer*4 readn
      character nameout*256,name*256,nametab*256,adum*256
      dimension p(1005,35,9,20,3,3)
      dimension prob(250,20,9),xlev(20,9),out(800000), nlev(9),iatten
     +(9,5),iconv(9,5)
      dimension nattn(9),period(9)
      dimension ifp(9)
      dimension b1(9),b2(9),b3(9),b4(9),b5(9),h(9),sigma(9)
      dimension gc1(9),gc2(9),gc3(9),gc4(9),gc5(9)
      dimension sc1(9),sc2(9),sc3(9),sc4(9),sc5(9),sc6(9),sh(9)
      dimension sd1(9),sd2(9),sd3(9),sd4(9),sd5(9),sd6(9),sd7(9)
      dimension ssigma1(9),ssigmacoef(9),ssigma2(9),smagsig(9)
      dimension cc1(9),cc2(9),cc3(9),cc4(9),cc5(9),ch(9)
      dimension cd1(9),cd2(9),cd3(9),cs1(9),cs2(9)
      dimension csigma1(9),csigmacoef(9),csigma2(9),cmagsig(9)
      dimension tlen(500),a(500),b(500),dmag(500)
      dimension itype(500),iftype(500),cmag(500),crate(500)
      dimension dip1(500),width1(500),depth0(500),npts(500)
      dimension dip2(500),width2(500)
      dimension thrust(9),sthrust(9)
      real*4 norm(9)
      dimension ab1(9),ab2(9),ab3(9),ab4(9)
      dimension gma(22,30,20),amag(20),rlog(20),f(20)
      dimension p1(9),p2(9),p3(9),p4(9),p5(9),p6(9),p7(9)
      dimension psigma(9),ph(9),absigma(9),abdepth(9)
      dimension xlogfac(9),bsigma(9),bdepth(9)
      dimension rate(50)
      dimension nrup(50),mm(50)
      dimension tarray(2)
      dimension wtdist(7,9),wt(7,9,2)
 
c      external fltdst
 
c      write(6,*) "enter name of input file"
      call getarg(1,name)
      open(unit=1,file=name)
c      write(6,*) "enter name of output file"
      write(6,*) "for sites: enter min lat, max lat, dlat"
      read(1,*) ymin, ymax, dy
      write(6,*) "for sites: enter min lon, max lon, dlon"
      read(1,*) xmin, xmax, dx
      write(6,*) xmin, xmax, dx
      nx= ceiling((xmax-xmin)/dx +1)
      ny= ceiling((ymax-ymin)/dy +1)
      write(6,*) nx,ny
      pi= 3.14159
c      read(5,*) idum
      nrec= nx*ny
cccccccccccccccccccccccccccccccccc
      write(6,*) "enter number of periods"
      read(1,*) nper
c---loop through periods
      do 40 ip=1,nper
         read(1,*) period(ip)
         write(6,*) "enter name of output file for this period"
         read(1,909) nameout
 909     format(a78)
         call openwx(ifp(ip),nameout)
         write(6,*) ip,ifp(ip)
         write(6,*) "enter number of atten. relations for this period"
         read(1,*) nattn(ip)
c--- loop through atten relations for that period
         do 20 ia=1,nattn(ip)
            write(6,*)
     +      "enter type of atten. relation, weight, mb to M conv."
            read(1,*) iatten(ip,ia),wt(ip,ia,1),wtdist(ip,ia),
     +      wt(ip,ia,2),iconv(ip,ia)
c---- Boore-Joyner type of attenuation relation
            if(iatten(ip,ia).eq.1) then
               write(6,*) "enter b1,b2,b3,b4,b5,h,sigma,thrust,norm"
               read(1,*) b1(ip),b2(ip),b3(ip),b4(ip),b5(ip),h(ip),sigma
     +         (ip), thrust(ip),norm(ip)
            endif
c--------  Geomatrix interface atten.
            if(iatten(ip,ia).eq.2) then
               write(6,*) "enter c1,c2,c3,c4,c5"
               read(1,*) gc1(ip),gc2(ip),gc3(ip),gc4(ip),gc5(ip)
            endif
c------- Sadigh type
            if(iatten(ip,ia).eq.3) then
               write(6,*) "enter c1,c2,c3,c4,c5,c6,h,thrust"
               read(1,*) sc1(ip),sc2(ip),sc3(ip),sc4(ip),sc5(ip), sc6
     +         (ip),sh(ip),sthrust(ip)
               read(1,*) sd1(ip),sd2(ip),sd3(ip),sd4(ip),sd5(ip),sd6(ip)
     +         ,sd7(ip)
               read(1,*) ssigma1(ip),ssigmacoef(ip),ssigma2(ip),smagsig
     +         (ip)
            endif
c------- Campbell type
            if(iatten(ip,ia).eq.4) then
               write(6,*) "enter c1,c2,c3,c4,c5,h,thrust"
               read(1,*) cc1(ip),cc2(ip),cc3(ip),cc4(ip),cc5(ip),ch(ip)
               read(1,*) cd1(ip),cd2(ip),cd3(ip)
               read(1,*) cs1(ip),cs2(ip)
               read(1,*) csigma1(ip),csigmacoef(ip),csigma2(ip),cmagsig
     +         (ip)
            endif
c------ Atkinson and Boore (95) look-up table
            if(iatten(ip,ia).eq.5) then
               write(6,*) "enter  coeficients, depth, sigma"
               read(1,*) ab1(ip),ab2(ip),ab3(ip),ab4(ip), abdepth(ip),
     +         absigma(ip)
               call gettab(gma,f,amag,rlog)
            endif
c------ Crouse attenuation form
            if(iatten(ip,ia).eq.6) then
               write(6,*) "enter p1,p2,p3,p4,p5,p6,p7,sigma,depth"
               read(1,*) p1(ip),p2(ip),p3(ip),p4(ip),p5(ip),p6(ip), p7
     +         (ip),psigma(ip),ph(ip)
            endif
c----- new Boore look-up table
            if(iatten(ip,ia).eq.7) then
               write(6,*) "enter file name of table"
               read(1,380) nametab
               write(6,*) "enter depth, sigma, log factor"
               read(1,*) bdepth(ip),bsigma(ip),xlogfac(ip)
               open(unit=15,file=nametab,status='old')
               read(15,380) adum
               do 10 idist=1,21
   10          read(15,*) xdum,(gma(imag,idist,ip),imag=1,20)
               close(15)
            endif
   20    continue
         totwt1 = 0.
         totwt2 = 0.
         do ia=1,nattn(ip)
            totwt1 = totwt1 + wt(ip,ia,1)
            totwt2 = totwt2 + wt(ip,ia,2)
         enddo
         if (totwt1.ne.1) write(*,*)
     &      "ERROR: wt1 is greater than zero and equal to", totwt1
         if (totwt2.ne.1) write(*,*)
     &      "ERROR: wt2 is greater than zero and equal to", totwt2

         write(6,*) "enter number of ground motion levels"
         read(1,*) nlev(ip)
         write(6,*) "enter ground motion levels"
         read(1,*) (xlev(k,ip),k=1,nlev(ip))
         do 30 k=1,nlev(ip)
   30    xlev(k,ip)= alog(xlev(k,ip))
   40 continue
ccccccccccccccccccccccccccccccccccc
      write(6,*) "enter increment dlen for source segment"
      read(1,*) dlen
      write(6,*) "enter distance increment, dmax"
      read(1,*) di,dmax
ccc-------read fault data
      do 50 ift=1,1000
         write(6,*) ift
         write(6,*)
     +   "enter 1 for char. 2 for G-R; 1 for SS, 2 for reverse"
         read(1,*,end=60) itype(ift),iftype(ift)
         if(itype(ift).eq.1) then
            write(6,*) "enter char. M and rate"
            read(1,*) cmag(ift),crate(ift)
         endif
         if(itype(ift).eq.2) then
            write(6,*) "enter a,b,min M,max M,dmag for freq mag curve"
            read (1,*) a(ift),b(ift),magmin(ift),magmax(ift),dmag(ift)
            nmag0(ift)= (magmax(ift)-magmin(ift))/dmag(ift) + 1.4
         endif
c******************************************************************
c      Replace code to read and setup model of subduction zone
c******************************************************************
c      write(6,*) "enter dip1, updip width, depth0 of lat, lons"
c      read(1,*) dip1(ift), width1(ift), depth0(ift)
c      write(6,*) "enter dip2, downdip width"
c      read(1,*) dip2(ift), width2(ift)
c      write(6,*) "enter number of segment points"
c      read(1,*) npts(ift)
c      write(6,*) "enter lat,lon for each point"
c      tlen(ift)= 0.
c      do 2 j=1,npts(ift)
c      read(1,*) y(j,ift),x(j,ift)
c      if(j.gt.1) then
c        sy= y(j-1,ift)
c        sx= x(j-1,ift)
c        ry= y(j,ift)
c        rx= x(j,ift)
c        if((sy.eq.ry).and.(sx.eq.rx)) then
c             delta=0.
c             az=0.
c             go to 1020
c        else
c           call delaz(sy,sx,ry,rx,delta,az,baz)
c        endif
c 1020   if(j.gt.2) xlen(j-1,ift)= delta +xlen(j-2,ift)
c        if(j.eq.2) xlen(1,ift)= delta
c        xaz(j-1,ift)= az
c        tlen(ift)= tlen(ift)+ delta
c        endif
c      write(6,*) xlen(j-1),delta,xaz(j-1)
c   2  continue
c***************************************************************
c     Insert new code to read and setup model
c***************************************************************
         call readmodel(ift,tlen)
 
         write(*,*) "tlen= ", tlen(ift)
 
 
c****************************************************************
c     End new code to read and setup model
c****************************************************************
         write(6,*) npts(ift),tlen(ift)
         coef= 3.14159/180.
c      dip1(ift) =dip1(ift)*coef
c      dip2(ift) =dip2(ift)*coef
   50 continue
   60 nft= ift-1
      write(6,*) "nft=",nft
ccccccccccccccccccccccccccccc
c---- write header --------------
      do 80 ip=1,nper
      headr%name(1)= name
      headr%period= period(ip)
      headr%nlev= nlev(ip)
         do 70 k=1,nlev(ip)
            headr%xlev(k)= exp(xlev(k,ip))
   70    continue
         ndata= 308
         call puthead(ifp(ip),headr,ndata,readn)
   80 continue
ccccccccccccccccccccccccccccccccccc
      ndist= dmax/di +1.5
      nmag= 35
      do 90 ii=1,ndist
         do 90 m=1,nmag
            do 90 ip=1,nper
               do 90 k=1,nlev(ip)
                  do 90 ifkey=1,3
                     do 90 iatype=1,3
   90 p(ii,m,ip,k,ifkey,iatype)=0.
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=35 dmag=0.1
c-- loop through distances
      do 170 ii=1,ndist
c      write(6,*) "dist counter=",ii
         dist0= (ii-1)*di
         if(dist0.gt.dmax) go to 170
c-- loop through magnitudes
         do 160 m=1,nmag
            xmag0= 6.0+(m-1)*0.1
c--- loop through periods
            do 150 ip=1,nper
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
               do 140 ia=1,nattn(ip)
                  xmag= xmag0
                  if(iconv(ip,ia).eq.1) xmag= 2. 715-0.277*xmag0 +0.127
     +            *xmag0*xmag0
c-- following for Joyner Boore WUS
                  if(iatten(ip,ia).eq.1) then
                     sigmasq= sigma(ip)*sqrt(2.)/0.4343
                     dist= sqrt(dist0*dist0+h(ip)*h(ip))
                     gnd= b1(ip)+b2(ip)*(xmag-6.)+b3(ip)*(xmag-6.) **2
 
                     gnd= gnd+ b4(ip)*dist + b5(ip)*alog10(dist)
c--- convert from psv in cm/sec to psa in g
                     if(period(ip).ne.0.) gnd=gnd+alog10(2.*pi /(980.
     +               *period(ip)))
                     gnd2= gnd+ thrust(ip)
                     gnd3= gnd+ norm(ip)
                     gnd= gnd/0.4343
                     gnd2=gnd2/0.4343
                     gnd3= gnd3/0.4343
                  endif
c--- following for Geomatrix (1993) Youngs and Coppersmith variation
c---  depth set= 20 km
                  if(iatten(ip,ia).eq.2) then
                     sig= gc4(ip)+gc5(ip)*xmag
                     if(xmag.gt.8.) sig= gc4(ip)+gc5(ip)*8.
                     sigmasq= sig*sqrt(2.)
                     gnd= 0.3633+1.414*xmag+ gc1(ip)+gc2(ip) *((10.
     +               -xmag)**3)
                     arg= 0.554*xmag
                     gnd=gnd+gc3(ip)*alog(dist0+1.782*exp(arg))
                     gnd2= gnd
                     gnd22= gnd2
                     gnd3=gnd
                  endif
c--- following is Sadigh form
                  if(iatten(ip,ia).eq.3) then
                     if(xmag.lt.smagsig(ip)) sig= ssigma1(ip)-
     +               ssigmacoef(ip)*xmag
                     if(xmag.ge.smagsig(ip)) sig= ssigma2(ip)
                     sigmasq= sig*sqrt(2.)
                     if(xmag.lt.6.5) then
                        gnd= sc1(ip)+sc2(ip)*xmag+sc3(ip) *((8.5-xmag)
     +                  **2.5)
                        gnd= gnd+ sc4(ip)*alog(dist0+ exp(sc5(ip) +sc6
     +                  (ip)*xmag))
                     endif
                     if((xmag.ge.6.5).and.(xmag.lt.8.5)) then
                        gnd= sd1(ip)+sd2(ip)*xmag+sd3(ip) *((8.5-xmag)
     +                  **2.5)
                        gnd= gnd+ sd4(ip)*alog(dist0+exp(sd5(ip) +sd6
     +                  (ip)*xmag))
                     endif
                     if(xmag.ge.8.5) then
                        gnd= sd1(ip)+sd2(ip)*8.5
                        gnd= gnd+ sd4(ip)*alog(dist0+exp(sd5(ip) +sd6
     +                  (ip)*8.5))
                     endif
                     gnd= gnd+ sd7(ip)*alog(dist0+2.)
                     gnd2= gnd+sthrust(ip)
                     gnd3= gnd
                     if(xmag.eq.7.9) write(12,*) dist0,period(ip),
     +                exp(gnd2),
     +               exp(gnd22),xmag
  100                continue
                  endif
c--- following is Campbell form
                  if(iatten(ip,ia).eq.4) then
                     if(xmag.lt.cmagsig(ip)) csigma= csigma1(ip)-
     +               csigmacoef(ip)*xmag
                     if(xmag.ge.cmagsig(ip)) csigma= csigma2(ip)
                     sigmasq= csigma*sqrt(2.)
                     arg= dist0**2+ (cc4(ip)*exp(cc5(ip)*xmag))**2
                     arg= sqrt(arg)
                     gnd= cc1(ip)+cc2(ip)*xmag+cc3(ip)*alog(arg)
                     if(dist0.eq.0.) then
                        gnd=0
                        gnd2=0
                        gnd3=0
                        go to 110
                     endif
                     gnd= gnd+ cs1(ip)+ cs2(ip)*alog(dist0)
                     gnd3= gnd
c--- following for Campbell thrust fault
                     gnd2= gnd+ cd1(ip)+cd2(ip)*alog(dist0)+cd3 (ip)
     +               *xmag
  110                continue
                  endif
c--- following for Atkinson and Boore (95) look-up table
                  if(iatten(ip,ia).eq.5) then
                     sigmasq= absigma(ip)*sqrt(2.)/0.4343
                     dist= dist0
                     if(dist.lt.10.) dist=10.
                     if(period(ip).gt.0.) freq=1./period(ip)
                     if(period(ip).eq.0.) freq= 99.
c---------if M> 7.25 use quadratic formula
                     if(xmag.gt.7.25) then
                        gnd= ab1(ip)+ab2(ip)*(xmag-6.)+ab3(ip)
     +                  *((xmag-6)**2) -alog10(dist)- ab4(ip)*dist
                        go to 120
                     endif
                     call table(gma,xmag,dist,freq,gnd,f,amag, rlog)
 
cc- convert to fraction of g
  120                fac= -alog10(980.)
                     gnd= gnd+fac
                     gnd= gnd/0.4343
                     gnd2= gnd
                     gnd3= gnd
                  endif
c--------following for Crouse(1991) atten.
                  if(iatten(ip,ia).eq.6) then
                     sigmasq= psigma(ip)*sqrt(2.)
                     dist= dist0
                     gnd= p1(ip)+p2(ip)*xmag+p3(ip)*xmag*xmag+ p4 (ip)
     +               *alog(dist + p5(ip)*exp(p6(ip)*xmag))+p7 (ip)*ph
     +               (ip)
c--- convert from psv in cm/sec to psa in g
                     if(period(ip).ne.0.) gnd=gnd+alog(2.*pi/(980.
     +               *period(ip)))
                     gnd2=gnd
                     gnd3= gnd
                  endif
c---- following for new Boore look-up table
                  if(iatten(ip,ia).eq.7) then
                     sigmasq= bsigma(ip)*sqrt(2.)/0.4343
                     dist= dist0
                     if(dist.lt.10.) dist=10.
                     imag= (xmag-4.4)/0. 2+1
                     rdist= alog10(dist)
                     idist= (rdist-0.70)/0. 1+ 1
                     xm1= (imag-1)*0. 2+ 4.4
                     fracm= (xmag-xm1)/0.2
                     xd1= (idist-1)*0. 1+0.7
                     fracd= (rdist-xd1)/0.1
                     gm1= gma(imag,idist,ip)
                     gm2= gma(imag+1,idist,ip)
                     gm3= gma(imag,idist+1,ip)
                     gm4= gma(imag+1,idist+1,ip)
                     arl= gm1 + fracm*(gm2-gm1)
                     aru= gm3 + fracm*(gm4-gm3)
                     gnd= arl +fracd*(aru-arl)
                     gnd= gnd+ xlogfac(ip)
                     gnd= gnd/0.4343
                  endif
c---- loop thru ground motion levels
                  iatype= iatten(ip,ia)
                  do 130 k=1,nlev(ip)
                     temp1= (gnd- xlev(k,ip))/sigmasq
                     temp2= (gnd2- xlev(k,ip))/sigmasq
                     temp3= (gnd3- xlev(k,ip))/sigmasq
c                     p(ii,m,ip,k,1,iatype)= (erf(temp1)+1.)/2.
c                     p(ii,m,ip,k,2,iatype)= (erf(temp2)+1.)/2.
c                     p(ii,m,ip,k,3,iatype)= (erf(temp3)+1.)/2.
      p(ii,m,ip,k,1,iatype)= (eerf(temp1)+1.)/2.
c      write(*,*) p(ii,m,ip,k,1,iatype), temp1
      p(ii,m,ip,k,1,iatype)= (p(ii,m,ip,k,1,iatype)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,1,iatype).lt.0.) p(ii,m,ip,k,1,iatype)=0.
      p(ii,m,ip,k,2,iatype)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,iatype)= (p(ii,m,ip,k,2,iatype)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,2,iatype).lt.0.) p(ii,m,ip,k,2,iatype)=0.
      p(ii,m,ip,k,3,iatype)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,iatype)= (p(ii,m,ip,k,3,iatype)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,3,iatype).lt.0.) p(ii,m,ip,k,3,iatype)=0.
  102 continue
  130             continue
  140          continue
  150       continue
  160    continue
  170 continue
ccccccccccccccccccccccccccccccccccccccccccc
      dum=dtime(tarray)
      write(*,*) 'Finished setting up erf matrix'
      write(*,*) tarray
ccccccccccccccccccccccccccccccccccc
      icnt=1
      do 180 j=1,250 
         do 180 ip=1,nper
            do 180 k=1,nlev(ip)
  180 prob(j,k,ip)=0.
ccccccccccccccccccccccccccccccccccc
c---Assume only one fault
      ift =1
c---Set up parameters of rupture zones
      do 185 m=1,nmag0(ift)
          xmag = magmin(ift)+(m-1)*dmag(ift)
c------find rupture length from Geomatrix interface eqs relation
          ruplen= (xmag-4.94)/1.39
          ruplen= 10**ruplen
          nrup(m)= (tlen(ift)-ruplen)/dlen +1.
          write(6,*) "nrup(m),ruplen,tlen",nrup(m),ruplen,tlen(ift)
          if(nrup(m).gt.1000) then
                write(6,*) "RUPTURES GREATER THAN 1000"
                stop
          endif
          if(ruplen.ge.tlen(ift)) then
		nrup(m)=1
                ruplen=tlen(ift)
          endif
          rate(m)= a(ift) - b(ift)*xmag
          rate(m)= 10.**rate(m)
      write(6,*) "rate(m)=",rate(m)
c       read(5,*) idum
          rate(m)= rate(m)/float(nrup(m))
          mm(m)= (xmag-6.)/0.1+1.5
          do 182 irup=1,nrup(m)
               xmove= dlen*(irup-1)
               s1 = xmove
               s2 = xmove + ruplen
               call rupturesetup(m,irup,s1,s2)
 182      continue
 185   continue
c
c---Here's the guts
c
c---loop through receiver sites 250 at a time
      i = 0
      nloop = (nrec-1) / 250+ 1
      do 370 iloop=1,nloop
 
c     loop through floating rupture zones
c     find minimum distance dmin2 to receiver for each rupture zone
c     nrup(m) rupture zones per magnitude bin m
c     move rupture zones by dlen horizontally along fault
 
c*****calculate distances for 250 receivers from all sources*********
        
         mloop = 250
 
         do 310 jloop=1,mloop
                  icnt = jloop
c                       i is the index of the receiver
                  i = (iloop-1)*250 + jloop
                  if(i.gt.nrec) then
                     i = i - 1
                     icnt = jloop -1          
                     go to 320
                  endif
 
                  iy= (i-1)/nx
                  ix= i-1-iy*nx
                  rx= xmin+(ix)*dx
                  ry= ymax-(iy)*dy
            do 300 m=1,nmag0(ift)
               do 290 irup=1,nrup(m) 
c-----------------------------------------------------------------------
c  dmin is minimum distance from fault to receiver
                  call qkdist(m,irup,rx,ry,distmin,fltdepmin, distminxy)
 
                  dmin2(m,irup,1) = distminxy
c                  write(*,*) "m,irup,rx,ry,distmin",m,irup,rx,ry,distmin
c        dmin2(m,irup,2) = fltdepmin
c                  dmin2(m,irup,2) = sqrt(distminxy *distminxy +
c     +            400.)
                  dmin2(m,irup,2) = distmin
                  dmin2(m,irup,3) = distmin
                  dmin2(m,irup,4) = distmin
 
c*****calculate hazard for 250 receivers from all sources*********
 
                           do 280 ip=1,nper
                           do 270 ia=1,nattn(ip)
                              iatype= iatten(ip,ia)
                              iatype0= iatype
                              if(iatype.gt.4) iatype=3
                              ii= dmin2(m,irup,iatype)/di +1.5
                              if(ii.gt.ndist) go to 270
         weight= wt(ip,ia,1)
c       if(weight.lt.0.and.ia.eq.1)
c    &     write(6,*) ia,ii,weight,dmin2(m,irup,iatype),wtdist(ip,ia),
c    &     wt(ip,ia,2),wt(ip,ia,1)
         if(dmin2(m,irup,iatype).gt.wtdist(ip,ia)) weight=wt(ip,ia,2)
         wt2= wt(ip,ia,2)
         wt1= wt(ip,ia,1)
         dist0= wtdist(ip,ia)
         dist= dmin2(m,irup,iatype)
         taper= 15.
         if((dist.ge.dist0-taper).and.(dist.le.dist0+taper)) 
     &     weight= wt1 + (wt2-wt1)*(dist-dist0+taper)/(2.*taper)
c        if(weight.lt.0) write(6,*) ia,ii,weight
                              do 260 k=1,nlev(ip)
c        write(*,*) ii,k,p(ii,mm(m),ip,k,iftype(ift),iatype0),
c    &             prob(icnt,k,ip),weight,rate(m)
  260                         prob(icnt,k,ip)= prob(icnt,k,ip)+ 
     +                        weight*rate(m)*p(ii,mm(m),ip,k,
     +                        iftype(ift),iatype0)
  270                      continue
  280                   continue
  290                continue
  300             continue
  310          continue
               go to 320
c---------------------------
c--- for characteristic event, tested
c                           if(itype(ift).eq.1) then
c                              xmag= cmag(ift)
c                              rate(1)= crate(ift)
c                              mm(1)= (xmag-6.)/0. 1+1.5
c                              do 310 ip=1,nper
c                                 do 310 ia=1,nattn(ip)
c                                    iatype= iatten(ip,ia)
c                                    iatype0= iatype
c                                    if(iatype.gt.4) iatype=3
c                                    ii= dmin(iatype)/di +1.5
c         weight= wt(ip,ia,1)
c         if(dmin(iatype).gt.wtdist(ip,ia)) weight=wt(ip,ia,2)
c                                    do 310 k=1,nlev(ip)
c  310                         prob(icnt,k,ip)= prob(icnt,k,ip)+ 
c     +                        weight*rate(1)*
c     +                        p(ii,mm(1),ip,k,iftype(ift),iatype0)
c                              go to 320
c                           endif
  320          continue
 
c*****write output for 250 receivers************************
c                  write(6,*) i,prob(icnt,1,1)
                  iend= 250
                  if(i.eq.nrec) iend=icnt
                  do 340 ip=1,nper
                     ndata= iend*nlev(ip)
                     do 330 ii=1,ndata
                        i2= ii-1
                        irec= i2/nlev(ip)
                        ilev= i2-irec*nlev(ip)
                        irec= irec+1
                        ilev= ilev+1
                        out(ii)= prob(irec,ilev,ip)
c                       write(*,*) ip,irec,ilev,out(ii)
  330                continue
c                     write(*,*) "hi"
                     call putbufx(ifp(ip),out,ndata,readn)
c                     write(*,*) "bye"
c                    readn= readn/4
c                     write(6,*) ndata,readn
  340             continue
                  icnt= 0
                  do 350 ii=1,250
                     do 350 ip=1,nper
                        do 350 k=1,nlev(ip)
  350             prob(ii,k,ip)= 0.
                  dum=dtime(tarray)
c                  write(*,*) 'Finished 250 receivers'
                  write(*,*) i, tarray
  360       continue
  370       continue
            dum=dtime(tarray)
            write(*,*) 'Finished loops over sources and receivers'
            write(*,*) tarray
            do 371 ip=1,nper
  371       close(ip)
  380 format(a256)
            end
c
 
c      gets corresponding log PSA values (cm/s**2) from interpolation
 
c      of data file ab94.tab.
 
c      Written by G.M. Atkinson, Jan. 1994
 
C----  modified by A. Frankel 10/95
      subroutine gettab(gma,f,xmag,rlog)
      dimension xmag(20),gma(20,30,20), rlog(30), f(20)
      character*60 header
 
c
c     Read AB94 simulated grd. motion table.
 
c
 
      open(3,file='ab94.tab',status='old')
 
      do 10 j = 1, 5
 
         read(3,60)header
         read(3,50) idum
 
 
 
   10 continue
 
 
      nf = 11
 
      nmag = 14
 
      nd = 18
 
      jdl = 0
 
      do 20 jf = 1, nf-2
 
         f(jf) = 10.**(-0. 3+ float(jf-1)*0.2)
 
   20 continue
 
      f(nf) = 89.
 
      f(nf-1) = 99.
 
c
 
      do 40 jm = 1, nmag
 
         read(3,50) idum
         read(3,50) idum
         read(3,70)xmag(jm)
 
         read(3,50) idum
 
         read(3,80)
         read(3,50) idum
         read(3,80)
         read(3,50) idum
 
         do 30 jd = 1, nd
 
            read(3,80)rlog(jd),(gma(jf,jd,jm), jf = 1,nf)
 
            read(3,50) idum
 
   30    continue
 
   40 continue
      close(3)
      write(6,*) gma(1,1,1)
      return
   50 format(i2)
   60 format(a)
   70 format(13x,f5.2)
   80 format(f6.2,2x,11f6.2)
      end
 
c***************************************************************
c*****************************************************************
 
      subroutine table(gma,amag,R,freq,amean,f,xmag,rlog)
 
      dimension xmag(20),gma(20,30,20), rlog(30), f(20)
 
      dimension avg(20)
 
c
 
c     Use interpolation to get amean for given freq(jfreq), amag, R.
 
c
      nf= 11
      nmag= 14
      nd= 18
 
      rl = alog10(R)
 
      jfl = 0
 
c     Find bounding jf values.
 
      do 10 j = 1,nf-3
 
         if (freq .ge. f(j) .and. freq .lt. f(j+1))jfl = j
 
   10 continue
 
      if (abs(freq-f(1)) .lt. 0.01) jfl = 1
 
      if (abs(freq-f(nf)) .lt. 0.01) jfl = nf
 
      if (abs(freq-f(nf-1)) .lt. 0.01) jfl = nf-1
 
      if (abs(freq-20.) .lt. 0.01)  jfl = nf-2
 
      if (jfl .eq. 0) write(*,*) ' ERROR. FREQ. OUTSIDE RANGE.'
 
      if (freq .gt. 20. .and. freq .lt. 89.) write(*,*)
 
     +' ERROR. FREQ. OUTSIDE RANGE'
 
      jfu = jfl + 1
 
      fracf = 0.
 
      if (freq .ge. 20.) go to 20
 
      fracf = (alog10 (freq) - alog10 (f(jfl)) ) / (alog10 (f(jfu)) -
 
     +alog10 (f(jfl)) )
 
   20 continue
 
c     Find bounding jm values.
 
      jml = 0
 
      do 30 j = 1, nmag-1
 
         if (amag .ge. xmag(j) .and. amag .lt. xmag(j+1)) jml = j
 
   30 continue
 
      if (amag .eq. xmag(nmag)) jml = nmag
 
      if (jml .eq. 0) write(*,*)' ERROR. MAGNITUDE OUTSIDE RANGE.'
 
      jmu = jml + 1
 
      fracm = (amag - xmag(jml)) / (xmag(jmu)-xmag(jml))
 
c     find bounding distance values.
 
      do 40 j = 1, nd-1
 
         if (rl .ge. rlog(j) .and. rl .lt. rlog(j+1)) jdl= j
 
   40 continue
 
      if (rl .eq. rlog(nd)) jdl = nd
 
      if (jdl .eq. 0) write(*,*)' ERROR. DISTANCE OUTSIDE RANGE.'
 
      jdu = jdl + 1
 
      fracd = (rl - rlog(jdl)) / (rlog(jdu)-rlog(jdl))
 
      do 50 j = jfl, jfu
 
         arl = gma(j,jdl,jml) + fracm * (gma(j,jdl,jmu)-gma(j,jdl,jml))
 
         aru = gma(j,jdu,jml) + fracm * (gma(j,jdu,jmu)-gma(j,jdu,jml))
 
c       arl is the amplitude for lower dist. bound, so arl>aru
 
         avg(j) = arl - fracd*(arl-aru)
 
   50 continue
 
      amean = avg(jfl) + fracf * (avg(jfu) - avg(jfl))
 
      return
 
      end
 
 
      subroutine readmodel(iflt,tlen)
c     yf,xf,zf(iflt,l,i) are the coordinates (latitude, longitude
c         and depth) of the i-th point on the l-th level
c         (top, hinge, bottom) of the iflt-th fault
      common /blk0/ xf(1,5,1000),yf(1,5,1000),zf(1,5,1000), npts
     +(10,5),dist(5),xfhld(5),yfhld(5),zfhld(5), shold(5), v0x,v0y, vnx,
     +vny, t0x(10,2),t0y(10,2),t0z(10,2), tnx(10,2),tny(10,2),tnz(10,2),
     +t0bx(10),t0by(10),t0bz(10), tnbx(10),tnby(10),tnbz(10)
 
      common /blk1/ rx(1,5,5000),ry(1,5,5000),rz(1,5,5000), rsx
     +(1,5,5000),rsy(1,5,5000),rsz(1,5,5000), send(10,5), vx
     +(1,3,5000),vy(1,3,5000), vz(1,3,5000)
 
 
      dimension tlen(500)
 
 
      parameter (coef=3.14159/180.,nlevels=2)
 
      external fltdst
      tol = 1.e-3
 
 
      do 20 l=1,nlevels
         read(1,*) npts(iflt,l)
         write(*,*) npts(iflt,l)
         do 10 i=1,npts(iflt,l)
            read(1,*) yf(iflt,l,i),xf(iflt,l,i),zf(iflt,l,i)
c             write(*,*) yf(iflt,l,i),xf(iflt,l,i),zf(iflt,l,i)
   10    continue
 
c         set up and resample
 
         write(*,*) 'resample'
         call resample(iflt,l,npts,xf,yf,zf)
   20 continue
 
      write(*,*) 'normvec'
 
      call normvec2(nlevels,tol)
c      tlen(iflt) = send(iflt,1)
       tlen(iflt)= send(iflt,1)
 
      return
      end
 
 
 
      subroutine rupturesetup(imag,jrup,s1,s2)
c
      type ruptparam_type
          sequence
          real s1
          real s2
          integer ist1
          integer ist2
          real xt1
          real yt1
          real zt1
          real xt2
          real yt2
          real zt2
          real sb1
          real sb2
          integer isb1
          integer isb2
          real xb1
          real yb1
          real zb1
          real xb2
          real yb2
          real zb2
          real v0x
          real v0y
          real vnx
          real vny
      end type ruptparam_type
      type(ruptparam_type) rupture(50,1001)
c
      common /blk0/ xf(1,5,1000),yf(1,5,1000),zf(1,5,1000), npts
     +(10,5),dist(5),xfhld(5),yfhld(5),zfhld(5), shold(5), v0x,v0y, vnx,
     +vny, t0x(10,2),t0y(10,2),t0z(10,2), tnx(10,2),tny(10,2),tnz(10,2),
     +t0bx(10),t0by(10),t0bz(10), tnbx(10),tnby(10),tnbz(10)
 
      common /blk1/ rx(1,5,5000),ry(1,5,5000),rz(1,5,5000), rsx
     +(1,5,5000),rsy(1,5,5000),rsz(1,5,5000), send(10,5), vx
     +(1,3,5000),vy(1,3,5000), vz(1,3,5000)
 
 
 
      common /blk2/ x,y,z,xflt,yflt,zflt,iflt,l
 
c     common /blk4/ rupture(50,1001)
      common /blk4/ rupture
 
      parameter (coef=3.14159/180.,nlevels=2)
 
      logical above
 
      external fltdst
 
c      write(*,*) 'rupturesetup'
      tol = 1.e-3
 
 
c     Get x,y,z coordinates of corners of rupture area corresponding to
c         s1,s2
c     s1,s2 are the distances along the top level of the rupture
      ist1 = int(s1) + 2
      ds1 = s1 - int(s1)
      xt1 = rx(iflt,1,ist1) + rsx(iflt,1,ist1)*ds1
      yt1 = ry(iflt,1,ist1) + rsy(iflt,1,ist1)*ds1
      zt1 = rz(iflt,1,ist1) + rsz(iflt,1,ist1)*ds1
      ist2 = int(s2) + 2
      ds2 = s2 - int(s2)
      xt2 = rx(iflt,1,ist2) + rsx(iflt,1,ist2)*ds2
      yt2 = ry(iflt,1,ist2) + rsy(iflt,1,ist2)*ds2
      zt2 = rz(iflt,1,ist2) + rsz(iflt,1,ist2)*ds2
 
c     sb1,sb2 are the distances along the bottom level of the rupture
      sb1 = s1*send(iflt,2)/send(iflt,1)
      sb2 = s2*send(iflt,2)/send(iflt,1)
      isb1 = int(sb1) + 2
      dsb1 = sb1 - int (sb1)
      xb1 = rx(iflt,2,isb1) + rsx(iflt,2,isb1)*dsb1
      yb1 = ry(iflt,2,isb1) + rsy(iflt,2,isb1)*dsb1
      zb1 = rz(iflt,2,isb1) + rsz(iflt,2,isb1)*dsb1
      isb2 = int(sb2) + 2
      dsb2 = sb2 - int(sb2)
      xb2 = rx(iflt,2,isb2) + rsx(iflt,2,isb2)*dsb2
      yb2 = ry(iflt,2,isb2) + rsy(iflt,2,isb2)*dsb2
      zb2 = rz(iflt,2,isb2) + rsz(iflt,2,isb2)*dsb2
 
 
c      write(*,*) xt1,yt1,zt1,xb1,yb1,zb1,xm1,ym1,zm1
c      write(*,*) xt2,yt2,zt2,xb2,yb2,zb2,xm2,ym2,zm2
 
c     Calculate inward pointing vectors perpendicular to end segments
 
 
      cosfac = cos(((yt1+yb1)/2.)*coef)
      v0x =  yb1-yt1
      v0y = -(xb1-xt1)*cosfac
      signfac = 1.
      signfac = sign(signfac, v0x*(xt2-xt1)*cosfac + v0y*(yt2-yt1) )
 
      vmag = sqrt(v0x**2 + v0y**2)
      v0x = v0x * signfac/vmag
      v0y = v0y * signfac/vmag
c      write(*,*) 'signfac,v0x,v0y=',
c     1            signfac,v0x,v0y
 
      cosfac = cos(((yt2+yb2)/2.)*coef)
      vnx =  yb2-yt2
      vny = -(xb2-xt2)*cosfac
      signfac = 1.
      signfac = sign(signfac, vnx*(xt1-xt2)*cosfac + vny*(yt1-yt2) )
 
      vmag = sqrt(vnx**2 + vny**2)
      vnx = vnx * signfac/vmag
      vny = vny * signfac/vmag
c      write(*,*) 'signfac,vnx,vny=',
c     1            signfac,vnx,vny
 
      rupture(imag,jrup)%s1 = s1
      rupture(imag,jrup)%s2 = s2
      rupture(imag,jrup)%ist1 = ist1
      rupture(imag,jrup)%ist2 = ist2
      rupture(imag,jrup)%xt1 = xt1
      rupture(imag,jrup)%yt1 = yt1
      rupture(imag,jrup)%zt1 = zt1
      rupture(imag,jrup)%xt2 = xt2
      rupture(imag,jrup)%yt2 = yt2
      rupture(imag,jrup)%zt2 = zt2
      rupture(imag,jrup)%sb1 = sb1
      rupture(imag,jrup)%sb2 = sb2
      rupture(imag,jrup)%isb1 = isb1
      rupture(imag,jrup)%isb2 = isb2
      rupture(imag,jrup)%xb1 = xb1
      rupture(imag,jrup)%yb1 = yb1
      rupture(imag,jrup)%zb1 = zb1
      rupture(imag,jrup)%xb2 = xb2
      rupture(imag,jrup)%yb2 = yb2
      rupture(imag,jrup)%zb2 = zb2
      rupture(imag,jrup)%v0x = v0x
      rupture(imag,jrup)%v0y = v0y
      rupture(imag,jrup)%vnx = vnx
      rupture(imag,jrup)%vny = vny
      return
      end
      subroutine qkdist(imag,jrup,xpt,ypt,distmin,fltdepmin,distminxy)
c
      type ruptparam_type
          sequence
          real s1
          real s2
          integer ist1
          integer ist2
          real xt1
          real yt1
          real zt1
          real xt2
          real yt2
          real zt2
          real sb1
          real sb2
          integer isb1
          integer isb2
          real xb1
          real yb1
          real zb1
          real xb2
          real yb2
          real zb2
          real v0x
          real v0y
          real vnx
          real vny
      end type ruptparam_type
      type(ruptparam_type) rupture(50,1001)
      common /blk0/ xf(1,5,1000),yf(1,5,1000),zf(1,5,1000), npts
     +(10,5),dist(5),xfhld(5),yfhld(5),zfhld(5), shold(5), v0x,v0y, vnx,
     +vny, t0x(10,2),t0y(10,2),t0z(10,2), tnx(10,2),tny(10,2),tnz(10,2),
     +t0bx(10),t0by(10),t0bz(10), tnbx(10),tnby(10),tnbz(10)
  
      common /blk1/ rx(1,5,5000),ry(1,5,5000),rz(1,5,5000), rsx
     +(1,5,5000),rsy(1,5,5000),rsz(1,5,5000), send(10,5), vx
     +(1,3,5000),vy(1,3,5000), vz(1,3,5000)
  
      common /blk2/ x,y,z,xflt,yflt,zflt,iflt,l
 
c     common /blk4/ rupture(50,1001)
      common /blk4/ rupture
      parameter (coef=3.14159/180.,nlevels=2)
 
      logical above, offend
 
      external fltdst
 
c      write(*,*) 'qkdist'
      tol = 1.e-3
      x = xpt
      y = ypt
      z = 0.
 
      s1 = rupture(imag,jrup)%s1 
      s2 = rupture(imag,jrup)%s2 
      ist1 = rupture(imag,jrup)%ist1  
      ist2 = rupture(imag,jrup)%ist2  
      xt1 = rupture(imag,jrup)%xt1 
      yt1 = rupture(imag,jrup)%yt1 
      zt1 = rupture(imag,jrup)%zt1  
      xt2 = rupture(imag,jrup)%xt2 
      yt2 = rupture(imag,jrup)%yt2 
      zt2 = rupture(imag,jrup)%zt2  
      sb1 = rupture(imag,jrup)%sb1  
      sb2 = rupture(imag,jrup)%sb2  
      isb1 = rupture(imag,jrup)%isb1 
      isb2 = rupture(imag,jrup)%isb2 
      xb1 = rupture(imag,jrup)%xb1  
      yb1 = rupture(imag,jrup)%yb1  
      zb1 = rupture(imag,jrup)%zb1  
      xb2 = rupture(imag,jrup)%xb2  
      yb2 = rupture(imag,jrup)%yb2  
      zb2 = rupture(imag,jrup)%zb2  
      v0x = rupture(imag,jrup)%v0x  
      v0y = rupture(imag,jrup)%v0y  
      vnx = rupture(imag,jrup)%vnx 
      vny = rupture(imag,jrup)%vny
c     Begin the analysis of where x,y,z is relative to the fault
 
c      write(*,*) 'x,y,s12,s2',x,y,s1,s2
 
c     The logical variable "above" will indicate whether the observation point
c     is directly above the the rupture surface, initially set to true
      above =.true.
      offend = .false.
 
c     Find the closest distance and coordinates of the closest points
c            of the fault level boundaries
 
      fltmindist = 1.e10
      do 20 l=1,nlevels 
 
         if(l.eq.1.) then
            as = s1
            cs = s2
         elseif(l.eq.2) then
            as = sb1
            cs = sb2
         else
            as = 0.
            cs = send(iflt,l)
         endif
         call mnbrak(as,bs,cs,fa,fb,fc,fltdst,tol)
c         write(*,*) 'as,bs,cs,fa,fb,fc  ',as,bs,cs,fa,fb,fc
         if(bs.eq.cs) then
            smin=cs
            dist2=fltdst(smin)
         else
            dist2=golden(as,bs,cs,fltdst,tol,smin)
c            write(*,*) 'as,bs,cs,smin  ',as,bs,cs,smin
c             write(*,*) 'yflt,xflt  ', yflt,xflt
         endif
         dist(l)=sqrt(dist2)
c         write(*,*) 'l,smin,dist(l) = ',l,smin,dist(l),yflt,xflt
         shold(l) = smin
         xfhld(l) = xflt
         yfhld(l) = yflt
         zfhld(l) = zflt
         if(dist(l).lt.fltmindist) then
            fltmindist = dist(l)
            ifltmin = l
         endif
c         write(*,*) 'l,yfhld(l)',l,yfhld(l)
   20 continue
 
c     Where is x,y,z relative to the top fault boundary?
      cosfac=cos((yfhld(1) + y)*coef/2.)
      xdiff1 = (x-xfhld(1))*cosfac*111.195
      ydiff1 = (y-yfhld(1))*111.195
      zdiff1 =  z-zfhld(1)
      is = int(shold(1))+1
      if(shold(1).eq.send(iflt,1)) is=is+1
      dotprod1 = xdiff1 * vx(iflt,1,is) + ydiff1 * vy(iflt,1,is) +
     +zdiff1 * vz(iflt,1,is)
 
c      write(*,*) 'dotprod1=',dotprod1
      dotprod01 = xdiff1 *vx(iflt,1,is) + ydiff1 * vy(iflt,1,is)
 
      if(dotprod01.le.0.) above = .false.
      if(dotprod1.lt.0) then
c        point on "updip" side of top fault boundary
         ireg=1
         distmin=dist(1)
         i fltdepmin = zfhld(1)
 
      else
c     Where is x,y,z relative to the bottom fault boundary?
         cosfac=cos((yfhld(2) + y)*coef/2.)
         xdiff2 = (x - xfhld(2))*cosfac*111.195
         ydiff2 = (y - yfhld(2))*111.195
         zdiff2 = z - zfhld(2)
         is = int(shold(2))+1
         if(shold(2).eq.send(iflt,2)) is=is+1
         dotprod2 = xdiff2 * vx(iflt,2,is) + ydiff2 * vy(iflt,2,is) +
     +   zdiff2 * vz(iflt,2,is)
 
c          write(*,*) 'dotprod2=', dotprod2
         dotprod02 = xdiff2 * vx(iflt,2,is) + ydiff2 * vy(iflt,2,is)
 
         if(dotprod02.ge.0.) above = .false.
         if(dotprod2.gt.0.) then
c             point on "downdip" side of bottom fault boundary
            ireg=2
            distmin=dist(2)
            fltdepmin=zfhld(2)
         else
            cosfac=cos((yt1 + y)*coef/2.)
            xdiff4 = (x-xt1)*cosfac*111.195
            ydiff4 = (y-yt1)*111.195
 
            dotprod4 = v0x*xdiff4 + v0y*ydiff4
            cosfac=cos((yt2 + y)*coef/2.)
            xdiff5 = (x-xt2)*cosfac*111.195
            ydiff5 = (y-yt2)*111.195
 
            dotprod5 = vnx*xdiff5 + vny*ydiff5
 
            if(dotprod5.lt.0.) above = .false.
            if(dotprod4.lt.0.) above = .false.
 
               if(dotprod4.le.0.) then
                  ireg = 4
               else
                  if(dotprod5.le.0.) then
                     ireg = 5
                  else
                     ireg = 3
                     cosfac=cos((yfhld(1) + yfhld(2))*coef/2.)
                     dxlev1 = (xfhld(2) - xfhld(1))*cosfac*111.195
                     dylev1 = (yfhld(2) - yfhld(1))*111.195
                     dzlev1 = zfhld(2) - zfhld(1)
                     isplus = int(shold(2)) + 5
                     if(isplus.gt.send(iflt,2)) isplus = isplus -10
                     dxlev2 = (rx(iflt,2,isplus) - xfhld(1))*cosfac
     +               *111.195
                     dylev2 = (ry(iflt,2,isplus) - yfhld(1))*111.195
                     dzlev2 = rz(iflt,2,isplus) - zfhld(1)
c                        write(*,*) x,y,z,xfhld(2),xfhld(3),
c     1                                   yfhld(2),yfhld(3)
 
c                        write(*,*) dxlev1,dylev1,dzlev1,
c     1                             dxlev2,dylev2,dzlev2,
c     2                             xdiff2,ydiff2,zdiff2
                     dxcros = dylev1*dzlev2-dzlev1*dylev2
                     dycros = -(dxlev1*dzlev2-dzlev1*dxlev2)
                     dzcros = dxlev1*dylev2-dylev1*dxlev2
c                        write(*,*) xdiff2,ydiff2,zdiff2,dxcros,dycros,
c     1                             dzcros
c                        write(*,*) dxlev1,dylev1,dzlev1,
c     1                             dxlev2,dylev2,dzlev2
                     crosmag = sqrt(dxcros**2+dycros**2+dzcros**2)
                     distmin = abs ( (dxcros*xdiff2 + dycros*ydiff2 +
     +               dzcros*zdiff2)) / crosmag
 
                     fltdepmin = distmin * dzcros / crosmag
 
                  endif
               endif
 
c     Calculate distance for off end points
            if(ireg.eq.4) then
               offend = .true.
               cosfac = cos((y + yt1)*coef/2.)
               ax = (x - xt1)*cosfac*111.195
               ay = (y - yt1)*111.195
               az = - zt1
               cosfac = cos((yb1 + yt1)*coef/2.)
               bx = (xb1 - xt1)*cosfac*111.195
               by = (yb1 - yt1)*111.195
               bz = zb1 - zt1
               bmag = sqrt(bx*bx + by*by + bz*bz)
               cx = ay*bz - az*by
               cy = -(ax*bz - az*bx)
               cz = ax*by - ay*bx
               distmin = sqrt(cx*cx + cy*cy + cz*cz)/bmag
               fltdepmin = zt1 + ( ax*bx+ay*by+az*bz)*(zb1-zt1)
     +         /(bmag*bmag)
               distminxy = abs(dotprod4) 
 
 
            elseif(ireg.eq.5) then
               offend = .true.
               cosfac = cos((y + yt2)*coef/2.)
               ax = (x - xt2)*cosfac*111.195
               ay = (y - yt2)*111.195
               az = - zt2
               cosfac = cos((yb2 + yt2)*coef/2.)
               bx = (xb2 - xt2)*cosfac*111.195
               by = (yb2 - yt2)*111.195
               bz = zb2 - zt2
               bmag = sqrt(bx*bx + by*by + bz*bz)
               cx = ay*bz - az*by
               cy = -(ax*bz - az*bx)
               cz = ax*by - ay*bx
c                  write(*,*) ax,ay,az,bx,by,bz,icx,cy,cz,bmag
               distmin = sqrt(cx*cx + cy*cy + cz*cz)/bmag
               fltdepmin = zt2 + ( ax*bx+ay*by+az*bz)*(zb2-zt2)
     +         /(bmag*bmag)
               distminxy = abs(dotprod5)
 
  
            endif
 
 
         endif
 
 
 
      endif
c     Final checks
c     Check for wierdness off ends
      if(fltmindist.lt.distmin) then
         distmin = fltmindist
         if(ifltmin.eq.1) ireg=1
         if(ifltmin.eq.2) ireg=2
      endif
 
      if(above.and.(ireg.eq.1.or.ireg.eq.2)) then
         cosfac=cos((yb1 + y)*coef/2.)
         xdiff4 = (x-xb1)*cosfac*111.195
         ydiff4 = (y-yb1)*111.195
 
         dotprod4 = v0x*xdiff4 + v0y*ydiff4
         cosfac=cos((yb2 + y)*coef/2.)
         xdiff5 = (x-xb2)*cosfac*111.195
         ydiff5 = (y-yb2)*111.195
 
         dotprod5 = vnx*xdiff5 + vny*ydiff5
 
         if(dotprod5.lt.0.) above = .false.
         if(dotprod4.lt.0.) above = .false.
      endif
      if(above) then
         distminxy = 0.
      else
         if(.not.offend)
     +   distminxy = sqrt(((x-xfhld(ifltmin))*cosfac)** 2+ (y-yfhld
     +   (ifltmin))**2) * 111.195
      endif
c      write(*,*) 'ireg,distmin,distminxy =', ireg,distmin,distminxy
c      write(9,*) y,x,ireg,distmin
 
   30 continue
c     goto 30
   40 return
 
      end
 
 
 
 
      subroutine resample(iflt,l,npts,xf,yf,zf)
 
      dimension xf(1,5,1000),yf(1,5,1000),zf(1,5,1000)
c     yf,xf,zf(iflt,l,i) are the coordinates (latitude, longitude
c         and depth) of the i-th point on the l-th level
c         (top, hinge, bottom) of the iflt-th fault
 
      dimension npts(10,5)
 
      dimension s(5000),xs(5000),ys(5000),zs(5000)
 
      common /blk1/ rx(1,5,5000),ry(1,5,5000),rz(1,5,5000), rsx
     +(1,5,5000),rsy(1,5,5000),rsz(1,5,5000), send(10,5), vx
     +(1,3,5000),vy(1,3,5000), vz(1,3,5000)
 
 
 
 
      parameter (coef=3.14159/180.)
 
      s(1)= 0.
      do 10 i=2,npts(iflt,l)
         ydiff = yf(iflt,l,i)-yf(iflt,l,i-1)
         cosfac=cos((yf(iflt,l,i-1)+.5*ydiff)*coef)
 
         s(i) = s(i-1) + 111.195* sqrt((cosfac*(xf(iflt,l,i)-xf
     +   (iflt,l,i-1)))** 2+ ydiff**2)
 
         xs(i) = (xf(iflt,l,i)-xf(iflt,l,i-1))/ (s(i)-s(i-1))
 
         ys(i) = (yf(iflt,l,i)-yf(iflt,l,i-1))/ (s(i)-s(i-1))
 
         zs(i) = (zf(iflt,l,i)-zf(iflt,l,i-1))/ (s(i)-s(i-1))
 
c          write(*,*) i,xf(iflt,l,i),yf(iflt,l,i),s(i),xs(i),ys(i)
 
   10 continue
      send(iflt,l) = s(npts(iflt,l))
 
c     Resample faults at points 1 km along trace
      rx(iflt,l,1)= xf(iflt,l,1)
      ry(iflt,l,1)= yf(iflt,l,1)
      rz(iflt,l,1)= zf(iflt,l,1)
      intv=2
      nptsr = int(s(npts(iflt,l))) + 2
c      write(*,*) 's(npts(iflt,l)),nptsr= ',
c     1            s(npts(iflt,l)),nptsr
      do 40 i=2,nptsr-1
         s1=i-1
         do 20 j=intv,npts(iflt,l)
            if(s1.le.s(j)) then
               intv=j
               goto 30
            endif
   20    continue
   30    ds= s1-s(intv-1)
         rx(iflt,l,i)= xf(iflt,l,intv-1) + ds*xs(intv)
         ry(iflt,l,i)= yf(iflt,l,intv-1) + ds*ys(intv)
         rz(iflt,l,i)= zf(iflt,l,intv-1) + ds*zs(intv)
         rsx(iflt,l,i)= rx(iflt,l,i) - rx(iflt,l,i-1)
         rsy(iflt,l,i)= ry(iflt,l,i) - ry(iflt,l,i-1)
         rsz(iflt,l,i)= rz(iflt,l,i) - rz(iflt,l,i-1)
c        write(*,*) s1,intv,xf(iflt,l,intv-1),xs(intv),
c     1             yf(iflt,l,intv-1),
c     1          ys(intv),rx(iflt,l,i),ry(iflt,l,i)
   40 continue
 
      rx(iflt,l,nptsr)=xf(iflt,l,npts(iflt,l))
      ry(iflt,l,nptsr)=yf(iflt,l,npts(iflt,l))
      rz(iflt,l,nptsr)=zf(iflt,l,npts(iflt,l))
      ds= s(npts(iflt,l)) - aint(s(npts(iflt,l)))
      rsx(iflt,l,nptsr)= (rx(iflt,l,nptsr)- rx(iflt,l,nptsr-1))/ds
 
      rsy(iflt,l,nptsr)= (ry(iflt,l,nptsr)- ry(iflt,l,nptsr-1))/ds
 
      rsz(iflt,l,nptsr)= (rz(iflt,l,nptsr)- rz(iflt,l,nptsr-1))/ds
 
      return
      end
 
 
 
 
      function fltdst(s)
c     On return fltdst is the squared distance from x,y,z to
c     a point on the fault, xflt(s),yflt(s),zflt(s)
c     s is arc length along the fault
 
      common /blk1/ rx(1,5,5000),ry(1,5,5000),rz(1,5,5000), rsx
     +(1,5,5000),rsy(1,5,5000),rsz(1,5,5000), send(10,5), vx
     +(1,3,5000),vy(1,3,5000), vz(1,3,5000)
 
 
 
      common /blk2/ x,y,z,xflt,yflt,zflt,iflt,l
      parameter (coef=3.14159/180.,distsq=111.195**2)
 
 
c      write(*,*) 's= ',s
      intv=int(s)+2
   10 ds=s-aint(s)
c      write(*,*) intv,ds
      xflt = rx (iflt,l,intv-1) + ds * rsx(iflt,l,intv)
      yflt = ry (iflt,l,intv-1) + ds * rsy(iflt,l,intv)
      zflt = rz (iflt,l,intv-1) + ds * rsz(iflt,l,intv)
 
      ydiff = y-yflt
      cosfac = cos((yflt+.5*ydiff)*coef)
      fltdst = (((x-xflt)*cosfac)**2 + ydiff**2 )*distsq + (z-zflt)**2
c      write(*,*) 'xflt,yflt,zflt,flddst= ', xflt,yflt,zflt,fltdst
      return
      end
 
 
      subroutine normvec2(nlevel,tol)
c     Find unit vectors normal to top of fault pointing in down
c          dip direction
c     At ends of fault, make vector parallel to end segment
      common /blk1/ rx(1,5,5000),ry(1,5,5000),rz(1,5,5000), rsx
     +(1,5,5000),rsy(1,5,5000),rsz(1,5,5000), send(10,5), vx
     +(1,3,5000),vy(1,3,5000), vz(1,3,5000)
 
 
      common /blk2/ x,y,z,xflt,yflt,zflt,iflt,lf
 
      parameter (coef=3.14159/180.)
 
      external fltdst
 
      iflt = 1
 
c     lv will refer to the level along which you seek the normal vectors
c     lf will refer to other level used to define the dip (i.e.the level
c          below if lv=1 or 2, above if lv =3)
 
 
      do 20 lv=1,nlevel
         if(lv.lt.nlevel) then
            lf = lv + 1
            signfac = 1.
            ls = 1
         else
            lf = lv - 1
            signfac = -1.
            ls = 2
         endif
c       write(*,*) iflt,int(send(iflt,lv))+2
         do 10 i=1,int(send(iflt,lv))+2
 
            x = rx(iflt,lv,i)
            y = ry(iflt,lv,i)
            z = rz(iflt,lv,i)
 
            as = 0.
            cs = send(iflt,lf)
 
c         write(*,*) lv,lf
c         write(*,*) 'x,y,z', x,y,z
            call mnbrak(as,bs,cs,fa,fb,fc,fltdst,tol)
c         write(*,*) 'as,bs,cs,fa,fb,fc  ',as,bs,cs,fa,fb,fc
            if(bs.eq.cs) then
               smin=cs
               dist2=fltdst(smin)
            else
               dist2=golden(as,bs,cs,fltdst,tol,smin)
c            write(*,*) 'as,bs,cs,smin  ',as,bs,cs,smin
c            write(*,*) 'yflt,xflt  ', yflt,xflt
            endif
c         write(*,*) 'x,y,z,xflt,yflt,zflt',x,y,z,xflt,yflt,zflt
            cosfac = cos(ry(iflt,lv,i)*coef)
 
            vx(iflt,lv,i) = signfac * (xflt - x) *cosfac * 111.195
            vy(iflt,lv,i) = signfac * (yflt - y) * 111.195
            vz(iflt,lv,i) = signfac * (zflt - z)
            vmag =sqrt(vx(iflt,lv,i)** 2+ vy(iflt,lv,i)** 2+ vz
     +      (iflt,lv,i)**2)
c     1                  iflt,l,i,rx(iflt,l,i),rx(iflt,l,i-1),
c     2                  ry(iflt,l,i),ry(iflt,l,i-1),
c     3                  vx(iflt,l,i),vy(iflt,l,i)
            vx(iflt,lv,i) = vx(iflt,lv,i)/vmag
            vy(iflt,lv,i) = vy(iflt,lv,i)/vmag
            vz(iflt,lv,i) = vz(iflt,lv,i)/vmag
 
 
 
c           write(*,*) 'vx(iflt,lv,i),vy(iflt,lv,i)',
c     1                vx(iflt,lv,i),vy(iflt,lv,i)
 
 
   10    continue
 
 
c      do 55 i=1,int(send(iflt,lv))+2
c         write(*,*) vx(iflt,lv,i),vy(iflt,lv,i),vz(iflt,lv,i)
c 55   continue
 
   20 continue
      return
      end
 
 
 
      function zbrent(func,x1,x2,tol)
      parameter (itmax=100,eps=3.e-8)
c      write(*,*) "zbrent"
      a=x1
      b=x2
c      write(*,*) "a,b", a,b
      fa=func(a)
      fb=func(b)
c      write(*,*) "fa,fb", fa,fb
      if(fb*fa.gt.0.) write(*,*) 'Root must be bracketed for zbrent.'
      fc=fb
      do 10 iter=1,itmax
         if(fb*fc.gt.0.) then
            c=a
            fc=fa
            d=b-a
            e=d
         endif
         if(abs(fc).lt.abs(fb)) then
            a=b
            b=c
            c=a
            fa=fb
            fb=fc
            fc=fa
         endif
         tol1=2.*eps*abs(b)+0.5*tol
         xm=.5*(c-b)
         if(abs(xm).le.tol1 .or. fb.eq.0.) then
            zbrent=b
            return
         endif
         if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
            s=fb/fa
            if(a.eq.c) then
               p=2.*xm*s
               q=1.-s
            else
               q=fa/fc
               r=fb/fc
               p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.))
               q=(q-1.)*(r-1.)*(s-1.)
            endif
            if(p.gt.0.) q=-q
            p=abs(p)
            if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then
               e=d
               d=p/q
            else
               d=xm
               e=d
            endif
         else
            d=xm
            e=d
         endif
         a=b
         fa=fb
         if(abs(d) .gt. tol1) then
            b=b+d
         else
            b=b+sign(tol1,xm)
         endif
         fb=func(b)
   10 continue
      write(*,*) 'zbrent exceeding maximum iterations.'
      zbrent=b
      return
      end
 
  
      SUBROUTINE MNBRAK(AX,BX,CX,FA,FB,FC,FUNC,tol)
c     Modified from Numerical Recipes
c     Want to establish whether minimum is within interval [ax,cx],
c     or at end point.  If minimum is within interval, then ax,bx,cx
c     will bracket minimum
      PARAMETER (GOLD=0.618034, GLIMIT=100., TINY=1.E-20)
      FA=FUNC(AX)
      FC=FUNC(CX)
      IF(FC.GT.FA)THEN
         DUM=AX
         AX=CX
         CX=DUM
         DUM=FC
         FC=FA
         FA=DUM
      ENDIF
   10 BX=AX+GOLD*(CX-AX)
      if(abs(cx-bx).le.tol) then
         bx=cx
         go to 30
      endif
      FB=FUNC(BX)
   20 IF(FB.LE.FC)THEN
         go to 30
c        Minimum is bracketed by ax,bx,cx
      else
c        Replace ax with bx, and choose new bx
         ax=bx
         fa=fb
         go to 10
      endif
   30 return
      end
 
 
 
 
      FUNCTION GOLDEN(AX,BX,CX,F,TOL,XMIN)
      PARAMETER (R=.61803399,C=.38196602)
      X0=AX
      X3=CX
      IF(ABS(CX-BX).GT.ABS(BX-AX))THEN
         X1=BX
         X2=BX+C*(CX-BX)
      ELSE
         X2=BX
         X1=BX-C*(BX-AX)
      ENDIF
      F1=F(X1)
      F2=F(X2)
   10 IF(ABS(X3-X0).GT.TOL*(ABS(X1)+ABS(X2)))THEN
         IF(F2.LT.F1)THEN
            X0=X1
            X1=X2
            X2=R*X1+C*X3
            F0=F1
            F1=F2
            F2=F(X2)
         ELSE
            X3=X2
            X2=X1
            X1=R*X2+C*X0
            F3=F2
            F2=F1
            F1=F(X1)
         ENDIF
         GOTO 10
      ENDIF
      IF(F1.LT.F2)THEN
         GOLDEN=F1
         XMIN=X1
      ELSE
         GOLDEN=F2
         XMIN=X2
      ENDIF
      RETURN
      END
 
      FUNCTION BRENT(AX,BX,CX,F,TOL,XMIN)
c
c     Numerical Recipes, 1st Edition
c
      PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0E-10)
      A=MIN(AX,CX)
      B=MAX(AX,CX)
      V=BX
      W=V
      X=V
      E=0.
      FX=F(X)
      FV=FX
      FW=FX
      DO 30 ITER=1,ITMAX
         XM=0.5*(A+B)
         TOL1=TOL*ABS(X)+ZEPS
         TOL2=2.*TOL1
         IF(ABS(X-XM).LE.(TOL2-.5*(B-A))) GOTO 40
         IF(ABS(E).GT.TOL1) THEN
            R=(X-W)*(FX-FV)
            Q=(X-V)*(FX-FW)
            P=(X-V)*Q-(X-W)*R
            Q=2.*(Q-R)
            IF(Q.GT.0.) P=-P
            Q=ABS(Q)
            ETEMP=E
            E=D
            IF(ABS(P).GE.ABS(.5*Q*ETEMP).OR.P.LE.Q*(A-X).OR. P.GE.Q
     +      *(B-X)) GOTO 10
            D=P/Q
            U=X+D
            IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=SIGN(TOL1,XM-X)
            GOTO 20
         ENDIF
   10    IF(X.GE.XM) THEN
            E=A-X
         ELSE
            E=B-X
         ENDIF
         D=CGOLD*E
   20    IF(ABS(D).GE.TOL1) THEN
            U=X+D
         ELSE
            U=X+SIGN(TOL1,D)
         ENDIF
         FU=F(U)
         IF(FU.LE.FX) THEN
            IF(U.GE.X) THEN
               A=X
            ELSE
               B=X
            ENDIF
            V=W
            FV=FW
            W=X
            FW=FX
            X=U
            FX=FU
         ELSE
            IF(U.LT.X) THEN
               A=U
            ELSE
               B=U
            ENDIF
            IF(FU.LE.FW .OR. W.EQ.X) THEN
               V=W
               FV=FW
               W=U
               FW=FU
            ELSE IF(FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN
               V=U
               FV=FU
            ENDIF
         ENDIF
   30 CONTINUE
      write(*,*) 'Brent exceed maximum iterations.'
   40 XMIN=X
      BRENT=FX
      RETURN
      END
