c--- program hazFXv7; 8/16/05; uses subroutines for attenuation relations 
c----  written by Art Frankel, with mods by Steve Harmsen
c--- includes test for 90 degree dip to speed up code
c---- averages strike of fault segment
c---- revised hanging wall test for AS97, including 22.5 degree wings
c-- can do up to 7 attenuation relations for each period 
c---- Added variable probability along fault 8/16/05 OSB
c---- added getMota 10/16/02
c---- revised Campbell2002
c---- added hanging wall term for Campbell 2002
c---- hanging wall term, separate reverse from thrust faulting
c---- fixed hanging wall for floating rupture
c---- changed to dmove increment for floating rupture 9/30/02
c---- added getBJF97 9/30/02
c--- changed length to WC all 7/1/02
c--- changed to only epistemic uncertainty for GR faults
c----- changed to Mmin= 5.8  nmag= 55
c---- added clamp for CEUS   max gnd motions
c---- added Mchar epistemic uncertainty 10/16/01
c---- added Mchar aleatory uncertainty 10/16/01
c---- checks to see if distribution goes below M6.5 for GR
c---- or below 6.0 for char.
c---- if so, ignores uncertainties
cc----- truncates at 3 sigma
c--- with clipping for Toro and Boore new table
c--- calculates mean annual frequency of exceedance for dipping faulta
c--- written by Art Frankel 4/95, modified 10/95
c--- this version handles multiple periods and gnd motion relations
c--- to compile: cc -c iosubs.c
c                 f77 -o hazFX hazFX.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
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
c---- p matrix is based on either hypocentral distance or
c---   BJF distance for BJF, Toro, and Somerville CEUS
      common / atten / p, phang, xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
      dimension x(500,50),y(500,50),xaz(500,50),xlen(500,50)
      dimension fdist(2000,100,6),nrup(100),dmin2(20,2000,6)
      dimension dmin(6),dminr(6),nmag0(500),torate(55),itest(500)
      dimension ale(100), ale2(100),ratenew(500,50,50),xmo2(500,50)
      dimension xmo2b(2000),xmo2c(20,2000),ab(2000),ac(20,2000)
      dimension ratenewb(2000,50),ratenewc(20,2000,50)
      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(50),magmax(50), mmax, mmax2
      integer*4 readn, ndata
      character nameout*256,name*256,adum*256
      dimension dmbranch(10),wtbranch(10), testhw3(2000,50)
      dimension prob(1000,20,7),out(800000),
     &  iatten(7,7),wt(7,7,2),wtdist(7,7)
      dimension nattn(7),period(7)
      dimension ifpn(7)
      dimension ifp(7)
      dimension tlen(50),a(500,50),b(50),dmag(50)
      dimension itype(50),iftype(50),cmag(50),crate(50)
c     dimension itype(500),iftype(500),cmag(500),crate(500)
      dimension dip(50),dip0(50),width(50),depth0(50),npts(50)
c      write(6,*) "enter name of input file"
 900  format(a30)
      call getarg(1,name)
      open(unit=1,file=name,status='old')
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
      nx= ceiling((xmax-xmin)/dx +1)
      ny= ceiling((ymax-ymin)/dy +1)
      write(6,*) nx,ny
      pi= 3.14159
      nrec= nx*ny
cccccccccccccccccccccccccccccccccc
      write(6,*) "enter distance increment, dmax"
      read(1,*) di,dmax
      write(6,*) "enter number of periods"
      read(1,*) nper
      ndist= dmax/di +10
      nmag= 55
      do 301 ii=1,ndist
      do 301 m=1,nmag
      do 301 ip=1,nper
      do 301 k=1,20
      do 301 ifkey=1,4
      do 301 iatype=1,7
 301  p(ii,m,ip,k,ifkey,iatype)=0.
      angtest= 22.5*3.14159/180.
      angtest= tan(angtest)
c---loop through periods
      do 700 ip=1,nper
      read(1,*) period(ip), sigmanf, distnf
      per= 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,*) ifp(ip)
      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 414 k=1,nlev(ip)
 414  xlev(k,ip)= alog(xlev(k,ip))
c-----------
      write(6,*) "enter number of atten. relations for this period"
      read(1,*) nattn(ip)
c--- loop through atten relations for that period
      do 701 ia=1,nattn(ip)
      write(6,*) "enter type of atten. relation, weight1, wtdist, 
     &  weight2 , mb to M conv."
      read(1,*) iatten(ip,ia),wt(ip,ia,1),wtdist(ip,ia),wt(ip,ia,2),
     &   iconv(ip,ia)
      if(iatten(ip,ia).eq.1) 
     &  call getBJF93(ip,per,ia,ndist,di,nmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.2) 
     &  call getToro(ip,per,ia,ndist,di,nmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.3) 
     & call getSadigh(ip,per,ia,ndist,di,nmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.4) 
     & call getCamp(ip,per,ia,ndist,di,nmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.5) 
     & call getAB95(ip,per,ia,ndist,di,nmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.6) 
     & call getFEA(ip,per,ia,ndist,di,nmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.7) 
     & call getSomer(ip,per,ia,ndist,di,nmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.8) 
     & call getAS(ip,per,ia,ndist,di,nmag,sigmanf,distnf)    
      if(iatten(ip,ia).eq.9)
     & call getCamp2000(ip,per,ia,ndist,di,nmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.10)
     & call getCampCEUS(ip,per,ia,ndist,di,nmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.11) 
     &  call getBJF97(ip,per,ia,ndist,di,nmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.12)
     &  call getMota(ip,per,ia,ndist,di,nmag,sigmanf,distnf)
 701  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
 700  continue
ccccccccccccccccccccccccccccccccccc
      write(6,*) "enter increment dlen and dmove for floating rup"
      read(1,*) dlen, dmove
      write(6,*) dlen, dmove
      write(6,*) "enter number of branches for Mchar logic tree"
      read(1,*) nbranch
      write(6,*) "enter dM of each branch"
      read(1,*) (dmbranch(i),i=1,nbranch)
      write(6,*) "enter weight of each branch"
      read(1,*) (wtbranch(i),i=1,nbranch)
      write(6,*) "enter std dev and width for Mchar aleatory uncer"
      read(1,*) sdal,mwid
      if(sdal.gt.0.) then
        denom=0.
        do 670 i= -mwid, mwid
          ale(i+mwid+1)= exp(-(i*0.05)*(i*0.05)/(2.*sdal*sdal))
          denom= denom+ ale(i+mwid+1)*10.**(1.5*i*0.05)
          write(6,*) ale(i+mwid+1)
 670    continue
        do 671 i=1,2*mwid+1
  671     ale2(i)= ale(i)/denom
        write(6,*) denom
c       read(5,*) idum
      endif
ccc-------read fault data
      do 1 ift=1,1000
        write(6,*) ift
        write(6,*) "enter 1 for char. 2 for G-R; 1 for SS, 2 for rev"
        read(1,*,end=999) itype(ift),iftype(ift)
        dmag0= dmag(ift)
        if(sdal.eq.0.) itest(ift)=1
        if(itype(ift).eq.1) then
          write(6,*) "enter char. M, rate"
          read(1,*) cmag(ift),crate(ift)
          itest(ift)=0
          test= cmag(ift)+dmbranch(1)-mwid*0.05
          if(test.lt.5.8) then
            itest(ift)=1
            write(6,*) "test.lt.6.",ift,itest(ift)
c           read(5,*) idum
          endif
        endif
        if(itype(ift).eq.2) then
          write(6,*) "b,min M,max M,dmag"
          read(1,*) b(ift),magmin(ift),magmax(ift),dmag(ift)
          write(6,*) b(ift),magmin(ift),magmax(ift),dmag(ift)
c---- center magnitude bins
          if(magmin(ift).ne.magmax(ift)) then
            magmin(ift)= magmin(ift)+dmag(ift)/2.
            magmax(ift)= magmax(ift)-dmag(ift)/2.
          endif
          nmag0(ift)= (magmax(ift)-magmin(ift))/dmag(ift) + 1.4
          write(6,*) "*** nmag0, ift ***",nmag0(ift),ift
          itest(ift)=0
          test= magmax(ift)+dmbranch(1)
          if((test.lt.6.5).and.(nmag0(ift).gt.1)) then
            itest(ift)=1
            write(6,*) "test.lt.6.5", ift, itest(ift)
c           read(5,*) idum
          endif
          if(nmag0(ift).eq.1) then
            test= magmax(ift)+dmbranch(1)-mwid*0.05
            if(test.lt.5.8) itest(ift)=1
          endif
        endif
        write(6,*) "enter dip, downdip width, depth0 of fault plane"
        read(1,*) dip(ift), width(ift), depth0(ift)
        write(6,*) "enter number of segment points"
        read(1,*) npts(ift)
        write(6,*) "enter lat,lon for each point"
        tlen(ift)= 0.
        do 2 j=1,npts(ift)
          if(itype(ift).eq.1) then
c           write(6,*) "enter y,x"
            read(1,*) y(j,ift),x(j,ift)
          endif
          if(itype(ift).eq.2) then
c           write(6,*) "enter y,x,a"
            read(1,*) y(j,ift),x(j,ift),a(j,ift)
c---- calculate moment rate
            xmorate=0.
            do 600 i=1,nmag0(ift)
              xmag= magmin(ift)+ dmag(ift)*(i-1)
              xmorate= xmorate+10.**(a(j,ift)-b(ift)*xmag+1.5*xmag+9.05)
 600        continue
            xmo2(j,ift)= xmorate
c           write(6,*) xmorate
c----- find rates for epistemic uncertainty
            if(nmag0(ift).gt.1) then
              do 603 ilt= 1, nbranch
                mmax= magmax(ift)+ dmbranch(ilt)
                sum2= 0.
                nmag= (mmax-magmin(ift))/dmag(ift) + 1.4
                sum= 0.
                do 602 m=1,nmag
                  xmag= magmin(ift)+ (m-1)*dmag(ift)
  602             sum= sum + 10.**(-b(ift)*xmag+1.5*xmag+9.05)
                ratenew(j,ift,ilt) = alog10(xmorate/sum)
                write(6,*) ratenew(j,ift,ilt),a(j,ift)
c               write(6,*) "sum2=",sum2
c               read(5,*) idum
  603         continue
            endif
          endif
c         read(1,*) y(j,ift),x(j,ift)
          if(j.gt.1) then
            sy= y(j-1,ift)
            sx= x(j-1,ift)
            ry= y(j,ift)
            rx= x(j,ift)
            if((sy.eq.ry).and.(sx.eq.rx)) then
              delta=0.
              az=0.
              go to 1020
            else
              call delaz(sy,sx,ry,rx,delta,az,baz)
            endif
 1020       if(j.gt.2) xlen(j-1,ift)= delta +xlen(j-2,ift)
            if(j.eq.2) xlen(1,ift)= delta
            xaz(j-1,ift)= az
            tlen(ift)= tlen(ift)+ delta
          endif
c         write(6,*) xlen(j-1),delta,xaz(j-1)
   2    continue
c----- next line added 10/18 based on Harmsen comment
        xaz(npts(ift),ift)= xaz(npts(ift)-1,ift)
        write(6,*) npts(ift),tlen(ift)
        if(tlen(ift).le.0.) then
          write (6,*) "tlen le 0",ift
          read(5,*) idum
        endif
        coef= 3.14159/180.
        dip0(ift)= dip(ift)
        dip(ift) =dip(ift)*coef
   1  continue
 999  nft= ift-1
      write(6,*) "nft=",nft
ccccccccccccccccccccccccccccc
c---- write header --------------
      do 201 ip=1,nper
        headr%name(1)= name
        headr%period= period(ip)
        headr%nlev= nlev(ip)
        do 702 k=1,nlev(ip)
 702      headr%xlev(k)= exp(xlev(k,ip))
        ndata= 308
        call puthead(ifp(ip),headr,ndata,readn)
c       write(6,900) headr%name(1)
        write(6,*) readn
 201  continue
ccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccc
      icnt=1
      do 30 j=1,1000
        do 30 ip=1,nper
          do 30 k=1,nlev(ip)
  30        prob(j,k,ip)=0.
c---Here's the guts
c
c---loop through receiver sites
      do 100 i=1,nrec
c     do 100 i=50000,nrec
        iy= (i-1)/nx 
        ix= i-1-iy*nx
        rx= xmin+(ix)*dx
        ry= ymax-(iy)*dy
c       if(i.gt.50830) write (6,*) "rec",i
        do 850 ift=1,nft
c         if(i.gt.50830) write (6,*) "ift",ift
c--- loop through faults
          xdiff= abs(rx-x(1,ift))
          ydiff= abs(ry-y(1,ift))
          xdmax= (dmax+1000.)/(111.11*cos(ymax*coef))
          ydmax= (dmax+1000.)/111.11
          if(xdiff.gt.xdmax) go to 850
          if(ydiff.gt.ydmax) go to 850
c      write(6,*) "prob",isz,i,prob(i-1,1)
          do 35 iatype=1,5
 35         dmin(iatype)=10000.
          nlenx= tlen(ift)/dlen +1
          nleny= width(ift)/dlen +1
          if(dip0(ift).eq.90.) nleny=1
c----- new for AS hanging wall test 
          width2= (50.-depth0(ift))/(sin(abs(dip(ift))))
          nleny2= width2/dlen +1
c--------------
          nseg= npts(ift)-1
c--- calculate distance to site for points (ix,iy) on fault plane
c---loop through increments along strike (ix) of each segment
          do 101 ix=1,nlenx
c           if(i.gt.50830) write (6,*) "ix",ix
            delta= (ix-1)*dlen
            do 500 iseg=1,nseg
              if(delta.lt.xlen(iseg,ift)) go to 501
 500        continue
 501        continue
            if(iseg.ge.2) delta= delta-xlen(iseg-1,ift)
            az= xaz(iseg,ift)
            sx= x(iseg,ift)
            sy= y(iseg,ift)
c      write(6,*) "ix, iseg, az", ix, iseg, az
            call back(delta,az,sy,sx,dfy,dfx)
            if(ix.eq.1) azprev= az
            if(abs(azprev-az).gt.10.) az= 0.5*(azprev+az)
            if(dip(ift).lt.0.) az2= az-90.
            if(dip(ift).gt.0.) az2= az+90.
            azprev= az
            if(itype(ift).eq.2) then
              xmo2b(ix) = xmo2(iseg,ift)
              ab(ix) = a(iseg,ift)
              do 703 ilt= 1, nbranch
 703            ratenewb(ix,ilt) = ratenew(iseg,ift,ilt)
            endif
c---loop through increments iy downdip
            do 111 iy=1,nleny
              deltaz= (iy-1)*dlen
              depth= depth0(ift)+ deltaz*sin(abs(dip(ift)))
              deltaz= deltaz*cos(dip(ift))
c---find surface point above downdip portion of fault
              call back(deltaz,az2,dfy,dfx,dfy2,dfx2)
c--- find distance from each source point on fault to receiver; fdist
              call delaz(dfy2,dfx2,ry,rx,dist,az0,baz)
              fdist(ix,iy,1)= dist
              fdist(ix,iy,2)= dist
              fdist(ix,iy,3)= sqrt(dist*dist+depth*depth)
              fdist(ix,iy,4)= fdist(ix,iy,3)
              fdist(ix,iy,5)= dist
c--- for Campbell, if depth less than seismogenic depth, ignore this part 
c--- of fault by making distance 1000 km
              if(depth.lt.3.) fdist(ix,iy,4)= 1000.
              if(dip0(ift).eq.90.) then
                if(depth0(ift).le.3.) fdist(ix,iy,4)= 
     &            sqrt(dist*dist+9.)
                if(depth0(ift).gt.3.) fdist(ix,iy,4)= 
     &            sqrt(dist**2 + (depth0(ift))**2)
              endif
c---- added 6/13/02
              if(dist.gt.(tlen(ift)+dmax)) go to 850
              do 36 iatype=1,4
  36            if(fdist(ix,iy,iatype).lt.dmin(iatype)) 
     &            dmin(iatype)=fdist(ix,iy,iatype)
              if(fdist(ix,iy,5).lt.dmin(5)) then
                dmin(5)= fdist(ix,iy,5)
                testhw= 2.
              endif
 111        continue
c----- new loop for AS97 hanging wall test
            if(dip0(ift).eq.90.) go to 101
            if(nleny2.le.(nleny+1)) go to 101
            if(ix.eq.1) azold= az
            if(abs(azold-az).gt.40.) then
              do 1114 iy=nleny+1,nleny2
1114            fdist(ix,iy,5)= 10000.
              go to 101
            endif
            do 1112 iy=nleny+1,nleny2
              deltaz= (iy-1)*dlen
              depth= depth0(ift)+ deltaz*sin(abs(dip(ift)))
              deltaz= deltaz*cos(dip(ift))
c--- find surface point above downdip portion of fault
              call back(deltaz,az2,dfy,dfx,dfy2,dfx2)
c--- find distance from each source point on fault to receiver; fdist
              call delaz(dfy2,dfx2,ry,rx,dist,az0,baz)
              fdist(ix,iy,5)= dist
              if(fdist(ix,iy,5).lt.dmin(5)) then
                dmin(5)=fdist(ix,iy,5)
                dhy= (iy-(nleny+1))*dlen*cos(abs(dip(ift)))
                testhw= dhy*angtest +2.
              endif
 1112       continue
c------- end of new loop
            azold= az
 101      continue
c         write(6,*) ry,rx,testhw,dhy,angtest,dmin(5)
c  dmin is minimum distance from fault to receiver 
          if(dmin(1).gt.dmax) go to 850 
c         if(i.gt.1) write (6,*) "line 420 rec",i
c         write (6,*) "ift",ift
          if(itype(ift).eq.2) then
c----  loop through floating rupture zones along fault
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
            dmag0= dmag(ift)
c           write (6,*) "nmag0",nmag0(ift)
            do 222 m=1,nmag0(ift)+2
              xmag= magmin(ift)+(m-1)*dmag(ift)
c------find rupture length from Wells Coppersmith relation 
c        if(iftype(ift).eq.1) ruplen=  (xmag-5.16)/1.12
c        if(iftype(ift).eq.2) ruplen=  (xmag-5.00)/1.22
c        if(iftype(ift).eq.3) ruplen=  (xmag-4.86)/1.32
c        if(iftype(ift).eq.1) ruplen=  -3.55+0.74*xmag
c        if(iftype(ift).eq.2) ruplen=  -2.86+0.63*xmag
c        if(iftype(ift).eq.3) ruplen=  -3.22+0.69*xmag
c---- changed 5/29/02
c        if(iftype(ift).eq.3) ruplen=  -2.01+0.50*xmag
c---- changed 7/1/02
              ruplen= -3.22+ 0.69*xmag
              ruplen= 10**ruplen
c--- changed to dmove move for floating rupture
              nrup(m)= (tlen(ift)-ruplen)/dmove +1.
              write(6,*) "nrup",nrup(m),"dmove",dmove
              if(nrup(m).gt.2000) then
                dmove = (tlen(ift)-ruplen)/1999
                nrup(m)= (tlen(ift)-ruplen)/dmove +1.
                write(6,*) "Number of ruptures greater than 2000"
                write(6,*) "dmove changed to",dmove
                write(6,*) "nrup",nrup(m)
              endif
              if(ruplen.ge.tlen(ift)) then
                nrup(m)=1
                ruplen= tlen(ift)
              endif 
              nxx= ruplen/dlen
              nyy= width(ift)/dlen
              if(dip0(ift).eq.90.) nyy=0
              nyy2= width2/dlen 
              do 302 irup=1,nrup(m)
                testhw3(irup,m)=2.
                do 37 iatype=1,5
  37              dminr(iatype)=10000.
                xmove= dmove*(irup-1)
                ix0= xmove/dlen+1
                do 304 ix= ix0,ix0+nxx
                  if (ix.gt.nlenx) go to 304
                  do 303 iy= 1,nyy+1
                    do 303 iatype=1,5
                      dist= fdist(ix,iy,iatype)
                      if(dist.lt.dminr(iatype)) dminr(iatype)=
     &                  dist
  303             continue
c----- new hanging wall test for AS97
                  xmo2c(m,irup) = xmo2b(ix)
                  ac(m,irup) = ab(ix)
                  do 704 ilt= 1, nbranch
 704                ratenewc(m,irup,ilt) = ratenewb(ix,ilt)
                  if(dip0(ift).eq.90.) go to 304
                  if(nyy2.le.(nyy+2)) go to 304
                  do 3031 iy=nyy+2, nyy2
                    dist= fdist(ix,iy,5)
                    if(dist.lt.dminr(5)) then
                      dminr(5)= dist
                      dhy= (iy-(nyy+2))*dlen*cos(abs(dip(ift)))
                      testhw3(irup,m)= dhy*angtest+ 2.
                    endif
 3031             continue
c-----------------------
  304           continue
                do 39 iatype=1,5
 39               dmin2(m,irup,iatype)=dminr(iatype)
  302         continue
  222       continue
          endif
c         if(i.gt.50833) write (6,*) "line 490 rec",i
c
c-----------------------------------------
c----  GR with nmag0>1 with Mmax uncertainties
          if((itype(ift).eq.2).and.(nmag0(ift).gt.1).and.
     &      (itest(ift).eq.0)) then
c           do 431 mm=1,55
c 431         torate(mm)= 0.
            do 403 ilt=1,nbranch
c             write(6,*) "ilt=",ilt
              mmax= magmax(ift)+ dmbranch(ilt)
              totmo= 0.
              nmag= (mmax- magmin(ift))/dmag(ift) + 1.4
c-- loop through magnitudes
c             write(6,*) "nmag=",nmag,"nmag0(ift)",nmag0(ift)
              do 285 m=1,nmag
                xmag= magmin(ift)+(m-1)*dmag(ift)
                mm= (xmag-5.8)/0.05 +1.5
c--- loop through floating rupture zones
                do 284 irup=1,nrup(m)
                  rate= ratenewc(m,irup,ilt) - b(ift)*xmag
                  rate= 10.**rate
                  totmo= totmo+ rate*(10.**(1.5*xmag+ 9.05))
c                 torate(mm)= torate(mm)+ rate*ale(idis+mwid+1)
                  rate= rate/float(nrup(m))
                  do 283 ip=1,nper
                    do 282 ia=1,nattn(ip)
                      iftype2= iftype(ift)
                      iatype= iatten(ip,ia)
                      iatype0= iatype
                      if(iatype.gt.4) iatype=3
c--------- Campbell new dist = Camp old dist
                      if(iatype0.eq.9) iatype=4
c------ Somerville dist = BJF dist ???
                      if(iatype0.eq.7) iatype=1
c----- for BJF97
                      if(iatype0.eq.11) iatype=1
                      ii= dmin2(m,irup,iatype)/di +1.5
                      weight= wt(ip,ia,1)
                      if(dmin2(m,irup,1).gt.wtdist(ip,ia)) 
     &                  weight=wt(ip,ia,2)
                      if(ii.gt.ndist) go to 282
c------- for AS hanging wall term
                      if((iatten(ip,ia).eq.8).and.(abs(dip(ift)).
     &                  ge.0.524)) then
                        if((dmin2(m,irup,5).lt.testhw3(irup,m)).and.
     &                    (iftype(ift).eq.2)) then
                          iftype2=4
                          itest2=1
c                         if(m.eq.1) then
c                           write(16,*) ry,rx,irup,m, testhw3(irup,m)
c                         endif
                        endif
                      endif
                      if((iatten(ip,ia).eq.8).and.(abs(dip(ift)).lt.
     &                  0.524)) then
                        if((dmin2(m,irup,1).lt.2.0).and.(iftype(ift)
     &                    .eq.2)) iftype2=4
                      endif
                      if((iatten(ip,ia).eq.9).and.(iftype(ift).eq.2)
     &                  .and.(abs(dip(ift)).gt.0.7854)) iftype2=4
ccccccc  hanging wall stuff for Campbell 2002
                      if((iatten(ip,ia).eq.9).and.(xmag.ge.5.5).and.
     &                  (dmin2(m,irup,1).le.5.)) then
                        ijb= dmin2(m,irup,1)/di +1.5
                        if(iftype2.eq.2) then
c                         write(6,*) xmag,ii,ijb,phang(ii,ijb,mm,ip,2,1)
c                         read(5,*) idum
                          do 3281 k=1,nlev(ip)
 3281                       prob(icnt,k,ip)= prob(icnt,k,ip)+ 
     &                        wtbranch(ilt)*weight*rate*phang(ii,ijb
     &                        ,mm,ip,k,1)
                          go to 282
                        endif
                        if(iftype2.eq.4) then
                          do 3282 k=1,nlev(ip)
 3282                       prob(icnt,k,ip)= prob(icnt,k,ip)+ 
     &                        wtbranch(ilt)*weight*rate*phang(ii,ijb
     &                        ,mm,ip,k,2)
                          go to 282
                        endif
                      endif
ccccccccc
                      do 281 k=1,nlev(ip)
 281                    prob(icnt,k,ip)= prob(icnt,k,ip)+ 
     &                    wtbranch(ilt)*weight*rate*p(ii,mm,ip,k,
     &                    iftype2,ia)
 282                continue
 283              continue
 284            continue
 285          continue
c             write(6,*) "totmo, xmo", ift,ilt,totmo,xmo2(ift)
c             read(5,*) idum
c             if(i.eq.1) then
c             do 432 mm=1,55
c             xmag= 6.0+ (mm-1)*0.05
c 432         write(4,*) ift, xmag, torate(mm)
c             endif
 403        continue
            go to 850
          endif
c         if(i.gt.50833) write (6,*) "line 590 rec",i
ccccccccccccccccccccccccccccc
ccc---- GR with nmag0=1, one magnitude floated, with uncertainties,
c       used for some faults outside of CA
ccc--- for very long faults when Mmmax set to 7.5
          if((itype(ift).eq.2).and.(nmag0(ift).eq.1).and.
     &      (itest(ift).eq.0)) then
            totmo=0.
            do 1403 ilt=1,nbranch
              xmag= magmax(ift)+ dmbranch(ilt)
              xmo= 10.**(1.5*xmag+9.05)
              do 1404 idis= -mwid, mwid
                xmag2= xmag+ idis*0.05
                mm= (xmag2-5.8)/0.05 +1.5
c--- loop through floating rupture zones
                do 1284 irup=1,nrup(1)
                  rate= xmo2c(1,irup)/xmo
                  rate= rate/float(nrup(1))
c                 if(i.gt.50833) write(6,*) "ale2 ind",idis+mwid+1
                  totmo= totmo+ ale2(idis+mwid+1)*rate*10.**(
     &              1.5*xmag2+ 9.05)
                  do 1283 ip=1,nper
                    do 1282 ia=1,nattn(ip)
                      iftype2= iftype(ift)
                      iatype= iatten(ip,ia)
                      iatype0= iatype
                      if(iatype.gt.4) iatype=3
c--------- Campbell new dist = Camp old dist
                      if(iatype0.eq.9) iatype=4
c------ Somerville dist = BJF dist ???
                      if(iatype0.eq.7) iatype=1
c----- for BJF97
                      if(iatype0.eq.11) iatype=1
                      ii= dmin2(1,irup,iatype)/di +1.5
                      weight= wt(ip,ia,1)
                      if(dmin2(1,irup,1).gt.wtdist(ip,ia)) 
     &                  weight=wt(ip,ia,2)
                      if(ii.gt.ndist) go to 1282
c------- for AS hanging wall term
                      if((iatten(ip,ia).eq.8).and.
     &                  (abs(dip(ift)).ge.0.524)) then
                        if((dmin2(1,irup,5).lt.testhw3(irup,1)).and.
     &                    (iftype(ift).eq.2)) iftype2=4
                      endif
                      if((iatten(ip,ia).eq.8).and.
     &                  (abs(dip(ift)).lt.0.524)) then
                        if((dmin2(1,irup,1).lt.2.0).and.
     &                    (iftype(ift).eq.2)) iftype2=4
                      endif
                      if((iatten(ip,ia).eq.9).and.(iftype(ift).eq.2)
     &                  .and.(abs(dip(ift)).gt.0.7854)) iftype2=4
ccccccc  hanging wall stuff for Campbell 2002
                      if((iatten(ip,ia).eq.9).and.(xmag.ge.5.5).and.
     &                  (dmin2(1,irup,1).le.5.)) then
                        ijb= dmin2(1,irup,1)/di +1.5
                        if(iftype2.eq.2) then
                          do 4281 k=1,nlev(ip)
 4281                       prob(icnt,k,ip)= prob(icnt,k,ip)+
     &                        ale2(idis+mwid+1)*wtbranch(ilt)*
     &                        weight*rate*phang(ii,ijb,mm,ip,k,1)
                          go to 1282
                        endif
                        if(iftype2.eq.4) then
                          do 4282 k=1,nlev(ip)
 4282                       prob(icnt,k,ip)= prob(icnt,k,ip)+
     &                        ale2(idis+mwid+1)*wtbranch(ilt)*
     &                        weight*rate*phang(ii,ijb,mm,ip,k,2)
                          go to 1282
                        endif
                      endif
ccccccccc
                      do 1281 k=1,nlev(ip)
1281                    prob(icnt,k,ip)= prob(icnt,k,ip)+
     &                    ale2(idis+mwid+1)*wtbranch(ilt)*
     &                    weight*rate*p(ii,mm,ip,k,iftype2,ia)
1282                continue
1283              continue
1284            continue
1404          continue
1403        continue
c           write(6,*) totmo,xmo2(ift)
c           read(5,*) idum
            go to 850
          endif
c---------------------
c------ all GR without uncertainties
          if((itype(ift).eq.2).and.(itest(ift).eq.1)) then
c-- loop through magnitudes
            iftype2= iftype(ift)
            do 2285 m=1,nmag0(ift)
              xmag= magmin(ift)+(m-1)*dmag(ift)
              mm= (xmag-5.8)/0.05 +1.5
c--- loop through floating rupture zones
              do 2284 irup=1,nrup(m)
                rate= ac(m,irup) - b(ift)*xmag
                rate= 10.**rate
                rate= rate/float(nrup(m))
                do 2283 ip=1,nper
                  do 2282 ia=1,nattn(ip)
                    iftype2= iftype(ift)
                    iatype= iatten(ip,ia)
                    iatype0= iatype
                    if(iatype.gt.4) iatype=3
ccc  Camp new dist = Camp old dist
                    if(iatype0.eq.9) iatype=4
cccc Somerville dist = BJF dist ???
                    if(iatype0.eq.7) iatype=1
c----- for BJF97
                    if(iatype0.eq.11) iatype=1
                    ii= dmin2(m,irup,iatype)/di +1.5
                    weight= wt(ip,ia,1)
                    if(dmin2(m,irup,1).gt.wtdist(ip,ia))
     &                weight=wt(ip,ia,2)
                    if(ii.gt.ndist) go to 2282
c---- for AS hanging wall term
                    if((iatten(ip,ia).eq.8).and.(abs(dip(ift)).ge.
     &                0.524)) then
                      if((dmin2(m,irup,5).lt.testhw3(irup,m)).and.
     &                  (iftype(ift).eq.2)) iftype2=4 
                    endif
                    if((iatten(ip,ia).eq.8).and.(abs(dip(ift)).lt.
     &                0.524)) then
                      if((dmin2(m,irup,1).lt.2.0).and.(iftype(ift).
     &                  eq.2)) iftype2=4
                    endif
                    if((iatten(ip,ia).eq.9).and.(iftype(ift).eq.2).and.
     &                (abs(dip(ift)).gt.0.7854)) iftype2=4
ccccccc  hanging wall stuff for Campbell 2002
                    if((iatten(ip,ia).eq.9).and.(xmag.ge.5.5).and.
     &                (dmin2(m,irup,1).le.5.)) then
                      ijb= dmin2(m,irup,1)/di +1.5
                      if(iftype2.eq.2) then
                        do 5281 k=1,nlev(ip)
 5281                     prob(icnt,k,ip)= prob(icnt,k,ip)+ 
     &                      weight*rate*phang(ii,ijb,mm,ip,k,1)
                        go to 2282
                      endif
                      if(iftype2.eq.4) then
                        do 5282 k=1,nlev(ip)
 5282                     prob(icnt,k,ip)= prob(icnt,k,ip)+ 
     &                      weight*rate*phang(ii,ijb,mm,ip,k,2)
                        go to 2282
                      endif
                    endif
ccccccccc
                    do 2281 k=1,nlev(ip)
2281                  prob(icnt,k,ip)= prob(icnt,k,ip)+
     &                  weight*rate*p(ii,mm,ip,k,iftype2,ia)
2282              continue
2283            continue
2284          continue
2285        continue
            go to 850
          endif
c---------------------------
c--- for characteristic event, tested
c---- characteristic with uncertainties
          if((itype(ift).eq.1).and.(itest(ift).eq.0)) then
            xmag= cmag(ift)
            xmo= 1.5*xmag +9.05
            xmo= 10.**xmo
            rate= crate(ift)
            xmorate= xmo *rate
            do 203 ilt=1,nbranch
              xmag= cmag(ift) + dmbranch(ilt)
              xmo= 1.5*xmag+9.05
              xmo= 10.**xmo
              rate= xmorate/xmo
              totmo =0.
              do 244 idis= -mwid, mwid
                xmag2= xmag + idis*0.05
                totmo= totmo+ rate*ale2(idis+mwid+1)*(10.**(1.5*xmag2
     &            +9.05))
                mm= (xmag2-5.8)/0.05 +1.5
                do 214 ip=1,nper
                  do 204 ia=1,nattn(ip)
                    iftype2= iftype(ift)
                    iatype= iatten(ip,ia)
                    iatype0=iatype
                    if(iatype.gt.4) iatype=3
ccc  Camp new dist = Camp old dist
                    if(iatype0.eq.9) iatype=4
cccc Somerville dist = BJF dist ???
                    if(iatype0.eq.7) iatype=1
c----- for BJF97
                    if(iatype0.eq.11) iatype=1
                    ii= dmin(iatype)/di +1.5 
                    weight= wt(ip,ia,1)
                    if(dmin(1).gt.wtdist(ip,ia)) weight=wt(ip,ia,2)
c------- for AS hanging wall term
                    if((iatten(ip,ia).eq.8).and.(abs(dip(ift)).ge.
     &                0.524)) then
                      if((dmin(5).le.testhw).and.(iftype(ift).eq.
     &                  2)) then
                        iftype2=4
c                       write(16,*) ry,rx, itest2
                      endif
                    endif
                    if((iatten(ip,ia).eq.8).and.(abs(dip(ift)).lt.
     &                0.524)) then
                      if((dmin(1).le.2.).and.(iftype(ift).eq.2)) then
                        iftype2=4
                      endif
                    endif
c                   write(14,*) p(ii,mm,ip,15,4,ia),p(ii,mm,ip,15,2,ia)
c                   write(14,*) dmin(1),dmin(iatype),xmag,ii,mm,ia
c                   read(5,*) idum
                    if((iatten(ip,ia).eq.9).and.(iftype(ift).eq.2).and.
     &                (abs(dip(ift)).gt.0.7854)) iftype2=4
ccccccc  hanging wall stuff for Campbell 2002
                    if((iatten(ip,ia).eq.9).and.(xmag.ge.5.5).and.
     &                (dmin(1).le.5.)) then
                      ijb= dmin(1)/di +1.5
                      if(iftype2.eq.2) then
                        do 6281 k=1,nlev(ip)
 6281                     prob(icnt,k,ip)= prob(icnt,k,ip)+
     &                      ale2(idis+mwid+1)*wtbranch(ilt)*weight*
     &                      rate*phang(ii,ijb,mm,ip,k,1)
                        go to 204
                      endif
                      if(iftype2.eq.4) then
                        do 6282 k=1,nlev(ip)
 6282                     prob(icnt,k,ip)= prob(icnt,k,ip)+
     &                      ale2(idis+mwid+1)*wtbranch(ilt)*weight*
     &                      rate*phang(ii,ijb,mm,ip,k,2)
                        go to 204
                      endif
                    endif
ccccccccc
                    do 202 k=1,nlev(ip)
 202                  prob(icnt,k,ip)= prob(icnt,k,ip)+
     &                  ale2(idis+mwid+1)*wtbranch(ilt)*weight*
     &                  rate*p(ii,mm,ip,k,iftype2,ia)
 204              continue
 214            continue
 244          continue
c             write(6,*) "totmo, xmorate",ift,ilt,totmo, xmorate
c             read(5,*) idum
 203        continue
          endif
ccccccccccccccccc
cc-- characteristic without uncertainties
          if((itype(ift).eq.1).and.(itest(ift).eq.1)) then
            xmag= cmag(ift)
            rate= crate(ift)
            iftype2= iftype(ift)
            mm= (xmag-5.8)/0.05 +1.5
            do 3203 ip=1,nper
              do 3203 ia=1,nattn(ip)
                iftype2= iftype(ift)
                iatype= iatten(ip,ia)
                iatype0=iatype
                if(iatype.gt.4) iatype=3
ccc    Camp new dist = Camp old dist
                if(iatype0.eq.9) iatype=4
cccc Somerville dist = BJF dist ???
                if(iatype0.eq.7) iatype=1
c----- for BJF97
                if(iatype0.eq.11) iatype=1
                ii= dmin(iatype)/di +1.5 
                weight= wt(ip,ia,1)
                if(dmin(1).gt.wtdist(ip,ia)) weight=wt(ip,ia,2)
c---- AS hanging wall term
                if((iatten(ip,ia).eq.8).and.(abs(dip(ift)).ge.0.524))
     &            then
                  if((dmin(5).le.testhw).and.(iftype(ift).eq.2)) then
                    iftype2=4
c                   write(16,*) ry,rx, itest2
                  endif
                endif
                if((iatten(ip,ia).eq.8).and.(abs(dip(ift)).lt.0.524))
     &            then
                  if((dmin(1).le.2.).and.(iftype(ift).eq.2)) then
                    iftype2=4
                  endif
                endif
                if((iatten(ip,ia).eq.9).and.(iftype(ift).eq.2).and.
     &            (abs(dip(ift)).gt.0.7854)) iftype2=4
ccccccc  hanging wall stuff for Campbell 2002
                if((iatten(ip,ia).eq.9).and.(xmag.ge.5.5).and.
     &            (dmin(1).le.5.)) then
                  ijb= dmin(1)/di +1.5
                  if(iftype2.eq.2) then
                    do 7281 k=1,nlev(ip)
 7281                 prob(icnt,k,ip)= prob(icnt,k,ip)+ 
     &                  weight*rate*phang(ii,ijb,mm,ip,k,1)
                    go to 3203
                  endif
                  if(iftype2.eq.4) then
                    do 7282 k=1,nlev(ip)
 7282                 prob(icnt,k,ip)= prob(icnt,k,ip)+ 
     &                  weight*rate*phang(ii,ijb,mm,ip,k,2)
                    go to 3203
                  endif
                endif
ccccccccc
                do 3202 k=1,nlev(ip)
 3202             prob(icnt,k,ip)= prob(icnt,k,ip)+
     &              weight*rate*p(ii,mm,ip,k,iftype2,ia)
 3203       continue
            go to 850
          endif
ccccccccccccccccccc
 850    continue
c--- write output files every 1000 sites
        if((icnt.eq.1000).or.(i.eq.nrec)) then
c         write(6,*) i,prob(icnt,1,1)
          iend= 1000
          if(i.eq.nrec) iend=icnt
          do 115 ip=1,nper
            ndata= iend*nlev(ip)
            do 50 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)
   50       continue
            call putbufx(ifp(ip),out,ndata,readn)
c           readn= readn/4
c           write(6,*) ndata,readn
 115      continue
          icnt= 0
          do 230 ii=1,1000
            do 230 ip=1,nper
              do 230 k=1,nlev(ip)
 230            prob(ii,k,ip)= 0.
        endif
        icnt= icnt+1
  100 continue
      end
c
      subroutine delaz(sorlat,sorlon,stnlat,stnlon,delta,az,baz)
      coef= 3.14159/180.
      xlat= sorlat*coef
      xlon= sorlon*coef
      st0= cos(xlat)
      ct0= sin(xlat) 
      phi0= xlon
      xlat= stnlat*coef
      xlon= stnlon*coef
      ct1= sin(xlat)
      st1= cos(xlat)
      sdlon= sin(xlon-phi0)
      cdlon= cos(xlon-phi0)
      cdelt= st0*st1*cdlon+ct0*ct1
      x= st0*ct1-st1*ct0*cdlon
      y= st1*sdlon
      sdelt= sqrt(x*x+y*y)
      delta= atan2(sdelt,cdelt)
      delta= delta/coef
      az= atan2(y,x)
      az= az/coef
      x= st1*ct0-st0*ct1*cdlon
      y= -sdlon*st0
      baz= atan2(y,x)
      baz= baz/coef
      delta= delta*111.2
      return
      end         
      subroutine back(delta0,az,sorlat,sorlon,stnlat,stnlon)
      coef= 3.14159/180.
      delta= delta0/111.11
      delta= delta*coef
      sdelt= sin(delta)
      cdelt= cos(delta)
      xlat= sorlat*coef
      xlon= sorlon*coef
      az2= az*coef
      st0= cos(xlat)
      ct0= sin(xlat)
      phi0= xlon
      cz0= cos(az2)
      ct1= st0*sdelt*cz0+ct0*cdelt
      x= st0*cdelt-ct0*sdelt*cz0
      y= sdelt*sin(az2)
      st1= sqrt(x*x+y*y)
      dlon= atan2(y,x)
      stnlat= atan2(ct1,st1)
      stnlat= stnlat/coef
      stnlon= phi0+dlon
      stnlon= stnlon/coef
      return
      end

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(22,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')
 900  format(i2)

      do 20 j = 1, 5

        read(3,10)header
        read(3,900) idum

10      format(a)


20    continue


      nf = 11

      nmag = 14

      nd = 18

      jdl = 0

      do 30 jf = 1, nf-2

        f(jf) = 10.**(-0.3 + float(jf-1)*0.2)

30    continue

      f(nf) = 89.

      f(nf-1) = 99.

c

      do 50 jm = 1, nmag

        read(3,900) idum
        read(3,900) idum
        read(3,35)xmag(jm)

35      format(13x,f5.2)
        read(3,900) idum

        read(3,45)
        read(3,900) idum
        read(3,45)
        read(3,900) idum

        do 40 jd = 1, nd

          read(3,45)rlog(jd),(gma(jf,jd,jm), jf = 1,nf)

          read(3,900) idum
45        format(f6.2,2x,11f6.2)

40      continue

50    continue
      close(3)
c      write(6,*) gma(1,1,1)
      return
      end

c***************************************************************
c*****************************************************************

      subroutine table(gma,amag,R,freq,amean,f,xmag,rlog)

      dimension xmag(20),gma(22,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(6,*) ' ERROR. FREQ. OUTSIDE RANGE.'

      if (freq .gt. 20. .and. freq .lt. 89.)

     *   write(6,*)' ERROR. FREQ. OUTSIDE RANGE'

      jfu = jfl + 1

      fracf = 0.

      if (freq .ge. 20.) go to 15

      fracf = (alog10 (freq) - alog10 (f(jfl)) ) /

     *        (alog10 (f(jfu)) - alog10 (f(jfl)) )

15    continue

c     Find bounding jm values.

      jml = 0

      do 20 j = 1, nmag-1

         if (amag .ge. xmag(j) .and. amag .lt. xmag(j+1)) jml = j

20    continue

      if (amag .eq. xmag(nmag)) jml = nmag

      if (jml .eq. 0) write(6,*)' ERROR. MAGNITUDE OUTSIDE RANGE.'

      jmu = jml + 1

      fracm = (amag - xmag(jml)) / (xmag(jmu)-xmag(jml))

c     find bounding distance values.
      jdl= 0

      do 30 j = 1, nd-1

        if (rl .ge. rlog(j) .and. rl .lt. rlog(j+1)) jdl= j

30    continue

      if (rl .eq. rlog(nd)) jdl = nd

      if (jdl .eq. 0) write(6,*)' ERROR. DISTANCE OUTSIDE RANGE.'

      jdu = jdl + 1

      fracd = (rl - rlog(jdl)) / (rlog(jdu)-rlog(jdl))

      do 40 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)

40    continue

      amean = avg(jfl) + fracf * (avg(jfu) - avg(jfl))

      return

      end
cccccccccccccccccccccccccccccc
      subroutine getBJF93(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
      common / atten / p, phang, xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
        write(6,*) "enter b1,b2,b3,b4,b5,h,sigma,thrust,norm"
        read(1,*) b1,b2,b3,b4,b5,h,sigma,thrust,xnorm
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=45 dmag=0.05
       pi= 4.*atan(1.)
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
        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(xmag.gt.7.7) xmag=7.7
          sig= sigma/0.4343
          if(dist0.lt.distnf) sig= sig+ sigmanf
          sigmasq= sig*sqrt(2.)
          dist= sqrt(dist0*dist0+h*h)
          gnd= b1+b2*(xmag-6.)+b3*(xmag-6.)**2
          gnd= gnd+ b4*dist + b5*alog10(dist)
c--- convert from psv in cm/sec to psa in g
          if(period.ne.0.) gnd=gnd+alog10(2.*pi/(980.*period))
          gnd2= gnd+ thrust
          gnd3= gnd+ xnorm
          gnd= gnd/0.4343
          gnd2=gnd2/0.4343
          gnd3= gnd3/0.4343
      do 102 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
  102 continue
  104 continue
  103 continue
      return
      end
cccccccccccccccc
      subroutine getToro(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
      common / atten / p, phang, xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
        write(6,*) "enter b1,b2,b3,b4,b5,b6,h,sigma"
        read(1,*) tc1,tc2,tc3,tc4,tc5,tc6,
     &    th,tsigma,clamp
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=45 dmag=0.05
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
        xmag= xmag0
        if(iconv(ip,ia).eq.1) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
          sigma= tsigma/0.4343
          if(dist0.lt.distnf) sigma= sigma+ sigmanf
          sigmasq= sigma*sqrt(2.)
          dist= sqrt(dist0*dist0+th*th)
          gnd= tc1+tc2*(xmag-6.)+ tc3*((xmag-6.)**2)
          gnd= gnd-tc4*alog(dist)-tc6*dist
          factor= alog(dist/100.)
          if(factor.gt.0.) gnd= gnd-(tc5-tc4)*factor
c---following is for clipping gnd motions: 1.5g PGA, 3.0g for rest 
          if((period.eq.0.).and.(gnd.gt.0.405)) gnd=0.405
          if((period.eq.0.3).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.2).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.1).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.5).and.(gnd.gt.1.099)) gnd=1.099
      test= exp(gnd + 3.*sigmasq/sqrt(2.))
      if((clamp.lt.test).and.(clamp.gt.0.)) clamp2= alog(clamp)
      if(clamp.ge.test) clamp2= gnd+ 3.*sigmasq/sqrt(2.)
      if(clamp.eq.0.) clamp2= gnd+ 3.*sigmasq/sqrt(2.)
      tempgt3= (gnd- clamp2)/sigmasq
      probgt3= (eerf(tempgt3)+1.)/2.
      do 105 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
  105 continue
  104 continue
  103 continue
      return
      end
cccccccccccccccccc
      subroutine getSadigh(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
      common / atten / p, phang, xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
         write(6,*) "enter c1,c2,c3,c4,c5,c6,h,thrust"
         read(1,*) sc1,sc2,sc3,sc4,sc5,
     &     sc6,sh,sthrust
         read(1,*) sd1,sd2,sd3,sd4,sd5,sd6,sd7
         read(1,*) ssigma1,ssigmacoef,ssigma2,smagsig
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=45 dmag=0.05
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
        xmag= xmag0
        if(iconv(ip,ia).eq.1) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
          if(xmag.lt.smagsig) sig= ssigma1- ssigmacoef*xmag
          if(xmag.ge.smagsig) sig= ssigma2
          if(dist0.lt.distnf) sig= sig+ sigmanf
          sigmasq= sig*sqrt(2.)
          if(xmag.lt.6.5) then
          gnd= sc1+sc2*xmag+sc3*((8.5-xmag)**2.5)
          gnd= gnd+ sc4*alog(dist0+ exp(sc5+sc6*xmag))
            endif
          if(xmag.ge.6.5) then
          gnd= sd1+sd2*xmag+sd3*((8.5-xmag)**2.5)
          gnd= gnd+ sd4*alog(dist0+exp(sd5+sd6*xmag))
            endif
          gnd= gnd+sd7*alog(dist0+2.)
          gnd2= gnd+sthrust
          gnd3= gnd
      do 102 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
  102 continue
  104 continue
  103 continue
      return
      end
cccccccccccccccccccccccc
      subroutine getCamp(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
      common / atten / p, phang, xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
         write(6,*) "enter c1,c2,c3,c4,c5,h,thrust"
         read(1,*) cc1,cc2,cc3,cc4,cc5,ch
         read(1,*) cd1,cd2,cd3
         read(1,*) cs1,cs2
         read(1,*) csigma1,csigmacoef,csigma2,cmagsig
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=45 dmag=0.05
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
        xmag= xmag0
        if(iconv(ip,ia).eq.1) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
          if(xmag.lt.cmagsig) csigma= csigma1- 
     &     csigmacoef*xmag
          if(xmag.ge.cmagsig) csigma= csigma2
          if(dist0.lt.distnf) csigma= csigma+ sigmanf
          sigmasq= csigma*sqrt(2.)
          arg=  dist0**2+ (cc4*exp(cc5*xmag))**2
          arg= sqrt(arg)
          gnd= cc1+cc2*xmag+cc3*alog(arg)
          if(dist0.eq.0.) then
             gnd=0
             gnd2=0
             gnd3=0
             go to 707
             endif
          gnd= gnd+ cs1+ cs2*alog(dist0)
          gnd3= gnd
c--- following for Campbell thrust fault
          gnd2= gnd+ cd1+cd2*alog(dist0)+cd3*xmag
  707     continue
      do 102 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
  102 continue
  104 continue
  103 continue
      return
      end
cccccccccccccccccccccc
      subroutine getAB95(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
      common / atten / p, phang, xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
      dimension gma(22,30,20),amag(20),rlog(20),f(20)
         write(6,*) "enter bcfac, coeficients, depth, sigma"
         read(1,*) bcfac, ab1,ab2,ab3,ab4,
     &            abdepth,absigma,clamp
         call gettab(gma,f,amag,rlog)
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=45 dmag=0.05
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
      dist= dist0
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
        xmag= xmag0
c--- following for Atkinson and Boore (95) look-up table
        if(iconv(ip,ia).eq.1) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
          sigma= absigma/0.4343
          if(dist0.lt.distnf) sigma= sigma+ sigmanf
          sigmasq= sigma*sqrt(2.)
cccccccccc  if fault depth0 term ne 0 then this depth term not needed
c          dist= sqrt(dist0*dist0+ abdepth*abdepth)
          if(dist.lt.10.) dist=10.
          if(period.gt.0.) freq=1./period
          if(period.eq.0.) freq= 99.
c---------if M> 7.25 use quadratic formula
          if((xmag.gt.7.25).or.(dist.gt.501.)) then
             gnd= ab1+ab2*(xmag-6.)+ab3*((xmag-6)**2)
     &            -alog10(dist)- ab4*dist
             go to 402
             endif
          call table(gma,xmag,dist,freq,gnd,f,amag,rlog)
cc- convert to fraction of g
  402     fac= -alog10(980.)
          gnd= gnd+fac + alog10(bcfac)
          gnd= gnd/0.4343
c----------- clipping
          if((period.eq.0.).and.(gnd.gt.0.405)) gnd=0.405
          if((period.eq.0.3).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.2).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.1).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.5).and.(gnd.gt.1.099)) gnd=1.099
          gnd2= gnd
          gnd3= gnd
      test= exp(gnd + 3.*sigmasq/sqrt(2.))
      if((clamp.lt.test).and.(clamp.gt.0.)) clamp2= alog(clamp)
      if(clamp.ge.test) clamp2= gnd+ 3.*sigmasq/sqrt(2.)
      if(clamp.eq.0.) clamp2= gnd+ 3.*sigmasq/sqrt(2.)
      tempgt3= (gnd- clamp2)/sigmasq
      probgt3= (eerf(tempgt3)+1.)/2.
      do 105 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
  105 continue
  104 continue
  103 continue
      return
      end
ccccccccccccccccccccccc
      subroutine getFEA(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
      common / atten / p, phang, xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
      dimension tabdist(21), gma(22,30,20)
      character nametab*30
         write(6,*) "enter file name of table"
         read(1,900) nametab
 900  format(a256)
         write(6,*) "enter depth, sigma, log factor"
         read(1,*) bdepth,bsigma,xlogfac,clamp
         open(unit=15,file=nametab,status='old')
         read(15,900) adum
         do 80 idist=1,21
   80    read(15,*) tabdist(idist),(gma(imag,idist,ip),imag=1,20)
         close(15)
c---- following for new Boore look-up table
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=45 dmag=0.05
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
        xmag= xmag0
        if(iconv(ip,ia).eq.1) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
          sigma= bsigma/0.4343
          if(dist0.lt.distnf) sigma= sigma+ sigmanf
          sigmasq= sigma*sqrt(2.)
          dist= dist0
ccccccccc if fault depth0 ne 0. then this depth term not needed 
c          dist= sqrt(dist0*dist0+bdepth*bdepth)
          if(dist.lt.10.) dist=10.
          imag= (xmag-4.4)/0.2 +1
          rdist= alog10(dist)
          idist= (rdist-tabdist(1))/0.1 + 1
          xm1= (imag-1)*0.2 + 4.4
          fracm= (xmag-xm1)/0.2
          xd1= (idist-1)*0.1 +tabdist(1)
          fracd= (rdist-xd1)/0.1
          idist1= idist+1
          if(idist.gt.21) idist=21
          if(idist1.gt.21) idist1=21
          gm1= gma(imag,idist,ip)
          gm2= gma(imag+1,idist,ip)
          gm3= gma(imag,idist1,ip)
          gm4= gma(imag+1,idist1,ip)
          arl= gm1 + fracm*(gm2-gm1)
          aru= gm3 + fracm*(gm4-gm3)
          gnd= arl +fracd*(aru-arl)
          gnd= gnd+ xlogfac
          gnd= gnd/0.4343
c          if(dist0.gt.950.) then
c            taper= (1001.-dist)/50.
c            if(taper.le.0.01) taper=.01
c            taper= log(taper)
c            gnd=gnd+taper
c            endif
c--- following is for clipping 
          if((period.eq.0.).and.(gnd.gt.0.405)) gnd=0.405
          if((period.eq.0.3).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.2).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.1).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.5).and.(gnd.gt.1.099)) gnd=1.099
      test= exp(gnd + 3.*sigmasq/sqrt(2.))
      if((clamp.lt.test).and.(clamp.gt.0.)) clamp2= alog(clamp)
      if(clamp.ge.test) clamp2= gnd+ 3.*sigmasq/sqrt(2.)
      if(clamp.eq.0.) clamp2= gnd+ 3.*sigmasq/sqrt(2.)
      tempgt3= (gnd- clamp2)/sigmasq
      probgt3= (eerf(tempgt3)+1.)/2.
      do 105 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
  105 continue
  104 continue
  103 continue
      return
      end
ccccccccccccccccccccccccc
      subroutine getSomer(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
c---- Somerville et al. (2001) for CEUS
      common / atten / p, phang, xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
        read(1,*) a1, a2, a3, a4, a5, a6, a7, sig0, clamp
      dist1= sqrt(50.*50.+ 6.*6.)
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=45 dmag=0.05
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
      sig= sig0
      if(dist0.lt.distnf) sig= sig0+ sigmanf
      sigmasq= sig*sqrt(2.)
      dist= sqrt(dist0*dist0+ 6.*6.)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
      xmag= xmag0
      if(iconv(ip,ia).eq.1) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
      if(dist0.lt.50.) then
      gnd= a1 + a2*(xmag-6.4) + a3*alog(dist)+a4*(xmag-6.4)*alog(dist)
     &  + a5*dist0 + a7*(8.5-xmag)*(8.5-xmag)
        endif
      if(dist0.ge.50.) then
      gnd= a1 + a2*(xmag-6.4) + a3*alog(dist1)+a4*(xmag-6.4)*alog(dist)
     & +a5*dist0 + a6*(alog(dist)-alog(dist1))+a7*(8.5-xmag)*(8.5-xmag)
      endif
c      write(16,*) period,dist0,xmag,exp(gnd)
c--- following is for clipping 
          if((period.eq.0.).and.(gnd.gt.0.405)) gnd=0.405
          if((period.eq.0.3).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.2).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.1).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.5).and.(gnd.gt.1.099)) gnd=1.099
      gnd2= gnd
      gnd3= gnd
      test= exp(gnd + 3.*sigmasq/sqrt(2.))
      if((clamp.lt.test).and.(clamp.gt.0.)) clamp2= alog(clamp)
      if(clamp.ge.test) clamp2= gnd+ 3.*sigmasq/sqrt(2.)
      if(clamp.eq.0.) clamp2= gnd+ 3.*sigmasq/sqrt(2.)
      tempgt3= (gnd- clamp2)/sigmasq
      probgt3= (eerf(tempgt3)+1.)/2.
      do 105 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
  105 continue
  104 continue
  103 continue
      return
      end
ccccccccccccccccccccccccc
      subroutine getAS(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
      common / atten / p, phang, xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
         write(6,*) "enter a1,a2,a3,a4,a5,a6"
         read(1,*) as1,as2,as3,as4,as5,as6
         write(6,*) "enter a9,a12,a13,c1,c4,b5,b6"
         read(1,*) as9,as12,as13,asc1,asc4,b5,b6
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=45 dmag=0.05
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
        xmag= xmag0
        if(iconv(ip,ia).eq.1) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
        if(xmag.lt.7.0) sig= b5- b6*(xmag-5.0)
        if(xmag.ge.7.0) sig= b5- 2.*b6 
        if(dist0.lt.distnf) sig= sig+ sigmanf
        sigmasq= sig*sqrt(2.)
          r= sqrt(dist0*dist0+ asc4*asc4)
          if(xmag.le.asc1) then
          gnd= as1+as2*(xmag-asc1)+as12*((8.5-xmag)**2)
     &         + (as3+ as13*(xmag-asc1))*alog(r)
          endif
          if(xmag.gt.asc1) then
          gnd= as1+as4*(xmag-asc1)+as12*((8.5-xmag)**2)
     &         + (as3+ as13*(xmag-asc1))*alog(r)
          endif
          if(xmag.le.5.8) fltfac= as5
          if((xmag.gt.5.8).and.(xmag.lt.asc1)) fltfac= as5
     &        +(as6-as5)*(xmag-5.8)/(asc1-5.8)
          if(xmag.ge.asc1) fltfac= as6
          gnd2= gnd+ fltfac
          gnd3= gnd
c------- for hanging wall use iftype2=4, gnd4, temp4
             if(xmag.lt.6.5) fac1= xmag-5.5
             if(xmag.ge.6.5) fac1= 1.
             if(dist0.lt.4.) fac2=0.
             if((dist0.ge.4.).and.(dist0.lt.8.)) 
     &          fac2= as9*(dist0-4.)/4.
             if((dist0.ge.8.).and.(dist0.lt.18.)) fac2= as9
             if((dist0.ge.18.).and.(dist0.lt.25.)) 
     &          fac2= as9*(1.-(dist0-18.)/7.)
             if(dist0.ge.25.) fac2= 0.
             gnd4= gnd2+ fac1*fac2
c      write(14,*) period,dist0,xmag,exp(gnd)
      do 102 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      temp4= (gnd4- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
      p(ii,m,ip,k,4,ia)= (eerf(temp4)+1.)/2.
      p(ii,m,ip,k,4,ia)= (p(ii,m,ip,k,4,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,4,ia).lt.0.) p(ii,m,ip,k,4,ia)=0.
  102 continue
  104 continue
  103 continue
      return
      end
ccccccccccccccccccccccccccccc
      subroutine getCamp2000(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
c  Campbell and Bozorgnia (2000)
      common / atten / p, phang, xlev, nlev, iconv
      dimension phang(30,6,55,7,20,2)
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      read(1,*) c1,c2,c3,c4,c5,c6,c7,c8,c9
      read(1,*) c10, c11, csite, c15
ccc   sigma??
      read(1,*) csigma1,csigmacoef,csigma2,cmagsig
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=45 dmag=0.05
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
      dist= dist0
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
        xmag= xmag0
        if(iconv(ip,ia).eq.1) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
      gnd= c1 + c2*xmag + c3*(8.5-xmag)*(8.5-xmag)
      arg0= c8*xmag + c9*(8.5-xmag)*(8.5-xmag)
      arg= dist*dist + ((c5 + 0.5*c6 +0.5*c7)*exp(arg0))**2
      arg= sqrt(arg)
      gnd = gnd + c4*alog(arg) + csite
c----  gnd2 for thrust faults
      gnd2= gnd+ c11
c----- gnd3 for normal faults
      gnd3= gnd
c--- gnd4 for reverse faults
      gnd4= gnd+ c10
          if(xmag.lt.cmagsig) csigma= csigma1- 
     &     csigmacoef*xmag
          if(xmag.ge.cmagsig) csigma= csigma2
          if(dist0.lt.distnf) csigma= csigma+ sigmanf
      sigmasq= csigma*sqrt(2.)
      do 102 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      temp4= (gnd4- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
      p(ii,m,ip,k,4,ia)= (eerf(temp4)+1.)/2.
      p(ii,m,ip,k,4,ia)= (p(ii,m,ip,k,4,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,4,ia).lt.0.) p(ii,m,ip,k,4,ia)=0.
  102 continue
c      if(dist.le.30.) then
      if(ii.le.30) then
      do 202 ijb=1,6
      rjb= (ijb-1)
      fhw= (5.-rjb)/5.
      if(xmag0.le.6.5) fhw= fhw*(xmag0-5.5)
      if(dist.lt.8.) fhw= fhw*c15*dist/8.
      if(dist.ge.8.) fhw= fhw*c15
c---- gnd5 for thrust faults hanging wall
      gnd5= gnd2 + fhw
c---- gnd6 for reverse faults hanging wall
      gnd6= gnd4 + fhw
c      if(dist.ge.3.) write(12,*) ip,xmag0,dist,rjb,exp(gnd),exp(gnd2),
c     &  exp(gnd3),exp(gnd4),exp(gnd5),exp(gnd6)
      do 203 k=1,nlev(ip)
      temp2h= (gnd5- xlev(k,ip))/sigmasq
      temp4h= (gnd6- xlev(k,ip))/sigmasq
      phang(ii,ijb,m,ip,k,1)= (eerf(temp2h)+1.)/2.
      phang(ii,ijb,m,ip,k,1)= (phang(ii,ijb,m,ip,k,1)-(1.3499e-3))
     &   /0.99865
      if(phang(ii,ijb,m,ip,k,1).lt.0.) phang(ii,ijb,m,ip,k,1)=0.
      phang(ii,ijb,m,ip,k,2)= (eerf(temp4h)+1.)/2.
      phang(ii,ijb,m,ip,k,2)= (phang(ii,ijb,m,ip,k,2)-(1.3499e-3))
     &  /0.99865
      if(phang(ii,ijb,m,ip,k,2).lt.0.) phang(ii,ijb,m,ip,k,2)=0.
  203 continue
  202 continue
      endif
  104 continue
  103 continue
      return
      end
cccccccccccccc
      subroutine getCampCEUS(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
c  Campbell CEUS 2001
      common / atten / p, phang,  xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
      read(1,*) c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,clamp
ccc   sigma??
      read(1,*) csigma1,csigmacoef,csigma2,cmagsig
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=45 dmag=0.05
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
      dist= dist0
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
        xmag= xmag0
        if(iconv(ip,ia).eq.1) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
      gnd= c1 + c2*xmag + c3*(8.5-xmag)*(8.5-xmag)
      arg= dist*dist + (c5*exp(c6*xmag))**2
      arg= sqrt(arg)
      fac=0.
      if(dist.gt.70.) fac= c7*(alog(dist)- alog(70.))
      if(dist.gt.130.) fac= fac+ c8*(alog(dist)-alog(130.))
      gnd = gnd + c4*alog(arg) + fac +(c9+c10*xmag)*dist
c      write(12,*) period,dist0,xmag,exp(gnd)
c--- following is for clipping 
          if((period.eq.0.).and.(gnd.gt.0.405)) gnd=0.405
          if((period.eq.0.3).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.2).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.1).and.(gnd.gt.1.099)) gnd=1.099
          if((period.eq.0.5).and.(gnd.gt.1.099)) gnd=1.099
      gnd2= gnd
      gnd3= gnd
          if(xmag.lt.cmagsig) csigma= csigma1+ 
     &     csigmacoef*xmag
          if(xmag.ge.cmagsig) csigma= csigma2
          if(dist0.lt.distnf) csigma= csigma+ sigmanf
      sigmasq= csigma*sqrt(2.)
      test= exp(gnd + 3.*sigmasq/sqrt(2.))
      if((clamp.lt.test).and.(clamp.gt.0.)) clamp2= alog(clamp)
      if(clamp.ge.test) clamp2= gnd+ 3.*sigmasq/sqrt(2.)
      if(clamp.eq.0.) clamp2= gnd+ 3.*sigmasq/sqrt(2.)
      tempgt3= (gnd- clamp2)/sigmasq
      probgt3= (eerf(tempgt3)+1.)/2.
      do 105 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-probgt3)/(1.-probgt3)
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
  105 continue
  104 continue
  103 continue
      return
      end
ccccccc
      subroutine getBJF97(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
      common / atten / p, phang, xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
        write(6,*) "enter b1,b2,b3,b4,b5,h,sigma,thrust,norm"
        read(1,*) b1,b2,b3,b4,b5,h,sigma,thrust,xnorm
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c--- Mmin=6.0 nmag=45 dmag=0.05
       pi= 4.*atan(1.)
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
        xmag= xmag0
        if(iconv(ip,ia).eq.1) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
c-- following for Joyner Boore WUS
          sig= sigma
          if(dist0.lt.distnf) sig= sig+ sigmanf
          sigmasq= sig*sqrt(2.)
          dist= sqrt(dist0*dist0+h*h)
          gnd= b1+b2*(xmag-6.)+b3*(xmag-6.)**2
          gnd= gnd+ b4*dist + b5*alog(dist)
          gnd2= gnd+ thrust
          gnd3= gnd+ xnorm
      do 102 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
  102 continue
  104 continue
  103 continue
      return
      end
cccccccc
      subroutine getMota(ip,period,ia,ndist,di,nmag,sigmanf,distnf)
      common / atten / p, phang, xlev, nlev, iconv
      dimension p(210,55,7,20,4,7),xlev(20,7),nlev(7),iconv(7,7)
      dimension phang(30,6,55,7,20,2)
c--  for Motazedian and Atkinson Puerto Rico relations
      write(6,*) "enter c1,c2,c3,c4,sigma"
      read(1,*) c1,c2,c3,c4,sigma
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di
        sig= sigma
        if(dist0.lt.distnf) sig= sig+ sigmanf
        sigmasq= sig*sqrt(2.)/0.4343
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= 5.8+(m-1)*0.05
c--- loop through atten. relations for each period
c-- gnd for SS; gnd2 for thrust; gnd3 for normal
        xmag= xmag0
        if(iconv(ip,ia).eq.1) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
        delta= -7.333 + 2.333*xmag
        dist= sqrt(dist0*dist0 + delta*delta)
        gnd= c1+c2*(xmag-6.)+c3*((xmag-6.)**2) + c4*dist
        if(dist.le.75.) fac= (-1.88+0.14*xmag)*alog10(dist)
        if((dist.gt.75.).and.(dist.le.100.))
     &  fac= (-1.88+0.14*xmag)*alog10(75.)
        if(dist.gt.100.) 
     &  fac= (-1.88+0.14*xmag)*alog10(75.) -0.5*alog10(dist/100.)
        gnd= gnd+ fac - alog10(980.)
        gnd= gnd/0.4343
c      write(15,*) period,dist0,xmag,exp(gnd)
        gnd2= gnd
        gnd3= gnd
      do 102 k=1,nlev(ip)
      temp1= (gnd- xlev(k,ip))/sigmasq
      temp2= (gnd2- xlev(k,ip))/sigmasq
      temp3= (gnd3- xlev(k,ip))/sigmasq
      p(ii,m,ip,k,1,ia)= (eerf(temp1)+1.)/2.
      p(ii,m,ip,k,1,ia)= (p(ii,m,ip,k,1,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,1,ia).lt.0.) p(ii,m,ip,k,1,ia)=0.
      p(ii,m,ip,k,2,ia)= (eerf(temp2)+1.)/2.
      p(ii,m,ip,k,2,ia)= (p(ii,m,ip,k,2,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,2,ia).lt.0.) p(ii,m,ip,k,2,ia)=0.
      p(ii,m,ip,k,3,ia)= (eerf(temp3)+1.)/2.
      p(ii,m,ip,k,3,ia)= (p(ii,m,ip,k,3,ia)-(1.3499e-3))/0.99865
      if(p(ii,m,ip,k,3,ia).lt.0.) p(ii,m,ip,k,3,ia)=0.
  102 continue
  104 continue
  103 continue
      return
  333 end
