c
c *** Simulated Acceleration on Engineering Bedock (2003, Apr.)***
c
      PROGRAM mwave2003
C
      PARAMETER(mm=33000)
      dimension time(mm),ddy(mm)
      dimension acc(mm),vel(mm),dis(mm),resa(mm,0:20),samp(mm)
      COMPLEX Fds(mm),Fds1(mm)
      CHARACTER INFILE*20
C
C  Prompt for input RECON'ed file
C
      write(*,*)
      write(*,*)'Simulation of Acceleration Wave on Engineering Bedrock'
	write(*,*)'  Coded by Y. Hisada (Kogakuin University)' 
	write(*,*)'  Based on "Manual for Seismic Inspection and Retrofitt 
	*ing (2003)"'
	write(*,*)
c
      WRITE(*,1010) 
 1010 FORMAT(1X,'Enter Acceleration Wave Data (ex: uni-mex85.acc)')
      READ(*,1020) INFILE
 1020 FORMAT(A20)
c
      write(*,*)
      write(*,*)'Choose Spectrum: Manual 2003 (=0) or Kokuji 1461 (=1)'
	read(*,*)iko
	if(iko.eq.1)then
	  write(*,*)'Choose Level 1 (Safety Level    : Enter=1)'
	  write(*,*)', or   Level 2 (Functional Level: Enter=2)'
        read(*,*)L12
	else
	  write(*,*)
	  write(*,*)'Choose Zone (A1 - B3: All Safety Level):'
	  write(*,*)
	  write(*,*)' Zone A1: Rock Site & South Kanto  (Enter:1)'
	  write(*,*)' Zone B1: Rock Site & North Kanto  (Enter:2)'
	  write(*,*)' Zone A2: Basin Edge & South Kanto (Enter:3)'
	  write(*,*)' Zone B2: Basin Edge & North Kanto (Enter:4)'
	  write(*,*)' Zone A3: Basin Site & South Kanto (Enter:5)'
	  write(*,*)' Zone B3: Basin Site & North Kanto (Enter:6)'
	  write(*,*)
	  write(*,*)'Enter Zone (A1=1,B1=2,A2=3,B2=4,A3=5,or B3=6)'
        read(*,*)izone
c
        write(*,*)
       write(*,*)'Consider Amplification of Surface Layer? (Y=0 or N=1)'
	  read(*,*)isr
	  SVs=0.0
	  STk=0.0
	  if(isr.eq.0)then
          write(*,*)
          write(*,*)'  Choose Average Vs of Surface Layer (m/s)'
	    write(*,*)'    Vs < 100 (Enter=1)'
	    write(*,*)'    Vs = 150 (Enter=2)'
	    write(*,*)'    Vs = 200 (Enter=3)'
 	    write(*,*)'    Vs > 250 (Enter=4)'
          read(*,*)iVs
          write(*,*)
          write(*,*)'  Choose Thickness of Surface Layer (m)'
	    write(*,*)'    H < 10 (Enter=1)'
	    write(*,*)'    H = 15 (Enter=2)'
	    write(*,*)'    H = 20 (Enter=3)'
 	    write(*,*)'    H = 25 (Enter=4)'
 	    write(*,*)'    H > 30 (Enter=5)'
          read(*,*)iTk
	  end if
c
	end if
c
      write(*,*)
   10 write(*,*)'Enter Number of Iteration Numbers (up to 20)'
      read(*,*)nit
      if(nit.gt.20)goto 10
C      
C *** Read input data file  ***
C
      OPEN(11,FILE=INFILE)
        i=1
   20   read(11,*,end=21)time(i),ddy(i)
        i=i+1
        goto 20
   21 CLOSE(11)
c
      NT=i-1
      dt=time(2)-time(1)
      write(*,*)'NT and DT=',NT,DT
c
C *** Find the Closest Number of 2**N ***
C
      NF=2
      do I=2,15
        NF=2*NF
        IF(NF.GE.NT)GOTO 30
      end do
      WRITE(*,*)' NT>2**15 !!!!!'
      STOP
   30 IF1=I
      WRITE(*,*)
      WRITE(*,*)' IF1,NF=',IF1,NF
c
C *** FOURIER TRANSFORM of the Original Wave ***
C
      do I=1,NT
        Fds(I)=CMPLX(ddy(I),0.0)
      end do
      do I=NT+1,NF
        Fds(I)=(0.0,0.0)
      end do
C
      CALL WVFAST(NF,Fds,mm,-1)
      DUR=dFLOAT(NF)*DT
      DOM=6.283185307/DUR
      NOM=NF/2+1
      EOM=DOM*NOM
      WRITE(*,*)
      WRITE(*,*)' DUR,DOM,NOM,EOM=',DUR,DOM,NOM,EOM
      WRITE(*,*)' Highest Frequency (hz)=',EOM/6.2831853070d0
      do I=1,NF
        Fds(I)=Fds(I)*DT
      end do
c
c *** Response Spectrum of Kokuji No.1461 (2000) ***
c
      h=0.05
      if(iko.eq.1)then

        do j=2,nom
	    samp(j)=1.0
        end do

        OM1=6.2831853070/0.64
        OM2=6.2831853070/0.16
        do j=2,nom
          OMEGA=FLOAT(j-1)*DOM
          period=6.283185307/omega
          resa(j,0)=800.0
          IF(OMEGA.LT.OM1)resa(j,0)=512.0/period
          IF(OMEGA.GT.OM2)resa(j,0)=period*3000.0+320.0
          if(L12.eq.2)resa(j,0)=resa(j,0)/5.0
        end do
c
      else
c
c *** Response Spectrum by Manual 2003 at Various Zones ***
c
c      * Calculation of Amplification Factor of Surface Layer *
c
	  if(isr.eq.1)then
c
          do j=2,nom
	      samp(j)=1.0
	    end do
c
	  else
c
c       * Zone 1 (Rock Site) *
          if((izone.eq.1).or.(izone.eq.2))then
            do j=2,nom
              OMEGA=FLOAT(j-1)*DOM
              period=6.283185307/omega
	        call samp1(period,iVs,iTk,amp1)
	        samp(j)=amp1
	      end do
c
	    else
c
c         * Zone 2 (Basin Edge Site) *
            if((izone.eq.3).or.(izone.eq.4))then
              do j=2,nom
                OMEGA=FLOAT(j-1)*DOM
                period=6.283185307/omega
	          call samp2(period,iVs,iTk,amp2)
	          samp(j)=amp2
	        end do
	      else
c         * Zone 3 (Basin Site) *
              do j=2,nom
                OMEGA=FLOAT(j-1)*DOM
                period=6.283185307/omega
	          call samp3(period,iVs,iTk,amp3)
	          samp(j)=amp3
	        end do
	      end if
c
	    end if
c
	  end if
c
	  if(izone.eq.1)then
          OM1=6.2831853070/0.5
          OM2=6.2831853070/0.08
          do j=2,nom
            OMEGA=FLOAT(j-1)*DOM
            period=6.283185307/omega
            resa(j,0)=800.0
            IF(OMEGA.LT.OM1)resa(j,0)=400.0/period
            IF(OMEGA.GT.OM2)resa(j,0)=period*6000.0+320.0
            resa(j,0)=resa(j,0)*samp(j)
          end do
	  else
	    if(izone.eq.2)then
            OM1=6.2831853070/0.5
            OM2=6.2831853070/0.08
            do j=2,nom
              OMEGA=FLOAT(j-1)*DOM
              period=6.283185307/omega
              resa(j,0)=800.0
              IF(OMEGA.LT.OM1)resa(j,0)=400.0/period
              IF(OMEGA.GT.OM2)resa(j,0)=period*6000.0+320.0
              resa(j,0)=resa(j,0)*samp(j)*0.9
            end do
	    else
	      if(izone.eq.3)then
              OM1=6.2831853070/0.64
              OM2=6.2831853070/0.16
              do j=2,nom
                OMEGA=FLOAT(j-1)*DOM
                period=6.283185307/omega
                resa(j,0)=800.0
                IF(OMEGA.LT.OM1)resa(j,0)=560.0/period**0.8
                IF(OMEGA.GT.OM2)resa(j,0)=period*3000.0+320.0
                resa(j,0)=resa(j,0)*samp(j)
              end do
	      else
	        if(izone.eq.4)then
                OM1=6.2831853070/0.64
                OM2=6.2831853070/0.16
                do j=2,nom
                  OMEGA=FLOAT(j-1)*DOM
                  period=6.283185307/omega
                  resa(j,0)=800.0
                  IF(OMEGA.LT.OM1)resa(j,0)=560.0/period**0.8
                  IF(OMEGA.GT.OM2)resa(j,0)=period*3000.0+320.0
                  resa(j,0)=resa(j,0)*samp(j)*0.9
                end do
	        else
	          if(izone.eq.5)then
                  OM1=6.2831853070/0.64
                  OM2=6.2831853070/0.16
                  do j=2,nom
                    OMEGA=FLOAT(j-1)*DOM
                    period=6.283185307/omega
                    resa(j,0)=800.0
                    IF(OMEGA.LT.OM1)resa(j,0)=585.0/period**0.7
                    IF(OMEGA.GT.OM2)resa(j,0)=period*3000.0+320.0
                    resa(j,0)=resa(j,0)*samp(j)
                  end do
	          else
                  OM1=6.2831853070/0.64
                  OM2=6.2831853070/0.16
                  do j=2,nom
                    OMEGA=FLOAT(j-1)*DOM
                    period=6.283185307/omega
                    resa(j,0)=800.0
                    IF(OMEGA.LT.OM1)resa(j,0)=585.0/period**0.7
                    IF(OMEGA.GT.OM2)resa(j,0)=period*3000.0+320.0
                    resa(j,0)=resa(j,0)*samp(j)*0.9
                  end do
	          end if
	        end if
	      end if
	    endif
        end if
	endif

C
C *** Modification of Original Acceleration ***
c
      write(*,*)
      do 100 iit=1,nit
	  err=0.0
        Fds(1)=(0.0,0.0)
        Fds1(1)=(0.0,0.0)
        do j=2,nom
          OMEGA=FLOAT(j-1)*DOM
          period=6.283185307/omega
          call resp(h,omega,dt,nt,ddy,acc,vel,dis,mm,sa,sv,sd)
          resa(j,iit)=sa
	    err=err+(resa(j,iit)-resa(j,0))/resa(j,0)
          Fds(j)=Fds(j)*resa(j,0)/resa(j,iit)
          Fds1(j)=Fds(j)
        end do
	  err=err/nom*100.0
        write(*,*)' Iteration Number =',iit,', Error (%) =',err
C
        do I=2,NF/2
           Fds1(NF+2-I)=CMPLX(REAL(Fds1(I)),-IMAG(Fds1(I)))
        end do
c
        CALL WVFAST(NF,Fds1,mm,1)
        do I = 1,NF
           ddy(I)=REAL(Fds1(I))/DUR
           t=dt*(i-1)
           if(t.lt.4.0)ddy(I)=ddy(I)*t*t/16.0
        end do
c
  100 CONTINUE
C
c *** output ***
c
      write(*,*)
      write(*,*)'Response Spectra are stored in "respa.csv"'
      write(*,*)'Data Order:'
      write(*,*)'Period, Amp., Manual or Kokuji, Manual or Kokuji * Amp,
     * Resp 1,2,..'
      OPEN(11,FILE='respa.csv',STATUS='unknown')
      do j=nom,2,-1
        OMEGA=FLOAT(j-1)*DOM
        period=6.283185307/omega
        call resp(h,omega,dt,nt,ddy,acc,vel,dis,mm,sa,sv,sd)
        write(11,2000)period,samp(j),resa(j,0)/samp(j)
	*               ,(resa(j,i),i=0,nit),sa
      end do
      CLOSE(11)
 2000 format(f8.3,',',23(f10.3,','))
c
      WRITE(*,2010) 
 2010 FORMAT(/1X,'Enter Output File Name of Simulation Wave (***.csv)')
      READ(*,2020) INFILE
 2020 FORMAT(A20)
      OPEN(12,FILE=INFILE,STATUS='unknown')
      do i=1,nt
      write(12,2100)dt*(i-1),ddy(I)
      end do
 2100 format(f10.4,f10.3)
      CLOSE(12)
c
        STOP
	END
C
C *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE WVFAST(N,A,MOM,IND)
C
C *** N   : NUMBER OF FOURIER COEFF. OR WAVE ( = 2**I )
C *** A(I): FOURIER COEFFICIENTS ( N/2+1 ) OR WAVE ( REAL PART )
C    ( FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0 )
C *** MOM : N IN MAIN DIMMENSION
C         : IF mm > 2**13(=8192), CAN NOT USE !!!
C *** IND : = -1 : TRANS.,   = 1 : INVERSE TRANS.
C
      COMPLEX A(MOM),TEMP,THETA
C
C *** CALUCURATION OF FFT
C
      PAI=3.141592654
      J=1
      DO 140 I=1,N
        IF(I.GE.J) GO TO 110
        TEMP=A(J)
        A(J)=A(I)
        A(I)=TEMP
  110   M=N/2
  120   IF(J.LE.M) GO TO 130
        J=J-M
        M=M/2
        IF(M.GE.2) GO TO 120
  130   J=J+M
  140 CONTINUE
      KMAX=1
  150 IF(KMAX.GE.N) RETURN
      ISTEP=KMAX*2
      DO 170 K=1,KMAX
        THETA=CMPLX(0.0,PAI*FLOAT(IND*(K-1))/FLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=A(J)*CEXP(THETA)
          A(J)=A(I)-TEMP
          A(I)=A(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END
c
c
c	Subroutine for Response of Single-Degree of System
c
	subroutine resp(h,w,dt,nn,ddy,acc,vel,dis,nd,sa,sv,sd)
c
	dimension ddy(nd),acc(nd),vel(nd),dis(nd)
c
	w2=w*w
	hw=h*w
	wd=w*sqrt(1.-h*h)
	wdt=wd*dt
	e=exp(-hw*dt)
	cwdt=cos(wdt)
	swdt=sin(wdt)
	a11= e*(cwdt+hW*swdt/wd)
	a12= e*swdt/wd
	a21=-e*w2*swdt/wd
	a22= e*(cwdt-hw*swdt/wd)
	ss=-hw*swdt-wd*cwdt
	cc=-hw*cwdt+wd*swdt
	s1=(e*ss+wd)/w2
	c1=(e*cc+hw)/w2
	s2=(e*dt*ss+hw*s1+wd*c1)/w2
	c2=(e*dt*cc+hw*c1-wd*s1)/w2
	s3=dt*s1-s2
	c3=dt*c1-c2
	b11=-s2/wdt
	b12=-s3/wdt
	b21=(hw*s2-wd*c2)/wdt
	b22=(hw*s3-wd*c3)/wdt
	acc(1)=2.*h*w*ddy(1)*dt
	vel(1)=-ddy(1)*dt
	dis(1)=0.
	dx=vel(1)
	x=0
	sa=0
	sv=0	
	sd=0
c
	do 110 m=2,nn
	dxf=dx
	xf=x
	ddym=ddy(m)
	ddyf=ddy(m-1)
	x = a12*dxf+a11*xf+b12*ddym+b11*ddyf
	dx= a22*dxf+a21*xf+b22*ddym+b21*ddyf
	ddx=-2.*hw*dx-w2*x
	acc(m)=ddx
	vel(m)=dx
	dis(m)=x
	sa=amax1(sa,abs(ddx))
	sv=amax1(sv,abs(dx))
	sd=amax1(sd,abs(x))
  110   continue
	return
	end
c
c
c	Subroutine for Amplification of Surface Layer (Zone 1: Rock Site)
c
	subroutine samp1(period,iVs,iTk,samp)
c
c   * VS = 100m/s *
      if(iVs.eq.1)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.463
	    T1=0.797
	    T2=0.957
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.48*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.3
	        else
	          samp=1.2*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.655
	    T1=1.126
	    T2=1.406
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.30*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.3
	        else
	          samp=2.4*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 20 m *
	  if(iTk.ge.3)then
	    T0=0.961
	    T1=1.592
	    T2=2.084
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.10*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.2
	        else
	          samp=4.5*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
c   * VS = 150m/s *
      if(iVs.eq.2)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.228
	    T1=0.378
	    T2=0.504
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.85*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.2
	        else
	          samp=0.35*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.413
	    T1=0.658
	    T2=0.855
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.54*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.1
	        else
	          samp=0.83*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 20 m *
	  if(iTk.eq.3)then
	    T0=0.520
	    T1=0.829
	    T2=1.050
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.42*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.1
	        else
	          samp=1.20*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 25 m *
	  if(iTk.eq.4)then
	    T0=0.667
	    T1=1.064
	    T2=1.394
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.29*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.1
	        else
	          samp=2.00*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 30 m *
	  if(iTk.eq.5)then
	    T0=0.808
	    T1=1.237
	    T2=1.664
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.19*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.0
	        else
	          samp=2.50*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
c   * VS = 200m/s *
      if(iVs.eq.3)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.155
	    T1=0.232
	    T2=0.332
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.89*period
	      else
	        if(period.lt.T2)then
	          samp=1.8
	        else
	          samp=0.11*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.251
	    T1=0.397
	    T2=0.563
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.68*period
	      else
	        if(period.lt.T2)then
	          samp=1.9
	        else
	          samp=0.32*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 20 m *
	  if(iTk.eq.3)then
	    T0=0.338
	    T1=0.535
	    T2=0.761
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.55*period
	      else
	        if(period.lt.T2)then
	          samp=1.9
	        else
	          samp=0.55*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 25 m *
	  if(iTk.eq.4)then
	    T0=0.478
	    T1=0.717
	    T2=1.014
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.40*period
	      else
	        if(period.lt.T2)then
	          samp=1.8
	        else
	          samp=0.82*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 30 m *
	  if(iTk.eq.5)then
	    T0=0.561
	    T1=0.842
	    T2=1.206
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.33*period
	      else
	        if(period.lt.T2)then
	          samp=1.8
	        else
	          samp=1.12*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
c   * VS = 250m/s *
      if(iVs.eq.4)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.116
	    T1=0.171
	    T2=0.212
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.92*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.7
	        else
	          samp=0.05*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.194
	    T1=0.267
	    T2=0.388
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.72*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.12*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 20 m *
	  if(iTk.eq.3)then
	    T0=0.278
	    T1=0.382
	    T2=0.583
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.58*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.24*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 25 m *
	  if(iTk.eq.4)then
	    T0=0.368
	    T1=0.506
	    T2=0.752
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.47*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.37*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 30 m *
	  if(iTk.eq.5)then
	    T0=0.463
	    T1=0.638
	    T2=0.919
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.38*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.52*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
	return
	end
c
c
c	Subroutine for Amplification of Surface Layer (Zone 2: Basin Edge Site)
c
	subroutine samp2(period,iVs,iTk,samp)
c
c   * VS = 100m/s *
      if(iVs.eq.1)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.463
	    T1=0.797
	    T2=0.965
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.48*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.3
	        else
	          samp=1.22*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.734
	    T1=1.217
	    T2=1.569
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.24*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.2
	        else
	          samp=2.7*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 20 m *
	  if(iTk.ge.3)then
	    T0=1.058
	    T1=1.686
	    T2=2.319
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.05*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.1
	        else
	          samp=5.0*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
c   * VS = 150m/s *
      if(iVs.eq.2)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.241
	    T1=0.385
	    T2=0.570
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.82*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.1
	        else
	          samp=0.40*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.429
	    T1=0.657
	    T2=0.902
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.52*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.0
	        else
	          samp=0.83*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 20 m *
	  if(iTk.eq.3)then
	    T0=0.561
	    T1=0.859
	    T2=1.220
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.38*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.0
	        else
	          samp=1.43*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 25 m *
	  if(iTk.eq.4)then
	    T0=0.749
	    T1=1.146
	    T2=1.569
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.23*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.0
	        else
	          samp=2.25*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 30 m *
	  if(iTk.eq.5)then
	    T0=1.018
	    T1=1.493
	    T2=2.110
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.07*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=1.9
	        else
	          samp=3.45*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
c   * VS = 200m/s *
      if(iVs.eq.3)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.158
	    T1=0.237
	    T2=0.332
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.88*period
	      else
	        if(period.lt.T2)then
	          samp=1.8
	        else
	          samp=0.11*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.263
	    T1=0.394
	    T2=0.601
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.66*period
	      else
	        if(period.lt.T2)then
	          samp=1.8
	        else
	          samp=0.32*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 20 m *
	  if(iTk.eq.3)then
	    T0=0.407
	    T1=0.610
	    T2=0.844
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.47*period
	      else
	        if(period.lt.T2)then
	          samp=1.8
	        else
	          samp=0.59*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 25 m *
	  if(iTk.eq.4)then
	    T0=0.512
	    T1=0.768
	    T2=1.054
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.37*period
	      else
	        if(period.lt.T2)then
	          samp=1.8
	        else
	          samp=0.88*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 30 m *
	  if(iTk.eq.5)then
	    T0=0.601
	    T1=0.902
	    T2=1.287
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.30*period
	      else
	        if(period.lt.T2)then
	          samp=1.8
	        else
	          samp=1.26*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
c   * VS = 250m/s *
      if(iVs.eq.4)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.113
	    T1=0.167
	    T2=0.212
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.93*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.7
	        else
	          samp=0.05*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.194
	    T1=0.267
	    T2=0.388
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.72*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.12*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 20 m *
	  if(iTk.eq.3)then
	    T0=0.292
	    T1=0.402
	    T2=0.611
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.56*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.26*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 25 m *
	  if(iTk.eq.4)then
	    T0=0.387
	    T1=0.533
	    T2=0.788
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.45*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.40*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 30 m *
	  if(iTk.eq.5)then
	    T0=0.488
	    T1=0.671
	    T2=0.950
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.36*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.55*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
	return
	end

c
c
c	Subroutine for Amplification of Surface Layer (Zone 3: Basin Site)
c
	subroutine samp3(period,iVs,iTk,samp)
c
c   * VS = 100m/s *
      if(iVs.eq.1)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.445
	    T1=0.782
	    T2=0.961
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.49*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.3
	        else
	          samp=1.21*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.793
	    T1=1.264
	    T2=1.826
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.20*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.1
	        else
	          samp=3.25*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 20 m *
	  if(iTk.ge.3)then
	    T0=1.058
	    T1=1.686
	    T2=2.345
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.05*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.1
	        else
	          samp=5.10*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
c   * VS = 150m/s *
      if(iVs.eq.2)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.241
	    T1=0.385
	    T2=0.570
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.82*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.1
	        else
	          samp=0.40*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.446
	    T1=0.683
	    T2=0.943
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.50*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.0
	        else
	          samp=0.90*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 20 m *
	  if(iTk.eq.3)then
	    T0=0.572
	    T1=0.876
	    T2=1.253
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.37*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=2.0
	        else
	          samp=1.50*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 25 m *
	  if(iTk.eq.4)then
	    T0=0.793
	    T1=1.163
	    T2=1.643
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.20*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=1.9
	        else
	          samp=2.20*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 30 m *
	  if(iTk.eq.5)then
	    T0=1.078
	    T1=1.581
	    T2=2.127
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.04*period**1.2
	      else
	        if(period.lt.T2)then
	          samp=1.9
	        else
	          samp=3.50*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
c   * VS = 200m/s *
      if(iVs.eq.3)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.162
	    T1=0.256
	    T2=0.326
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.87*period
	      else
	        if(period.lt.T2)then
	          samp=1.9
	        else
	          samp=0.12*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.269
	    T1=0.425
	    T2=0.582
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.65*period
	      else
	        if(period.lt.T2)then
	          samp=1.9
	        else
	          samp=0.34*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 20 m *
	  if(iTk.eq.3)then
	    T0=0.416
	    T1=0.659
	    T2=0.791
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.46*period
	      else
	        if(period.lt.T2)then
	          samp=1.9
	        else
	          samp=0.59*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 25 m *
	  if(iTk.eq.4)then
	    T0=0.512
	    T1=0.811
	    T2=1.030
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.37*period
	      else
	        if(period.lt.T2)then
	          samp=1.9
	        else
	          samp=0.95*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 30 m *
	  if(iTk.eq.5)then
	    T0=0.601
	    T1=0.952
	    T2=1.278
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.30*period
	      else
	        if(period.lt.T2)then
	          samp=1.9
	        else
	          samp=1.40*period**(-1.8)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
c   * VS = 250m/s *
      if(iVs.eq.4)then
c
c     * H = 10 m *
	  if(iTk.eq.1)then
	    T0=0.116
	    T1=0.171
	    T2=0.212
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.92*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.7
	        else
	          samp=0.05*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 15 m *
	  if(iTk.eq.2)then
	    T0=0.194
	    T1=0.267
	    T2=0.388
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.72*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.12*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 20 m *
	  if(iTk.eq.3)then
	    T0=0.300
	    T1=0.413
	    T2=0.639
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.55*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.28*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H = 25 m *
	  if(iTk.eq.4)then
	    T0=0.418
	    T1=0.576
	    T2=0.833
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.42*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.44*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
c     * H > 30 m *
	  if(iTk.eq.5)then
	    T0=0.500
	    T1=0.689
	    T2=0.980
	    if(period.lt.T0)then
	      samp=1.2
	    else
	      if(period.lt.T1)then
	        samp=10**0.35*period**0.9
	      else
	        if(period.lt.T2)then
	          samp=1.6
	        else
	          samp=0.58*period**(-1.7)+1
	        end if
	      end if
	    end if
	  end if
c
      end if
c
	return
	end
