      program GMPEtable
c  Control file for this program is GMPEtable.ctl
c  This program reads GMPE table specified in GMPEtable.ctl
c  e.g. (AB06RevA_Rcd.dat, A08RevA_Rjb.dat, P11A_Rcd.dat)
c  the tables contain predicted median ground motions in log10  cgs units
c  Also reads a list of magnitudes, distances, frequencies for which GMPEs will be evaluated.
c  Writes them out in convenient format(to GMPEtable.out) to import to excel or other spreadsheet and plot
c
c  Provided to USGS by Gail Atkinson, March 17, 2011
c  Modified for GSC Aug. 2011 to include option for ENA or WNA fault sizes
c
c  modified Feb. 2012: 1) to include option to read a list of freq-dep. factors
c   these log factors are added to the GMPE to create a "higher" or "lower" GMPE.
c    Note the factors must be given for the same freq. as requested for the output table.
c    There is also an option for a distance dep. additive factor (as well as freq. dep;  can use either or both)
c    The dist-dep. factor is independent of freq., given by constant and slope, and bounded low and high.
c
c    2) Also increase all array sizes for input tables to 60
c
c    3) Also add in the option for inslab area factors, using Strasser et al., 2010 SRL;
c     - the option for the inslab area fault size is invoked using fault size factor = -1.
c
c     modification to algorithm in geta subroutine so that it does not bomb in some compilers due to tiny diff.-Aug. 2012
      common xmag(60), gma(60,60,60), rlog(60), f(60), itype
      common nf, nd, nm
      common jfreq, freq(60)
      dimension FDEPTH(100), ampfac(60)
      dimension yout(60) 
      character*80 header
      character*30 fname
   
      open(5,file='GMPEtable.ctl')
      open(6,file='GMPEtable.out',status='unknown')
      read(5,3) header
      write(6,3) header
c     Read filename with the GMPE table data to be read.
      read(5,3) fname
c     Read the desired 1st magnitude, increment, and number of mags for table
      read(5,*) xm1, xmdelta, nmtable
c     Read same parameters for epicentral distance, but in log10 units.
c     Note we are specifying epicentral distance and will do approx. conversion from table's metrics.
      read(5,*) d1, ddelta, ndtable
c     Specify frequencies: an error will occur if not within table range.
      read(5,*) nftable
      read(5,*) (freq(i), i=1, nftable)
	 INGS=1
c     read focal depth for a single source, have set INGS=1
      read(5,*) FDEPTH(INGS)
c     Use isat=1 if we want to include recommended dist saturation for AB06 (0 otherwise)
      read(5,*) isat
c     Fault size factor relative to Wells&Coppersmith; use 1.0 for WNA, 0.6 for ENA
c     Use -1.  for factors of Strasser et al. for subduction zones
      read(5,*) sizefac
c     read logfactors to make bigger/smaller GMPE table if desired. (factors read but ignored if noptfac=0)
      read(5,*) noptfac
	read(5,*) (ampfac(i),i=1,nftable)
c     read factors to specify dist-dep. added sigma to GMPE table if desired.  
      read(5,*) noptfacD
      if (noptfacD .gt. 0) read(5,*) facDc, facDb, 
     *                                       facDmin, facDmax	
	write(*,*) ' have read GMPEtable.ctl'

   
c***************************************************************
c     Read grd. motion table.  2nd last col. is PGA(f=99.), last is PGV(f=89.)
c     itype tells us the distance metric of the table  (all magnitudes are moment M)
c
      open(3,file=fname)
	read(3,3)header
3     format(a)
	read(3,*) nm, nd, nf, itype
	write(*,9) fname, itype
9     format(' GMfile, itype(0=Repi,1=Rhypo,2=Rcd,3=Rjb): ',/, a30,i5)
      if (itype .eq. 0) write(6,*)
     *    ' NO DISTANCE CONVERSIONS DONE FROM INPUT Rtable!'
      if (itype .gt. 0) write(6,*)
     * ' FINITE-FAULT DISTANCE CONVERSIONS TO Repi APPLIED FROM Rtable!'

	read(3,*) (f(jf), jf=1,nf)
c	write(*,*) (f(jf), jf=1,nf)
c
      do 5000 jm = 1, nm
	  read(3,*)xmag(jm)
	  do 4000 jd = 1, nd
	    read(3,*)rlog(jd),(gma(jf,jd,jm), jf = 1,nf)
c	    write(*,*) xmag(jm), rlog(jd)
4000    continue
5000  continue
c****************************************************************
      write(6,9) fname, itype
      write(6,10) FDEPTH(INGS), isat
10    format(' Assumed focal depth = ', f8.1, ' isat:(1 for saturation)'
     *  ,i5)
      write(6,12) sizefac
12    format(' Fault size factor relative to W&C(-1 for inslab)= ',f5.2)
      if (noptfac .gt. 0) write(6,13) (ampfac(i),i=1,nftable)
13    format(' Additive logfactors:', 60f12.3)
      if (noptfacD .gt. 0) write(6,14) facDc, facDb, facDmin, facDmax
14    format(' Add Dist-factors (const, slope, minval, maxval): ', f6.3,
     *                                                  f8.5, 2f6.3) 
      write(6,*)'   Mw     Rtable   Repi  PSA(cgs) for freq:'
      write(6,20) (freq(i), i=1, nftable)
20    format(24x, 60f12.3)

c     Now loop over mag, dist, freq. to make table.
c     Note this is similar to loop inside PSHA.
          
      do 3000 im= 1, nmtable
        amag= xm1 + float(im-1)*xmdelta
        
        do 2000 id= 1, ndtable
c         recall that distance was specified in logs, so convert
          REPIS = 10**(d1 + float(id-1)*ddelta)
c	    write(*,*) id, repis
          
          do 1000 jfreq=1, nftable

c**   In the following, I call the interpolation subroutine.
c*******************************************************************************************
C**   This approach matches PSHA where software is based on a point source epicentre
c     where the stepping out is being done in epicentral distance(REPIS) 
c     FDEPTH(INGS) is the focal depth for the source zone considered
c     I assume calculations are being done for some specified frequency (freq(jfreq))
c***********************************************************************************
c     Here we first do a conversion from REPIS to the appropriate dist. measure
c     Wells and Coppersmith, 1984 fault length and width are:
c
      if (sizefac .gt. 0.) then
        faultl= 10**(-2.44 + 0.59*amag)
	  faultw= 10**(-1.01 + 0.32*amag)

c***  We can scale these by fault size factors (Somerville finds 0.6 for ENA)
        faultl=faultl*sizefac
	  faultw=faultw*sizefac

C**    then find the top of the fault, assuming a vertical fault
	  Dtop=FDEPTH(INGS) - 0.5*faultw
	else
c**     if negative sizefac is specified, used the Strasser formula for in-slab source at depth	
        Dtop = FDEPTH(INGS)
        area= 10**(-3.225 +0.89*amag)
        faultl= sqrt(area/3.141592)
c**     we divide this by 0.3 because we mult by 0.3 below 
        faultl = faultl/0.3
      endif
c***************************************************************************************
c**   This section is optional if you wish to impose a minimum
c**   top of rupture depth (Ztor) to force near-source saturation.
c**   I use this only for AB06' (but it usually doesn't matter much and could be skipped)

      if (isat .eq. 1) then
        Ztor= 21.-2.5*amag
       	if (Dtop .lt. Ztor) Dtop = Ztor
  	endif
c***************************************************************************************

	if (Dtop .lt. 0.) Dtop = 0.
c**   make approximation that surface distance is epicentral dist reduced by 0.3 *fault length
c**   this assumes epicentre in mid-fault, with random fault orientation
c**    (such that reduction varies from 0 to L/2 depending on azimuth; avg(by cos(theta)L/2)is 0.3L)
c**   I have tested this algorithm (for AB06) vs. using equiv. Rhypo gmpe, at 10km depth, for
c     moderate seismicity environment in Quebec, and it gives ~same answer within about 5% at 10^-4.
	Dsurf = REPIS - 0.3*faultl
	if (Dsurf .lt. 1.) Dsurf = 1.

c**   Select appropriate distance measure
      if (itype .eq. 0) R = REPIS
	if (itype .eq. 1)R = SQRT(REPIS*REPIS + FDEPTH(INGS)*FDEPTH(INGS)) 
	if (itype .eq. 2) R = sqrt(Dsurf*Dsurf + Dtop*Dtop)
	if (itype .eq. 3) R = Dsurf

c**************************************************************************
c     Call the interpolation subroutine with magnitude(amag), the appropriate distance metric(R).
c     It will return the median ground motion in log10 units (amean).
c     We are assuming we are in a loop over frequency (freq), and doing index jfreq
c     This can be changed if running a single frequency
c     We are placing the value freq(jfreq) within the table array of frequencies f(nf)

      call geta(amag,R,amean)
c      write(*,*) im, id, jfreq, amag, R, amean

c*************************************************************************
c     The returned amplitude at this freq. is amean.  Add log factors as applicable
      if (noptfac .gt. 0) amean = amean + ampfac(jfreq)
c      we need to add the factor ampfac at freq. specified in array freq(), if noptfac>0

      if (noptfacD .gt. 0) then
c      we need to add a distance-dependent factor.
       addDfac = facDc + facDb*R
       if (addDfac .lt. facDmin) addDfac = facDmin
       if (addDfac .gt. facDmax) addDfac = facDmax
       amean = amean + addDfac
      endif

c    Store the returned ground motion value for printing and do next freq.
c    Convert from logs for printing
      yout(jfreq) = 10.**amean
c	write(*,*) yout(jfreq)
      
1000  continue
c**   Done all frequencies for this magnitude, distance combo. Print.

      write(6,1010) amag, R, REPIS,(yout(jf), jf=1, nftable)
1010  format(3f8.2, 60(1pe12.3))

2000  continue
3000  continue
      stop
      end
      
c************************************************************
      subroutine geta(amag,r,amean)
c     Gets ground motion value (log PSA) from table 

      common xmag(60), gma(60,60,60), rlog(60), f(60), itype
      common nf, nd, nm
      common jfreq, freq(60)

      dimension avg(60)
c     
c     Use interpolation to get amean for given freq(jfreq), amag, R (hazard looping values).
c     Note that freq(jfreq) is the defined block of freq. for hazard calcs (common block)
c     GMPE tables give values(gma) for freq array f(nf), for distances rlog(nd), magnitudes xmag(nm)
c
      rl = alog10(R)
      jfl = 0
      jdflag = 0
c     Find bounding jf values.
      do 10 j = 1,nf-3
	if (freq(jfreq) .ge. f(j) .and. freq(jfreq) .lt. f(j+1))jfl = j
10    continue
c     check near the edges (in case match of frequencies not exact)
      if (abs(freq(jfreq)-f(1)) .lt. 0.01) jfl = 1
      if (abs(freq(jfreq)-f(nf)) .lt. 0.01) jfl = nf
      if (abs(freq(jfreq)-f(nf-1)) .lt. 0.01) jfl = nf-1
      if (abs(freq(jfreq)-f(nf-2)) .lt. 0.01)  jfl = nf-2

      if (jfl .eq. 0) write(*,*) ' ERROR. FREQ. NOT FOUND.'
      if (freq(jfreq) .gt. f(nf-2) .and. freq(jfreq) .lt. 89.)
     *   write(*,*)' ERROR. FREQ. OUTSIDE RANGE'
      jfu = jfl + 1
      fracf = 0.
      if (abs(freq(jfreq)-f(jfl)) .lt. 0.01) then
c**     no interpolation needed, fracf=0.
	  go to 15
	else
c**     interpolate if needed (redefine fracf)
        fracf = (alog10 (freq(jfreq)) - alog10 (f(jfl)) ) /
     *        (alog10 (f(jfu)) - alog10 (f(jfl)) )
	endif

c	write(*,*) jfl, jfu, fracf
c**   fracf gives interpolation fraction in frequency
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**   fracm give interpolation fraction in magnitude
c     find bounding distance values.
      if (rl .lt. rlog(1)) write(*,*)' ERROR. DIST. < RANGE '
	  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 (rl .gt. rlog(nd)) then
c*       dist. beyond table maximum.  Attenuate ground motion by 1/R from last point in table.
         rscale= rl - rlog(nd)
         jdflag= 1
         jdl = nd
	  endif

      if (jdl .eq. 0) write(*,*)' ERROR. CANNOT FIND DIST'
      jdu = jdl + 1
      fracd = (rl - rlog(jdl)) / (rlog(jdu)-rlog(jdl))
c**   fracd gives interpolation fraction in distance
      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))
      if (jdflag .eq. 1) amean = amean - rscale
      return
      end

