C ***** New code for AG Subduction model **** c ---------------------------------------------------------------------- subroutine AG2020 ( mag, evType, rRup, vs30, z25, ztor, region, mu, sigma, phi, tau, 1 rockPGA, specT, period1, iFlag, ACadjfac, MSAS_flag, 2 regName1, sigmaSS, phiSS, phiS2S, iBasin) C Model Version: December 11, 2020 C PGA (T=0.0) set equal to T=0.01 coefficients implicit none integer MAXPER, NPER, ACadjfac real sigma, sigmaSS parameter (MAXPER=25) character*80 RegName(10), RegName1 integer i, j, k, region, i9, iParam, count1, count2, iflag, MSAS_flag real mag, rRup, ZTOR, evType real rockPGA, period1, period(MAXPER), specT real mu, r real VLin(MAXPER), b_soil(MAXPER), c1, vs real vs30, sum1, angle1, tapera_T real n, c, T1, T2 real vs1, z25, Z25ref, ln_z25_prime c integer iprint real c1_slab(8), c1_inter(MAXPER), depthLimit, c4 real term1, term1a, term1b, term1c, term2, term3, term4, 1 term5, term6, term7, term6a real part(45), NL_soil, mu1 real a1, a2(MAXPER), a3, a4, a5, 1 a6(MAXPER), a7(MAXPER), a8(MAXPER), a9, a10(MAXPER), 1 a11(MAXPER), a12(MAXPER), a13(MAXPER), a14(MAXPER), a15, 1 a16(MAXPER), a17, a18(MAXPER), a19, a20(MAXPER), 1 a21(MAXPER), a22(MAXPER), a23(MAXPER), a24(MAXPER), a25(MAXPER), 1 a26(MAXPER), a27(MAXPER), a28(MAXPER), a29(MAXPER), a30(MAXPER), 1 a31(MAXPER), a32(MAXPER), a33(MAXPER), a34(MAXPER), a35(MAXPER), 1 a36(MAXPER), a37(MAXPER), a38, a39(MAXPER), a40, 1 a41(MAXPER), a42, a43, a44, a45 real AKFac(MAXPER), CasFac(MAXPER), AKFacT, CasFacT real c1_inter_T, vLin_T, b_soil_T, temp1 real d1(MAXPER), d2(MAXPER), d1_T, d2_T real alpha1, alpha2, f2, f3, phi, tau, phi1, phi2, phi3 real d0, d3, d4, d5, d6 real a1Global(MAXPER), a1G_T real T1_phi2, T2_Phi2, T3_phi2, T4_phi2, d3_phi2, d4_phi2, d5_phi2 real T1_phi3, T2_Phi3, T3_phi3, T4_phi3, d3_phi3, d4_phi3, d5_phi3 real alpha_phi2, alpha_phi3, phi1_SQ, A_phi2 , A_phi3 real phi_LIN, phi_LIN_PGA, tau_LIN, tau_LIN_PGA real phi_B, phi_B_PGA, phi_amp, partial_f_PGA, phiSQ_NL, tauSQ_NL, PGA_1000 real rho_W_T, rho_B_T, rho_W(MAXPER), rho_B(MAXPER) real phi1_SQ_PGA, F2_PGA, f3_PGA integer iPer real f_mag, f_site, f_slab, f_ztor, f_basin, eq_3_1 real theta(45) real phi_S2S, phi1_SQ_100, phi3_SQ_100, A_phi2_100 real phiSS, phiSS_B, phiSS_B_PGA, phi1_sq_PGA_100, 1 phiSS_LIN, phiSS_LIN_PGA, phiSS_SQ_NL, phi_S2S_PGA real phiS2S_G1 (MAXPER), phiS2S_G2 (MAXPER), phiS2S_G3 (MAXPER) real phiS2S, phiS2S_G1T, phiS2S_G2T, phiS2S_G3T real f2_phi2, f2_phi2_PGA, phi2_SQ, phi2_SQ_PGA real f2_phi3, f2_phi3_PGA, phi3_SQ, phi3_SQ_PGA integer iBasin c slab c1: Alaska, CAS, CenAm, Japan, NewZealand, SouthAm, Taiwan, global data c1_slab / 7.9, 7.1, 7.4, 7.6, 8.0, 7.5, 7.7, 7.5 / data RegName / 'Alaska', 'Cascadia', 'Central_America', 'Japan', 1 'New_Zealand', 'South_America', 'Taiwan', 'Global', 'Alaska-adjusted', 2 'Cascadia-Adjusted' / C Updated Coefficients 12/11/20 data period / 0.0, 0.01, 0.02, 0.03, 0.05, 0.075, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 1 0.5, 0.6, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.5, 10.0 / data vLin / 865.1, 865.1, 865.1, 907.8, 1053.5, 1085.7, 1032.5, 877.6, 748.2, 654.3, 587.1, 1 503.0, 456.6, 430.3, 410.5, 400.0, 400.0, 400.0, 400.0, 400.0, 400.0, 400.0, 2 400.0, 400.0, 400.0 / data b_soil / -1.186, -1.186, -1.219, -1.273, -1.346, -1.471, -1.624, -1.931, -2.188, 1 -2.381, -2.518, -2.657, -2.669, -2.599, -2.401, -1.955, -1.025, -0.299, 2 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 / data c1_inter / 8.2, 8.2, 8.2, 8.2, 8.2, 8.2, 8.2, 8.2, 8.2, 8.2, 8.2, 8.2, 8.2, 8.2, 1 8.15, 8.1, 8.05, 8.0, 7.95, 7.9, 7.85, 7.8, 7.8, 7.8, 7.8 / data a1Global / 4.596, 4.596, 4.678, 4.773, 5.029, 5.334, 5.455, 1 5.376, 4.936, 4.636, 4.423, 4.124, 3.838, 2 3.562, 3.152, 2.544, 1.636, 1.076, 0.658, 3 0.424, 0.093, -0.145, -0.32, -0.556, -0.860 / data a2 / -1.45, -1.45, -1.45, -1.45, -1.45, -1.45, -1.45, -1.425, -1.335, -1.275, 1 -1.231, -1.165, -1.115, -1.071, -1.02, -0.95, -0.86, -0.82, -0.798, -0.793, 2 -0.793, -0.793, -0.793, -0.793, -0.793 / data a6 / -0.0043, -0.0043, -0.0043, -0.0044, -0.0046, -0.0047, -0.0048, -0.0047, 1 -0.0045, -0.0043, -0.0042, -0.004, -0.0037, -0.0035, -0.0032, -0.0029, 2 -0.0026, -0.0024, -0.0022, -0.0021, -0.002, -0.002, -0.002, -0.002, -0.002 / data a7 / 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 1 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.13, 2.985, 2.818, 2.682, 2.515, 2.3 / data a8 / 0.044, 0.044, 0.044, 0.044, 0.044, 0.044, 0.044, 0.044, 0.043, 0.042, 0.041, 1 0.04, 0.039, 0.038, 0.037, 0.035, 0.034, 0.032, 0.031, 0.03, 0.029, 0.028, 2 0.027, 0.026, 0.025 / data a10 / 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 1 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.21, 3.13, 2.985, 2.818, 2.682, 2.515, 2.3 / data a11 / 0.007, 0.007, 0.007, 0.007, 0.007, 0.007, 0.007, 0.007, 0.0062, 0.0056, 1 0.0051, 0.0043, 0.0037, 0.0033, 0.0027, 0.0019, 0.0008, 0.0, 0.0, 2 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 / data a12 / 0.9, 0.9, 1.008, 1.127, 1.333, 1.565, 1.679, 1.853, 2.022, 2.181, 2.281, 1 2.379, 2.339, 2.217, 1.946, 1.416, 0.394, -0.417, -0.725, -0.695, -0.638, 2 -0.597, -0.561, -0.53, -0.486 / data a13 / 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.002, -0.007, -0.011, 1 -0.015, -0.021, -0.028, -0.041, -0.05, -0.057, -0.065, -0.077, -0.088, 2 -0.098, -0.11, -0.127 / data a14 / -0.46, -0.46, -0.46, -0.46, -0.46, -0.46, -0.46, -0.46, -0.46, -0.46, 1 -0.46, -0.47, -0.48, -0.49, -0.50, -0.51, -0.52, -0.53, -0.54, -0.54, 2 -0.54, -0.54, -0.54, -0.54, -0.54 / data a16 / 0.090, 0.090, 0.090, 0.090, 0.090, 0.090, 0.090, 0.090, 0.084, 0.080, 1 0.078, 0.075, 0.072, 0.070, 0.067, 0.063, 0.059, 0.059, 0.060, 0.059, 2 0.050, 0.043, 0.038, 0.032, 0.024 / data a18 / -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.2, -0.186, -0.15, -0.14, -0.12, 1 -0.1, -0.08, -0.06, -0.047, -0.035, -0.018, -0.01, -0.005, 0.0, 0.0, 2 0.0, 0.0, 0.0, 0.0 / data a20 / 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.055, -0.105, -0.134, -0.15, 1 -0.15, -0.15, -0.15, -0.15, -0.15, -0.13, -0.11, -0.095, -0.085, 2 -0.073, -0.065, -0.06, -0.055, -0.045 / data a21 / 0.04, 0.04, 0.04, 0.04, 0.04, 0.06, 0.1, 0.135, 0.17, 0.17, 0.17, 0.17, 1 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17 / data a22 / 0.04, 0.04, 0.04, 0.04, 0.04, 0.06, 0.1, 0.135, 0.17, 0.17, 0.17, 0.17, 1 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17, 0.17 / data a23 / 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.069, 0.14, 0.164, 0.19, 0.206, 0.22, 1 0.225, 0.217, 0.185, 0.083, 0.045, 0.026, 0.035, 0.053, 0.072, 2 0.086, 0.115, 0.151 / data a24 / 0.0015, 0.0015, 0.0015, 0.0015, 0.0011, 0.0011, 0.0012, 0.0013, 0.0013, 1 0.0013, 0.0014, 0.0015, 0.0015, 0.0015, 0.0014, 0.0013, 0.0014, 0.0015, 2 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014 / data a25 / 0.0007, 0.0007, 0.0006, 0.0006, 0.0006, 0.0004, 0.0003, -0.0002, -0.0007, 1 -0.0009, -0.001, -0.001, -0.0011, -0.0012, -0.0011, -0.0008, -0.0004, 0.0002, 2 0.0004, 0.0007, 0.001, 0.0013, 0.0015, 0.0017, 0.0017 / data a26 / 0.0036, 0.0036, 0.0036, 0.0037, 0.0039, 0.0039, 0.0039, 0.0037, 0.0031, 1 0.0027, 0.002, 0.0013, 0.0009, 0.0006, 0.0003, 0.0001, -0.0001, 0.0, 2 0.0, 0.0003, 0.0007, 0.0014, 0.0015, 0.0015, 0.0015 / data a27 / -0.0004, -0.0004, -0.0005, -0.0007, -0.0009, -0.0009, -0.0008, -0.0009, 1 -0.001, -0.0011, -0.0009, -0.0007, -0.0007, -0.0007, -0.0007, -0.0008, 2 -0.0008, -0.0007, -0.0007, -0.0007, -0.0006, -0.0004, -0.0003, -0.0002, -0.0001 / data a28 / 0.0025, 0.0025, 0.0025, 0.0025, 0.0026, 0.0026, 0.0026, 0.0022, 0.0018, 1 0.00155, 0.0014, 0.0011, 0.0008, 0.0006, 0.0004, 0.0002, 0.0001, 0.0002, 2 0.0002, 0.0004, 0.0006, 0.0008, 0.0011, 0.0014, 0.0017 / data a29 / 0.0006, 0.0006, 0.0005, 0.0005, 0.0004, 0.0003, 0.0003, 0.0001, -0.0001, 1 -0.0003, -0.0002, 0.0, 0.0002, 0.0002, 0.0002, 0.0001, 0.0, 0.0, -0.00015, 2 -0.0002, -0.0002, -0.0001, 0.0, 0.0001, 0.0002 / data a30 / 0.0033, 0.0033, 0.0033, 0.0034, 0.0036, 0.0037, 0.0038, 0.0037, 0.0035, 1 0.0033, 0.0032, 0.003, 0.0027, 0.0025, 0.0022, 0.0019, 0.0016, 0.0014, 2 0.0012, 0.0011, 0.001, 0.001, 0.001, 0.001, 0.001 / data a31 / 3.7783, 3.7783, 3.8281, 3.8933, 4.2867, 4.5940, 4.7077, 4.6065, 1 4.1866, 3.8515, 3.5783, 3.2493, 2.9818, 2.7784, 2.4780, 1.9252, 2 0.9924, 0.4676, 0.0579, -0.1391, -0.3030, -0.4094, -0.5010, 3 -0.6209, -0.6221 / data a32 / 3.3468, 3.3468, 3.4401, 3.5087, 3.6553, 3.9799, 4.1312, 4.2737, 1 3.9650, 3.6821, 3.5415, 3.3256, 3.1334, 2.9215, 2.5380, 1.9626, 2 1.3568, 0.8180, 0.4389, 0.1046, -0.1597, -0.2063, -0.3223, 3 -0.4223, -0.5909 / data a33 / 3.8025, 3.8025, 3.9053, 4.0189, 4.2952, 4.5464, 4.6138, 4.5290, 1 4.1656, 3.9147, 3.7846, 3.5702, 3.3552, 3.0922, 2.6572, 2.1459, 2 1.3499, 0.8148, 0.3979, 0.1046, -0.2324, -0.5722, -0.8631, -1.1773, -1.4070 / data a34 / 5.0361, 5.0361, 5.1375, 5.2699, 5.6157, 6.0204, 6.1625, 5.9614, 1 5.3920, 5.0117, 4.7057, 4.2896, 3.9322, 3.6149, 3.1785, 2.5722, 2 1.6499, 1.0658, 0.6310, 0.3882, 0.0164, -0.2802, -0.4822, -0.7566, -1.0870 / data a35 / 4.6272, 4.6272, 4.6958, 4.7809, 5.0211, 5.3474, 5.5065, 5.5180, 5.1668, 1 4.8744, 4.6544, 4.3660, 4.0779, 3.8146, 3.4391, 2.8056, 1.8546, 1.3020, 2 0.8017, 0.5958, 0.3522, 0.1874, -0.1243, -0.3316, -0.6783 / data a36 / 4.8044, 4.8044, 4.8943, 5.0028, 5.2819, 5.6123, 5.7668, 5.7313, 5.2943, 1 5.0058, 4.7588, 4.3789, 4.0394, 3.7366, 3.2930, 2.6475, 1.6842, 1.1002, 2 0.6737, 0.4126, 0.0097, -0.2715, -0.4591, -0.6822, -0.9173 / data a37 / 3.5669, 3.5669, 3.6425, 3.7063, 3.9184, 4.2207, 4.3536, 4.3664, 1 4.0169, 3.7590, 3.5914, 3.3704, 3.1564, 2.9584, 2.6556, 2.0667, 2 1.3316, 0.7607, 0.3648, 0.1688, -0.0323, -0.1516, -0.2217, 3 -0.3338, -0.5441 / data a39 / 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.101, 0.184, 0.315, 0.416, 1 0.499, 0.600, 0.731, 0.748, 0.761, 0.770, 0.778, 0.790, 0.799, 1 0.807, 0.817, 0.829 / data a41 / -0.029, -0.029, -0.024, -0.034, -0.061, -0.076, -0.049, -0.026, -0.011, -0.009, 1 0.005, 0.040, 0.097, 0.145, 0.197, 0.269, 0.347, 0.384, 0.397, 0.404, 2 0.397, 0.378, 0.358, 0.333, 0.281 / data d1 / 0.325, 0.325, 0.325, 0.325, 0.325, 0.325, 0.325, 0.325, 0.325, 0.325, 1 0.325, 0.325, 0.325, 0.325, 0.325, 0.325, 0.312, 0.302, 0.295, 0.289, 2 0.280, 0.273, 0.267, 0.259, 0.250 / data d2 / 0.137, 0.137, 0.137, 0.137, 0.137, 0.137, 0.137, 0.137, 0.137, 0.137, 1 0.137, 0.137, 0.137, 0.137, 0.137, 0.137, 0.113, 0.096, 0.082, 0.072, 2 0.055, 0.041, 0.030, 0.017, 0.000 / c corrected AK and CAS factors 07/20/21 data AKfac / 0.487, 0.487, 0.519, 0.543, 0.435, 0.410, 0.397, 0.428, 0.442, 0.494, 2 0.565, 0.625, 0.634, 0.581, 0.497, 0.469, 0.509, 0.478, 0.492, 0.470, 4 0.336, 0.228, 0.151, 0.051, -0.251 / data CasFac / 0.828, 0.828, 0.825, 0.834, 0.895, 0.863, 0.842, 0.737, 0.746, 0.796, 1 0.782, 0.768, 0.728, 0.701, 0.685, 0.642, 0.325, 0.257, 0.211, 0.296, 2 0.232, 0.034, -0.037, -0.178, -0.313 / data rho_W / 1.00, 1.00, 0.99, 0.99, 0.97, 0.95, 0.92, 0.90, 0.87, 1 0.84, 0.82, 0.74, 0.66, 0.59, 0.50, 0.41, 0.33, 0.30, 2 0.27, 0.25, 0.22, 0.19, 0.17, 0.14, 0.10 / data rho_B / 1.00, 1.00, 0.99, 0.99, 0.985, 0.98, 0.97, 0.96, 0.94, 1 0.93, 0.91, 0.86, 0.80, 0.78, 0.73, 0.69, 0.62, 0.56, 2 0.52, 0.495, 0.43, 0.40, 0.37, 0.32, 0.28 / data PhiS2S_G1/ 0.396, 0.396, 0.396, 0.396, 0.396, 0.396, 0.396, 0.396, 0.396, 1 0.396, 0.396, 0.396, 0.396, 0.396, 0.396, 0.396, 0.379, 2 0.366, 0.356, 0.348, 0.335, 0.324, 0.314, 0.301, 0.286 / data phiS2S_G2/ 0.396, 0.396, 0.396, 0.396, 0.467, 0.516, 0.516, 0.516, 0.516, 1 0.501, 0.488, 0.468, 0.451, 0.438, 0.420, 0.396, 0.379, 2 0.366, 0.356, 0.348, 0.335, 0.324, 0.314, 0.301, 0.286 / data phiS2S_G3/ 0.545, 0.545, 0.545, 0.545, 0.644, 0.713, 0.713, 0.647, 0.596, 1 0.539, 0.488, 0.468, 0.451, 0.438, 0.420, 0.396, 0.379, 2 0.366, 0.356, 0.348, 0.335, 0.324, 0.314, 0.301, 0.286 / c set period-independent terms depthLimit = 200. c4 = 10. n = 1.18 c = 1.88 theta(1) = 0. theta(3) = 0.1 theta(4) = 0.73 theta(5) = 0. theta(9) = 0.4 theta(15) = -0.1 theta(17) = 0. theta(19) = 0. theta(38) = 0. theta(40) = 0. theta(42) = 0. theta(43) = 0. theta(44) = 0. theta(45) = 0.34 nper = 25 C First check for the PGA case if (specT .eq. 0.0) then write (*,'( 2x,''using Pga'')') vlin_T = vLin(1) b_soil_T = b_soil(1) c1 = c1_inter(1) c1_inter_T = c1_inter(1) if ( region .eq. 8 ) then theta(1) = a1Global(1) else theta(1) = 0. endif theta(2) = a2(1) theta(6) = a6(1) theta(7) = a7(1) theta(8) = a8(1) theta(10) = a10(1) theta(11) = a11(1) theta(12) = a12(1) theta(13) = a13(1) theta(14) = a14(1) theta(16) = a16(1) theta(18) = a18(1) theta(20) = a20(1) theta(21) = a21(1) theta(22) = a22(1) theta(23) = a23(1) theta(24) = a24(1) theta(25) = a25(1) theta(26) = a26(1) theta(27) = a27(1) theta(28) = a28(1) theta(29) = a29(1) theta(30) = a30(1) theta(31) = a31(1) theta(32) = a32(1) theta(33) = a33(1) theta(34) = a34(1) theta(35) = a35(1) theta(36) = a36(1) theta(37) = a37(1) theta(39) = a39(1) theta(41) = a41(1) d1_T = d1(1) d2_T = d2(1) AKfacT = AKFac(1) CasfacT = CasFac(1) rho_B_T = rho_B(1) rho_W_T = rho_W(1) phiS2S_G1T = phiS2S_G1(1) phiS2S_G2T = phiS2S_G2(1) phiS2S_G3T = phiS2S_G3(1) goto 1011 endif C For other periods, loop over the spectral period range of the attenuation relationship. do i=2,nper-1 if (specT .ge. period(i) .and. specT .le. period(i+1) ) then count1 = i count2 = i+1 goto 1020 endif enddo C Selected spectral period is outside range defined by attenuaton model. write (*,*) write (*,*) 'AG (2020) Horizontal' write (*,*) 'attenuation model is not defined for a ' write (*,*) ' spectral period of: ' write (*,'(a10,f10.5)') ' Period = ',specT write (*,*) 'This spectral period is outside the defined' write (*,*) 'period range in the code or beyond the range' write (*,*) 'of spectral periods for interpolation.' write (*,*) 'Please check the input file.' write (*,*) stop 99 C Interpolate the coefficients for the requested spectral period. 1020 T1 = period(count1) T2 = period(count2) call s24_interp (T1,T2,b_soil(count1),b_soil(count2), + specT,b_soil_T,iflag) call s24_interp (T1,T2,vlin(count1),vlin(count2), + specT,vlin_T,iflag) call s24_interp (T1,T2,c1_inter(count1),c1_inter(count2), + specT,c1_inter_T,iflag) if ( region .eq. 8 ) then call S24_interp (T1,T2,a1Global(count1),a1Global(count2),specT,theta(1),iflag) else theta(1) = 0. endif call S24_interp (T1,T2,a2(count1),a2(count2), specT,theta(2),iflag) call S24_interp (T1,T2,a6(count1),a6(count2), specT,theta(6),iflag) call S24_interp (T1,T2,a7(count1),a7(count2), specT,theta(7),iflag) call S24_interp (T1,T2,a8(count1),a8(count2), specT,theta(8),iflag) call S24_interp (T1,T2,a10(count1),a10(count2), specT,theta(10),iflag) call S24_interp (T1,T2,a11(count1),a11(count2), specT,theta(11),iflag) call S24_interp (T1,T2,a12(count1),a12(count2), specT,theta(12),iflag) call S24_interp (T1,T2,a13(count1),a13(count2), specT,theta(13),iflag) call S24_interp (T1,T2,a14(count1),a14(count2), specT,theta(14),iflag) call S24_interp (T1,T2,a16(count1),a16(count2), specT,theta(16),iflag) call S24_interp (T1,T2,a18(count1),a18(count2), specT,theta(18),iflag) call S24_interp (T1,T2,a20(count1),a20(count2), specT,theta(20),iflag) call S24_interp (T1,T2,a21(count1),a21(count2), specT,theta(21),iflag) call S24_interp (T1,T2,a22(count1),a22(count2), specT,theta(22),iflag) call S24_interp (T1,T2,a23(count1),a23(count2), specT,theta(23),iflag) call S24_interp (T1,T2,a24(count1),a24(count2), specT,theta(24),iflag) call S24_interp (T1,T2,a25(count1),a25(count2), specT,theta(25),iflag) call S24_interp (T1,T2,a26(count1),a26(count2), specT,theta(26),iflag) call S24_interp (T1,T2,a27(count1),a27(count2), specT,theta(27),iflag) call S24_interp (T1,T2,a28(count1),a28(count2), specT,theta(28),iflag) call S24_interp (T1,T2,a29(count1),a29(count2), specT,theta(29),iflag) call S24_interp (T1,T2,a30(count1),a30(count2), specT,theta(30),iflag) call S24_interp (T1,T2,a31(count1),a31(count2), specT,theta(31),iflag) call S24_interp (T1,T2,a32(count1),a32(count2), specT,theta(32),iflag) call S24_interp (T1,T2,a33(count1),a33(count2), specT,theta(33),iflag) call S24_interp (T1,T2,a34(count1),a34(count2), specT,theta(34),iflag) call S24_interp (T1,T2,a35(count1),a35(count2), specT,theta(35),iflag) call S24_interp (T1,T2,a36(count1),a36(count2), specT,theta(36),iflag) call S24_interp (T1,T2,a37(count1),a37(count2), specT,theta(37),iflag) call S24_interp (T1,T2,a39(count1),a39(count2), specT,theta(39),iflag) call S24_interp (T1,T2,a41(count1),a41(count2), specT,theta(41),iflag) call S24_interp (T1,T2,d1(count1),d1(count2), specT,d1_T,iflag) call S24_interp (T1,T2,d2(count1),d2(count2), specT,d2_T,iflag) call S24_interp (T1,T2,rho_W(count1),rho_W(count2), specT,rho_W_T,iflag) call S24_interp (T1,T2,rho_B(count1),rho_B(count2), specT,rho_B_T,iflag) call S24_interp (T1,T2,Casfac(count1),Casfac(count2), specT,CASfacT,iflag) call S24_interp (T1,T2,AKfac(count1),AKfac(count2), specT,AKfacT,iflag) call S24_interp (T1,T2,phiS2S_G1(count1),phiS2S_G1(count2), specT,phiS2S_G1T,iflag) call S24_interp (T1,T2,phiS2S_G2(count1),phiS2S_G2(count2), specT,phiS2S_G2T,iflag) call S24_interp (T1,T2,phiS2S_G3(count1),phiS2S_G3(count2), specT,phiS2S_G3T,iflag) 1011 period1 = specT c This code is taken from the regression program, so it shows what was run c Set c1_slab if (evType .eq. 0. ) then c1 = c1_inter_T else c1 = c1_slab(region) endif r = rrup + c4*exp((mag-6.)*theta(9)) c Set break in slope of soil term VS1 = 1000. if ( vs30 .gt. VS1 ) then vs = VS1 else vs = vs30 endif c set fixed non-linear soil if ( vs .lt. vLin_T) then nl_soil = ( - b_soil_T*alog(c+rockPGA) 1 + b_soil_T*alog(rockPga+c*(vs/vLin_T)**(n)) ) else nl_soil = (b_soil_T*n*alog(vs/vLIN_T) ) endif c Set partial derivatives with respect to coefficients part(1) = 1. part(2) = alog(r) part(3) = (mag - 7.) * alog(r) if ( mag .le. c1 ) then part(4) = mag - c1 + (c1_slab(region)-7.5) * evType part(5) = 0. else part(4) = (c1_slab(region)-7.5) * evType part(5) = mag - c1 endif part(6) = rRup part(7) = 0. part(8) = 0. part(9) = 0. if ( region .eq. 3 .or. region .eq. 4 .or. region .eq. 6 ) then part(7) = 0. part(10) = evType else part(7) = evType part(10) = 0. endif if ( ztor .lt. depthLimit ) then if ( ztor .lt. 50. ) then part(8) = (ztor-50.) * evType part(11) = 0. else part(8) = 0. part(11) = evType * (ztor - 50.) endif else part(11) = evType * (depthLimit - 50.) endif part(12) = alog(vs/vLin_T) part(13) = (10. - mag )**2 part(14) = alog(r) * evType c class1 / class2 flag part(15) = float(MSAS_flag) if ( region .eq. 7 ) then part(16) = alog(r) else part(16) = 0. endif c Initialize regional terms do i9=17,45 part(i9) = 0. enddo c Region-dependent VS30 terms if ( region .eq. 1 ) then part(17) = alog(vs/vLin_T) elseif ( region .eq. 2 ) then part(18) = alog(vs/vLin_T) elseif ( region .eq. 3 ) then part(19) = alog(vs/vLin_T) elseif ( region .eq. 4 ) then part(20) = alog(vs/vLin_T) elseif ( region .eq. 5 ) then part(21) = alog(vs/vLin_T) elseif ( region .eq. 6 ) then part(22) = alog(vs/vLin_T) elseif ( region .eq. 7 ) then part(23) = alog(vs/vLin_T) endif c Region-dependent R terms (interface) if ( region .eq. 1 ) then part(24) = rRup elseif ( region .eq. 2 ) then part(25) = rRup elseif ( region .eq. 3 ) then part(26) = rRup elseif ( region .eq. 4 ) then part(27) = rRup elseif ( region .eq. 5 ) then part(28) = rRup elseif ( region .eq. 6 ) then part(29) = rRup elseif ( region .eq. 7 ) then part(30) = rRup endif c Region-dependent Constant terms if ( region .eq. 1 ) then part(31) = 1. elseif ( region .eq. 2 ) then part(32) = 1. elseif ( region .eq. 3 ) then part(33) = 1. elseif ( region .eq. 4 ) then part(34) = 1. elseif ( region .eq. 5 ) then part(35) = 1. elseif ( region .eq. 6 ) then part(36) = 1. elseif ( region .eq. 7 ) then part(37) = 1. endif c set z2.5_ref (meters) c Calc reference z25 for given vs30 for Cascadia(Region=2) or Japan (Region=4) if ( Region .eq. 2 ) then call z25_interp ( vs30, 200., 570., 8.52, 7.6, z25ref ) elseif ( Region .eq. 4 ) then call z25_interp ( vs30, 170., 800., 7.3, 4.1, z25ref ) else z25Ref = -99. endif c Check if basin is included c (Basin depth is not included for the PGARock calculation) if ( iBasin .eq. 0 ) then z25ref = -99. endif c Z2.5 scaling c Note: z2.5 passed into AG20 is in km, but the model uses meters if ( z25 .ge. 0. .and. z25ref .ge. 0.) then ln_z25_prime = alog ( ( z25*1000.0 + 50.) /(z25ref+50.) ) if ( z25 .gt. 0. ) then if ( region .eq. 1 ) then part(38) = ln_z25_prime elseif ( region .eq. 2 ) then if ( ln_z25_prime .gt. 0. ) then part(39) = ln_z25_prime-0. else part(39) = 0. endif elseif ( region .eq. 3 ) then part(40) = ln_z25_prime elseif ( region .eq. 4 ) then if ( ln_z25_prime .gt. -2. ) then part(41) = ln_z25_prime else part(41) = -2. endif elseif ( region .eq. 5 ) then part(42) = ln_z25_prime elseif ( region .eq. 6 ) then part(43) = ln_z25_prime elseif ( region .eq. 7 ) then part(44) = ln_z25_prime endif endif else ln_z25_prime = 0. endif c Add mag scaling difference for interface and slab part(45) = (c1_slab(region)-7.5) * evType c if ( evType .eq. 1 .and. mag .le. c1 ) then if ( evType .eq. 1 .and. mag .le. c1) then part(45) = mag - c1 + part(45) endif c compute median mu = NL_soil do iParam=1,45 mu = mu + theta(iParam) * part(iParam) enddo C Apply Alaska or Cascadia adjustment if requested. if (Region .eq. 1 ) then mu = mu + ACadjfac * AKfacT elseif (Region .eq. 2 ) then mu = mu + ACadjfac * CASfacT endif c Set regions name regName1 = regName(region) if ( region .eq. 1 .and. ACadjfac .eq. 1 ) then regName1 = regName(9) elseif ( region .eq. 2 .and. ACadjfac .eq. 1 ) then regName1 = regName(10) endif c ----------- This part of the code is to allow checking with the equations as written in PEER report and spreadsheet ----------- C remove this for hazard runs to speed up calculations c eq 3.3 if ( mag .lt. c1 ) then f_mag = (theta(4)+evType*theta(45))*(mag-c1) + theta(13)*(10-mag)**2 else f_mag = theta(5)*(mag-c1)+theta(13)*(10.-mag)**2 endif c eq 3.5 f_slab = theta(10)+(theta(4)+theta(45))*(C1-7.5) + theta(14)*alog(r) c eq 3.6 if ( ztor .lt. 50. ) then f_ztor = theta(8)*(ztor-50) * evtype elseif ( ztor .lt. 200. ) then f_ztor = theta(11)*(ztor-50.)*evtype else f_ztor = theta(11)*150. * evtype endif c Add regional terms if ( region .eq. 1 ) then theta(1) = theta(1) + theta(31) theta(6) = theta(6) + theta(24) theta(12) = theta(12) + theta(17) elseif ( region .eq. 2 ) then theta(1) = theta(1) + theta(32) theta(6) = theta(6) + theta(25) theta(12) = theta(12) + theta(18) elseif ( region .eq. 3 ) then theta(1) = theta(1) + theta(33) theta(6) = theta(6) + theta(26) theta(12) = theta(12) + theta(19) elseif ( region .eq. 4 ) then theta(1) = theta(1) + theta(34) theta(6) = theta(6) + theta(27) theta(12) = theta(12) + theta(20) elseif ( region .eq. 5 ) then theta(1) = theta(1) + theta(35) theta(6) = theta(6) + theta(28) theta(12) = theta(12) + theta(21) elseif ( region .eq. 6 ) then theta(1) = theta(1) + theta(36) theta(6) = theta(6) + theta(29) theta(12) = theta(12) + theta(22) elseif ( region .eq. 7 ) then theta(1) = theta(1) + theta(37) theta(6) = theta(6) + theta(30) theta(12) = theta(12) + theta(23) theta(2) = theta(2) + theta(16) endif c eq 3.7 f_site = nl_soil + theta(12)*alog(Vs/vlin_T) c Add adjustment for Alaska or Cascadia if ( region .eq. 1 ) then theta(1) = theta(1) + ACadjfac * AKfacT elseif (region .eq. 2 ) then theta(1) = theta(1) + ACadjfac * CasfacT endif c Basin terms, initialize to zero f_basin = 0. c Basin terms for Cascadia c eq 3.11 if ( region .eq. 2 .and. ln_z25_prime .gt. 0. ) then f_basin = theta(39)*ln_z25_prime endif c Basin terms for Japam c eq 3.10 if ( region .eq. 4 ) then if ( ln_z25_prime .gt. -2. ) then f_basin = theta(41)*ln_z25_prime else f_basin = theta(41)*(-2.) endif endif c Eq 3.1 eq_3_1 = theta(1) + (theta(2)+theta(3)*(Mag-7))*alog(R) + 1 theta(6)*rrup + f_mag + f_ztor + f_site + f_slab*evType 1 + f_basin write (21,1000) region, mag, evType, ZTOR, Rrup, 1 vs30, z25, rockPGA, specT, 1 mu, eq_3_1, f_mag, f_ztor, f_site, f_slab*evType, f_basin 1000 format (i5,f8.2,f8.0,4f8.2,2f9.4, 12f10.4) c ------------------------------------------------------ c Sigma models c set fixed coefficents d0 = 0.47 T1_phi2 = 0.03 T2_phi2 = 0.075 T3_phi2 = 0.20 T4_phi2 = 1.0 T1_phi3 = 0.03 T2_phi3 = 0.075 T3_phi3 = 0.10 T4_phi3 = 0.3 d3_phi2 = 0.109 d4_phi2 = 0.062 d5_phi2 = 0.470 d3_phi3 = 0.242 d4_phi3 = 0.0 d5_phi3 = 0.0 alpha_phi3 = 0.42 c tau model tau_LIN = d0 tau_LIN_PGA = d0 c phi1 model if ( rrup .lt. 150. ) then phi1_SQ = d1_T phi1_SQ_PGA = d1(1) elseif ( Rrup .lt. 450. ) then phi1_SQ = d1_T + d2_T * (Rrup-150)/300. phi1_SQ_PGA = d1(1) + d2(1) * (Rrup-150)/300. else phi1_SQ = d1_T + d2_T phi1_SQ_PGA = d1(1) + d2(1) endif c phi2 model C Set alpha for phi2 (eq 5.6) if (rrup .le. 250. ) then alpha_phi2 = 1. else if ( rrup .lt. 450. ) then alpha_phi2 = 1 - 0.0036*(Rrup-250.) else alpha_phi2 = 0.28 endif c f2 term for Phi2 (eq 5.3) if ( specT .lt. T1_phi2 ) then f2_phi2 = 1 - alpha_phi2 elseif ( specT .lt. T2_phi2 ) then f2_phi2 = 1 - alpha_phi2 *alog(specT/T2_phi2)/alog(T1_phi2/T2_phi2) elseif ( specT .lt. T3_phi2 ) then f2_phi2 = 1. elseif ( specT .lt. T4_phi2 ) then f2_phi2 = alog(specT/T4_phi2)/alog(T3_phi2/T4_phi2) else f2_phi2 = 0. endif f2_phi2_PGA = 1 - alpha_phi2 c A_phi term (eq 5.4) if ( Rrup .lt. 225. ) then A_phi2 = d3_phi2 elseif ( Rrup .lt. 450. ) then A_phi2 = d3_phi2 + d4_phi2 * (Rrup-225)/225. + d5_phi2 * ((Rrup-225)/225.)**2 else A_phi2 = d3_phi2 + d4_phi2 + d5_phi2 endif c compute phi2^2 for SpecT and for PGA phi2_SQ = A_phi2 * f2_phi2 phi2_SQ_PGA = A_phi2*f2_phi2_PGA c phi3 model c (same for as eq 5.3, but different coeff) alpha_Phi3 = 0.42 if ( specT .lt. T1_phi3 ) then f2_phi3 = 1 - alpha_phi3 elseif ( specT .lt. T2_phi3 ) then f2_phi3 = 1 - alpha_phi3 *alog(specT/T2_phi3)/alog(T1_phi3/T2_phi3) elseif ( specT .lt. T3_phi3 ) then f2_phi3 = 1. elseif ( specT .lt. T4_phi3 ) then f2_phi3 = alog(specT/T4_phi3)/alog(T3_phi3/T4_phi3) else f2_phi3 = 0. endif A_phi3 = d3_phi3 f2_phi3_PGA = 1 - alpha_phi3 c compute phi3^2 for SpecT and for PGA phi3_SQ = A_phi3 * f2_phi3 phi3_SQ_PGA = A_phi3*f2_phi3_PGA c compute region-specific phi_LIN and PhiSS_LIN if ( region .eq. 3) then phi_LIN = sqrt( phi1_SQ + phi2_SQ ) phi_LIN_PGA = sqrt( phi1_SQ_PGA + A_phi2*f2_phi2_PGA ) PhiSS_LIN = sqrt( phi_LIN**2 - phiS2S_G2T**2 ) PhiSS_LIN_PGA = sqrt( phi_LIN_PGA**2 - phiS2S_G2(1)**2 ) PhiS2S = phiS2S_G2T elseif ( region .eq. 4 .or. region .eq. 6 ) then phi_LIN = sqrt( phi1_SQ + phi2_SQ + phi3_SQ) phi_LIN_PGA = sqrt( phi1_SQ_PGA + phi2_SQ_PGA + phi3_SQ_PGA) PhiSS_LIN = sqrt( phi_LIN**2 - phiS2S_G3T**2 ) PhiSS_LIN_PGA = sqrt( phi_LIN_PGA**2 - phiS2S_G3(1)**2 ) PhiS2S = phiS2S_G3T else phi_LIN = sqrt(phi1_SQ) phi_LIN_PGA = sqrt(phi1_SQ_PGA) PhiSS_LIN = sqrt( phi_LIN**2 - phiS2S_G1T**2 ) PhiSS_LIN_PGA = sqrt( phi_LIN_PGA**2 - phiS2S_G1(1)**2 ) PhiS2S = phiS2S_G1T endif c NL site effects on phi and tau c eq 5.9 If ( vs .ge. VLIN_T ) then c No NL effects phi = phi_LIN phiSS = phiSS_LIN tau = tau_LIN sigma = sqrt( phi**2 + tau**2) sigmaSS = sqrt( phiSS**2 + tau**2) goto 200 else partial_f_PGA = b_soil_T * rockPGA * 1 ( -1./(rockPGA+c) + 1. / (rockPGA+c*(Vs/VLIN_T)**n) ) endif c eq 5.7 phi_AMP = 0.3 phi_B_PGA = sqrt( phi_LIN_PGA**2 - phi_AMP**2) phi_B = sqrt( phi_LIN**2 - phi_AMP**2) c eq 5.8 phiSQ_NL = phi_LIN**2 + partial_f_PGA**2 * phi_B_PGA**2 1 + 2. * partial_f_PGA * phi_B_PGA * phi_B * rho_W_T c eq 5.10 tauSQ_NL = tau_LIN**2 + partial_f_PGA**2 * tau_LIN_PGA**2 1 + 2. * partial_f_PGA * tau_LIN_PGA * tau_LIN * rho_B_T phi = sqrt(phiSQ_NL) tau = sqrt(tauSQ_NL) sigma = sqrt (phiSQ_NL + tauSQ_NL) c NL effects on phiSS and sigmaSS c eq 5.7 phi_AMP = 0.3 phiSS_B_PGA = sqrt( phiSS_LIN_PGA**2 - phi_AMP**2) phiSS_B = sqrt( phiSS_LIN**2 - phi_AMP**2) c eq 5.8 phiSS_SQ_NL = phiSS_LIN**2 + partial_f_PGA**2 * phiSS_B_PGA**2 1 + 2. * partial_f_PGA * phiSS_B_PGA * phiSS_B * rho_W_T phiSS = sqrt(phiSS_SQ_NL) sigmaSS = sqrt (phiSS_SQ_NL + tauSQ_NL) 200 continue write (22,1001) region, Rrup, vs30, rockPGA, specT, 1 tau_LIN, sqrt(phi1_SQ), sqrt(phi2_SQ), sqrt(phi3_SQ), phi_LIN, 2 phiS2S, phiSS_LIN, sqrt(phiSQ_NL), sqrt(tauSQ_NL), sigma, phiSS, sigmaSS, partial_f_PGA, 3 A_phi2, f2_phi2, A_phi3, f2_phi3 1001 format (i5,1x,f8.2,f8.0,2f8.2,20f10.4) return end c ------------------------------------------------------------------- subroutine z25_interp ( vs30, x1, x2, y1, y2, y) real vs30, x1, x2, y1, y2, y, x, x1Log, x2Log x = alog(vs30) x1Log = alog(x1) x2Log = alog(x2) if ( x .lt. x1Log ) then y = y1 elseif ( x .gt. x2Log ) then y = y2 else y = (x-x1Log)/(x2Log-x1Log) * (y2-y1) + y1 endif y = exp(y) return end