-
Notifications
You must be signed in to change notification settings - Fork 11
/
hazgridXnga13l.f
12216 lines (11838 loc) · 551 KB
/
hazgridXnga13l.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
c--- hazgridXnga13l.f for USGS PSHA runs, Last changed 09/21/ 2016. Long header version.
c insure that Campbell CEUS GMM (index 10) is called with negative ipia when vs30 > 1200 m/s. Mod sept 21 2016 SH
c
c nov 2,2015: repair a little bug with "oktogo" variable so that 2-s SA can be
c computed with Frankel and others
c nov 5, 2015: repair initializations for "getToro" so that 2-s SA works with this GMM.
c In Toro, the period index for 2-s SA is 7.
c Oct 20, include the Atkinson_2015 Induced seismicity GMM "a15ind" with index 51.
c Nov 18, 2015. Use dtor(1) to model source depth in A15ind. Previously, set at 5 km.
c Aug 31 2015: include the Hawaii terms to BA08 GMM ( use index 50 for Hawaii GMM.)
c change the unit=15 to unit=95 because of a conflict with station file open.
c
c Crustal-source Zhao now allowed June 24 2015. Use index 49 for Zhao GMM (29 is BSSA)
c June 2015: include 0.4s for ABslab (getABsub) calcs.
c May 27 2015: changed a few coeffs for soil response in ASK2013 to agree with Spectra and Peer reports
c May 15 2015: add 2.5s coeffs to ASK2013, CB13, CY2013, Idriss13, and AB03(slab)
c April 7, 2015: Correct c(0.4s) in Zhao et al. to 0.01 as in Table 4. Previous
c value was from an email from Zhao. His email set c(0.4s) to 0.006.
c
c Feb 13, 2015: Define a clamp for the 14 periods in the three newer Atkinson&Pezeshk
c CEUS GMMs. Remove the "iq" index from these subroutines.
c Note: Units of PGV are cm/s in the Atkinson CEUS GMMs. Provisional
c CEUS PGV clamp is set at 400 cm/s. This value should be discussed
c Feb 10, 2015: Use a coarser test to accept 3.33 hz as equiv. to 0.3-s in newer
c CEUS Atkinson GMMs.
c Dec 9 2014: add 0.15 and 0.25 s coeffs. to getABsub
c Oct 30 2014: icampceus array repair for pga.
c Feb 27: add coeffs for T=4 and 5 s to the ABsub routine for in-slab events.
c Feb 20: add hard-rock site term for Zhao et al relation. This change does not
c affect the BC rock calcs, but will affect B- to B+ rock for intraplate source.
c Also extended number of periods available to 22 (from 11 in previous vers) Zhao subr.
c
c ABsub: added 1.5s coeffs Jan 24 2014 SH. see AB03slab.1p5s.f in my Srcf dir.
c ABsub: add 4 and 5 s coeffs (extrapolations) Feb 27
c jan31 2014: include 1.5s for the three "GailTable" tables
c Jan 24, 2014: getAB06: Add coeffs for 1.5s jan 24 2014 SH. See AB06.1p5s.f for the details.
c Jan 24, 2014: getCampCEUS: include interpo. coeffs for 1.5s spectral period see "campCEUS.1p5s.f" for details
c jan 31, 2014: getSomer: add 1.5s coeffs. (CEUS relation used for finite src)
c Jan 24, 2014: getSilva: include interpo. coeffs for 1.5s spectral period see "silva.1p5s.f" for details
c Jan 23, 2014: getToro: include interpo. coeffs for 1.5s spectral periodsee "toro.1p5s.f" for details.
c Jan 22 2014: Correct Idriss sigma_aleatory to conform to Eq Spectra article
c Jan 21 2014: Correct NAAsub (BCHydro) sigma to 0.74.
c Jan 21, 2014: NAAsub: Central DelC has been corrected to -0.3 for all periods
c intraplate source. More work will be needed when we get info from K Addo.
c Jan 7, 2014: slight modification of getNAAsub for PGArock. from Peter Powers email.
c
c Jan 3 2014: Apply same high-limit (40 hz) median clamp in AB06',A08',Pez11.
c 12/16/2013: Update ASK2013 a1 and vlin for some short periods. from Sanaz R email
c OCT 1, 2013: use Mcap of 7.5 on s.d. computation in Idriss2013.
c Oct 18 2013: standardize Rrup for Idriss too. Lowers dip-slip hazard on hanging wall
c compared to previous. Footwall hazard same. Strike-slip same.
c Aug 29 2013: standardize Rrup and Rx to OpenSHA from P.Powers notes. This mod
c affects ASK13, CB13 and CY13. It does not affect other GMMs.
c 9/10/2013: In CB13, zbot is initialized for the first time. Several unused subroutines removed.
c 8/28/2013: CB13: always calculate phi_lnY(23) (a PGA quant) early in subroutine. It is needed
c for all spectral periods sigma.
c 8/27/2013: Include an mmin matrix option. Previously code just had an mmax distribution. This option
c is invoked if maxmat=-2. not finished
c 8/27/2013: update c0 vector in CB13 model.
c 8/21/2013: Inline coeffs in CB13. Do not need coeff. file any longer. Standardize
c Zhyp in CB13.
c 8/13/2013: z1_rock = -1 ASK13.
c 8/07/2013: common/epistemic/ now same size everywhere.
c 8/02/2013: The correlation coeff vector rho was updated in CB13. It is now
c stored as a data vector rather than being read in.
c 8/01: New CY2013 and ASK2013 models. Both use inferred Vs30 as the default condition
c BSSA NGAW(2): clin and Vclin updated July 2013.
c 7/09/2013: Include tapered GR distribution option. Invoked if magref >=6.5.
c Magref was previously unused. Now, it is m_c of the tapered GR (Pareto distribution)
c QA 7/10/2013: put mZ_tor inside M loop in CY2013.
c CB13_SPEC Modified May 22 2013 to make depth to 2.5 km/s rock 0.398 km for hardrock model.
c 5/21/2013: correct c11 definition in gk13v2.
c 5/20/2013: Add May 18 BSSA. Correct definition of z1km used in GK13 and AS08.
c 5/17/2013: correct some e5 coeffs in bssa2013drv. Make cDPP = 0 in CY2013. (no directivity)
c 5/16/2013: Update AS2013 NGA-W. Use the W&C M(A) to infer a fault width, as in CB13.
c 5/15/2013: update CY NGA-W to May version. 5/16: correct small typo from Chiou email
c 5/13/2013: Add anelastic attenuation in CB13, for R>80 km. Bozorgnia email May 2013.
c 5/10/2013: Further update hypocenter in CB13. USE W&C M(A) to estimate flt width
c
c 5/09/2013: Increase the R index to 310 to allow di=1 and Rmax of 300 km. Steering committee
c recommendation is to step up Rmax to 300 km to capture afault influence in Arizona.
c Put the hypocenter 1/2 way downdip in CB13 virtual faults. Was 3/4 down previously.
c
c 5/08/2013: Deep seismicity can be treated with a staircase distribution of depths
c These depths and the longitudes where they change are input in the dtor line.
c There can be up to three deep seismicity levels.
c 4/16/2013: BSSA13 running for all 107 periods incl PGV (index 1).
c 4/12/2013: CY2013 updated. Corrected a line about HW_taper3 in AS-2013 routine (their only update).
c 4/11/2013: Update Idriss NGA-W GMPE coeffs to latest available (emailed 4/10)
c 4/10/2013: correct logic for fixed-strike source and agrid with header (iflt=20)
c Increase number of maximum-magnitude zones to 15 (controlled by nzonex parameter)
c 4/09/2013: modify CB13 and GK13. Modify Basin and Range Q for GK13 to 205
c 4/04/2013: asum first dim up dimension to 40 to accomodate big range of magnitudes.
c 4/02/2013: added a definition for di2 = di/2 that was needed (PC Liu Discovery). also small
c repair of a line controlling logic for variable Vs30 input model.
c 4/01/2013: Add GK13. Calif Q is 156.6 here (new). Same index as GK12.
c 3/20/2013: clean up declarations. dont use dimension without specifying type. This
c step is a disambiguation where subroutines specify "implicit none". Some array items
c were replaced with scalar items in calls to some subroutines. Work with Morgan 3/20
c 3/19/2013: Revise BSSA to the March version. Has large GMPE Std. dev. Add sdi option to BSSA13.
c 3/152013: Include AS08 index 24.
c 3/06/2013: incl Mar 2013 version of CY NGA(2) GMPE- this one removed and the april one replaced.
c 3/05/2013: include Feb 2013 version of AS NGA(2) GMPE.
c 2/25/2013: include Feb 2013 version of CB NGA(2) GMPE. work in progress.
c 2/22/2013: add option to compute inelastic spectral displacement (Units: cm). Fcn of dy.
c To invoke option, specify dy as arg2 and the value of dy as arg3 (for example 10.0 cm).
c 1/18/2013: corrections to sigma in Idriss GMPE. Loop corr in BSSA preparation.
c
c 1/09/2013: allow user to specify aleatory sigma with arg2=fixsigma and arg3= value
c Example using fixed sigma option: hazgridXnga13l CAmapC_21.in fixsigma 0.5
c 1/09/2013: improve the asum() accumulation based on P. Powers observation
c about "banding" of event-rates at certain latitudes. New vers. has no banding.
c The effect of this change is visible when comparing with OpenSHA but otherwise too subtle.
c 12/28/2012: add BSSA12 index 29
c 12/27/2012: add Idriss2013 index 38. This one has linear siteamp function.
c 12/26/2012: add CY2012 index 33 (removed Mar 2013) Use CY2013.
c 12/24/2012: add AS2012 model index 34. As in CB12, you can vary dip of virtual flts
c by using icode values in the range 0 to 9, or 0 to -9 to put all sites on footwall side
c DEC 20 2012: iconv becomes icode to add options associated with WUS gmpes as well as CEUS gmpes
c add Campbell-Bozorgnia 2012 GMPE for wus. Icode for CB12 set up to vary dip of random sources:
c icode=0 90 =>d; 1=>80 d, 2=>70 d, and so on. Index 32.
c Add GK12 model with basin effect. index 39. Q_s is 435 everywhere in the initial model setup.
c this version has the long header records (896 instead of 308 byte)
c 11/16/2012: add NAAsub corresponding to the BCHYDRO GMPE of 2010. Index=31. For intraplate sources.
c we dont intend to use NAAsub for subduction sources in this code (see hazSUBXnga.test)
c GetGeom : use BA nonlinear siteamp, from AF hazgridXGT.f.
c GetABsub: do not modify. Use the original 2003 formulation of siteamp. Some corrections to getABsub
c were discovered by Pengsheng and fixed in this code Jan 9, 2013. SHarmsen.
c
c more comments in previous versions.
c Some history of code mods:
c 7/22/2010: hazgridXnga13l with double precision for the summation step prob(,,,)=prob()+blah blah
c Idea is that the small distant source rate may be omitted during summation due to word-size.
c Testing shows d.p. accumulation is a very minor improvement. But does not impact run time so we
c implement it here.
c 10/31/2012: Add AB06', A08', Pez11 CEUS GMPEs initially set up for A-rock or BC-rock only. The
c A->BC conversion uses Gail Atkinson's factors very different from Frankel's factors.
c
c 7/19/2010: Add Geomatrix subduction (previously just inslab). iatten=14 for subd.
c 7/13/2010: Add Zhao interface and intraplate GMPEs. These have slightly different
c sigma_t and considerably different medians.
c 7/14/2010: add ABsub subduction (previously only inslab). gmpe 16 for Cascadia,
c 17 for global
c 2/03/2010: Add option for fixed strike with strike direction a fcn of site location.
c To invoke this new option, make iflt -2
c The field of strikes is in a binary file that replaces the strike direction & is written
c on that line of the input file. Use iosubs to read and write this file.
c The trial WUS strike directions were written by azimuth_to_pole.f. An Euler pole
c was put at 45N, 116W for trial runs. Some places (PACNW, So CAL) were modified
c 10/13/2009: correct TP05 c1(0.3 s) for BC rock.
c Oct 9 2009: put M(m) into an array xmw() where M is the moment mag for index m.
c Conversion should be done just once, and then subsequently should access xmw(m) rather
c than recomputing the CEUS mb->M .
c Jan 7 2009: add PGV coeffs to CY-NGA relation. PGV available in Mar 2008 Eq spectra.
c Jan 9 2009: zero out the avghwcy array when entering routine. necessary on PC runs
c May 8 2009: correct some D-soil coeffs, c6, in getABsub (10hz, 3.33, 2 hz)
c May 13 2009: add E-soil coeffs, c7, to getABsub. Include B rock, prev. BC only
c June 17 2009: correct tau_nl in getCY2007H
c Previous version had tau a decreasing fcn of distance. This version is more accurate.
c June 18 2009: use current xmag to define the magnitude index, m_ind, in rjbmean array
c This mod assumes that xmag, or M, at the time it is looked at is Moment Magnitude, and
c that the calculation of rjbmean was for a M(SRL) or M(A) relation.
c
c Dec 18 2008: vulnerability in output file names for the .m and .p files is finally repaired.
c Nov 20 2008: add Zhao et al. atten. for inslab and interface source. Need variety of models at LP Sa, but
c AB03 is only good to 3 s. Zhao's inslab goes to 5 s Sa.
c Nov 19 2008: increase set of periods available with Geomatrix inslab attn.
c mod Oct 22 2008. M (saturation) limit at 8.0 (AB03, BSSA v93 #4, p 1709)
c August 2008: Mmax may be treated as a distribution for the cases iflt=3 and iflt=4.
c these are CEUS or CENA flags indicating mblg is the native magnitude and
c if iflt=3 mblg will be convertd to moment mag. using an Arch Johnston conversion;
c if iflt=4 mblg will be convertd to moment mag. using Atknson-Boore 95 conversion formula.
c wtmj_cra, wtmj_ext
c wtmab_cra, wtmab_ext complementary distribution vectors.
c Vector sizes increased to allow up to mb 8.05 Feb 1 2012. Largest mb corresponds
c to a moment mag approaching 9.
c These vector values are defined in data statements below (or f95 equiv. of data statements).
c These vectors can be modified to change the epistemic range of Mmax in the CEUS (CENA).
c The values are relative to the local magnitude (mb) and its conversion to moment
c magnitude. NSHMP and this code use two conversions, Johnston and AB95.
c These are designated "j" and "ab" in above vectors. The first distribution value corresponds
c to mb 5.05. The second corresponds to mb 5.15, and nth corresponds to 5.05+ (n-1)* 0.1. 0.1 is the
c expected dm in the below code. If mmin or dm changes, you would have to change the distributions above.
c The upper tail endpoint that has been tested corresponds to Mmax 7.7 (moment magnitude).
c "cra" defines mmax distribution in USGS Craton. "ext" defines mmax distribution in USGS Extended Margin.
c logical arrays craton and margin are read in. The arrays are defined over the CEUS standard map grid.
c If i is the geographic location index, craton(i)=.true. if i is in the craton. margin(i)=.true. if i
c is in the extended margin. Both are false if i is in an undefined location (offshore).
c Undefined results
c in a default offshore Mmax which is the magmax read in in the input file. This is a scalar Mmax (no distribution).
c This option also requires that you read in an mmax array which defines the absolute upper limit of Mmax at
c each grid point (maxmat = 1 must be selected). If you set maxmat to 0, the input magmax will be applied at
c all geographic locations (no craton / margin distinction).
c maxmat has additional option in 2013: maxmat>1 means read in a model of Mmax regions. In each Mmax
c region, a distribution of Mmax values is read in. The hazard is calculated based on the CCD of Mmax
c in each of these regions (up to 10 such zones are allowed)
c
c sept 30 2008: define dy inside getBooreNGA (prior omitted)
c October 26 2007: add Oct version of CY NGA relation. Subroutine CY2007H.
c modified sept 12, 2007: does not call erf() within intraplate subroutines.
c This seems to be required for Windows PC gfortran. Other subroutines will
c also have to be modified similarly (such as CEUS models).
c Nov 26, 2007: fix a bug in getAB06 for the non-default stress factor case such as 200 bar
c June 4 2008: deaggregate at a site
c May 21 2008: further consolidate CEUS background runs, from 16 to 2. The
c remaining 2 are for Johnston mb to M and AB mb to M, respectively.
c The mmax branches are included as weights. Two binary logical files craton and margin are
c read in iff abs(iflt)= 3 or 4.
c --- iflt <= 0. Uses only point src.
c --- iflt=-3 uses point src but also uses Johnston mb to Mw New sept 2008
c --- iflt=-4 uses point src but also uses AB mb to Mw
c April 30 2008: Read in precomputed array of mean distances to randomly oriented
c faults. This array is called rjbmean and is indexed by mag and distance. Reading
c in this array from file 'meanrjb.bin' can save a few seconds of CPU time
c per CEUS background-hazard run (16 of them). Indexes ilat and ilon no longer needed.
c Feb 1 2013: the mean distance from src S needed to be increased to allow
c Mmax of 8.1 for the CEUS SSC logic tree branches. The new mean distance
c array is called rjbmean.bin.Aeast and has 26 mag bins instead of 16.
c
c aug 19, 2008: clamp is consistent for short period defined from 0.5 > T > 0.02 s.
c Apr 10 2008: clamp in getFEA is now used. Was commented out prior to today. (T. Cao noticed this)
c feb 5 2008: add 3s coeffs. to absub and Geomatrix atten models.
c March 15 2008: Hardrock for AB06 always called if Vs30>=2000 m/s and can be called if Vs30>1500 m/s.
c March 15 2008: Hardrock coeffs used for Frankel, Toro, TP05, Silva if Vs30>=1500 regardless of whether negative
c iatten() index is used. This is meant to assist users. However, there is a debate about what constitutes
c A-class rock in the CEUS versus A in the WUS. May need an AB class in nhbd of 1500 m/s.
c March 20 2008: BooreNGA308 subroutine replaces BooreNGA407: Only effect is to
c nonlinear calculations and pga_nl. Does not affect 760 m/s or BC rock calcs.
c This update should be used for the SoilD and SoilC PSHA calcs.
c April 8 2008: Toro 2-hz and 3.33-hz c1 coeffs checked and changed.
c july 6 2008: add 2.5 and 25 hz to getSilva and getTP05. For the latter,
c use cubic spline interp from Numerical Recipes. This seems to do a good job
c putting the median between the nearest available period medians over Rcd range<1000 km
c
c Dec 3, 2007: Toro CEUS finite fault corr. uses Mw to make the correction.
c Sept 20, 2007: new mmax rule: if the mmax matrix indicator is -1, use the
c minimum of the scalar Mmax and the matrix Mmax for that source. The purpose
c of this mod is to perform a rate calculation for M<6.5 earthquakes which could
c have M<6.0 in a few places, such as the Creeping Section of the SAF. SH 9/20/2007
c
c Greater Ztor allowed aug 15 2007 to accomodate very deep Benioff zones
c revised Apr 10, 2007, ABsub. for C and D site classes. SHarmsen.
c rev feb 13 Art Frankel, corrects AB06 for the BC boundary
c revised june 11 2007 based on erratum in bssa for june 2007
c --- Revised Mar 27, 2007 SHarmsen, USGS. [email protected] 303 273 8567;
c This code implements nga relations; Kanno, new 2006 CEUS models, and previously used
c attenuation models.
c --- Oct 2007 rev: update B&A NGA relation, which now has 23 spectral periods (see Eq Spectra Mar 2008)
c === Jan 25 2007 rev: 3-branch gnd, to increase epistemic uncert where needed.
c --- Important, this 3-branch has been included in 4 NGA relations only as of jan 25 07
c --- These are A&B (21), C&B (22), C&Y (23), and Idriss (26). Dont use A&S2005
c
c 3-branch gnd is an all-or-none proposition as of jan 31. If last period doesnt
c ask for it, none will get it. THis is a coding error. But not fixed.
c === Oct revision: uses array constructors which replace previous "data " statements that were
c used to initialize arrays but cannot be used with some F95 compilers
c --- Dec revision: add effects of buried blind thrust and normal if present. Extra
c --- branching is implied if wtnormal or wtrev > 0. If everything is StrikeSlip, these
c --- additional steps are not done.
c --- Dec revision: use mean (rjb) instead of sample (rjb) i.e., replace spinning fault with mean
c --- distance from random-orientation fault, uses Numerical Recipes el2 routine
c ---- For details, see getmeanrjb subroutine. These mean distances are read
c ---- in from a precomputed file in hazgridXnga13l. Saves a few sec
c --- Future directions: finite reverse and normal slip faults need to be coded up.
c --- Simplest model might be a circular surface with center at each source point.
c
c --- Compile with f95 or gfortran. On Solaris, use -e for extended line length. link to iosubs.o
c Try this on sun using Solaris 10 fortran 95:
c f95 hazgridXnga13l.f -o hazgridXnga13l -O iosubs_128.o -e -ftrap=%none
c the -fast flag has been known to fail on hazgridXnga13l runs...
c Try this on PCs with gfortran:
c gfortran hazgridXnga13l.f iosubs.o -ffixed-line-length-none -static -o hazgridXnga13l.exe -finit-local-zero -ffpe-trap=
c
c the flag -finit-local-zero initializes uninitialized variables (arrays) to 0 and
c logicals to .false. This flag was brought in to veersions of gfortran after 2007, and is not a standard.
c
c --- f95 man pages say to compile with -ftrap=%none if -fast flag is used.
c ---- For listing use flag -Xlist
c --- Further notes:
c
c Feb 9 2007: can have up to 8 atten models per spectral period.
c this limit used to be 7. However, we need 8 for
c NM & Charleston Seismic Zones in the 2007 update. Steve H.
c
c --- Computation approach is that of Frankel in this version: compute
c the cumu. prob. exceedance over R and M for each atten and add the
c weight*pr to a multi-dimension matrix, which is here called pr (actually mean
c freq. of exceedance).
c These multi-dimensions correspond to distance, magnitude, ground-motions sampled,
c depth-to-top-of-rupture,
c and spectral period. Depth-to-top is a new prediction variable. Dtor requires a
c separate dimension because median motion can vary with uncertainty in Dtor.
c To some extent, atten models of 2002 had some sensitivity to dtor (not all however).
c This depth was previously fixed at 5 km for WUS gridded. 5 or 10 for CEUS gridded.
c Now, however, dtor can vary.
c This version of the gridded hazard program
c also precomputes the complementary normal prob and stores in p(), a 1-dim array
c
c -- iatten=21 boore-atkinson nga updated to the Apr 2007 version. SH. Note:
c B&A have 7.5 and 10-s SA predictions but 10s had to be modified for normal faulting.
c -- iatten=22 campbell-bozorgnia nga updated to the 3-2008 vers,
c the CB update includes peak displacement, a novelty. Sigma for
c random horizontal component is the default now.
c === iatten=32 CB 2012 for vertical dipping faults. (testing dec 20 2012)
c -- iatten=23 chiou-youngs nga vers 3-2008. (earlier 2006 version is also available)
c -- iatten=24 abrahamson-silva partially set up mar 06 (this relation will probably change).
c -- iatten=25 idriss pga oct 2005.
c -- iatten=26 Kanno et al. BSSA 2006. This model has large aleatory sigma for
c --- all spectral periods, about 50% larger than NGA relations above.
c Next,
c Older atten. relations. Some are for fixed site conditions and some for
c Vs-30 dependent site conditions. CEUS fixed site is HR or FR; WUS FR or soil.
c --- added iatten= 1 Spudich ea, 2000. From BJF93. Has siteamp from BJF97.
c --- added iatten= 2 toro ceus BC rock
c --- added iatten= -2 toro ceus hard rock
c --- added iatten= 3 Sadigh et al ( rock-site coeffs.& eqn) nov 22 2005.rock
c --- added iatten= -3 Sadigh et al (soils-site coeffs.&eqn) in prep aug06
c iatten = 4 AB06 BC Atkinson and B00re 2006 (added Nov 14 2006)
c iatten = -4 AB06 hardrock. There is a siteamp that is added to hardrock median;
c however, it is 0 (in logspace) for vs30=760.
c --- iatten == 20 :: AB06 with 200bar stress, siteamp. HR coeffs used if vs30>1000
c --- add iatten= 5 AB94 ceus (New, previously had same 6 index as fea. )
c --- add iatten= -5 AB94 HRceus (New, previously had same 6 index as fea.)
c --- added iatten= 6 Frankel ea BC rock, ceus
c --- add iatten= -6 FEA HRceus (New, previously had same 6 index as fea.)
c
c --- added iatten= 7 Somerville ceus. BCrock. Note: Somerville is used
c for the finite-fault portion of gridded hazard. Used with Charleston
c --- added iatten= -7 Somerville ceus. hardrock.
c--- added iatten= 8 Abrahamson-Silva 1997. rock. july 25 2006
c --- added iatten= 9 Campbell and Bozorgnia 2003. rock. july 25 2006
c --- added iatten= -9 Campbell and Bozorgnia 2003. D soil. future 2006
c --- added iatten= 10 Campbell CEUS BC or firmrock 2003. july 25 2006
c --- added iatten= -10 Campbell CEUS A or hardrock 2003. aug 2006
c--- added iatten= 11 BJF 1997. All Vs30 allowed, like NGA relations. Mech dependent. july 26 2006.
c --- added iatten= 12 AB intraslab seismicity Puget Sound region B or BC-rock condition.
c --- repaired rpga calc feb 12 2007. Modify c3 dtor to min(100,dtor)
c --- added iatten= -12 AB intraslab seismicity Puget Sound region C D E-soil condition
c --- added iatten= 18 AB intraslab seismicity world data BC-rock condition
c --- added iatten= -18 AB intraslab seismicity world data region C D E-soil condition
c --- added iatten= 13 Geomatrix slab seismicity rock, 1997 srl. july 25 2006
c --- added iatten= -13 Geomatrix slab seismicity soil, 1997 srl. july 25 2006
c --- added iatten= 14 Geomatrix subduction. Added July 19 2010 SH. Uses BA siteamp
c for soil Vs30 and rock Vs30 .ne. Vref (760 m/s)
c --- added iatten= 15 Silva 2002 added jan 31 2007. hr or bc only
c === iatten = 19 Tavakoli and Pezeshk 2005 added nov 14 2006.
c ---- TP05 has hardrock coef c1 used if vs30>900 m/s (no negative index though).
c - - - iatten = 27 Zhao et al for deep inslab. Can use for LP. Nov 20 2008.
c - - - iatten = 28 Zhao et al for deep interface. Can use for LP to 5 sc Tmax.
c GMPE index 28 became active July 12 2010. Zhao interface (may be used for
c America Samoa Oceanic Sources NSHMP PSHA work. To be determined later).
c --- index 29 BSSA 2013 April.
c - - - iatten = 41 Motazettian and Atkinson (changed from index 14 Jul 2010)
c - - - Motazettian and Atkinson has siteamp from BJF97
c - - - iatten = 31 BCHydro for inslab. added Nov 2012. This one has nonln siteamp.
c ---- iatten = 32 CB13 march
c ---- iatten = 33 CY13 march
c ----- iatten= 34 AS13 march
c --- iatten=35 A08'
c --- iatten=36 AB06'
c --- iatten=37 Pez11
c --- iatten=38 Idriss Apr 2013 (this one has been updated- STD. Deviation, however, has not).
c --- iatten=39 GK12 with continuous basin response.
c --- iatten = 50 Atkinson Hawaiian Islands based on BA08
c --- iatten = 51 Atkinson Induced Earthquake GMM works with Vs30 of 760 m/s or near B/C boundary
c ---
c --- SOme New Features (compared to 2002 update versions):
c --- code works for a large grid of sites or a small set (<=30) of sites.
c --- Features of hazgridXv31.f have been included here so that only
c --- one code is necessary.
c --- Input file has a new first line : number of sites, or nrec.
c --- Make nrec 0 to get the grid features of hazgridXv3.
c --- Make nrec 1 to 30 to specify set of sites <lat(i),long(i),i=1,...,nrec>
c --- Output file names are specified in input file, one per period.
c
c --- New Feature:
c --- Source box is now independent of site box and is specified same way.
c
c --- More new features...:
c---- now includes siteamp w/ some older attenuation relations i.e., variable VS30.
c --- How siteamp is handled is an open question as of July 2006. But using the
c --- simple BJF model for initial try. Many relations have nonlinear siteamp fcn
c
c -- Geotechnical input : Vs30 (which can be a geographic array), H=depth to Vs2500 m/s.
c --- To do: array of Vs30 will require a separ. Vs30 dim. on pr(...). Code is triggered
c --- to read Vs30 array if the initial input value of Vs30 is set to 0.0.
c --- Relatively low Vs30 produces nonlinear siteamp in all NGA relations (except idriss).
c --- Geotechnical information is communicated in common /geotec/
c
c --- Depth_to_top: now is a distribution with 1 to 3 depths
c --- Enter wt as two distributions, for M<6.5 and for M>=6.5:
c ntor, (dtor(k),wtor(k),wtor65(k),k=1,ntor) with 1 <= ntor <=3. wtor65 for M>=6.5
c wtor is weight distribution form M<6.5, wtor65 is weight distribution for M>=6.5
c ---We might want to model a shallower average top-of-rupture for larger events. Example
c --- 3 2. .3 .6 6. .4 .3 9. .3 .1
c --- says "3 dtors, 2 km, 6 km, and 9 km, resp. For smaller source, roughly uniform
c --- distribution of depths. For M>=6.5, shift center of distribution to shallower depths."
c --- Testing this distribution concept Mar 2006. For most sites,
c ---plausible variation of dtor seems to have little significance for Pex> o(10^-6)
c --- The distance metric is r_jb.
c --- Variable dtor is also of potential interest for intraslab seismic hazard.
c --- Dtor variability could be linked to geographic position, but not programmed.
c Assumes that r_cd = sqrt(r_jb**2+dtor**2) in all instances. SH Mar 2006.
c
c -- Additional source sense-of-slip input: wtss, wtrev, and wtnormal. Please
c dont mix reverse and normal but it is fine to mix ss+rev or ss+normal in any combination
c In 2002 style-of-slip variability was controlled in individual attenuation
c model inputs. Here it is a global variable, in common /mech/
c
c --- Older atten. model subroutines have been brought up to NGA style for
c --- standard 7 periods and BC rock site conditions. A few have variable vs30
c --- modeling capability (BJF, Spudich2000, Motazetian). CEUS models have A-rock.
c Older notes (also see Frankel codes for more notes)
c--- with clipping for Toro and new Boore tables
c--- calculates mean annual rates of exceedances for gridded a-values
c--- can use b-value and Mmax matrix
c--- this version can use finite faults centered on grid cells for M>=6.5
c--- randomizes strike of faults
c--- this version can use multiple attenuation functions at different periods
c--- choice of attenuation relation: Joyner-Boore, Toro, Sadigh, Campbell, etc.
c--- From hazgridXnga written by Art Frankel
c
c--- to run: hazgridXnga13l inputfile > log.file
c--- try: hazgridXnga13l WUSmapC.in > WUSmapC.log
c --- SDI option: hazgridXnga13l input file dy 1.0 > log file
c ---- this SDI option has been built in to more recent WUS GMPEs and to Intraplate GMPEs
c ---- SDI has not been built into CEUS GMPEs or to older WUS GMPEs as of Mar 21 2013 SH.
c
c--- output files have hazard curve for each site concatenated
c--- one output file per period (unless extra epistemic uncert, which adds .p and .m files)
c --- Messages are written to have more useful diagnostic information. Look at log file if
c you are having problems.
c--- Ground motion levels should be in units of g, except for PGV. Code checks for >12 g
c --- and stops if this is found. For PGV there is a related check, for PGV max < 20 cm/s.
c--- period=0 indicates PGA. period = -1 indicates PGV
parameter ( d2r=0.0174533,pi=3.14159265,fourpisq=39.4784175)
parameter (nlmx=20,npmx=8,nsrcmx=275000,nzonex=15)
real*8 emax/3.0/,sqrt2/1.4142135623/
real, dimension(2) :: mcut,dcut,tarray
real, dimension(50) :: a_fac
real, dimension(nzonex) :: mwmaxx !mwmaxx=the absolutely largest mag in a zone. new 4/04/2013.
c You can raise p dim & make some minor code changes to improve accuracy (currently
c p() has about 4 decimal place accuracy which is likely to be good enough)
common/prob/p(25005),plim,dp2 !table of complementary normal probab.
common/mech/wtss,wtrev,wtnormal
c there are several expressions that are sensitive to fault dip. What to do about this for
c gridded sources?
common/dipinf_90/dipang,cosDELTA,cdipsq,cyhwfac,cbhwfac
common/dipinf_50/dipang2,cosDELTA2,cdip2sq,cyhwfac2,cbhwfac2
common/geotec/vs30,dbasin
common/depth_rup/ntor,dtor,wtor,wtor65
common/gail3/freq
c rjbmean = mean rjb distance (km) from random oriented src
c Sizeof of rjbmean will work for Rmax of 1000 km. Mmin 6.05 Mmax 8.55 (feb 1 2013)
real rjbmean(26,1001)
common/ipindx/iperba(10),ipgeom(10)
common / atten / pr, xlev, nlev, icode, wt, wtdist
common/ia08/ia08
c add SDI-related common block sdi feb 22 2013
common/sdi/sdi,dy_sdi,fac_sde
real, dimension(8) :: fac_sde
real pr(310,38,20,8,3,3),xlev(20,8),wt(8,10,2),wtdist(8,10)
real ylev(20,8),m_c
c m_c is the optional corner magnitude of a tapered GR distribution.
integer nlev(8),icode(8,10)
real, dimension(3):: dtor,wtor,wtor65,w_edge
real, dimension(0:9):: dipbck
real, dimension (3,3,10) :: gnd_ep
real, dimension(16,40,8,3,3) :: hbin,ebar,rbar,mbar
real, dimension(40,nzonex):: mmax_ccd !new 4/02/2013 complementary cum. distribution of mmax
c (first dimension) in zones 1 to 10 (2nd dimension)
common/epistemic/nfi,e_wind,gnd_ep,mcut,dcut
common/deagg/deagg
c epsilon_0 pre-calcs. will occur if and only if deagg = .true.
c separate arrays are kept for each atten. model. epsilon band contribs may be computed
c later, during the summing-up stage.
common/e0_wus/e0_wus(310,31,8,3,3)
common/fix_sigma/fix_sigma,SIGMA_FX
LOGICAL fix_sigma,sdi/.false./,pnw_deep/.false./
c
c If multiple depths (3 or less) are run, the below arrays will work. However, typical
c case for WUS intrslab is just one depth at 50 km. May need to remove this dim.
c to reduce memory demand.
common/e0_sub/e0_sub(310,31,8,3)
c
c 3 intraplate models e0s will be combined in this array as a fcn of distance, M, period,
c ad depth of source-zone top.
c
c e0_ceus not saving a depth of rupture dim, not sens. to this. last dim is ip (period index)
common/e0_ceus/e0_ceus(310,31,8)
common/ceus_sig/lceus_sigma,ceus_sigma
real, dimension (107):: perbssa13
real, dimension(10) :: pertoro !added Nov 5 2015
real, dimension(40,nzonex):: mwmax,wtmw,wt_zone
common/cb13p/Percb13
common/gnome/name
integer m_ind
real, dimension(14):: a11fr,freq
c wtmj_cra = wt applied to rate when using Johnston mb to M for src in stable craton
c wtmj_ext = wt applied to M when using Johnston mb to M for src in extended margin
c wtmab_cra, wtmab_ext = ditto for Atkinson-Boore mb to M
real, dimension(40) :: wtmj_cra,wtmj_ext,wtmab_cra,wtmab_ext,wt_cra,wt_marg
real, dimension(40) :: wt_mask,xmw
c CEUS, 2008:wt_mask is assigned to one of the above depending on current site location and
c current mb --> Mw in input
c CEUS, 2013: wt_mask is assigned to one of the wtmw() vectors, the one corresponding to zone(k)
real v30(100000) !possible array of site-specific vs30 new mar 2006.
real maga,lowmag,himag,ri,minlon,maxlon,dlon,minlat,maxlat,dlat
type header_type
character*128 name(6)
real period
integer nlev
real xlev(20)
real extra(10)
end type header_type
type(header_type) :: headr,hd
type h_type
character*128 name(6)
real bval !bvalue >=0 usually
integer icum !icum=1 if cumulative rate of exceed; 0 if incremental
real atmp(20) !some compilers complain during the "gethead" step if not array
c real maga !the magnitude at which a is computed, typically 0
c real dmag !the delta-M interval used (usually 0.1)
c real lowmag !minimum M for which this rate should be used
c real cmag !the magnitude for which a recurrence interval is available
c real himag !maximum M for which this should be used
c real ri !the recurrence interval years
c real excess(14) !items without a current use.
c real s,minlon,maxlon,dlon,minlat,maxlat,dlat,en
real extra(10)
end type h_type
type(h_type) :: hd_a
real magmin,magmax,magref,sigmanf,distnf
real ymax/-100./,ymin/100./,dy/0.1/ !latitude of sites?
c real, dimension(13) :: pcut
c integer, dimension(13) :: icut
integer readn,iq_ka,iq_as,jabs,vs30_class/0/
logical finite,grid,isok/.false./,m_zones/.false./,taperGR/.false./ ! grid=.true. if stations form a regular grid
c m_zones is an indicator that magnitude zones are (are not) active
logical lceus_sigma/.false./,wus/.false./,ceus/.false./,slab/.false./
logical byeca,byesoc,byeext,byepug,hawaii/.false./,v30a,override_vs/.false./,l_mmax,hardrock,useRy0
c override_vs becomes true if for deaggregation work user inputs a vs30 on command line.
logical deagg/.false./,ss,rev,normal,obl,okabs,okgeo,okzhao,oktogo,e_wind(8)
c craton and margin determine whether a source is in craton, margin, or neither
logical, dimension(200000):: craton,margin
c oktogo is a check that period is in set of 7 available for 2003 PSHA work.
c similar logical variable for AB03 is okabs. Similar for geomatrix is okgeo.
c The pre-NGA atten. models are only called if oktogo = .true.
real slat(32),slong(32) !station coordinates if using list option
character*12 sname(32) !station names? might be useful
character*72 progname,pname*36
character*8 date,time*10,zone*5
character*80 namein,nameout,name,name3,name4,namea
c character cmnts2skip(200)*80
character*12 pithy(3),adum
dimension xlen2(50),perabs(16),perabu(10)
integer, dimension(8):: ia08
c 7/22/2010: promote hazard curve to double precision for better accumulation of small contributors
real*8, dimension (1000,20,8,3) :: prob
c the sum will be put back into single precision array called "out" prior to writing.
real, dimension(40000) :: out
integer, dimension(8,10):: iatten
integer, dimension(8,10):: irab
real, dimension(1000):: xlim,xwide
real*8 dp,pr0,prl,prlr,prr
real prd(106),camper(24),perb(23),perka(0:37),tpper(16)
real, dimension(13):: pergeo
real, dimension(22) :: NAAper
real, dimension(23) :: perId12
real, dimension(24) :: Percb13,peras13
real arat, aratemx/0.0/,Mtaper
real z1_ref, z1_refr !z1 reference values for ASK13.
real, dimension(27):: abper, abfrq
c above spectral period vectors are associated with various NGA and other
c atten models. perka corresponds to Kanno et al. added Nov 8 2006.
dimension aperiod(108),ival(8)
integer, dimension (npmx) :: nattn,iper,iperb,iperab,ipertp,
+ ipa15, isilva,ipcb13,ipbssa,icampCEUS, itoro
integer, dimension(npmx,3):: ifp
c real, dimension (16,20,10,npmx,5) :: prob5
real, dimension(0:35) :: pdgk
real, dimension (25) :: percy13 !5/2013
real, dimension(npmx) :: perx,period,safix
real, dimension(npmx+2) :: perx_pl
real, dimension(12):: perCampCEUS,per_camp,pdSilva,pda15 !add 1.5s jan 24 2014.
real, dimension(22):: perzhao
c some arrays for BSSA NGAW model
c real, dimension(5) :: dumb
c safix is for fixed-SA or fixed PGA runs, usually with deaggregation.
dimension asum(40,1000,3),arate(nsrcmx,40)
c nov 14 2007: add wtgrid matrix for use if Mtaper > 5 (agrid tapering magnitude)
real, dimension(nsrcmx) :: a,b,mmax,wtgrid,fltstrk
integer, dimension(nsrcmx) :: mzone !zone indicator
c predefined functions f and tmo:
f(u,x,y,z)=(u/x)**y *exp (-(x-u)/z) !complementary pareto distribution new 4/2012
tmo(x)=10.**(1.5*x+9.05) !moment as fcn of mag x (N-m)
c Spectral period -1 is reserved for pgv
c pcut=(/0.0062,0.0228,0.0668,0.1587,0.3085,0.5,0.6915,.8413,.9332,
c + 0.9772,0.9938,1.,1./)
c above cuts at e=-2.5,-2,-1.5,-1,-0.5,0,.5,1,1.5,2,2.5, respectively
c for possible deagg work.
c special GMM for induced seismicity, Atkinson 2015. Has PGV...
pda15=(/5.0,3.0,2.0,1.0,0.5,0.3,0.2,0.1,0.05,0.03,0.01,-1.0/) !atkinson15
c percy13 changed in the may 2013 update: they omit pgd and pgv this time.
percy13= (/ 0.0100, 0.0200, 0.0300, 0.0400, 0.0500,
1 0.0750, 0.1000, 0.1200, 0.1500, 0.1700,
1 0.2000, 0.2500, 0.3000, 0.4000, 0.5000,
1 0.7500, 1.0000, 1.5000, 2.0000, 2.5,3.0000,
1 4.0000, 5.0000, 7.5000,10.0000/)
c peras13 23 periods plus a -1 at the tail.
peras13= (/ 0.00, 0.02, 0.03, 0.05, 0.075, 0.1, 0.15, 0.2, 0.25,
1 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2., 2.5, 3., 4., 5., 6., 7.5, 10.,-1. /)
perId12 = (/0.01,0.02,0.03,0.04,0.05,0.075,0.1,0.15,0.2,0.25,0.3,0.4,0.5,
+ 0.75,1.,1.5,2.,2.5,3.,4.,5.,7.5,10./)
c Idriss 2012 GMPE periods in perId12. Add 2.5 s May 15 2015. SH.
c perabs: add 0.15 and 0.25 s Dec 9 2014. SH. Incl 0.4s june 2015
perabs = (/0.,0.2,1.0,0.1, 0.150, 0.250,0.3,0.4,0.5,0.75,1.5,2.0,2.5,3.,4.,5./)
perabu = (/0.,0.2,1.0,0.1,0.3,0.4,0.5,0.75,2.0,3./)
perx_pl= (/0.0,0.2,1.0,0.1,0.3,0.5,1.5,2.0,0.04,0.4/) !add 1.5s jan 24 2014.
c pergeo Geomatrix inslab periods available Nov 19 2008.
pergeo= (/0.,0.2,1.0,0.1,0.3,0.5,2.0,.4,0.75,1.5,3.,4.,5./)
perx= (/0.,0.2,1.0,0.1,0.3,0.5,1.5,2.0/)
perCampCEUS = (/0.01,0.2,1.0,0.1,0.3,0.4,0.5,1.5,2.0,.03,.04,.05/)
c pdgk = period set for GK12 model.
pdgk= (/0.,0.01,0.02,0.03,0.04,0.06,0.08,0.1,0.12,0.14,
& 0.16,0.18,0.20,0.22,0.24,0.27,0.30,0.33,
& 0.36,0.4,0.46,0.5,0.6,0.75,0.85,1.0,1.5,
& 2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0/)
c perb = Boore-Atkinson NGA 4/2007 period set, -1 = pgv. Now 23. 10 s is longest in April.
perb= (/-1.000, 0.000, 0.010, 0.020, 0.030, 0.050, 0.075, 0.100,
+ 0.150, 0.200, 0.250, 0.300, 0.400, 0.500, 0.750, 1.000,
+ 1.500, 2.000, 3.000, 4.000, 5.0, 7.5, 10.0/)
aperiod= (/ 0.0, -1.0, -2.0, 0.01, 0.02, 0.022, 0.025, 0.029,
1 0.03, 0.032, 0.035, 0.036, 0.04, 0.042, 0.044,
2 0.045, 0.046, 0.048, 0.05, 0.055, 0.06, 0.065,
3 0.067, 0.07, 0.075, 0.08, 0.085, 0.09, 0.095, 0.1,
4 0.11, 0.12, 0.13, 0.133, 0.14, 0.15, 0.16, 0.17,
5 0.18, 0.19, 0.2, 0.22, 0.24, 0.25, 0.26, 0.28, 0.29,
6 0.3, 0.32, 0.34, 0.35, 0.36, 0.38, 0.4, 0.42, 0.44,
7 0.45, 0.46, 0.48, 0.5, 0.55, 0.6, 0.65, 0.667, 0.7,
8 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4,
9 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.5, 2.6, 2.8,
1 3.0, 3.2, 3.4, 3.5, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8,
1 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0 /)
c modified the .3155 to .3, modified the .0996 to 0.1 sh mar 16. mod .3968 to 0.4
c june 30 2007.
perbssa13=(/-1.000000, 0.000000, 0.010000, 0.020000, 0.022000, 0.025000, 0.029000,
+ 0.030000, 0.032000, 0.035000, 0.036000, 0.040000, 0.042000, 0.044000, 0.045000, 0.046000, 0.048000,
+ 0.050000, 0.055000, 0.060000, 0.065000, 0.067000, 0.070000, 0.075000, 0.080000, 0.085000, 0.090000,
+ 0.095000, 0.100000, 0.110000, 0.120000, 0.130000, 0.133000, 0.140000, 0.150000, 0.160000, 0.170000,
+ 0.180000, 0.190000, 0.200000, 0.220000, 0.240000, 0.250000, 0.260000, 0.280000, 0.290000, 0.300000,
+ 0.320000, 0.340000, 0.350000, 0.360000, 0.380000, 0.400000, 0.420000, 0.440000, 0.450000, 0.460000,
+ 0.480000, 0.500000, 0.550000, 0.600000, 0.650000, 0.667000, 0.700000, 0.750000, 0.800000, 0.850000,
+ 0.900000, 0.950000, 1.000000, 1.100000, 1.200000, 1.300000, 1.400000, 1.500000, 1.600000, 1.700000,
+ 1.800000, 1.900000, 2.000000, 2.200000, 2.400000, 2.500000, 2.600000, 2.800000, 3.000000, 3.200000,
+ 3.400000, 3.500000, 3.600000, 3.800000, 4.000000, 4.200000, 4.400000, 4.600000, 4.800000, 5.000000,
+ 5.500000, 6.000000, 6.500000, 7.000000, 7.500000, 8.000000, 8.500000, 9.000000, 9.50,10.0/)
c PerCB13 = period set for the CB13 NGAW(2) GMM.
PerCB13=(/0.01,0.02,0.03,0.05,0.075,0.1,0.15,0.2,0.25,0.3,0.4,
+ 0.5,0.75,1.,1.5,2.,2.5, 3.,4.,5.,7.5,10.,0.,-1./)
c per_camp available spectral periods for CampCEUS (2003). PGA is 0.0 s here.
per_camp = (/0.00,0.2,1.0,0.1,0.3,0.4,0.5,1.5,2.0,.03,.04,.05/)
abper = (/5.0000, 4.0000, 3.1250, 2.5000, 2.0000, 1.5873, 1.5,1.2500, 1.0000,
1 0.7937, 0.6289, 0.5000, 0.4, 0.3, 0.2506, 0.2000, 0.1580,
1 0.1255, 0.1, 0.0791, 0.0629, 0.0499, 0.0396, 0.0315, 0.0250,
1 0.0000, -1.000/)
c BCHydro periods add nov 2012
NAAper =(/0.00, 0.050, 0.075,0.100, 0.150,0.200,0.250,0.300,0.400,
+ 0.500, 0.600, 0.750, 1.000, 1.500, 2.000, 2.500, 3.000, 4.000, 5.000, 6.0, 7.500, 10.0/)
c
c Silva Gregor and Darragh 2002 CEUS periods available as of july 1 2008.
pdSilva=(/0.,0.04,0.05,0.1,0.2,0.3,0.4,0.5,1.,1.5,2.,5./)
c Tavakoli periods 0 = pga. added 0.4 s june 30 2008 (interpolated)
tpper = (/0.00e+00,.04,5.00e-02,8.00e-02,1.00e-01,1.50e-01,2.00e-01,
1 0.3,0.40,0.5,
1 7.50e-01,1.00e+00,1.50e+00,2.00e+00,3.00e+00,4.00e+00/)
c a11fr frequencies for CEUS 2011 GMPES: (2.5 hz not available)
a11fr =(/0.20, 0.333, 0.50, 0.6667, 1.00, 2.00, 3.33, 5.00,10.00,20.,33.00,50.00,99.00,89.00/)
c available periods for CY as of jan 2009. pga=0.0, pgv is -1.0 here
prd=(/0.0,0.020,0.022,0.025,0.029,0.030,0.032,0.035,0.036,0.040,0.042,0.044,0.045,0.046,
10.048,0.050,0.055,0.060,0.065,0.067,0.070,0.075,0.080,0.085,0.090,0.095,0.100,0.110,0.120,
10.130,0.133,0.140,0.150,0.160,0.170,0.180,0.190,0.200,0.220,0.240,0.250,0.260,0.280,0.290,
10.300,0.320,0.340,0.350,0.360,0.380,0.400,0.420,0.440,0.450,0.460,0.480,0.500,0.550,0.600,
10.650,0.667,0.700,0.750,0.800,0.850,0.900,0.950,1.000,1.100,1.200,1.300,1.400,1.500,1.600,
11.700,1.800,1.900,2.000,2.200,2.400,2.500,2.600,2.800,3.000,3.200,3.400,3.500,3.600,3.800,
14.000,4.200,4.400,4.600,4.800,5.000,5.500,6.000,6.500,7.000,7.500,8.000,8.500,9.000,9.500,
110.0,-1./)
c Toro 97 periods, with 1.5s added...
pertoro = (/0.,0.2,1.0,0.1,0.3,0.5,2.0,0.04,0.4,1.5/)
c available periods for CB as of Mar 2008. pga=0.0 here. Displacement per is -2
camper=(/0.010,0.020,0.030,0.050,0.075,0.100,0.150,0.200,0.250,0.300,0.400,0.500,0.750,
+ 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 7.5,10.0, 0.0,-1.0,-2.0/)
perka =(/0.,0.05,0.06,0.07,0.08,0.09,0.10,0.11,0.12,0.13,0.15,0.17,0.20,0.22,
+0.25,0.30,0.35,0.40,0.45,0.50,0.60,0.70,0.80,0.90,1.00,1.10,1.20,
+1.30,1.50,1.70,2.00,2.20,2.50,3.00,3.50,4.00,4.50,5.00/)
c perabs: period set for ab slab-zone (deep) eqs.
c ab06 frequencies, these don't seem to be extremely close to 1/T
abfrq = (/2.00e-01,2.50e-01,3.20e-01,4.00e-01,5.00e-01,6.30e-01,0.667,8.00e-01,1.00e+00,
1 1.26e+00,1.59e+00,2.00e+00,2.52e+00,3.17e+00,3.99e+00,5.03e+00,6.33e+00,
1 7.97e+00,1.00e+01,1.26e+01,1.59e+01,2.00e+01,2.52e+01,3.18e+01,4.00e+01,
1 0.00e+00,-1.00e+00/)
c modified perzhao to have 22
PerZhao = (/0.01,0.05,0.10,0.15,0.20,0.25,0.30,0.40,
& 0.50,0.60,0.70,0.75,0.80,0.90,1.0,1.25,1.5,2.0,2.5,3.0,4.0,5.0/)
c use 2008 USGS logic tree for Mmax for the wts for mbLg from 5.05 to 7.45. We may
c read in this distribution but this is safe for envisioned purpose: web site server for
c site-specific deagg.
dipbck =(/90.0,80.,70.,60.,50.,45.,40.,35.,30.,25.0/)
c the size of below arrays was increased to 40 because of potentially greater Mmax under consideration.
wtmj_cra=(/1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,
+0.9,0.7,0.2,0.2,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./)
wtmj_ext=(/1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,
+1.,1.,1.,0.9,0.7,0.7,0.2,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./)
wtmab_cra=(/1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,
+1.,1.,0.9,0.9,0.7,0.2,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./)
wtmab_ext=(/1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,
+1.,1.,1.,1.,1.,1.,0.9,0.7,0.2,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./)
pithy = (/'Using Median','Median+EpUnc','Median-EpUnc'/)
coef= 3.14159/180.
c write(6,*) "enter name of input file"
900 format(a)
if(iargc().lt.1)stop'Usage hazgridXnga13l inputfile > logfile'
if(iargc().eq.3)then
call getarg(2,adum)
if(adum.eq.'fixsigma'.or.adum.eq.'FIXSIGMA')then
call getarg(3,adum)
read(adum,'(f6.3)')sigma_fx
fix_sigma=.true.
print *,'This run will fix sigma at ',sigma_fx,' for all WUS GMPEs all sp. periods'
elseif(adum.eq.'dy'.or.adum.eq.'dY'.or.adum.eq.'DY')then
sdi=.true.
call getarg(3,adum)
read(adum,'(f6.3)')dy_sdi
print *,'The code will compute inelastic displ spectra with dy(cm) ',dy_sdi
else
print *,'dont know what to do with ',adum
endif
elseif(iargc().gt.4)then
deagg=.true.
call getarg(2,adum)
read(adum,'(f7.4)')rlatd
ind=index(adum,'.')
jnd=index(adum,' ')-1
if(ind.eq.0)adum=adum(1:jnd)//'.'
read(adum,'(f8.4)')rlatd
call getarg(3,adum)
ind=index(adum,'.')
jnd=index(adum,' ')-1
if(ind.eq.0)adum=adum(1:jnd)//'.'
read(adum,'(f9.4)')rlond
c rlatd rlond are the coordinates of the deagg-analysis site. This location overrides the stuff
c in the input file.
write(6,*)'#HazgridXnga13l: deagg site location: ',rlatd,rlond
call getarg(4,adum)
read(adum,'(i1)')npd
do i=5,npd +4
call getarg(i,adum)
c adum could be sa(g) or pgv (cm/s). need flexi format
ind=index(adum,'.')
jnd=index(adum,' ')-1
override_vs= .false.
if(ind.eq.0)adum=adum(1:jnd)//'.'
read(adum,'(f8.4)')safix(i-4)
write(6,*)'#command line sa ',safix(i-4)
enddo
if(iargc().gt.npd+4)then
call getarg(npd+5,adum)
ind=index(adum,'.')
jnd=index(adum,' ')-1
if(ind.eq.0)adum=adum(1:jnd)//'.'
read(adum,'(f6.1)')vs30d
override_vs=.true.
write(6,*)'#vs30 will be ',vs30d,' m/s for this run. From command line'
else
override_vs=.false.
endif
rbar=0.
mbar=0.
ebar=0.
hbin=0.0
c prob5=0.
endif
call getarg(1,namein)
inquire(file=namein,exist=isok)
if(isok)then
open(unit=1,file=namein,status='old')
else
write(6,*)'File not found: ',namein
stop 'Put in working dir and retry'
endif
call date_and_time(date,time,zone,ival)
write (6,61)date,time,zone,namein
61 format('hazgridXnga13l (09/30/2016) log file. Pgm run on ',a,' at ',a,1x,a,/,
+ '# Control file:',a)
call getarg(0,progname)
ind=index(progname,' ')
if(ind.gt.30)then
pname=progname(ind-36:ind-1)
else
pname=progname(1:ind)
endif
c Initialize truncated normal array Pex, store in p().
c The indep. variable is a real*8 to reduce discretization error.
prl=0.5*(derf(emax/sqrt2) + 1.0)
prlr=1.0/prl
prr=1.0-prl
plim=-emax/sqrt2
pr0=plim
ii=1
dp=0.00004*(3.3-plim)
dp2=1.0/dp
p(1)=1.e-9
do i=2,25001
pr0=pr0+dp
p(i)=((derf(pr0)+1.)*0.5-prr)*prlr
c if(p(i).ge.pcut(ii))then
c icut(ii)=i
c ii=ii+1
c endif
enddo
p(25002)=1.0
c write(6,*)'epsilon indexes: '
c do ii=1,13
c write(6,*)ii,icut(ii)
c enddo
mcut(1)=6.0
mcut(2)=7.0
dcut(1)=10.
dcut(2)=30.
pr=0.
a_fac = 1. !truncated GR or tapered GR? If truncated, a_fac=1
c call initialize()
c End initializing the truncated normal PEx() array=p()
c Made the indep. variable a real*8 to reduce discretization error. However,
c you have to be careful that your erf can take a real*8. If not, spr0=pr0 assignment
c where spr0 is single precision should be tried.
v30a = .false.
distnf = 0.0; sigmanf = 0.0
wt_mask=1.0
craton=.false.
margin=.false.
write(6,580)'Enter a zero for grid of sites 1 to 30 for list: '
read(1,*,err=2106)nrec
write(6,*)nrec
if(nrec.eq.0)then
grid=.true.
c write(6,*) "for sites: enter min lat, max lat, dlat"
read(1,*) ymin, ymax, dy
if(.not.deagg)then
latmin=nint(ymin)
ylatmin=float(latmin)
write(6,*)'Receiver latitude range ',ymin,ymax,dy
endif
c write(6,*) "for sites: enter min lon, max lon, dlon"
read(1,*) xmin, xmax, dx
if(.not.deagg)then
write(6,*)' & Longitude range ',xmin,xmax,dx
nx= nint((xmax-xmin)/dx) +1
ny= nint((ymax-ymin)/dy) +1
nxbog=(xmax-xmin)/dx+1 !the old way, can yield bogus estimate of nx.
nybog=(ymax-ymin)/dy+1
write(6,*) nx,ny,' old calc: ',nxbog,nybog
write(6,*)'Grid_of_sites hazcurves underway'
nrec= nx*ny
endif !if not deagg
elseif(nrec.lt.33)then
grid=.false.
write(6,*)'Program will compute gridded haz at ',nrec,' sites'
dx =0.1
dy =0.1 !need defaults
do i=1,nrec
c write(6,*)'Enter station lat,long (dec deg) and name(1 word): '
580 format(a,$)
read(1,*)slat(i),slong(i),sname(i)
write(6,*)slat(i),slong(i),' ',sname(i)
ymin=min(ymin,slat(i))
ymax=max(ymax,slat(i))
if(slat(i).lt.-88.)stop 'invalid station location. Antarctica?'
enddo
latmin=nint(ymin)
ylatmin=ymin
latmax=nint(ymax)
else
write(6,*)'Code expected first line of input to be nrec'
write(6,*)'Valid nrec are 0, 1, ..., 29,30. Just read in ',nrec
stop 'Please correct input file.'
endif
if(deagg)then
c replace whatever was in the file with the command line location
nrec=1
grid=.false.
slat(1)=rlatd
ymax=rlatd; ymin=rlatd
slong(1)=rlond
sname(1)='DEAG1'
ylatmin=rlatd
latmin=nint(rlatd); latmax=nint(rlatd)
nx=1
ny=1
dmbin = 0.1 !try fine dM for initial deagg
endif !deagg is true?
byeext=index(namein,'EXTmap').gt.0
byeca=index(namein,'CAmap').gt.0.or.index(namein,'creepmap').gt.0
byesoc=index(namein,'brawmap').gt.0
byepug=index(namein,'pugetmap').gt.0
normal=byeext
c ss=.not.normal !temporary testing: no reverse or oblique.
c *** NEW 11/05 **** Enter soil Vs30 condition ******NEW*******
write(6,*)'Softrock has Vs<2500 m/s in below question'
write(6,*)"For sites, enter Vs30(m/s). Z25 will be computed internally"
read(1,*)vs30 !,dbasin
dbasin = exp(7.089 - 1.144*alog(vs30)) !New Campbell default depth Z25 may 22 2013
if(override_vs)vs30=vs30d
c use Chiou-Youngs 10-2007 default depth to 1 km/s rock. Z1 Units: m.
Z1cal = exp(-7.15/4 * log(((VS30/1000.)**4 + .57094**4)/(1.360**4 + .57094**4)))
c Norm Abrahamson's CA z1 reference (eq 18). z1_ref is in units km.
z1_ref = exp ( -7.67/4. * alog( (Vs30**4 + 610.**4)/(1360.**4+610.**4) ) ) / 1000.
z1_refr=exp ( -7.67/4. * alog( (1180.**4 + 610.**4)/(1360.**4+610.**4) )) / 1000.
c z1_refr added 8/13/2013. Z1 for hard rock. This value is .0028 km or 2.8 m
c deltaZ1=0.0 !dont know use 0. from guidance in CY doc.
z1=z1cal !CY2013 function used until we know better for wus...
z1km= z1_ref !for AS need units km. Redefined as z1_ref 8/27/2013. From P.Powers
c from B Chiou email of Apr 15 2013.
deltaZ1 = Z1cal -
1 exp(-7.15/4 *
1 log(((VS30/1000.)**4 + .57094**4)/(1.360**4 + .57094**4)))
write(6,*)' Vs30 (m/s), Z1 (m) and depth of basin (km): ',vs30,Z1,dbasin,deltaZ1
if(vs30.lt.90..and.vs30.gt.0.)then
write(6,*)'Vs30 = ',vs30,'. This looks unreasonable.'
stop'hazgridXnga13l: Please check input file just before distance incr,dmax'
elseif(vs30.eq.0.)then
write(6,580)'Enter the name of the binary vs30 array: '
read (1,900)name
inquire(file=name,exist=isok)
if(.not.isok)stop'Vs30 array not in working directory'
call openr(name)
write(6,580)'Enter ymin, ymax, dy (lat, d) of Vs30 array: '
read (1,*)vymin,vymax,vdy
write(6,580)'Enter xmin, xmax, dx (long, d) of Vs30 array: '
read (1,*)vxmin,vxmax,vdx
write(6,580)'Enter a default vs30 in case site not inbounds: '
read (1,*)vs30dflt
if(vymin.gt.vymax)then
y0=vymin
vymin=vymax
vymax=y0
c fix it if nec.
endif
nvx= nint((vxmax-vxmin)/vdx) +1
nvy= nint((vymax-vymin)/vdy) +1
c nint insures portability among different computers
v30a=.true.
nv30=nvx*nvy
c 30 characters
c nhd=308
c 128 characters
nhd=896
call gethead(hd,nhd,nread)
if(hd%extra(2).ne.vxmin)stop'mismatch vs30'
call getbuf2(v30,nv30,nvread)
write(6,*)'For vs30 expect ',nv30,' got ',nvread
write(6,*)'**** Warning on site variability: ******************'
write(6,*)'If using ChiouYoungs NGA, Z1 does not vary with Vs30'
endif
c up to three depth to top of rupture values may be enterred to define uncertainty
c in this important variable. Units km.
write(6,*)'Separate weights for dtor | M<6.5 or M>=6.5 in below question'
write(6,*)"Enter ndtor, (dtor(k),wt_65-(k),wt_65+),k=1,..ndtor<=3"
read (1,*)ntor,(dtor(k),wtor(k),wtor65(k),k=1,ntor)
write(6,*)ntor, dtor(1),dtor(2),dtor(3),' km was input'
c check reasonableness of distribution
wtor6=0.0
wtor7=0.0
tormin=10000.
do k=1,ntor
wtor6=wtor6+wtor(k)
wtor7=wtor7+wtor65(k)
tormin=min(tormin,dtor(k))
enddo
if (abs(wtor6-1.0).gt.0.001.and.tormin.lt.32.)then
write(6,*)'sum of dtor weights not equal to 1. Please check input'
stop 'and retry with improved weights'
elseif(abs(wtor6-float(ntor)).le.0.001.and.tormin.ge.32.)then
write(6,*)'For deep seismicity all focal-depth weights are 1'
pnw_deep = .true. !special study if tormin>32 and 3 depths
w_edge=wtor65 !the locations of steps are input in wtor65 in this special case
elseif(tormin.ge.32.)then
write(6,*)'For deep seismicity all focal-depth weights should be 1'
stop'Reset these to 1.0 and retry'
endif
if (abs(wtor7-1.0).gt.0.001.and.tormin.lt.32.)then
write(6,*)'sum of M6.5+ dtor Wt not equal to 1. Please check input'
stop 'and retry with improved weights'
endif
c large tormin could be associated with deep Benioff zone. For crustal
c earthquakes, this check isn't very interesting.
if(tormin.lt.0..or.(tormin.gt.202..and.xmin.lt.-100..and.ymin.gt.20.))stop'Top of rupture distribution