c version aug 10, 2007 c scenario: program to interpolate hazard results for a specific probability c of exceedance at a user specified location. c The input files c for this consist of psa and pga in units of percent g. The c input files are in the following format: c longitude, lat, 0.5hz, 1.0hz, 2hz, 3.33hz, 5hz, 6.67hz, 13hz, pga. c for 1247 sites on a 1/8 x 1/8 degree grid. c example input file is in directory scdot/results/2percent_soil.dat c c An input file corresponding to 2%,5%, 7% and 10& probability of exceedance, for c both soil (geologically realistic site condition) and hard-rock basement outcrop c must exist in the same directory as the executable. Also, a file with the lat, lon c and thickness of coastal plain sediments in meters for the very same 1247 sites c must exist in the same directory. These files are the following c************************************************************************************ c NECESSARY INPUT FILES WHICH MUST BE IN THE SAME DIRECTORY AS THE EXECUTABLE c sites.dat ! sediment thickness in meters. c scenario-s2, scenario-r2 ! results for 2%, soil and rock c scenario-s5, scenario-r5 ! 5% c scenario-s7, scenario-r7 c scenario-s10, scenario-r10 !10% c scenario-sites ! sediment thickness in meters c c these files are just renamed versions of c 2percent_rock.dat 2percent_soil.dat ! hazard results for 2% probability in 50Y c 5percent_rock.dat 5percent_soil.dat c 7percent_rock.dat 7percent_soil.dat c 10percent_rock.dat 10percent_soil.dat c sites.dat c*********************************************************************************** c c the user is prompted to determine which file is read by the program. c c the user has the option to generate a scaled time series to match either c the whole uhs or a single frenquency or pga value. c c this option requres determination of the thickness of coastal plain sediments at the c chosen site location. The input file "sites.dat" must exist in the same directory as c the executable. c c Basic input time c series is defined on the basis of a user specified c modal magnitude and distance. The stochastic model of Boore and c Atkinson (1985) is used (w sq, hanks fmax, etc.). The SCDOT c geological structure (1 or 2 layer model) is used to modify the c rock motion, for the user defined site location. c c The time series can be modified to match a uniform c hazard spectrum using phase invariant spectral matching. c Or, the option exists to scale to either pga or a particular c psa oscillator response value (single frequency). c c c c the length of the time series is hardwired to 4096 points. c c martin chapman, december 2002,original c modified Feb. 2006 c USE MSFLIB DIMENSION FREQ(16384),DATA(32768) INTEGER IPOW2(14) DIMENSION RR(16384),WIN(16384),A(16384) COMPLEX SOIL,ROCK,RESULT character*1 yesno1,yesno2,yesno3,yesno4 character*66 filein character*66 filein2 character*66 filename real target(1000),rtarfrq(1000),uhs(8),ff(99),tt(4096),hazeq(6) real shake(4096),bsngl(4096),asngl(4096),sdtemp(1000),osfreq(8) real sv(100),sbeta(99),sh(99) character dirname*20 integer dnlen=0 data osfreq/0.5,1.0,2.0,3.33,5.0,6.67,13.0,-1./ c write(*,*)' ' write(*,*)' WELCOME TO PROGRAM SCENARIO_PC' write(*,*)' DESIGNED FOR SCDOT' write(*,*)' MARTIN CHAPMAN - VIRGINIA TECH' write(*,*)' ' write(*,*)' ' c write(*,*)'DIRECTORY FOR OUTPUT FILE LOCATION' write(*,*)'Enter name of new folder for placement of output' write(*,*)' files produced by this run.' write(*,*)' (Be sure that no folder by this name already exists.)' write(*,*)' (Folder name must be less than 16 characters long.)' read(*,*)dirname c determine the length of the folder name just entered do 6000 i=1,20 if(dirname(i:i).eq.' ') then dnlen=i-1 goto 5999 end if 6000 continue 5999 if(.not.makedirqq(dirname)) then c 5999 if(system('mkdir '//dirname).ne.0) then write(*,*)'Could not create directory -- does it already exist?' write(*,*)'Press Enter to terminate program.' read(*,*) stop end if write(*,*) write(*,*)' ALL OUTPUT FROM THIS RUN WILL BE PLACED IN: ',dirname write(*,*) c c summary of output filename=dirname(1:dnlen)//'/summary.txt' open(unit=3,name=filename,type='new',err=9000) write(3,*)'THIS FILE CONTAINS THE RESULTS FROM ONE EXECUTION' write(3,*)' OF PROGRAM scenario_pc (Martin Chapman, 2006).' write(3,*) write(3,*)'THE NAME OF THE DIRECTORY CONTAINING THIS FILE' write(3,*)' AND ALL ASSOCIATED OUTPUT FILES IS: ',dirname write(3,*)' ' write(3,*)' ' c write(*,*)' ' write(*,*)'SITE LOCATION' write(*,*)'Enter the latitude and longitude of the site' write(*,*)' of interest, in decimal degrees.' write(*,*)' Acceptable Latitude Range: 31.75 - 35.25 (N)' write(*,*)' Acceptable Longitude Range: 78.25 - 83.50 (W)' write(*,*)' Example: 33.231 81.2348' read(*,*)slat,slon slon=-abs(slon) write(*,*)' ' c c Check that lat lon are within appropriate bounds c if(slat.lt.31.75.or.slat.gt.35.25.or.slon.lt.-83.5.or.slon.gt.-78.25)then write(*,*)'Latitude and/or Longitude value(s) are out of range!' write(*,*)'Press Enter to terminate program' read(*,*) stop endif c c interpolate for the uniform hazard spectrum c write(*,*)' ' write(*,*)'PROBABILITY OF EXCEEDENCE' write(*,*)'There are 4 choices for probability of exceedance' write(*,*)' ' write(*,*)' 1: 2% probability of exceedance for 50 year exposure' write(*,*)' 2: 5% probability of exceedance for 50 year exposure' write(*,*)' 3: 7% probability of exceedance for 50 year exposure' write(*,*)' 4: 10%probability of exceedance for 50 year exposure' 5998 write(*,*)' ' write(*,*)'Enter 1, 2, 3, or 4, according to choice:' read(*,*)iprob if(iprob.ne.1.and.iprob.ne.2.and.iprob.ne.3.and.iprob.ne.4) goto 5998 write(*,*)' ' write(*,*)' ' write(*,*)'SITE CONDITIONS' write(*,*)'Hazard results are available for two conditions' write(*,*)' ' write(*,*)' 1: Geologically Realistic' write(*,*)' 2: Hard-rock Basement outcrop' write(*,*)' ' 5997 write(*,*)'Enter 1 or 2, according to choice:' read(*,*)isoil if(isoil.ne.1.and.isoil.ne.2)goto 5997 write(*,*)' ' c if(iprob.eq.1.and.isoil.eq.1)then filein='scenario-data/scenario-s2' !2percent_soil.dat' write(*,*)' ' write(*,*)' 2% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(*,*)' FOR GEOLOGICALLY REALISTIC SITE CONDITION' write(3,*)' 2% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(3,*)' FOR GEOLOGICALLY REALISTIC SITE CONDITION' write(3,*)' ' endif c if(iprob.eq.2.and.isoil.eq.1)then filein='scenario-data/scenario-s5' !5percent_soil.dat' write(*,*)' ' write(*,*)' 5% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(*,*)' FOR GEOLOGICALLY REALISTIC SITE CONDITION' write(3,*)' 5% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(3,*)' FOR GEOLOGICALLY REALISTIC SITE CONDITION' write(3,*)' ' endif c if(iprob.eq.3.and.isoil.eq.1)then filein='scenario-data/scenario-s7' !7percent_soil.dat' write(*,*)' ' write(*,*)' 7% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(*,*)' FOR GEOLOGICALLY REALISTIC SITE CONDITION' write(3,*)' 7% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(3,*)' FOR GEOLOGICALLY REALISTIC SITE CONDITION' write(3,*)' ' endif c if(iprob.eq.4.and.isoil.eq.1)then filein='scenario-data/scenario-s10' !10percent_soil.dat' write(*,*) write(*,*)' 10% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(*,*)' FOR GEOLOGICALLY REALISTIC SITE CONDITION' write(3,*)' 10% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(3,*)' FOR GEOLOGICALLY REALISTIC SITE CONDITION' write(3,*)' ' endif c c if(iprob.eq.1.and.isoil.eq.2)then filein='scenario-data/scenario-r2' !2percent_rock.dat' write(*,*)' ' write(*,*)' 2% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(*,*)' FOR BASEMENT OUTCROP' write(3,*)' 2% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(3,*)' FOR BASEMENT OUTCROP' write(3,*)' ' endif c c if(iprob.eq.2.and.isoil.eq.2)then filein='scenario-data/scenario-r5' !5percent_rock.dat' write(*,*)' ' write(*,*)' 5% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(*,*)' FOR BASEMENT OUTCROP' write(3,*)' 5% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(3,*)' FOR BASEMENT OUTCROP' write(3,*)' ' endif c if(iprob.eq.3.and.isoil.eq.2)then filein='scenario-data/scenario-r7' !7percent_rock.dat' write(*,*)' ' write(*,*)' 7% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(*,*)' FOR BASEMENT OUTCROP' write(3,*)' 7% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(3,*)' FOR BASEMENT OUTCROP' write(3,*)' ' endif c if(iprob.eq.4.and.isoil.eq.2)then filein='scenario-data/scenario-r10' !10percent_rock.dat' write(*,*)' ' write(*,*)' 10% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(*,*)' FOR BASEMENT OUTCROP' write(3,*)' 10% PROBABILITY OF EXCEEDANCE (For 50 year Exposure)' write(3,*)' FOR BASEMENT OUTCROP' write(3,*)' ' endif c open(unit=1,name=filein,type='old',err=9001) c call interpolate(slat,slon,filein,uhs) close(unit=1) ! make sure interpolate reads from filein above c call interpdepth(slat,slon,sedth) close(unit=1) ! make sure interpdepth reads from scenario-sites write(*,'(a33,f7.1)')'Thickness of sediments, meters: ',sedth write(3,'(a33,f7.1)')'Thickness of sediments, meters: ',sedth c*************************************************************************** c Modification 2/18/2006 to allow for Langston's result. c The orginal hazard calculations for isoil=1 (geological realistic) are c based on the assumption that Q=32 in the first layer. The following code changes allow the user to c still select that value for designing the uniform hazard spectrum, c but the default is Q=100 in the first layer. The result is realized by c c 1) first removing the site response effect calculated for Q=32, (only done for isoil=1) c 2) applying the site response for Q=100, unless the default is changed: then the c original response is used. c c hazard results for rock outcrop conditions are unchanged. The modified code c uses the same subroutines for calculating the site response used in workingversion.f c in the scdot /martin_giles/hazard/scdot directory on miles. c if(isoil.eq.1)then write(*,*)' The default value for Q in the Coastal Plain sediments is: Q=100' write(*,*)'Do you want to change the default value to Q=32 [y/n]?' read(*,'(a1)')yesno4 if(yesno4.eq.'Y'.or.yesno4.eq.'y')then qsed=32. else qsed=100. endif nl=2 ! number of layers sv(1)=0.7 ! velocity of first layer, km/s sv(2)=2.5 ! velocity of second layer km/s sh(1)=sedth/1000. ! first layer thickness, km sh(2)=0.25 vhalf=3.5 ! half space velocity km/s sbeta(2)=1./(2.*600.) ! Q in the second layer c first, remove the site response effect assumed for the hazard calculations, for isoil.eq.1 sbeta(1)=1./(2.*32.) do 333 i=1,8 if(osfreq(i).ne.-1.0)call layerresp(osfreq(i),nl,sv,sbeta,sh,vhalf,sitres) if(osfreq(i).eq.-1.0)call pgaratio1(sh(1),sitres) uhs(i)=uhs(i)/sitres 333 continue c c now apply the site response for the correct value of Q c sbeta(1)=1./(2.*qsed) do 334 i=1,8 if(osfreq(i).ne.-1.0) &call layerresp(osfreq(i),nl,sv,sbeta,sh,vhalf,sitres) if(osfreq(i).eq.-1.0.and.qsed.eq.32.)call pgaratio1(sh(1),sitres) if(osfreq(i).eq.-1.0.and.qsed.eq.100.)call pgaratio2(sh(1),sitres) uhs(i)=uhs(i)*sitres 334 continue endif c******************************************************************* write(*,*)' PSA and PGA as Percentage of g' write(*,103) write(3,*) write(3,*)' PSA and PGA as Percentage of g' write(3,103) 103 format(t4,'0.5Hz',t14,'1.0Hz',t24,'2.0Hz',t34,'3.3Hz',t44,'5Hz',t54,'6.7Hz',t64,'13Hz',t74,'PGA') write(*,102)(uhs(i),i=1,8) write(*,*) write(3,102)(uhs(i),i=1,8) write(3,*) 102 format(8(f9.5,1x)) if(iprob.eq.1)then !only have deagg eq for 2% and 10% write(*,*)' ' write(3,*)' ' write(*,'(a62)')'Interpolated results from USGS Deaggregation 2002' write(3,'(a62)')'Interpolated results from USGS Deaggregation 2002' write(*,1001)'Freq.','R(mean) km','mag(mean)','eps0(mean)','R(modal) km','mag(modal)','eps0(modal)' write(3,1001)'Freq.','R(mean) km','mag(mean)','eps0(mean)','R(modal) km','mag(modal)','eps0(modal)' filein2='scenario-data/scenario-pga2pe' open(unit=4,name=filein2,type='old',err=9001) call interphazard(slat,slon,filein2,hazeq) write(*,1002)' PGA ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) write(3,1002)' PGA ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) close(unit=4) filein2='scenario-data/scenario-5hz2pe' open(unit=4,name=filein2,type='old',err=9001) call interphazard(slat,slon,filein2,hazeq) write(*,1002)'5 Hz ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) write(3,1002)'5 Hz ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) close(unit=4) filein2='scenario-data/scenario-1hz2pe' open(unit=4,name=filein2,type='old',err=9001) call interphazard(slat,slon,filein2,hazeq) write(*,1002)'1 Hz ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) write(3,1002)'1 Hz ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) close(unit=4) 1001 format(a6,a11,1x,a11,1x,a12,1x,a11,1x,a11,1x,a12) 1002 format(a6,f11.1,1x,f11.2,1x,f12.2,1x,f11.1,1x,f11.2,1x,f12.2) endif if(iprob.eq.4)then !only have deagg eq for 2% and 10% write(*,*)' ' write(3,*)' ' write(*,'(a62)')'Interpolated results from USGS Deaggregation 2002' write(3,'(a62)')'Interpolated results from USGS Deaggregation 2002' write(*,1003)'Freq.','R(mean) km','mag(mean)','eps0(mean)','R(modal) km','mag(modal)','eps0(modal)' write(3,1003)'Freq.','R(mean) km','mag(mean)','eps0(mean)','R(modal) km','mag(modal)','eps0(modal)' filein2='scenario-data/scenario-pga10pe' open(unit=4,name=filein2,type='old',err=9001) call interphazard(slat,slon,filein2,hazeq) write(*,1004)' PGA ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) write(3,1004)' PGA ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) close(unit=4) filein2='scenario-data/scenario-5hz10pe' open(unit=4,name=filein2,type='old',err=9001) call interphazard(slat,slon,filein2,hazeq) write(*,1004)'5 Hz ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) write(3,1004)'5 Hz ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) close(unit=4) filein2='scenario-data/scenario-1hz10pe' open(unit=4,name=filein2,type='old',err=9001) call interphazard(slat,slon,filein2,hazeq) write(*,1004)'1 Hz ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) write(3,1004)'1 Hz ',hazeq(1),hazeq(2),hazeq(3),hazeq(4),hazeq(5),hazeq(6) close(unit=4) 1003 format(a6,a11,1x,a11,1x,a12,1x,a11,1x,a11,1x,a12) 1004 format(a6,f11.1,1x,f11.2,1x,f12.2,1x,f11.1,1x,f11.2,1x,f12.2) write(3,*)' ' endif c c c above gives the interpolated values of PSA uniform hazard spectrum, c from 0.5 to 13 hz and pga. c option exists to scale to a single psa or pga value.. this c require specifiying thickness of sediments c c have to interpolate to determine the depth thickness of c sediments at the site (not really require if we match the entire c uhs, but it is important for case where we match only one frequency c or pga.,so we go ahead and do, regardless. c c sedth is the interpolated thickness of the sediments. If the c interpolated value is less than 1m, it is set to zero,below. c c the interpolated uniform hazard psa spectrum, plus pga c units are in percent g c c write(*,*)' ' write(*,*)' ' write(*,*)'Generate a time series [y/n]?' read(*,*)yesno1 if(yesno1.eq.'Y'.or.yesno1.eq.'y')then c print to summary file that time series will be generated write(3,*)' ' write(3,*)' ' write(3,*)'TIME SERIES GENERATED BY THIS RUN' write(3,*)' ' c c Handle the case for the basement outcrop. This is a mod from the original version: c here the weathered rock layer is NOT included in the time series. The unscaled version c for isoil=1 will be the true stochastic simulation for rock outcrop conditions. The scaled c version will match the uniform hazard spectrum for hard rock outcrop according to the c attenuation models in the hazard calculations if(isoil.eq.2)then ! basement outcrop conditions sedth=0. else c c First, decide if user wants to change the default thickness of sediments. c this option is important, because he may be wanting to get a rock c time series simulation: e.g. he may have read in the hazard results for c basement outcrop, or he may simply be wanting to make a time series regardless c of hazard. write(*,*)' ' write(*,*)' ' write(*,*)'!! IMPORTANT DECISION FOR TIME SERIES GENERATION:' write(*,*)' Do you want to change the default value of' write(*,*)' sediment thickness (given above) to some other' write(*,*)' value for the time series simulations?' write(*,*)' This won`t make much difference if you choose' write(*,*)' to scale the time series to match the entire' write(*,*)' uniform hazard spectrum (below); but otherwise,' write(*,*)' the choice can have a significant effect on the' write(*,*)' time series.' write(*,*)' ' write(*,*)' Change the default sediment thickness? [y/n]' yesno3='N' read(*,*)yesno3 if(yesno3.eq.'Y'.or.yesno3.eq.'y')then write(*,*)' ' write(*,*)'Enter choice for sediment thickness, in meters' read(*,*)sedth endif c endif c c c c write(*,*)' ' write(*,*)'Do you want to match the entire uniform' write(*,*)'hazard spectrum [y/n]?' read(*,*)yesno2 c if(yesno2.eq.'Y'.or.yesno2.eq.'y')then c write(3,*)'The time series defined below is scaled using a' write(3,*)' phase-invariant method to match the entire' write(3,*)' uniform hazard PSA spectrum given above.' write(3,*)' ' c c define the target spectrum, on the basis of the uniform hazard c spectrum interpolated above. Note that the high frequency part c of the target spectrum is defined using pga value, and assigning it c a frequency of 40 Hz. c do 4491 i=1,8 target(i+1)=uhs(i) 4491 continue target(1)=target(2)*(0.1/0.5)**2 ! target at 0.1 Hz c based on value at 0.5 Hz, assuming c w**-2 decay at lower frequencies target(10)=target(9) ! pga ntarg=10 rtarfrq(1)=.1 rtarfrq(2)=0.5 rtarfrq(3)=1. rtarfrq(4)=2. rtarfrq(5)=3.33 rtarfrq(6)=5. rtarfrq(7)=6.7 rtarfrq(8)=13. rtarfrq(9)=40. rtarfrq(10)=50. c c write out the target psa spectrum in cgs units c filename=dirname(1:dnlen)//'/targetpsa.dat' open(unit=2,name=filename,type='new',err=9000) do 3391 i=1,10 write(2,*)rtarfrq(i),target(i)*981./100. 3391 continue close(unit=2) c c do the same thing in term of period in seconds c filename=dirname(1:dnlen)//'/targetpsa_t.dat' open(unit=2,name=filename,type='new',err=9000) ix=11 do 5391 i=1,10 ix=ix-1 write(2,*)1./rtarfrq(ix),target(ix)*981./100. 5391 continue close(unit=2) endif ! re line 72 c ITYPE=0 write(*,*) n=4096 ! hardwired. C DO 4017 I=1,12 IPOW2(I)=NINT(2.0**I) ! MAXIMUM NUMBER OF POINTS ALLOWED IS 4096 4017 CONTINUE C DO 4018 I=1,11 IF(N.GT.IPOW2(I))NN=IPOW2(I+1) 4018 CONTINUE DO 4019 I=1,12 IF(N.EQ.IPOW2(I))NN=IPOW2(I) 4019 CONTINUE N=NN ! ROUND (UP) TO THE NEAREST POWER OF 2, LE 16384. write(*,*) WRITE(*,1033)N WRITE(3,1033)N 1033 FORMAT(' Number of points used in time series is: ',I5) C write(*,*) write(*,*)'The recommended time step is 0.01 sec.' write(*,*)'Enter time step: ' read(*,*)dt write(*,*) write(3,'(a18,1x,f4.2)')' Time step used:',dt write(*,*)'The simulation requires a random number seed.' write(*,*)'Enter seed: (any integer)' read(*,*)idum write(*,*) if(idum.gt.0)idum=-idum if(idum.eq.0)idum=-1 ISEED=IDUM write(3,*)'Seed used: ',ISEED write(3,*)' ' c c option to use the atkinson and boore 1995 source model c iab95=0 c write(*,*)' ' write(*,*)'SPECIFY THE PARAMETERS OF THE SCENARIO EARTHQUAKE.' if(iprob.eq.1.or.iprob.eq.4)then write(*,*)' ' write(*,*)'The choice of magnitude and distance will affect the ' write(*,*)' duration of the synthetic event generated by this' write(*,*)' program. The recommended event for most cases is' write(*,*)' the 1 Hz modal event from the USGS deaggregation' write(*,*)' web site. (See above table of USGS deagg results.)' write(*,*)' NOTE: The 1 Hz modal event is not the appropriate' write(*,*)' choice for every application! You may need' write(*,*)' to choose a different event depending on' write(*,*)' your application.' write(*,*)' ' write(*,1005)' The 1 Hz modal magnitude for the chosen site is:',hazeq(5) write(*,1006)' The 1 Hz Modal epicentral distance (km) is:',hazeq(4) 1005 format(a51,f6.2) 1006 format(a51,f6.1) write(*,*)' ' write(*,*)'From the above information or from a deaggregation' write(*,*)' plot, enter the magnitude and epicentral distance' write(*,*)' for the chosen site.' else write(*,*)' ' write(*,*)'Enter the magnitude and epicentral distance from a' write(*,*)' deaggregation plot.' write(*,*)' ' end if write(*,*)'Moment magnitude:' read(*,*)catmag rmomag=catmag write(3,*)' ' write(3,*)'SCENARIO EVENT' write(3,*)' ' write(3,'(a26,f6.2)')' Moment magnitude:',rmomag rmo=(rmomag+10.7)*1.5 rmo=10.**rmo ! moment in dyne-cm write(*,*)' ' write(*,*)'Epicentral distance:' read(*,*)epdis write(3,'(a26,f6.1)')' Epicentral distance (km): ',epdis write(3,*)' ' c if(isoil.eq.1)then ! for geologically realistic conditions write(3,*)' Motions are based on either a 1 or 2 layer geological structure' write(3,*)' overlying a basement half space.' write(3,*)' The first layer represents a sedimentary layer with' write(3,'(a37,f5.1,a17)')' Vs=700 m/s, density 2.0 gm/cc and Q=',qsed,' The thickness of' write(3,*)' the first layer is by default that assumed in the hazard' write(3,*)' calculations. This can be changed by the user.' write(3,*)' The second layer represents weathered rock (Vs=2.5 km/sec).' write(3,'(a61,f7.2,a2)')' Sediment layer thickness used for time series simulations: ',sedth,'m' write(3,*)' ' write(3,*)' ' write(3,*)' thickness(meters) Vs(m/sec) density(gm/cc)',' Q' write(3,*)' -------------------------------------------------------------' if(sedth.ne.0.0)then write(3,*)' ' write(3,*)' ' write(3,'(t6,f7.2,t27,i3.3,t46,f3.1,9x,f5.1)')sedth,700,2.0,qsed write(3,*)' ' write(3,*)' ' write(3,*)' -------------------------------------------------------------' endif write(3,*)' ' write(3,*)' 250 2500 2.5 600' write(3,*)' ' write(3,*)' -------------------------------------------------------------' write(3,*)' ' write(3,*)' INF 3500 2.6 680f^0.36' write(3,*)' ' write(3,*)' ' write(3,*)' ' endif c if(isoil.eq.2)then ! for basement outcrop conditions write(3,*)' Motions are for basement outcrop conditions ' write(3,*)' ' write(3,*)' thickness(meters) Vs(m/sec) density(gm/cc)', &' Q' write(3,*)' -------------------------------------------------------------' write(3,*)' ' write(3,*)' INF 3500 2.6 680f^0.36' write(3,*)' ' write(3,*)' ' write(3,*)' ' endif c c c BETA=3.6 RHO=2.6 SD=100. ! use old w**2 model c c c IKAPPA=0 WM=40. WM=2.*3.141592654*WM C IAKI=0 C q=680. rqe=0.36 c DEPTH=10. DIST=SQRT(EPDIS**2+DEPTH**2) c RSTAR=100. c FC=(SD/RMO)**(1./3.) FC=FC*BETA*4.9E06 ! CORNER FREQUENCY (HZ) (brune) C IAMP=0 C C GET THEM RANDOM NUMBERS C DO 114 I=1,N RR(I)=GASDEV(IDUM) 114 CONTINUE C C ADJUST VARIANCE TO BE EQUAL TO BANDWIDTH: C NOTE: JUST WHAT THIS BANDWIDTH SHOULD BE IS A MYSTERY: C CURRENTLY USE 100 HZ IF DT IS .O1 SEC, FOR EXAMPLE C DO 1114 I=1,N RR(I)=RR(I)/SQRT(DT) ! CAUSES AMP. SPECTRUM OF RR TO BE 1, ON AVERAGE. 1114 CONTINUE C C CHECK TO SEE THAT MEAN IS ZERO AND THAT THE VARIANCE IS 1/DT. C CALL STAT2(RR,N) C C C cfreq=fc ! brune (hanks and Mcguire) CALL WINDOW(DT,N,WIN,cfreq,EPDIS) C C CHECK THAT WINDOW IS NORMALIZED TO HAVE A SQUARED INTEGRAL EQUAL TO 1 C SUM=0. DO 5555 I=1,N if(win(i).lt.1.e-10)win(i)=0.0 ! get rid of ieee underflow/inexact. SUM=SUM+((WIN(I)**2)*DT) 5555 CONTINUE C DO 115 I=1,N RR(I)=RR(I)*WIN(I) 115 CONTINUE C C C NOTE: AMPLITUDE SPECTRUM OF RR HERE HAS AN AVERAGE VALUE OF 0.904, C DETERMINED FROM 19 TRIAL RUNS. I DONT TRY TO CORRECT FOR THIS C HERE (BUT SEE BELOW). THE REASON FOR THIS IS NOT UNDERSTOOD. C C TRANSFORM TO FREQUENCY DOMAIN C J=0 DO 1155 I=1,N*2,2 J=J+1 DATA(I)=RR(J) DATA(I+1)=0.0 1155 CONTINUE CALL FOUR1(DATA,N,1) ! FORWARD FFT DO 1156 I=1,N*2 1156 DATA(I)=DATA(I)*DT ! FOR CONSISTANCY WITH JAS. C C NOTE THAT THE FFT USES ONLY REAL ARITHMETIC AND THEREFORE NEEDS C TWICE AS MANY ARRAY ELEMENTS AS REAL TIME SAMPLES. THE REAL PART OF C THE FT IS CONTAINED IN THE ODD ELEMENTS OF THE DATA ARRAY. ALSO, THE C FFT DOES NOT MULTIPLY BY DT ON THE FORWARD TRANSFOMATION, AND DOES NOT C DIVIDE BY DT*N ON THE INVERSE TRANSFORM. N MUST BE A POWER OF 2. C C DETERMINE THE FOURIER FREQUENCIES CONSISTENT WITH THE TRANSFORM C DO 116 I=1,N/2+1 ! FOR 0 UP TO AND INCLUDING THE NYQUIST FREQ. FREQ(I)=FLOAT(I-1)/(DT*FLOAT(N)) ! = (I-1)*DF 116 CONTINUE J=N/2+1 ! INDEX OF THE NYQUIST FREQ. DO 1162 I=N/2+2,N ! FOR -NYQUIST+DF TO -DF HERTZ. J=J-1 FREQ(I)=-FREQ(J) ! NOTE THAT THE FREQS ARE SYMETRIC ABOUT NYQUIST. 1162 CONTINUE C C C NOW MULTIPLY BY THE MODEL EARTHQUAKE SOURCE AND PATH, ETC C WC=2.*3.141592654*FC r=0.55 ! default for radiaion pattern FS=2. ! FREE SURFACE AMPLIFICATION PRTITN=.707106781 ! SQUARE ROOT OF 2 (BOORE'S CHOICE) BETA=BETA*100000. ! CONVERT OF CM/SEC C=R*FS*PRTITN/(4.*3.141592654*RHO*(BETA**3.)) C C DETERMINE NUMBER OF POSITIVE FREQUENCIES C NPF=N/2+1 ! THE NUMBER OF POSITIVE FREQEUNCIES, COUNTING 0.0 HZ C C C RENORMALIZE THE AMPLITUDE SPECTRUM OF THE WINDOWED RANDOM SEQUENCE C TO ENSURE THAT IT'S AVERAGE VALUE IS 1.0. C SUM=0.0 DO 234 I=1,N*2,2 SUM=SUM+SQRT(DATA(I)**2+DATA(I+1)**2) 234 CONTINUE SUM=SUM/FLOAT(N) DO 235 I=1,N*2,2 DATA(I)=DATA(I)/SUM DATA(I+1)=DATA(I+1)/SUM 235 CONTINUE C C THE AMPLITUDE SPECTRUM OF THE WINDOWED, RANDOM SEQUENCE C SHOULD NOW HAVE AN AVERAGE MODULUS VALUE OF 1.0 C C NOW EVALUATE MODEL SPECTRUM AT THE ABOVE FREQUECIES C DO 44 I=1,NPF ! ONLY POSITIVE FREQUENCIES FIRST W=2.*3.141592654*FREQ(I) S=(W/WC)**2 S=S+1. S=1./S ! DISPLACEMENT AMPLITUDE SPECTRUM: NOT ACCELERATION C SPECTRUM AS PER BOORE (1983). THIS NECESSARY IN ORDER C TO USE THE HILBERT TRANSFORM BELOW. C c c P=(W/WM)**8. ! BOORE'S CHOICE FOR HIGH FREQUENCY FALL-OFF P=P+1. P=1./SQRT(P) C IF(DIST.LT.RSTAR)A(I)=C*RMO*S*P/(DIST*100000.) IF(DIST.GE.RSTAR)A(I)=C*RMO*S*P/SQRT(DIST*100000.*RSTAR*100000.) ! SEE BOORE & ATKINSON c c IF(FREQ(I).NE.0.)THEN QF=Q*(FREQ(I)**RQE) A(I)=A(I)*EXP(-W*DIST*100000./(2.*QF*BETA)) ENDIF IF(FREQ(I).EQ.0.)A(I)=A(I) ! SIMPLY NO ANELASTIC ATTENUATION C 44 CONTINUE C C NOW HANDLE NEGATIVE FREQUENCIES C J=NPF DO 45 I=NPF+1,N ! REALLY N/2+2 TO N J=J-1 A(I)=A(J) ! JUST MAKE IT SYMETRICAL ABOUT THE NYQUIST. 45 CONTINUE C C WE'VE DEFINED THE MODULUS OF A MODEL DISPLACEMENT SPECTRUM ABOVE, IN ARRAY A. C NOW WE NEED A COMPLETE CAUSAL FUNCTION WITH THE SAME MODULUS. C USE HILBERT TRANSFORM OF LOG MODULUS TO GET A MINIMUM PHASE. C C C GET NAT LOG OF A C DO 451 I=1,N A(I)=ALOG(A(I)) ! HAVE TO HAVE NON-ZERO AMPLITUDES: C THEREFORE THE USE OF DISPLACEMENT RATHER THAN C ACCELERATION. 451 CONTINUE C C GET MINIMUM PHASE C CALL HILAPHA(A,RR,N) ! RR WILL HOLD THE PHASE SPECTRUM C C MAKE THE CAUSAL FUNCTION FOR THE "MODEL" OR "TARGET" C DO 452 I=1,N A(I)=EXP(A(I)) ! RECOVER THE MODULUS 4455 WIN(I)=A(I)*COS(RR(I)) ! REAL PART OF TARGET SPECTRUM RR(I)=A(I)*SIN(RR(I)) ! IMAGINARY PART OF TARGET SPECTRUM 452 CONTINUE C C NOW CONVERT THE MODEL FROM CAUSAL DISPLACMENT TO C CAUSAL ACCELERATION C DO 453 I=1,N ! MULTIPLY BY (-IW)**2 W=2.*3.141592654*FREQ(I) WIN(I)=-WIN(I)*W*W RR(I)=-RR(I)*W*W 453 CONTINUE C C MULTIPLY FT OF NORMALIZED TIME SERIES BY MODEL ACCELERATION SPECTRUM C J=0 DO 100 I=1,N*2,2 J=J+1 SOIL=CMPLX(DATA(I),DATA(I+1)) ROCK=CMPLX(WIN(J),RR(J)) SOIL=SOIL*ROCK DATA(I)=REAL(SOIL) DATA(I+1)=AIMAG(SOIL) 100 CONTINUE C IF(ITYPE.NE.0)THEN ! MAKE A REAL VELOCITY OR DISPLACEMENT RECORD C (CORRECT PHASE) C W/R TO THE ACCELERATION RECORD. THIS IS NEW C FEATURE WITH THIS VERSION (MARCH, 1994). J=0 DO 101 I=1,N*2,2 J=J+1 IF(I.NE.1)THEN ROCK=CMPLX(0.0,-2.0*3.141592654*FREQ(J)) IF(ITYPE.EQ.3)ROCK=ROCK*ROCK SOIL=CMPLX(DATA(I),DATA(I+1)) SOIL=SOIL/ROCK ! SIMPLY DIVIDING BY -IW OR (-IW)**2 HERE. DATA(I)=REAL(SOIL) DATA(I+1)=AIMAG(SOIL) ENDIF IF(I.EQ.1)THEN ! ZERO FREQUENCY DATA(I)=0.0 DATA(I+1)=0.0 ENDIF 101 CONTINUE ENDIF c c get the proper site response C if(isoil.eq.1)then ! for "geologically reasonable" condtions c note that in case where sedth=0., the second layer c in the model (the weathered rock layer, is present. c if isoil=2, "basement outcrop", no site response is c included call newshakesub(qsed,sedth,freq,asngl,bsngl) icc=0 DO 742 I=1,N*2,2 icc=icc+1 SOIL=CMPLX(asngl(icc),bsngl(icc)) ROCK=CMPLX(DATA(I),DATA(I+1)) RESULT=SOIL*ROCK DATA(I)=REAL(RESULT) DATA(I+1)=AIMAG(RESULT) 742 CONTINUE endif C C C INVERT BACK TO TIME DOMAIN C CALL FOUR1(DATA,N,-1) J=0 DO 4193 I=1,N*2,2 J=J+1 RR(J)=DATA(I)/(DT*FLOAT(N)) 4193 CONTINUE C C WRITE RESULTS C write(*,*)' ' write(*,*)'Enter start time of seismogram (sec):' read(*,*)TSTART filename=dirname(1:dnlen)//'/unscaledtime.dat' OPEN(UNIT=2,NAME=filename,TYPE='new',err=9000) DO 3033 I=1,N T=((I-1)*DT)+TSTART tt(i)=t WRITE(2,*)T,RR(I) 3033 CONTINUE CLOSE (UNIT=2) write(*,*)' ' write(*,*)' ' write(*,*)'The unscaled time series is contained' write(*,*)' in the file "unscaledtime.dat". The' write(*,*)' units are cm/sec^2 versus seconds.' write(*,*)' ' write(3,*)'The unscaled time series is contained' write(3,*)' in the file "unscaledtime.dat". The' write(3,*)' units are cm/sec^2 versus seconds.' write(3,*)' ' c c write out the psa spectrum of the unscaled time series at 99 c frequencies between 0.1 and 50 hz c filename=dirname(1:dnlen)//'/unscaledpsa.dat' open(unit=2,name=filename,type='new',err=9000) df=(alog(50.)-alog(.1))/99. damp=.05 do 4920 i=1,99 ff(i)=alog(0.1)+(i-1)*df ff(i)=exp(ff(i)) call rspectsub(tt,rr,n,ff(i),.05,sdd) sdtemp(i)=sdd write(2,*)ff(i),sdd*(2.*3.141592654*ff(i))**2 4920 continue close(unit=2) c c do the same thing, with the x axis being period instead of hertz c ix=100 filename=dirname(1:dnlen)//'/unscaledpsa_t.dat' open(unit=2,name=filename,type='new',err=9000) do 6921 i=1,99 ix=ix-1 write(2,*)1./ff(ix),sdtemp(ix)*(2.*3.141592654*ff(ix))**2 6921 continue close(unit=2) write(*,*)'The PSA spectrum for the unscaled time series is' write(*,*)' contained in the file "unscaledpsa.dat". The' write(*,*)' units are cm/sec^2 versus frequency in Hz.' write(*,*)' ' write(3,*)'The PSA spectrum for the unscaled time series is' write(3,*)' contained in the file "unscaledpsa.dat". The' write(3,*)' units are cm/sec^2 versus frequency in Hz.' write(3,*)' ' write(*,*)'The same spectrum with x axis in terms of period' write(*,*)' in seconds is in file "unscaledpsa_t.dat"' write(*,*)' ' write(3,*)'The same spectrum with x axis in terms of period' write(3,*)' in seconds is in file "unscaledpsa_t.dat"' write(3,*)' ' c C c do the spectral scaling c if(yesno2.eq.'Y'.or.yesno2.eq.'y')then c call rat(rr,n,dt,target,ntarg,rtarfrq) ! match entire uhs goto 5567 c else write(*,*)' ' write(*,*)'You can either scale to match an arbitrary PGA value' write(*,*)'or you can scale to match a PSA response spectrum ordinate.' write(*,*)' ' write(*,*)'Enter 1 if you want to scale to PGA: 2 otherwise' read(*,*)iscale c if(iscale.eq.1)then ! scale to PGA c write(3,*)' ' write(3,*)'The time series for this scenario is scaled to' write(3,*)' match a user specified peak acceleration value of' xpga=0. do 5689 izz=1,n if(abs(rr(izz)).gt.xpga)then xpga=abs(rr(izz)) xxpga=rr(izz) endif 5689 continue xpga=xpga*100./981. xxpga=xxpga*100./981. write(*,*)' ' write(*,'(a39,f9.4,a10)')'The PGA of the unscaled time series is:',xxpga,' percent g' write(*,*)' ' write(*,*)'Enter value of PGA for the scaled time series,' write(*,*)'in units of PERCENT g' read(*,*)xscale write(3,'(a6,f8.3,a3)')' PGA =',xscale,' %g' write(3,*)' ' if(xscale.lt.0.0.and.xxpga.lt.0.0)sign=1.0 if(xscale.lt.0.0.and.xxpga.ge.0.0)sign=-1.0 if(xscale.ge.0.0.and.xxpga.ge.0.0)sign=1.0 if(xscale.ge.0.0.and.xxpga.lt.0.0)sign=-1.0 xscale=sign*abs(xscale)/xpga c else ! scale to some PSA value c write(*,*)' ' write(*,*)'Enter the oscillator frequency (in Hertz) at' write(*,*)'which you want to scale the scenario time series.' read(*,*)xfreq write(3,*)' ' write(3,*)'The time series for this scenario is scaled' write(3,*)' independently of frequency to match a user' write(3,*)' specified PSA value of' call rspectsub(tt,rr,n,xfreq,0.05,xsd) write(*,*)' ' write(*,'(a12,f6.3,a24)')'The PSA at ',xfreq,' Hz of the unscaled time' write(*,'(a12,f8.4,a10)')'series is: ',xsd*(2.*3.14159*xfreq)**2*100./981.,' percent g' write(*,*)' ' write(*,*)' ' write(*,*)'Enter the value of PSA for the scaled time series,' write(*,*)'in units of PERCENT g' read(*,*)xscale write(3,'(t2,f8.4,a7,1x,f7.3,a3)')xscale,' %g at ',xfreq,' Hz' write(3,*)' ' xscale=xscale/(xsd*(2.*3.141592654*xfreq)**2*100./981.) endif c do 1122 i=1,n rr(i)=rr(i)*xscale 1122 continue c endif ! re line 459 c c 5567 continue c write out the modified time series c filename=dirname(1:dnlen)//'/scaledtime.dat' open(unit=2,name=filename,type='new',err=9000) do 3233 i=1,n t=(i-1)*dt+tstart write(2,*)t,rr(i) 3233 continue close(unit=2) write(*,*)' ' write(*,*)'The scaled time series is contained in' write(*,*)' the file "scaledtime.dat." The units' write(*,*)' are cm/sec^2.' write(*,*)' ' write(3,*)' ' write(3,*)'The scaled time series is contained in' write(3,*)' the file "scaledtime.dat." The units' write(3,*)' are cm/sec^2.' write(3,*)' ' c c write out the psa spectrum of the scaled time series c filename=dirname(1:dnlen)//'/scaledpsa.dat' open(unit=2,name=filename,type='new',err=9000) do 4921 i=1,99 call rspectsub(tt,rr,n,ff(i),.05,sdd) sdtemp(i)=sdd write(2,*)ff(i),sdd*(2.*3.141592654*ff(i))**2 4921 continue close(unit=2) write(*,*)'The PSA spectrum for the scaled time series is' write(*,*)' contained in the file "scaledpsa.dat". The units' write(*,*)' are cm/sec^2 versus frequency in Hz.' write(*,*)' ' write(3,*)'The PSA spectrum for the scaled time series is' write(3,*)' contained in the file "scaledpsa.dat". The units' write(3,*)' are cm/sec^2 versus frequency in Hz.' write(3,*)' ' c c do the same thing in terms of period c filename=dirname(1:dnlen)//'/scaledpsa_t.dat' open(unit=2,name=filename,type='new',err=9000) ix=100 do 5921 i=1,99 ix=ix-1 write(2,*)1./ff(ix),sdtemp(ix)*(2.*3.141592654*ff(ix))**2 5921 continue close(unit=2) write(*,*)'The same spectrum with x axis in terms of' write(*,*)'period in seconds is in file "scaledpsa_t.dat"' write(*,*)' ' write(3,*)'The same spectrum with x axis in terms of' write(3,*)'period in seconds is in file "scaledpsa_t.dat"' write(3,*)' ' C C CALCULATE SOME SUMMARY STUFF C C MAX AND MIN VALUES OF TIME SERIES C RMAX=-9.E20 RMIN=9.E20 DO 1021 I=1,N IF(RR(I).LE.RMIN)RMIN=RR(I) IF(RR(I).GE.RMAX)RMAX=RR(I) 1021 CONTINUE C C INTEGRAL OF SQUARED TIME SERIES C RINT=0. DO 1022 I=1,N RINT=RINT+((RR(I)**2)*DT) 1022 CONTINUE C C RMS AMPLITUDE C RMS=SQRT(RINT/(DT*FLOAT(N))) C C DURATION C RLOW=0.05*RINT RHI=0.95*RINT SUM=0. DO 721 I=1,N SUM=SUM+((RR(I)**2)*DT) IF(SUM.LE.RLOW)T1=(I-1)*DT IF(SUM.LE.RHI)T2=(I-1)*DT 721 CONTINUE DUR=T2-T1 C C WRITE(3,*)' ' WRITE(3,*)' ' WRITE(3,*)' ' WRITE(3,*)'SUMMARY OF STOCHASTIC MODELING:' WRITE(3,*)' ' IF(ITYPE.EQ.1)WRITE(3,746) 746 FORMAT(' VELOCITY RECORD') IF(ITYPE.EQ.0)WRITE(3,747) 747 FORMAT(' ACCELERATION RECORD') IF(ITYPE.EQ.3)WRITE(3,748) 748 FORMAT(' DISPLACEMENT RECORD') if(iab95.eq.1)write(3,330) 330 format(' Atkinson and Boore 1995 source model') if(iab95.ne.1)write(3,331) 331 format(' brune source model') WRITE(3,800)RMIN,RMAX 800 FORMAT(' PEAK AMPLITUDES: ',1X,E13.6,1X,E13.6) WRITE(3,801)RINT 801 FORMAT(' INTEGRAL OF SQUARED AMPLITUDES: ',E13.6) WRITE(3,802)RMS 802 FORMAT(' RMS AMPLITUDE: ',E13.6) WRITE(3,803)T1,T2,DUR 803 FORMAT( ' T1: ',E13.6,1X,'T2: ',E13.6,1X,'DURATION: ',E13.6) WRITE(3,804)N,DT,ISEED 804 FORMAT(' N: ',I5,1X,'DT: ',F5.3,1X,'SEED: ',I8) WRITE(3,805)RMO,BETA,RHO 805 FORMAT(' MOMENT: ',E13.6,1X,'VELOCITY: ',E13.6,1X,'DENSITY: ',E13.6) C IF(IKAPPA.EQ.0)THEN if(iab95.ne.1)WRITE(3,806)SD,WM if(iab95.eq.1)write(3,332)fa,fb,wm 332 format(' Fa (hz): ',e13.6,1x,' Fb (hz): ',e13.6,1x,' fmax: ',e13.6,'(rad/sec)') 806 FORMAT(' STRESS DROP: ',E13.6,1X,'FMAX (RAD/SEC): ',E13.6) ELSE if(iab95.ne.1)WRITE(3,8066)SD,RKAPPA if(iab95.eq.1)write(3,3333)fa,fb,rkappa 3333 format(' Fa (hz): ',e13.6,1x,' fb (hz): ',e13.6,1x,' Kappa: ',e13.6) 8066 FORMAT(' STRESS DROP: ',E13.6,1X,'KAPPA: ',E13.6) ENDIF C WRITE(3,807)Q,RQE 807 FORMAT(' Q: ',E13.6,1X,'EXPONENT: ',E13.6) WRITE(3,808)EPDIS,RSTAR,DEPTH 808 FORMAT(' EP. DIS.: ',E13.6,1X,'CROSSOVER: ',E13.6,1X,'FOCAL DEPTH: ',E13.6) CLOSE (UNIT=3) write(*,*)' ' write(*,*)' ' write(*,*)'SHAKE91 FILES' write(*,*)'Do you want the time series files written' write(*,*)' into SHAKE91 input format? [y/n]' yesno1='n' read(*,*)yesno1 c if(yesno1.eq.'Y'.or.yesno1.eq.'y')then filename=dirname(1:dnlen)//'/unscaledtime.dat' open(unit=1,name=filename,type='old',err=9001) do 624 i=1,n read(1,*)tt(i),shake(i) 624 continue close(unit=1) filename=dirname(1:dnlen)//'/unscaledshake.dat' call mcc2shake(n,shake,tt,dt,filename) c filename=dirname(1:dnlen)//'/scaledtime.dat' open(unit=1,name=filename,type='old',err=9001) do 625 i=1,n read(1,*)tt(i),shake(i) 625 continue close(unit=1) filename=dirname(1:dnlen)//'/scaledshake.dat' call mcc2shake(n,shake,tt,dt,filename) write(*,*) write(*,*)'The SHAKE compatible input files are' write(*,*)' "unscaledshake.dat" and "scaledshake.dat"' write(*,*) endif c endif write(*,*) write(*,*)'Analysis Complete' write(*,*) WRITE(*,*)'Summary results are contained in file "summary.txt"' write(*,*)' ' write(*,*)'Press enter to terminate program' read(*,*) continue STOP 9000 write(*,*)'I/O ERROR: Cannot open an input file.' write(*,*)'Make sure all input files are in "scenario-data".' write(*,*)'Press Enter to terminate program.' read(*,*) STOP 9001 write(*,*)'I/O ERROR: Cannot open an output file.' write(*,*)'Make sure output files do not already exist.' write(*,*)'Press Enter to terminate program.' read(*,*) STOP END C C C ******** THIS IS THE END OF THE MAIN ********* C C C FUNCTION gasdev(idum) c this is the version of gasdev that I got from my microsoft pc fortran numerical recipes c library. I had to clean up the line endings, but it is the same as the version c printed in the Press et al f77 (blue) book. It, combined with the updated version of c ran1 works on linux. The older version did not. INTEGER idum REAL gasdev CU USES ran1 INTEGER iset REAL fac,gset,rsq,v1,v2,ran1 SAVE iset,gset DATA iset/0/ if (iset.eq.0) then 1 v1=2.*ran1(idum)-1. v2=2.*ran1(idum)-1. rsq=v1**2+v2**2 if(rsq.ge.1..or.rsq.eq.0.)goto 1 fac=sqrt(-2.*log(rsq)/rsq) gset=v1*fac gasdev=v2*fac iset=1 else gasdev=gset iset=0 endif return END FUNCTION ran1(idum) c updated version of ran1, from f77 version of press et al. It c works on linux. the old version does not. INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV REAL ran1,AM,EPS,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS) INTEGER j,k,iv(NTAB),iy SAVE iv,iy DATA iv /NTAB*0/, iy /0/ if (idum.le.0.or.iy.eq.0) then idum=max(-idum,1) do 11 j=NTAB+8,1,-1 k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM j=1+iy/NDIV iy=iv(j) iv(j)=idum ran1=min(AM*iy,RNMX) return END C C SUBROUTINE WINDOW(DT,N,W,FC,DIST) DIMENSION W(16384) C C DEVELOPES THE TIME DOMAIN WINDOW USED BY BOORE, C BSSA VOL 73, NO. 6, PP. 1865-1894 C AND INCORPORATES DISTANCE DEPENDENT EFFECT ON DURATION: C SEE BOORE AND ATKINSON, BSSA V. 77 PP. 445. TD=1./FC ! SOURCE DURATION c TD=TD+(0.05*DIST) ! SEE BOORE AND ATKINSON bssa, 1987 c this is from Herrmann, and has c been replaced: atkinson and boore c if(dist.lt.10.0)dd=0.0 ! at& boore bssa 1995. if(dist.ge.10.0.and.dist.le.70.)dd=-1.6+.16*dist if(dist.gt.70.0.and.dist.le.130.)dd=11.7-.03*dist if(dist.gt.130.)dd=2.6+.04*dist td=td+dd E= 0.2 ! BOORE'S CHOICE FOR EPSILON RNU=0.05 ! BOORE'S CHOICE FOR NU B=ALOG(E)-1. B=E*B B=1.+B B=-E*ALOG(RNU)/B ! BOORE'S b TW=2.*TD C=B/(E*TW) ! BOORE'S c TWOBP1=(2.*B)+1. GAM=GAMMLN(TWOBP1) GAM=EXP(GAM) A=(2.*C)**TWOBP1 A=A/GAM A=SQRT(A) ! BOORE'S a DO 1 I=1,N T=FLOAT(I-1)*DT if(c*t.gt.70.)then ! this avoids exponential underflow. w(i)=0.0 else W(I)=A*(T**B)*EXP(-C*T) endif 1 CONTINUE RETURN END FUNCTION GAMMLN(XX) REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, * -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/ DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/ X=XX-ONE TMP=X+FPF TMP=(X+HALF)*LOG(TMP)-TMP SER=ONE DO 11 J=1,6 X=X+ONE SER=SER+COF(J)/X 11 CONTINUE GAMMLN=TMP+LOG(STP*SER) RETURN END SUBROUTINE STAT2(RR,N) DIMENSION RR(16384) SUM=0. DO 1 I=1,N SUM=SUM+RR(I) 1 CONTINUE RMEAN=SUM/FLOAT(N) SUM=0. DO 2 I=1,N SUM=SUM+(RR(I)-RMEAN)**2 2 CONTINUE VAR=SUM/FLOAT(N) c WRITE(*,100)RMEAN,VAR 100 FORMAT(' MEAN: ',E13.7,' VARIANCE: ',E13.7) RETURN END C C************************************************************************** C SUBROUTINE ROCKAMP(F,ANS) C C THIS IMPLIMENTS AN AMPLIFICATION FACTOR DERIVED BY BOORE 1986, BSSA 76 C PP 43 FOR ROCK SITES IN WESTERN NORTH AMERICA. ITS USE IS ALSO DISCUSSED C IN JOYNER AND BOORE 1988 ASCE PAPER. I FIT HIS AMPLIFICATION VALUES IN C LOG A VS LOG FREQ WITH A STRAIGHT LINE FOR IMPLIMENTATION. BOORE INTERPOLATES C BUT I DONT THINK THE DIFFERENCE IS SIGNIFICANT C IF(F.GE.0.1341.AND.F.LE.10.)AMPF=1.5066*(F**0.204) IF(F.LT.0.1341)AMPF=1. IF(F.GT.10.)AMPF=2.409894 ANS=ANS*AMPF RETURN END C************************************************************************** SUBROUTINE HILAPHA(LNMOD,PHASE,NN) COMPLEX P(16384),Z REAL LNMOD(1),PHASE(1),DATA(32768) C COMPUTES PHASE FROM LOG AMP SPECTRUM USING FREQUENCY DOMAIN REALIZATION OF C HILBERT TRANSFROM. C C LNMOD CONTAINS THE NATURAL LOG OF THE MODULUS OF THE FOURIER TRANSFORM OF C SOME TIME SERIES IN WRAP AROUND ORDER: ZERO TO NYQUIST, -NYQUIST TO ZERO. C NN IS THE SAME AS IN FOUR1: THE ACTUAL NUMBER OF TIME SERIES DATA POINTS. C C PUT THE LOG OF THE MODULUS IN THE REAL PART OF DATA ARRAY (ODD ELEMENTS) C DO 55 I=1,32684 55 DATA(I)=0. IC=0 DO 3 I=1,NN*2,2 IC=IC+1 DATA(I)=LNMOD(IC) 3 DATA(I+1)=0. C C TRANSFORM C CALL FOUR1(DATA,NN,-1) DO 4 I=1,NN*2 4 DATA(I)=DATA(I)/NN C IC=0 DO 5 I=1,NN*2,2 IC=IC+1 5 P(IC)=CMPLX(DATA(I),DATA(I+1)) Z=(0.0,1.0) DO 6 I=1,NN/2 ! FOR POSITIVE FREQS 6 P(I)=-P(I)*Z DO 7 I=NN/2+2,NN ! NEG FREQS 7 P(I)=P(I)*Z P(1)=(0.0,0.0) ! ZERO FREQ P(NN/2+1)=(0.0,0.0) ! NYQUIST C IC=0 DO 8 I=1,NN*2,2 IC=IC+1 DATA(I)=REAL(P(IC)) 8 DATA(I+1)=AIMAG(P(IC)) C CALL FOUR1(DATA,NN,1) C IC=0 DO 10 I=1,NN*2,2 IC=IC+1 10 PHASE(IC)=DATA(I) ! PHASE C RETURN END subroutine rat(sig,nsig,dt,tarspect,nfreq,rfreq) c c specially modified version of myrattle. c c Does spectral matching. Given a response spectrum, c this program will modify an input c time series signal to match it. c c The input time series must be less than 4096 points c The spectrum to be matched must have less than 1000 points, and c be in mcc format, frequency (hz), spectral amplitude. c real signal(4096),rfreq(1000),tarspect(1000),x(1000),y(1000) real y2(1000),t(4096),psv(4097),psa(4097),omega(4097) real target(4097),ratio(4097),modt(4097),modamp(4097) complex ans(4097),mod(4097) real sig(16384) do 203 i=1,nsig signal(i)=sig(i) ! this is because of the calling program array size. t(i)=(i-1)*dt 203 continue c c itype=1 damp=.05 c write(*,*)' ' write(*,*)'*****************************************************' write(*,*)'ITERATIVE SPECTRAL MATCHING' write(*,*)'Set the number of iterations to match target spectrum.' write(*,*)'5 iterations usually do a good job.' write(*,*) read(*,*)nit c c convert from percent g to cm/sec**2 c do 44 i=1,nfreq tarspect(i)=tarspect(I)*981./100. 44 continue c c c fit it with a cubic spline C DO 35 I=1,nfreq X(I)=ALOG(rfreq(I)) Y(I)=ALOG(tarspect(I)) 35 CONTINUE C diff1=(y(2)-y(1))/(x(2)-x(1)) diffn=(y(nfreq)-y(nfreq-1))/(x(nfreq)-x(nfreq-1)) call splineX(x,y,nfreq,diff1,diffn,y2) c do 57 izzq=1,nit write(*,*)'iteration: ',izzq c c get response spectrum of input time series at each fourier c frequency, in the band of the target response spectrum c c first, get the time series into the frequency domain. c call ftfwdpad(nsig,dt,signal,ans,omega,nposfreq) c c routine takes times series in array signal (which is length c narray), pads it with zeros to lenght 4096*2, and then c flips to the frequency domain. It creates a complex array ans c of length 4097 that holds the positive frequency fourier c transform of the padded time series, and also returns in c array omega, the corresponding RADIAN frequencies. It also c returns the value 4097(number of pos frequencies in padded c transform) in variable nfreq) c narray must of coarse, be less than or equal to 4096, to c avoid c wrap-around. After manipulating the answer, use c subroutine c ftinvpad to flip to the time domain, and get the time c series c of unpadded length narray. c c ic=0 do 10 i=1,nposfreq frq=omega(i)/(2.*3.141592654) if(frq.ge.rfreq(1).and.frq.le.rfreq(nfreq))then ic=ic+1 call rspectsub(t,signal,nsig,frq,damp, &sd) psv(ic)=omega(i)*sd psa(ic)=psv(ic)*omega(i) c c c determine ratio of the response spectrum of the input time c series divided by the response spectrum of the spline fit to c the target spectrum, evaluated at the fourier frequencies of c the input time series, between the maximum and minimum freq. of c the target response spectrum c val=alog(frq) call splint(x,y,y2,nfreq,val,result) C target(ic)=exp(result) if(itype.eq.1)ratio(ic)=psa(ic)/target(ic) if(itype.eq.2)ratio(ic)=psv(ic)/target(ic) endif 10 continue c c c now multiply the fourier amplitude spectrum of the input c time series by the inverse of the ratio above, in the freqency c band of the target specrum. Keep the phase the same as the c input time series. Outside the band of the target spectrum, c we simply scale according to the ratio(1) and ratio(ic). c icc=0 do 20 i=1,nposfreq frq=omega(i)/(2.*3.141592654) if(i.ne.1)then phase=atan2(aimag(ans(I)),real(ans(I))) else phase=0. endif c if(frq.lt.rfreq(1))then modamp(i)=cabs(ans(i))/ratio(1) endif c if(frq.ge.rfreq(1).and.frq.le.rfreq(nfreq))then icc=icc+1 modamp(i)=cabs(ans(i))/ratio(icc) endif c if(frq.gt.rfreq(nfreq))then modamp(i)=cabs(ans(i))/ratio(ic) endif c c rea=modamp(i)*cos(phase) rimag=modamp(i)*sin(phase) mod(i)=cmplx(rea,rimag) ! fourier transform of modified c 20 continue c call ftinvpad(nsig,dt,mod,modt) c c re-set the signal array for another iteration do 329 i=1,nsig signal(i)=modt(i) 329 continue c 57 continue cc do 492 i=1,nsig sig(i)=signal(i) 492 continue return END c c SUBROUTINE SPLINEX(X,Y,N,YP1,YPN,Y2) PARAMETER (NMAX=100) DIMENSION X(N),Y(N),Y2(N),U(NMAX) IF (YP1.GT..99E30) THEN Y2(1)=0. U(1)=0. ELSE Y2(1)=-0.5 U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) ENDIF DO 11 I=2,N-1 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2. Y2(I)=(SIG-1.)/P U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 11 CONTINUE IF (YPN.GT..99E30) THEN QN=0. UN=0. ELSE QN=0.5 UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) ENDIF Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.) DO 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE RETURN END SUBROUTINE SPLINTX(XA,YA,Y2A,N,X,Y) DIMENSION XA(N),YA(N),Y2A(N) KLO=1 KHI=N 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF(XA(K).GT.X)THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF H=XA(KHI)-XA(KLO) IF (H.EQ.0.) PAUSE 'Bad XA input.' A=(XA(KHI)-X)/H B=(X-XA(KLO))/H Y=A*YA(KLO)+B*YA(KHI)+ * ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6. RETURN END SUBROUTINE splint(xa,ya,y2a,n,x,y) INTEGER n REAL x,y,xa(n),y2a(n),ya(n) INTEGER k,khi,klo REAL a,b,h klo=1 khi=n 1 if (khi-klo.gt.1) then k=(khi+klo)/2 if(xa(k).gt.x)then khi=k else klo=k endif goto 1 endif h=xa(khi)-xa(klo) if (h.eq.0.) pause 'bad xa input in splint' a=(xa(khi)-x)/h b=(x-xa(klo))/h y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h** *2)/6. return END SUBROUTINE FOUR1(DATA,NN,ISIGN) REAL*8 WR,WI,WPR,WPI,WTEMP,THETA DIMENSION DATA(2*NN) N=2*NN J=1 DO 11 I=1,N,2 IF(J.GT.I)THEN TEMPR=DATA(J) TEMPI=DATA(J+1) DATA(J)=DATA(I) DATA(J+1)=DATA(I+1) DATA(I)=TEMPR DATA(I+1)=TEMPI ENDIF M=N/2 1 IF ((M.GE.2).AND.(J.GT.M)) THEN J=J-M M=M/2 GO TO 1 ENDIF J=J+M 11 CONTINUE MMAX=2 2 IF (N.GT.MMAX) THEN ISTEP=2*MMAX THETA=6.28318530717959D0/(ISIGN*MMAX) WPR=-2.D0*DSIN(0.5D0*THETA)**2 WPI=DSIN(THETA) WR=1.D0 WI=0.D0 DO 13 M=1,MMAX,2 DO 12 I=M,N,ISTEP J=I+MMAX TEMPR=SNGL(WR)*DATA(J)-SNGL(WI)*DATA(J+1) TEMPI=SNGL(WR)*DATA(J+1)+SNGL(WI)*DATA(J) DATA(J)=DATA(I)-TEMPR DATA(J+1)=DATA(I+1)-TEMPI DATA(I)=DATA(I)+TEMPR DATA(I+1)=DATA(I+1)+TEMPI 12 CONTINUE WTEMP=WR WR=WR*WPR-WI*WPI+WR WI=WI*WPR+WTEMP*WPI+WI 13 CONTINUE MMAX=ISTEP GO TO 2 ENDIF RETURN END c subroutine costape(n,data,percent,dt) dimension data(n) c puts a cosine taper on the end of the data array c percent is the percentage of the array that is UNTAPERED c the program is smart enough to figure out if you sent down c the percentage as a decimal fraction or as a percentage c if(percent.gt.1.1)percent=percent/100.0 c See c tmax=dt*n tstart=percent*tmax ! time of first point to begin tapering omega=(acos(-1.0)/2.)/(tmax-tstart) do 1 i=1,n if(i*dt.gt.tstart)then win=cos(omega*(i*dt-tstart)) if(win.lt.0.0)win=0.0 data(i)=data(I)*win endif 1 continue return end c subroutine decimate(xdt,x,nx,ans) c takes the data array x, filters and decimates it to c one half the number of data elements. c xdt: the time step of the x array, provided by main calling routine c x: the oversampled data array with time step xdt, provided by main c nx: the dimension of the oversampled array x, provided by main. c ans: the filtered and decimated array,must be dimensioned correctly in main c ie., must be dimensioned as nx/2 in main routine. nx must be a c power of 2. c real x(nx),dum(32768),ans(nx/2) complex comp c if(nx.gt.16384)write(*,*)'nx is too big in subroutine decimate' c do 1 i=1,32768 dum(i)=0.0 1 continue c ic=0 do 2 i=1,nx*2,2 ic=ic+1 dum(i)=x(ic) 2 continue c nn=32768/2 ! total number of padded real time series points. rnyq=1./(2.*xdt) ! the current nyquist frequency df=1./(nn*xdt) ! the current frequency increment fc=0.7*rnyq/2. ! the low pass corner frequency. It will be 70% of the c nyquist of the decimated time series. call four1(dum,nn,1) c ic=0 do 3 i=1,(nn/2+1)*2,2 ! filter the positive freqs. ic=ic+1 comp=cmplx(dum(i),dum(i+1)) freq=(ic-1)*df call butter(comp,freq,fc,1,fc,xdt,4) ! 4 pole lowpass at fc hz. dum(i)=real(comp)*xdt dum(i+1)=aimag(comp)*xdt 3 continue c c now generate the spectrum from -nyquist to 0-df Hz. c ic=((nn/2+1)*2)-1 ! index of the real part of the nyquist do 7 i=((nn/2+1)*2)+1,nn*2,2 ic=ic-2 dum(i)=dum(ic) dum(i+1)=-dum(ic+1) 7 continue c call four1(dum,nn,-1) c decimate and return c ic=0 do 8 i=1,nx*2,4 ic=ic+1 ans(ic)=dum(i)*df 8 continue if(ic*2.ne.nx)write(*,*)'decimation error' return end c C SUBROUTINE BUTTER(X,FREQ,FC,ITYPE,FCHI,DT,NPOLE) C C BUTTERWORTH FILTER WITH VARIABLE CORNER FREQUENCY. C THIS WAS DERIVED FROM A BW LOW PASS PROTOTYPE WITH A CORNER OF 0.001 HZ, C AND UNIT SAMPLE INTERVAL. IT USES PRE-WARPING OF THE FREQUENCY AXIS C AND THE BILINEAR S TO Z TRANSFORMATION. C SEE BOSE, N.K. (1985), DIGITAL FILTERS, PP 163,167 FOR BUTTERWORTH FILTERS. C C ITYPE=1 LOW PASS C ITYPE=2 HI PASS C ITYPE=3 BAND PASS C ITYPE=4 BAND REJECT C C FC IS THE CUTOFF FREQ FOR BOTH HI AND LOW PASS: ALSO LOW FREQ CUTOFF C FOR BAND PASS, LOW FREQUENCY CORNER ON BAND REJECT. C FCHI IS THE HI FREQ CUTOFF FOR BANDPASS, HI FREQUENCY CORNER ON BAND REJECT. C C NPOLE IS NUMBER OF POLES: 1,2,..10: GIVES THE SLOPE ON LOG AMP VS LOG F PLOT. C C PHASE IS CONSISTENT WITH AKI & RICHARDS FOURIER TRANSFORM CONVENTION: C IT INCREASES WITH FREQUENCY. FOR THE OTHER CONVENTION, DELETE THE C X=CONJG(H) AND REPLACE WITH X=H. C C MCC (LATEST MOD: 3/91) C COMPLEX II,X,H,ZLP,Z,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,C II=CMPLX(0.,1.) PI2=2.*ACOS(-1.0) W=PI2*FREQ*DT ! THIS NORMALIZES THE FREQUENCY SCALE TO THE PROTOTYPE ZLP=CEXP(II*W) IF(W.GT.PI2/2.0)write(*,*)'FREQUENCIES ARE ABOVE THE NYQUIST' WC=0.001*PI2 ! CUTOFF OF ORIGINAL LOW PASS FILTER C IF(ITYPE.EQ.2)THEN ! HI PASS FILTER WCHP=PI2*FC*DT A=-COS((WCHP+WC)/2.0) A=A/(COS((WCHP-WC)/2.0)) Z=-((1.0/ZLP)+A)/(1.0+A/ZLP) Z=1.0/Z ENDIF C IF(ITYPE.EQ.1)THEN ! LOW PASS WCLP=PI2*FC*DT A=SIN((WC-WCLP)/2.) A=A/SIN((WC+WCLP)/2.) Z=((1./ZLP)-A)/(1.-A/ZLP) Z=1./Z ENDIF C IF(ITYPE.EQ.3)THEN ! BAND PASS WCL=PI2*FC*DT WCH=PI2*FCHI*DT A=COS((WCH+WCL)/2.) A=A/COS((WCH-WCL)/2.) B=TAN((WCH-WCL)/2.) B=(1./B)*TAN(WC/2.) Z=1./(ZLP**2)-(2.*A*B)/((B+1.)*ZLP)+(B-1.)/(B+1.) Z=Z/((B-1.)/((B+1.)*ZLP*ZLP)-2.*A*B/((B+1.)*ZLP)+1.) Z=-1./Z ENDIF C IF(ITYPE.EQ.4)THEN ! BAND REJECT WCL=PI2*FC*DT WCH=PI2*FCHI*DT A=COS((WCH+WCL)/2.) A=A/COS((WCH-WCL)/2.) B=TAN((WCH-WCL)/2.)*TAN(WC/2.) Z=1./(ZLP**2)-(2.*A)/((1.+B)*ZLP)+(1.-B)/(1.+B) Z=Z/((1.-B)/((1.+B)*ZLP*ZLP)-(2.*A)/((1.+B)*ZLP)+1.) Z=1./Z ENDIF C IF(NPOLE.EQ.1)THEN H=3.131765E-03*(Z+1.0)/(Z-0.9937364715881357) ENDIF C IF(NPOLE.EQ.2)THEN P1=CMPLX(0.995571461339282,4.4232020327230668E-03) P2=CMPLX(0.995571461339282,-4.4232020327230668E-03) C=(2.0**2)/((1.0-P1)*(1.0-P2)) C=1.0/C H=C*((Z+1.0)**2)/((Z-P1)*(Z-P2)) ENDIF C IF(NPOLE.EQ.3)THEN P1=CMPLX(0.9968485892806065,5.4243213534056843E-03) P2=CMPLX(0.9937364715881357,0.0) P3=CMPLX(0.9968485892806065,-5.4243213534056843E-03) C=(2.0**3)/((1.0-P1)*(1.0-P2)*(1.0-P3)) C=1.0/C H=C*((Z+1.0)**3)/((Z-P1)*(Z-P2)*(Z-P3)) ENDIF C IF(NPOLE.EQ.4)THEN P1=CMPLX(0.9975816205264062,5.7909443190584446E-03) P2=CMPLX(0.9942090084129309,2.3905782210500918E-03) P3=CMPLX(0.9942090084129309,-2.3905782210500917E-03) P4=CMPLX(0.9975816205264062,-5.7909443190584445E-03) C=(2.0**4)/((1.0-P1)*(1.0-P2)*(1.0-P3)*(1.0-P4)) C=1.0/C H=C*((Z+1.0)**4)/((Z-P1)*(Z-P2)*(Z-P3)*(Z-P4)) ENDIF C IF(NPOLE.EQ.5)THEN P1=CMPLX(0.9980424632487176,5.9640455141251982E-03) P2=CMPLX(0.9949228981409190,3.6744616215257588E-03) P3=CMPLX(0.9937364712450427,0.0) P4=CMPLX(0.9949228981409190,-3.6744616215257590E-03) P5=CMPLX(0.9980424632487176,-5.9640455141251983E-03) C=(2.0**5)/((1.0-P1)*(1.0-P2)*(1.0-P3)*(1.0-P4)*(1.0-P5)) C=1.0/C H=C*((Z+1.0)**5)/((Z-P1)*(Z-P2)*(Z-P3)*(Z-P4)*(Z-P5)) ENDIF C IF(NPOLE.EQ.6)THEN P1=CMPLX(0.9983567357277001,6.0591978630835325E-03) P2=CMPLX(0.9955571458898052,4.4232022746910896E-03) P3=CMPLX(0.9939479398032704,1.6163874602061155E-03) P4=CMPLX(0.9939479398032704,-1.6163874602061153E-03) P5=CMPLX(0.9955571458898052,-4.4232022746910895E-03) P6=CMPLX(0.9983567357277001,-6.0591978630835325E-03) C=(2.0**6)/((1.0-P1)*(1.0-P2)*(1.0-P3)*(1.0-P4)*(1.0-P5)* &(1.0-P6)) C=1.0/C H=C*((Z+1.0)**6)/((Z-P1)*(Z-P2)*(Z-P3)*(Z-P4)*(Z-P5)*(Z-P6)) ENDIF C IF(NPOLE.EQ.7)THEN P1=CMPLX(0.9985841093252858,6.1170602668453920E-03) P2=CMPLX(0.9960781482436489,4.8931910342132977E-03) P3=CMPLX(0.9943513201895502,2.7108084639421203E-03) P4=CMPLX(0.9937364712450427,0.0) P5=CMPLX(0.9943513201895502,-2.7108084639421202E-03) P6=CMPLX(0.9960781482436489,-4.8931910342132977E-03) P7=CMPLX(0.9985841093252858,-6.1170602668453921E-03) C=(2.0**7)/((1.0-P1)*(1.0-P2)*(1.0-P3)*(1.0-P4)*(1.0-P5)* &(1.0-P6)*(1.0-P7)) C=1.0/C H=C*((Z+1.0)**7)/((Z-P1)*(Z-P2)*(Z-P3)*(Z-P4)*(Z-P5)*(Z-P6)* &(Z-P7)) ENDIF C IF(NPOLE.EQ.8)THEN P1=CMPLX(0.9987560050814114,6.1548708872084610E-03) P2=CMPLX(0.9965017443899614,5.2060705493843435E-03) P3=CMPLX(0.9947832707928261,3.4725862864172033E-03) P4=CMPLX(0.9938557091182976,1.2182731326107853E-03) P5=CMPLX(0.9938557091182976,-1.2182731326107851E-03) P6=CMPLX(0.9947832707928261,-3.4725862864172032E-03) P7=CMPLX(0.9965017443899614,-5.2060705493843433E-03) P8=CMPLX(0.9987560050814114,-6.1548708872084608E-03) C=(2.0**8)/((1.0-P1)*(1.0-P2)*(1.0-P3)*(1.0-P4)*(1.0-P5)* &(1.0-P6)*(1.0-P7)*(1.0-P8)) C=1.0/C H=C*((Z+1.0)**8)/((Z-P1)*(Z-P2)*(Z-P3)*(Z-P4)*(Z-P5)*(Z-P6)* &(Z-P7)*(Z-P8)) ENDIF C IF(NPOLE.EQ.9)THEN P1=CMPLX(0.9988904149228698,6.1809454234714143E-03) P2=CMPLX(0.9968485891069061,5.4243216505241300E-03) P3=CMPLX(0.9951902432714839,4.0193813261956157E-03) P4=CMPLX(0.9941108076665426,2.1363484214415723E-03) P5=CMPLX(0.9937364712450427,0.0) P6=CMPLX(0.9941108076665426,-2.1363484214415725E-03) P7=CMPLX(0.9951902432714839,-4.0193813261956156E-03) P8=CMPLX(0.9968485891069061,-5.4243216505241299E-03) P9=CMPLX(0.9988904149228698,-6.1809454234714143E-03) C=(2.0**9)/((1.0-P1)*(1.0-P2)*(1.0-P3)*(1.0-P4)*(1.0-P5)* &(1.0-P6)*(1.0-P7)*(1.0-P8)*(1.0-P9)) C=1.0/C H=C*((Z+1.0)**9)/((Z-P1)*(Z-P2)*(Z-P3)*(Z-P4)*(Z-P5)*(Z-P6)* &(Z-P7)*(Z-P8)*(Z-P9)) ENDIF C IF(NPOLE.EQ.10)THEN P1=CMPLX(0.9989983450694978,6.1996946500145186E-03) P2=CMPLX(0.9971359427383281,5.5823988071140633E-03) P3=CMPLX(0.9955571458898052,4.4232022746910894E-03) P4=CMPLX(0.9944132149472717,2.8366075603253747E-03) P5=CMPLX(0.9938128685505996,9.7683826023020200E-04) P6=CMPLX(0.9938128685505996,-9.7683826023020184E-04) P7=CMPLX(0.9944132149472717,-2.8366075603253746E-03) P8=CMPLX(0.9955571458898052,-4.4232022746910895E-03) P9=CMPLX(0.9971359427383281,-5.5823988071140637E-03) P10=CMPLX(0.9989983450694978,-6.1996946500145186E-03) C=(2.0**10)/((1.0-P1)*(1.0-P2)*(1.0-P3)*(1.0-P4)*(1.0-P5)* &(1.0-P6)*(1.0-P7)*(1.0-P8)*(1.0-P9)*(1.0-P10)) C=1.0/C H=C*((Z+1.0)**10)/((Z-P1)*(Z-P2)*(Z-P3)*(Z-P4)*(Z-P5)*(Z-P6)* &(Z-P7)*(Z-P8)*(Z-P9)*(Z-P10)) ENDIF C X=X*CONJG(H) RETURN END c c subroutine ftfwdpad(narray,dt,pulse,ans,omega,nfreq) c routine takes times series in array pulse (which is length c narray), pads it with zeros to lenght 4096*2, and then c flips to the frequency domain. It creates a complex array ans c of length 4097 that holds the positive frequency fourier c transform of the padded time series, and also returns in c array omega, the corresponding RADIAN frequencies. It also c returns the value 4097(number of pos frequencies in padded c transform) in variable nfreq) c narray must of coarse, be less than or equal to 4096, to avoid c wrap-around. After manipulating the answer, use subroutine c ftinvpad to flip to the time domain, and get the time series c of unpadded length narray. c complex ans(4096*2/2+1) real pulse(narray),dum(4096*2*2),omega(4096*2/2+1) if(narray.gt.4096)write(*,*)'error 1 in ftpadforward' nn=4096*2 ! number of real time samples in time series c four1 operates using nn*2 do 1 i=1,nn*2 1 dum(i)=0. ic=0 do 2 i=1,narray*2,2 ic=ic+1 dum(i)=pulse(ic) 2 continue call four1(dum,nn,1) do 3 i=1,nn*2 dum(i)=dum(i)*dt 3 continue c ic=0 df=1./(nn*dt) do 4 i=1,(nn/2+1)*2,2 ic=ic+1 ans(ic)=cmplx(dum(i),dum(i+1)) omega(ic)=(ic-1)*df*2.*3.141592654 4 continue nfreq=nn/2 +1 return end subroutine ftinvpad(narray,dt,ftpad,ans) c routine takes the complex array ftpad, which c represents the positive frequency fourier transform of a time c series that has been padded to a lenght of 4096*2, by c subroutine ftfwdpad, and returns c the inverse transform of the original, unpadded pulse which c is lenght narray. Narray MUST be less than 4096 for these c routines to work correctly, and to avoid wrap-around. c complex ftpad(4096*2/2+1) real ans(narray),dum(4096*2*2) if(narray.gt.4096)write(*,*)'error 1 in ftpadforward' nn=4096*2 ! number of real time samples in padded time series c four1 operates using nn*2 c create the frequency domain representation in wrap-around c order, in array dum c nposf=nn/2+1 ! number of positive frequencies ic=0 do 1 i=1,nposf*2,2 ! create the array for 0 to nyquist ic=ic+1 dum(i)=real(ftpad(ic)) dum(i+1)=aimag(ftpad(ic)) 1 continue c ic=(nposf*2)-1 ! index of the real part at the nyquist do 2 i=(nposf*2)+1,nn*2,2 ! create array for neg. freqs ic=ic-2 dum(i)=dum(ic) dum(i+1)=-dum(ic+1) 2 continue c c flip back to the time domain and return the answer c in the unpadded array ans. c df=1./(nn*dt) ic=0 call four1(dum,nn,-1) do 3 i=1,narray*2,2 ic=ic+1 ans(ic)=dum(i)*df 3 continue c return end subroutine ftmul(a,b,n) c multiplies two frequency domain functions together. c note that the b array gets overwritten with the product complex a(32768),b(32768) do 1 i=1,n b(i)=a(I)*b(I) 1 continue return end c subroutine timeshift(ang,n,h,vp,vs,den,itype,shift) c figures the time shift in the haskell matrix method for the c arrive of the direct S wave or direct P at the surface. c note that in the case of incident SV, it aint easy to determine the c time of the converted P that arrives ahead of the direct S, c so I dont do it here. The answers here are only for the Direct c S wave when either sv or sh are incident, or the direct p when c p is incident. c c n is the number of layers + 1. c the logic is consistent with subroutine qlayer. c c the timeshift is the difference in time between travel time of c the direct ray through the stack, and the horzontal time of the c wave across the top of the half space. c real h(200),vp(200),vs(200),den(200) c hor=0. ray=0. if(itype.ge.2)then ! incident SH or SV c=vs(n)/sin(ang) do 2 i=1,n-1 hlen=h(I)*tan(asin(vs(I)/c)) rlen=sqrt(h(i)**2+hlen**2) ray=ray+rlen/vs(i) hor=hor+hlen 2 continue shift=ray-hor/c go to 20 endif c c=vp(n)/sin(ang) ! incident p do 22 i=1,n-1 hlen=h(I)*tan(asin(vp(I)/c)) rlen=sqrt(h(i)**2+hlen**2) ray=ray+rlen/vp(i) hor=hor+hlen 22 continue shift=ray-hor/c 20 continue return end c c subroutine rotate(fp,fn,narray,strike) real fp(narray),fn(narray) c rotate: rotates two fault-referenced components to north and c south, east and west. c fp,fn: fault parallel and fault normal components. fault para. c motion is positive in the direction of fault strike: c fn motion is positive in the horizontal direction of c fault dip. They get rotated into NS and EW, respectively. c c narray: array size c strike: fault strike, in degrees c pi=3.141592654 s=strike*pi/180. do 1 i=1,narray rns=fp(i)*cos(s)+fn(i)*cos(pi/2.+s) ew=fp(I)*sin(s)+fn(I)*sin(pi/2.+s) fp(i)=rns fn(I)=ew 1 continue return end c subroutine rspectsub(stt,sdata,npts,sosfrq,zh,dmax) IMPLICIT REAL*8 (A-H,O-Y) DOUBLE PRECISION DATA(65536),TT(65536),OSFRQ DOUBLE PRECISION DIS(65536),VEL(65536),ABAC(65536) DOUBLE PRECISION RELAC(65536),XBAR(2),ABAR(2),A(2,2),B(2,2) REAL*4 dmax,vmax,abmax,stt(65536),sdata(65536),sosfrq,zh real psa real relmax,psv c DO 4721 I=1,npts data(i)=dble(sdata(i)) tt(i)=dble(stt(i)) 4721 CONTINUE C C READ THE ACCELERATION TIME SERIES, INTO DATA ARRAY C C N=npts ! NUMBER OF DATA POINTS DT=TT(2)-TT(1) C OSFRQ=DBLE(sOSFRQ) H=DBLE(ZH) PI2=8.0D00*DATAN(1.0D00) C C WS=PI2*OSFRQ C C SET UP THE A AND B MATRICES C EXPO=DEXP(-H*WS*DT) SQR=DSQRT(1.-H*H) ARG=WS*SQR*DT A(1,1)=EXPO*((H/SQR)*DSIN(ARG)+DCOS(ARG)) A(1,2)=(EXPO/(WS*SQR))*DSIN(ARG) A(2,1)=(-WS/SQR)*EXPO*DSIN(ARG) A(2,2)=EXPO*(DCOS(ARG)-(H/SQR)*DSIN(ARG)) P1=((2.*H*H)-1.)/(WS*WS*DT) P2=((2.*H)/(WS*WS*WS*DT))+1./(WS*WS) P3=(2.*H)/(WS*WS*WS*DT) B(1,1)=EXPO*((P1+(H/WS))*DSIN(ARG)/(WS*SQR)+P2*DCOS(ARG))-P3 B(1,2)=-EXPO*(P1*DSIN(ARG)/(WS*SQR)+P3*DCOS(ARG)) &-(1./(WS*WS))+P3 B(2,1)=EXPO*((P1+(H/WS))*(DCOS(ARG)-(H/SQR)*DSIN(ARG)) &-P2*(WS*SQR*DSIN(ARG)+H*WS*DCOS(ARG))) &+(1./(WS*WS*DT)) B(2,2)=-EXPO*(P1*(DCOS(ARG)-(H/SQR)*DSIN(ARG)) &-P3*(WS*SQR*DSIN(ARG)+H*WS*DCOS(ARG))) &-(1./(WS*WS*DT)) C XBAR(1)=0. XBAR(2)=0. ABAR(1)=0. WSSQ=WS*WS HWSTWO=2.*H*WS C DO 5 J=1,N ! FOR EACH OS. FREQUENCY, INTEGRATE THE RESPONSE ABAR(2)=DATA(J) DIS(J)=(A(1,1)*XBAR(1)+A(1,2)*XBAR(2)) &+(B(1,1)*ABAR(1)+B(1,2)*ABAR(2)) VEL(J)=(A(2,1)*XBAR(1)+A(2,2)*XBAR(2)) &+(B(2,1)*ABAR(1)+B(2,2)*ABAR(2)) XBAR(1)=DIS(J) XBAR(2)=VEL(J) ABAR(1)=DATA(J) ABAC(J)=-(HWSTWO*VEL(J)+WSSQ*DIS(J)) RELAC(J)=ABAC(J)-DATA(J) C C ARRAYS DIS AND VEL HOLD THE RELATIVE DISPLACEMENTS AND VELOCITIES; C ABAC HOLDS THE ABSOLUTE ACCELERATION, RELAC HOLDS THE RELATIVE C ACCELERATION C 5 CONTINUE C C NOW FIND MAXIMUM ABSOLUTE VALUES C RELMAX=0. ABMAX=0. VMAX=0. DMAX=0. DO 8 J=1,N IF(DABS(DIS(J)).GE.DMAX)THEN ITSD=J DMAX=sngl(DABS(DIS(J))) ! SD ENDIF IF(DABS(VEL(J)).GE.VMAX)THEN VMAX=sngl(DABS(VEL(J))) ! SV ITSV=J ENDIF IF(DABS(ABAC(J)).GE.ABMAX)THEN ABMAX=sngl(DABS(ABAC(J))) ! SA ITSA=J ENDIF IF(DABS(RELAC(J)).GE.RELMAX)THEN RELMAX=sngl(DABS(RELAC(J))) ! RA ITRA=J ENDIF 8 CONTINUE C C PSV=WS*DMAX ! PSEUDO VELOCITY PSA=WS*PSV ! PSEUDO ACCELERATION TSD=(ITSD-1)*DT ! TIME OF SD TSV=(ITSV-1)*DT ! TIME OF SV TSA=(ITSA-1)*DT ! TIME OF SA TRA=(ITRA-1)*DT ! TIME OF RA C C return END SUBROUTINE TRAPAZ(SIGNAL,DT,N) implicit real*8 (a-h,o-z) DOUBLE PRECISION SIGNAL(65536),DUMMY1(65536) C C INTEGRATES SIGNAL USING THE TRAPAZOID RULE C DO 1 I=1,N DUMMY1(I)=SIGNAL(I) 1 continue SIGNAL(1)=DT*(0.5*DUMMY1(1)) DO 62 J=1,N-1 JP1=J+1 SIGNAL(JP1)=SIGNAL(J)+DT*(0.5*(DUMMY1(J)+DUMMY1(JP1))) 62 CONTINUE RETURN END subroutine interpolate(slat,slon,filein,ans) c Does Bicubic Spline interpolation of a z(longitude,latitude) array. c Designed to work using an array Z derived from the 1247 site hazard c calculation for South Carolina (SCDOT). c Based on Numerical Recipes (Press et al) pp 101. c c Generate the ya=z array. It is derived from the 1247 line *.dat c file containing longitude (negaitive, in degress), latitude (deg) c and motion values for 0.5hz-pga (8 values). c typically the name of the file is something like 2percent_soil.dat c c the main routine passes down the site latitude and longitude c and the name of a file (filein) containing c all the results for a specific probability of exceedance (e.g., c 2percent_soil. The routine returns the interpolated values for c location slat,slon. c real rlon(43),rlat(29),ya(29,43),y2a(29,43) real ya1(29,43),ya2(29,43),ya3(29,43),ya4(29,43),ya5(29,43) real ya6(29,43),ya7(29,43),ya8(29,43) real ans(8) character*66 filein m=29 ! number of latitude values in file (rows) n=43 ! number of longitude values in file (columns) do 1 i=1,m ! loop on latitude do 2 j=1,n ! loop on longitude read(1,*)rlon(j),rlat(i),ya1(i,j),ya2(i,j),ya3(i,j),ya4(i,j) &,ya5(i,j),ya6(i,j),ya7(i,j),ya8(i,j) 2 continue 1 continue close(unit=1) write(3,*)' ' c slon=-abs(slon) c do 44 ik=1,8 ! loop over oscillator frequency and pga c do 45 i=1,m do 46 j=1,n if(ik.eq.1)ya(i,j)=ya1(i,j) ! 0.5hz if(ik.eq.2)ya(i,j)=ya2(i,j) if(ik.eq.3)ya(i,j)=ya3(i,j) if(ik.eq.4)ya(i,j)=ya4(i,j) if(ik.eq.5)ya(i,j)=ya5(i,j) if(ik.eq.6)ya(i,j)=ya6(i,j) if(ik.eq.7)ya(i,j)=ya7(i,j) if(ik.eq.8)ya(i,j)=ya8(i,j) ! pga 46 continue 45 continue call splie2(rlat,rlon,ya,m,n,y2a) c c y2a is array with second-derivatives of ya c c call splin2(rlat,rlon,ya,y2a,m,n,slat,slon,an) ans(ik)=an 44 continue c c c identfiy the nearest grid point and calculate the distance and c azimuth to it, from the site of interest. c dmin=100000. do 11 i=1,m-1 ! loop on latitude do 22 j=1,n-1 ! loop on longitude c if(slat.eq.rlat(i).and.slon.eq.rlon(j))then dmin=0.0 ilat=i ilon=j goto 5393 endif if(slat.ne.rlat(i).and.slon.ne.rlon(j))then call disaz(slat,slon,rlat(i),rlon(j),d,az) if(d.lt.dmin)then dmin=d ilat=i ilon=j endif endif 22 continue 11 continue 5393 dmin=dmin*111.199 write(*,*) write(*,*)'RESULTS OF INTERPOLATION' write(*,*) write(*,107)'Site Location:',slat,'N',-slon,'W' write(3,*)' RESULTS OF INTERPOLATION' write(3,*) write(3,107)'Site Location:',slat,'N',-slon,'W' 107 format(a20,1x,f7.4,1x,a1,1x,f7.4,1x,a1) write(*,101)'Nearest Grid Point:',rlat(ilat), &'N',-rlon(ilon),'W', 'Distance From Site:',dmin, 'Km' write(3,101)'Nearest Grid Point:',rlat(ilat), &'N',-rlon(ilon),'W', 'Distance From Site:',dmin, 'Km' 101 format(a20,1x,f7.4,1x,a1,1x,f7.4,1x,a1,1x,a19,1x,f6.2,1x,a2) write(*,*) 20 return end C Begin subroutine interphazard inserted by Jake Beale (Nov 2005) C This subroutine is a modified copy of subroutine interpolate subroutine interphazard(slat,slon,filein2,ans) c Does Bicubic Spline interpolation of a z(longitude,latitude) array. c Designed to work using an array Z derived from the 1247 site hazard c calculation for South Carolina (SCDOT). c Based on Numerical Recipes (Press et al) pp 101. c c Generate the ya=z array. It is derived from the 1247 line *.dat c file containing longitude (negaitive, in degress), latitude (deg) c and motion values for 0.5hz-pga (8 values). c typically the name of the file is something like 2percent_soil.dat c c the main routine passes down the site latitude and longitude c and the name of a file (filein) containing c all the results for a specific scenario earthquake (e.g., c scenario-eq2. The routine returns the interpolated values for c location slat,slon. c real rlon(43),rlat(29),ya(29,43),y2a(29,43) real ya1(29,43),ya2(29,43),ya3(29,43),ya4(29,43),ya5(29,43) real ya6(29,43) real ans(6) !contains mean and modal eq mag, dist, eps character*66 filein2 m=29 ! number of latitude values in file (rows) n=43 ! number of longitude values in file (columns) do 1 i=1,m ! loop on latitude do 2 j=1,n ! loop on longitude read(4,*)rlat(i),rlon(j),ya1(i,j),ya2(i,j),ya3(i,j),ya4(i,j) &,ya5(i,j),ya6(i,j) 2 continue 1 continue close(unit=4) c slon=-abs(slon) c do 44 ik=1,6 ! loop over oscillator frequency and pga c do 45 i=1,m do 46 j=1,n if(ik.eq.1)ya(i,j)=ya1(i,j) !mean dist if(ik.eq.2)ya(i,j)=ya2(i,j) !mean mag if(ik.eq.3)ya(i,j)=ya3(i,j) !mean eps if(ik.eq.4)ya(i,j)=ya4(i,j) !mod dist if(ik.eq.5)ya(i,j)=ya5(i,j) !mod mag if(ik.eq.6)ya(i,j)=ya6(i,j) !mod eps 46 continue 45 continue call splie2(rlat,rlon,ya,m,n,y2a) c c y2a is array with second-derivatives of ya c c call splin2(rlat,rlon,ya,y2a,m,n,slat,slon,an) ans(ik)=an 44 continue 20 return end C End of subroutine interphazard inserted by Jake Beale (Nov 2005) SUBROUTINE SPLIE2(X1A,X2A,YA,M,N,Y2A) c PARAMETER (NN=100) PARAMETER (NN=1300) DIMENSION X1A(M),X2A(N),YA(M,N),Y2A(M,N),YTMP(NN),Y2TMP(NN) DO 13 J=1,M DO 11 K=1,N YTMP(K)=YA(J,K) 11 CONTINUE CALL SPLINE(X2A,YTMP,N,1.E30,1.E30,Y2TMP) DO 12 K=1,N Y2A(J,K)=Y2TMP(K) 12 CONTINUE 13 CONTINUE RETURN END SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2) c PARAMETER (NMAX=100) PARAMETER (NMAX=1300) DIMENSION X(N),Y(N),Y2(N),U(NMAX) IF (YP1.GT..99E30) THEN Y2(1)=0. U(1)=0. ELSE Y2(1)=-0.5 U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) ENDIF DO 11 I=2,N-1 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2. Y2(I)=(SIG-1.)/P U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 11 CONTINUE IF (YPN.GT..99E30) THEN QN=0. UN=0. ELSE QN=0.5 UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) ENDIF Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.) DO 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE RETURN END SUBROUTINE SPLIN2(X1A,X2A,YA,Y2A,M,N,X1,X2,Y) c PARAMETER (NN=100) PARAMETER (NN=1300) DIMENSION X1A(M),X2A(N),YA(M,N),Y2A(M,N),YTMP(NN),Y2TMP(NN),YYTMP( *NN) DO 12 J=1,M DO 11 K=1,N YTMP(K)=YA(J,K) Y2TMP(K)=Y2A(J,K) 11 CONTINUE CALL SPLINT(X2A,YTMP,Y2TMP,N,X2,YYTMP(J)) 12 CONTINUE CALL SPLINE(X1A,YYTMP,M,1.E30,1.E30,Y2TMP) CALL SPLINT(X1A,YYTMP,Y2TMP,M,X1,Y) RETURN END c subroutine disaz(oelat,oelon,oslat,oslon,odel,oaz) implicit real*8 (a-h,p-z) real*8 k,kp c disaz: calculates distance in degrees and azimuth in degrees for two c points on spherical earth, according to bullens exact equations. c the azimuth az is reckoned from epicenter to station, east of north. c elat and slon are geographic latitudes, in decimal degrees. c elon and slon are east longitude, in decimal degrees: (west longitudes c are to be sent down as negatives. c elat and elon refer to the "epicenter", whereas slat and slon refer to c the "station". c NOTE: calling program provides the epicenter and station coordinates in c SINGLE PRECISION, and the azimuth and distance are returned in c single precision. c elat=dble(oelat) elon=dble(oelon) slat=dble(oslat) slon=dble(oslon) if(elon.lt.0.0)elon=360.0+elon if(slon.lt.0.0)slon=360.0+slon c c get into radian mode: c rmin1=-1.00 pi=dacos(rmin1) elatrad=elat*pi/180.0 elonrad=elon*pi/180.0 slatrad=slat*pi/180.0 slonrad=slon*pi/180.0 c c calculate geocentric colatitudes: c scl=datan(.993277*dtan(slatrad)) ecl=datan(.993277*dtan(elatrad)) scl=pi/2.-scl ecl=pi/2.-ecl c c use bullens variable naming scheme c theta=ecl thetap=scl phi=elonrad phip=slonrad a=dsin(theta)*dcos(phi) ap=dsin(thetap)*dcos(phip) b=dsin(theta)*dsin(phi) bp=dsin(thetap)*dsin(phip) c=dcos(theta) cp=dcos(thetap) d=dsin(phi) dp=dsin(phip) e=-dcos(phi) ep=-dcos(phip) g=dcos(theta)*dcos(phi) gp=dcos(thetap)*dcos(phip) h=dcos(theta)*dsin(phi) hp=dcos(thetap)*dsin(phip) k=-dsin(theta) kp=-dsin(thetap) c cosdel=a*ap+b*bp+c*cp sinz=(((ap-d)**2+(bp-e)**2+cp**2)-2.)/(2.*dsin(dacos(cosdel))) cosz=(((ap-g)**2+(bp-h)**2+(cp-k)**2)-2.)/(2.*dsin(dacos(cosdel))) c c z=datan2(sinz,cosz)*180.0/pi if(z.lt.0.)z=360.+z oaz=sngl(z) del=dacos(cosdel)*180./pi odel=sngl(del) return end C********************************************** function cmn2m(xmn) c c converts from mblg (Mn) to moment magnitude c using the atkinson and boore SRL 1987 relation c xmn=mblg,cmn2m=moment mag: also can use johnston relation: c watch out for comments! c c cmn2m=2.689-.252*xmn+.127*xmn*xmn ! atkinson and boore (srl) c cmn2m=2.715-.277*xmn+.127*xmn*xmn ! atkinson and boore bssa c cmn2m=5.4867-1.433*xmn+.2533*xmn*xmn ! arch cmn2m=3.45-.473*xmn+.145*xmn*xmn ! frankel et al 1996 (arch) c open file report on maps. return end C********************************************** subroutine mcc2shake(n,a,t,dt,fileout) c mcc2shake: writes out file in shake91 input format c to unit 3, acceleration values in g. c maximum number of points is 4096 character*80 fileout real a(4096),t(4096) open(unit=3,name=fileout,type='new',err=9002) do 555 i=1,n a(i)=a(i)/981. 555 continue dt=t(2)-t(1) nlines=n/8 c c write(3,100)dt,n 100 format(1x,'sample interval (sec)',1x,f6.4,' npts ',i7) k=-7 do 1 i=1,nlines k=k+8 write(3,88)(a(j),j=k,k+7),i 88 format(8(f9.6),t73,i7) 1 continue c close(unit=3) return 9002 write(*,*)'I/O ERROR: Cannot open an output file.',fileout write(*,*)'Make sure output files do not already exist.' write(*,*)'Press Enter to terminate program.' read(*,*) stop end c subroutine interpdepth(slat,slon,ans) c Does Bicubic Spline interpolation of a z(longitude,latitude) array. c Designed to work using an array Z derived from the 1247 site hazard c calculation for South Carolina (SCDOT). c Based on Numerical Recipes (Press et al) pp 101. c c Generate the ya=z array. It is derived from the 1247 line scenario-sites c file containing longitude (negaitive, in degress), latitude (deg) c and thickness of sediments, in meters C c NOTE : file scenario-sites has to exist in the same directory as this routine c real rlon(43),rlat(29),ya(29,43),y2a(29,43) character*80 filein filein='scenario-data/scenario-sites' ! just renamed version of old sites.dat' open(unit=1,name=filein,type='old',err=9003) isites=1 m=29 ! number of latitude values in file (rows) n=43 ! number of longitude values in file (columns) do 1 i=1,m ! loop on latitude do 2 j=1,n ! loop on longitude if(isites.eq.1)read(1,*)rlat(i),rlon(j),ya(i,j) 2 continue 1 continue close(unit=1) c call splie2(rlat,rlon,ya,m,n,y2a) c c y2a is array with second-derivatives of ya c slon=-abs(slon) c call splin2(rlat,rlon,ya,y2a,m,n,slat,slon,ans) c if(ans.lt.1.)ans=0. 20 return 9003 write(*,*)'I/O ERROR: Cannot open an input file. Exiting.' write(*,*)'Make sure all input files are in "scenario-data".' write(*,*)'Press Enter to terminate program.' read(*,*) stop end c subroutine newshakesub(qsed,sedth,freqs,asngl,bsngl) c c derived from newshake to do the SCDOT boore and joyner approximation. c c modified 2/2006 to allow either qsed=32 (as in the original version of the c code) or qsed=100, as per Langston's BSSA 2006 result c note that this routine is only called for the isoil=1 case (geologically c realistic). (This is different from original version). c c calling routine must provide the freqs c array with the 4096 frequency values, c the thickness of the sedimentary layer, sedth (zero is ok) c and the value of qsed, which is the value of Q in the sedimentary layer c c program returns the real and imaginery parts of the b and j 1991 c 1/4 simplification transfer function in arrays asngl and bsngl, which c must be dimensioned 4096 in the calling routine C*************************************************************************** C IOPT=3 IS COMPLETELY DIFFERENT !!!!!!!!!!!!!!!!!!!!!!!!!!!!! C*************************************************************************** C IOPT 3 IS NOT, REPEAT NOT, A RIGOROUS SOLUTION OF THE EOM. C IT IS THE JOYNER AND BOORE (1981), BOORE AND JOYNER (1991), BSSA C VO1 81, NO. 6 PP 2167-2185 APPROXIMATE SOLUTION FOR VERTICAL C INCIDENCE, FOR THE TRANSFER FUNCTION OF SOIL SURFACE OVER BASEMENT C OUTCROP. C C Note that the dimensions of measurement for density do not have to c be consistent with velocity and thickness (the latter two have to c be consistent) C c c note the this version makes minimum phase transfer function for c iopt 3. C REAL*8 RHO(100),V(100),BETA(100),H(99) REAL*8 FFREQ(32768),A(32768),B(32768) REAL*8 SUMT,SUMD,KAPPA,TBOTLAY(100),DBOTLAY(100),TTARG,Z,VEFF REAL*8 TH,T,XIMAG,XREAL real*4 asngl(4096),bsngl(4096),freqs(4096) do 4304 i=1,4096 ffreq(i)=dble(freqs(i)) 4304 continue C if(sedth.eq.0.)nl=1 if(sedth.ne.0.)nl=2 c if(nl.eq.1)then rho(1)=2.5 v(1)=2.5 h(1)=0.25 beta(1)=1./(2.*600.) ! attenuation in weathered rock layer v(2)=3.5 rho(2)=2.6 endif c if(nl.eq.2)then v(1)=.7 rho(1)=2. beta(1)=1./(2.*qsed) !qsed will be either 32 or 100, from main h(1)=dble(sedth)/1000. ! thickness in kilometers v(2)=2.5 rho(2)=2.5 h(2)=.25 beta(2)=1./(2.*600.) v(3)=3.5 rho(3)=2.6 endif C C IOPT=3 NF=4096 ! TOTAL NUMBER OF FREQUENCIES ipos=nf/2+1 ! number of positive freqs. C C IF(IOPT.EQ.3)THEN ! SIMPLIFIED METHOD J&B DO 6111 I=1,NF SUMT=0. SUMD=0. DO 6112 J=1,NL SUMT=SUMT+H(J)/V(J) SUMD=SUMD+H(J) TBOTLAY(J)=SUMT ! TIME TO BOTTOM OF LAYER DBOTLAY(J)=SUMD ! DEPTH TO BOTTOM OF LAYER 6112 CONTINUE IF(FFREQ(I).NE.0.0)TTARG=DABS(1./(4.*FFREQ(I))) IF(FFREQ(I).EQ.0.0)GOTO 600 DO 6113 J=1,NL IF(TTARG.LE.TBOTLAY(J))THEN IF(J.EQ.1)Z=TTARG*V(1) IF(J.NE.1)Z=(TTARG-TBOTLAY(J-1))*V(J)+DBOTLAY(J-1) VEFF=Z/TTARG GOTO 1278 ENDIF 6113 CONTINUE 1278 IF(TTARG.GT.TBOTLAY(NL))THEN Z=DBOTLAY(NL)+(TTARG-TBOTLAY(NL))*V(NL+1) VEFF=Z/TTARG ENDIF 600 IF(FFREQ(I).EQ.0.0)VEFF=V(NL+1) A(I)=DSQRT(V(NL+1)/VEFF) ! AMPLIFICATION 6111 CONTINUE C C NOW GET KAPPA C KAPPA=0. DO 6114 I=1,NL KAPPA=KAPPA+H(I)*2.*BETA(I)/V(I) 6114 CONTINUE TH=DBOTLAY(NL) T=TBOTLAY(NL) ANGI=0. C C MAKE COMPLETE TRANSFER FUNCTION C FIRST GET PHASE FROM LOG MODULUS OF TRANSFER FUNCTION C DO 6115 I=1,NF A(I)=A(I)*DEXP(-3.141592654*KAPPA*DABS(FFREQ(I))) A(I)=DLOG(A(I)) ! NAT LOG OF MODULUS 6115 CONTINUE C C GET MINIMUM PHASE FROM HILBERT TRANSFORM OF LOG MODULUS C CALL HILAPHAD(A,B,NF) ! MINIMUM PHASE IS IN ARRAY B C C MAKE CAUSAL T.F. C DO 4124 I=1,NF A(I)=DEXP(A(I)) XREAL=A(I)*DCOS(B(I)) XIMAG=A(I)*DSIN(B(I)) A(I)=XREAL B(I)=XIMAG asngl(i)=sngl(a(i)) bsngl(i)=sngl(b(i)) 4124 CONTINUE C ENDIF C END OF J&B SIMPLIFIED OPTION C return END c SUBROUTINE HILAPHAD(LNMOD,PHASE,NN) c this version unchanged from routine in newshake.f c uses double precision from the calling routine COMPLEX P(32768),Z REAL*8 LNMOD(1),PHASE(1) DIMENSION DATA(65536) C COMPUTES PHASE FROM LOG AMP SPECTRUM USING FREQUENCY DOMAIN REALIZATION OF C HILBERT TRANSFROM. C C LNMOD CONTAINS THE NATURAL LOG OF THE MODULUS OF THE FOURIER TRANSFORM OF C SOME TIME SERIES IN WRAP AROUND ORDER: ZERO TO NYQUIST, -NYQUIST TO ZERO. C NN IS THE SAME AS IN FOUR1: THE ACTUAL NUMBER OF TIME SERIES DATA POINTS. C C PUT THE LOG OF THE MODULUS IN THE REAL PART OF DATA ARRAY (ODD ELEMENTS) C DO 55 I=1,NN*2 55 DATA(I)=0. IC=0 DO 3 I=1,NN*2,2 IC=IC+1 DATA(I)=SNGL(LNMOD(IC)) 3 DATA(I+1)=0. C C TRANSFORM C CALL FOUR1(DATA,NN,-1) DO 4 I=1,NN*2 4 DATA(I)=DATA(I)/NN C IC=0 DO 5 I=1,NN*2,2 IC=IC+1 5 P(IC)=CMPLX(DATA(I),DATA(I+1)) Z=(0.0,1.0) DO 6 I=1,NN/2 ! FOR POSITIVE FREQS 6 P(I)=-P(I)*Z DO 7 I=NN/2+2,NN ! NEG FREQS 7 P(I)=P(I)*Z P(1)=(0.0,0.0) ! ZERO FREQ P(NN/2+1)=(0.0,0.0) ! NYQUIST C IC=0 DO 8 I=1,NN*2,2 IC=IC+1 DATA(I)=REAL(P(IC)) 8 DATA(I+1)=AIMAG(P(IC)) C CALL FOUR1(DATA,NN,1) C IC=0 DO 10 I=1,NN*2,2 IC=IC+1 10 PHASE(IC)=DBLE(DATA(I)) ! PHASE C RETURN END C subroutine layerresp(sfreq,nl,sv,sbeta,sh,vhalf,ans) REAL*8 V(100),BETA(99),H(99) REAL*8 A,SUMT,SUMD,KAPPA,TBOTLAY(100),DBOTLAY(100),TTARG,Z,VEFF REAL*8 PI2,FREQ,RONE real sfreq,sv(100),sbeta(99),sh(99),vhalf,ans C c sfreq = frequency in Hertz. c nl = number of layers, c sv = array with layer shear wave velocities c sbeta = array with 1/2Q in each layer c sh = array with layer thickness c vhalf = shear velcity in half space (Q is infinite) c ans = modulus of transfer function C c program allows up to 99 layers. c temp: c standard input for SCDOT project c nl=2 c sv(1)=.7 ! 700 m/sec c sv(2)=2.5 ! 2.5 km/sec c sbeta(1)=1./(2.*32.) ! shallow layer damping, Q=32 c sbeta(2)=1./(2.*600.) ! deep weathering layer, Q=600 c sh(1) is determined on input by the main routine and sent down in call to sub c sh(2)=0.25 ! second layer thickness in km. c vhalf=3.5 ! 3.5 km/sec half space c sfreq is determined on input by the main routine and sent down in call to sub C c standard input for no site response correction: c n1=1 c sv(1)=1. c sh(1)=1. c sbeta(1)=0. c vhalf=1. C c NOTE that this is not an exact solution and is simply the JB quarter wavelength c approximation C RONE=1.0 pi2=DABS(8.*DATAN(RONE)) freq=dble(sfreq) sv(nl+1)=vhalf ! puts the half space velocity in the nl+1 array position. v(nl+1)=dble(sv(nl+1)) C C C SUMT=0. SUMD=0. DO 6112 J=1,NL v(j)=dble(sv(j)) h(j)=dble(sh(j)) beta(j)=dble(sbeta(j)) SUMT=SUMT+H(J)/V(J) SUMD=SUMD+H(J) TBOTLAY(J)=SUMT ! TIME TO BOTTOM OF LAYER DBOTLAY(J)=SUMD ! DEPTH TO BOTTOM OF LAYER 6112 CONTINUE IF(FREQ.NE.0.0)TTARG=DABS(1./(4.*FREQ)) IF(FREQ.EQ.0.0)GOTO 600 DO 6113 J=1,NL IF(TTARG.LE.TBOTLAY(J))THEN IF(J.EQ.1)Z=TTARG*V(1) IF(J.NE.1)Z=(TTARG-TBOTLAY(J-1))*V(J)+DBOTLAY(J-1) VEFF=Z/TTARG GOTO 1278 ENDIF 6113 CONTINUE 1278 IF(TTARG.GT.TBOTLAY(NL))THEN Z=DBOTLAY(NL)+(TTARG-TBOTLAY(NL))*V(NL+1) VEFF=Z/TTARG ENDIF 600 IF(FREQ.EQ.0.0)VEFF=V(NL+1) A=DSQRT(V(NL+1)/VEFF) ! AMPLIFICATION C C NOW GET KAPPA C KAPPA=0. DO 6114 I=1,NL KAPPA=KAPPA+H(I)*2.*BETA(I)/V(I) 6114 CONTINUE tstar=kappa C A=A*DEXP(-3.141592654*KAPPA*DABS(FREQ)) ans=sngl(a) c C return END c subroutine pgaratio1(thick,ratio) real ratio1(6),thick1(6),ratio2(12),thick2(12) data ratio1/1.18,1.35,1.58,1.8,1.95,2.0/ data thick1/0.,.002,.004,.006,.008,.01/ data ratio2/2.0,1.58,1.32,1.13,.98,.849,.734,.592, &.486,.4,.3,.2/ data thick2/.01,.1,.2,.3,.4,.5,.6,.8,1.,1.4,1.8, &2.2/ c c The velocity model consists of 2 layers over half space: c layer 1 velocity =.7 km/sec, Q=32 c layer 2 velocity =2.5 km/sec Q=600., thickness .25 km. c halfspace velocity 3.5 km/sec c c the values tabulated as ratio1 and ratio2 are the pga amplification c factors as functions of layer 1 thickness. This was determined by c comparing the soil rock time series simulations using transfer c function derived from jb quarter wavelength approximation. c c ratio1 contains the ratio of soil/rock pga for layer thickness c 0 to .01 km (0 to 10 meters). c c ratio2 contains the ratio of soil/rock pga for layer thickness c .01 to 2.2 km (10 to 2200 m). c if(thick.le.0.01)then ! less than 10 meters call lininterp(6,thick1,ratio1,thick,ratio) else call lininterp(12,thick2,ratio2,thick,ratio) endif return end c subroutine lininterp(n,x,y,xval,ans) c c does linear interpolation between two points. c calling routine sends down arrays y and x, and c number of points in array n, and xvalue at which c y value is needed. c real x(1000),y(1000) c c fool proof it c if(xval.lt.x(1))write(*,*)'Error 1 in subroutine lininterp' if(xval.gt.x(n))write(*,*)'Error 2 in subroutine lininterp' if(xval.lt.x(1))write(2,*)'Error 1 in subroutine lininterp' if(xval.gt.x(n))write(2,*)'Error 2 in subroutine lininterp' do 1 i=1,n-1 if(xval.ge.x(i).and.xval.le.x(i+1))index=i 1 continue rise=y(index+1)-y(index) run=x(index+1)-x(index) ans=y(index)+(xval-x(index))*rise/run if(xval.eq.x(n))ans=y(n) return end c subroutine pgaratio2(thick,ratio) real ratio1(6),thick1(6),ratio2(12),thick2(12) data ratio1/1.17,1.3,1.49,1.78,1.9,2.0/ data thick1/0.,.002,.004,.006,.008,.01/ data ratio2/2.0,2.04,1.93,1.8,1.66,1.58,1.5,1.36,1.22,.96,.86,.76/ data thick2/.01,.1,.2,.3,.4,.5,.6,.8,1.,1.4,1.8,2.2/ c c The velocity model consists of 2 layers over half space: c layer 1 velocity =.7 km/sec, Q=100, from langston et al. bssa 2006 c layer 2 velocity =2.5 km/sec Q=600., thickness .25 km. c halfspace velocity 3.5 km/sec c c the values tabulated as ratio1 and ratio2 are the pga amplification c factors as functions of layer 1 thickness. This was determined by c comparing the soil rock time series simulations using transfer c function derived from jb quarter wavelength approximation. c c ratio1 contains the ratio of soil/rock pga for layer thickness c 0 to .01 km (0 to 10 meters). c c ratio2 contains the ratio of soil/rock pga for layer thickness c .01 to 2.2 km (10 to 2200 m). c if(thick.le.0.01)then ! less than 10 meters call lininterp(6,thick1,ratio1,thick,ratio) else call lininterp(12,thick2,ratio2,thick,ratio) endif return end c