Function pH_nbs2total(S, T_C, pH_NBS)

' CONVERT pH (ON NBS SCALE) TO pH on TOTAL SCALE
'
' This code was adapted from the work of others:
'     CO2SYS.exe was originally written by Ernie Lewis and Doug Wallace.
'     It was converted to Matlab by Richard E. Zeebe & Dieter A. Wolf-Gladrow.
'     Parts were extracted and modified by Rachel M. Dunk.
'
' I am indebted to all of them for their initial development of this code.
'
' FIRST: convert pH_NBS to pH_SWS
'  aH = 10^(-pH_NBS) = fH*H_SWS
'  fH = activity coefficient of H+ ion & includes liquid junction effects 
'  The fit used in CO2SYS is valid for salinity 20-40 (from Takahashi et al,
'   Chapter 3 in GEOSECS Pacific Expedition,v.3, 1982, p. 80).
' For future ref, CO2SYS assumes fH is independent of P

' Convert T to degree K
    T = T_C + 273.15

' Convert pH to H+ concentration
    aH = 10 ^ -pH_NBS
    fH = 1.2948 - .002036 * T
    fH = fH + (.0004607 - 1.475E-06 * T) * S^2
    H  = aH / fH

' On NBS scale, aH is in mol/dm3_H2O, whereas H_SWS is in mol/kg_SW
'   to convert units we need the density of SW at given S,T
' Call DENSITY.bas (author=ETP), calculates rho(S,T,P=1atm) in g/cm3
'   NOTE: in this function T is in deg C
'   Equation of State is from Millero & Poisson (1981) DSR V28: 625-629.

    rho_SW = DENSITY(S, T_C)

    pH_SWS = -log(H / rho_SW) / log(10)

'SECOND: convert pH_SWS to pH_total
' CONVERT pH (on SEAWATER SCALE) TO pH on TOTAL SCALE

    H_SWS = 10^(- pH_SWS)
    Hi = 10^(- (pH_SWS + 0.01))
    F_T = 0.00007 * S/35
    lnK_F = 874/T - 9.68 + 0.111*S^0.5
    K_F = exp(lnK_F)
    HF = F_T / (1 + K_F/Hi)
    H_T = H_SWS - HF

    pH_nbs2total = - log(H_T) / log(10)

End Function