c--- program hazgridXv3; last revised 10/17/02; 
c---- now handles 7 attenuation relations
c--- added getMota
c----- revised Atkinson and Boore intraslab
c--- changed index for Geomatrix (Youngs et al.) intraslab to 13
c---  added choice of mblg to Mw conversion
c---  revised Campbell 2000
c---  added BJF97
c--- iflt=3 and 4 uses mblg to Mw to calculate fault length
c----- option for AS for thrust or random
c---- fixed AB sub
c---- added Campbell CEUS
c---- added clamp for CEUS max gnd motion
c----   added distance dependent weighting
c--------added input depth for Geomatrix
c----  changed routine to find nearest distance
c---  intraslab
c--- with clipping for Toro and new Boore tables
c--- calculates mean annual rates of exceedances for gridded a-values
c--- can use b-value and Mmax matrix
c--- this version can use finite faults centered on grid cells for M>=6.5
c--- randomizes strike of faults
c--- this version can use multiple attenuation functions at different periods
c--- choice of attenuation relation: Joyner-Boore, Toro, Sadigh, Campbell, etc.
c--- written by Art Frankel  6/94, modified 5/95, 7/95, 10/95
c--- to run: hazgridX inputfile
c--- try: hazgridX WUSmap1.in
c--- output files have hazard curve for each site concatenated
c--- one output file per period
c--- gnd motion levels should be in units of g
c--- period=0 indicates PGA
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     &  wtdist(7,7)
      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,magmax,mmax(200000),magref
      integer*4 readn
      character namein*256,nameout*256,name*256,name3*256,name4*256
      dimension xlen2(50)
      dimension prob(1000,20,7),out(800000),
     &  iatten(7,7)
      dimension nattn(7),period(7)
      dimension ifp(7)
      dimension a(200000),asum(30,1000),b(200000),arate(200000,40)
      write(6,*) "enter name of input file"
 900  format(a256)
      call getarg(1,namein)
      open(unit=1,file=namein,status='old')
      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
      nrec= nx*ny
      write(6,*) "enter distance increment, max. distance"
      read(1,*) di,dmax
      pi= 3.14159
      write(6,*) "enter increments dlat,dlon for source cell"
      read(1,*) dsy, dsx
      write(6,*) "ymin=",ymin
      nsx= ceiling((xmax-xmin)/dsx +1)
      nsy= ceiling((ymax-ymin)/dsy +1)
      nsrc= nsx*nsy
      do 32 i=1,nsrc
      do 32 j=1,40
  32  arate(i,j)=0.
      write(6,*) "enter bval,magmin,magmax,dmag,magref"
      read(1,*) bval,magmin,magmax,dmag,magref
c--- iflt=1 uses finite faults for M>6.0
c--- iflt=2 uses finite faults and fixes strike
c--- iflt=3 uses finite faults with Johnston mblg to Mw conv.
c---- iflt=4 uses finite faults with Boore-Atkinson nblg to Mw conv.
c--- ibmat=1 uses b-value matrix
c--- maxmat=1 uses Mmax matrix
c-- set each to zero if you don't want these
      write(6,*) "enter iflt,ibmat,maxmat"
      read (1,*) iflt,ibmat,maxmat
      if(ibmat.eq.1) then
        write(6,*) "enter name of b-value file"
        read(1,900) name3
        call openr(name3)
        call getbuf2(b,nsrc,readn)
        write(6,*) "ymin=",ymin
        write(6,*) nsrc,readn
        endif
      if(maxmat.eq.1) then
        write(6,*) "enter name of Mmax file"
        read(1,900) name4
c       write(6,*) "Filename is: ",name4
        call openr(name4)
        call getbuf2(mmax,nsrc,readn)
        write(6,*) nsrc,readn
        endif
      dmag2= dmag/2.
      fac= alog10(10.**(bval*dmag2)-10.**(-bval*dmag2))
      nmagmax=0
      write(6,*) "enter name of a-value file"
      read(1,900) name
      call openr(name)
      call getbuf2(a,nsrc,readn)
      write(6,*) nsrc,readn
      write(6,*) "enter number of years for agrid file, conv to incr"
      read(1,*) cyr,incr
c--- convert from cumulative a-value to incremental a-value if necessary
      if (incr.eq.1) then
      do 112 j=1,nsrc
      if(a(j).ne.0.) then
         a(j)= alog10(a(j))+fac
         a(j)= 10.**a(j)
         endif
 112  continue
      endif
c----------- new line for fixing fault strike if iflt=2
      if(iflt.eq.2) read(1,*) fltstr
c--- center magnitude bins
      if(magmax.ne.magmin) magmin= magmin+dmag2
c-------------------------
c---- set up rate matrix (arate) for each source cell and mag increment
      do 90 j=1,nsrc
      if(a(j).eq.0.) go to 90
      if((ibmat.eq.1).and.(b(j).eq.0.)) b(j)=bval
      if(ibmat.eq.0) b(j)=bval
c---- some changes 5/29 for negative Mmax
c      if((maxmat.eq.1).and.(mmax(j).eq.0.)) mmax(j)=magmax
      if((maxmat.eq.1).and.(mmax(j).lt.0.)) a(j)=1.e-10
      if((maxmat.eq.1).and.(mmax(j).le.0.)) mmax(j)=magmax
      if(ibmat.eq.1)
     &  fac2= alog10(10.**(b(j)*dmag2)-10.**(-b(j)*dmag2))
      if(maxmat.eq.0) nmag= (magmax-dmag2-magmin)/dmag +1.4
      if((maxmat.eq.0).and.(magmax.eq.magmin)) nmag=1
      if(maxmat.eq.1) nmag= (mmax(j)-dmag2-magmin)/dmag +1.4
      if(nmag.gt.nmagmax) nmagmax=nmag
      do 91 m=1,nmag
      xmag= magmin+(m-1)*dmag
      if(ibmat.eq.0) rate= alog10(a(j)/cyr)-bval*xmag
      if(ibmat.eq.1) then
           rate= alog10(a(j)/cyr)-b(j)*xmag
           endif
      arate(j,m)= 10.**rate
  91  continue
  90  continue
      coef= 3.14159/180.
      ndist= dmax/di +0.5
      nmag= nmagmax
cccccccccccccccccccccccc
ccccccccccccccccccccccccc
      write(6,*) "enter number of periods"
      read(1,*) nper
      do 301 ii=1,ndist
      do 301 m=1,nmagmax
      do 301 ip=1,nper
      do 301 k=1,nlev(ip)
 301  p(ii,m,ip,k)=0.
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 403 k=1,nlev(ip)
 403  xlev(k,ip)= alog(xlev(k,ip))
      ndist= dmax/di + 0.5
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,magmin,dmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.2) 
     &  call getToro(ip,per,ia,ndist,di,nmag,magmin,dmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.3) 
     & call getSadigh(ip,per,ia,ndist,di,nmag,magmin,dmag,
     &  sigmanf,distnf)
      if(iatten(ip,ia).eq.4) 
     & call getCamp(ip,per,ia,ndist,di,nmag,magmin,dmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.5) 
     & call getAB95(ip,per,ia,ndist,di,nmag,magmin,dmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.6) 
     & call getFEA(ip,per,ia,ndist,di,nmag,magmin,dmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.7) 
     & call getSomer(ip,per,ia,ndist,di,nmag,magmin,dmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.8) 
     & call getAS(ip,per,ia,ndist,di,nmag,magmin,dmag,sigmanf,distnf)    
      if(iatten(ip,ia).eq.9) 
     & call getCamp2000(ip,per,ia,ndist,di,nmag,magmin,dmag,
     &  sigmanf,distnf)    
      if(iatten(ip,ia).eq.10)
     & call getCampCEUS(ip,per,ia,ndist,di,nmag,magmin,dmag,
     &  sigmanf,distnf)
      if(iatten(ip,ia).eq.11) 
     & call getBJF97(ip,per,ia,ndist,di,nmag,magmin,dmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.12)
     & call getABsub(ip,per,ia,ndist,di,nmag,magmin,dmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.13)
     & call getGeom(ip,per,ia,ndist,di,nmag,magmin,dmag,sigmanf,distnf)
      if(iatten(ip,ia).eq.14)
     & call getMota(ip,per,ia,ndist,di,nmag,magmin,dmag,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
c---- write header --------------
      do 291 ip=1,nper
      headr%name(1)= namein
      headr%name(2)= 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
      headr%extra(1)= cyr
      call puthead(ifp(ip),headr,ndata,readn)
 291  continue
c
c---loop through receiver sites
      do 30 j=1,1000
      do 30 ip=1,nper
      do 30 k=1,nlev(ip)
 30   prob(j,k,ip)= 0.
c----- uses Wells and Coppersmith relation to get length from moment mag.
c         arg= (xmag-5.08)/1.16
       do 37 i=1,nmagmax
         xmag= magmin+(i-1)*dmag
c---- Johnston 96 conversion
         if(iflt.eq.3) xmag=  1.14 + 0.24*xmag+0.0933*xmag*xmag
c----Boore Atkinson 87 conversion
         if(iflt.eq.4) xmag=  2.715 - 0.277*xmag+0.127*xmag*xmag
         arg= -3.22+0.69*xmag
         xlen= 10**arg
         xlen2(i)= xlen/2.
  37     continue
c---Here's the guts
      icnt=1
      do 100 i=1,nrec
      do 31 ii=1,1000
      do 31 jj=1,30
 31   asum(jj,ii)= 0
      iy= (i-1)/nx 
      ix= i-1-iy*nx
c      write(6,*) i,ix,iy
c      read(5,*) idum
      rx= xmin+(ix)*dx
      ry= ymax-(iy)*dy
c-- loop through source cells
c-- bin rate by distance and magnitude (asum)
      do 101 j=1,nsrc
      if(a(j).eq.0.) go to 101
      isy= (j-1)/nsx
      isx= j-1-isy*nsx
      sx= xmin+(isx)*dsx
      sy= ymax-(isy)*dsy
      dlen= 1.
      xdiff= abs(sx-rx)
      ydiff= abs(sy-ry)
      xdmax= dmax/(111.11*cos(ymax*coef))
      ydmax= dmax/111.11
      if((xdiff.gt.xdmax).or.(ydiff.gt.ydmax)) go to 101
      call delaz(sy,sx,ry,rx,dist,az,baz)
      if(dist.gt.dmax) go to 101
      idist= dist/di +1
      do 111 m=1,nmagmax
      xmag= magmin+(m-1)*dmag
      if((xmag.lt.6.0).or.(iflt.eq.0)) then
         asum(m,idist)= asum(m,idist) + arate(j,m)
         go to 111
         endif
      if((xmag.ge.6.0).and.(iflt.ge.1)) then
         if(iflt.eq.2) then
           azrand= fltstr
           go to 678
           endif
         rnd= rand(0)
c---- randomize strike of fault
         rnd= 2.*(rnd-0.5)
         azrand=  rnd*90.
  678    continue
c---- find closest distance to fault
         daz= coef*(az-azrand)
         sdist= abs(dist*sin(daz))
         cdist= abs(dist*cos(daz))
         if(cdist.le.xlen2(m)) then
          dmin2= sdist
         else
         dmin2= sqrt(sdist*sdist +(cdist-xlen2(m))**2)
         endif
         if(dmin2.gt.dmax) go to 111
         idist= dmin2/di +1
         asum(m,idist)= asum(m,idist) + arate(j,m)
         endif
 111  continue
 101  continue
c
      do 202 ii=1,ndist
      do 202 m=1,nmagmax
      do 202 ip=1,nper
      do 202 k= 1,nlev(ip)
 202  prob(icnt,k,ip)= prob(icnt,k,ip)+p(ii,m,ip,k)*asum(m,ii)
cc-- output hazard curves every 1000 sites
      if((icnt.eq.1000).or.(i.eq.nrec)) then
        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)
        readn= readn/4
        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)
      if((sorlat.eq.stnlat).and.(sorlon.eq.stnlon)) then
          delta=0.
          az=0.
          return
          endif
      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(20,30,20), rlog(30), f(20)
      character*60 header

c
c     Read AB94 simulated grd. motion table.

c

      open(14,file='ab94.tab',status='old')
 900  format(i2)

      do 20 j = 1, 5

        read(14,10)header
        read(14,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(14,900) idum
        read(14,900) idum
        read(14,35)xmag(jm)

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

        read(14,45)
        read(14,900) idum
        read(14,45)
        read(14,900) idum

        do 40 jd = 1, nd

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

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

40      continue

50    continue
      close(14)
      return
      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 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(*,*)' ERROR. MAGNITUDE OUTSIDE RANGE.'

      jmu = jml + 1

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

c     find bounding distance values.

      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(*,*)' 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
cccccccccccccccccccc
      subroutine getBJF93(ip,period,ia,ndist,di,nmag,
     & magmin,dmag,sigmanf,distnf)
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
        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
      pi = 4.*atan(1.)
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0+h*h)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin + (m-1)*dmag
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= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
c-- following for Joyner Boore WUS
          sig= sigma/0.4343
          if(dist0.lt.distnf) sig= sig+ sigmanf
          sigmasq= sig*sqrt(2.)
          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
c       if(period.eq.1.) write(12,*) ia, period, xmag,  dist0, 10.**gnd
          if(period.ne.0.) gnd=gnd+alog10(2.*pi/(980.*period))
          gnd= gnd/0.4343
      do 102 k=1,nlev(ip)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1- 1.3499e-3)/0.99865
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  102 continue
  104 continue
  103 continue
      return
      end
cccccccccccccccc
      subroutine getToro(ip,period,ia,ndist,di,nmag,
     &   magmin,dmag,sigmanf,distnf)
      real magmin
      common / atten / p, xlev, nlev, iconv, wt,wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7)
        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-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0+th*th)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin + (m-1)*dmag
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= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) 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.)
          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.75g 0.3, 3.75g 0.2 
          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)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1-probgt3)/(1.-probgt3)
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  105 continue
  104 continue
  103 continue
      return
      end
cccccccccccccccccc
      subroutine getSadigh(ip,period,ia,ndist,di,nmag,
     &  magmin,dmag,sigmanf,distnf)
      real magmin
      common / atten / p, xlev, nlev, iconv, wt,wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7)
         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-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0+ sh*sh)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin + (m-1)*dmag
        xmag= xmag0
        if(iconv(ip,ia).eq.1) xmag= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) 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(dist+ 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(dist+exp(sd5+sd6*xmag))
            endif
          gnd= gnd+sd7*alog(dist0+2.)
      do 102 k=1,nlev(ip)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1- 1.3499e-3)/0.99865
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  102 continue
  104 continue
  103 continue
      return
      end
cccccccccccccccccccccccc
      subroutine getCamp(ip,period,ia,ndist,di,nmag,
     c  magmin,dmag,sigmanf,distnf)
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
         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-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0+ ch*ch)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin + (m-1)*dmag
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= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) 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=  dist**2+ (cc4*exp(cc5*xmag))**2
          arg= sqrt(arg)
          gnd= cc1+cc2*xmag+cc3*alog(arg)
          if(dist.eq.0.) then
             gnd=0
             gnd2=0
             gnd3=0
             go to 707
             endif
          gnd= gnd+ cs1+ cs2*alog(dist)
c--- following for Campbell thrust fault
          gnd= gnd+ cd1+cd2*alog(dist0)+cd3*xmag
  707     continue
      do 102 k=1,nlev(ip)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1- 1.3499e-3)/0.99865
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  102 continue
  104 continue
  103 continue
      return
      end
cccccccccccccccccccccc
      subroutine getAB95(ip,period,ia,ndist,di,nmag,
     &   magmin,dmag,sigmanf,distnf)
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
      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-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0+ abdepth*abdepth)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin + (m-1)*dmag
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= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) 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.)
          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      write(12,*) period,dist,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
      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)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1-probgt3)/(1.-probgt3)
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  105 continue
  104 continue
  103 continue
      return
      end
ccccccccccccccccccccccc
      subroutine getFEA(ip,period,ia,ndist,di,nmag,
     &   magmin,dmag,sigmanf,distnf)
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
      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-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0+bdepth*bdepth)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin + (m-1)*dmag
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= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) 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.)
c          dist= dist0
          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
c          write(19,*) ip,xmag,imag,dist,idist
          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)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1-probgt3)/(1.-probgt3)
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  105 continue
  104 continue
  103 continue
      return
      end
ccccccccccccccccccccccccc
      subroutine getSomer(ip,period,ia,ndist,di,nmag,
     & magmin,dmag,sigmanf,distnf)
c---- Somerville et al (2001) for CEUS
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
        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-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      sig= sig0
      if(dist0.lt.distnf) sig= sig0+ sigmanf
      sigmasq= sig*sqrt(2.)
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0+ 6.*6.)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin + (m-1)*dmag
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= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) 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--- 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)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1-probgt3)/(1.-probgt3)
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  105 continue
  104 continue
  103 continue
      return
      end
ccccccccccccccccccccccccc
      subroutine getAS(ip,period,ia,ndist,di,nmag,
     &  magmin,dmag,sigmanf,distnf)
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
         write(6,*) "enter a1,a2,a3,a4,a5,a6,h,ithrust"
         read(1,*) as1,as2,as3,as4,as5,as6,h,ithrust
         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-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist1= sqrt(dist0*dist0 + h*h)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin + (m-1)*dmag
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= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) 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(dist1*dist1+ 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
          if(ithrust.eq.1) gnd= gnd+ fltfac
          if(ithrust.eq.2) gnd= gnd + 0.5*fltfac
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
      do 102 k=1,nlev(ip)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1- 1.3499e-3)/0.99865
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  102 continue
  104 continue
  103 continue
      return
      end
ccccccccccccccccccccccccccccc
      subroutine getCamp2000(ip,period,ia,ndist,di,nmag,
     &     magmin,dmag,sigmanf,distnf)
c  Campbell and Bozorgnia (2000)
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
      read(1,*) c1,c2,c3,c4,c5,c6,c7,c8,c9,h
      read(1,*) c10,c11,csite,ithrust
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-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0 + h*h)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin +(m-1)*dmag
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= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) 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---- below for reverse-thrust fault types
      if(ithrust.eq.1) gnd= gnd+ 0.5*c10 + 0.5*c11
c---- below for random fault type
      if(ithrust.eq.2) gnd= gnd+ 0.25*c10 + 0.25*c11
c      write(12,*) period, xmag, dist0, exp(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.)
      do 102 k=1,nlev(ip)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1- 1.3499e-3)/0.99865
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  102 continue
  104 continue
  103 continue
      return
      end
cccccccccccccccc
      subroutine getGeom(ip,period,ia,ndist,di,nmag,
     &     magmin,dmag,sigmanf,distnf)
c  Geomatrix (Youngs et al. intraslab)
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
        write(6,*) "enter c1,c2,c3,c4,c5,depth"
        read(1,*) gc1,gc2,gc3,gc4,gc5,gch
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0+ gch*gch)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin +(m-1)*dmag
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= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
          sig= gc4+gc5*xmag
          sigmasq= sig*sqrt(2.)
          gnd= 0.2418+0.3846+0.00607*gch+1.414*xmag+ 
     &   gc1+gc2*((10.-xmag)**3) 
          arg= 0.554*xmag
          gnd=gnd+gc3*alog(dist+1.7818*exp(arg))
      do 102 k=1,nlev(ip)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1- 1.3499e-3)/0.99865
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  102 continue
  104 continue
  103 continue
      return
      end
ccccccccccccccc
      subroutine getABsub(ip,period,ia,ndist,di,nmag,
     &     magmin,dmag,sigmanf,distnf)
c  Atkinson and Boore subduction zone intraslab
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
        write(6,*) "enter c1,c2,c3,c4,c5,sig,depth"
        read(1,*) c1,c2,c3,c4,c5,sig,depth
          sigmasq= sig*sqrt(2.)/0.4343
c set up erf matrix p as ftn of dist,mag,period,level,flt type,atten type
c-- loop through distances
      if(period.ne.0.) freq= 1./period
      if(period.eq.0.) freq= 100.
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0+ depth*depth)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin +(m-1)*dmag
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= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) xmag= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
          delta= 0.00724*(10.**(0.507*xmag))
          dist2= sqrt(dist*dist + delta*delta)
          g= 10.**(0.301-0.01*xmag)
          gnd= c1+c2*xmag+c3*depth+c4*dist2
     &         -g*alog10(dist2)
c--- calculate rock PGA for C site amp
c--- Cascadia inslab
       rc1= -0.25
       rc2= 0.6909
       rc3= 0.01130
       rc4= -0.00202
       rpga= rc1+ rc2*xmag + rc3*depth+ rc4*dist2- g*alog10(dist2)
       rpga= 10.**rpga
       if((rpga.le.100.).or.(freq.le.1.)) sl=1.
       if((rpga.gt.100.).and.(rpga.lt.500.).and.(freq.gt.1.).and.
     &   (freq.lt.2.)) sl= 1.-(freq-1.)*(rpga-100)/400.
       if((rpga.ge.500.).and.(freq.gt.1.).and.(freq.lt.2.))
     &   sl= 1.-(freq-1.)
       if((rpga.gt.100.).and.(rpga.lt.500.).and.(freq.ge.2.))
     &   sl= 1.-(rpga-100.)/400.
       if((rpga.ge.500.).and.(freq.ge.2.)) sl= 0.
c-----
c---   take log ave of B (rock) and C site
          gnd= gnd + (sl*c5)/2. - alog10(980.) 
          gnd= gnd/0.4343
c      write(15,*) period, xmag, dist, exp(gnd), rpga, sl
      do 102 k=1,nlev(ip)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1- 1.3499e-3)/0.99865
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  102 continue
  104 continue
  103 continue
      return
      end
cccccccc
      subroutine getCampCEUS(ip,period,ia,ndist,di,nmag,
     &     magmin,dmag,sigmanf,distnf)
c----- Campbell 2001 CEUS
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
      read(1,*) c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,h,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 + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0 + h*h)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin+(m-1)*dmag
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= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
        if(iconv(ip,ia).eq.2) 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
          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)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1-probgt3)/(1.-probgt3)
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  105 continue
  104 continue
  103 continue
      return
      end
cccccc
      subroutine getBJF97(ip,period,ia,ndist,di,nmag,
     & magmin,dmag,sigmanf,distnf)
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
        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
      pi = 4.*atan(1.)
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
      dist= sqrt(dist0*dist0+h*h)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin + (m-1)*dmag
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) then
           xmag1= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
           xmag2= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
           xmag= 0.5*(xmag1+xmag2)
           endif
c-- following for Joyner Boore WUS
          sig= sigma
          if(dist0.lt.distnf) sig= sig+ sigmanf
          sigmasq= sig*sqrt(2.)
          gnd= b1+b2*(xmag-6.)+b3*(xmag-6.)**2
          gnd= gnd+ b4*dist + b5*alog(dist)
      do 102 k=1,nlev(ip)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1- 1.3499e-3)/0.99865
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  102 continue
  104 continue
  103 continue
      return
      end
cccccc
      subroutine getMota(ip,period,ia,ndist,di,nmag,
     & magmin,dmag,sigmanf,distnf)
      real magmin
      common / atten / p, xlev, nlev, iconv, wt, wtdist
      dimension p(210,30,7,20),xlev(20,7),nlev(7),iconv(7,7),wt(7,7,2),
     & wtdist(7,7) 
c--  for Motazedian and Atkinson Puerto Rico relations
      write(6,*) "enter c1,c2,c3,c4,sigma,h"
      read(1,*) c1,c2,c3,c4,sigma,h
c-- loop through distances
      do 103 ii=1,ndist
      dist0= (ii-1)*di + di/2.
      dist0= sqrt(dist0*dist0 + h*h)
        sig= sigma
        if(dist0.lt.distnf) sig= sig+ sigmanf
        sigmasq= sig*sqrt(2.)/0.4343
      weight= wt(ip,ia,1)
      if(dist0.gt.wtdist(ip,ia)) weight= wt(ip,ia,2)
c-- loop through magnitudes
      do 104 m=1,nmag
      xmag0= magmin + (m-1)*dmag
        xmag= xmag0
        if(iconv(ip,ia).eq.1) then
           xmag1= 1.14 +0.24*xmag0+0.0933*xmag0*xmag0
           xmag2= 2.715 -0.277*xmag0+0.127*xmag0*xmag0
           xmag= 0.5*(xmag1+xmag2)
           endif
        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)
      do 102 k=1,nlev(ip)
      temp= (gnd- xlev(k,ip))/sigmasq
      temp1= (eerf(temp)+1.)/2.
      temp1= (temp1- 1.3499e-3)/0.99865
      if(temp1.lt.0.) temp1=0.
      p(ii,m,ip,k)= p(ii,m,ip,k) + weight*temp1
  102 continue
  104 continue
  103 continue
      return
      end
