diff --git a/process/pfcoil.py b/process/pfcoil.py
index a31651891d..63451837cb 100644
--- a/process/pfcoil.py
+++ b/process/pfcoil.py
@@ -603,7 +603,7 @@ def pfcoil(self):
                 i = i + 1
 
         # Calculate peak field, allowable current density, resistive
-        # power losses and volumes and weights for each PF coil
+        # power losses and volumes and weights for each PF coil, index i
         i = 0
         it = 0
         pfv.powpfres = 0.0e0
@@ -620,8 +620,8 @@ def pfcoil(self):
                         i + 1, iii + 1, it
                     )  # returns bpf, bpf2
 
-                # Allowable current density (for superconducting coils)
-
+                # Issue 1871.  MDK
+                # Allowable current density (for superconducting coils) for each coil, index i
                 if pfv.ipfres == 0:
                     bmax = max(abs(pfv.bpf[i]), abs(pf.bpf2[i]))
                     pfv.rjpfalw[i], jstrand, jsc, tmarg = self.superconpf(
@@ -2905,11 +2905,11 @@ def deltaj_hijc_rebco(temperature):
             jcritstr = jcritsc * (1.0e0 - fcu)
 
             # The CS coil current at EOF
-            ioheof = bv.hmax * pfv.ohhghf * bv.ohcth * 2.0 * pfv.coheof
+            # ioheof = bv.hmax * pfv.ohhghf * bv.ohcth * 2.0 * pfv.coheof
             # The CS coil current/copper area calculation for quench protection
             # Copper area = (area of coil - area of steel)*(1- void fraction)*
             # (fraction of copper in strands)
-            rcv.copperaoh_m2 = ioheof / (pfv.awpoh * (1.0 - pfv.vfohc) * pfv.fcuohsu)
+            # rcv.copperaoh_m2 = ioheof / (pfv.awpoh * (1.0 - pfv.vfohc) * pfv.fcuohsu)
 
         elif isumat == 7:
             # Durham Ginzburg-Landau critical surface model for Nb-Ti
@@ -2919,7 +2919,7 @@ def deltaj_hijc_rebco(temperature):
             jcritstr = jcritsc * (1.0e0 - fcu)
 
             # The CS coil current at EOF
-            ioheof = bv.hmax * pfv.ohhghf * bv.ohcth * 2.0 * pfv.coheof
+            # ioheof = bv.hmax * pfv.ohhghf * bv.ohcth * 2.0 * pfv.coheof
 
         elif isumat == 8:
             # Durham Ginzburg-Landau critical surface model for REBCO
@@ -2930,9 +2930,9 @@ def deltaj_hijc_rebco(temperature):
             jcritstr = jcritsc * (1.0e0 - fcu)
 
             # The CS coil current at EOF
-            ioheof = bv.hmax * pfv.ohhghf * bv.ohcth * 2.0 * pfv.coheof
+            # ioheof = bv.hmax * pfv.ohhghf * bv.ohcth * 2.0 * pfv.coheof
             # The CS coil current/copper area calculation for quench protection
-            rcv.copperaoh_m2 = ioheof / (pfv.awpoh * (1.0 - pfv.vfohc) * pfv.fcuohsu)
+            # rcv.copperaoh_m2 = ioheof / (pfv.awpoh * (1.0 - pfv.vfohc) * pfv.fcuohsu)
 
         elif isumat == 9:
             # Hazelton experimental data + Zhai conceptual model for REBCO
@@ -2945,15 +2945,25 @@ def deltaj_hijc_rebco(temperature):
             jcritstr = jcritsc * (1.0e0 - fcu)
 
             # The CS coil current at EOF
-            ioheof = bv.hmax * pfv.ohhghf * bv.ohcth * 2.0 * pfv.coheof
+            # ioheof = bv.hmax * pfv.ohhghf * bv.ohcth * 2.0 * pfv.coheof
             # The CS coil current/copper area calculation for quench protection
-            rcv.copperaoh_m2 = ioheof / (pfv.awpoh * (1.0 - pfv.vfohc) * pfv.fcuohsu)
+            # rcv.copperaoh_m2 = ioheof / (pfv.awpoh * (1.0 - pfv.vfohc) * pfv.fcuohsu)
 
         else:
             # Error condition
             eh.idiag[0] = isumat
             eh.report_error(156)
 
+            # Issue 1871 MDK. The CS calculation has been removed from the isumat option list,
+            # and only calculated if the CS properties are needed.
+            if bv.iohcl != 0:
+                # CS coil current at EOF
+                ioheof = bv.hmax * pfv.ohhghf * bv.ohcth * 2.0 * pfv.coheof
+                # CS coil current/copper area calculation for quench protection
+                rcv.copperaoh_m2 = ioheof / (
+                    pfv.awpoh * (1.0 - pfv.vfohc) * pfv.fcuohsu
+                )
+
         # Critical current density in winding pack
         jcritwp = jcritstr * (1.0e0 - fhe)
         jstrand = jwp / (1.0e0 - fhe)