-
Notifications
You must be signed in to change notification settings - Fork 2
/
bgc_isocap.m
64 lines (51 loc) · 2.31 KB
/
bgc_isocap.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
function [Q,Qx,Psi] = bgc_isocap(pco2,dic,temp,sal)
%bgc_isocap Calculate isocapnic quotient from input conditions
% Requires functions from https://github.com/mvdh7/oceancarb-constants
% === INPUTS ==============================================================
% pco2: seawater CO2 partial pressure (pCO2) / microatm
% dic: dissolved inorganic carbon (DIC) / micromol/kg
% temp: seawater temperature / degrees C
% sal: practical salinity /
% === OUTPUTS =============================================================
% Q: isocapnic quotient /
% Qx: approximation of isocapnic quotient /
% Psi: Phi for calcification (i.e. RA_RC = [-2 -1]) /
% =========================================================================
% --- Written by Matthew P. Humphreys -------------------------------------
% v1: last updated 2017-06-29
% Citation: Humphreys MP, Daniels CJ, Wolf-Gladrow DA, Tyrrell T, &
% Achterberg EP (2017): "On the influence of marine
% biogeochemical processes over CO2 exchange between the
% atmosphere and ocean", Marine Chemistry, in press,
% doi:10.1016/j.marchem.2017.12.006
% =========================================================================
% Absolute zero
azero = 273.15; % degC
% Evaluate coefficients
[~,pkHSO4] = cD90(temp+azero,sal);
TSO4 = cMR66(sal);
kHSO4 = 10.^-pkHSO4;
t2f = log10(1 + TSO4./kHSO4);
k0 = cW74(temp,sal);
[pKstarT_C1,pKstarT_C2] = cLDK00(temp+azero,sal);
k1 = 10.^-(pKstarT_C1 + t2f);
k2 = 10.^-(pKstarT_C2 + t2f);
tB = cLKB10(sal) * 1e-6;
pKstarT_w = cDSC07(temp+azero,sal);
kw = 10.^-(pKstarT_w + t2f);
pKstarT_B = cD90(temp+azero,sal);
kB = 10.^-(pKstarT_B + t2f); % mol/kg
% Aqueous CO2(aq) concentration
co2aq = pco2.*k0*1e-6;
% Hydrogen ion concentration
h = (-co2aq.*k1 - sqrt((co2aq.*k1).^2 ...
- 4*(co2aq - dic*1e-6).*co2aq.*k1.*k2)) ./ (2*(co2aq - dic*1e-6));
% Isocap slope <Q>
Q = ((k1.*co2aq.*h + 4*k1.*k2.*co2aq + kw.*h + h.^3) .* (kB + h).^2 ...
+ kB.*tB.*h.^3) ./ (k1.*co2aq .* (2*k2 + h) .* (kB + h).^2);
% Approximate isocap slope <Qx>
k2_k0k1 = k2 ./ (k0 .* k1);
Qx = 1 + 2 * k2_k0k1 .* dic ./ pco2;
% Psi (Frankignoulle et al., 1994) <Psi>
Psi = -1 + 2./Q;
end %function bgc_isocap