mocsy
2.0
Fortran 95 routines to model ocean carbonate system thermodynamics

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.
The first step, before importing mocsy in python, is to make a shared object module mocsy.so using the utility f2py. This fortrantopython 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
make mocsy.so
If successful, output from that command will end in "Removing build directory ..." and the new file mocsy.so will have been produced. If unsuccesful, test to see that you have gfortran installed by typing
f2py c helpfcompiler
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.
After a successful build, the sharedobject file mocsy.so will be found in the local directory. It is then imported into a python script like any other python library with
import mocsy
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).
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 1D vectors (not multidimensional arrays); the conversion, both ways, is illustrated in the script mocsy_GLODAP.py 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:
python mocsy_GLODAP.py
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, insitu density, potential temperature, and insitu 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), insitu 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 airsea CO2 flux & surfaceocean 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.e6, dic=2000.e6, 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 bestpractices 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*1028e6, dic=2000*1028e6, 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), insitu 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.0e3 po4 = po4 * 1.0e3 # Convert Alk and DIC units (umol/kg) to model units (mol/m3) dic = dic * rhois * 1e6 alk = alk * rhois * 1e6 # 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 * 1e6)
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 = "DICAlkP_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') lat.fill(45.) # 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 = "DICAlkPvaryfrommocsy.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") outfile.close()
4) Example with netCDF files for input and output:
MODEL example: See the python script mocsy_CMIP5.py to compute carbonate chemistry variables from CMIP5 model input.
DATA example: See the python script mocsy_GLODAP.py to compute carbonate chemistry variables from GLODAP DIC and Alk along with T, S, PO43 and SiO2 from the World Ocean Atlas (2009)