Function pH_total2nbs(S, T_C, pH_total)
' CONVERT pH (ON TOTAL SCALE) TO pH on NBS 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 T_C to degree K
T = T_C + 273.15
' SECOND: convert pH_total to pH_sws
' (CONVERT pH on TOTAL SCALE TO pH on SEAWATER SCALE)
H_T = 10^(- pH_total)
Hi = 10^(- (pH_total - 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_SWS = H_T + HF
' THIRD: convert pH_SWS to pH_NBS
' (CONVERT pH on SEAWATER SCALE TO pH on NBS SCALE)
'
' This is a two-step process: first adjust concentration units,
' then convert conc to activity...
'
' 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)
H_L = H_SWS * rho_SW
' Now, we can convert H+ concentration to activity, aH
' aH = 10^(-pH_NBS) = fH*H_L
' 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 H+ concentration to pH_NBS
fH = 1.2948 - .002036 * T
fH = fH + (.0004607 - 1.475E-06 * T) * S^2
aH = H_L * fH
pH_total2nbs = - log(aH) / log(10)
End Function