diff --git a/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java b/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java
new file mode 100644
index 0000000000000000000000000000000000000000..85de2d625d02ff1bd9edaec84125076b96675cc7
--- /dev/null
+++ b/src/main/java/gov/usgs/earthquake/nshmp/gmm/ParkerEtAl_2020.java
@@ -0,0 +1,274 @@
+package gov.usgs.earthquake.nshmp.gmm;
+
+import static gov.usgs.earthquake.nshmp.gmm.GmmInput.Field.MW;
+import static gov.usgs.earthquake.nshmp.gmm.GmmInput.Field.RRUP;
+import static gov.usgs.earthquake.nshmp.gmm.GmmInput.Field.VS30;
+import static gov.usgs.earthquake.nshmp.gmm.GmmInput.Field.ZHYP;
+import static java.lang.Math.exp;
+import static java.lang.Math.log;
+import static java.lang.Math.min;
+import static java.lang.Math.pow;
+
+import java.util.Map;
+
+import com.google.common.collect.Range;
+
+import gov.usgs.earthquake.nshmp.Maths;
+import gov.usgs.earthquake.nshmp.gmm.GmmInput.Constraints;
+
+/**
+ * Implementation of the PEER NGA-Subduction GMM, Parker et al., 2020 *
+ * <p><b>Reference:</b> Parker, G.A, Stewart, J.P., Boore, D.M., Atkinson, G.M.,
+ * Hassani, B. (2020) NGA-Subduction global ground-motion models with regional
+ * adjustment factors. Report no. 2020/03. Berkeley, CA: PEER, 131 p.
+ *
+ * <p><b>Component:</b> TODO
+ *
+ * @author U.S. Geological Survey
+ */
+
+public abstract class ParkerEtAl_2020 implements GroundMotionModel {
+
+ static final String NAME = "Parker et al. (2020) NGA-Subduction";
+
+ static final Constraints CONSTRAINTS_INTERFACE = Constraints.builder()
+ .set(MW, Range.closed(4.5, 9.5))
+ .set(RRUP, Range.closed(20.0, 1000.0))
+ .set(ZHYP, Range.closed(0.0, 40.0))
+ .set(VS30, Range.closed(150.0, 2000.0))
+ .build();
+ static final Constraints CONSTRAINTS_INTRASLAB = Constraints.builder()
+ .set(MW, Range.closed(4.5, 8.5))
+ .set(RRUP, Range.closed(35.0, 1000.0))
+ .set(ZHYP, Range.closed(0.0, 200.0))
+ .set(VS30, Range.closed(150.0, 2000.0))
+ .build();
+
+ static final CoefficientContainer COEFFS = new CoefficientContainer("TBD.csv");
+
+ private static final double V1 = 270.0; // m/s
+ private static final double VREF = 760.0; // m/s
+ private static final double DB1 = 20.0; // km
+ private static final double DB2 = 67.0; // km
+
+ private static final class Coefficients {
+ final Imt imt;
+ final double c0, c1, b4, a0, c4, c5, c6, d, m, db, v2, s1, s2, f4, f5, e1, e2, e3;
+
+ Coefficients(Imt imt, CoefficientContainer cc) {
+ this.imt = imt;
+ Map<String, Double> coeffs = cc.get(imt);
+ c0 = 0.0;
+ c1 = 0.0;
+ b4 = 0.0;
+ a0 = 0.0;
+ c4 = 0.0;
+ c5 = 0.0;
+ c6 = 0.0;
+ d = 0.0;
+ m = 0.0;
+ db = 0.0;
+ v2 = 0.0;
+ s1 = 0.0; // s1 = s2 except for Japan and Taiwan
+ s2 = 0.0;
+ f4 = 0.0;
+ f5 = 0.0;
+ e1 = 0.0;
+ e2 = 0.0;
+ e3 = 0.0;
+ }
+ }
+
+ private final Coefficients coeffs;
+ private final Coefficients coeffsPGA;
+
+ ParkerEtAl_2020(final Imt imt) {
+ coeffs = new Coefficients(imt, COEFFS);
+ coeffsPGA = new Coefficients(imt, COEFFS);
+ }
+
+ @Override
+ public final ScalarGroundMotion calc(final GmmInput in) {
+ /*
+ * Equation 4.1 μ_lny = c0 + Fp + Fm + Fd + Fs
+ */
+
+ double pgaRef =
+ exp(calcMean(coeffsPGA, isSlab(), false, getMc(), 0.0, in.Mw, in.rRup, in.zHyp, in.vs30,
+ in.z2p5));
+
+ double μ = calcMean(coeffs, isSlab(), useBasinTerm(), getMc(), pgaRef, in.Mw, in.rRup, in.zHyp,
+ in.vs30, in.z2p5);
+
+ double σ = 0.0;
+ return DefaultScalarGroundMotion.create(μ, σ);
+ }
+
+ abstract boolean isSlab();
+
+ // Note: mc, corner magnitude, is a regional saturation magnitude
+ abstract double getMc();
+
+ abstract boolean useBasinTerm();
+
+ /* Equation 4.1 */
+ private static final double calcMean(
+ final Coefficients c,
+ final boolean slab,
+ final boolean useBasinTerm,
+ final double mc,
+ final double pgaRef,
+ final double Mw,
+ final double rRup,
+ final double zHyp,
+ final double vs30,
+ final double z2p5) {
+
+ /* Equation 4.7: Fs = Flin + Fnl + Fb */
+ double siteTerm = getSiteTermLinear(c, vs30) + getSiteTermNonLinear(c, pgaRef, vs30) +
+ getSiteBasinTerm(c, useBasinTerm, z2p5);
+
+ double h = nearSourceSaturationTerm(Mw, mc, slab);
+
+ double μ = c.c0 +
+ getPathTerm(c, rRup, Mw, h) +
+ getMagnitudeTerm(c, mc, Mw) +
+ getDepthTerm(c, slab, zHyp) +
+ siteTerm;
+
+ return 0.0;
+ }
+
+ /* Equations 4.2 and 4.3 */
+ private static final double getPathTerm(final Coefficients c, double rRup, double Mw, double h) {
+ double R = Maths.hypot(rRup, h);
+ return c.c1 * log(R) + c.b4 * Mw * log(R) + c.a0 * R;
+ }
+
+ private static final double nearSourceSaturationTerm(double Mw, double mc, boolean slab) {
+ if (slab) {
+ /* Equation 4.4b */
+ return (Mw <= mc)
+ ? pow(10, (1.050 / (mc - 4.0)) * (Mw - mc) + 1.544)
+ : 35;
+ } else {
+ /* Equation 4.4a */
+ return pow(10.0, -0.82 + 0.252 * Mw);
+ }
+ }
+
+ /* Equation 4.5 */
+ private static final double getMagnitudeTerm(final Coefficients c, double mc, double Mw) {
+ double dM = Mw - mc;
+ return (Mw <= mc)
+ ? c.c4 * dM + c.c5 * dM * dM
+ : c.c6 * dM;
+ }
+
+ /* Equation 4.6 */
+ private static final double getDepthTerm(final Coefficients c, boolean slab, double zHyp) {
+ if (!slab) {
+ return 0.0;
+ }
+ // Zhyp (see section 4.3.3) can be replaced with mean depth (Zhyp_bar)
+ // that depends on Ztor, fault width, and dip angle
+ double dTerm;
+ if (zHyp < DB1) { // "less than"
+ dTerm = c.m * (DB1 - DB2) + c.d;
+ } else if (zHyp > DB2) {
+ dTerm = c.d;
+ } else { // "db1 <= zHyp" possible typo? (written as "db1 < xHyp")
+ dTerm = c.m * (zHyp - DB2) + c.d;
+ }
+ return dTerm;
+ }
+
+ /* Equation 4.8 */
+ private static final double getSiteTermLinear(final Coefficients c, double vs30) {
+ double Flin;
+ if (vs30 <= V1) {
+ Flin = c.s1 * log(vs30 / V1) + c.s2 * log(V1 / VREF);
+ } else if (vs30 > c.v2) {
+ Flin = c.s2 * log(c.v2 / VREF);
+ } else {
+ Flin = c.s2 * log(vs30 / VREF);
+ }
+ return Flin;
+ }
+
+ /* Equation 4.9 */
+ private static final double getSiteTermNonLinear(final Coefficients c, double pgaRef,
+ double vs30) {
+ double f2 = c.f4 * (exp(c.f5 * (min(vs30, 760) - 200)) - exp(c.f5 * (760 - 200)));
+ return 0.0 + f2 * log((pgaRef + 0.05) / 0.05);
+ }
+
+ /* Equations 4.11 - 4.13 */
+ private static final double getSiteBasinTerm(final Coefficients c, boolean basin, double z2p5) {
+ if (!basin) {
+ return 0.0;
+ } else {
+ // TODO: implement for Cascadia - this term is 0 except in Japan and
+ // Cascadia
+ return 0.0;
+ }
+ }
+
+ /*
+ * Global interface model
+ */
+ static class Interface extends ParkerEtAl_2020 {
+ static final String NAME = ParkerEtAl_2020.NAME + " : Interface";
+
+ Interface(Imt imt) {
+ super(imt);
+ }
+
+ @Override
+ final boolean isSlab() {
+ return false;
+ }
+
+ @Override
+ final boolean useBasinTerm() {
+ return false;
+ }
+
+ @Override
+ final double getMc() {
+ return 7.9;
+ }
+
+ }
+
+ /*
+ * Global intraslab model
+ */
+ static class Slab extends ParkerEtAl_2020 {
+ static final String NAME = ParkerEtAl_2020.NAME + " : Slab";
+
+ Slab(Imt imt) {
+ super(imt);
+ }
+
+ @Override
+ final boolean isSlab() {
+ return true;
+ }
+
+ @Override
+ final boolean useBasinTerm() {
+ return false;
+ }
+
+ @Override
+ final double getMc() {
+ return 7.6;
+ }
+ }
+
+ public static void main(String[] args) {
+
+ }
+
+}