Examples of propagating uncertainties in mocsy

James Orr - 11 November 2018

LSCE/IPSL, CEA-CNRS-UVSQ, Gif-sur-Yvette, France

1. Setup

1.1 Download and build `mocsy` (Do this only once)

Put code into ./code/mocsy subdirectory and build mocsy.so for import into python

The build process might take a couple of minutes, depending our your computer



In [ ]:

    
%%bash
pwd
mkdir code
cd code
git clone https://github.com/jamesorr/mocsy.git
cd mocsy
make
pwd

1.2 Import standard python libraries



In [1]:

    
## Preliminaries
import numpy as np
import pandas as pd
import os, sys

1.3 Import mocsy (after specifying location of its shared object file)

The first uncommented line below is the subdirectory where the make command was executed (where mocsy.so is located)



In [2]:

    
# Comment out 1st line below, and Uncomment 2nd line below (if you have run the cell above under section 1.1)
mocsy_dir = "/homel/orr/Software/fortran/mocsy"
#mocsy_dir = "./code/mocsy"

sys.path.append (mocsy_dir)
import mocsy

2. Compute derived variables using mocsy's `vars` routine

For starters, begin with a simple example of how to call (in python) mocsy's most used routine, vars

Some basic documentation



In [3]:

    
print mocsy.mvars.vars.__doc__









    



ph,pco2,fco2,co2,hco3,co3,omegaa,omegac,betad,rhosw,p,tempis = vars(temp,sal,alk,dic,sil,phos,patm,depth,lat,optcon,optt,optp,[optb,optk1k2,optkf,optgas,opts,lon,verbose])

Wrapper for ``vars``.

Parameters
----------
temp : input rank-1 array('d') with bounds (n)
sal : input rank-1 array('d') with bounds (n)
alk : input rank-1 array('d') with bounds (n)
dic : input rank-1 array('d') with bounds (n)
sil : input rank-1 array('d') with bounds (n)
phos : input rank-1 array('d') with bounds (n)
patm : input rank-1 array('d') with bounds (n)
depth : input rank-1 array('d') with bounds (n)
lat : input rank-1 array('d') with bounds (n)
optcon : input string(len=6)
optt : input string(len=7)
optp : input string(len=2)

Other Parameters
----------------
optb : input string(len=3), optional
    Default: 'l10'
optk1k2 : input string(len=3), optional
    Default: 'l'
optkf : input string(len=2), optional
    Default: 'pf'
optgas : input string(len=7), optional
    Default: 'pinsitu'
opts : input string(len=4), optional
    Default: 'sprc'
lon : input rank-1 array('d') with bounds (n)
verbose : input int

Returns
-------
ph : rank-1 array('d') with bounds (n)
pco2 : rank-1 array('d') with bounds (n)
fco2 : rank-1 array('d') with bounds (n)
co2 : rank-1 array('d') with bounds (n)
hco3 : rank-1 array('d') with bounds (n)
co3 : rank-1 array('d') with bounds (n)
omegaa : rank-1 array('d') with bounds (n)
omegac : rank-1 array('d') with bounds (n)
betad : rank-1 array('d') with bounds (n)
rhosw : rank-1 array('d') with bounds (n)
p : rank-1 array('d') with bounds (n)
tempis : rank-1 array('d') with bounds (n)

Specify input variables and options



In [4]:

    
# Options:
optCON  = 'mol/kg' 
optT    = 'Tinsitu'
optP    = 'm'      
#optB    = 'l10'
optB    = 'u74'
optK1K2 = 'l'
optKf   = 'dg'

# Standard 6 input variables
temp   = 18.0  
sal    = 35.0  
alk    = 2300.e-6
dic    = 2000.e-6
#phos   = 2.0e-6   
#sil   = 60.0e-6 
phos   = 0.0e-6   
sil    = 0.0e-6 

# Other standard input (depth, atm pressure, latitude)
depth  = 0.
Patm   = 1.0        
lat    = 0.

Call the `vars` routine



In [5]:

    
ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis =    mocsy.mvars.vars(
temp, sal, alk, dic, sil, phos, Patm, depth, lat,
optCON, optT, optP, optb=optB, optk1k2=optK1K2, optkf=optKf    )

Print output



In [6]:

    
print ph, pco2, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis









    



[ 8.15244557] [ 298.14547549] [ 297.10842043] [  1.01871371e-05] [ 0.00177925] [ 0.00021056] [ 3.25176155] [ 5.02857045] [ 9.68872071] [ 1025.27206239] [ 0.] [ 18.]

3. Example use of `errors` routine



In [7]:

    
print mocsy.merrors.__doc__









    



eh,epco2,efco2,eco2,ehco3,eco3,eomegaa,eomegac = errors(temp,sal,alk,dic,sil,phos,patm,depth,lat,etemp,esal,ealk,edic,esil,ephos,optcon,optt,optp,[optb,optk1k2,optkf,optgas,opts,lon,r,epk,ebt])

Wrapper for ``errors``.

Parameters
----------
temp : input rank-1 array('d') with bounds (n)
sal : input rank-1 array('d') with bounds (n)
alk : input rank-1 array('d') with bounds (n)
dic : input rank-1 array('d') with bounds (n)
sil : input rank-1 array('d') with bounds (n)
phos : input rank-1 array('d') with bounds (n)
patm : input rank-1 array('d') with bounds (n)
depth : input rank-1 array('d') with bounds (n)
lat : input rank-1 array('d') with bounds (n)
etemp : input rank-1 array('d') with bounds (n)
esal : input rank-1 array('d') with bounds (n)
ealk : input rank-1 array('d') with bounds (n)
edic : input rank-1 array('d') with bounds (n)
esil : input rank-1 array('d') with bounds (n)
ephos : input rank-1 array('d') with bounds (n)
optcon : input string(len=6)
optt : input string(len=7)
optp : input string(len=2)

Other Parameters
----------------
optb : input string(len=3), optional
    Default: 'u74'
optk1k2 : input string(len=3), optional
    Default: 'l'
optkf : input string(len=2), optional
    Default: 'pf'
optgas : input string(len=7), optional
    Default: 'pinsitu'
opts : input string(len=7), optional
    Default: 'sprc'
lon : input rank-1 array('d') with bounds (n)
r : input float, optional
    Default: 0.0
epk : input rank-1 array('d') with bounds (7), optional
    Default: (0.002,0.0075,0.015,0.01,0.01,0.02,0.02)
ebt : input float, optional
    Default: 0.02

Returns
-------
eh : rank-1 array('d') with bounds (n)
epco2 : rank-1 array('d') with bounds (n)
efco2 : rank-1 array('d') with bounds (n)
eco2 : rank-1 array('d') with bounds (n)
ehco3 : rank-1 array('d') with bounds (n)
eco3 : rank-1 array('d') with bounds (n)
eomegaa : rank-1 array('d') with bounds (n)
eomegac : rank-1 array('d') with bounds (n)

Define errors



In [9]:

    
#etemp, esal = 0.01, 0.01
etemp, esal = 0.0, 0.0
ealk, edic = 2e-6, 2e-6
esil = 5e-6
ephos = 0.1e-6

Basic call to `errors` routine



In [10]:

    
[eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC] = \
mocsy.merrors(
temp, sal, alk, dic, sil, phos, Patm, depth, lat,
etemp, esal, ealk, edic, esil, ephos,
optCON, optT, optP, optb=optB, optk1k2=optK1K2, optkf=optKf, optgas='Pinsitu', ebt=0.01)
print eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC









    



[  1.88235588e-10] [ 9.42437971] [ 9.39158749] [  3.18579388e-07] [  4.31598697e-06] [  3.08021391e-06] [ 0.15712288] [ 0.24297706]



In [14]:

    
[eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC] = \
mocsy.merrors(
temp, sal, alk, dic, sil, phos, Patm, depth, lat,
etemp, esal, ealk, edic, esil, ephos,
optCON, optT, optP, optb=optB, optk1k2=optK1K2, optkf=optKf, optgas='Pinsitu')
print eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC









    



[  1.92787864e-10] [ 9.62209854] [ 9.58861836] [  3.25406441e-07] [  4.44296389e-06] [  3.27782684e-06] [ 0.15807355] [ 0.2444472]

Call `errors` specifying most all of the arguments, including the optional ones for just the errors routine (r, epk, and ebt).

In the cell below, the results would be identical if those 3 options were not present, because the values listed are the defaults in the errors routine:

r = 0.0 (no correlation between uncertainties in input pair (ealk and edic).
epk = a 7-member vector for the uncertainties in the equil. constants (K0, K1, K2, Kb, Kw, Ka, Kc) in pK units
ebt = the uncertainty in total boron, a relative fractional error [i.e., ebt = 0.02 is a 2% error, the default]

3.1 Errors, specifying standard r, epK, eBt



In [15]:

    
[eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC] = \
mocsy.merrors(
temp, sal, alk, dic, sil, phos, Patm, depth, lat,
etemp, esal, ealk, edic, esil, ephos,
optCON, optT, optP, optb=optB, optk1k2=optK1K2, optkf=optKf, 
r=0.0, epk=np.array([0.002,0.0075,0.015,0.01,0.01,0.02,0.02]), ebt=0.02)

print eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC









    



[  1.92787864e-10] [ 9.62209854] [ 9.58861836] [  3.25406441e-07] [  4.44296389e-06] [  3.27782684e-06] [ 0.15807355] [ 0.2444472]

3.2 Errors, assuming defaults for r, epK, eBt



In [16]:

    
[eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC] = \
mocsy.merrors(
temp, sal, alk, dic, sil, phos, Patm, depth, lat,
etemp, esal, ealk, edic, esil, ephos,
optCON, optT, optP, optb=optB, optk1k2=optK1K2, optkf=optKf)

print eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC









    



[  1.92787864e-10] [ 9.62209854] [ 9.58861836] [  3.25406441e-07] [  4.44296389e-06] [  3.27782684e-06] [ 0.15807355] [ 0.2444472]

3.3 Errors, specifying 0 for r, epK, eBt



In [17]:

    
[eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC] = \
mocsy.merrors(
temp, sal, alk, dic, sil, phos, Patm, depth, lat,
etemp, esal, ealk, edic, esil, ephos,
optCON, optT, optP, optb=optB, optk1k2=optK1K2, optkf=optKf, 
epk=np.array([0.000,0.000,0.00,0.00,0.00,0.00,0.00]), ebt=0.00)

print eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC









    



[  7.51501817e-11] [ 3.72766291] [ 3.71469247] [  1.27367921e-07] [  3.43481808e-06] [  1.87214484e-06] [ 0.02891188] [ 0.04470975]

3.3 Errors, specifying r=1.0



In [18]:

    
[eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC] = \
mocsy.merrors(
temp, sal, alk, dic, sil, phos, Patm, depth, lat,
etemp, esal, ealk, edic, esil, ephos,
optCON, optT, optP, optb=optB, optk1k2=optK1K2, optkf=optKf, r=1.0)

print eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC









    



[  1.77609844e-10] [ 8.88676459] [ 8.85584301] [  2.99999693e-07] [  3.40300633e-06] [  2.69160507e-06] [ 0.15541128] [ 0.24033022]