Repository: oscarbranson/cbsyst

Issue:
Implement different pH scales
- [x] Total
- [x] SWS
- [x] Free
- [x] NBS

Discussion:

Patch:
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 5524ea6..f26c8bb 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -45,16 +45,4 @@ Currently using the unittest module, and testing via ``setup.py test``.
 
 Current unittests check internal consistency of functions against stable-state reference values, and compare the output of Csys against reference carbon speciation data. 
 
-*Do not change existing tests without good reason*.
-
-## Project Structure
-
-```
-cbsyst/
-  |--- boron_fns.py  : Functions for calculating relating to B speciation / isotopes.
-  |--- carbon_fns.py  : Functions for calculating C speciation.
-  |--- cbsyst.py  : User-facing functions used for calculating seawater chemistry.
-  |--- helpers.py  : General functions that are used elsewhere.
-  |--- MyAMI_V2.py  : Functions for K0, K1, K2, KB, KW, KS, KspA, KspC and [Mg] and [Ca] corrections.
-  |--- non_MyAMI_constants.py  : Functions to calculate any constants that are not handled by MyAMI.
-```
\ No newline at end of file
+*Do not change existing tests without good reason*.
\ No newline at end of file
diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100644
index 0000000..f076406
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,2 @@
+graft cbsyst/test_data/Lueker2000
+graft cbsyst/test_data/GLODAP_data
\ No newline at end of file
diff --git a/README.md b/README.md
index dc42ad5..930efa9 100644
--- a/README.md
+++ b/README.md
@@ -1,11 +1,10 @@
+# cbsyst
+
 <div align="right">
   <a href="https://travis-ci.org/oscarbranson/cbsyst"><img src="https://travis-ci.org/oscarbranson/cbsyst.svg?branch=master" alt="TravisCI build" height="18"></a>
   <a href="https://badge.fury.io/py/cbsyst"><img src="https://badge.fury.io/py/cbsyst.svg" alt="PyPI version" height="18"></a>
-  <a href="https://anaconda.org/conda-forge/cbsyst"> <img src="https://anaconda.org/conda-forge/cbsyst/badges/version.svg" ald="conda-forge version" height="18"/></a>
 </div>
 
-# cbsyst
-
 **A Python module for calculating seawater carbon and boron chemistry.** 
 
 This will be particularly useful for anyone thinking about oceans in the distant past, when Mg and Ca concentrations were different. I use [Mathis Hain's MyAMI model](http://www.mathis-hain.net/resources/Hain_et_al_2015_GBC.pdf) to adjust speciation constants for Mg and Ca concentration.
@@ -20,7 +19,7 @@ This will be particularly useful for anyone thinking about oceans in the distant
 - [ ] [Compare to CO2SYS](https://github.com/oscarbranson/cbsyst/issues/6), a la [Orr et al (2015)](http://www.biogeosciences.net/12/1483/2015/bg-12-1483-2015.pdf)?
 
 If anyone wants to help with any of this, please do contribute!
-A full list of bite-sized tasks that need doing is available on the [Issues](https://github.com/oscarbranson/cbsyst/issues) page.
+A full list of bite-sized tasks that need doing is available in the [TODO list]().
 
 ## Acknowledgement
 The development of `cbsyst` has been greatly aided by [CO2SYS](http://cdiac.ornl.gov/oceans/co2rprt.html), and the [Matlab conversion of CO2SYS](http://cdiac.ornl.gov/ftp/oceans/co2sys/).
@@ -81,15 +80,12 @@ Delta values can be provided as an input, and are given as an output.
 **Requires Python 3.5+**. 
 Does *not* work in 2.7. Sorry.
 
-### PyPi
 ```bash
 pip install cbsyst
 ```
 
-### Conda-Forge
-```bash
-conda install cbsyst -c conda-forge
-```
+## Build Status
+
 
 ## Example Usage
 
diff --git a/cbsyst/MyAMI_V2.py b/cbsyst/MyAMI_V2.py
index 6b5286c..d887f87 100644
--- a/cbsyst/MyAMI_V2.py
+++ b/cbsyst/MyAMI_V2.py
@@ -1266,7 +1266,7 @@ start_params = {'K0': np.array([-60.2409, 93.4517, 23.3585, 0.023517,
                 'KB': np.array([148.0248, 137.1942, 1.62142, -8966.90,
                                 -2890.53, -77.942, 1.728, -0.0996,
                                 -24.4344, -25.085, -0.2474, 0.053105]),
-                'KW': np.array([148.9652, -13847.26, -23.6521, 118.67,
+                'KW': np.array([148.9802, -13847.26, -23.6521, 118.67,
                                 -5.977, 1.0495, -0.01615]),
                 'KspC': np.array([-171.9065, -0.077993, 2839.319, 71.595, -0.77712,
                                   0.0028426, 178.34, -0.07711, 0.0041249]),
@@ -1591,7 +1591,7 @@ def MyAMI_K_calc(TempC=25., Sal=35., Ca=0.0102821, Mg=0.0528171, P=None, param_d
     return Ks
 
 
-def MyAMI_K_calc_multi(T=25., S=35., Ca=0.0102821, Mg=0.0528171, P=None):
+def MyAMI_K_calc_multi(TempC=25., Sal=35., Ca=0.0102821, Mg=0.0528171, P=None):
     """
     Calculate MyAMI equilibrium constants for multiple T, S and Mg and Ca conditions.
 
@@ -1617,8 +1617,8 @@ def MyAMI_K_calc_multi(T=25., S=35., Ca=0.0102821, Mg=0.0528171, P=None):
     """
     # package data in a bunch of 1D arrays.
     d = Bunch()
-    d.T = np.array(T, ndmin=1)
-    d.S = np.array(S, ndmin=1)
+    d.T = np.array(TempC, ndmin=1)
+    d.S = np.array(Sal, ndmin=1)
     d.Ca = np.array(Ca, ndmin=1)
     d.Mg = np.array(Mg, ndmin=1)
     d.P = np.array(P, ndmin=1)
diff --git a/cbsyst/boron_fns.py b/cbsyst/boron_fns.py
index 61939b8..cd8790e 100644
--- a/cbsyst/boron_fns.py
+++ b/cbsyst/boron_fns.py
@@ -1,5 +1,5 @@
 import numpy as np
-from cbsyst.helpers import ch
+from cbsyst.helpers import ch, cp, Bunch
 
 
 # B conc fns
@@ -160,3 +160,39 @@ def R11_2_A11(R11):
     Convert Ratio to Abundance notation.
     """
     return R11 / (1 + R11)
+
+
+def calc_B_species(pHtot=None, BT=None, BO3=None, BO4=None, Ks=None):
+    # B system calculations
+    if pHtot is not None and BT is not None:
+        H = ch(pHtot)
+    elif BT is not None and BO3 is not None:
+        H = BT_BO3(BT, BO3, Ks)
+    elif BT is not None and BO4 is not None:
+        H = BT_BO4(BT, BO4, Ks)
+    elif BO3 is not None and BO4 is not None:
+        BT = BO3 + BO3
+        H = BT_BO3(BT, BO3, Ks)
+    elif pHtot is not None and BO3 is not None:
+        H = ch(pHtot)
+        BT = pH_BO3(pHtot, BO3, Ks)
+    elif pHtot is not None and BO4 is not None:
+        H = ch(pHtot)
+        BT = pH_BO4(pHtot, BO4, Ks)
+
+    # The above makes sure that BT and H are known,
+    # this next bit calculates all the missing species
+    # from BT and H.
+
+    if BO3 is None:
+        BO3 = cBO3(BT, H, Ks)
+    if BO4 is None:
+        BO4 = cBO4(BT, H, Ks)
+    if pHtot is None:
+        pHtot = np.array(cp(H), ndmin=1)
+
+    return Bunch({'pHtot': pHtot,
+                  'H': H,
+                  'BT': BT,
+                  'BO3': BO3,
+                  'BO4': BO4})
diff --git a/cbsyst/carbon_fns.py b/cbsyst/carbon_fns.py
index 2fda38c..13093b2 100644
--- a/cbsyst/carbon_fns.py
+++ b/cbsyst/carbon_fns.py
@@ -1,6 +1,6 @@
 import scipy.optimize as opt
 import numpy as np
-from cbsyst.helpers import ch, noms, cast_array, maxL
+from cbsyst.helpers import ch, noms, cast_array, maxL, calc_pH_scales, Bunch, cp
 # from cbsyst.boron_fns import cBO4
 
 
@@ -353,9 +353,9 @@ def TA_DIC(TA, DIC, BT, TP, TSi, TS, TF, Ks):
 #     Returns H
 #     """
 #     TA, DIC, BT = noms(TA, DIC, BT)  # get nominal values of inputs
-#     par = cast_array(TA, DIC, BT, Ks.K1, Ks.K2, Ks.KB, Ks.KW)  # cast parameters into array
+#     = cast_array(TA, DIC, BT, Ks.K1, Ks.K2, Ks.KB, Ks.KW)  # cast meters into array
 
-#     return np.apply_along_axis(_zero_wrapper, 0, par, fn=zero_TA_DIC)
+#     return np.apply_along_axis(_zero_wrapper, 0,  fn=zero_TA_DIC)
 
 
 def zero_TA_DIC(h, TA, DIC, BT, K1, K2, KB, KW):
@@ -503,3 +503,156 @@ def fCO2_to_pCO2(fCO2, Tc):
 
     return fCO2 / np.exp(P * (B + 2 * delta) / RT)
 
+
+def calc_C_species(pHtot=None, DIC=None, CO2=None,
+                   HCO3=None, CO3=None, TA=None,
+                   fCO2=None, pCO2=None,
+                   T_in=None, BT=None, TP=0, TSi=0,
+                   TS=0, TF=0, Ks=None):
+    """
+    Calculate all carbon species from minimal input.
+    """
+
+    # if fCO2 is given but CO2 is not, calculate CO2
+    if CO2 is None:
+        if fCO2 is not None:
+            CO2 = fCO2_to_CO2(fCO2, Ks)
+        elif pCO2 is not None:
+            CO2 = fCO2_to_CO2(pCO2_to_fCO2(pCO2, T_in), Ks)
+
+    # Carbon System Calculations (from Zeebe & Wolf-Gladrow, Appendix B)
+    # 1. CO2 and pH
+    if CO2 is not None and pHtot is not None:
+        H = ch(pHtot)
+        DIC = CO2_pH(CO2, pHtot, Ks)
+    # 2. CO2 and HCO3
+    elif CO2 is not None and HCO3 is not None:
+        H = CO2_HCO3(CO2, HCO3, Ks)
+        DIC = CO2_pH(CO2, cp(H), Ks)
+    # 3. CO2 and CO3
+    elif CO2 is not None and CO3 is not None:
+        H = CO2_CO3(CO2, CO3, Ks)
+        DIC = CO2_pH(CO2, cp(H), Ks)
+    # 4. CO2 and TA
+    elif CO2 is not None and TA is not None:
+        # unit conversion because OH and H wrapped
+        # up in TA fns - all need to be in same units.
+        pHtot = CO2_TA(CO2=CO2,
+                       TA=TA,
+                       BT=BT,
+                       TP=TP,
+                       TSi=TSi,
+                       TS=TS,
+                       TF=TF,
+                       Ks=Ks)
+        H = ch(pHtot)
+        DIC = CO2_pH(CO2, pHtot, Ks)
+    # 5. CO2 and DIC
+    elif CO2 is not None and DIC is not None:
+        H = CO2_DIC(CO2, DIC, Ks)
+    # 6. pHtot and HCO3
+    elif pHtot is not None and HCO3 is not None:
+        H = ch(pHtot)
+        DIC = pH_HCO3(pHtot, HCO3, Ks)
+    # 7. pHtot and CO3
+    elif pHtot is not None and CO3 is not None:
+        H = ch(pHtot)
+        DIC = pH_CO3(pHtot, CO3, Ks)
+    # 8. pHtot and TA
+    elif pHtot is not None and TA is not None:
+        H = ch(pHtot)
+        DIC = pH_TA(pH=pHtot,
+                    TA=TA,
+                    BT=BT,
+                    TP=TP,
+                    TSi=TSi,
+                    TS=TS,
+                    TF=TF,
+                    Ks=Ks)
+    # 9. pHtot and DIC
+    elif pHtot is not None and DIC is not None:
+        H = ch(pHtot)
+    # 10. HCO3 and CO3
+    elif HCO3 is not None and CO3 is not None:
+        H = HCO3_CO3(HCO3, CO3, Ks)
+        DIC = pH_CO3(cp(H), CO3, Ks)
+    # 11. HCO3 and TA
+    elif HCO3 is not None and TA is not None:
+        Warning('Nutrient alkalinity not implemented for this input combination.\nCalculations use only C and B alkalinity.')
+        H = HCO3_TA(HCO3,
+                    TA,
+                    BT,
+                    Ks)
+        DIC = pH_HCO3(cp(H), HCO3, Ks)
+    # 12. HCO3 amd DIC
+    elif HCO3 is not None and DIC is not None:
+        H = HCO3_DIC(HCO3, DIC, Ks)
+    # 13. CO3 and TA
+    elif CO3 is not None and TA is not None:
+        Warning('Nutrient alkalinity not implemented for this input combination.\nCalculations use only C and B alkalinity.')
+        H = CO3_TA(CO3,
+                   TA,
+                   BT,
+                   Ks)
+        DIC = pH_CO3(cp(H), CO3, Ks)
+    # 14. CO3 and DIC
+    elif CO3 is not None and DIC is not None:
+        H = CO3_DIC(CO3, DIC, Ks)
+    # 15. TA and DIC
+    elif TA is not None and DIC is not None:
+        pHtot = TA_DIC(TA=TA,
+                       DIC=DIC,
+                       BT=BT,
+                       TP=TP,
+                       TSi=TSi,
+                       TS=TS,
+                       TF=TF,
+                       Ks=Ks)
+        H = ch(pHtot)
+
+    # The above makes sure that DIC and H are known,
+    # this next bit calculates all the missing species
+    # from DIC and H.
+    if CO2 is None:
+        CO2 = cCO2(H, DIC, Ks)
+    if fCO2 is None:
+        fCO2 = CO2_to_fCO2(CO2, Ks)
+    if pCO2 is None:
+        pCO2 = fCO2_to_pCO2(fCO2, T_in)
+    if HCO3 is None:
+        HCO3 = cHCO3(H, DIC, Ks)
+    if CO3 is None:
+        CO3 = cCO3(H, DIC, Ks)
+    # Calculate all elements of Alkalinity
+    (TA, CAlk, BAlk,
+     PAlk, SiAlk, OH,
+     Hfree, HSO4, HF) = cTA(H=H,
+                            DIC=DIC,
+                            BT=BT,
+                            TP=TP,
+                            TSi=TSi,
+                            TS=TS,
+                            TF=TF,
+                            Ks=Ks, mode='multi')
+
+    # if pH not calced yet, calculate on all scales.
+    if pHtot is None:
+        pHtot = np.array(cp(H), ndmin=1)
+
+    return Bunch({'pHtot': pHtot,
+                  'TA': TA,
+                  'DIC': DIC,
+                  'CO2': CO2,
+                  'H': H,
+                  'HCO3': HCO3,
+                  'fCO2': fCO2,
+                  'pCO2': pCO2,
+                  'CO3': CO3,
+                  'CAlk': CAlk,
+                  'BAlk': BAlk,
+                  'PAlk': PAlk,
+                  'SiAlk': SiAlk,
+                  'OH': OH,
+                  'Hfree': Hfree,
+                  'HSO4': HSO4,
+                  'HF': HF})
\ No newline at end of file
diff --git a/cbsyst/cbsyst.py b/cbsyst/cbsyst.py
index 31c7cc6..f7a06b1 100644
--- a/cbsyst/cbsyst.py
+++ b/cbsyst/cbsyst.py
@@ -9,48 +9,48 @@ from cbsyst.non_MyAMI_constants import *
 
 # Helper functions
 # ----------------
-def get_Ks(ps):
+def calc_Ks(T, S, P, Mg, Ca, TS, TF, Ks=None):
     """
     Helper function to calculate Ks.
 
-    If ps.Ks is a dict, those Ks are used
-    transparrently, with no pressure modification.
+    If Ks is a dict, those Ks are used
+    transparrently (i.e. no pressure modification).
     """
-    if isinstance(ps.Ks, dict):
-        Ks = Bunch(ps.Ks)
+    if isinstance(Ks, dict):
+        Ks = Bunch(Ks)
     else:
-        if maxL(ps.Mg, ps.Ca) == 1:
-            if ps.Mg is None:
-                ps.Mg = 0.0528171
-            if ps.Ca is None:
-                ps.Ca = 0.0102821
-            Ks = MyAMI_K_calc(TempC=ps.T, Sal=ps.S,
-                              Mg=ps.Mg, Ca=ps.Ca, P=ps.P)
+        if maxL(Mg, Ca) == 1:
+            if Mg is None:
+                Mg = 0.0528171
+            if Ca is None:
+                Ca = 0.0102821
+            Ks = MyAMI_K_calc(TempC=T, Sal=S, P=P,
+                              Mg=Mg, Ca=Ca)
         else:
             # if only Ca or Mg provided, fill in other with modern
-            if ps.Mg is None:
-                ps.Mg = 0.0528171
-            if ps.Ca is None:
-                ps.Ca = 0.0102821
+            if Mg is None:
+                Mg = 0.0528171
+            if Ca is None:
+                Ca = 0.0102821
             # calculate Ca and Mg specific Ks
-            Ks = MyAMI_K_calc_multi(TempC=ps.T, Sal=ps.S,
-                                    Ca=ps.Ca, Mg=ps.Mg, P=ps.P)
+            Ks = MyAMI_K_calc_multi(TempC=T, Sal=S, P=P,
+                                    Ca=Ca, Mg=Mg)
 
-    # non-MyAMI Constants
-    Ks.update(calc_KPs(ps.T, ps.S, ps.P))
-    Ks.update(calc_KF(ps.T, ps.S, ps.P))
-    Ks.update(calc_KSi(ps.T, ps.S, ps.P))
+        # non-MyAMI Constants
+        Ks.update(calc_KPs(T, S, P))
+        Ks.update(calc_KF(T, S, P))
+        Ks.update(calc_KSi(T, S, P))
 
-    # pH conversions to total scale.
-    #   - KPn are all on SWS
-    #   - KSi is on SWS
-    #   - MyAMI KW is on SWS... DOES THIS MATTER?
+        # pH conversions to total scale.
+        #   - KP1, KP2, KP3 are all on SWS
+        #   - KSi is on SWS
+        #   - MyAMI KW is on SWS... DOES THIS MATTER?
 
-    SWStoTOT = (1 + ps.TS / Ks.KSO4) / (1 + ps.TS / Ks.KSO4 + ps.TF / Ks.KF)
-    # FREEtoTOT = 1 + ps.TS / Ks.KSO4
-    conv = ['KP1', 'KP2', 'KP3', 'KSi', 'KW']
-    for c in conv:
-        Ks[c] *= SWStoTOT
+        SWStoTOT = (1 + TS / Ks.KSO4) / (1 + TS / Ks.KSO4 + TF / Ks.KF)
+        # FREEtoTOT = 1 + 'T_' + mode]S / Ks.KSO4
+        conv = ['KP1', 'KP2', 'KP3', 'KSi', 'KW']
+        for c in conv:
+            Ks[c] *= SWStoTOT
 
     return Ks
 
@@ -61,9 +61,10 @@ def Csys(pHtot=None, DIC=None, CO2=None,
          HCO3=None, CO3=None, TA=None,
          fCO2=None, pCO2=None,
          BT=None, Ca=None, Mg=None,
-         T=25., S=35., P=None,
+         T_in=25., S_in=35., P_in=None,
+         T_out=None, S_out=None, P_out=None,
          TP=0., TSi=0.,
-         pHsws=None, pHfree=None,
+         pHsws=None, pHfree=None, pHNBS=None,
          Ks=None, pdict=None, unit='umol'):
     """
     Calculate the carbon chemistry of seawater from a minimal parameter set.
@@ -97,12 +98,20 @@ def Csys(pHtot=None, DIC=None, CO2=None,
     Ca, Mg : arra-like
         The [Ca] and [Mg] of the seawater, in mol / kg.
         Used in calculating MyAMI constants.
-    T, S : array-like
-        Temperature in Celcius and Salinity in PSU.
-        Used in calculating MyAMI constants.
-    P : array-like
-        Pressure in Bar.
-        Used in calculating MyAMI constants.
+    T_in, S_in : array-like
+        Temperature in Celcius and Salinity in PSU that the
+        measurements were conducted under.
+        Used in calculating constants.
+    P_in : array-like
+        Pressure in Bar that the measurements were conducted under.
+        Used in pressure-correcting constants.
+    T_out, S_out : array-like
+        Temperature in Celcius and Salinity in PSU of the desired
+        output conditions.
+        Used in calculating constants.
+    P_in : array-like
+        Pressure in Bar of the desired output conditions.
+        Used in pressure-correcting constants.
     unit : str
         Concentration units of C and B parameters (all must be in
         the same units).
@@ -150,164 +159,70 @@ def Csys(pHtot=None, DIC=None, CO2=None,
 
     # Conserved seawater chemistry
     if 'TS' not in ps:
-        ps.TS = calc_TS(ps.S)
+        ps.TS = calc_TS(ps.S_in)
     if 'TF' not in ps:
-        ps.TF = calc_TF(ps.S)
+        ps.TF = calc_TF(ps.S_in)
     if ps.BT is None:
-        ps.BT = calc_TB(ps.S)
-    elif isinstance(BT, (int, float)):
-        ps.BT = ps.BT * ps.S / 35.
-
-    # Calculate Ks
-    ps.Ks = get_Ks(ps)
-
-    # Calculate pH scales (does nothing if no pH given)
-    ps.update(calc_pH_scales(ps.pHtot, ps.pHfree, ps.pHsws,
-                             ps.TS, ps.TF, ps.Ks))
-
-    # if fCO2 is given but CO2 is not, calculate CO2
-    if ps.CO2 is None:
-        if ps.fCO2 is not None:
-            ps.CO2 = fCO2_to_CO2(ps.fCO2, ps.Ks)
-        elif ps.pCO2 is not None:
-            ps.CO2 = fCO2_to_CO2(pCO2_to_fCO2(ps.pCO2, ps.T), ps.Ks)
-
-    # Carbon System Calculations (from Zeebe & Wolf-Gladrow, Appendix B)
-    # 1. CO2 and pH
-    if ps.CO2 is not None and ps.pHtot is not None:
-        ps.H = ch(ps.pHtot)
-        ps.DIC = CO2_pH(ps.CO2, ps.pHtot, ps.Ks)
-    # 2. ps.CO2 and ps.HCO3
-    elif ps.CO2 is not None and ps.HCO3 is not None:
-        ps.H = CO2_HCO3(ps.CO2, ps.HCO3, ps.Ks)
-        ps.DIC = CO2_pH(ps.CO2, cp(ps.H), ps.Ks)
-    # 3. ps.CO2 and ps.CO3
-    elif ps.CO2 is not None and ps.CO3 is not None:
-        ps.H = CO2_CO3(ps.CO2, ps.CO3, ps.Ks)
-        ps.DIC = CO2_pH(ps.CO2, cp(ps.H), ps.Ks)
-    # 4. ps.CO2 and ps.TA
-    elif ps.CO2 is not None and ps.TA is not None:
-        # unit conversion because OH and H wrapped
-        # up in TA fns - all need to be in same units.
-        ps.pHtot = CO2_TA(CO2=ps.CO2,
-                          TA=ps.TA,
-                          BT=ps.BT,
-                          TP=ps.TP,
-                          TSi=ps.TSi,
-                          TS=ps.TS,
-                          TF=ps.TF,
-                          Ks=ps.Ks)
-        ps.H = ch(ps.pHtot)
-        ps.DIC = CO2_pH(ps.CO2, ps.pHtot, ps.Ks)
-    # 5. ps.CO2 and ps.DIC
-    elif ps.CO2 is not None and ps.DIC is not None:
-        ps.H = CO2_DIC(ps.CO2, ps.DIC, ps.Ks)
-    # 6. ps.pHtot and ps.HCO3
-    elif ps.pHtot is not None and ps.HCO3 is not None:
-        ps.H = ch(ps.pHtot)
-        ps.DIC = pH_HCO3(ps.pHtot, ps.HCO3, ps.Ks)
-    # 7. ps.pHtot and ps.CO3
-    elif ps.pHtot is not None and ps.CO3 is not None:
-        ps.H = ch(ps.pHtot)
-        ps.DIC = pH_CO3(ps.pHtot, ps.CO3, ps.Ks)
-    # 8. ps.pHtot and ps.TA
-    elif ps.pHtot is not None and ps.TA is not None:
-        ps.H = ch(ps.pHtot)
-        ps.DIC = pH_TA(pH=ps.pHtot,
-                       TA=ps.TA,
-                       BT=ps.BT,
-                       TP=ps.TP,
-                       TSi=ps.TSi,
-                       TS=ps.TS,
-                       TF=ps.TF,
-                       Ks=ps.Ks)
-    # 9. ps.pHtot and ps.DIC
-    elif ps.pHtot is not None and ps.DIC is not None:
-        ps.H = ch(ps.pHtot)
-    # 10. ps.HCO3 and ps.CO3
-    elif ps.HCO3 is not None and ps.CO3 is not None:
-        ps.H = HCO3_CO3(ps.HCO3, ps.CO3, ps.Ks)
-        ps.DIC = pH_CO3(cp(ps.H), ps.CO3, ps.Ks)
-    # 11. ps.HCO3 and ps.TA
-    elif ps.HCO3 is not None and ps.TA is not None:
-        Warning('Nutrient alkalinity not implemented for this input combination.\nCalculations use only C and B alkalinity.')
-        ps.H = HCO3_TA(ps.HCO3,
-                       ps.TA,
-                       ps.BT,
-                       ps.Ks)
-        ps.DIC = pH_HCO3(cp(ps.H), ps.HCO3, ps.Ks)
-    # 12. ps.HCO3 amd ps.DIC
-    elif ps.HCO3 is not None and ps.DIC is not None:
-        ps.H = HCO3_DIC(ps.HCO3, ps.DIC, ps.Ks)
-    # 13. ps.CO3 and ps.TA
-    elif ps.CO3 is not None and ps.TA is not None:
-        Warning('Nutrient alkalinity not implemented for this input combination.\nCalculations use only C and B alkalinity.')
-        ps.H = CO3_TA(ps.CO3,
-                      ps.TA,
-                      ps.BT,
-                      ps.Ks)
-        ps.DIC = pH_CO3(cp(ps.H), ps.CO3, ps.Ks)
-    # 14. ps.CO3 and ps.DIC
-    elif ps.CO3 is not None and ps.DIC is not None:
-        ps.H = CO3_DIC(ps.CO3, ps.DIC, ps.Ks)
-    # 15. ps.TA and ps.DIC
-    elif ps.TA is not None and ps.DIC is not None:
-        ps.pHtot = TA_DIC(TA=ps.TA,
-                          DIC=ps.DIC,
-                          BT=ps.BT,
-                          TP=ps.TP,
-                          TSi=ps.TSi,
-                          TS=ps.TS,
-                          TF=ps.TF,
-                          Ks=ps.Ks)
-        ps.H = ch(ps.pHtot)
-
-    # The above makes sure that DIC and H are known,
-    # this next bit calculates all the missing species
-    # from DIC and H.
-    if ps.CO2 is None:
-        ps.CO2 = cCO2(ps.H, ps.DIC, ps.Ks)
-    if ps.fCO2 is None:
-        ps.fCO2 = CO2_to_fCO2(ps.CO2, ps.Ks)
-    if ps.pCO2 is None:
-        ps.pCO2 = fCO2_to_pCO2(ps.fCO2, ps.T)
-    if ps.HCO3 is None:
-        ps.HCO3 = cHCO3(ps.H, ps.DIC, ps.Ks)
-    if ps.CO3 is None:
-        ps.CO3 = cCO3(ps.H, ps.DIC, ps.Ks)
-    # Always calculate elements of alkalinity
-    try:
-        # necessary for use with CBsyst in special cases
-        # where BT is not known before Csys is run.
-        (ps.TA, ps.CAlk, ps.BAlk,
-         ps.PAlk, ps.SiAlk, ps.OH,
-         ps.Hfree, ps.HSO4, ps.HF) = cTA(H=ps.H,
-                                         DIC=ps.DIC,
-                                         BT=ps.BT,
-                                         TP=ps.TP,
-                                         TSi=ps.TSi,
-                                         TS=ps.TS,
-                                         TF=ps.TF,
-                                         Ks=ps.Ks, mode='multi')
-    except TypeError:
-        pass
-    if ps.pHtot is None:
-        ps.pHtot = np.array(cp(ps.H), ndmin=1)
-        # Calculate other pH scales
-        ps.update(calc_pH_scales(ps.pHtot, ps.pHfree, ps.pHsws,
-                                 ps.TS, ps.TF, ps.Ks))
-
-    # clean up for output
+        ps.BT = calc_TB(ps.S_in)
+    # elif isinstance(BT, (int, float)):
+    #     ps.BT = ps.BT * ps.S_in / 35.
+
+    # Calculate Ks at input conditions
+    ps.Ks = calc_Ks(ps.T_in, ps.S_in, ps.P_in,
+                    ps.Mg, ps.Ca, ps.TS, ps.TF, ps.Ks)
+
+    # Calculate pH scales at input conditions (does nothing if no pH given)
+    ps.update(calc_pH_scales(ps.pHtot, ps.pHfree, ps.pHsws, ps.pHNBS,
+                             ps.TS, ps.TF, ps.T_in + 273.15, ps.S_in, ps.Ks))
+
+    # calculate C system at input conditions
+    ps.update(calc_C_species(pHtot=ps.pHtot, DIC=ps.DIC, CO2=ps.CO2,
+                             HCO3=ps.HCO3, CO3=ps.CO3, TA=ps.TA,
+                             fCO2=ps.fCO2, pCO2=ps.pCO2,
+                             T_in=ps.T_in, BT=ps.BT, TP=ps.TP, TSi=ps.TSi,
+                             TS=ps.TS, TF=ps.TF, Ks=ps.Ks))
+
+    # clean up output
     for k in ['BT', 'CO2', 'CO3', 'Ca', 'DIC', 'H',
-              'HCO3', 'Mg', 'S', 'T', 'TA',
+              'HCO3', 'Mg', 'S_in', 'T_in', 'TA',
               'CAlk', 'PAlk', 'SiAlk', 'OH']:
         if not isinstance(ps[k], np.ndarray):
             # convert all outputs to (min) 1D numpy arrays.
             ps[k] = np.array(ps[k], ndmin=1)
     if ps.unit != 1:
-        for p in upar + ['CAlk', 'BAlk', 'PAlk', 'SiAlk', 'OH', 'HSO4', 'HF', 'Hfree']:
+        for p in upar + ['CAlk', 'BAlk', 'PAlk', 'SiAlk',
+                         'OH', 'HSO4', 'HF', 'Hfree']:
             ps[p] *= ps.unit  # convert back to input units
 
+    # Calculate Output Conditions
+    # ===========================
+    if ps.T_out is not None or ps.S_out is not None or ps.P_out is not None:
+        if ps.T_out is None:
+            ps.T_out = ps.T_in
+        if ps.S_out is None:
+            ps.S_out = ps.S_in
+        if ps.P_out is None:
+            ps.P_out = ps.P_in
+        # assumes conserved alkalinity and DIC
+        out_cond = Csys(TA=ps.TA, DIC=ps.DIC, 
+                        T_in=ps.T_out,
+                        S_in=ps.S_out,
+                        P_in=ps.P_out,
+                        unit=ps.unit)
+        # Calculate pH scales at output conditions (does nothing if no pH given)
+        out_cond.update(calc_pH_scales(out_cond.pHtot, out_cond.pHfree, out_cond.pHsws, out_cond.pHNBS,
+                                       out_cond.TS, out_cond.TF, out_cond.T_in + 273.15, out_cond.S_in, out_cond.Ks))
+
+        # rename parameters in output conditions
+        outputs = ['BAlk', 'BT', 'CAlk', 'CO2', 'CO3',
+                   'DIC', 'H', 'HCO3', 'HF',
+                   'HSO4', 'Hfree', 'Ks', 'OH',
+                   'PAlk', 'SiAlk', 'TA', 'TF', 'TP',
+                   'TS', 'TSi', 'fCO2', 'pCO2',
+                   'pHfree', 'pHsws', 'pHtot', 'pHNBS']
+
+        ps.update({k + '_out': out_cond[k] for k in outputs})
+
     # remove some superfluous outputs
     rem = ['pdict']
     for r in rem:
@@ -323,9 +238,10 @@ def Bsys(pHtot=None, BT=None, BO3=None, BO4=None,
          ABT=None, ABO3=None, ABO4=None,
          dBT=None, dBO3=None, dBO4=None,
          alphaB=None,
-         T=25., S=35., P=None,
+         T_in=25., S_in=35., P_in=None,
+         T_out=None, S_out=None, P_out=None,
          Ca=None, Mg=None,
-         pHsws=None, pHfree=None,
+         pHsws=None, pHfree=None, pHNBS=None,
          Ks=None, pdict=None):
     """
     Calculate the boron chemistry of seawater from a minimal parameter set.
@@ -382,52 +298,33 @@ def Bsys(pHtot=None, BT=None, BO3=None, BO4=None,
 
     # Conserved seawater chemistry
     if 'TS' not in ps:
-        ps.TS = calc_TS(ps.S)
+        ps.TS = calc_TS(ps.S_in)
     if 'TF' not in ps:
-        ps.TF = calc_TF(ps.S)
+        ps.TF = calc_TF(ps.S_in)
 
-    # if neither Ca nor Mg provided, use default Ks
-    ps.Ks = get_Ks(ps)
+    # Calculate Ks
+    ps.Ks = calc_Ks(ps.T_in, ps.S_in, ps.P_in,
+                    ps.Mg, ps.Ca, ps.TS, ps.TF, ps.Ks)
 
     # Calculate pH scales (does nothing if none pH given)
-    ps.update(calc_pH_scales(ps.pHtot, ps.pHfree, ps.pHsws,
-                             ps.TS, ps.TF, ps.Ks))
+    ps.update(calc_pH_scales(ps.pHtot, ps.pHfree, ps.pHsws, ps.pHNBS,
+                             ps.TS, ps.TF, ps.T_in + 273.15, ps.S_in, ps.Ks))
 
-    # B system calculations
-    if ps.pHtot is not None and ps.BT is not None:
-        ps.H = ch(ps.pHtot)
-    elif ps.BT is not None and ps.BO3 is not None:
-        ps.H = BT_BO3(ps.BT, ps.BO3, ps.Ks)
-    elif ps.BT is not None and ps.BO4 is not None:
-        ps.H = BT_BO4(ps.BT, ps.BO4, ps.Ks)
-    elif ps.BO3 is not None and ps.BO4 is not None:
-        ps.BT = ps.BO3 + ps.BO3
-        ps.H = BT_BO3(ps.BT, ps.BO3, ps.Ks)
-    elif ps.pHtot is not None and ps.BO3 is not None:
-        ps.H = ch(ps.pHtot)
-        ps.BT = pH_BO3(ps.pHtot, ps.BO3, ps.Ks)
-    elif ps.pHtot is not None and ps.BO4 is not None:
-        ps.H = ch(ps.pHtot)
-        ps.BT = pH_BO4(ps.pHtot, ps.BO4, ps.Ks)
-
-    # The above makes sure that BT and H are known,
-    # this next bit calculates all the missing species
-    # from BT and H.
-    if ps.BO3 is None:
-        ps.BO3 = cBO3(ps.BT, ps.H, ps.Ks)
-    if ps.BO4 is None:
-        ps.BO4 = cBO4(ps.BT, ps.H, ps.Ks)
+    ps.update(calc_B_species(pHtot=ps.pHtot, BT=ps.BT, BO3=ps.BO3, BO4=ps.BO4, Ks=ps.Ks))
+
+    # If pH not calced yet, calculate on all scales
     if ps.pHtot is None:
         ps.pHtot = np.array(cp(ps.H), ndmin=1)
         # Calculate other pH scales
-        ps.update(calc_pH_scales(ps.pHtot, ps.pHfree, ps.pHsws,
-                                 ps.TS, ps.TF, ps.Ks))
+        ps.update(calc_pH_scales(ps.pHtot, ps.pHfree, ps.pHsws, ps.pHNBS,
+                                 ps.TS, ps.TF, ps.T_in + 273.15, ps.S_in, ps.Ks))
 
+    # If any isotope parameter specified, calculate the isotope systen.
     if NnotNone(ps.ABT, ps.ABO3, ps.ABO4, ps.dBT, ps.dBO3, ps.dBO4) != 0:
         ps.update(ABsys(pdict=ps))
 
     for k in ['BT', 'H', 'BO3', 'BO4',
-              'Ca', 'Mg', 'S', 'T']:
+              'Ca', 'Mg', 'S_in', 'T_in']:
         # convert all outputs to (min) 1D numpy arrays.
         if not isinstance(ps[k], np.ndarray):
             # convert all outputs to (min) 1D numpy arrays.
@@ -448,9 +345,9 @@ def ABsys(pHtot=None,
           ABT=None, ABO3=None, ABO4=None,
           dBT=None, dBO3=None, dBO4=None,
           alphaB=None,
-          T=25., S=35., P=None,
+          T_in=25., S_in=35., P_in=None,
           Ca=None, Mg=None,
-          pHsws=None, pHfree=None,
+          pHsws=None, pHfree=None, pHNBS=None,
           Ks=None, pdict=None):
     """
     Calculate the boron isotope chemistry of seawater from a minimal parameter set.
@@ -512,19 +409,19 @@ def ABsys(pHtot=None,
     if isinstance(pdict, dict):
         ps.update(pdict)
 
-
     # Conserved seawater chemistry
     if 'TS' not in ps:
-        ps.TS = calc_TS(ps.S)
+        ps.TS = calc_TS(ps.S_in)
     if 'TF' not in ps:
-        ps.TF = calc_TF(ps.S)
+        ps.TF = calc_TF(ps.S_in)
 
-    # if neither Ca nor Mg provided, use default Ks
-    ps.Ks = get_Ks(ps)
+    # Calculate Ks
+    ps.Ks = calc_Ks(ps.T_in, ps.S_in, ps.P_in,
+                    ps.Mg, ps.Ca, ps.TS, ps.TF, ps.Ks)
 
     # Calculate pH scales (does nothing if none pH given)
-    ps.update(calc_pH_scales(ps.pHtot, ps.pHfree, ps.pHsws,
-                             ps.TS, ps.TF, ps.Ks))
+    ps.update(calc_pH_scales(ps.pHtot, ps.pHfree, ps.pHsws, ps.pHNBS,
+                             ps.TS, ps.TF, ps.T_in + 273.15, ps.S_in, ps.Ks))
 
     # if deltas provided, calculate corresponding As
     if ps.dBT is not None:
@@ -580,9 +477,10 @@ def CBsys(pHtot=None, DIC=None, CO2=None, HCO3=None, CO3=None, TA=None, fCO2=Non
           BT=None, BO3=None, BO4=None,
           ABT=None, ABO3=None, ABO4=None, dBT=None, dBO3=None, dBO4=None,
           alphaB=None,
-          T=25., S=35., P=None,
+          T_in=25., S_in=35., P_in=None,
+          T_out=None, S_out=None, P_out=None,
           Ca=None, Mg=None, TP=0., TSi=0.,
-          pHsws=None, pHfree=None,
+          pHsws=None, pHfree=None, pHNBS=None,
           Ks=None, pdict=None, unit='umol'):
     """
     Calculate carbon, boron and boron isotope chemistry of seawater from a minimal parameter set.
@@ -659,7 +557,6 @@ def CBsys(pHtot=None, DIC=None, CO2=None, HCO3=None, CO3=None, TA=None, fCO2=Non
     -------
     dict(/Bunch) containing all calculated parameters.
     """
-
     # Bunch inputs
     ps = Bunch(locals())
     if isinstance(pdict, dict):
@@ -683,109 +580,107 @@ def CBsys(pHtot=None, DIC=None, CO2=None, HCO3=None, CO3=None, TA=None, fCO2=Non
         if ps[p] is not None:
             ps[p] = np.divide(ps[p], ps.unit)  # convert to molar
 
-    # reassign unit, so conversions aren't repeated by Csys
+    # reassign unit, convert back at end
     orig_unit = ps.unit
     ps.unit = 1.
 
-    # determin max lengths
-    kexcl = ['Ks', 'pdict', 'unit']
-    ks = [k for k in ps.keys() if k not in kexcl]
-    L = maxL(*[ps[k] for k in ks])
-    # make inputs same length
-    for k in ks:
-        if ps[k] is not None:
-            if isinstance(ps[k], (int, float)):
-                ps[k] = np.full(L, ps[k])
-
     # Conserved seawater chemistry
     if 'TS' not in ps:
-        ps.TS = calc_TS(ps.S)
+        ps.TS = calc_TS(ps.S_in)
     if 'TF' not in ps:
-        ps.TF = calc_TF(ps.S)
+        ps.TF = calc_TF(ps.S_in)
 
     # Calculate Ks
-    ps.Ks = get_Ks(ps)
+    ps.Ks = calc_Ks(ps.T_in, ps.S_in, ps.P_in,
+                    ps.Mg, ps.Ca, ps.TS, ps.TF, ps.Ks)
 
     # Calculate pH scales (does nothing if none pH given)
-    ps.update(calc_pH_scales(ps.pHtot, ps.pHfree, ps.pHsws,
-                             ps.TS, ps.TF, ps.Ks))
+    ps.update(calc_pH_scales(ps.pHtot, ps.pHfree, ps.pHsws, ps.pHNBS,
+                             ps.TS, ps.TF, ps.T_in + 273.15, ps.S_in, ps.Ks))
 
     # if fCO2 is given but CO2 is not, calculate CO2
     if ps.CO2 is None:
         if ps.fCO2 is not None:
             ps.CO2 = fCO2_to_CO2(ps.fCO2, ps.Ks)
         elif ps.pCO2 is not None:
-            ps.CO2 = fCO2_to_CO2(pCO2_to_fCO2(ps.pCO2, ps.T), ps.Ks)
+            ps.CO2 = fCO2_to_CO2(pCO2_to_fCO2(ps.pCO2, ps.T_in), ps.Ks)
 
     # if no B info provided, assume modern conc.
     nBspec = NnotNone(ps.BT, ps.BO3, ps.BO4)
     if nBspec == 0:
-        ps.BT = calc_TB(ps.S)
+        ps.BT = calc_TB(ps.S_in)
     elif isinstance(BT, (int, float)):
-        ps.BT = ps.BT * ps.S / 35.
-
-    # This section works out the order that things should be calculated in.
+        ps.BT = ps.BT * ps.S_in / 35.
+    # count number of not None C parameters
+    nCspec = NnotNone(ps.DIC, ps.CO2, ps.HCO3, ps.CO3)  # used below
+
+    # if pH is given, it's easy
+    if ps.pHtot is not None or nBspec == 2:
+        ps.update(calc_B_species(pHtot=ps.pHtot, BT=ps.BT, BO3=ps.BO3, BO4=ps.BO4, Ks=ps.Ks))
+        ps.update(calc_C_species(pHtot=ps.pHtot, DIC=ps.DIC, CO2=ps.CO2,
+                                 HCO3=ps.HCO3, CO3=ps.CO3, TA=ps.TA,
+                                 fCO2=ps.fCO2, pCO2=ps.pCO2,
+                                 T_in=ps.T_in, BT=ps.BT, TP=ps.TP, TSi=ps.TSi,
+                                 TS=ps.TS, TF=ps.TF, Ks=ps.Ks))
+    # if not, this section works out the order that things should be calculated in.
     # Special case: if pH is missing, must have:
-    #   a) two C
-    #   b) two B
+    #   a) two C or one C and both TA and BT
+    #   b) two B (above)
     #   c) one pH-dependent B, one pH-dependent C... But that's cray...
     #      (c not implemented!)
-    if ps.pHtot is None:
-        nCspec = NnotNone(ps.DIC, ps.CO2, ps.HCO3, ps.CO3)
-        # a) if there are 2 C species, or one C species and TA and BT
-        if ((nCspec == 2) | ((nCspec == 1) & (NnotNone(ps.TA, ps.BT) == 2))):
-            ps.update(Csys(pdict=ps))  # calculate C first
-            ps.update(Bsys(pdict=ps))  # then B
-            # Note on the peculiar syntax here:
-            #  ps is a dict of parameters, where
-            #  everything that needs to be calculated
-            #  is None.
-            #  We give this to the [N]sys function
-            #  as pdict, which passes the paramters from THIS
-            #  function to [N]sys.
-            #  As the output of Csys is also a dict (Bunch)
-            #  with exactly the same form as ps, we can then
-            #  use the .update attribute of ps to update
-            #  all the paramters that were calculated by
-            #  Csyst.
-            #  Thus, all calculation is incremental, working
-            #  with the same parameter set. As dicts are
-            #  mutable, this has the added benefit of the
-            #  parameters only being stored in memory once.
-            if ps.TA is None:
-                (ps.TA, ps.CAlk, ps.BAlk,
-                 ps.PAlk, ps.SiAlk, ps.OH,
-                 ps.Hfree, ps.HSO4, ps.HF) = cTA(H=ps.H,
-                                                 DIC=ps.DIC,
-                                                 BT=ps.BT,
-                                                 TP=ps.TP,
-                                                 TSi=ps.TSi,
-                                                 TS=ps.TS,
-                                                 TF=ps.TF,
-                                                 Ks=ps.Ks, mode='multi')
-                # necessary becayse TA in Csys fails if there's no BT
-        # b) if there are 2 B species
-        elif nBspec == 2:
-            ps.update(Bsys(pdict=ps))  # calculate B first
-            ps.update(Csys(pdict=ps))  # then C
-        else:  # if neither condition is met, throw an error
-            raise ValueError(("Impossible! You haven't provided enough parameters.\n" +
-                              "If you don't know pH, you must provide either:\n" +
-                              "  - Two of [DIC, CO2, HCO3, CO3], and one of [BT, BO3, BO4]\n" +
-                              "  - One of [DIC, CO2, HCO3, CO3], and TA and BT\n" +
-                              "  - Two of [BT, BO3, BO4] and one of [DIC, CO2, HCO3, CO3]"))
-
-    else:  # if we DO have pH, it's dead easy!
-        ps.update(Bsys(pdict=ps))  # calculate B first
-        ps.update(Csys(pdict=ps))  # then C
+    elif ((nCspec == 2) | ((nCspec == 1) & (NnotNone(ps.TA, ps.BT) == 2))):  # case A
+        ps.update(calc_C_species(pHtot=ps.pHtot, DIC=ps.DIC, CO2=ps.CO2,
+                                 HCO3=ps.HCO3, CO3=ps.CO3, TA=ps.TA,
+                                 fCO2=ps.fCO2, pCO2=ps.pCO2,
+                                 T_in=ps.T_in, BT=ps.BT, TP=ps.TP, TSi=ps.TSi,
+                                 TS=ps.TS, TF=ps.TF, Ks=ps.Ks))
+        ps.update(calc_B_species(pHtot=ps.pHtot, BT=ps.BT, BO3=ps.BO3, BO4=ps.BO4, Ks=ps.Ks))
+    # elif nBspec == 2:  # case B -- moved up
+    #     ps.update(calc_B_species(pHtot=ps.pHtot, BT=ps.BT, BO3=ps.BO3, BO4=ps.BO4, Ks=ps.Ks))
+    #     ps.update(calc_C_species(pHtot=ps.pHtot, DIC=ps.DIC, CO2=ps.CO2,
+    #                              HCO3=ps.HCO3, CO3=ps.CO3, TA=ps.TA,
+    #                              fCO2=ps.fCO2, pCO2=ps.pCO2,
+    #                              T_in=ps.T_in, BT=ps.BT, TP=ps.TP, TSi=ps.TSi,
+    #                              TS=ps.TS, TF=ps.TF, Ks=ps.Ks))  # then C
+    else:  # if neither condition is met, throw an error
+        raise ValueError(("Impossible! You haven't provided enough parameters.\n" +
+                          "If you don't know pH, you must provide either:\n" +
+                          "  - Two of [DIC, CO2, HCO3, CO3], and one of [BT, BO3, BO4]\n" +
+                          "  - One of [DIC, CO2, HCO3, CO3], and TA and BT\n" +
+                          "  - Two of [BT, BO3, BO4] and one of [DIC, CO2, HCO3, CO3]"))
 
     for p in upar + ['CAlk', 'BAlk', 'PAlk', 'SiAlk', 'OH', 'HSO4', 'HF', 'Hfree']:
         ps[p] *= orig_unit  # convert back to input units
 
+    # Recursive approach to calculate output params.
+    # if output conditions specified, calculate outputs.
+    if ps.T_out is not None or ps.S_out is not None or ps.P_out is not None:
+        if ps.T_out is None:
+            ps.T_out = ps.T_in
+        if ps.S_out is None:
+            ps.S_out = ps.S_in
+        if ps.P_out is None:
+            ps.P_out = ps.P_in
+        # assumes conserved alkalinity
+        out_cond = CBsys(TA=ps.TA, DIC=ps.DIC, BT=ps.BT, T_in=ps.T_out,
+                         S_in=ps.S_out, P_in=ps.P_out, unit=ps.unit)
+        # Calculate pH scales (does nothing if no pH given)
+        out_cond.update(calc_pH_scales(out_cond.pHtot, out_cond.pHfree, out_cond.pHsws, out_cond.pHNBS,
+                                       out_cond.TS, out_cond.TF, out_cond.T_in + 273.15, out_cond.S_in, out_cond.Ks))
+        # rename parameters in output conditions
+        outputs = ['BAlk', 'BT', 'CAlk', 'CO2', 'CO3',
+                   'DIC', 'H', 'HCO3', 'HF',
+                   'HSO4', 'Hfree', 'Ks', 'OH',
+                   'PAlk', 'SiAlk', 'TA', 'TF', 'TP',
+                   'TS', 'TSi', 'fCO2', 'pCO2',
+                   'pHfree', 'pHsws', 'pHtot', 'pHNBS', 'BO3', 'BO4',
+                   'ABO3', 'ABO4', 'dBO3', 'dBO4']
+
+        ps.update({k + '_out': out_cond[k] for k in outputs})
+
         # remove some superfluous outputs
     rem = ['pdict', 'unit']
     for r in rem:
         if r in ps:
             del ps[r]
-
     return ps
diff --git a/cbsyst/helpers.py b/cbsyst/helpers.py
index 90b4fb0..6a9d815 100644
--- a/cbsyst/helpers.py
+++ b/cbsyst/helpers.py
@@ -248,14 +248,27 @@ def calc_TF(Sal):
     return (a / b) * (Sal / c)  # mol/kg-SW
 
 
+# def calc_TB(Sal):
+#     """
+#     Calculate total Boron
+
+#     Lee, Kim, Byrne, Millero, Feely, Yong-Ming Liu. 2010.
+#     Geochimica Et Cosmochimica Acta 74 (6): 1801-1811
+#     """
+#     a, b = (0.0004326, 35.)
+#     return a * Sal / b
+
+
 def calc_TB(Sal):
     """
     Calculate total Boron
-
-    Lee, Kim, Byrne, Millero, Feely, Yong-Ming Liu. 2010.
-    Geochimica Et Cosmochimica Acta 74 (6): 1801-1811
+    
+    Directly from CO2SYS:
+    Uppstrom, L., Deep-Sea Research 21:161-162, 1974:
+    this is 0.000416 * Sal/35. = 0.0000119 * Sal
+    TB(FF) = (0.000232 / 10.811) * (Sal / 1.80655) in mol/kg-SW
     """
-    a, b = (0.0004326, 35.)
+    a, b = (0.0004157, 35.)
     return a * Sal / b
 
 
@@ -269,28 +282,36 @@ def calc_fH(TempK, Sal):
 
 
 # Convert between pH scales
-def calc_pH_scales(pHtot, pHfree, pHsws, TS, TF, Ks):
+def calc_pH_scales(pHtot, pHfree, pHsws, pHNBS, TS, TF, TempK, Sal, Ks):
     """
     Calculate pH on all scales, given one.
     """
 
     # check if any pH scale is given.
-    npH = NnotNone(pHfree, pHsws, pHtot)
+    npH = NnotNone(pHfree, pHsws, pHtot, pHNBS)
 
     if npH == 1:
         # pH scale conversions
         FREEtoTOT = -np.log10((1 + TS / Ks.KSO4))
         SWStoTOT = -np.log10((1 + TS / Ks.KSO4) /
                              (1 + TS / Ks.KSO4 + TF / Ks.KF))
+        fH = calc_fH(TempK, Sal)
 
         if pHtot is not None:
             return {'pHfree': pHtot - FREEtoTOT,
-                    'pHsws': pHtot - SWStoTOT}
+                    'pHsws': pHtot - SWStoTOT,
+                    'pHNBS': pHtot - SWStoTOT - np.log10(fH)}
         elif pHsws is not None:
             return {'pHfree': pHsws + SWStoTOT - FREEtoTOT,
-                    'pHtot': pHsws + SWStoTOT}
+                    'pHtot': pHsws + SWStoTOT,
+                    'pHNBS': pHsws - np.log10(fH)}
         elif pHfree is not None:
             return {'pHsws': pHfree + FREEtoTOT - SWStoTOT,
-                    'pHtot': pHfree + FREEtoTOT}
+                    'pHtot': pHfree + FREEtoTOT,
+                    'pHNBS': pHfree + FREEtoTOT - SWStoTOT - np.log10(fH)}
+        elif pHNBS is not None:
+            return {'pHsws': pHNBS + np.log10(fH),
+                    'pHtot': pHNBS + np.log10(fH) + SWStoTOT,
+                    'pHfree': pHNBS + np.log10(fH) + SWStoTOT - FREEtoTOT}
     else:
         return {}
diff --git a/setup.py b/setup.py
index 3be183c..8508ee5 100644
--- a/setup.py
+++ b/setup.py
@@ -25,6 +25,5 @@ setup(name='cbsyst',
                         'pandas',
                         'uncertainties',
                         'tqdm'],
-      package_data={'cbsyst': ['test_data/GLODAP_data/*.py',
-                               'test_data/Lueker2000/*.csv']},
+      package_data={'cbsyst': ['test_data/*']},
       zip_safe=True)


Test patch:
diff --git a/cbsyst/test_data/Orr2015/plot_CO2SYS_comparison.py b/cbsyst/test_data/Orr2015/plot_CO2SYS_comparison.py
new file mode 100644
index 0000000..1c049d7
--- /dev/null
+++ b/cbsyst/test_data/Orr2015/plot_CO2SYS_comparison.py
@@ -0,0 +1,18 @@
+"""
+TODO: Script to compare CO2SYS.m output to cbsyst output.
+Comparison similar to Orr et al (2015)
+
+Parameters
+----------
+raw_data : str
+    Path to GLODAP-WOA comparison data.
+co2sys_data : str
+    Path to output from CO2SYS.m
+output_dir : str
+    Output directory for plots
+
+Returns
+-------
+Prints a range of statistics comparing CO2SYS and cbsyst.
+Creates a variety of plots with statistics in output_dir.
+"""
diff --git a/test_cbsyst.py b/test_cbsyst.py
index 06216fe..4e464e8 100644
--- a/test_cbsyst.py
+++ b/test_cbsyst.py
@@ -323,21 +323,21 @@ class ReferenceDataTestCase(unittest.TestCase):
 
         # Csys calculations
         # TA from pH and DIC
-        cTA = Csys(pHtot=pH, DIC=DIC, BT=BT, S=S)
+        cTA = Csys(pHtot=pH, DIC=DIC, BT=BT, S_in=S)
         # Calculate % differences from measured
         dTA = (100 * (TA - cTA.TA) / TA)
 
         self.assertLess(max(abs(dTA)), 0.2, msg='TA from DIC and pH')
 
         # pH from TA and DIC
-        cpH = Csys(DIC=DIC, TA=TA, BT=BT, S=S)
+        cpH = Csys(DIC=DIC, TA=TA, BT=BT, S_in=S)
         # Calculate % differences from measured
         dpH = (100 * (pH - cpH.pHtot) / pH)
 
         self.assertLess(max(abs(dpH)), 0.2, msg='pH from TA and DIC')
 
         # DIC from pH and TA
-        cDIC = Csys(pHtot=pH, TA=TA, BT=BT, S=S)
+        cDIC = Csys(pHtot=pH, TA=TA, BT=BT, S_in=S)
         # Calculate % differences from measured
         dDIC = (100 * (DIC - cDIC.DIC) / DIC)
 
@@ -365,21 +365,21 @@ class ReferenceDataTestCase(unittest.TestCase):
 
         # Csys calculations
         # TA from pH and DIC
-        cTA = CBsys(pHtot=pH, DIC=DIC, BT=BT, S=S)
+        cTA = CBsys(pHtot=pH, DIC=DIC, BT=BT, S_in=S)
         # Calculate % differences from measured
         dTA = (100 * (TA - cTA.TA) / TA)
 
         self.assertLess(max(abs(dTA)), 0.2, msg='TA from DIC and pH')
 
         # pH from TA and DIC
-        cpH = CBsys(DIC=DIC, TA=TA, BT=BT, S=S)
+        cpH = CBsys(DIC=DIC, TA=TA, BT=BT, S_in=S)
         # Calculate % differences from measured
         dpH = (100 * (pH - cpH.pHtot) / pH)
 
         self.assertLess(max(abs(dpH)), 0.2, msg='pH from TA and DIC')
 
         # DIC from pH and TA
-        cDIC = CBsys(pHtot=pH, TA=TA, BT=BT, S=S)
+        cDIC = CBsys(pHtot=pH, TA=TA, BT=BT, S_in=S)
         # Calculate % differences from measured
         dDIC = (100 * (DIC - cDIC.DIC) / DIC)
 
@@ -395,7 +395,7 @@ class ReferenceDataTestCase(unittest.TestCase):
 
         # Calculate using cbsys
         # TA from DIC and fCO2
-        cTA = Csys(DIC=ld.DIC.values, fCO2=ld.fCO2.values, T=ld.Temp.values, S=ld.Sal.values)
+        cTA = Csys(DIC=ld.DIC.values, fCO2=ld.fCO2.values, T_in=ld.Temp.values, S_in=ld.Sal.values)
         dTA = ld.TA - cTA.TA
         dTA_median = np.median(dTA)
         dTA_pc95 = np.percentile(dTA, [2.5, 97.5])
@@ -403,7 +403,7 @@ class ReferenceDataTestCase(unittest.TestCase):
         self.assertTrue(all(abs(dTA_pc95 - dTA_median) <= 16), msg='TA 95% Conf <= 16')
 
         # fCO2 from TA and DIC
-        cfCO2 = Csys(TA=ld.TA.values, DIC=ld.DIC.values, T=ld.Temp.values, S=ld.Sal.values)
+        cfCO2 = Csys(TA=ld.TA.values, DIC=ld.DIC.values, T_in=ld.Temp.values, S_in=ld.Sal.values)
         dfCO2 = ld.fCO2 - cfCO2.fCO2
         dfCO2_median = np.median(dfCO2)
         # dfCO2_pc95 = np.percentile(dfCO2, [2.5, 97.5])
@@ -414,7 +414,7 @@ class ReferenceDataTestCase(unittest.TestCase):
         # self.assertTrue(all(abs(dfCO2_pc95) <= 70), msg='fCO2 95% Conc <= 70')
 
         # DIC from TA and fCO2
-        cDIC = Csys(TA=ld.TA.values, fCO2=ld.fCO2.values, T=ld.Temp.values, S=ld.Sal.values)
+        cDIC = Csys(TA=ld.TA.values, fCO2=ld.fCO2.values, T_in=ld.Temp.values, S_in=ld.Sal.values)
         dDIC = ld.DIC - cDIC.DIC
         dDIC_median = np.median(dDIC)
         dDIC_pc95 = np.percentile(dDIC, [2.5, 97.5])
@@ -431,7 +431,7 @@ class ReferenceDataTestCase(unittest.TestCase):
 
         # Calculate using cbsys
         # TA from DIC and fCO2
-        cTA = CBsys(DIC=ld.DIC.values, fCO2=ld.fCO2.values, T=ld.Temp.values, S=ld.Sal.values)
+        cTA = CBsys(DIC=ld.DIC.values, fCO2=ld.fCO2.values, T_in=ld.Temp.values, S_in=ld.Sal.values)
         dTA = ld.TA - cTA.TA
         dTA_median = np.median(dTA)
         dTA_pc95 = np.percentile(dTA, [2.5, 97.5])
@@ -439,7 +439,7 @@ class ReferenceDataTestCase(unittest.TestCase):
         self.assertTrue(all(abs(dTA_pc95 - dTA_median) <= 16), msg='TA 95% Conf <= 16')
 
         # fCO2 from TA and DIC
-        cfCO2 = CBsys(TA=ld.TA.values, DIC=ld.DIC.values, T=ld.Temp.values, S=ld.Sal.values)
+        cfCO2 = CBsys(TA=ld.TA.values, DIC=ld.DIC.values, T_in=ld.Temp.values, S_in=ld.Sal.values)
         dfCO2 = ld.fCO2 - cfCO2.fCO2
         dfCO2_median = np.median(dfCO2)
         # dfCO2_pc95 = np.percentile(dfCO2, [2.5, 97.5])
@@ -450,7 +450,7 @@ class ReferenceDataTestCase(unittest.TestCase):
         # self.assertTrue(all(abs(dfCO2_pc95) <= 70), msg='fCO2 95% Conc <= 70')
 
         # DIC from TA and fCO2
-        cDIC = CBsys(TA=ld.TA.values, fCO2=ld.fCO2.values, T=ld.Temp.values, S=ld.Sal.values)
+        cDIC = CBsys(TA=ld.TA.values, fCO2=ld.fCO2.values, T_in=ld.Temp.values, S_in=ld.Sal.values)
         dDIC = ld.DIC - cDIC.DIC
         dDIC_median = np.median(dDIC)
         dDIC_pc95 = np.percentile(dDIC, [2.5, 97.5])
@@ -477,8 +477,8 @@ class ReferenceDataTestCase(unittest.TestCase):
         gd = gd.loc[gd.cruise != 270]
 
         # calculate pH from TA and DIC
-        cpH = Csys(TA=gd.talk, DIC=gd.tco2, T=gd.temperature, S=gd.salinity,
-                   P=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
+        cpH = Csys(TA=gd.talk, DIC=gd.tco2, T_in=gd.temperature, S_in=gd.salinity,
+                   P_in=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
         pH_resid = gd.phtsinsitutp - cpH.pHtot
         pH_median = np.median(pH_resid)
         pH_pc95 = np.percentile(pH_resid, [2.5, 97.5])
@@ -487,8 +487,8 @@ class ReferenceDataTestCase(unittest.TestCase):
         self.assertTrue(all(abs(pH_pc95) <= 0.05), msg='pH 95% Conf <= 0.05')
 
         # calculate TA from pH and DIC
-        cTA = Csys(pHtot=gd.phtsinsitutp, DIC=gd.tco2, T=gd.temperature, S=gd.salinity,
-                   P=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
+        cTA = Csys(pHtot=gd.phtsinsitutp, DIC=gd.tco2, T_in=gd.temperature, S_in=gd.salinity,
+                   P_in=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
         TA_resid = gd.talk - cTA.TA
         TA_median = np.median(TA_resid)
         TA_pc95 = np.percentile(TA_resid, [2.5, 97.5])
@@ -497,8 +497,8 @@ class ReferenceDataTestCase(unittest.TestCase):
         self.assertTrue(all(abs(TA_pc95) < 13), msg='TA 95% Conf <= 15')
 
         # calculate DIC from TA and pH
-        cDIC = Csys(pHtot=gd.phtsinsitutp, TA=gd.talk, T=gd.temperature, S=gd.salinity,
-                    P=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
+        cDIC = Csys(pHtot=gd.phtsinsitutp, TA=gd.talk, T_in=gd.temperature, S_in=gd.salinity,
+                    P_in=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
         DIC_resid = gd.tco2 - cDIC.DIC
         DIC_median = np.median(DIC_resid)
         DIC_pc95 = np.percentile(DIC_resid, [2.5, 97.5])
@@ -526,8 +526,8 @@ class ReferenceDataTestCase(unittest.TestCase):
         gd = gd.loc[gd.cruise != 270]
 
         # calculate pH from TA and DIC
-        cpH = CBsys(TA=gd.talk, DIC=gd.tco2, T=gd.temperature, S=gd.salinity,
-                    P=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
+        cpH = CBsys(TA=gd.talk, DIC=gd.tco2, T_in=gd.temperature, S_in=gd.salinity,
+                    P_in=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
         pH_resid = gd.phtsinsitutp - cpH.pHtot
         pH_median = np.median(pH_resid)
         pH_pc95 = np.percentile(pH_resid, [2.5, 97.5])
@@ -536,8 +536,8 @@ class ReferenceDataTestCase(unittest.TestCase):
         self.assertTrue(all(abs(pH_pc95) <= 0.05), msg='pH 95% Conf <= 0.05')
 
         # calculate TA from pH and DIC
-        cTA = CBsys(pHtot=gd.phtsinsitutp, DIC=gd.tco2, T=gd.temperature, S=gd.salinity,
-                    P=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
+        cTA = CBsys(pHtot=gd.phtsinsitutp, DIC=gd.tco2, T_in=gd.temperature, S_in=gd.salinity,
+                    P_in=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
         TA_resid = gd.talk - cTA.TA
         TA_median = np.median(TA_resid)
         TA_pc95 = np.percentile(TA_resid, [2.5, 97.5])
@@ -546,8 +546,8 @@ class ReferenceDataTestCase(unittest.TestCase):
         self.assertTrue(all(abs(TA_pc95) < 13), msg='TA 95% Conf <= 15')
 
         # calculate DIC from TA and pH
-        cDIC = CBsys(pHtot=gd.phtsinsitutp, TA=gd.talk, T=gd.temperature, S=gd.salinity,
-                     P=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
+        cDIC = CBsys(pHtot=gd.phtsinsitutp, TA=gd.talk, T_in=gd.temperature, S_in=gd.salinity,
+                     P_in=gd.pressure, TP=gd.phosphate, TSi=gd.silicate, BT=415.7)
         DIC_resid = gd.tco2 - cDIC.DIC
         DIC_median = np.median(DIC_resid)
         DIC_pc95 = np.percentile(DIC_resid, [2.5, 97.5])