program EXSIM !! release version 1.0, October 10, 2005, extensively modified by D. Boore c For more detalis see Motazedian,Dariush and Gail M. Atkinson, (2005) c "Stochastic Finite-Fault Modeling Based on a Dynamic Corner Frequency", c Bulletin of the Seismological Society of America,95,995-1010. c c This EXSIM is designed to work for a list of observation sites for a c user-defined fault geometry and location. Reads the input parameters and c outputs acceleration time series for the first trial, and the average Fourier c Spectra and PSA over trials. * !!! NOTE: If the PSA and/or acceleration time series of each realization is * desired, suitable modifications need to be made to the program *!Control file for program exsim_dmb_working *! Revision of program involving a change in the control file on this date: * 05/28/09 *!Title * Runs for comparing EXSIM and SMSIM *!Write acc, psa, husid files for each site? * N *!MW, Stress, flag (0=fmax; 1=kappa), fmax or kappa * 7.0 140.0 1 0.005 *!lat and lon of upper edge of fault * 0.0 0.0 *!strike,dip, depth of fault * 0.0 50.0 0.0 *!fault type (S=strikeslip; R=reverse; N=normal; U=undifferentiated) *! (Only used if Wells and Coppersmith is used to obtain FL and FW). * R *!fault length and width, dl, dw, stress_ref *!Note: Force program to use Wells and Coppersmith for FL and/or FW if either entry = 0.0 * 0 0 1.5 1.5 70.0 ! used for WC scaling (need an entry as a placeholder even if not used_ *!y (=vrup/beta) * 0.8 *!hypo location in along fault and down dip distance from the fault *!reference point (an upper corner)(-1.0, -1.0 for a random location); *!number of iterations over hypocenter (need an entry, but only used if *!either of the first two values are -1.0, indicating a random location) * -1.0 -1.0 10 * 0 0 10 *!Enter type of risetime (1=original, 2=1/f0) * 2 *!tpadl, tpadt, delta t * 20.0 20.0 0.002 *!beta , rho * 3.7 2.8 *!gsprd: r_ref, nsegs, (rlow(i), a_s, b_s, m_s(i)) (Usually set r_ref = 1.0 km) * 1.0 * 3 * 1.0 -1.3 0.0 6.5 * 70.0 +0.2 0.0 6.5 * 140.0 -0.5 0.0 6.5 *!q: fr1, Qr1, s1, ft1, ft2, fr2, qr2, s2, c_q * 1.0 1000.0 0.0 1.4242 1.4242 1.0 893.0 0.32 3.7 *!trilinear duration and properties (rmin,rd1,rd2,durmin,b1,b2,b3) * 10. 70. 130. 0.0 0.16 -0.03 0.04 *!Type of window: 1 for Saragoni-Hart taper windows, 0 for tapered boxcar window * 0 *!low-cut filter corner, nslope (0 ==> no filter) * 0.05 8 *! %damping of response spectra * 5.0 *!# of f and Min and Max F for response spectra * 100 0.1 100. *!no. of frequencies for summary output (10 max): * 2 *!frequency (-1.0, 99.0 for pgv, pga): * 0.5 5.0 *!Output file names stem: * M7.0_60sub_10_10_pulse50_orig_rt_normal *!Name of crustal amplification file: * ab06_hr_crustal_amps.txt *!No. of frequencies in the crustal amplification file * 5 *!Name of site amplification file: * siteAmplification.txt *!No. of frequencies in the site amplification file (0 to skip) * 0 *!DynamicFlag (0=no), PulsingPercent * 1 50.0 *!iflagscalefactor (1=vel^2; 2=acc^2; 3=asymptotic acc^2 (dmb)) * 2 *!iflagfas_avg (1=arithmetic; 2=geometric, 3=rms: USE 3!) * 3 *!iflagpsa_avg (1=arithmetic; 2=geometric: USE 2!, 3=rms) * 2 *!deterministic flag,gama,nu,t0, impulse peak * 0 1.0 90.0 4.0 10. *!iseed, # of trials * 309 10 *!Number of Sites, coord flag (1=lat,long; 2=R,Az; 3=N,E) * 2 3 *!Coordinates of each site * 0.0 02.5 * 0.0 160.0 *!Move the sites to the midpoint and end of the surface projection *!of the upper edge of the fault if isitecoordflag = 3 and *!siteLocation(1) = 0 (center) and siteLocation(2) = 0 (end), respectively *!(this only makes sense if the strike of the fault = 0.0)? * Y *!islipweight=-1 -> unity slip for all subfaults,=0 -> specify slips below, =1 -> random weights * -1 *!Matrix of slip weights * 1.0 1.0 1.0 1.0 1.0 * 1.1 1.3 1.4 1.1 1.0 * 1.1 1.4 1.5 1.2 1.2 * 1.2 1.5 1.6 2.2 1.3 * Dates: 08/17/08 - Modifications by D. M. Boore * 08/17/08 - Corrected typo "ISlipWeigth" and "falg", and allow all * unit slip weights if flag -1. * 08/18/08 - Read in y (=vrup/beta), write R in PSA_FA output * Revise column headers of PSA_FA output. * Lots of changes to output, including writing * separate PSA_FA files for each site (to make it * easier to import into a graphics program). * Allow the site coordinates to be entered as * isitecoordflag coords * 1 lat,long * 2 R, Az * 3 N, E * 08/19/08 - Reverse order of fault input, and allow use of Wells * and Coppersmith if FL=0 or FW=0, put magnitude input before * FL, FW (because the later needs amag). Also get fault * type (needed for Wells and Coppermsith). * 08/20/08 - Change averaging algorithm (log avg for PSA, E(FAS^2)^1/2 for FAS) * Return E(FAS^2)^1/2, not log of the value, out of the * binning routine SAMPLE, called by compute FACCN. * * Add input variable iflagscalefactor to determine whether to use * scale factor based on the following: * iflagscalefactor integrand of integral * 1 velocity spectrum squared * 2 acceleration spectrum squared * 3 D. Boore's factor based on acc. sq. high-frequency asymptotes * Add input variable to determine if use arithmetic or * geometric average for PSA * 08/21/08 - Added scalefactor (h) to computeStochasticWave output and print out * the scalefactor and corner frequency for each subfault for * the first simulation. Also include an option for tapering the * scaling factor. * 08/22/08 - Write out the time series for the first realization only. * 08/25/08 - Read in qmin * In checkFFTpoints, increased upper index so that 2**n=32768. * Also, doubled dimensions of wave, totalWave, and subwave. * In many subroutines replaced dimensions of arrays passed * in through the parameter list with "(*)". * 08/27/08 - Added dur_sub to common /sub/, in order to print it out at some point. * 11/06/08 - Add option to compute RMS of PSA, and also use same options for computing the average of FAS * (for study purposes). * 11/07/08 - Changed way that averages are computed. Now the accumulated sum is computed * in computeAverageResponse and computeAverageFA and this is used in the main program * to compute an average. But a modification was needed to allow * for negative values of log10(y) for geometric means. I set the average to -9.99 * if accum_fas = 0.0 or accum_psa = 0.0. * 11/08/08 - I may be having problems with the rms average of FAS, so compute this differently now. * 11/15/08 - Increase maximum number of nl, nw to 200. * 11/25/08 - I am concerned about the calls to ran1. Note that the distributed version of * exsim uses urand, not ran1. Apparently I susbtituted ran1 for urand, but I think * I made a mistake. If used to * generate a series of "random" numbers, ran1 should be called the * first time with a negative integer (idum) as a seed. * Ran1 performs some initializations, including assigning * values to Saved variables iv and iy, as well as changing * idum. In exsim ran1 is called for a number of things, * possibly including locateHypocentreRandomly, * createRandomWeight, computeStochasticWave, and * computing a random delay for the time series from * each subsource. I would think that the idum, iv, * and iy values would need to be in sync with one another, * but this will not be the case if ran1 is called with * a different value for idum that is not the same as in * previous calls. In the previous version of exsim, iseed * is used as idum in the calls in all calls except for * the call for computing the random delay. But iseed was * not constrained to be less than 0.0 initially, so the ran1 * initialization was probably not done. And also there was * one call to ran1 that used islipweight for idum, which * could have values less than, equal, to, or grater than 0.0 * (in my runs it was lsess than 0.0). When less than 0.0, as * it was in all of my runs, the iv and iy would have been * initialzed and would have been used in the next call to * ran1, but the idum would not have been the reset value * of islipweight. I do not know the consequences of all of * this. I have changed exsim to contrain iseed to be less * than 0 initially and have used iseed in all calls to ran1. * 11/27/08 - Added duration calculation * 11/28/08 - Use 1/f0 for risetime * 11/28/08 - Change input file to have a comment line before the name of the scale * factor file. * 11/30/08 - Move computation of random slipweights * if Islipweight .eq. 1.0 from main program to getInputParameters * (a more logical place because the slipweights array is filled * in the same place for all values of ISlipWeight). * 11/30/08 - Move "findSubFaultMoment" out of loop over isite * 12/01/08 - Use "if then" in loop over subsources in * finding NumberofActiveSubs * 12/01/08 - Remove "NumberofActivesSubs" from the arguments passed * to subroutines computeStochasticWave and sourceSpectra, as * the variable is not used. * 12/01/08 - Added writing of istart, istop, nptsTotalWave, maxPointsOfWave * in WritePar (I'm trying to figure out if findindex is needed and if * code setting istop is correct). * 11/28/08 - Use 1/firstElementF0 for risetime * 12/01/08 - Use original exsim risetime * 12/01/08 - Add option to choose one of three ways of obtaining risetime * 12/06/08 - Allow iteration of random hypocenters for each station * 12/08/08 - Read initial seed. * 12/09/08 - Add scaling_coefficient, following Dariush's suggestion. * 12/09/08 - Increase length of output file names. * 12/17/08 - Move computations that do not depend on the loop over * nsims out of the loop. This required defining f0, etc as * arrays (one value per subsource). This is computationally efficient. * Another advantage is that * I can calculate the maximum duration for the total wave before doing * the subsource wave calculation, and therefore I can use dynamic * allocation for the time series array dimension. This required * changing common sub not to include arrays of variables that * depend on the subsource, passing these things through argument lists. * I renamed common /sub/ to common /params/. I also removed dur_sub * from this common block. * I also compute dur_sub outside of computeStochasticWave and pass the * value into that subroutine through the argument list. * 12/17/08 - Major revision involving the time series (use tpads and front and back, do not truncate * the motions from each subsource; use dynamic allocations). NOTE: * Using Mavroeidis and Papageorgiou has not been tested, and it will probably not work * because I have not taken into account the added npadl. I comment out a section of code * before the call to M&P that may have to do with determining the index for which the pulse * starts. It would be easy to fix this later using t_arrive(i,j). * 12/22/08 - If nl = nw = 1, over-ride the computation of a random hypocenter and * set n_hypocenters, i0, j0 = 1, 1, 1 * 12/31/08 - Write sqrt(nl*nw)*fas to help in deciding if that factor can correct for * the incoherent summation. * 12/31/08 - replace "np" for the taper correction at low frequencies by "ftaper_f0", * the ratio of ftaper and f0 read from the input file. Also, there is no * need for "taper_scalefactor", because ftaper_f0 .eq. 0.0 implies * no tapering toward low frequencies. * 01/04/09 - Apply a new taper that attempts to correct for the incoherent summation, * the decay of the subsource spectra for frequencies less than the subsource * corner frequency (following Frankel, 1995), and the improper long-period spectral * level. * 01/06/09 - Remove "scaling_coefficient" from input and program. * - Put computation of c4taper, fc4taper outside of the sourcespectra function. * - Revise format of input, using a header line before each line containing input. * - Include a low-cut filter. * - Compute average pga and pgv. * 01/10/09 - Increase max number of frequencies to 500, without changing anything else. * 01/13/09 - "apply_taper" was read in from the parameter file, but it has not been used since * I added the Frankel-like scale factor (on 01/04/09). Thus I have formally removed it from this * version of the program. * - Add SD, PSV to FAS, PSA output file and change order of FAS and PSA. * 02/06/09 - Compute and write out husid time series and 95%-5% duration * - Change input to be closer to smsim input * - Read only output file stem name * 02/14/09 - Correct error in setting i0, j0 if n_hypocenters = 1 and nsites > 1 * 02/16/09 - Change input, use stress scaling to adjust Wells and Coppersmith FL, FW * - Compute RJB, RCD, print all distances for each site. * - Include an option to move the sites to positions relative to * the midpoint of the fault or to the tip of the fault * if move_site .eq. .true. and isitecoordflag .eq. 2 (r and epi) * (reset isitecoordflag to 3 in this case) * or if move_site .eq. .true. and isitecoordflag .eq. 3 and one of the * site locations = 0.0 * - Moved writePar before loops, and just after call to getInputParameters. * 02/17/09 - Modify findDistanceAndAzimuth subroutine * - Modify findSubfaultDistance subroutine * - Change FaultStrikeDeg to FaultStrike * - Delete convertDegreeToRadian subroutine * - Delete findDistance subroutine * - Delete shortestSubfaultDistance subroutine * - Read in dl, dw, and compute nl, nw * 02/18/09 - Replace calls to get_time with a call to the system routine date_time * 02/23/09 - Hardwire in choice of 75% rather than 95% in determining the duration. * 75% is recommended by Ou and Herrmann (SRL) and is consistent with a few * Husid plots I made. * 05/28/09 - Specify a hypocenter location by distances along the fault and down the dip. * - Add period column to psa_fs output real wave(:), totalWave(:), vel_total(:), husid(:) allocatable :: wave, totalWave, vel_total, husid dimension slipWeights(200,200), : slat(100),slon(100) : ,dist(100), : accum_fas(500), accum_psa(500), : averagelogPSA(500),averagePSA(500), : averagesqFA(500),averageFA(500),nnFA(500) dimension weightedMoment(200,200) dimension freq(500),psa(500),fa(500),fh(500),hd1(500),hd2(500) real subfaultMoment(200,200), f0(200,200), risetime(200,200) integer NumberOfActiveSubs(200,200) real subfaultDistance(200,200), actualSlip(200,200), : dur_sub(200,200), dur_path(200,200), : t_arrive(200,200), delay(200,200), freq_out(10), : freq4intrp(500), psa4intrp(500) real dur_75_05, arias real alogFASsum(500), alogPSAsum(500) integer nFASsum(500), nPSAsum(500) real avgavgFAS(500), avgavgPSA(500), avgavgPSA_out(10) character f4head(10)*5 real amag, r_ref, rlow(10), a_s(10), b_s(10), m_s(10) real stutter dimension afreq1(500), amp1(500), afreq2(500), amp2(500) dimension inormal(10), xrandom1(10), xrandom2(10), * xrmin(10), xrmax(10) dimension siteLocation(300,2) character fresp1*30,fresp2*30, InputFileName*30 character f_stem*120 character fpar*120, facc*120, fpsa*120, f_h_f0*120, fhusid*120, : f_times*120, f_dur_75_05*120, f_dist_psa*120 logical f_exist, rmv_trend, specify_length, specify_width, : move_site, write_site_files character fault_type*10 character datx*8, time_start*10, time_stop*10 character*30 fparrandom common /par/ iFFTsub, dt, npadl, npadt common /seed/ iseed common /fnames/ fresp1,fresp2 common /params/ : rho,beta,iKapa,fmax, : fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q, : rmin,rd1,rd2,durmin,b1,b2,b3, : iwind,namp1,namp2,totalMoment, : c4taper, fc4taper, : nsubs, flocut, nslope, pi, twopi integer soFarNumberOfSubs integer nptsTotalWave c nptsTotalWave is the total number of points in the final timeseries. nu_ctl = 99 ioPar = nu_ctl - 1 ioPSA = nu_ctl - 2 ioHUS = nu_ctl - 3 ioACC = nu_ctl - 4 ioresp1 = nu_ctl - 5 ioresp2 = nu_ctl - 6 ioa = nu_ctl - 7 nu_h = nu_ctl - 8 nu_i0j0 = nu_ctl - 9 nu_times = nu_ctl - 10 nu_dur_75_05 = nu_ctl - 11 nu_dist_psa = nu_ctl - 20 * Standard Fortran 90 intrinsic Subroutine DATE_AND_TIME call DATE_AND_TIME( datx, time_start ) * Date is returned as 'CCYYMMDD' * Time is returned as 'hhmmss.sss' pi=4.0*atan(1.0) twopi = 2.0 * pi f_exist = .false. do while (.not. f_exist) InputFileName = ' ' write(*, '(a\)') : ' Enter name of input parameter file '// : '(cr = exsim_dmb.params): ' read(*, '(a)') InputFileName if (InputFileName(1:4) .eq. ' ') : InputFileName = 'exsim_dmb.params' call trim_c(InputFileName, nc_f_in) inquire(file=InputFileName(1:nc_f_in), exist=f_exist) if (.not. f_exist) then write(*,'(a)') ' ******* FILE DOES NOT EXIST ******* ' end if end do open(unit=nu_ctl, file=InputFileName(1:nc_f_in), status='unknown') print *, ' Control file: '//InputFileName(1:nc_f_in) c Read (and echo) the input parameters that are common to all grid points. call getInputParameters(nu_ctl, : FaultStrike,FaultDip,h, : FaultLat,FaultLon, move_site, write_site_files, : siteLocation,numberOfSites,isitecoordflag, : fault_type, : FaultLength, FaultWidth, : dl, dw, nl,nw, nsubs, : specify_length, specify_width, stress_ref, : y, : hyp_loc_fl, hyp_loc_fw, i0_in,j0_in,n_hypocenters, : tpadl, tpadt, dt, beta, rho, amag, : stress,iKapa,fmax, : fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q, : r_ref, nsprd_segs, rlow, a_s, b_s, m_s, : rmin,rd1,rd2,durmin,b1,b2,b3, ! trilinear duration parameters : iwind, : flocut, nslope, : nfreq,freq1,freq2, : nfout, freq_out, : f_stem, : namp1,fresp1, namp2,fresp2, : iseed, nsims,damp, : Islipweight, slipWeights, : iRow,jColumn,pulsingPercent, : iPapaFlag,PapaGama,PapaNu,PapaT0,PapaPeak, : iDynamicFlag, : iflagscalefactor, : iflagfas_avg, iflagpsa_avg, : i_rise_time) close(nu_ctl) v=y*beta subFaultRadius=sqrt((dl*dw)/pi) !%%%%%%%%%% riseTime_original=subFaultRadius/v !%%%%%%%%%% call trim_c(f_stem, nc_f_stem) fpar = ' ' fpar = f_stem(1:nc_f_stem)//'_parameters.out' call trim_c(fpar, nc_fpar) open (ioPar,file=fpar(1:nc_fpar),status='unknown') call writePar(ioPar, : FaultStrike,FaultDip,h, : FaultLat, FaultLon, : siteLocation,numberOfSites, move_site, : fault_type, : FaultLength, FaultWidth,nl,nw,dl,dw, : specify_length, specify_width, stress_ref, : y, : hyp_loc_fl, hyp_loc_fw, i0_in,j0_in,n_hypocenters, : iseed, nsims, : amag, : stress, : pulsingPercent,iDynamicFlag, : i_rise_time, : iPapaFlag,PapaGama,PapaNu,PapaT0,PapaPeak, : iflagscalefactor, : iflagfas_avg, iflagpsa_avg, : tpadl, tpadt, : r_ref, nsprd_segs, rlow, a_s, b_s, m_s : ) close(ioPar) * Initialize q_f dummy = q_f_setup(fr1, qr1, s1, ft1, ft2, : fr2, qr2, s2, c_q) * Initialize gsprd_q_f: numsource = 1 dummy = gsprd_f_setup(r_ref,nsprd_segs,rlow, : a_s,b_s,m_s, : numsource,amag) npadl = tpadl/dt npadt = tpadt/dt f_h_f0 = ' ' f_h_f0 = f_stem(1:nc_f_stem)//'_h_f0.out' call trim_c(f_h_f0, nc_f_h_f0) open(nu_h, file=f_h_f0(1:nc_f_h_f0), status='unknown') write(nu_h,'(2x,a, : 1x,a, : 2x, a, 2x,a, : 3x, a, 3x,a, 1x,a, : 1x,a, : 1x,a, : 1x,a, : 5x,a, 1x,a, 9x,a, : 8x,a, : 1x,a, : 3x,a, 4x,a, 7x,a, : 1x,a, 4x,a, : 1x,a, : 2x,a)') : 'isite', : 'ihyp', : 'i0', 'j0', : 'i', 'j', 'nsubs', : 'Diff_T_arrive', : 'NoOfEffectiveSubfaults', : 'NumberOfActiveSubs', : 'f0main', '1stElmntF0', 'f0', : 'hij', : 'iflagscalefactor', : '1/f0main', '1/1stF0', '1/f0', : 'riset_orig', 'dur_sub', : 'subFaultM0', : 'Tot/subM0' f_dur_75_05 = ' ' f_dur_75_05 = f_stem(1:nc_f_stem)//'_dur_75_05.out' call trim_c(f_dur_75_05, nc_f_dur_75_05) open(nu_dur_75_05, file=f_dur_75_05(1:nc_f_dur_75_05), : status='unknown') write(nu_dur_75_05,'(2x,a, : 1x,a, : 2x, a, 2x,a, : 1x,a, : 1x,a)') : 'isite', : 'ihyp', : 'i0', 'j0', : 'dur_75_05', : 'arias(cm/s)' c Generate some logarithmically-spaced frequencies freq(1)=freq1 if (nfreq.ge.2) then frinc=alog(freq2/freq1)/(nfreq-1) do k=2,nfreq freq(k)=freq1*exp((k-1)*frinc) enddo endif c End of Generate some logarithmically-spaced frequencies iseed = -iabs(iseed) if(i0_in.eq.0.or.j0_in.eq.0) : call locateHypocentreRandomly(i0,j0,nl,nw) NoOfEffectiveSubfaults=float(nl)*(pulsingPercent/100.0) NoOfEffectiveSubfaults=NoOfEffectiveSubfaults/2. !!double side propogation if(NoOfEffectiveSubfaults.le.1) NoOfEffectiveSubfaults=1 totalMoment=10.**(1.5*amag+16.05) avgSubFaultMoment=totalMoment/float(nw*nl) firstElementF0= : 4.9e+6*beta*(stress/avgSubFaultMoment)**(1.0/3.0) F0main=4.9e+6*beta*(stress/totalMoment)**(1.0/3.0) call findSubFaultMoment(nl,nw,slipWeights : ,totalMoment,weightedMoment) f_dist_psa = ' ' f_dist_psa = f_stem(1:nc_f_stem)//'_distances_psa.out' call trim_c(f_dist_psa, nc_f_dist_psa) open(nu_dist_psa, file=f_dist_psa(1:nc_f_dist_psa), : status='unknown') f4head = '00000' do i = 1, nfout if (freq_out(i) .lt. 0.0 .or. freq_out(i) .ge. 10.0) then write(f4head(i)(1:5),'(f5.2)') freq_out(i) else write(f4head(i)(2:5),'(f4.2)') freq_out(i) end if end do write(nu_dist_psa,'(1x,a, : 1x,a, 1x,a, 1x,a, : 1x,a, 1x,a, : 5x,a, 3x,a, : 10(7x,a4, 3x,a8))') : 'isite', : 'sitecoord(1)', 'sitecoord(2)', 'isitecoordflag', : 'site_lat_degrees', 'site_lon_degrees', : 'd_jb', 'd_cd2f', : ('freq', 'psa'//f4head(i), i = 1, nfout) DO isite=1,numberOfSites ! Loop on Sites around the Fault SiteLat=siteLocation(isite,1) SiteLon= siteLocation(isite,2) print *, "Working on Site #",isite,SiteLat, SiteLon * Compute some distances: call site_coords_in_degrees( : FaultLat, FaultLon, : SiteLat, SiteLon, isitecoordflag, : site_lat_degrees, site_lon_degrees ) h_min_c = 3.0 ! Campbell depth to seismogenic region w1 = 0.0 w2 = faultwidth s1 = 0.0 s2 = faultlength call dist_3df( : site_lat_degrees, site_lon_degrees, : FaultLat, FaultLon, h, FaultStrike, FaultDip, : w1, w2, s1, s2, : h_min_c, d_jb, az_jb, d_cd2f, az_cd2f, d_c, az_c, : d_sta_n, d_sta_e, irgn_cd2f, irgn_c, irgn_jb) nnFA=0 totalDuration=0 fi2=0 call findDistanceAndAzimuth(FaultLat,FaultLon,SiteLat, : SiteLon,R,fi2,isitecoordflag) c Note that R, azm fi2 are the distance, azm. w.r.t. origin (not epi), c Note also I must input faultstrike in degrees and its converted to rad. alogFASsum = 0.0 alogPSAsum = 0.0 nFASsum = 0 nPSAsum = 0 alogPGAsum = 0.0 alogPGVsum = 0.0 DO ihypo = 1, n_hypocenters if (n_hypocenters .gt. 1) then call locateHypocentreRandomly(i0,j0,nl,nw) else if (i0_in. eq. 0 .or. j0_in .eq. 0) then call locateHypocentreRandomly(i0,j0,nl,nw) else i0 = i0_in j0 = j0_in end if nsubfaults = 0 t_arrive_min=10000. t_end_max=0. risetime_max = 0.0 DO i=1,nl DO j=1,nw nsubfaults = nsubfaults + 1 subFaultMoment(i,j)=weightedMoment(i,j) if(iDynamicFlag.eq.0) then NumberOfActiveSubs(i,j)=1 else NumberOfActiveSubs(i,j)=NumberOfPulsingSubs(i0,j0,i,j : ,nl,nw,NoOfEffectiveSubfaults) if(NumberOfActiveSubs(i,j).eq.0) then NumberOfActiveSubs(i,j)=1 end if end if f0(i,j)= : firstElementF0*(NumberOfActiveSubs(i,j)**(-1.0/3.0)) if (i_rise_time .eq. 1) then risetime(i,j) = riseTime_original else if (i_rise_time .eq. 2) then riseTime(i,j) = 1.0/f0(i,j) else print *,' ERROR: invalid value of i_rise_time, = ', : i_rise_time stop end if if(risetime(i,j) .gt. risetime_max) then risetime_max = risetime(i,j) end if actualSlip(i,j)=1.e-22*subFaultMoment(i,j)/ : (rho*beta**2.*dl*dw) subfaultDistance(i,j)= : findSubfaultDistance(R,h,FaultStrike, : fi2,FaultDip,dl,dw,i,j) c calculate duration of subsource time history at distance c subfaultDistance call dur_path_cmp(subfaultDistance(i,j),dur_path(i,j)) ! computes path part of duration dur_sub(i,j) = dur_path(i,j) + risetime(i,j) if ((i-i0) .ne. 0 .or. (j-j0) .ne. 0) then delay(i,j)=sqrt((dl*(i-i0))**2.+(dw*(j-j0))**2.)/v else delay(i,j)=0. ! delay from hypocenter to itself end if t_arrive(i,j) = delay(i,j) + subfaultDistance(i,j)/beta if (t_arrive(i,j) .lt. t_arrive_min) : t_arrive_min = t_arrive(i,j) t_end= t_arrive(i,j) + dur_sub(i,j) if (t_end .gt. t_end_max) then imax = i jmax = j t_end_max = t_end end if END DO ! loop over nw to calculate quantities not dependent on nsims END DO ! loop over nl to calculate quantities not dependent on nsims nptsTotalWave=(t_end_max - t_arrive_min + tpadl + tpadt + : risetime_max)/dt ! need risetime_max because the stutter ! delay can shift the subsource time series ! corresponding to t_end_max by an amount up to ! risetime of that subsource; to be safe I ! use risetimemax rather than risetime(imax,jmax), ! where imax, jmax are the indices of the subsource ! corresponding to t_end_max. signnpw2 = +1.0 call get_npw2(nptsTotalWave,signnpw2,npw2TotalWave) if (write_site_files) then if (ihypo .eq. 1) then f_times = ' ' f_times = f_stem(1:nc_f_stem)//'_times_s000.out' write(f_times(nc_f_stem+9:nc_f_stem+11),'(i3.3)') isite call trim_c(f_times, nc_f_times) open (nu_times,file=f_times(1:nc_f_times), : status='unknown') write(nu_times,'(1x,a, 2x,a, 2x,a, 3x,a, 3x,a, : 1x,a, 1x,a, : 1x,a, 1x,a, 1x,a, : 1x,a, 1x,a, : 1x,a, : 1x,a, 1x,a, 1x,a, : 1x,a, : 1x,a, : 1x,a, 1x,a)') : 'itr', 'i0', 'j0', 'i', 'j', : 't_arrive(i0,j0)', 't_arrive_min', : 't_arrive(i,j)', 'risetime(i,j)', 'risetime_max', : 'dur_path(i,j)', 'dur_sub(i,j)', : 't_arrive(i,j)+dur_sub(i,j)', : 'imx', 'jmx', 't_arrive(imax,jmax)', : 'dur_sub(imax,jmax)', : 't_arrive(imax,jmax)+dur_sub(imax,jmax)', : 'risetime_max', 't_end_max' do i = 1, nl do j = 1, nw write(nu_times,'(5(1x,i3), : 8x,f8.1, 5x,f8.1, : 6x,f8.1, 6x,f8.1, 5x,f8.1,, : 5x,f9.2, 4x,f9.2, : 19x,f8.1, : 2(1x,i3), 12x,f8.1, : 11x,f8.1, : 31x,f8.1, : 5x,f8.1, 2x,f8.1)') : ihypo, i0, j0, i, j, : t_arrive(i0, j0), t_arrive_min, : t_arrive(i,j), risetime(i,j), risetime_max, : dur_path(i,j), dur_sub(i,j), : t_arrive(i,j)+dur_sub(i,j), : imax, jmax, t_arrive(imax,jmax), : dur_sub(imax,jmax), : t_arrive(imax,jmax)+dur_sub(imax,jmax), : risetime_max, t_end_max end do end do close(nu_times) end if end if HypoDistance=findSubfaultDistance(R,h,FaultStrike, : fi2,FaultDip,dl,dw,i0,j0) c now do actual generation and summation of subfault time series *** ! Initialize arrays (index over frequency) containing cumulative sums for averages averagelogPSA=0.0 averagePSA=0.0 averageFA=0.0 accum_fas = 0.0 accum_psa = 0.0 nnFA = 0 accum_logpga = 0.0 accum_logpgv = 0.0 accum_dur_75_05 = 0.0 accum_arias = 0.0 DO isim=1,nsims print *, ' Site isite, Hypocenter iteration ihypo, '// : 'Simulation number isim = ', isite, ihypo, isim allocate (totalWave(npw2TotalWave)) allocate (husid(npw2TotalWave)) totalWave = 0.0 DO i=1,nl DO j=1,nw nsubsource = npadl + dur_sub(i,j)/dt + npadt signnpw2 = +1.0 call get_npw2(nsubsource,signnpw2,iFFTsub) allocate (wave(iFFTsub)) wave=0 call computeStochasticWave( : wave, : subFaultMoment(i,j), subfaultDistance(i,j), : F0Main, f0(i,j), risetime(i,j), dur_sub(i,j), : nl,nw, : iflagscalefactor, hij) ! dmb ! NOTE: The length of wave is iFFTsub; this is passed through common /par/ if (isim .eq. 1) then write(nu_h,'(4x,i3, : 2x,i3, : 1x,i3, 1x,i3, : 1x,i3, 1x,i3, 1x,i5, : 1p, 4x,e10.3, 0p, : 18x,i5, : 14x,i5, : 1p, 4(1x,e10.3), : 0p, 16x,i1, : 1p, 7(1x,e10.3))') : isite, : ihypo, : i0, j0, : i, j, nsubfaults, : t_arrive(i,j)-t_arrive(i0,j0), : NoOfEffectiveSubfaults, : NumberOfActiveSubs(i,j), : f0main, firstElementF0, f0(i,j), : hij, : iflagscalefactor, : 1.0/f0main, 1.0/firstElementF0, 1.0/f0(i,j), : risetime_original, dur_sub(i,j), : subFaultMoment(i,j), : TotalMoment/subFaultMoment(i,j) end if c subsource accelerogram is randomly delayed to c simulate complexity in slip process c subtract minimum delay to eliminate static time shift stutter = Ran1(iseed) * riseTime(i,j) ishift = ((t_arrive(i,j) - t_arrive_min) + stutter)/dt if (npadl+ishift+dur_sub(i,j)+npadt .gt. : npw2TotalWave) then print *,'npadl+ishift+dur_sub(i,j)+'// : 'npadt.gt.npw2TotalWave' print *, npadl+ishift+dur_sub(i,j)+npadt, : npw2TotalWave end if ! No need to call shift---just figure out the shifted index and use it ! to fill the totalWave vector. c add current accelerogram to total wave field nhigh = min(iFFTsub, npw2TotalWave - ishift) DO k = 1, nhigh totalWave(k+ishift)=totalWave(k+ishift)+wave(k) END DO deallocate (wave) END DO ! loop over nw END DO ! loop over nl cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c include Mavroeidis and Papageorgiou, 2003 approach if( int(iPapaFlag).eq.1)then call MavroPapa(totalWave,npw2TotalWave, * iPapaFlag,PapaGama,PapaNu,PapaT0,amag,PapaPeak) endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c ************************** end of summation ********************** c accelerogram and model parameters are saved only once * DEBUG print *,' iseed = ', iseed * DEBUG pga=findPeakValue(npw2TotalWave,totalWave) accum_logpga = accum_logpga + alog10(pga) allocate (vel_total(npw2TotalWave)) vel_total = 0.0 call acc2v(totalWave, npw2TotalWave, dt, .false., vel_total) pgv=findPeakValue(npw2TotalWave,vel_total) accum_logpgv = accum_logpgv + alog10(pgv) deallocate(vel_total) rmv_trend = .false. husid = 0.0 call accsqint(totalWave, npw2TotalWave, dt, rmv_trend, : husid) husid_end = husid(npw2TotalWave) do i = 1, npw2TotalWave husid(i) = husid(i)/husid_end end do call locate(husid,npw2TotalWave,0.05,j05) call locate(husid,npw2TotalWave,0.75,j75) accum_dur_75_05 = accum_dur_75_05 + : alog10(float(j75-j05)*dt) ! for geometric mean accum_arias = accum_arias + : alog10((pi/(2.0*981.0))*husid_end) ! for geometric mean if (write_site_files) then if (ihypo.eq.1 .and. isim .eq. 1) then facc = ' ' facc = f_stem(1:nc_f_stem)//'_acc_s000.out' write(facc(nc_f_stem+7:nc_f_stem+9),'(i3.3)') isite call trim_c(facc, nc_facc) open (ioAcc,file=facc(1:nc_facc),status='unknown') call writeAcc(ioAcc, totalWave, npw2TotalWave, pga, : isim,isite,HypoDistance,d_cd2f, fpar) close(ioAcc) fhusid = ' ' fhusid = f_stem(1:nc_f_stem)//'_husid_s000.out' write(fhusid(nc_f_stem+9:nc_f_stem+11),'(i3.3)') isite call trim_c(fhusid, nc_fhusid) open (ioHus,file=fhusid(1:nc_fhusid),status='unknown') dur_75_05_sim1 = float(j75-j05)*dt arias_sim1 = (pi/(2.0*981.0))*husid_end call writeHUS(ioHus, husid, npw2TotalWave, : arias_sim1, dur_75_05_sim1, : isim,isite,HypoDistance,d_cd2f,fpar) close(ioHUS) end if end if dur=float(npw2TotalWave-1)*dt-5.0 psa=0.0 call computeResponse(totalWave,dt,dur,damp,nfreq,freq,psa) call computeAverageResponse(accum_psa,psa,nfreq,nsims, : iflagpsa_avg) ! dmb FA=0.0 call computeFACCN(totalWave,npw2TotalWave,dt,nfreq,freq,FA) call computeAverageFA(accum_fas,FA,nfreq,nsims,nnFA, : iflagfas_avg) write (*,"(' *** finished with trial',i4)") isim deallocate(totalWave) deallocate(husid) END DO ! End loop over nsims averagePGA = 10.0**(accum_logpga/float(nsims)) alogPGAsum = alogPGAsum + alog10(averagePGA) averagePGV = 10.0**(accum_logpgv/float(nsims)) alogPGVsum = alogPGVsum + alog10(averagePGV) DO i = 1, nfreq if (iflagpsa_avg .eq. 2 .and. accum_psa(i) .eq. 0.0) then averagePSA(i) = -9.99 else if (iflagpsa_avg .eq. 1) then ! arithmetic averagePSA(i) = accum_psa(i)/float(nsims) else if (iflagpsa_avg .eq. 2) then ! geometric averagePSA(i) = 10.0**(accum_psa(i)/float(nsims)) else if (iflagpsa_avg .eq. 3) then ! rms averagePSA(i) = sqrt(accum_psa(i)/float(nsims)) else print *, : ' ERROR, iflagpsa_avg not 1, 2, or 3; = ', iflagpsa_avg stop end if end if if (iflagfas_avg .eq. 2 .and. accum_fas(i) .eq. 0.0) then averageFA(i) = -9.99 else if (iflagfas_avg .eq. 1) then ! arithmetic averageFA(i) = accum_fas(i)/float(nnFA(i)) else if (iflagfas_avg .eq. 2) then ! geometric averageFA(i) = 10.0**(accum_fas(i)/float(nnFA(i))) else if (iflagfas_avg .eq. 3) then ! rms averageFA(i) = sqrt(accum_fas(i)) else print *, : ' ERROR, iflagfas_avg not 1, 2, or 3; = ', iflagfas_avg stop end if end if if (averagePSA(i) .gt. 0.0) then alogPSAsum(i) = alogPSAsum(i) + alog10(averagePSA(i)) nPSAsum(i) = nPSAsum(i) + 1 end if if (averageFA(i) .gt. 0.0) then alogFASsum(i) = alogFASsum(i) + alog10(averageFA(i)) nFASsum(i) = nFASsum(i) + 1 end if END DO ! loop over nfreq dur_75_05 = 10.0**(accum_dur_75_05/float(nsims)) arias = 10.0**(accum_arias/float(nsims)) write(nu_dur_75_05,'(4x,i3, : 2x,i3, : 1x,i3, 1x,i3, : 1x,f9.2, : 1p, 1x,e10.3)') : isite, : ihypo, : i0, j0, : dur_75_05, : arias END DO !loop over random hypocenters * Compute average of average FAS and PSA avgavgPGA = 10.0**(alogPGAsum/float(n_hypocenters)) avgavgPGV = 10.0**(alogPGVsum/float(n_hypocenters)) DO i = 1, nfreq avgavgPSA(i) = 10.0**(alogPSAsum(i)/float(nPSAsum(i))) avgavgFAS(i) = 10.0**(alogFASsum(i)/float(nFASsum(i))) END DO ! end loop over nfreq nfreq4intrp = nfreq do i = 1, nfreq4intrp freq4intrp(i) = freq(i) psa4intrp(i) = avgavgPSA(i) end do do i = 1, nfout if (freq_out(i) .eq. -1.0) then ! pgv avgavgPSA_out(i) = avgavgPGV else if (freq_out(i) .ge. 99.0) then ! pga avgavgPSA_out(i) = avgavgPGA else avgavgPSA_out(i) = : yintrf( freq_out(i), freq4intrp, psa4intrp, nfreq4intrp) end if end do if (write_site_files) then fpsa = ' ' fpsa = f_stem(1:nc_f_stem)//'_psa_fs_s000.out' write(fpsa(nc_f_stem+10:nc_f_stem+12),'(i3.3)') isite call trim_c(fpsa, nc_fpsa) open (ioPsa,file=fpsa(1:nc_fpsa),status='unknown') call writePSA_FA(ioPSA, freq,avgavgPGA,avgavgPGV, : avgavgPSA,nfreq,nsims,damp,amag, : siteLat,siteLon,h,r,fi2*180./pi,d_cd2f,avgavgFAS, : isite, HypoDistance,fPar, nl, nw) close(ioPSA) end if write(nu_dist_psa, '( : 1x,i5, : 4x,f9.3, 4x,f9.3, 14x,i1, : 8x,f9.3, 8x,f9.3, : 1x,f8.2, 1x,f8.2, : 1p,10(1x,e10.3, 1x,e10.3))') : isite, : SiteLat, SiteLon, isitecoordflag, : site_lat_degrees, site_lon_degrees, : d_jb, d_cd2f, : (freq_out(i), avgavgPSA_out(i), i= 1, nfout) END DO ! End loop over sites print * * Standard Fortran 90 intrinsic Subroutine DATE_AND_TIME call DATE_AND_TIME( datx, time_stop ) * Date is returned as 'CCYYMMDD' * Time is returned as 'hhmmss.sss' call time_diff(time_start, time_stop, time_elapsed) print *, : ' Elapsed time (sec): ', time_elapsed close(nu_h) close(nu_dur_75_05) close(nu_dist_psa) STOP END cccccccccccccccccccccccc End of Main Program cccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !------------------- SUBPROGRAMS ----------------------------------------- cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real function amplf (f,afreq,amp,namp) c computes amplification value at frequency f by linear c interpolation between its adjoining values dimension amp(*),afreq(*) if (f.le.afreq(1)) then amplf=amp(1) return endif if (f.gt.afreq(namp)) then amplf=amp(namp) return endif do i=2,namp if (f.le.afreq(i)) then d=(amp(i)-amp(i-1))/(afreq(i)-afreq(i-1))*(f-afreq(i-1)) amplf=amp(i-1)+d return endif enddo end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine avg(n,s,avamp2) c calculates average squared amplitude spectrum complex s dimension s(*) sum=0. do j=1,n amp=cabs(s(j)) sum=sum+amp*amp enddo avamp2=sum/float(n) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine computeAverageFA(accum_fas,FA,nfreq,nsims,nn, : iflagfas_avg) c calculate average Fourier spectrum ! revised by dmb dimension accum_fas(*),FA(*), nn(500) do i=1,nfreq if(FA(i).gt.0)then nn(i)=nn(i)+1 if (iflagfas_avg .eq. 1) then ! arithmetic accum_fas(i)= : accum_fas(i) + FA(i) else if (iflagfas_avg .eq. 2) then ! geometric accum_fas(i)= : accum_fas(i) + alog10(FA(i)) else if (iflagfas_avg .eq. 3) then ! rms accum_fas(i)= : (accum_fas(i)*float(nn(i)-1) + FA(i)**2.0 )/float(nn(i)) ! else if (iflagfas_avg .eq. 3) then ! rms ! accum_fas(i)= ! : accum_fas(i) + FA(i)**2.0 else print *, : ' ERROR, iflagfas_avg not 1, 2, or 3; = ', iflagfas_avg stop end if endif enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine computeAverageResponse(accum_psa,psa,nfreq,nsims, : iflagpsa_avg) c calculate average response spectrum ! revised by dmb dimension accum_psa(*),psa(*) do i=1,nfreq if (iflagpsa_avg .eq. 1) then ! arithmetic accum_psa(i)= : accum_psa(i) + psa(i) else if (iflagpsa_avg .eq. 2) then ! geometric accum_psa(i)= : accum_psa(i) + alog10(psa(i)) else if (iflagpsa_avg .eq. 3) then ! rms accum_psa(i)= : accum_psa(i) + psa(i)**2.0 else print *, : ' ERROR, iflagpsa_avg not 1, 2, or 3; = ', iflagpsa_avg stop end if enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine computeFACCN(totalWave,npts,dt,nfreq,freq, FA) c c Compute Fourier acceleration spectrum for final timeseries. dimension totalWave(*), freq(500), FA(500) complex fft(:) real spectrum(:) allocatable :: fft, spectrum allocate(fft(npts), spectrum(npts)) fft = 0.0 c Copy totalWave into new array because we need to save original for later. do j=1, npts spectrum(j)=totalWave(j) enddo call makeItComplex(spectrum,fft,npts) c Take padded time series into freq. domain call fork(npts,fft,-1.) ncall = npts/2 +1 df = 1./(dt * float(npts)) c Scale spectrum and put value back into spectrum array. fft(ncall) = cmplx(real(fft(ncall)),0.) do jf = 1, ncall spectrum(jf) = dt* sqrt(float(npts)) * cabs(fft(jf)) enddo c Now sample the spectrum into the same freq. bins used for PSA output. call sample(spectrum, ncall, df, nfreq, freq, FA) c Array FA now contains the nfreq values of sampled faccn. c TotalWave input array is unchanged. deallocate(fft, spectrum) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine computeResponse(a,dt,dur,damp, Nfreq,freq,psa) c calculates response spectrum from acceleration time history * Dates: 11/23/08 - Some changes by DMB (compute pi, use rdcalcdp) dimension Freq(*),psa(*),a(*) real pi,freq1,freq2,damp,Maxtim c damp is % critical damping. freq2 is maximum freq. to consider. pi = 4.0*atan(1.0) c Convert % damping to fraction beta=damp*0.01 c Add 5 sec. to dur for shake-down maxtim=dur+5.0 n=ifix(dur/dt) ntime=ifix(maxtim/dt) do j=n+1,ntime a(j)=0. enddo c Call new response spectrum routine for each desired freq do k=1,nfreq omega=2.*pi*freq(k) call rdcalcdp(a, ntime, omega, beta, dt, sd) PSA(k)=sd*omega*omega enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine computeStochasticWave( : seism, : subFaultMoment, subfaultDistance, : F0main, f0, risetime, dur_sub, : nl,nw, : iflagscalefactor, scalingFactor) c calculates accelerogram from individual subfault by using c band-limited white noise method (Boore, BSSA, no.6/83). c An w-squared spectrum is assumed * Dates: 08/20/08 - Added scale factor flag so can choose which factor * to use without recompiling the program * (iflagscalefactor = 1, 2, 3, 4 = vel^2, acc^2, acc, DMB respectively). * 08/21/08 - Add scalingfactor (h) to output. * Do/do not apply the scalefactor taper, depending on the * value of taper_scalefactor. * Doubled dimension of s * 12/01/08 - Removed "NumberOfActiveSubs" from argument lists * 12/17/08 - Added f0 and risetime to argument list and removed iii, jjj, z. * Also, the program used dur rather than dur_sub in ndur, etc. * I changed this. I also pass dur_sub through the argument list * rather than compute it here. dimension seism(*) complex s(:) allocatable :: s common /par/ iFFTsub, dt, npadl, npadt common /params/ : rho,beta,iKapa,fmax, : fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q, : rmin,rd1,rd2,durmin,b1,b2,b3, : iwind,namp1,namp2,totalMoment, : c4taper, fc4taper, : nsubs, flocut, nslope, pi, twopi common /seed/ iseed ! Dave added this statement c specify point-source constants tmax=float(iFFTsub)*dt df=1.0/tmax taper=0.02 prtitn=1.0/sqrt(2.0) rtp=0.55 fs=2.0 c define spectral constant, for subfaultDistance=1 km const=prtitn*rtp*fs*(1.e-20)/(4.*pi*rho*beta**3.) nnyq=iFFTsub/2+1 if (dur_sub .lt. 0.0) then print *,'STOPPING!!! Negative duration. Check duration model' stop end if ndur=dur_sub/dt ntaper=int(taper*ndur) nstop=ndur+ntaper+ntaper if(nstop.gt.iFFTsub) then write (*,*) 'PROGRAM STOPS: subsource duration too long.' write (*,*) 'Options to fix the problem:' write (*,*) '1) Increase iFFTsub or dt, or' write (*,*) '2) Reduce distance to observation point, or' write (*,*) '3) Change parameters of duration model' stop endif allocate (s(iFFTsub)) ! iFFTsub is passed through the common block "par" do i=1,iFFTsub seism(i)=0.0 s(i)=cmplx(0.0,0.0) seism(i)=0.0 enddo c generate Gaussian white noise with zero mean, unit variance. Noise is c tapered. do i=1,nstop if (iwind .eq. 0) then call window(i,1,nstop,ntaper,wind) else if (iwind .eq. 1) then call wind2 (i,1,nstop,ntaper,dur_sub,0.2,0.2,wind) else print *,'STOPPING: invalid value of iwind (= ', iwind, ')' stop end if s(i+npadl)=wind*cmplx(ggnqf(iseed),0.) enddo c transform to frequency domain call fork (iFFTsub,s,-1.) c find scaling factor, H averageMoment=totalMoment/float(nsubs) if (iflagscalefactor .eq. 3) then ScalingFactor=sqrt(float(nsubs))*(F0main/f0)**2 !! David Boore else sum1=0 sum2=0 do i=1,nnyq f=(i-1)*df if (iKapa.eq.0) highc=1./sqrt(1.+(f/fmax)**8.) if (iKapa.eq.1) highc=exp(-pi*fmax*f) if (iflagscalefactor .eq. 1) then ! vel spectrum s1=const*totalMoment* : (2.*pi*f)**1.0*(1/(1.+(f/F0main)**2.))* : highc s2=const*averageMoment* : (2.*pi*f)**1.0*(1/(1.+(f/f0)**2.))* : highc else if (iflagscalefactor .eq. 2) then ! acc spectrum s1=const*totalMoment* : (2.*pi*f)**2.0*(1/(1.+(f/F0main)**2.))* : highc s2=const*averageMoment* : (2.*pi*f)**2.0*(1/(1.+(f/f0)**2.))* : highc else print *,' ERROR, iflagscalefactor = ', iflagscalefactor print *, '; not a legal value; QUITTING!!!' stop end if sum1=sum1+s1**2/float(nsubs) sum2=sum2+s2**2 enddo scalingFactor=sqrt(sum1/sum2) end if * Compute parameters used to taper spectrum toward low frequencies c4taper = sqrt(float(nsubs))/scalingFactor fc4taper = f0/sqrt(c4taper) c multiply spectrum by an w-squared spectral shape after normalizing to c square of unit spectral amplitude call avg (nnyq,s,avamp2) avamp=sqrt(avamp2) do i=1,nnyq f=(i-1)*df s(i)=sourceSpectra(f,subfaultDistance,const, : scalingFactor,subFaultMoment,f0) * : s(i)/avamp if (i.ne.1) s(iFFTsub+2-i)=conjg(s(i)) enddo s(nnyq)=cmplx(real(s(nnyq)),0.) c transform back to time domain call fork(iFFTsub,s,1.) afact=sqrt(float(iFFTsub))/tmax do i=1,iFFTsub seism(i)=afact*real(s(i)) enddo deallocate (s) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine createRandomWeight(nl,nw,slipWeights) common /seed/ iseed dimension slipWeights(200,200) do i=1,nl do j=1,nw c** draw random number 0 to 1. This was changed from c using =ggnqf(iseed)+1. Also impose nonzero weight. slipWeights(i,j)=Ran1(iseed) if (slipWeights(i,j).le.0.) slipWeights(i,j)=0.001 enddo enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine dur_path_cmp (subfaultDistance, dur_path) c computes duration dur due to path for a given subsource common /params/ : rho,beta,iKapa,fmax, : fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q, : rmin,rd1,rd2,durmin,b1,b2,b3, : iwind,namp1,namp2,totalMoment, : c4taper, fc4taper, : nsubs, flocut, nslope, pi, twopi if (subfaultDistance.le.rmin) : dur_path= durmin if (subfaultDistance.gt.rmin.and.subfaultDistance.le.rd1) : dur_path= durmin+b1*(subfaultDistance-rmin) if (subfaultDistance.gt.rd1.and.subfaultDistance.le.rd2) : dur_path= durmin+b1*(rd1-rmin)+b2*(subfaultDistance-rd1) if (subfaultDistance.gt.rd2) : dur_path= durmin+b1*(rd1-rmin)+b2*(rd2-rd1)+ : b3*(subfaultDistance-rd2) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine findDistanceAndAzimuth (FaultLat,FaultLon,SiteLat, * SiteLon,epi,azi,isitecoordflag) c calculates distance and azimuth between two c points using their latitudes and longitudes * Dates: 02/17/09 - renamed some variables to avoid possible confusion * between angles in degrees and radians. * The input and output angles are in degrees pi = 4.0*atan(1.0) d2r = pi/180.0 if (isitecoordflag .eq. 1) then ! lat,long re=6371. alat=(FaultLat+SiteLat)/2. alat_radians=d2r*alat dlat=SiteLat-FaultLat ! degrees dlon=SiteLon-FaultLon ! degrees r1 = d2r*re*(SiteLat-FaultLat) epi= d2r*re*sqrt(dlat**2.+cos(alat_radians)**2.*dlon**2.) if(r1.ge.epi)then azi_radians=pi else azi_radians=acos(r1/epi) endif if (dlon.le.0.) azi_radians=2.0*pi-azi_radians azi= azi_radians/d2r else if (isitecoordflag .eq. 2) then ! R,Az epi = sitelat azi = sitelon else ! assume cartesian coordinates xn = sitelat xe = sitelon epi = sqrt(xn**2 + xe**2) azi = 180.0*atan2(xe, xn)/pi if (azi .lt. 0.0) then azi = 360.0 + azi end if end if return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real function findPeakValue(nptsTotalWave,totalWave) c find absolute maximum acceleration in simulated time history dimension totalWave(*) peakValue=0. do i=1,nptsTotalWave if(abs(totalWave(i)).gt.peakValue) peakValue=abs(totalWave(i)) enddo findPeakValue=peakValue return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real function findSubfaultDistance (R,h,FaultStrike, * fi2,FaultDip,dl,dw,i,j) c computes distance subfaultDistance from center of subfault (i,j) c to observation point using formula (1) (Figure 1) * Note confusion in what dip means in the figure 1 of Beresnev and Atkinson. * What they label as "delta" is what was used in the original version of * this subroutine as "dip". Their figure has delta_1 as the true * fault dip, probably added as a result of a review. I have * replaced "dip" by 90_faultdip here to eliminate confusion (D. Boore, 17Feb09). pi = 4.0*atan(1.0) d2r = pi/180.0 a90_faultdip_radians = d2r*(90.0-FaultDip) phi2_strike_radians = d2r*(fi2-FaultStrike) t1=R*cos(phi2_strike_radians)-(2.*i-1)*dl/2. t2=R*sin(phi2_strike_radians)-(2.*j-1)*dw/2.* : sin(a90_faultdip_radians) t3=-h-(2.*j-1)*dw/2.*cos(a90_faultdip_radians) findSubfaultDistance=sqrt(t1**2.+t2**2.+t3**2.) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine findSubFaultMoment(nl,nw,slipWeights,totalMoment, * weightedMoment) dimension slipWeights(200,200) dimension weightedMoment(200,200) c calculate moments of subfaults. First calculate the total of all weights c (totalSlipWeights), then moment per unit weight (elem), then individual moments and c put them back into array slipWeights(i,j) c Seismic Moment = stress para* (L**3) totalSlipWeights=0. do i=1,nl do j=1,nw totalSlipWeights=totalSlipWeights+slipWeights(i,j) enddo enddo elem=totalMoment/totalSlipWeights do i=1,nl do j=1,nw weightedMoment(i,j)=slipWeights(i,j)*elem enddo enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine findTotalDuration(nl,nw,dl,dw,FaultStrike,fi2, : beta,i0,j0,dip,R,h,v,t_arrive_min, : risetime, subfaultdistance, dur_path, : totalDuration) c calculate miminum time delay (t_arrive_min) for subfault radiation to reach c observation point. Calculate estimated total duration (t_end_max-t_arrive_min) of finite- c fault accelerogram. Numbers i,j define current subfault (Figure 1) duration. * The original version added one second to the duration "for sufficiency", but it * also had the motion starting 2 sec before the actual first arrival, thus * adding 2 sec to the duration. I've kept these here, but in a future version * I will use tpads given in the input control file rather than the values of 2 and 1 sec. * Dates: 12/17/08 - Since compute risetime, dur_path, subfault distance before calling this, * pass, those arrays through the argument list rather than * doing the calculations again. * I made extensive changes to this subroutine, including * renaming variables ("delay" was used for the start time of * a subsource as well as the time of arrival at the station). real risetime(200,200), subfaultdistance(200,200), : dur_path(200,200) t_arrive_min=10000. t_end_max=0. DO i=1,nl DO j=1,nw if ((i-i0) .ne. 0 .or. (j-j0) .ne. 0) then delay=sqrt((dl*(i-i0))**2.+(dw*(j-j0))**2.)/v else delay=0. ! delay from hypocenter to itself end if t_arrive = delay + subfaultDistance(i,j)/beta if (t_arrive .lt. t_arrive_min) t_arrive_min = t_arrive t_end= t_arrive + dur_path(i,j) + risetime(i,j) if (t_end .gt. t_end_max) t_end_max = t_end END DO END DO totalDuration=t_end_max - (t_arrive_min - 2.0) + 1.0 c** the total duration includes delay times of subfaults and 2sec timeshift. return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine getInputParameters(nu_ctl, : FaultStrike,FaultDip,h, : FaultLat,FaultLon, move_site, write_site_files, : siteLocation,numberOfSites,isitecoordflag, : fault_type, : FaultLength, FaultWidth, : dl, dw, nl, nw, nsubs, : specify_length, specify_width, stress_ref, : y, : hyp_loc_fl, hyp_loc_fw, i0_in,j0_in,n_hypocenters, : tpadl, tpadt, dt, beta, rho, amag, : stress,iKapa,fmax, : fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q, : r_ref, nsprd_segs, rlow, a_s, b_s, m_s, : rmin,rd1,rd2,durmin,b1,b2,b3, ! trilinear duration parameters : iwind, : flocut, nslope, : nfreq,freq1,freq2, : nfout, freq_out, : f_stem, : namp1,fresp1, namp2,fresp2, : iseed, nsims,damp, : Islipweight, slipWeights, : iRow,jColumn,pulsingPercent, : iPapaFlag,PapaGama,PapaNu,PapaT0,PapaPeak, : iDynamicFlag, : iflagscalefactor, : iflagfas_avg, iflagpsa_avg, : i_rise_time) * Dates: 11/30/08 - Move computation of random slipweights * if Islipweight .eq. 1.0 from main program to here * (a more logical place). * 12/02/08 - Get i_rise_time, to determine what type of risetime is used real rlow(*), a_s(*), b_s(*), m_s(*), freq_out(*) dimension slipWeights(200,200), siteLocation(300,2) character f_stem*(*), : fresp1*30,fresp2*30,aline*60, fault_type*(*), : version_ctl*8, version_in*30 character cmnts2skip(50)*80, buf_in*10 logical f_exist, specify_length, specify_width, : move_site, write_site_files pi = 4.0*atan(1.0) d2r = pi/180.0 version_ctl = ' ' version_ctl = '05/28/09' call trim_c(version_ctl,nc_version_ctl) call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) version_in = ' ' read(nu_ctl,'(a)') version_in call trim_c(version_in,nc_version_in) if (version_ctl(1:nc_version_ctl) .ne. : version_in(1:nc_version_in)) then write(*,'(a)') : ' The control file has the wrong version number; STOP!' close(nu_ctl) stop end if call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) aline = ' ' read (nu_ctl,'(a)')aline call trim_c(aline, nc_aline) ! print *, ' Title: '//aline(1:nc_aline) call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) buf_in = ' ' read (nu_ctl,'(a)') buf_in call upstr(buf_in) call trim_c(buf_in, nc_buf_in) if (buf_in(1:1) .eq. 'Y') then write_site_files = .true. else write_site_files = .false. end if call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) amag,stress,iKapa,fmax ! if( iKapa.eq.1)then ! print *,' mag,stress,kappa = ', amag,stress,fmax ! else ! print *,' mag,stress,fmax = ', amag,stress,fmax ! endif call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) FaultLat, FaultLon ! print *,' FaultLat,FaultLon = ', FaultLat, FaultLon call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) FaultStrike, FaultDip, h ! print *,' FaultStrike,FaultDip(degree),Depth(km) = ', ! : FaultStrike,FaultDip,h call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) fault_type = ' ' read(nu_ctl,'(a)') fault_type call trim_c(fault_type, nc_fault_type) call upstr(fault_type(1:1)) call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) FaultLength, FaultWidth, dl, dw, stress_ref stress_factor = (stress_ref/stress)**(1.0/3.0) specify_length = .true. specify_width = .true. if (faultlength .eq. 0.0) then specify_length = .false. if(fault_type(1:1) .eq. 'S') then faultlength = 10.0**(-2.57+0.62*amag) else if(fault_type(1:1) .eq. 'R') then faultlength = 10.0**(-2.42+0.58*amag) else if(fault_type(1:1) .eq. 'N') then faultlength = 10.0**(-1.88+0.50*amag) else faultlength = 10.0**(-2.44+0.59*amag) end if ! print *, ' WC FL = ', faultlength faultlength = stress_factor * faultlength ! print *, ' WC FL, after scaling = ', faultlength end if if (faultwidth .eq. 0.0) then specify_width = .false. if(fault_type(1:1) .eq. 'S') then faultwidth = 10.0**(-0.76+0.27*amag) else if(fault_type(1:1) .eq. 'R') then faultwidth = 10.0**(-1.61+0.41*amag) else if(fault_type(1:1) .eq. 'N') then faultwidth = 10.0**(-1.14+0.35*amag) else faultwidth = 10.0**(-1.01+0.32*amag) end if ! print *, ' WC FW = ', faultwidth faultwidth = stress_factor * faultwidth ! print *, ' WC FW, after scaling = ', faultwidth end if nl=FaultLength/dl if (nl .lt. 1) nl = 1 nw=FaultWidth/dw if (nw .lt. 1) nw = 1 nsubs = nl * nw dl = FaultLength/float(nl) ! Need to reset dl and dw because nl and nw are integers dw = FaultWidth/float(nw) ! write(*,'(a,4f9.3)')'FaultLength,FaultWidth (km) :' ! : , FaultLength,FaultWidth ! write(*,'(a,2i8)') ! : 'No. of subfaults along length and width of fault plane:', ! : nl,nw if(nl.gt.200.or.nw.gt.200) then print *,"Burp!, you exceeded the : the array dimentions, (200,200)" stop end if call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) y ! print *,' y (=vrup/beta) = ', y call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) hyp_loc_fl, hyp_loc_fw, n_hypocenters if (nl .eq. 1 .and. nw .eq. 1) then ! a single subsource i0_in = 1 j0_in = 1 n_hypocenters = 1 else if (hyp_loc_fl .gt. 0.0 .and. hyp_loc_fw .ge. 0.0 ) then ! nonrandom hypocenter ! convert hyp_loc to subfault index i0_in = int(hyp_loc_fl/dl)+1 j0_in = int(hyp_loc_fw/dw)+1 n_hypocenters = 1 ! not a random hypocenter, force this else ! not a single subsource and random i0_in = 0 j0_in = 0 end if print *,' IN get_params: n_hypocenters, '// : 'hyp_loc_fl, hyp_loc_fw, dl, '// : 'dw, i0_in, j0_in =' print *, n_hypocenters, : hyp_loc_fl, hyp_loc_fw, dl, dw, i0_in, j0_in call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read(nu_ctl,*) i_rise_time call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) tpadl, tpadt, dt ! write(*,'(a,3(1x,f9.3))') ! : 'tpadl, tpadt, dt(sec) :', ! : tpadl, tpadt, dt call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) beta,rho ! write(*,'(a,2f9.3)')"beta,rho :" ! * , beta,rho * gsprd: call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read(nu_ctl, *) r_ref read(nu_ctl, *) nsprd_segs do i = 1, nsprd_segs read(nu_ctl, *) rlow(i), a_s(i), b_s(i), m_s(i) end do * Q: call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read(nu_ctl, *) fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q * Duration: call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) rmin,rd1,rd2,durmin,b1,b2,b3 ! write(*,'(a,7f8.2)')"Dur. info;rmin,rd1,rd2,durmin,b1,b2,b3:", ! * rmin,rd1,rd2,durmin,b1,b2,b3 call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) iwind if (iwind.eq.0) write (*,*)"Windowing info : * tapered boxcar" if (iwind.eq.1) write (*,*)"Windowing info : * Saragoni-Hart" call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read(nu_ctl,*) flocut, nslope call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) damp write(*,'(a, 1x,f8.3)') ! : ' damping: ', ! : damp call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) nfreq,freq1,freq2 ! write(*,'(a,i3,2f8.2)')"No.of freq,lowest and highest freq (Hz)for ! * PSA : ", nfreq,freq1,freq2 nfreq_max = 500 if (iabs(nfreq).gt.nfreq_max) then print *,' ' print *, ' nfreq = ', nfreq, : ' but it cannot exceed ', nfreq_max print *,' STOPPING!!!' stop endif call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read(nu_ctl,*) nfout call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read(nu_ctl,*) (freq_out(i),i=1,nfout) call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) f_stem = ' ' read (nu_ctl,'(a)') f_stem call trim_c(f_stem, nc_f_stem) call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) f_exist = .false. fresp1 = ' ' read (nu_ctl,'(a)') fresp1 call trim_c(fresp1,nc_fresp1) inquire(file=fresp1(1:nc_fresp1), exist=f_exist) if (.not. f_exist) then print *,' file '//fresp1(1:nc_fresp1)// : ' does not exist; QUITTING!!!' stop else ! write(*,'(2a)') ' Input file name for crustal'// ! : ' amplification :'//fresp1(1:nc_fresp1) end if call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) namp1 ! write(*,'(a,i5)') ' No. of lines in the file for crustal'// ! : ' amplification:', namp1 call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) f_exist = .false. fresp2 = ' ' read (nu_ctl,'(a)') fresp2 call trim_c(fresp2,nc_fresp2) inquire(file=fresp2(1:nc_fresp2), exist=f_exist) if (.not. f_exist) then print *,' file '//fresp2(1:nc_fresp2)// : ' does not exist; QUITTING!!!' stop else ! write(*,'(2a)') ' Input file name for'// ! : ' site amplification :'//fresp2(1:nc_fresp2) end if call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) namp2 ! write(*,'(a,i5)') ' No. of lines in the file for site'// ! : ' amplification:', namp2 call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*)iDynamicFlag,pulsingPercent ! write(*,'(a,i8,f8.0)')'Dynamic Flag,pulsing Percent: ', ! : iDynamicFlag,pulsingPercent call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read(nu_ctl,*) iflagscalefactor ! (1=vel^2; 2=acc^2; 3=asymptotic acc^2 (dmb)) * Get types of averages for FAS and PSA call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read(nu_ctl,*) iflagfas_avg call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read(nu_ctl,*) iflagpsa_avg call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) iPapaFlag,PapaGama,PapaNu,PapaT0,PapaPeak write(*,'(a,i2,4f8.3)')"Analytical flag and its Gama,Nu,T0,peak * :", iPapaFlag,PapaGama,PapaNu,PapaT0,PapaPeak call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) iseed, nsims write(*,'(a,1x,i6, 1x, i4)') ! : 'Iseed, No. of trials for PSA calculation: ', ! : iseed, nsims call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*)numberOfSites, isitecoordflag ! write(*,*)"Number of Sites Around the Fault :", ! * numberOfSites print *, ' Type of coordinates (1=lat,long; 2=R,Az; 3=N,E):', : isitecoordflag call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) do i=1,numberOfSites read (nu_ctl,*)siteLocation(i,1),siteLocation(i,2) print *,' SiteLat,SiteLon for Site #', : i,' = ', siteLocation(i,1),siteLocation(i,2) enddo call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) buf_in = ' ' read (nu_ctl,'(a)') buf_in call upstr(buf_in) call trim_c(buf_in, nc_buf_in) if (buf_in(1:1) .eq. 'Y') then move_site = .true. else move_site = .false. end if if (move_site .and. FaultStrike .ne. 0.0) then print *,' ERROR: Cannot request move site if'// : ' the fault strike .ne. 0.0; QUITTING' stop end if isitecoordflag_in = isitecoordflag Do i=1,numberOfSites if (move_site .and. isitecoordflag_in .eq. 1 ) then dum = 0.0 ! eventually might do something with this case else if (move_site .and. isitecoordflag_in .eq. 2 ) then isitecoordflag = 3 r = siteLocation(i,1) az = siteLocation(i,2) xn = r*cos(d2r*az) xe = r*sin(d2r*az) siteLocation(i,1) = faultLength/2.0 + xn siteLocation(i,2) = xe *DEBUG print *, : ' isitecoordflag_in .eq. 2: isite, siteLocation(i,1&2) = ' print *, : i, siteLocation(i,1), siteLocation(i,2) *DEBUG else if (move_site .and. isitecoordflag_in .eq. 3) then if (sitelocation(i,1) .eq. 0.0) then ! move to midpoint siteLocation(i,1) = faultLength/2.0 end if if (siteLocation(i,2) .eq. 0.0) then ! move to end siteLocation(i,1) = faultLength + siteLocation(i,1) end if end if enddo call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read (nu_ctl,*) Islipweight write(*,'(a,i8)')"Slip flag 1=Random, 0=yours * :",Islipweight if (Islipweight.eq.-1) then ! assign unity to all cells do i = 1, nl do j = 1, nw slipWeights(i,j) = 1.0 end do end do else if (Islipweight.eq.0) then call skipcmnt(nu_ctl,cmnts2skip,nc_cmnts2skip) read(nu_ctl,*)((slipWeights(i,j),i=1,nl),j=1,nw) write(*,'(/a)')" *** read slip distribution ***" else if (Islipweight.eq.1) then call createRandomWeight(nl,nw,slipWeights) else print *,' Islipweight = ', Islipweight print *, 'NOT A LEGAL VALUE; STOPPING' stop endif return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function ggnqf(iseed) c generates Gaussian white noise by approx. method of Ross (1985) usum=0. do j=1,12 usum=usum+Ran1(iseed) enddo ggnqf=usum-6. return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine locateHypocentreRandomly(i0,j0,nl,nw) common /seed/ iseed c** draw random hypocentre integer i0,j0,nl,nw i0= nint(Ran1(iseed)*nl) if (i0 .lt. 1) i0=1 if (i0 .gt. nl) i0=nl j0=nint(Ran1(iseed)*nw) if (j0 .lt. 1) j0=1 if (j0 .gt. nw) j0=nw return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine makeItComplex(x,cx,noOfSegPoints) dimension x(*) complex cx(*) do jjj = 1, noOfSegPoints cx(jjj) = cmplx(x(jjj),0.) enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine MavroPapa(totalWave,nptsTotalWave, * iPapaFlag,PapaGama,PapaNu,PapaT0,amag,PapaPeak) c include Mavroeidis and Papageorgiou, 2003 approch common /par/ iFFTsub,dt, npadl, npadt dimension totalWave(*),PapaPulseW(70000),papaW(70000) complex cTotalWave(70000),cPapaPulseW(70000),cPapaW(70000) complex cTemp(70000) integer nptsTotalWave pi=4.0*atan(1.0) peakValue=findPeakValue(nptsTotalWave,totalWave) papaA=PapaPeak papaTp=10**(-2.9+0.5*aMag) papaFp=1.0/papaTp Fp=papaFp samplingRate=1./dt initialnptsTotalWave=nptsTotalWave T0=papaT0 cccccc step1 calculate pusle wave in the time domain PapaPulseW=0.0 do i=1,nptsTotalWave t=(i-1)*dt if(t.ge.(T0-papaGama/2./Fp).and.t.le.(T0+papaGama/2./Fp))then c1=papaA*pi*Fp/papaGama c2=sin(2*pi*Fp*(t-T0)/papaGama)*cos(2*pi*Fp*(t-T0)+papaNu) c3=papaGama*sin(2*pi*Fp*(t-T0)+papaNu) c4=1+cos(2*pi*Fp*(t-T0)/papaGama) PapaPulseW(i)=c1*(c2+c3*c4) else PapaPulseW(i)=0.0 endif enddo cccccc step3-a take totalWave (stochastic) to the frequency domain call padding(TotalWave,nptsTotalWave) call makeItComplex(totalWave,cTotalWave,nptsTotalWave) call fork(nptsTotalWave,cTotalWave,-1.) cccccc step3-b take PapaPulseW to the frequency domain call padding(PapaPulseW,nptsTotalWave) call makeItComplex(PapaPulseW,cPapaPulseW,nptsTotalWave) call fork(nptsTotalWave,cPapaPulseW,-1.) cccccc step 4, keep the phase of stocastic,ccalculate abs(cTotalWave)-abs(cPapaPulseW) cccccc calculate real and imaginary part of the resulus do i=1,nptsTotalWave if(real(cTotalWave(i)).ge.0)then !!very importat when using atan phase=atan(aimag(cTotalWave(i))/real(cTotalWave(i))) else !!very importat when using atan phase=atan(aimag(cTotalWave(i))/real(cTotalWave(i)))+pi endif amp=abs(abs(cTotalWave(i))-abs(cPapaPulseW(i))) x=amp*cos(phase) y=amp*sin(phase) cTemp(i)=cmplx(x,y) enddo cccccc take care of complex conjugates cTemp(nptsTotalWave/2+1)=cmplx(1.,0.)* * cTemp(nptsTotalWave/2+1) do i=1, nptsTotalWave/2+1 cTemp(nptsTotalWave+2-i) = conjg(cTemp(i)) enddo cccccc go back to the time domain call fork(nptsTotalWave,cTemp,+1.) cccccc add them up do i=1,nptsTotalWave totalWave(i)=real(cTemp(i))+PapaPulseW(i) enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function NumberOfPulsingSubs(i0,j0,i,j,nl,nw * ,NoeffectiveSubfaults) c** this function determines how many subfaults are active at this time (i,j). * Dates: 12/01/08 - DMB added integer specification for Rmin, Rmax * and used iabs rather than abs. integer Rmin, Rmax, r,NumberOfPulsingSubs n=0 RMax=max(iabs(i-i0)+1,iabs(j-j0)+1) RMin=Rmax-NoeffectiveSubfaults if(RMin.lt.0)Rmin=0 do ii=1,nl do jj=1,nw r=max(iabs(ii-i0)+1,iabs(jj-j0)+1) if(r.gt.RMin.and.r.le.RMax)n=n+1 enddo enddo NumberOfPulsingSubs=n return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine padding(x,noOfSegPoints) dimension x(70000) do ii=1,100 ncheck=2**(ii) if (ncheck.ge.noOfSegPoints) go to 110 if(ncheck .gt. 70000) then write(*,*)' Burp !! Exceeds array dimensions.' stop endif enddo 110 continue do j = noOfSegPoints+1, ncheck x(j)=0. enddo noOfSegPoints=ncheck return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc * ------------------- BEGIN RAN1 -------------------------- FUNCTION ran1(idum) 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) Copr. 1986-92 Numerical Recipes Software $!6)$-"11j. * ------------------- END RAN1 -------------------------- cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine sample(x,npts,df,nfreq,freq,FA) dimension x(*), freq(500), FA(500) c** takes a large FFT array and samples for nfreq points, c in array FA. Uses avg (sq) FA of 5 around each freq. FA=0.0 do jf = 1, nfreq itarg= ifix(freq(jf)/df)+1 f=df*float(itarg-1) ! dmb replaced itarg with itarg-1 if (itarg .lt. 4) then FA(jf)=-9.999 ! dmb ! FA(jf)=9.999 else c find spectrum around itarg sum = 0.0 ! dmb do ii=itarg-2, itarg+2 sum=0.2*x(ii)*x(ii) + sum ! dmb enddo if(sum.eq.0) then FA(jf)= -9.999 ! dmb else FA(jf) = sqrt(sum) ! dmb end if endif ! dmb ! if(FA(jf).eq.0)FA(jf)=9.999 ! FA(jf) = alog10(sqrt(FA(jf))) ! endif enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine site_coords_in_degrees( : FaultLat, FaultLon, : SiteLat, SiteLon, isitecoordflag, : site_lat_degrees, site_lon_degrees ) pi = 4.0*atan(1.0) dtor = pi/180.0 ! If site coords not in degrees, convert if (isitecoordflag .eq. 1) then ! degrees site_lat_degrees = sitelat site_lon_degrees = sitelon else if (isitecoordflag .eq. 2) then ! R,Az r = sitelat az = sitelon xn = r*cos(dtor*az) xe = r*sin(dtor*az) call km2deg_f( xn, xe, FaultLat, FaultLon, : site_lat_degrees, site_lon_degrees ) else ! assume cartesian coordinates xn = sitelat xe = sitelon call km2deg_f( xn, xe, FaultLat, FaultLon, : site_lat_degrees, site_lon_degrees ) end if return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real function sourceSpectra (f,subfaultDistance, : const,scalingFactor,aMoment,cornerF) c Computes amplitude value of w-squared spectrum at c distance subfaultDistance. c Source terms correspond to subfaultDistance=1 km * Dates: 02/07/09 - Determine gspread and q using gsprd_f and q_f (these * need to be initialized and deallocated in the main program). dimension afreq1(500),afreq2(500),amp1(500),amp2(500) character fresp1*30,fresp2*30 save iflag_source_spectra common /params/ : rho,beta,iKapa,fmax, : fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q, : rmin,rd1,rd2,durmin,b1,b2,b3, : iwind,namp1,namp2,totalMoment, : c4taper, fc4taper, : nsubs, flocut, nslope, pi, twopi common /fnames/ fresp1,fresp2 data iflag_source_spectra/0/ c read site-response functions if (iflag_source_spectra .eq. 0) then if (namp1.ne.0) then nu_site_amp = 2 open (nu_site_amp,form='formatted',file=fresp1) do i=1,namp1 read (nu_site_amp,*) afreq1(i),amp1(i) enddo close (nu_site_amp) endif if (namp2.ne.0) then nu_site_amp = 2 open (nu_site_amp,form='formatted',file=fresp2) do i=1,namp2 read (nu_site_amp,*) afreq2(i),amp2(i) enddo close (nu_site_amp) endif iflag_source_spectra=1 endif w=twopi*f scalingFactor_hf_taper= scalingFactor * : c4taper*(1.0+(f/cornerF)**2.0)/(1.0+(f/fc4taper)**2.0) sa=const*aMoment*scalingFactor_hf_taper*w**2.0* : (1.0/(1.0+(f/cornerF)**2.0)) if (f .eq. 0.0) gamma=0. if (f .ne. 0.0) then Q=q_f(f) gamma=w/(2.*Q*beta) endif anelas=exp(-gamma*subfaultDistance) if (iKapa.eq.0) highc=1./sqrt(1.+(f/fmax)**8.) if (iKapa.eq.1) highc=exp(-pi*fmax*f) ga = gsprd_f(subfaultdistance) am1=1. am2=1. if (namp1.ne.0) am1=amplf(f,afreq1,amp1,namp1) if (namp2.ne.0) am2=amplf(f,afreq2,amp2,namp2) sourceSpectra=sa*ga*highc*anelas*am1*am2 * : buttrlcf( f, flocut, nslope) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine window (i,nstart,nstop,ntaper,wind) c applies cosine tapered window c unit amplitude assumed c Subroutine from Dave Boore. wind=0. if (i.lt.nstart.or.i.gt.nstop) return wind=1. if(i.ge.nstart+ntaper.and.i.le.nstop-ntaper) return pi=3.141593 c Value of wind goes from 0 to 1 from nstart to nstart+ntaper, c then from 1 to 0 from nstop-ntaper to nstop dum1=(nstop+nstart)/2. dum2=(nstop-nstart-ntaper)/2. wind=0.5*(1.-sin(pi*(abs(float(i)-dum1)-dum2)/float(ntaper))) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine wind2 (i,nstart,nstop,ntaper,tw,eps,eta,wind) c applies Saragoni and Hart (1974) window, with parameters c tw (dur), eps (fraction of dur to reach peak), and c eta (fraction of peak ampl. at tw) c The Saragoni and Hart window is applied between c nstart + ntaper, and nstop - ntaper. Outside these c bounds we do a quick cosine taper down to 0 c pi=3.141592654 wind=0. if(i.lt.nstart.or.i.gt.nstop) return wind=1. c Apply Sargoni and Hart window b=-eps*alog(eta)/(1.+eps*(alog(eps)-1.)) c=b/(eps*tw) a=(2.7182818/(eps*tw))**b if (i.ge.(nstart+ntaper).and.i.le.(nstop-ntaper)) then t=(tw/float(nstop-nstart-2*ntaper)) * * (float(i-nstart-ntaper+1)) wind=a*t**b*exp(-c*t) return else c Cos. taper goes from 0 to 1 from nstart to nstart+ntaper, from c 1 to 0 from nstop-ntaper to nstop. if (i.lt.(nstart+ntaper)) then t1=tw/float(nstop-nstart-2*ntaper) wf=a*t1**b*exp(-c*t1) wind=abs(sin((float(i-nstart)/float(ntaper))*pi/2.))*wf else wf=a*tw**b*exp(-c*tw) wind=abs(sin((float(nstop-i)/float(ntaper))*pi/2.))*wf endif return endif end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine writeAcc(ioAcc, totalWave, npts, peakValue, : isim, isite, HypoDistance, faultDistance, fpar) c writes acceleration time history into specified ascii file common /par/ iFFTsub,dt, npadl, npadt character fpar*(*) dimension totalWave(*) integer npts write (ioAcc,"('**********************************************')") write (ioAcc,"('******* Acceleration Time Series ************')") write (ioAcc,"('*** Site #',i4)") isite write (ioAcc,"('*** trial ',i4)") isim write (ioAcc,"('Input Parameters file = ',a)") fPar write (ioAcc,50) npts 50 format (i6,' samples') write (ioAcc,30) dt 30 format ('dt: ',f6.3,' sec') write (ioAcc,'(a)') 'data format: (1x,f8.3,1p,2x,e11.4)' write (ioAcc,10) peakValue 10 format ('maximum absolute acceleration:',f7.2) write(ioAcc,'("fault Dis.(km)=",f8.2)')faultDistance write(ioAcc,'("Hypo Dis.(km)=",f8.2)')HypoDistance write (ioAcc,'(2x,a, 1x,a)') 'time(s)', 'acc(cm/s**2)' do i=1,npts time=(i-1)*dt write (ioAcc,20) time,totalWave(i) enddo 20 format (1x,f8.3,1p,2x,e11.4) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine writeHUS( ioHus, husid, npts, : arias, dur_75_05, : isim,isite,HypoDistance,faultDistance,fpar) * writes husid time history, arias intensity, and 75%-5% duration * into specified ascii file common /par/ iFFTsub,dt, npadl, npadt character fpar*(*) dimension husid(*) integer npts write (iohus,"('**********************************************')") write (iohus,"('******* Husid Time Series ************')") write (iohus,"('*** Site #',i4)") isite write (iohus,"('*** trial ',i4)") isim write (iohus,"('Input Parameters file = ',a)") fPar write (iohus,50) npts 50 format (i6,' samples') write (iohus,30) dt 30 format ('dt: ',f6.3,' sec') write (iohus,'(a)') 'data format: (1x,f8.3,1p,2x,e11.4)' write (iohus,'(a,1x,1pe10.3)') ' Arias intensity = ', arias write (iohus,'(a,1x,1pe10.3)') ' 75%-5% duration = ', dur_75_05 write(iohus,'("fault Dis.(km)=",f8.2)')faultDistance write(iohus,'("Hypo Dis.(km)=",f8.2)')HypoDistance write (iohus,'(5x,a, 7x,a)') 'time', 'husid' do i=1,npts time=(i-1)*dt write (iohus,'(1x,f8.3,1p,1x,e11.4)') time,husid(i) enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine writePar(ioPar, : FaultStrike,FaultDip,h, : FaultLat, FaultLon, : siteLocation,numberOfSites, move_site, : fault_type, : FaultLength, FaultWidth,nl,nw,dl,dw, : specify_length, specify_width, stress_ref, : y, : hyp_loc_fl, hyp_loc_fw, i0_in,j0_in,n_hypocenters, : iseed, nsims, : amag, : stress, : pulsingPercent,iDynamicFlag, : i_rise_time, : iPapaFlag,PapaGama,PapaNu,PapaT0,PapaPeak, : iflagscalefactor, : iflagfas_avg, iflagpsa_avg, : tpadl, tpadt, : r_ref, nsprd_segs, rlow, a_s, b_s, m_s : ) c writes modeling parameters to specified ascii file * Dates: 12/01/08 - DMB added writing of istart, istop, nptsenvelope, maxPointsOfWave character typ*7,fresp1*30,fresp2*30 common /par/ iFFTsub, dt, npadl, npadt logical specify_length, specify_width, move_site character fault_type*(*) dimension siteLocation(300,2) real amag, r_ref, rlow(*), a_s(*), b_s(*), m_s(*) common /params/ : rho,beta,iKapa,fmax, : fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q, : rmin,rd1,rd2,durmin,b1,b2,b3, : iwind,namp1,namp2,totalMoment, : c4taper, fc4taper, : nsubs, flocut, nslope, pi, twopi common /fnames/ fresp1,fresp2 write(ioPar,'(" modeling parameters ")' ) write(ioPar,'("Fault Strike = ",f8.2)')FaultStrike write(ioPar,'("Fault dip = ",f8.2)')FaultDip write(ioPar,'("Fault depth to upper edge = ",f8.2)')h if(.not. specify_length) then write(ioPar,'(a,1x, f6.1)') : ' Fault length from Wells and Coppersmith for fault type '// : fault_type(1:1)//', using a reference stress of ', : stress_ref end if write(ioPar,'("Fault Length = ",f8.2)')FaultLength if(.not. specify_width) then write(ioPar,'(a,1x, f6.1)') : ' Fault width from Wells and Coppersmith for fault type '// : fault_type(1:1)//', using a reference stress of ', : stress_ref end if write(ioPar,'("Fault Width = ",f8.2)')FaultWidth write(ioPar,'("ratio of rupture to s-wave velocity = ",f8.2)') y write(ioPar,'("FaultLat = ",f8.2)')FaultLat write(ioPar,'("FaultLon = ",f8.2)')FaultLon write(ioPar,'("No.of subs along strike = ",i8)')nl write(ioPar,'("No.of subs along dip = ",i8)')nw write(ioPar,'("subfault length = ",f8.2)')dl write(ioPar,'("subfault width = ",f8.2)')dw write(ioPar,'("i_rise_time (1=orig,2=1/f0) = ",i2)') : i_rise_time write(ioPar,'("iseed, nsims = ",1x,i11, 1x,i4)') : iseed, nsims write(ioPar,'("-----------------------------------------------")') write(ioPar,'("input hypocenter at position = ", : 1p,2(1x,e10.3))') : hyp_loc_fl, hyp_loc_fw write(ioPar,'("input hypocenter at subfault = ",2i4)') : i0_in,j0_in write(ioPar,'("n_hypocenters = ",i4)') n_hypocenters write(ioPar,'("Mag. = ",f8.2)')amag write(ioPar,'("-----------------------------------------------")') write(ioPar,'("FFTpoints = ",i8)')iFFTsub write(ioPar,'("dt (sec) = ",f8.2)')dt write(ioPar,'("beta (km/s) = ",f8.2)')beta write(ioPar,'("density (rho), gr/cm3 = ",f8.2)')rho write(ioPar,'("pulsing Percentage = ",f8.2)')pulsingPercent write(ioPar,'(a,i1)') 'iflagscalefactor (1=vel^2; 2=acc^2; '// : '3=asymptotic acc^2 (dmb)) = ', : iflagscalefactor write(ioPar,'("flocut, nslope = ",1x, f6.3, 1x, i2)') : flocut, nslope write(ioPar,'("iflagfas_avg = ",i1)')iflagfas_avg write(ioPar,'("iflagpsa_avg = ",i1)')iflagpsa_avg write(ioPar,'("stress parameter (bars) = ",f8.2)')stress typ='fmax' if (iKapa.eq.1)then write(ioPar,'("kappa = ",f8.2)')fmax endif if (iKapa.eq.0)then write(ioPar,'("fmax = ",f8.2)')fmax endif write(ioPar,'("-----------------------------------------------")') write(ioPar,'(" Corner Frequency ")' ) if (iDynamicFlag.eq.1)then write(ioPar,'("Dynamic Corner Frequency Flag is ON = ",i4)') * iDynamicFlag else write(ioPar,'("Dynamic Corner Frequency Flag is OFF = ",i4)') * iDynamicFlag write(ioPar,'("Corner Frequency is static")') endif write(ioPar,'("-----------------------------------------------")') write(ioPar,'(a)') : ' fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q = ' write(ioPar,*) fr1, qr1, s1, ft1, ft2, fr2, qr2, s2, c_q write(ioPar,'("-----------------------------------------------")') write(ioPar,'("rMin in duration model = ",f8.2)')rmin write(ioPar,'("1st hinge in duration model = ",f8.2)')rd1 write(ioPar,'("2nd hinge in duration model = ",f8.2)')rd2 write(ioPar,'("1st slope in duration model = ",f8.2)')b1 write(ioPar,'("2nd slope in duration model = ",f8.2)')b2 write(ioPar,'("3rd slope in duration model = ",f8.2)')b3 write(ioPar,'("Min of duration = ",f8.2)')durmin write(ioPar,'("-----------------------------------------------")') write(ioPar,'(a)') : ' gspread: i, nsprd_segs, r_ref, rlow, a_s, b_s, m_s' do i = 1, nsprd_segs write(ioPar,*) : i, nsprd_segs, r_ref, rlow(i), a_s(i), b_s(i), m_s(i) end do write(ioPar,'("-----------------------------------------------")') if (iwind.eq.0) write (ioPar,120) if (iwind.eq.1) write (ioPar,130) 120 format ('window applied = tapered boxcar') 130 format ('window applied = Saragoni-Hart') write(ioPar,'("-----------------------------------------------")') if (namp1.eq.0.and.namp2.eq.0) * write (ioPar,'(a)') ' no local amplification' if (namp1.ne.0) write (ioPar,200) fresp1 if (namp2.ne.0) write (ioPar,200) fresp2 200 format ('amplification as in file ',a30) write(ioPar,'("-----------------------------------------------")') write(ioPar,'("-----------------------------------------------")') write(ioPar,'(" Analytical Pulse parameters ")' ) write(ioPar,'("-----------------------------------------------")') if (iPapaFlag.eq.1)then write(ioPar,'("Analytical Flag is ON = ",i4)')iPapaFlag write(ioPar,'("Analytical Gama = ",f8.3)')PapaGama write(ioPar,'("Analytical Nu = ",f8.3)')PapaNu write(ioPar,'("Analytical T0 = ",f8.3)')PapaT0 write(ioPar,'("Analytical Peak = ",f8.3)')PapaPeak else write(ioPar,'("Analytical Flag is OFF = ",i4)')iPapaFlag endif write(ioPar,*) write(ioPar,'(a,3(1x,f9.3))') : 'tpadl, tpadt, dt(sec) :', : tpadl, tpadt, dt do i = 1, numberOfSites write(ioPar,'(a, 1x,i3, 1x,a, 1x,f8.2, 1x, f8.2)') : 'For site', i, 'siteLocation coordinates 1&2 = ', : siteLocation(i,1), siteLocation(i,2) end do if (move_site) then write(ioPar,'(a)') : 'Site may have been moved to midpoint or end of '// : 'surface projection of upper edge of fault' end if return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine writePSA_FA(ioPSA, freq,averagePGA,averagePGV, : averagePSA,nfreq,nsims,damp,amag, : siteLat,siteLon,depth,r,azimuth,faultDistance,averageFA,isite, : HypoDistance,fPar, nl, nw) ! : HypoDistance,fPar) c writes simulated average PSA spectrum at nfreq frequencies c into specified ascii file. First column is frequency in Hz, c second column is PSA value in cm/s**2 * Dates: 12/31/08 - Write sqrt(nl*nw)*psa, sqrt(nl*nw)*fas * 01/06/09 - Write pgv, pga as first two values * 01/13/09 - Write sqrt(nl*nw)*psa, sqrt(nl*nw)*fas and add SD, PSV. * - Reverse order of FAS and PSA dimension freq(*),averagePSA(*),averageFA(*) character fname*20, fPar*(*) real omega, omega_sq pi=4.0*atan(1.0) twopi = 2.0 * pi sqrt_n = sqrt(float(nl*nw)) write(ioPsa,"('**********************************************')") write(ioPsa,"('******* Site #',i4)") isite write(ioPsa,40) nsims write(ioPsa,50) int(damp) write(ioPsa,20) nfreq amag=anint(amag*10)/10. write(ioPsa,'("Mag. = ",f8.2)')amag write(ioPsa,'("siteLat = ",f8.2)')siteLat write(ioPsa,'("siteLon = ",f8.2)')siteLon write(ioPsa,'("depth = ",f8.2)')depth write(ioPsa,'("R = ",f8.2)')R write(ioPsa,'("azimuth = ",f8.2)')azimuth write(ioPsa,'("fault Dis.(km)= ",f8.2)')faultDistance write(ioPsa,'("Hypo Dis.(km)= ",f8.2)')HypoDistance write(ioPsa,"('Input Parameters file = ',a)") fPar write (ioPsa,'(1x,a)') 'data format: '// : '(1x,f10.3, 1p, 1x,e11.4, 5x,e11.4,'// : ' 2x,e11.4, 1x,e11.4)' write(ioPSA, '(3x,a, 5x,a, : 1x,a, 1x,a, : 1x,a, 3x,a)') : 'freq(Hz)', 'per(s)', : 'FAS(cm/sec)', 'avgPSA(cm/s**2)', : 'avgPSV(cm/s)', 'avgSD(cm)' freq_dum = -1.0 write (ioPsa,'(1x,f10.3, 8x,a3, : 9x,a3, 1p, 13x,a3, : 2x,e11.4, 9x,a3)') : freq_dum, 'NaN', : 'NaN', 'Nan', : averagePGV, 'Nan' freq_dum = 999.9 per_dum = 0.0 write (ioPsa,'(1x,f10.3, 1x,f10.3, : 9x,a3, 1p, 5x,e11.4, : 10x,a3, 9x,a3)') : freq_dum, per_dum, : 'Nan', averagePGA, : 'NaN','NaN' do i=1, nfreq omega = twopi*freq(i) omega_sq = omega**2.0 per_dum=1.0/freq(i) if(averageFA(i).gt.0) then ! dmb write (ioPsa,'(1x,f10.3, 1x,f10.3, : 1p, 1x,e11.4, 5x,e11.4, : 2x,e11.4, 1x,e11.4)') : freq(i), per_dum, : averageFA(i), averagePSA(i), : averagePSA(i)/omega, averagePSA(i)/omega_sq else per_dum=1.0/freq(i) write (ioPsa,'(1x,f10.3, 1x,f10.3, : 9x,a3, 1p, 5x,e11.4, : 1x,e11.4, 1x,e11.4)') : freq(i), per_dum, : 'NaN', averagePSA(i), : averagePSA(i)/omega, averagePSA(i)/omega_sq endif enddo write (ioPsa,*) write (ioPsa,*) 20 format (i4,' samples') 40 format('Average PSA and Fourier spectrum over ',i5,' trials') 50 format ('damping: ',i3,'%') return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc include 'exsim_dmb_included_subprograms.for' ! include '\forprogs\accsqint.for' ! include '\forprogs\dcdt.for' ! include '\forprogs\locate.for' ! include '\forprogs\locate_d.for' ! include '\forprogs\zbrent.for' ! include '\forprogs\dist_3df.for' ! include '\forprogs\km2deg_f.for' ! include '\forprogs\yintrf.for' ! include '\smsim\gsprd_f.for' ! include '\smsim\gsprd_q_f.for' ! include '\smsim\gsprd_q_avg_f.for' ! include '\smsim\q_f.for'