
c************************************************************
      subroutine geta(freq,amag,R,amean)
c     Gets ground motion value (log PSA) from table 
c     Note that common variables include:
c     array of magnitudes in table (xmag), dist in logs (rlog), freq (f)
c     with ground motion values assigned to gma(jf,jd,jm)
c     So need to have read in f(), xmag(), rlog(), and gma()
c     
      common xmag(20), gma(20,30,20), rlog(30), f(20)
      common nf, nd, nm
      dimension avg(20)
c     
c     Use interpolation to get amean for given freq, amag, R.
c
      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-f(nf-2)) .lt. 0.01)  jfl = nf-2
      if (jfl .eq. 0) write(*,*) ' ERROR. FREQ. OUTSIDE RANGE.'
      if (freq .gt. f(nf-2) .and. freq .lt. 89.)
     *   write(*,*)' ERROR. FREQ. OUTSIDE RANGE'
      jfu = jfl + 1
      fracf = 0.
      if (freq .ge. f(nf-2)) 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, nm-1
	 if (amag .ge. xmag(j) .and. amag .lt. xmag(j+1)) jml = j
20    continue
      if (abs(amag - xmag(nm)) .lt. 0.01) jml = nm
      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.
      if (rl .lt. rlog(1)) write(*,*)' ERROR. DIST. < RANGE '
	  if (rl .gt. rlog(nd)) then
c*      dist. beyond 500 km.  Assign negligible grd. motion
        amean = -5.
		return
	  endif
	  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(*,*)' ERROR. CANNOT FIND DIST'
      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

