mocsy  2.0
Fortran 95 routines to model ocean carbonate system thermodynamics
 All Classes Namespaces Files Functions Variables Pages
Using mocsy in python

This page is for users that want to access the vast capabilities of python, e.g. reading & writing netCDF files easily.

The mocsy package, written in Fortran, is also designed to be called within a python script. Using mocsy from within python offers a simple and more flexible interface to handle input and output files and to integrate with other packages. For example, the python interface makes it easy to read and write netCDF files as well as to interface with R using rpy2.

Make it useable

The first step, before importing mocsy in python, is to make a shared object module using the utility f2py. This fortran-to-python lilbrary converter is part of the python's widely used NumPy library. Secondly, you will need a Fortran 95 compiler (e.g., GNU's gfortran). With both NumPy and gfortran installed, just go to the directory where you downloaded the mocsy routines, and type


If successful, output from that command will end in "Removing build directory ..." and the new file will have been produced. If unsuccesful, test to see that you have gfortran installed by typing

f2py -c --help-fcompiler

If –fcompiler=gnu95 does not show up just under the line "Fortram compilers found:", then try again to reinstall it. Otherwise, if you already have another Fortran 95 compiler installed, simply edit the makefile and change the name of the FC and F90 variables, accordingly. In addition, change the corresponding part of the last line –fcompiler=gnu95 to the name of your compiler, as listed under the "Fortran compilers found:" line from the above command.

Import it

After a successful build, the shared-object file will be found in the local directory. It is then imported into a python script like any other python library with

import mocsy

Documentation in python

Once mocsy is imported, you can see the basic routines that are available, with the standard python library documentation command:

print mocsy.__doc__

For details about the arguments of each routine, just add the name of the "module", e.g.,

print mocsy.mvars.__doc__

One could in addition add the name of the "routine", e.g.,

print mocsy.mvars.vars__doc__

but that is usually unnecessary. Typically the documentation is identical because there is but 1 routine per module, except for the gasx routine. More extensive documentation (on arguments, units, options, etc.) is available in this manual by clicking on the "Data File Types" tab and then the name of the module. The routines in python have the same names as in Fortran. The module names are the routine names preceded by a "m" (except for the gasx module, which contains several routines).

Use it

As an example, here is the call in python of the main routine vars

    pH, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, DENis, p, Tis = (
    mocsy.mvars (t, s, alk, dic, sio2, po4, Patm, depth, lat, optcon='mol/kg', optt='Tinsitu', optp='m', optb="u74", optk1k2='l', optkf="dg", optGAS="Ppot")

Input variables and options are on the right side of the equal sign, whereas output variables are on the left. All input and output arguments must be 1-D vectors (not multidimensional arrays); the conversion, both ways, is illustrated in the script as are other finer points for using mocsy in python. For instance, that script also uses the python library cdms to read input files (for DIC, ALk, T, S, \(\hbox{PO}_4^{-3}\), and \(\hbox{SiO}_2\)) and write the output files, all in netCDF format. Moreover, to save disk space for the user, input files are read remotely, using OpenDAP. This script is executed with the following command:


Other subroutines needed by mocsy.mvars are also available to the user. For example, to get the equilibrium constants that were used in the above call to mocsy.mvars, the command is

    Kh, K1, K2, Kb, Kw, Ks, Kf, Kspc, Kspa, K1p, K2p, K3p, Ksi, St, Ft, Bt = mocsy.mconstants (t, s, Patm, depth, lat, optt='Tinsitu', optp='m', optb="u74", optk1k2='l', optkf="dg", optGAS="Ppot")

Other routines are available to compute pressure, in-situ density, potential temperature, and in-situ temperature as well as to convert between pCO2, fCO2 and xCO2.

# Compute pressure (db) from depth (m) and latitude (degrees north) following Saunders (1981)
pfd = mocsy.mdepth2press(depth,lat)   

# Compute in situ density (kg/m3) from salinity (psu), in-situ temp (C), and hydrostatic pressure (db)
rhois = mocsy.mrhoinsitu(s, t, pfd)  

# Compute potential temperature (C) from salinity (psu), in situ temperature (C), and hydrostatic pressure (db)
Tpot = mocsy.mtpot(s, Tis, pfd, 0)     

# Compute in situ temperature (C) from salinity (psu), potential temperature (C), and hydrostatic pressure (db)
Tis = mocsy.mtis(s, Tpot, pfd, 0)

# Compute oceanic pCO2 (uatm) from arrays of fCO2 (uatm), in situ temp (C), atm pressure (atm), & hydrostatic pressure (db)
pCO2 = mocsy.mf2pCO2(fCO2, temp, Patm, p)

# Compute oceanic fCO2 (uatm) from arrays of pCO2 (uatm), in situ temp (C), atm pressure (atm), & hydrostatic pressure (db)
fCO2 = mocsy.mp2fCO2(pCO2, temp, Patm, p)

# Compute atmospheric xCO2 (ppm) from arrays of atmospheric pCO2 (uatm), surface ocean T (C), S (psu), & atm pressure (atm)
xCO2atm = mocsy.mgasx.pCO2atm2xCO2(pCO2atm, temp, salt, Patm)

# Compute atmospheric pCO2 (uatm) from arrays of atmospheric xCO2 (ppm), surface ocean T (C), S (psu), & atm pressure (atm)
pCO2atm = mocsy.mgasx.x2pCO2atm(xCO2, temp, salt, Patm)

# Compute vapor pressure of seawater (atm) following procedure from Weiss & Price (1980)
vpsw = mocsy.mgasx.vapress(temp, salt)

# Computes air-sea CO2 flux & surface-ocean carbonate system vars (see mocsy doc & OCMIP2 Abiotic HOWTO)
co2flux,co2ex,dpco2,ph,pco2,fco2,co2,hco3,co3,OmegaA,OmegaC,BetaD,rhoSW,p,tempis = \
mocsy.mgasx.flxco2(temp, sal, alk, dic, sil, phos, kw660, xco2, Patm, dz1,         \
optCON, optT, optP, optB, optK1K2, optKf, optGAS)


1) Simple scaler variables:

For DATA input: DIC and ALk in mol/kg, in situ temperature, pressure.

$ python
import mocsy
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')

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 \(K_1\) and \(K_2\):

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')

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').

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')

2) Simple arrays (numpy):

$ python
import numpy as np
import mocsy

#Make dummy array with 11 members
one = np.ones(11,dtype='float32')

sal   = 35*one
temp  = 2*one
depth = np.arange(0,11000,1000,dtype='float32') # units in 'm'

# 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)  

# 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

# 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

# Compute carbonate system variables, using 'model' options
pH,pco2,fco2,co2,hco3,co3,OmegaA,OmegaC,BetaD,DENis,p,Tis = \
mocsy.mvars (temp=temp, sal=sal, alk=alk, dic=dic, sil=sio2, phos=po4, patm=1, depth=depth, lat=lat, optcon='mol/m3', optt='Tpot', optp='m')

# Print out depth and carbonate ion concentration (converted to umol/kg)
print depth, co3 / (DENis * 1e-6)

3) Read input from a .csv (spreadsheet) file:

$ python
import numpy as np
import mocsy

# 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)
infilename = "DIC-Alk-P_vary_input.csv"
indata = np.loadtxt(infilename, delimiter=',', dtype='float32', skiprows=1)
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.zeros([n,], dtype='float32')

# Recover computed carbonate system variables
ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis = \
mocsy.mvars (temp, sal, alk,dic, sil, phos, depth, lat, optcon='mol/kg', optt='Tinsitu', optp='db', optb="u74", optk1k2='l', optkf="dg")

# Save results to another .csv file
results = np.dstack ((sal, temp, p, ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, dic, alk))
outfilename = "DIC-Alk-P-vary-from-mocsy.csv"
outfile = open(outfilename, "w")
outfile.write ("S, T, P, pH, pCO2, fCO2, CO2, HCO3, CO3, OmegaAragonite, OmegaCalcite, BetaD, DIC, ALK\n")
np.savetxt(outfile, results[0], delimiter=',', fmt="%1.6g")

4) Example with netCDF files for input and output:

MODEL example: See the python script to compute carbonate chemistry variables from CMIP5 model input.

DATA example: See the python script to compute carbonate chemistry variables from GLODAP DIC and Alk along with T, S, PO43- and SiO2 from the World Ocean Atlas (2009)