12/12 and 11/9/2015. Emilio Mayorga. Reproduce and extend Examples 1-3 from the mocys Python documentation page.
In [1]:
import mocsy
For DATA input: DIC and ALk in mol/kg, in situ temperature, pressure.
In [2]:
pH,pco2,fco2,co2,hco3,co3,OmegaA,OmegaC,BetaD,DENis,p,Tis = \
mocsy.mvars(temp=18, sal=35, alk=2300.e-6, dic=2000.e-6, sil=0, phos=0,
patm=1, depth=100, lat=0,
optcon='mol/kg', optt='Tinsitu', optp='db',
optb="u74", optk1k2='l', optkf="dg", optgas='Pinsitu')
print(pH,pco2,fco2,co2,hco3,co3,OmegaA,OmegaC,BetaD,DENis,p,Tis)
Compute the CONSTANTS with the same scalar input for S, T, and P (as above) but change to the newer options published since the best-practices guide: Lee et al. (2010) for total boron and Millero (2010) for K1 and K2:
In [3]:
Kh,K1,K2,Kb,Kw,Ks,Kf,Kspc,Kspa,K1p,K2p,K3p,Ksi,St,Ft,Bt = \
mocsy.mconstants(temp=18, sal=35, patm=1, depth=0,lat=0,
optt='Tinsitu', optp='m',
optb="l10", optk1k2='m10', optkf="dg", optgas='Ppot')
print(Kh,K1,K2,Kb,Kw,Ks,Kf,Kspc,Kspa,K1p,K2p,K3p,Ksi,St,Ft,Bt)
For MODEL input: DIC and Alk in mol/m3, potential temperature, and depth. Also do NOT specify last 4 options, thus using defaults (optb='l10', optk1k2='l', optkf='pf', optgas='Pinsitu').
In [4]:
pH,pco2,fco2,co2,hco3,co3,OmegaA,OmegaC,BetaD,DENis,p,Tis = \
mocsy.mvars(temp=18, sal=35, alk=2300*1028e-6, dic=2000*1028e-6,
sil=0, phos=0, patm=1, depth=100, lat=0,
optcon='mol/m3', optt='Tpot', optp='m')
print(pH,pco2,fco2,co2,hco3,co3,OmegaA,OmegaC,BetaD,DENis,p,Tis)
In [5]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
In [6]:
# Make dummy array with 11 members
one = np.ones(11, dtype='float32')
sal = 35*one
temp = 2*one
patm = 1*one
depth = np.arange(0, 11000, 1000, dtype='float32') # units in 'm'
In [7]:
# Compute in situ density (3 steps using 2 mocsy routines):
# a) Specify latitude = 60°S for depth to pressure conversion formula
lat = -60*one
# b) Compute pressure (db) from depth (m) and latitude following Saunders (1981)
pfd = mocsy.mdepth2press(depth,lat) # units in 'db'
# c) Compute in situ density (kg/m3) from salinity (psu), in-situ temp (C), and pressure (db)
rhois = mocsy.mrhoinsitu(sal, temp, pfd)
In [8]:
# Specify surface concentrations typical of 60°S (from GLODAP and WOA2009)
alk = 2295*one # umol/kg
dic = 2154*one # umol/kg
sio2 = 50.*one # umol/L
po4 = 1.8*one # umol/L
In [9]:
# Convert nutrient units(mol/L) to model units (mol/m3)
sio2 = sio2 * 1.0e-3
po4 = po4 * 1.0e-3
# Convert Alk and DIC units (umol/kg) to model units (mol/m3)
dic = dic * rhois * 1e-6
alk = alk * rhois * 1e-6
In [10]:
# Compute carbonate system variables, using 'model' options.
# Note that original example code had an error:
# patm argument was passed as scalar rather than array shaped like one
response_tup = mocsy.mvars(temp=temp, sal=sal, alk=alk, dic=dic,
sil=sio2, phos=po4, patm=patm, depth=depth, lat=lat,
optcon='mol/m3', optt='Tpot', optp='m')
pH,pco2,fco2,co2,hco3,co3,OmegaA,OmegaC,BetaD,DENis,p,Tis = response_tup
In [11]:
depth, pH, OmegaA
Out[11]:
In [12]:
data_for_pd = dict(depth=depth,
temp=temp, sal=sal, alk=alk, dic=dic, sil=sio2, phos=po4,
pH=pH, pCO2=pco2, fCO2=fco2, CO3=co3, OmegaA=OmegaA, OmegaC=OmegaC
)
In [13]:
data_df = pd.DataFrame(data_for_pd,
columns=['depth', 'temp', 'sal', 'alk', 'dic', 'sil', 'phos',
'pH', 'pCO2', 'fCO2', 'CO3', 'OmegaA', 'OmegaC']
)
data_df.set_index('depth', drop=False, inplace=True)
In [14]:
data_df
Out[14]:
In [15]:
data_df['OmegaA'].plot(title='Omega Aragonite (OmegaA) vs depth');
In [16]:
fig, ax1 = plt.subplots(figsize=(5,7))
ax1.plot(data_df['OmegaA'], data_df['depth'], lw=2, color="blue")
ax1.set_xlabel(r"$\Omega_A (fraction)$", fontsize=16, color="blue")
for label in ax1.get_xticklabels():
label.set_color("blue")
ax2 = ax1.twiny()
ax2.plot(data_df['pCO2'], data_df['depth'], lw=2, color="red")
ax2.set_xlabel(r"$pCO_2 (\mu atm)$", fontsize=16, color="red")
for label in ax2.get_xticklabels():
label.set_color("red")
ax1.invert_yaxis()
ax1.set_ylabel(r"$depth (m)$", fontsize=16);
Note that pCO2 here is calculated for "in-situ" conditions, a capability that apparently is fairly unique to mocsy. See
Orr, J. C. and Epitalon, J.-M.: Improved routines to model the ocean carbonate system: mocsy 2.0, Geosci. Model Dev., 8, 485-499, doi:10.5194/gmd-8-485-2015, 2015
In [17]:
# Extra Packages used.
# import numpy as np # Already imported above
# import pandas as pd # Already imported above
import urllib
import StringIO
In [18]:
# Read input .csv file from mocsy github repository
infileurl = "https://raw.githubusercontent.com/jamesorr/mocsy/master/DIC-Alk-P_vary_input.csv"
readfileurl = urllib.urlopen(infileurl).read()
infile = StringIO.StringIO(readfileurl)
In [19]:
# Read .csv input file with 7 columns: flag, Alk(mol/kg), DIC(mol/kg), salinity(psu),
# temp(C), press(db), PO43-(mol/kg), SiO2(mol/kg)
# infile = "DIC-Alk-P_vary_input.csv"
indata = np.loadtxt(infile, delimiter=',', dtype='float32', skiprows=1)
In [20]:
n = len(indata)
alk = indata[:,1]
dic = indata[:,2]
sal = indata[:,3]
temp = indata[:,4]
depth = indata[:,5] * 10 #Convert from bars to decibars
phos = indata[:,6]
sil = indata[:,7]
# Specify that latitude = 45°N for all samples (for depth to pressure conversion)
# note: "lat" is required input for mocsy, but results are only weakly sensitive
lat = np.full_like(depth, 45)
# Specify the atmospheric pressure for all samples (1 atm)
patm = np.full_like(depth, 1)
In [21]:
# Recover computed carbonate system variables
response_tup = mocsy.mvars(temp=temp, sal=sal, alk=alk, dic=dic,
sil=sil, phos=phos, patm=patm, depth=depth, lat=lat,
optcon='mol/kg', optt='Tinsitu', optp='db',
optb="u74", optk1k2='l', optkf="dg", optgas='Pzero')
pH, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis = response_tup
Notes for the above calculation:
optgas='Pzero', to expose the optgas option and match the pre-existing output on https://github.com/jamesorr/mocsy/blob/master/DIC-Alk-P-vary-from-mocsy.csv. From the mvars module in the Fortran code, regarding Pzero: "'zero order' fCO2 and pCO2 (typical approach, which is flawed) considers in situ T & only atm pressure (hydrostatic=0)"
In [22]:
data_for_pd = dict(depth=depth,
temp=temp, sal=sal, alk=alk, dic=dic, sil=sil, phos=phos,
pH=pH, pCO2=pco2, fCO2=fco2, CO3=co3, OmegaA=OmegaA, OmegaC=OmegaC
)
data_df = pd.DataFrame(data_for_pd,
columns=['depth', 'temp', 'sal', 'alk', 'dic', 'sil', 'phos',
'pH', 'pCO2', 'fCO2', 'CO3', 'OmegaA', 'OmegaC']
)
data_df.set_index('depth', drop=False, inplace=True)
data_df
Out[22]:
This output is a good match to the output csv file on the mocsy repository.
In [23]:
fig, ax1 = plt.subplots(figsize=(5,7))
ax1.plot(data_df['OmegaA'], data_df['depth'], lw=2, color="blue")
ax1.set_xlabel(r"$\Omega_A (fraction)$", fontsize=16, color="blue")
for label in ax1.get_xticklabels():
label.set_color("blue")
ax2 = ax1.twiny()
ax2.plot(data_df['pH'], data_df['depth'], lw=2, color="red")
ax2.set_xlabel(r"$pH (total scale)$", fontsize=16, color="red")
for label in ax2.get_xticklabels():
label.set_color("red")
ax1.invert_yaxis()
ax1.set_ylabel(r"$depth (m)$", fontsize=16);
In [ ]: