Corrections and modifications to SIMQKE1 (VAX version)
This note describes some corrections and modifications made to SIMQKE1 (as supplied by NISEE at UC Berkeley) by Edmund Booth (edmund@booth-
The note first lists the changes to the Fortran code made by Edmund Booth to correct apparent errors; for convenience, the changed parts of the code are listed. The note also describes some changes made to improve the usability of the program; the associated changes to the Fortran code are not however listed.
A: CORRECTIONS
A1: 'Baseline correction'
A 'baseline correction' is applied to make the velocity at the end of the record equal zero. The baseline correction in the original program appears to be based on a modified version of
Boyce, W.H. Treatment of earthquake accelerograms. BNZSEE 3(4) Page 163-172
There are two major problems with the adopted procedure:
|
·
|
The initial acceleration as generated does not equal zero.
|
|
·
|
The generated final velocity is also non-
|
In fact, the baseline correction method adopted originally appears overcomplex and inappropriate.
The procedure has been modified so that an acceleration pulse of the same shape as the chosen intensity envelope (ie figure 2 in the SIMQKE manual) is applied to bring the final velocity to zero. The peak value of this acceleration pulse, expressed as a ratio of the peak ground acceleration, is written to an output file. The correction has been changed to take place after the accelerations have been adjusted so that there is only one peak equal to the specified peak ground acceleration. This change means that the pga of the record may differ very slightly from the specified pga value. However, the error will not be greater than the normalised value of the peak acceleration pulse referred to above; values not exceeding 5% have been found so far.
The procedure has been modified so that an acceleration pulse of the same shape as the chosen intensity envelope (ie figure 2 in the SIMQKE manual) is applied to bring the final velocity to zero. The peak value of this acceleration pulse, expressed as a ratio of the peak ground acceleration, is written to an output file. The correction has been changed to take place after the accelerations have been adjusted so that there is only one peak equal to the specified peak ground acceleration. This change means that the pga of the record may differ very slightly from the specified pga value. However, the error will not be greater than the normalised value of the peak acceleration pulse referred to above; values not exceeding 5% have been found so far.
A2: Rise time for ICASE=1
The manual does not state what envelope is used if ICASE =1. It might be assumed that ICASE =1 corresponds to no envelope. Actually, this is not the case and an envelope is applied; between 0 and 2 seconds, a linear ramp up is applied, then after two seconds, the motions are not enveloped.
A3: A0 value for 'exponential' envelope (ICASE=3)
The variable A0, which sets the shape of the 'exponential' envelope option (ICASE=3) has been calculated internally to normalise the envelope to a maximum value of 1. The associated code is as follows.
C EXPONENTIAL INTENSITY ENVELOPE
3004 DO 3006 KK=2,NACCG
C!!!
C The next lines set AO for ICASE=3 such that the max intensity is 1
C!!!!
T00=(Log(BETAO) -
GO TO 3011
A4: Plateau time of 'compound' envelope(ICASE=4)
For ICASE=4, a 'compound' envelope is taken as bounding the generated acceleration time history. However, figure 2 of the user manual is inconsistent; the variable TLVL is defined in the figure as the elapsed time between the start of the motion and the end of the envelope plateau , whereas the formula in the figure takes TLVL as equal to T_LEVEL, the duration of the envelope. The code is inconsistent in that it assumes that TLVL equals both the elapsed time at the end of the plateau and the plateau duration, which is wrong unless the initial rise time equals zero. For the 'trapezoidal' envelope, TLVL is taken as the duration of the plateau, which is inconsistent with Figure 2.
The fix adopted has been to make TLVL consistently equal to the elapsed time (ie as defined in Figure 2). This necessitates the some changes to the Fortran code, as follows; the changes are identified by comment lines starting C!!!.
C DEFINE SLOPES OF ENVELOPE
C
IF (ICASE.GT.2) GO TO 6
IF(TRISE.GT..0) GO TO 33
TRISE=0.25*DUR
TLVL=0.
33 IF(ICASE.LE.1) GO TO 7
8 FTC1=1./TRISE
C
C!!!
C The original code seems to confuse the time length of the pateau and
C the time at the end of the plateau.
C TLVL has been taken as the elapsed time at the end of the plateau
C!!!!
C
FTC2=-
7 FTC1=0.5
C!!! FTC1=0.5 corresponds to a ramp up time of 2 secs for ICASE=1
C
FTC2=0.
6 WRITE(6,114) WA,TA,NFQ,WBEGIN,WEND
C TRAPEZOIDAL INTENSITY ENVELOPE
3003 IF(ICASE.LE.1) GO TO 18
TX=TRISE
GO TO 19
18 TX=2.
C!!! TX=2. corresponds to a ramp up time of 2 secs for ICASE=1
C
C DEFINE MAXIMUM HEIGHTS IN TERMS OF SLOPES
C
19 DO 16 KK=2,NACCG
TI=(KK-
FT=FTC1*TI
GO TO 16
C
C!!!
C The original code seems to confuse the time length of the pateau and
C the time at the end of the plateau.
C TLVL has been taken as the time at the end of the plateau
C!!!
C
15 IF(ICASE.LE.1) GO TO 28
IF((TI-
GO TO 16
29 FT=1.+(TI-
C COMPUTE ACCELERATION
C
16 ACCG(KK)=ACCG(KK)*FT
GO TO 3011
C
C EXPONENTIAL INTENSITY ENVELOPE
3004 DO 3006 KK=2,NACCG
TI=(KK-
GO TO 3011
C
C!!!
C The original code seems to confuse the time length of the pateau and
C the time at the end of the plateau.
C TLVL has been taken as the time at the end of the plateau
C!!!!
C
C COMPOUND INTENSITY ENVELOPE
3007 DO 3010 KK= 2,NACCG
TI=(KK-
FT=(TI/TRISE)**IPOW
GO TO 3010
3008 IF ((TI-
3009 FT=1.0
3010 ACCG(KK)=ACCG(KK)*FT
3011 CONTINUE
C
A5: Generation of multiple time histories
Currently, the VAX version of the code merely regenerates identical time histories where multiple histories are requested (ie where NPA>1 is chosen). The following changes have been made to the Fortran code to correct this.
IX=(IIX/2)*2+1
C!!! The next line has been added to get statistically independent records to be
C generated from the same target spectrum
PA(1)=random(ix)
C!!!
C
C LOOP OVER NPA, NUMBER OF ARTIFICIAL EARTHQUAKES DESIRED
DO 585 NTOTAL=1,NPA
C!!! The following line has been modified so that a new random number is used
C for successive phase angles for successive time histories
31 pa(i)=random(-1)
C!!!
C ACCELERATION COMPUTATIONS
C
8603 NACCG=DUR/DELT+1.000001
C for successive phase angles for successive time histories
31 pa(i)=random(-
C ACCELERATION COMPUTATIONS
C
8603 NACCG=DUR/DELT+1.000001
B: CHANGES TO INPUT FILES
The following changes have been made to the way the input file is read, in order to improve usability. The associated changes to the Fortran code are not listed here.
|
·
|
Lines 1,3,5,7,9,11 and 13 are treated as comment lines, and are ignored. This makes it much easier to understand (and edit) the input file.
|
|
·
|
Line 2 is an extra data line, consisting of up to 8 characters (1A8) describing the IDENT of the run (eg MYFILE). IDENT is also prompted for from the screen at the time of running, and must be typed in by the user.
|
|
·
|
The data file name must have the extension .DAT and must correspond to IDENT, eg MYFILE.DAT
|
|
·
|
The variable TLVL describing the plateau of the generated time history envelope has been revised to become the elapsed time at the end of the plateau rather than the duration of the plateau. (The original version is generally ambiguous in this respect and in error for the 'compound' envelope type).
|
|
·
|
Consistent units (eg m/s and m/s2) are assumed, rather than in/sec and g.
|
|
·
|
The variable A0, which sets the shape of the 'exponential' envelope option 3 has been calculated internally to normalise the envelope to a maximum value of 1. AO is therefore not required as input. This means that the format for line 6 becomes (I5,5F10.4,I5).
|
C: OTHER MODIFICATIONS
The program has also been modified in other ways as follows.
