diff --git a/README.md b/README.md
new file mode 100644
index 0000000..279ae54
--- /dev/null
+++ b/README.md
@@ -0,0 +1,78 @@
+
+# pCO2 box
+> Box model simulating carbon solubility and air-sea exchange in the ocean.
+
+
+This package includes all code neccessary to run simulation and generate figues for the paper "Controls on Surface Water Carbonate Chemistry along North American Ocean Margins" by Cai et al.
+
+## Background and boundary conditions
+To study the relationships that control surface ocean carbonate parameters and air-sea CO2 fluxes, we use a box model to representing an idealized mixed layer. The box model is simulating the carbonate systems as a homogenous surface-ocean water mass located either in the North-West Atlantic (NWA, 40.1-42.5°N, 69.2-67.5°W), South Atlantic Bight (28.3°N–29.7°N, 77.6°W–79.2°W), or the California Current System (CCS, 40.8–41.6°N, 124.5–125°W). The model is driven by realistic temperatures and salinities from the Mercator 1/12° data-assimilated General Circulation Model with a daily resolution in time[^1]. 
+
+Biological production is assessed from weekly averages of daily changes in satellite derived Chlorophyll for the year 2015 using the daily MODIS aqua 4-km Level 3 product (NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group. Moderate-resolution Imaging Spectroradiometer (MODIS) Aqua Chlorophyll Data; 2018 Reprocessing. NASA OB.DAAC, Greenbelt, MD, USA. doi: data/10.5067/AQUA/MODIS/L3M/CHL/2018. Accessed on 10/02/2019) or for the year 2016 for the west coast box. We use extract daily time series for each grid cell that falls within the NWA, SAB and CCS model regions, identify all pairs of consecutive days with valid data and convert changes in Chl to a carbon flux by using a fixed  C:Chl ratio of 60 (Ref63). The resulting weekly averages are interpolated to daily time series. The resulting change in DIC is significantly smaller than other sources and sinks in the current study. 
+
+Winds are prescribed to 7 m/s in the summer and 10.5 m/s in the winter, to be consistent with representative wind data from the NOAA/Seawinds blended wind data set62.  Winter mixing in NWA is simulated by adding 0.1 mmol m-3 DIC daily from October through February (over 155 days). The value is based on a scaling analysis of vertical DIC gradients and diffusivity estimates in the region based on observational in 2015 in the NWA region and 2016 in the CCS. Literature values of physical vertical diffusion combined with vertical profiles of DIC from the ECOA 2015 cruse suggests that diffusive transports of DIC is negligible during summer. We use constant winds to minimize noise and to isolate the effect by changes in solubility on the carbonate system. Atmospheric pCO2 is set to 395 μatm. All carbonate equilibria and concentrations are calculated using the CO2SYS package[^29].
+
+## Box model initiation and iteration
+The model is initiated with the carbonate system defined by pCO2 in equilibrium with atmospheric values and alkalinity defined by initial salinity. The carbonate system is defined from DIC and TA and re-adjusted to changes in temperature and salinity at each time step using the CO2SYS package. pCO2_t(aq) is calculated using CO2SYS as well. ΔDICBio and ΔDICVertical are prescribed while ΔpCO2_Air-Sea is calculated at each time step using the relationship ΔpCO2_Air-Sea = 0.24 * kappa * K0 * (pCO2_t(aq) - (pCO2_(atm)). K0 is the CO2 gas solubility[^2] and kappa is defined as kappa = 0.251 * W2 * (Sc/660)-0.5 where W is wind speed in m/s and Sc is the Schmidt number[^65,^66].  Total alkalinity is calculated from prescribed salinity[^67].  The model is initiated to an equilibrium where pCO2_t0(aq) = pCO2_(atm), TA_t0 is based on salinity and temperature at t=0, and DIC is calculated from CO2SYS using pCO2(aq),  TAt0, temperature and salinity at t = 0. The model is spun up for 180 days using prescribed physical and biological conditions. The model iterated in time using the iteration
+
+```
+DICt+1 = DICt + ΔpCO2_Air-Sea + ΔDICBio + ΔDICVertical.
+```
+
+where each timestep is specifically conducted with the steps
+
+```
+1. Calculate the carbonate system defined by temperature, salinity, DIC, and TA at t_n.
+2. Apply changes to the DIC concentration at t_n+1 by biological processes, 
+vertical mixing and air-sea exchange (based on pCO2 values at t_n).
+3. Recalculate the carbonate system using temperature, salinity, and TA at t_n,  and DIC at t_n+1.
+4. Record pH, pCO2, and Ωarag at t_n+1.
+```
+
+The reason to use temperature, salinity, and TA at t_n  when performing step 3 is to make sure that all carbon parameters are calculated using the same physical conditions. We also calculate DICeq, a property that is based on the same TA, salinity, and temperature as the model but with pCO2(aq) relaxed to pCO2(atm) which can be interpreted as the air-sea exchange being infinitely fast. Here for simplicity, we simply assume pCO2(atm)  equals the dry CO2 mole fraction. 
+
+
+## Installing / Getting started
+
+Install the model with the following steps. The code and approaches has been tested on macOS and Linux and require a functioning python3.6 installed.
+
+```shell
+python3 -m venv /tmp/pco2box
+source /tmp/pco2box/bin/activate
+curl -O https://rsg.pml.ac.uk/shared_files/brj/pco2box.zip
+unzip pco2box.zip
+cd pco2box
+pip install .
+```
+
+Run the model with
+
+```python
+import pco2box
+
+pco2box.plot_perturbation(pert="ca") # California Current condtitions
+pco2box.plot_perturbation(pert="ne") # North West Atlantic conditions
+pco2box.plot_tmat_wind_mld()
+```
+
+This should create a folder called `pCO2box_figs` with the figures used in the paper.
+
+## Links
+
+Even though this information can be found inside the project on machine-readable
+format like in a .json file, it's good to include a summary of most useful
+links to humans using your project. You can include links like:
+
+- Project homepage: https://your.github.com/awesome-project/
+- Repository: https://github.com/your/awesome-project/
+- Issue tracker: https://github.com/your/awesome-project/issues
+  - In case of sensitive bugs like security vulnerabilities, please contact
+    my@email.com directly instead of using issue tracker. We value your effort
+    to improve the security and privacy of this project!
+- Related projects:
+  - Your other project: https://github.com/your/other-project/
+  - Someone else's project: https://github.com/someones/awesome-project/
+
+
+## Licensing
+The code in this project is licensed under MIT license.
diff --git a/figs/pertubation_timeseries_ca.pdf b/figs/pertubation_timeseries_ca.pdf
new file mode 100644
index 0000000..6a94715
Binary files /dev/null and b/figs/pertubation_timeseries_ca.pdf differ
diff --git a/figs/pertubation_timeseries_ne.pdf b/figs/pertubation_timeseries_ne.pdf
new file mode 100644
index 0000000..780e3d4
Binary files /dev/null and b/figs/pertubation_timeseries_ne.pdf differ
diff --git a/oceanbox/__init__.py b/oceanbox/__init__.py
new file mode 100644
index 0000000..14cffea
--- /dev/null
+++ b/oceanbox/__init__.py
@@ -0,0 +1,635 @@
+import time
+import copy
+
+import numpy as np
+import pylab as pl
+import pandas as pd
+
+import seawater as sw
+from .wind2pv import wind2pv, schmidt
+from . import stload
+from . import reuerflux
+import co2sys
+
+try:
+    import click
+    HAS_CLICK = True
+except ImportError:
+    HAS_CLICK = False
+    
+def nodim(a,b):  return (a-b) / ((a+b)/2)
+def nodimb(a,b): return (a-b) / b
+
+def clean(a):
+    a[np.isinf(a)]=np.nan
+    return a[~isnan(a)]
+
+
+class BoxModel(object):
+    """1D boxmodel to calculate air-sea fluxes of O2 and CO2"""
+    def __init__(self, svec=365):
+        self.zerovec = np.squeeze(np.zeros(svec))
+        self.tvec = np.arange(len(self.zerovec[:,...]))
+        self.reset()
+
+    def reset(self):
+        self.setup_physics()
+        self.setup_o2ar()
+        self.setup_carb()
+        
+    def setup_physics(self):
+        initpar = {"temp":28.0, "salt":33.0, "nwnd":10.0, "mldp":30.0,
+                   "ncpm":0.0,  "ncpc":0.0,  "bubl":0.01}
+        for key in initpar:
+            setattr(self, key, self.zerovec + initpar[key])
+
+    def setup_o2ar(self):
+        """Setup needed parameters to model O2 and Ar"""
+        for key in ['o2armsk', 'o2fl', 'arfl']:
+            setattr(self, key, self.zerovec.copy())
+        self.o2ct = sw.satO2(self.salt, self.temp)
+        self.arct = sw.satAr(self.salt, self.temp)
+        self.o2st = sw.satO2(self.salt, self.temp)
+        self.arst = sw.satAr(self.salt, self.temp)
+        #self.o2pv =  wind2pv(self.temp, self.nwnd)
+    
+    def setup_carb(self, pco2_atm=None, pco2=None, talk=None):
+        """Setup needed parameters to model CO2"""
+        self.pco2_atm = self.zerovec+395.0 if pco2_atm is None else pco2_atm
+        self.pco2     = self.zerovec+395.0 if pco2     is None else pco2.copy()
+        if talk is None:
+            self.talk = self.zerovec + co2sys.calc_TAlk(self.salt, self.temp)
+        else:
+            self.talk = talk
+        co = co2sys.CarbonateSystem(
+            self.salt, self.temp, TA=self.talk, pCO2=self.pco2)
+        self.pH   = self.zerovec + co.pH
+        self.DIC  = self.zerovec + co.TC
+        self.OmAr = self.zerovec + co.OmAr
+
+    def load():
+        return stload.load('tpz_nwnd_ncpm.npz')
+
+    @property
+    def fvent(self):
+        return self.o2pv / self.mldp
+
+    @property
+    def dens(self):
+        return sw.dens0(self.salt, self.temp)
+
+    @property
+    def o2pv(self):
+        return wind2pv(self.temp, self.nwnd)
+
+    
+    def calc_CO2_flux(self, tpos=None, w14=True):
+        """Calculate CO2 air-sea flux"""
+        Sc = schmidt(self.temp[tpos], "co2")
+        dpCO2 = self.pco2[tpos] - self.pco2_atm[tpos]
+        k0 = co2sys.constants.calculate_K0(self.salt[tpos], self.temp[tpos])
+        if w14:
+            kappa = 0.251 * self.nwnd[tpos]**2 * (Sc/660)**-0.5
+        else:
+            kappa = 0.24 * 0.0280 * self.nwnd[tpos]**3 * (Sc/660)**-0.5
+        co2fl = kappa * k0 * dpCO2 # (mmol m^-2 d^-1)
+        return co2fl
+        
+    def calc_wwfl(self, wtlen=60):   
+        end = self.pv.shape[0]
+        wpv = self.pv.copy()
+        weights = np.ones(self.pv[wtlen:,...].shape)
+        wghtsum = np.ones(self.pv[wtlen:,...].shape)
+        for t in np.arange(1,wtlen):
+            weights = (weights * (1-self.fvent[wtlen-t+1:end-t+1,...]))
+            wghtsum += weights
+            wpv[wtlen:,...] += self.pv[wtlen-t:end-t,...] * weights
+        wpv[wtlen:,...] = wpv[wtlen:,...] / wghtsum
+        wpv[:wtlen,...] = np.nan
+        self.wpv = wpv
+        self.wwfl = wpv * (self.o2ar) * self.o2st/1000 * self.dens
+
+    def step_co2(self,tpos):
+        """Calculate the carbonate system at tpos+1"""
+        co = co2sys.CarbonateSystem(self.salt[tpos], self.temp[tpos],
+                                 TA=self.talk[tpos], TC=self.DIC[tpos])
+        self.pco2[tpos] = co.pCO2
+        #print("%03i, %2.2f, %s" % (tpos, self.temp[tpos], co))
+
+        co2fl = self.calc_CO2_flux(tpos=tpos) * self.dens[tpos] / 1000
+        self.DIC[tpos+1] = (self.DIC[tpos]
+                                - co2fl / self.mldp[tpos] -
+                            self.ncpc[tpos])
+        co = co2sys.CarbonateSystem(self.salt[tpos],    self.temp[tpos],
+                                 TA=self.talk[tpos], TC=self.DIC[tpos+1])
+        self.pH[tpos+1]   = co.pH
+        self.pco2[tpos+1] = co.pCO2
+        self.OmAr[tpos+1] = co.OmAr
+        
+    def step_o2(self, tpos):
+        # Independent vars
+        self.o2st[tpos] = sw.satO2(self.salt[tpos], self.temp[tpos])
+        o2fl = self.o2pv[tpos] * (self.o2ct[tpos]-self.o2st[tpos])
+        self.o2ct[tpos+1] = (self.o2ct[tpos] -
+            (o2fl - self.ncpm[tpos])/self.mldp[tpos])
+
+    """
+    def step_o2ar(self, tpos):
+        # Independent vars
+        self.o2ct[self.o2armsk==0] = self.o2st[self.o2armsk==0] 
+        nrm = self.o2st[t] / self.arst[t]
+        self.o2fl[t] = self.pv[t] * (self.o2ct[t]-self.o2st[t])
+        self.arfl[t] = self.pv[t] * (self.arct[t]-self.arst[t])        
+        self.o2ct[t+1] = (self.o2ct[t] -
+                        (self.o2fl[t] - self.bubl[0]*20 - self.ncpm[t])/self.mldp[t])
+        self.arct[t+1] = self.arct[t] - (self.arfl[t] - self.bubl[t])/self.mldp[t]
+        self.o2ar = (self.o2ct/self.o2st)/(self.arct/self.arst)*1-1
+    """
+
+
+    """
+    def calc_nrm(self, ncpcut=3):
+        if ncpcut >= 0:
+            self.ncpm[self.ncpm<=ncpcut] = np.nan
+        self.ncpm[self.o2ar<=0] = np.nan
+        self.nn10[np.isnan(self.ncpm)] = np.nan
+        if ncpcut >= 0:
+            self.wwfl[self.ncpm<=ncpcut] = np.nan
+        self.wwfl[self.o2ar<=0] = np.nan
+
+        self.nrm10 = nodimb(self.wwfl,self.nn10)
+        self.nrm01 = nodimb(self.wwfl,self.ncpm)
+        self.nrm10[self.nrm10==-1] = np.nan
+        self.nrm01[self.nrm01==-1] = np.nan
+        self.nrm10[self.nrm10>10] = np.nan
+        self.nrm01[self.nrm01>10] = np.nan
+        self.nrm10[np.isinf(self.nrm10)] = np.nan
+        self.nrm01[np.isinf(self.nrm01)] = np.nan
+    """
+
+
+    def temp_drop(self):
+        self.setup_physics()
+        self.setup_o2ar()
+        self.setup_carb()
+        self.temp[100:] = self.temp[100:] - 3
+        self.run()
+        
+
+    def run(self, label=None):
+
+        #self.salt[100:] = self.salt[100:] - 2
+        if HAS_CLICK:
+            label= "Run model" if label is None else label
+            with click.progressbar(self.tvec[:-1], label=label) as bar:
+                for tpos in bar:
+                    self.step_co2(tpos)
+                    self.step_o2(tpos)
+        else:
+            for tpos in self.tvec[:-1]:
+                self.step_co2(tpos)
+                self.step_o2(tpos)
+
+    def to_csv(self, filename, extra_cols={}):
+        if type(extra_cols) is str:
+            extra_cols = {extra_cols:getattr(self, extra_cols)}
+        elif type(extra_cols) is list:
+            extra_cols = {key:getattr(self, key) for key in extra_cols}
+
+        df = pd.DataFrame(
+            {**{"salt":self.salt, "temp":self.temp, "wind":self.nwnd,
+                "pco2_atm":self.pco2_atm,
+                "pco2":self.pco2, "talk":self.talk, "pH":self.pH,
+                "DIC":self.DIC, "OmAr":self.OmAr}, **extra_cols})
+        df.to_csv(filename)
+
+
+    def plot_temps(self):
+        def run(dt):
+            self.reset()
+            self.temp[100:] = self.temp[100:] - dt
+            self.run()
+            pl.plot(self.DIC, label="%i $\degree$C drop" % dt, alpha=0.5)
+        pl.clf()
+        run(1)
+        run(2)
+        run(3)
+        pl.xlabel("Time (days)", fontsize=12)
+        pl.ylabel("DIC (mmol/m$^3$)", fontsize=12)
+        pl.tick_params(axis='both', which='major', labelsize=12)
+        pl.legend()
+        pl.title("Windspeed = 10 m/s")
+
+    def plot_winds(self):
+        def run(dt):
+            self.reset()
+            self.temp[100:] = self.temp[100:] - 3
+            self.nwnd[:]    = dt
+            self.run()
+            pl.plot(self.DIC, label="%i m/s wind" % dt, alpha=0.5)
+        pl.clf()
+        run(5)
+        run(10)
+        run(15)
+        run(20)
+        pl.xlabel("Time (days)", fontsize=12)
+        pl.ylabel("DIC (mmol/m$^3$)", fontsize=12)
+        pl.tick_params(axis='both', which='major', labelsize=12)
+        pl.legend()
+        pl.title("Temperature drop = 3 $\degree$C")
+
+
+    def plot_hurricane(self):
+        self.reset()
+        self.temp[100:103] = self.temp[100:103] - 3
+        self.nwnd[100] = 35
+        self.run()
+        pl.clf()
+        #pl.plot(self.DIC, alpha=0.5)
+        #pl.ylabel("DIC (mmol/m$^3$)", fontsize=12)
+        #pl.ylim(1900,1925)
+        #
+        #pl.plot(self.pco2, alpha=0.5)
+        #pl.ylabel("pCO$_2$ ($\mu$Atm)", fontsize=12)
+        #
+        pl.plot(self.OmAr, alpha=0.5)
+        pl.ylabel("$\Omega_{Ar}$", fontsize=12)
+        
+        pl.xlabel("Time (days)", fontsize=12)
+        pl.tick_params(axis='both', which='major', labelsize=12)
+        pl.title("Wind 35 m/s, Temperature drop = 3 $\degree$C")
+
+
+        
+    def plot(self, var1, var2):
+
+        ax1 = pl.gca()
+        ax2 = ax1.twinx()
+
+        p1, = ax2.plot(self.DIC,  'r-', label="DIC")
+        p2, = ax2.plot(self.pco2, 'g-', label="pCO$_2$")
+        ax1.set_xlabel("Time (days)")
+        ax1.set_ylabel("DIC")
+        ax2.set_ylabel("pCO$_2$")
+        """
+        lines = [p1, p2, p3]
+        host.legend(lines, [l.get_label() for l in lines])
+
+        for ax in [par1, par2]:
+            ax.set_frame_on(True)
+            ax.patch.set_visible(False)
+
+        plt.setp(ax.spines.values(), visible=False)
+        ax.spines["right"].set_visible(True)
+        """
+        ax1.yaxis.label.set_color(p1.get_color())
+        ax2.yaxis.label.set_color(p2.get_color())
+        #par1.spines["right"].set_edgecolor(p2.get_color())
+        #par2.spines["right"].set_edgecolor(p3.get_color())
+        ax1.tick_params(axis='y', colors=p1.get_color())
+        ax2.tick_params(axis='y', colors=p2.get_color())
+
+        
+
+#####################################
+
+    
+def run_model(st):
+    for t in np.arange(0,364):
+        # Independent vars
+        self.o2ct[self.o2armsk==0] = self.o2st[self.o2armsk==0] 
+        nrm = self.o2st[t]/self.arst[t]
+        self.o2fl[t] = self.pv[t] * (self.o2ct[t]-self.o2st[t])
+        self.arfl[t] = self.pv[t] * (self.arct[t]-self.arst[t])        
+        self.o2ct[t+1] = (self.o2ct[t] -
+            (self.o2fl[t] - self.bubl[0]*20 - self.ncpm[t])/self.mldp[t])
+        self.arct[t+1] = (self.arct[t] -
+            (self.arfl[t] - self.bubl[t])/self.mldp[t])
+    self.o2ar = (self.o2ct/self.o2st)/(self.arct/self.arst)*1-1
+    return st
+
+
+#####################################
+
+    
+def steady():    
+    st = setup()
+    st = setup_o2ar(st)
+    st = run_model(st)
+    plot1D(st)
+
+def pulse():
+    st = setup()
+    self.nwnd[100:110] = 15
+    self.ncpm[100:110] = 50
+    st = setup_o2ar(st)
+    st = run_model(st)
+    pl.figure(1)
+    plot1D(st)
+    pl.figure(2)
+    plotReuer(st)
+    return st
+
+def growth_season(tpp,im,tl=180,i1=None,i2=None,j1=None,w=False):
+    if w:
+        tp = tpp
+        st = setup(tp.nwnd[:,10:12,10:12].shape)
+        self.nwnd[...] = w
+        self.ncpm = np.zeros(self.nwnd.shape) *0.01
+    else:
+        tp = subset(tpp,im)
+        st = setup(tp.nwnd[:,:im,:].shape)
+        self.nwnd = tp.nwnd[:,...].copy()
+        self.ncpm = np.zeros(self.nwnd.shape) *0.01
+    xvec = np.arange(tl)
+    tsin = np.zeros(self.nwnd.shape[0]) *0.01
+    tsin[90:90+tl,...] = np.sin(xvec/((len(xvec)-1)/np.pi))*40
+    self.ncpm[:,...] = tsin[:,np.newaxis,np.newaxis]
+    self.ncpf = self.ncpm.copy()
+
+    st = setup_o2ar(st)
+    st = run_model(st)
+    st = calc_wwfl(st)
+    st = calc_nn10(st)
+    st = calc_nrm(st,ncpcut=-1)
+    return st
+
+def legtext(txt,nrm):
+    return ('%s  avg=%.2f   std=%.2f' %
+            (txt,np.mean(nrm[~np.isnan(nrm)]),
+             np.std(nrm[~np.isnan(nrm)])) )
+
+def season_hist(tp,im):
+    pl.close(1)
+    F = pl.figure(1,(11,8))
+    slvec =  [30,90,150,210]
+    pl.clf()
+    legvec = []
+    for n,s in enumerate(slvec):
+        st = growth_season(tp,im=im,tl=s)
+        pl.hist(self.nrm10[~np.isnan(self.nrm10)],np.arange(-2,2,0.02),
+                histtype='step')
+        legvec.append(legtext('%i days' % s,self.nrm10))
+    pl.legend(legvec)
+    pl.title('Misfit for different lenghts of Growth Seasons')
+    pl.ylabel('Number of data points')
+    pl.xlabel('Misfit ( (bioflux-NCP)/ncp )')
+    pl.savefig('figs/toymodel/seasonhiself.png')
+
+def season_wind(tp,im,i1=10,j1=10):
+    slvec =  [30,90,150,210]
+    pl.close(1)
+    F = pl.figure(1,(11,8))
+    w = 10
+    for n,s in enumerate(slvec):
+        print(n)
+        st = growth_season(tp,im,tl=s,w=w)
+        pl.subplot(4,1,n+1)
+        pl.plot(self.ncpf[:,0,0])
+        pl.plot(self.nn10[:,0,0])
+        pl.plot(self.wwfl[:,0,0])
+        pl.xticks(np.arange(0,400,50),[])
+        pl.yticks([10,30])
+        pl.title('Growth period = %03i days' % s)
+        pl.legend ( ('NCP instant','NCP 10d avg', 'Bioflux'))
+        pl.xlim(0,365)
+        pl.ylabel('mmol/m2')
+    pl.xticks(np.arange(0,400,50),np.arange(0,400,50))
+    pl.xlabel('Days')
+    pl.suptitle('Constant wind (%i m/s)' % w)
+    pl.savefig('figs/toymodel/seasonwind_%03i.png' % w)
+    
+    slvec = np.arange(30,250,5)
+    wvec = [5,10,15]
+
+    pl.close(2)
+    F = pl.figure(2,(11,8))
+    legvec = []
+    for w in wvec:
+        msvec = []
+        for s in slvec:
+            st = growth_season(tp,im,tl=s,w=w)
+            msvec.append(np.mean(abs(self.nrm10[~np.isnan(self.nrm10)])))
+        pl.plot(slvec,msvec,'.')
+        legvec.append('%02i m/s' % w)
+    pl.legend(legvec)
+    pl.ylim(0,1)
+    pl.ylabel('Average misfit during growth season, sign removed.')
+    pl.xlabel('Growing season (days)')
+    pl.suptitle('Constant wind speed vs length of growth season')
+    pl.savefig('figs/toymodel/seasonwind_all.png')
+  
+
+def season_plot(tp,i1=10,j1=10):
+    slvec =  [30,90,150,210]
+    pl.figure(2)
+    pl.clf()
+    for n,s in enumerate(slvec):
+        print(n)
+        st = growth_season(tp,tl=s,i1=i1,j1=j1)
+        pl.figure(2)
+        pl.subplot(4,1,n+1)
+        pl.plot(self.nn10)
+        pl.plot(self.wwfl)
+        pl.xticks(np.arange(0,400,50),[])
+        pl.yticks([10,30])
+        pl.title('Growth period = %03i days' % s)
+    pl.xticks(np.arange(0,400,50),np.arange(0,400,50))
+    pl.xlim(0,365)
+    
+    slvec = np.arange(30,250,5)
+    msvec = []
+    for s in slvec:
+        st = growth_season(tp,tl=s,i1=i1,j1=j1)
+        msvec.append(np.mean(abs(self.msft)))
+
+
+    pl.legend ( ('NCP 10d avg', 'Bioflux'))
+    pl.figure(3)
+    #pl.clf()
+    pl.plot(slvec,msvec,'.')
+    #pl.xlim(20,200)
+    pl.ylim(0,1)
+
+    pl.ylabel('Average misfit during growing season')
+    pl.xlabel('Growing sweason (days)')
+    pl.figure(2)
+    pl.savefig('figs/seasonplot_%03i_%03i_2.png' % (i1,j1) )
+    pl.figure(3)
+    pl.savefig('figs/seasonplot_%03i_%03i_3.png' % (i1,j1) )
+
+def season_plot_2D(tpp,im=30):
+    class tp: pass
+    tpvars = ['llat', 'llon', 'mldp','ncpm','nwnd', 'temp']    
+    for v in tpvars:
+        tp.__dict__[v] = tpp.__dict__[v][...,:im,:].copy()
+    slvec =  [30,90,150,210]
+    pl.figure(2)
+    pl.clf()
+    for n,s in enumerate(slvec):
+        st = growth_season(tp,tl=s)
+        pl.figure(n+1)
+        ncphist(st,self.msft)
+
+def realwind(tp):
+    st = setup(im,jm)
+
+    self.nwnd = tp.nwnd[:,40,10]
+    #self.ncpm = tp.ncpm[:,40,10]
+    #self.temp = tp.temp[:,40,10]
+
+    self.o2ct[0] = 200
+    self.arct[0] = 11
+    self.o2st = sw.satO2(self.salt, self.temp)
+    self.arst = sw.satAr(self.salt, self.temp)
+    self.pv   = wind2pv(self.temp,self.nwnd)
+    
+    st = run_model(st)
+    pl.figure(1)
+    plot1D(st)
+    pl.figure(2)
+    plotReuer(st)
+    return st
+
+def realwind2D(tpp,im,realvars=['nwnd','ncpm']):
+    tp = subset(tpp,im)
+    st = setup(tp.nwnd[:,:im,:].shape)
+    if 'nwnd' in realvars:
+        self.nwnd = tp.nwnd[:,:im,:].copy()
+    if 'ncpm' in realvars:
+        self.ncpm = tp.ncpm[:,:im,:].copy()
+        self.ncpm[np.isnan(self.ncpm)] = 0
+    if 'mldp' in realvars:
+        self.mldp = tp.mldp[:,:im,:].copy()
+    if 'temp' in realvars:
+        self.temp = tp.temp[:,:im,:].copy()
+    if 'o2armsk' in realvars:
+        self.o2armsk = tp.o2ar[:,:im,:].copy()
+        self.o2armsk[self.o2armsk>0] = 1
+        self.o2armsk[self.o2armsk<0] = 0
+    st = setup_o2ar(st)
+    st = run_model(st)
+    st = calc_wwfl(st)
+    st = calc_nn10(st)
+    st = calc_nrm(st)
+    return st
+
+def hist_realwind2D(tp,im,realvars=['nwnd','ncpm']):
+    st = realwind2D(tp,im,realvars=realvars)
+
+    pl.clf()
+    nvec = self.nrm10[~np.isnan(self.nrm10)]
+    h = pl.hist(nvec,np.arange(-2,2,0.02),histtype='step')
+    nvec = self.nrm01[~np.isnan(self.nrm01)]
+    h = pl.hist(nvec,np.arange(-2,2,0.02),histtype='step')
+
+    if 1==0:
+        tp = calc_nrm(tp)
+        self.nrm01[np.isnan(tp.nrm01)] = np.nan
+        nvec = tp.nrm10[~np.isnan(tp.nrm10)]
+        h = pl.hist(nvec,np.arange(-2,2,0.02),histtype='step')
+        nvec = tp.nrm01[~np.isnan(tp.nrm01)]
+        h = pl.hist(nvec,np.arange(-2,2,0.02),histtype='step')
+    pl.legend( ('NCP 10-day mean','NCP instant','full10d','fullinst') )
+    return st
+
+
+def ncphist(st):
+
+    self.nrm10[np.isinf(self.nrm10)]=np.nan
+    nvec= [0,10,20,30,40,50]
+    leglist = []
+    def hist(nrm):
+        return pl.hist(nrm[~np.isnan(nrm)],np.arange(-2,2,0.02),
+                       histtype='step')
+    hist(self.nrm10)
+    legliself.append(legtext('All data',self.nrm10))
+    for n in nvec:
+        nrm = self.nrm10.copy()
+        nrm[self.nn10>n+10] = np.nan
+        nrm[self.nn10<n]    = np.nan 
+        hist(nrm)
+        legliself.append(legtext('%02i-%02i' % (n,n+10),nrm))
+    pl.legend (leglist)
+    pl.title('Misfit for different NCP values')
+    pl.ylabel('Number of data points')
+    pl.xlabel('Misfit ( (bioflux-NCP)/ncp )')
+    pl.savefig('figs/toymodel/ncphiself.png')
+
+
+def plot1D(st):
+    pl.clf()
+    pl.subplot(4,1,1)
+    pl.plot(self.o2ct/self.o2st*100 - 100)
+    pl.plot(self.arct/self.arst*100 - 100)
+    pl.ylabel('Percent')
+    lg = pl.legend(('O2','Ar'))
+    lg.get_frame().set_fill(False)
+    pl.title('Saturation')
+    pl.ylim(-10,10)
+    pl.xticks([])
+
+    pl.subplot(4,1,2)
+    pl.plot(self.o2fl)
+    pl.plot(self.arfl*20)
+    lg = pl.legend(('O2','Ar'))
+    lg.get_frame().set_fill(False)
+    pl.title('Flux')
+    pl.ylabel('mmol / d')
+    pl.ylim(0,100)
+    pl.xticks([])
+    
+    pl.subplot(4,1,3)
+    pl.plot(self.o2ar*100)
+    pl.title('O2Ar Saturation')
+    pl.ylim(-10,10)
+    pl.ylabel('Percent')
+    pl.xticks([])
+    
+    pl.subplot(4,1,4)
+    pl.plot(self.o2ar * self.pv*self.o2st)
+    pl.plot(self.ncpm)
+    pl.title('O2bio Flux')
+    pl.ylabel('mmol / d')
+    pl.ylim(0,100)
+
+def plotReuer(st, wtlen=60):
+    st = calc_wwfl(st)
+    pl.clf()
+    pl.subplot(4,1,1)
+    pl.plot(self.nwnd)
+    pl.plot(self.pv)
+    lg = pl.legend(('Wind','PV'))
+    lg.get_frame().set_fill(False)
+    pl.title('Wind (m/s) & PV (m/d)')
+    pl.ylim(0,20)
+    pl.xticks([])
+    pl.xlim(0,400)
+    
+    pl.subplot(4,1,2)
+    pl.plot(self.fvent)
+    pl.title('F-vent')
+    pl.ylabel('m^2 / d')
+    pl.ylim(0,1)
+    pl.xlim(0,400)
+    
+    pl.subplot(4,1,3)
+    pl.plot(self.wpv)
+    pl.title('Weighted PV')
+    pl.ylabel('m / d')
+    pl.ylim(0,20)
+    #pl.ylabel('Percent')
+    pl.xlim(0,400)
+    
+    pl.subplot(4,1,4)
+    pl.plot(self.wwfl)
+    pl.plot(self.ncpm)
+
+    pl.title('O2bio Flux')
+    pl.ylabel('mmol / d')
+    pl.ylim(0,100)
+    pl.xlim(0,400)
+
+
+#md = BoxModel()
+#md.run()             
diff --git a/oceanbox/reuerflux.py b/oceanbox/reuerflux.py
new file mode 100644
index 0000000..3de3db4
--- /dev/null
+++ b/oceanbox/reuerflux.py
@@ -0,0 +1,89 @@
+import numpy as np
+
+import seawater as sw
+from . import wind2pv
+
+
+def pv2wpv(temp, mld, wind, wtlen=60, salt=False, dens=None):
+    """Calculate bio-o2flux using weighted wind """
+    dens = sw.dens0(salt, temp) if dens is None else dens
+    pv    = wind2pv(temp,wind)
+    fvent = pv / mld
+    fvent[np.isnan(fvent)] = 0
+    pv[np.isnan(pv)] = 0
+    if wtlen == 0:
+        wpv = pv
+    else:
+        end=pv.shape[0]
+        wpv = pv.copy()
+        weights = np.ones(pv[wtlen:,...].shape)
+        wghtsum = np.ones(pv[wtlen:,...].shape)
+        for t in np.arange(1,wtlen):
+            weights = (weights * (1-fvent[wtlen-t+1:end-t+1,...]))
+            wghtsum += weights
+            wpv[wtlen:,...] += pv[wtlen-t:end-t,...] * weights
+        wpv[wtlen:,...] = wpv[wtlen:,...] / wghtsum
+    wpv[:wtlen,...] = np.nan
+    return wpv
+
+def wpvweight(temp, mld, wind, wtlen=60, salt=False, dens=False):
+    """Calculate bio-o2flux using weighted wind """
+    if type(dens) == bool : dens = sw.dens0(salt, temp)
+    pv    = wind2pv(temp,wind)
+    fvent = pv / mld
+    fvent[np.isnan(fvent)] = 0
+    pv[np.isnan(pv)] = 0
+    if wtlen == 0:
+        wpv = pv
+    else:
+        end=pv.shape[0]
+        wpv = field.copy()
+        weights = np.ones(pv[wtlen:,...].shape)
+        wghtsum = np.ones(pv[wtlen:,...].shape)
+        for t in np.arange(1,wtlen):
+            weights = (weights * (1-fvent[wtlen-t+1:end-t+1,...]))
+            wghtsum += weights
+            wpv[wtlen:,...] += pv[wtlen-t:end-t,...] * weights
+        wpv[wtlen:,...] = wpv[wtlen:,...] / wghtsum
+    wpv[:wtlen,...] = np.nan
+    return wpv[:,0,...]
+
+def test_range(patm):
+    if patm.min() < 500 or patm.min() > 1500: 
+        print("Atmospheric pressure on in range")
+        raise
+
+def pv2fl(o2ar,temp,mld,wind,wtlen=60,
+          salt=False, o2st=False, dens=False, patm=False):
+    """Calculate bio-o2flux using weighted wind """
+    if type(o2st) == bool : o2st = sw.satO2(salt, temp)
+    if type(patm) == np.ndarray:
+        print("Using atm press in reuerflux")
+        test_range(patm)
+        o2st=o2st * (patm/1013.25)
+    if type(dens) == bool : dens = sw.dens0(salt, temp)
+    pv    = wind2pv(temp,wind)
+    fvent = pv / mld
+    fvent[np.isnan(fvent)] = 0
+    pv[np.isnan(pv)] = 0
+    if wtlen == 0:
+        wpv = pv
+    else:
+        end=pv.shape[0]
+        wpv = pv.copy()
+        weights = np.ones(pv[wtlen:,...].shape)
+        wghtsum = np.ones(pv[wtlen:,...].shape)
+        for t in np.arange(1,wtlen):
+            weights = (weights * (1-fvent[wtlen-t+1:end-t+1,...]))
+            wghtsum += weights
+            wpv[wtlen:,...] += pv[wtlen-t:end-t,...] * weights
+        wpv[wtlen:,...] = wpv[wtlen:,...] / wghtsum
+    wpv[:wtlen,...] = np.nan
+    return (wpv[:,np.newaxis] * (o2ar/100) * o2st/1000 * dens)[:,0,...]
+
+
+def bottle2o2ar(o2ar_air, temp, salt):
+
+    air = 20.946/0.9340
+    o2ar_sat = (sw.satO2(salt, temp) / sw.satAr(salt, temp)/air-1) * 1000
+    return o2ar_air / o2ar_sat - 1
diff --git a/oceanbox/stload.py b/oceanbox/stload.py
new file mode 100644
index 0000000..a1e0daf
--- /dev/null
+++ b/oceanbox/stload.py
@@ -0,0 +1,75 @@
+
+import os 
+import numpy as np
+
+def save(S,name):
+
+
+    savedict ={}
+
+    for at in S.__dict__:
+        if type(S.__dict__[at]) == str:
+            savedict[at] = S.__dict__[at]
+        elif type(S.__dict__[at]) == list:
+            savedict[at] = S.__dict__[at]
+        elif type(S.__dict__[at]) == int:
+            savedict[at] = S.__dict__[at]
+        elif type(S.__dict__[at]) == float:
+            savedict[at] = S.__dict__[at]
+        elif type(S.__dict__[at]) == np.ndarray:
+            savedict[at] = S.__dict__[at]
+        elif type(S.__dict__[at]) == np.ma.core.MaskedArray:
+            savedict[at] = S.__dict__[at].data
+            savedict[at + '_mask'] = S.__dict__[at].mask
+    
+    np.savez(name,**savedict)
+
+def load(name,S=-999, t1=None, t2=None, k1=None, k2=None,
+         i1=None,i2=None,j1=None,j2=None):
+
+    if S == -999:
+        class S:pass
+    tl = np.load(name)
+
+    pref = os.path.split(name)[1][:-4]
+    try:
+        if not '__tm__' in tl.files:
+            f = np.load(pref + "_ijk.npz")
+            __tm__ = f['tm'].item()
+            __km__ = f['km'].item()
+            __im__ = f['im'].item()
+            __jm__ = f['jm'].item()   
+            
+        if not '__shape__' in tl.files:
+            __shape__ = np.load(pref + "_shape.npz")["shp"].item()
+        load = sliceload
+    except IOError:
+        def load(var): return tl[var]
+
+    def sliceload(var):
+        if __shape__[var] == (__tm__,__km__,__im__,__jm__):
+            return tl[f][t1:t2,k1:k2,i1:i2,j1:j2]
+        elif __shape__[var] == (__tm__,__im__,__jm__):
+            return tl[f][t1:t2,i1:i2,j1:j2]
+        elif __shape__[var] == (__km__,__im__,__jm__):
+            print(var + " kij")
+            return tl[f][k1:k2,i1:i2,j1:j2]
+        elif __shape__[var] == (__im__,__jm__):
+            return tl[f][i1:i2,j1:j2]
+        elif  len(__shape__[var]) == 1:
+            if __shape__[var] == (__jm__):
+                return tl[f][j1:j2]
+            elif __shape__[var] == (__im__):
+                return tl[f][i1:i2]
+        return tl[f]
+        
+    for f in tl.files:
+        if f + '_mask' in tl.files:
+            field = load(f)
+            mask =  load(f + '_mask')
+            setattr(S, f, np.ma.masked_array(field,mask))
+        if '_mask' in f:
+            pass
+        else:
+           setattr(S, f, load(f))
+    return S
diff --git a/oceanbox/wind2pv.py b/oceanbox/wind2pv.py
new file mode 100644
index 0000000..29189df
--- /dev/null
+++ b/oceanbox/wind2pv.py
@@ -0,0 +1,39 @@
+import warnings
+import numpy as np
+
+def wind2pv(temp, wind, model='swee', gas="o2"):
+    temp = np.array(temp)
+    wind = np.array(wind)
+    if temp.ndim > wind.ndim:
+        temp = temp[:,0,...]
+    assert not wind[~np.isnan(wind)].min() < 0, "Negative wind speeds."
+    if any([s in model.lower() for s in ['wann','w92']]):
+        return _wind2pv_wann1992(temp, wind, gas)
+    elif  any([s in model.lower() for s in ['swee','sw07']]):
+        return _wind2pv_sweeney2007(temp, wind, gas)
+
+def _wind2pv_wann1992(temp, wind, gas="o2"):
+    sc = schmidt(temp, gas)
+    return (0.31 * wind**2 * (sc/660)**(-0.5)) * (24./100) 
+
+def _wind2pv_sweeney2007(temp, wind, gas="o2"):
+    sc = schmidt(temp, gas)
+    return (0.27 * wind**2 * (sc/660)**(-0.5)) * (24./100) 
+
+#def _wind2pv_
+
+
+def schmidt(temp, gas="o2"):
+    """Calculate the Schmidt number for different gases
+    
+    Formulation from Saramiento, Chemical Ocenaography, p 85. 
+    Based on Wilke and Chang (1955) via Wanninkhof (1992).
+    """
+    if gas =="o2":
+        Sc = {"A":1953.4, "B":128.00, "C":3.9918, "D":0.050091}
+    elif gas == "ar":
+        Sc = {"A":1909.1, "B":125.09, "C":3.9012, "D":0.048953}
+    elif gas == "co2":
+        Sc = {"A":2073.1, "B":125.62, "C":3.6276, "D":0.043219}
+    return Sc["A"] - Sc["B"]*temp + Sc["C"]*temp**2 - Sc["D"]*temp**3
+
diff --git a/pco2box/.DS_Store b/pco2box/.DS_Store
new file mode 100644
index 0000000..e3b4ef2
Binary files /dev/null and b/pco2box/.DS_Store differ
diff --git a/pco2box/__init__.py b/pco2box/__init__.py
new file mode 100644
index 0000000..a112690
--- /dev/null
+++ b/pco2box/__init__.py
@@ -0,0 +1,10 @@
+"""Model to study the solubility and gasexchange in a generalized ocean box """
+import os
+
+from .timeseries import plot_perturbation
+from .recovery import plot_tmat_wind_mld
+
+
+FIGDIR = os.path.abspath("./pCO2box_figs")
+if not os.path.exists(FIGDIR):
+    os.makedirs(FIGDIR)
diff --git a/pco2box/figures.py b/pco2box/figures.py
new file mode 100644
index 0000000..680b39b
--- /dev/null
+++ b/pco2box/figures.py
@@ -0,0 +1,85 @@
+
+import numpy as np
+import pylab as pl
+import seaborn as sns
+#import cartopy.crs as ccrs
+#crs = ccrs.PlateCarree()
+from matplotlib.dates import MonthLocator, WeekdayLocator, DateFormatter
+
+jdvec = np.arange(735599, 736481)[:730]
+months = MonthLocator(range(1, 13), bymonthday=1, interval=3)
+monthsFmt = DateFormatter("%b")
+
+#def east_map():
+#    pl.clf()
+#    fig = pl.figure()
+#    ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal(
+#        central_longitude=-70, central_latitude=45))
+#    ax.set_extent([-81, -58, 23, 47], crs=ccrs.PlateCarree())
+#    ax.coastlines('10m',zorder=3, color="0.5")
+#
+#    latvec = np.load("indata/merc_uscoast_latlon.npz")["lat"][[300,400,460]]
+#    lonvec = np.load("indata/merc_uscoast_latlon.npz")["lon"][[710,770,840]]
+#    pl.scatter(lonvec, latvec, transform=crs)
+#    transform = crs._as_mpl_transform(pl.gca())
+#    for txt,lon,lat in zip(["SAB","MAB", "NWA"], lonvec, latvec):
+#        pl.annotate(txt, xy=(lon,lat), xycoords=transform)
+#    pl.savefig("boxmodel_eastwest_east_map.png", bbox_inches="tight")
+
+def fldplot(md, fldname, ax=None, **kwargs):
+    """Plot list from model"""
+    if ax is None:
+        ax = pl.gca()
+    for args in {"lw":1, "alpha":0.5}.items():
+        kwargs[args[0]] = kwargs.get(args[0], args[1])
+        
+    p = ax.plot(jdvec,getattr(md, fldname), **kwargs)
+    ax.xaxis_date()
+    ax.xaxis.set_major_formatter(monthsFmt) 
+    ax.xaxis.set_major_locator(months) 
+    ax.grid(True)
+    ax.set_xlim("2015-01-01","2017-01-01")
+    return p[0].get_color()
+
+def plot_timeseries(md, **kwargs):
+    tsp=3
+    sns.set_style("whitegrid")
+    clrs =sns.color_palette()
+    if not "axes" in kwargs:
+        fig,axes = pl.subplots(tsp,1, sharex=True, num=1,figsize=(6,6))
+    else:
+        axes = kwargs.pop("axes")
+    pl.subplots_adjust(hspace=0.03)
+
+    def line(ax, fldname, ylabel, ylim=(None,None),axpos="left", yticks=None,
+                 **kwargs):
+        if axpos == "right":
+            kwargs["c"] = clrs[1]
+            ax = pl.twinx(ax=ax)
+            ax.yaxis.label.set_color(clrs[1])
+            ax.tick_params(axis='y', colors=clrs[1])
+        else:
+            kwargs["c"] = clrs[0]
+            ax.yaxis.label.set_color(clrs[0])
+            ax.tick_params(axis='y', colors=clrs[0])
+        fldplot(md, fldname, ax=ax, **kwargs)
+        if fldname =="DIC":
+            kwargs["label"] = None
+            fldplot(md, "DICeq", ax, ls=":", **kwargs)
+        ax.set_ylabel(ylabel)
+        ax.set_ylim(*ylim)
+        if yticks is not None:
+            ax.set_yticks(yticks)
+
+    args = [axes[0], "temp", "Temperature $\degree$C", [0,30]]
+    line(*args, yticks=[7.5,15,22.5], **kwargs)
+    args = [axes[0], "salt", "Salinity ", [32,36]]
+    line(*args, axpos="right", yticks=[33,34,35], **kwargs)
+    args = [axes[1], "pco2", "pCO$_2$ $\mu{}atm$", [300,1100]]
+    line(*args, yticks=[500,700,900], **kwargs)
+    args = [axes[1], "pH", "pH", [7.5,8.5], ]
+    line(*args, axpos="right", yticks=[7.75, 8.0, 8.25], **kwargs)
+    args = [axes[2], "DIC",  "DIC mmol m$^{-3}$", [1900,2300]]
+    line(*args, yticks=[2000,2100,2200], **kwargs) 
+    args = [axes[2], "OmAr",  "$\Omega$", [0.6, 3.0]]   
+    line(*args, axpos="right", yticks=[1.2, 1.8, 2.4],**kwargs)
diff --git a/pco2box/indata/merc_uscoast_latlon.npz b/pco2box/indata/merc_uscoast_latlon.npz
new file mode 100644
index 0000000..2e1528e
Binary files /dev/null and b/pco2box/indata/merc_uscoast_latlon.npz differ
diff --git a/pco2box/indata/mercator_tseries.h5 b/pco2box/indata/mercator_tseries.h5
new file mode 100644
index 0000000..9b82ef1
Binary files /dev/null and b/pco2box/indata/mercator_tseries.h5 differ
diff --git a/pco2box/indata/mercator_us_coasts_tseries.npz b/pco2box/indata/mercator_us_coasts_tseries.npz
new file mode 100644
index 0000000..d712081
Binary files /dev/null and b/pco2box/indata/mercator_us_coasts_tseries.npz differ
diff --git a/pco2box/model.py b/pco2box/model.py
new file mode 100644
index 0000000..65dd1d4
--- /dev/null
+++ b/pco2box/model.py
@@ -0,0 +1,261 @@
+
+import os
+
+import numpy as np
+import pylab as pl
+import pandas as pd
+from scipy.stats import linregress
+ 
+import oceanbox
+import co2sys
+
+from .figures import plot_timeseries
+
+DATADIR = os.path.join(os.path.dirname(__file__), "indata")
+
+def calc_talk(salt, pref):
+    """Calculate Total Alkalinity
+    
+    Return Total Alkalinity based on parametrizations for different parts
+    of the US East Coast. Based on personal communications with Wei-Jun Cai
+    and [REF]
+    
+    Parameters
+    ----------
+    salt : array_like
+        Salinity used to calculate TALk
+    pref : string
+        Part of the coast [se=South Atlantic Bight, 
+                           me=Mid Atlantic Bight, 
+                           ne=North Atlantic outside Gulf of Maine]
+
+    Returns
+    -------
+        TAlk : array_like
+            Calculated Total Alkalinity
+    """
+    if pref == "se":
+        talk = 48.7 * salt + 608.8
+    if pref == "me":
+        talk = 46.6 * salt + 670.6
+        #md.salt = md.salt - 1
+    if pref == "ne":
+        talk = 39.1 * salt + 932.7
+    return talk
+
+def setup_model(pref="ne", wind=7, mldp=30, wintermixing=True, winterwind=True):
+    """Setup Box model
+    
+    Initiate the boxmodel class with data from the indata/mercator_tseries.h5 
+    file. Wind and Mixed Layer Depths (MLDs) can be prescribed to fixed values 
+    and winter conditions can be turned on or off.
+
+    Parameters
+    ----------
+    pref : string, optional
+        Part of the coast [se=South Atlantic Bight, 
+                           me=Mid Atlantic Bight, 
+                           ne=North Atlantic outside Gulf of Maine]
+                           ca=California Coastal Current]
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    wintermixing : bool, optional
+        Simulate wintermixing by removing 0.1 mM DIC daily during winter
+    winterwind : bool, optional
+        Simulate winter condtions by increasing winds 1.5 times during winter
+    """
+    df = pd.read_hdf(os.path.join(DATADIR, "mercator_tseries.h5"))
+    md = oceanbox.BoxModel(svec=730)
+    md.nwnd = np.zeros((365)) + wind
+    if winterwind:
+        md.nwnd[:60]  = wind*1.5
+        md.nwnd[270:] = wind*1.5
+    md.nwnd = np.hstack((md.nwnd, md.nwnd))  
+    md.mldp[:] = mldp
+
+    def loop_years(vec):
+        slope = np.linspace(0,(vec[0] - vec[364]), 365)
+        vec = vec + slope
+        return np.hstack((vec, vec))
+
+    md.temp = loop_years(df[pref + "temp"][:365])
+    md.salt = loop_years(df[pref + "salt"][:365])
+    try:
+        md.ncpc = loop_years(df[pref + "ncpm"][:365])
+    except KeyError:
+        md.ncpc = loop_years(np.zeros(365))
+    if wintermixing:
+        md.ncpc[270:] = md.ncpc[270:] - 0.1
+        md.ncpc[:60]  = md.ncpc[:60]  - 0.1
+    return md
+
+def run_model(pref, reg, temp=None, salt=None, wind=7, mldp=30,
+        deepmld1={"pert":0, "tpos":165},
+        deepmld2={"pert":0, "tpos":553}, 
+        uppwell1={"pert":0, "tpos":553, "pdays":10},
+        uppwell2={"pert":0, "tpos":553, "pdays":10}):
+    """Run the boxmodel in pertubation mode
+
+    Setup and run the box model with the ability to apply an pertubation.
+
+    pref : string
+        Part of the coast [se=South Atlantic Bight, 
+                           me=Mid Atlantic Bight, 
+                           ne=North Atlantic outside Gulf of Maine]
+                           ca=California Coastal Current]
+    reg : string
+        Region name to be used in title and file name
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    temp : float, optional
+        Prescribed fixed temperature
+    salt : float, optional
+        Prescribed fixed salinity
+    deepmld1 : dict, optional
+        Dict describing a mixed layer deepening pertubation
+    deepmld2 : dict, optional
+        Dict describing a mixed layer deepening pertubation
+    uppwell1 : dict, optional
+        Dict describing a upwelling pertubation
+    uppwell2 : dict, optional
+        Dict describing a upwelling pertubation
+    """
+    md = setup_model(pref=pref, mldp=mldp, wind=wind)
+    if temp is not None:
+        md.temp[:] = temp
+    if salt is not None:
+        md.salt[:] = salt
+
+    def mix_deepmld(deepmld):
+        md.ncpc[deepmld["tpos"]] = md.ncpc[deepmld["tpos"]] - deepmld["pert"]
+        md.mldp[deepmld["tpos"]:deepmld["tpos"]+10] = mldp*2
+        ncpslice = slice(deepmld["tpos"]+5, deepmld["tpos"]+15)
+
+    def mix_upwell(uppwell):
+        tpos = uppwell["tpos"]
+        pdays = uppwell["pdays"]
+        md.ncpc[tpos:tpos+pdays+1] -= uppwell["pert"]/pdays
+
+    if deepmld1["pert"] != 0:
+        mix_deepmld(deepmld1)
+    if deepmld2["pert"] != 0:
+        mix_deepmld(deepmld2)
+    if uppwell1["pert"] != 0:
+        mix_upwell(uppwell1)
+    if uppwell2["pert"] != 0:
+        mix_upwell(uppwell2)
+
+    talk = calc_talk(md.salt, "ne")
+    co = co2sys.CarbonateSystem(md.salt, md.temp, TA=talk, pCO2=talk*0+395)
+    md.DICeq = co.TC            
+    md.setup_carb()
+    md.run(label="Run %s" % reg)
+    return md
+
+def perturbation(wind=7, mldp=30, temp=12, salt=33, mpos=165, pdays=40):
+    """Run the model twice, once with perturbation once without
+    
+    Parameters
+    ----------
+    pref : string, optional
+        Part of the coast [se=South Atlantic Bight, 
+                           me=Mid Atlantic Bight, 
+                           ne=North Atlantic outside Gulf of Maine]
+                           ca=California Coastal Current]
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    temp : float, optional
+        Prescribed fixed temperature
+    salt : float, optional
+        Prescribed fixed salinity
+    mpos : int, optional
+        Day of year of pertubation
+    pdays : int, optional
+        Length of upwelling pertubation
+
+    Returns
+    -------
+    mds : model class instance
+        Background case
+    mpd : model class instance
+        Pertubation case
+    """
+    kwargs = dict(pref="ne", reg="NWA",
+                  temp=temp, salt=salt, wind=wind, mldp=mldp,
+                  deepmld1={"pert":0, "tpos":mpos})
+    mds = run_model(**kwargs)
+    kwargs["deepmld1"] = {"pert":17*3, "tpos":mpos}
+    mdp = run_model(**kwargs)
+    return mds,mdp
+
+def plot_perturbation(wind=7, mldp=30, pref="ca"):
+    """Run the model in pertubation mode and plot the results"""
+    pl.clf()
+    fig,axes = pl.subplots(3, 1, sharex=True, num=1,figsize=(6,6))
+    model_kws = dict(pref=pref, reg="pert",
+                temp=None, salt=None, wind=wind, mldp=mldp,
+                deepmld1={"pert":0, "tpos":165},
+                deepmld2={"pert":0, "tpos":553},
+                uppwell1={"pert":0, "tpos":165, "pdays":5},
+                uppwell2={"pert":0, "tpos":553, "pdays":5})
+    md = run_model(**model_kws)
+    plot_timeseries(md, axes=axes, alpha=0.5) 
+    if pref == "ca":
+        preftxt = "CCS"
+        model_kws["uppwell1"]["pert"] = 82.5
+        model_kws["uppwell2"]["pert"] = 165
+    else:
+        preftxt = "NWA"
+        model_kws["deepmld1"]["pert"] = 17
+        model_kws["deepmld2"]["pert"] = 34
+    md = run_model(**model_kws)
+    plot_timeseries(md, axes=axes, alpha=1)
+    pl.suptitle(
+        f"Perturbations, temp and salt for {preftxt}, wind:{wind}m/s, mld:{mldp}m")
+    pl.savefig(f"figs/pertubation_timeseries_{pref}.pdf")
+
+def to_csv(wind=7, mldp=30, pref="ca"):
+    """Save model results to csv
+    
+    Parameters
+    ----------
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    pref : string, optional
+        Part of the coast [ne=North Atlantic outside Gulf of Maine]
+                           ca=California Coastal Current]    
+    """
+    ptext = "CCS" if pref == "ca" else "NWA"
+    if pref == "ca":
+        preftxt = "CCS"
+        model_kws = dict(pref=pref, reg="pert",
+                     temp=None, salt=None, wind=wind, mldp=mldp,
+                     deepmld1={"pert":0, "tpos":165},
+                     deepmld2={"pert":0, "tpos":553},
+                     uppwell1={"pert":82.5, "tpos":165, "pdays":5},
+                     uppwell2={"pert":165,  "tpos":553, "pdays":5})
+    else:
+        preftxt = "NWA"
+        model_kws = dict(pref=pref, reg="pert",
+                     temp=None, salt=None, wind=wind, mldp=mldp,
+                     deepmld1={"pert":17, "tpos":165},
+                     deepmld2={"pert":34, "tpos":553},
+                     uppwell1={"pert":0,  "tpos":165, "pdays":5},
+                     uppwell2={"pert":0,  "tpos":553, "pdays":5})
+
+    md = run_model(**model_kws)
+    md.to_csv(f"boxmodel_baseline_{ptext}__wind_{wind:02}__mldp_{mldp:03}.csv")
+    for key in ["deepmld1", "deepmld1", "uppwell1", "uppwell2"]:
+        model_kws[key] = {"pert":0, "tpos":165}
+    md = run_model(**model_kws)
+    md.to_csv(f"boxmodel_pertubation_{ptext}__wind_{wind:02}__mldp_{mldp:03}.csv")
+
+ 
diff --git a/pco2box/model.py~ b/pco2box/model.py~
new file mode 100644
index 0000000..fbc573b
--- /dev/null
+++ b/pco2box/model.py~
@@ -0,0 +1,257 @@
+
+import numpy as np
+import pylab as pl
+import pandas as pd
+from scipy.stats import linregress
+ 
+import oceanbox
+import co2sys
+
+from figures import plot_timeseries
+
+def calc_talk(salt, pref):
+    """Calculate Total Alkalinity
+    
+    Return Total Alkalinity based on parametrizations for different parts
+    of the US East Coast. Based on personal communications with Wei-Jun Cai
+    and [REF]
+    
+    Parameters
+    ----------
+    salt : array_like
+        Salinity used to calculate TALk
+    pref : string
+        Part of the coast [se=South Atlantic Bight, 
+                           me=Mid Atlantic Bight, 
+                           ne=North Atlantic outside Gulf of Maine]
+
+    Returns
+    -------
+        TAlk : array_like
+            Calculated Total Alkalinity
+    """
+    if pref == "se":
+        talk = 48.7 * salt + 608.8
+    if pref == "me":
+        talk = 46.6 * salt + 670.6
+        #md.salt = md.salt - 1
+    if pref == "ne":
+        talk = 39.1 * salt + 932.7
+    return talk
+
+def setup_model(pref="ne", wind=7, mldp=30, wintermixing=True, winterwind=True):
+    """Setup Box model
+    
+    Initiate the boxmodel class with data from the indata/mercator_tseries.h5 
+    file. Wind and Mixed Layer Depths (MLDs) can be prescribed to fixed values 
+    and winter conditions can be turned on or off.
+
+    Parameters
+    ----------
+    pref : string, optional
+        Part of the coast [se=South Atlantic Bight, 
+                           me=Mid Atlantic Bight, 
+                           ne=North Atlantic outside Gulf of Maine]
+                           ca=California Coastal Current]
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    wintermixing : bool, optional
+        Simulate wintermixing by removing 0.1 mM DIC daily during winter
+    winterwind : bool, optional
+        Simulate winter condtions by increasing winds 1.5 times during winter
+    """
+    df = pd.read_hdf("indata/mercator_tseries.h5")
+    md = oceanbox.BoxModel(svec=730)
+    md.nwnd = np.zeros((365)) + wind
+    if winterwind:
+        md.nwnd[:60]  = wind*1.5
+        md.nwnd[270:] = wind*1.5
+    md.nwnd = np.hstack((md.nwnd, md.nwnd))  
+    md.mldp[:] = mldp
+
+    def loop_years(vec):
+        slope = np.linspace(0,(vec[0] - vec[364]), 365)
+        vec = vec + slope
+        return np.hstack((vec, vec))
+
+    md.temp = loop_years(df[pref + "temp"][:365])
+    md.salt = loop_years(df[pref + "salt"][:365])
+    try:
+        md.ncpc = loop_years(df[pref + "ncpm"][:365])
+    except KeyError:
+        md.ncpc = loop_years(np.zeros(365))
+    if wintermixing:
+        md.ncpc[270:] = md.ncpc[270:] - 0.1
+        md.ncpc[:60]  = md.ncpc[:60]  - 0.1
+    return md
+
+def run_model(pref, reg, temp=None, salt=None, wind=7, mldp=30,
+        deepmld1={"pert":0, "tpos":165},
+        deepmld2={"pert":0, "tpos":553}, 
+        uppwell1={"pert":0, "tpos":553, "pdays":10},
+        uppwell2={"pert":0, "tpos":553, "pdays":10}):
+    """Run the boxmodel in pertubation mode
+
+    Setup and run the box model with the ability to apply an pertubation.
+
+    pref : string
+        Part of the coast [se=South Atlantic Bight, 
+                           me=Mid Atlantic Bight, 
+                           ne=North Atlantic outside Gulf of Maine]
+                           ca=California Coastal Current]
+    reg : string
+        Region name to be used in title and file name
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    temp : float, optional
+        Prescribed fixed temperature
+    salt : float, optional
+        Prescribed fixed salinity
+    deepmld1 : dict, optional
+        Dict describing a mixed layer deepening pertubation
+    deepmld2 : dict, optional
+        Dict describing a mixed layer deepening pertubation
+    uppwell1 : dict, optional
+        Dict describing a upwelling pertubation
+    uppwell2 : dict, optional
+        Dict describing a upwelling pertubation
+    """
+    md = setup_model(pref=pref, mldp=mldp, wind=wind)
+    if temp is not None:
+        md.temp[:] = temp
+    if salt is not None:
+        md.salt[:] = salt
+
+    def mix_deepmld(deepmld):
+        md.ncpc[deepmld["tpos"]] = md.ncpc[deepmld["tpos"]] - deepmld["pert"]
+        md.mldp[deepmld["tpos"]:deepmld["tpos"]+10] = mldp*2
+        ncpslice = slice(deepmld["tpos"]+5, deepmld["tpos"]+15)
+
+    def mix_upwell(uppwell):
+        tpos = uppwell["tpos"]
+        pdays = uppwell["pdays"]
+        md.ncpc[tpos:tpos+pdays+1] -= uppwell["pert"]/pdays
+
+    if deepmld1["pert"] != 0:
+        mix_deepmld(deepmld1)
+    if deepmld2["pert"] != 0:
+        mix_deepmld(deepmld2)
+    if uppwell1["pert"] != 0:
+        mix_upwell(uppwell1)
+    if uppwell2["pert"] != 0:
+        mix_upwell(uppwell2)
+
+    talk = calc_talk(md.salt, "ne")
+    co = co2sys.CarbonateSystem(md.salt, md.temp, TA=talk, pCO2=talk*0+395)
+    md.DICeq = co.TC            
+    md.setup_carb()
+    md.run(label="Run %s" % reg)
+    return md
+
+def perturbation(wind=7, mldp=30, temp=12, salt=33, mpos=165, pdays=40):
+    """Run the model twice, once with perturbation once without
+    
+    Parameters
+    ----------
+    pref : string, optional
+        Part of the coast [se=South Atlantic Bight, 
+                           me=Mid Atlantic Bight, 
+                           ne=North Atlantic outside Gulf of Maine]
+                           ca=California Coastal Current]
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    temp : float, optional
+        Prescribed fixed temperature
+    salt : float, optional
+        Prescribed fixed salinity
+    mpos : int, optional
+        Day of year of pertubation
+    pdays : int, optional
+        Length of upwelling pertubation
+
+    Returns
+    -------
+    mds : model class instance
+        Background case
+    mpd : model class instance
+        Pertubation case
+    """
+    kwargs = dict(pref="ne", reg="NWA",
+                  temp=temp, salt=salt, wind=wind, mldp=mldp,
+                  deepmld1={"pert":0, "tpos":mpos})
+    mds = run_model(**kwargs)
+    kwargs["deepmld1"] = {"pert":17*3, "tpos":mpos}
+    mdp = run_model(**kwargs)
+    return mds,mdp
+
+def plot_perturbation(wind=7, mldp=30, pref="ca"):
+    """Run the model in pertubation mode and plot the results"""
+    pl.clf()
+    fig,axes = pl.subplots(3, 1, sharex=True, num=1,figsize=(6,6))
+    model_kws = dict(pref=pref, reg="pert",
+                temp=None, salt=None, wind=wind, mldp=mldp,
+                deepmld1={"pert":0, "tpos":165},
+                deepmld2={"pert":0, "tpos":553},
+                uppwell1={"pert":0, "tpos":165, "pdays":5},
+                uppwell2={"pert":0, "tpos":553, "pdays":5})
+    md = run_model(**model_kws)
+    plot_timeseries(md, axes=axes, alpha=0.5) 
+    if pref == "ca":
+        preftxt = "CCS"
+        model_kws["uppwell1"]["pert"] = 82.5
+        model_kws["uppwell2"]["pert"] = 165
+    else:
+        preftxt = "NWA"
+        model_kws["deepmld1"]["pert"] = 17
+        model_kws["deepmld2"]["pert"] = 34
+    md = run_model(**model_kws)
+    plot_timeseries(md, axes=axes, alpha=1)
+    pl.suptitle(
+        f"Perturbations, temp and salt for {preftxt}, wind:{wind}m/s, mld:{mldp}m")
+    pl.savefig(f"figs/pertubation_timeseries_{pref}.pdf")
+
+def to_csv(wind=7, mldp=30, pref="ca"):
+    """Save model results to csv
+    
+    Parameters
+    ----------
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    pref : string, optional
+        Part of the coast [ne=North Atlantic outside Gulf of Maine]
+                           ca=California Coastal Current]    
+    """
+    ptext = "CCS" if pref == "ca" else "NWA"
+    if pref == "ca":
+        preftxt = "CCS"
+        model_kws = dict(pref=pref, reg="pert",
+                     temp=None, salt=None, wind=wind, mldp=mldp,
+                     deepmld1={"pert":0, "tpos":165},
+                     deepmld2={"pert":0, "tpos":553},
+                     uppwell1={"pert":82.5, "tpos":165, "pdays":5},
+                     uppwell2={"pert":165,  "tpos":553, "pdays":5})
+    else:
+        preftxt = "NWA"
+        model_kws = dict(pref=pref, reg="pert",
+                     temp=None, salt=None, wind=wind, mldp=mldp,
+                     deepmld1={"pert":17, "tpos":165},
+                     deepmld2={"pert":34, "tpos":553},
+                     uppwell1={"pert":0,  "tpos":165, "pdays":5},
+                     uppwell2={"pert":0,  "tpos":553, "pdays":5})
+
+    md = run_model(**model_kws)
+    md.to_csv(f"boxmodel_baseline_{ptext}__wind_{wind:02}__mldp_{mldp:03}.csv")
+    for key in ["deepmld1", "deepmld1", "uppwell1", "uppwell2"]:
+        model_kws[key] = {"pert":0, "tpos":165}
+    md = run_model(**model_kws)
+    md.to_csv(f"boxmodel_pertubation_{ptext}__wind_{wind:02}__mldp_{mldp:03}.csv")
+
+ 
\ No newline at end of file
diff --git a/pco2box/recovery.py b/pco2box/recovery.py
new file mode 100644
index 0000000..e2ca5f3
--- /dev/null
+++ b/pco2box/recovery.py
@@ -0,0 +1,101 @@
+import os
+
+import numpy as np
+import pylab as pl
+ 
+from .model import perturbation
+from .figures import plot_timeseries
+
+FIGDIR = os.path.abspath("./pCO2box_figs")
+
+
+def recovery_time(wind=7, mldp=30, temp=12, salt=33, pfunc=perturbation):
+    """Calculate recovery time for a specific mld and windspeed
+    
+    Run the model with and without a pertubation calculate the time it
+    takes for the difference to be less than 1/e of the initial difference.
+
+    Parameters
+    ----------
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    temp : float, optional
+        Prescribed fixed temperature
+    salt : float, optional
+        Prescribed fixed salinity
+    pfunc : function, optional
+        Function to calculate pertubation and baseline
+
+    Returns
+    -------
+    time : int
+        Time in days to return to 1/e difference between pertubation
+        and baseline timeseries.
+    """
+    mds,mdp = pfunc(wind=wind, mldp=mldp, temp=temp, salt=salt)
+    pvec = (mdp.DIC-mds.DIC)[mpos+1:]/(mdp.DIC-mds.DIC)[mpos+1] 
+    pvec[pvec<1/np.exp(1)] = 0 
+    return np.nonzero(pvec==0)[0].min()
+
+def perturbation_recovery_matrix(parname="DIC"):
+    """Compute array with recovery timescales for different winds and mld's
+    
+    Returns an array with recovery timescales based on different 
+    Mixed Layer Depths (MLDs) and windspeeds defined by the matrices
+
+        mmat,wmat = np.mgrid[10:100:5,1:20:2]
+
+    The boxmodel is run with fixed a prescribed temperature of temp 12°C and a
+    prescribed salinity of 33 PSU. The pertubation is set to occur after
+    165 day. This values has no effect on the results since T and S are fixed.
+
+    use the function 'plot_tmat_wind_mld' to visualize the results.
+
+    Parameters
+    ----------
+    parname : string, optional
+        Name of carbon species to base the recovery (default is DIC)
+
+    Returns
+    -------
+    recovery_time_matrix : ndarray
+        matrix with wind as xaxis and MLD as y-axis
+    """
+    mpos = 165 # Time of pertubation. Has no effect since T and S are fixed
+    mmat,wmat = np.mgrid[10:100:5,1:20:2]
+    tmat = mmat * np.nan
+    for m,mldp in enumerate(mmat[:,0]):
+        for n,wind in enumerate(wmat[0,:]):
+            mds,mdp = perturbation(wind=wind, mldp=mldp, mpos=mpos)
+            tvecs = getattr(mds, parname)
+            tvecp = getattr(mdp, parname)
+            pvec = (tvecp-tvecs)[mpos+1:]/(tvecp-tvecs)[mpos+1]             
+            pvec[pvec<1/np.exp(1)] = 0 
+            try:
+              tmat[m,n] = np.nonzero(pvec==0)[0].min()
+            except:
+                print(mldp,wind)
+                pass
+    return tmat
+
+def plot_tmat_wind_mld(tmat=None, parname="DIC", color="k", text=""):
+    """Plot tmat array"""
+    mmat,wmat = np.mgrid[10:100:5,1:20:2]
+    if tmat is None:
+        tmat = perturbation_recovery_matrix(parname=parname)
+    CS = pl.contour(wmat, mmat, tmat,
+                    [10,20,30,40,50,60,80,100,150,200],
+                    colors=color)
+    pl.gca().clabel(CS, inline=1, fontsize=10, fmt="%i")
+    CS.collections[0].set_label(text)
+    pl.xlim(3,18)
+    pl.xlabel("Wind speed (m/s)")
+    pl.ylabel("Mixed Layer Depth (m)")
+    pl.title("Time to recover from 51 mmol DIC perturbation (days)")
+    pl.savefig(os.path.join(FIGDIR, f"recovery_matrix.pdf"))
+
+    return CS
+
+#In [87]: plot([0.01, 0.3, 0.4, 0.5, 1, 1.5,2, 3, 4, 5], [30,28,27,26,22,18,14, 12, 10, 9], ".-")                      
diff --git a/pco2box/recovery.py~ b/pco2box/recovery.py~
new file mode 100644
index 0000000..eddd160
--- /dev/null
+++ b/pco2box/recovery.py~
@@ -0,0 +1,96 @@
+import numpy as np
+import pylab as pl
+ 
+from model import perturbation
+from figures import plot_timeseries
+
+
+
+def recovery_time(wind=7, mldp=30, temp=12, salt=33, pfunc=perturbation):
+    """Calculate recovery time for a specific mld and windspeed
+    
+    Run the model with and without a pertubation calculate the time it
+    takes for the difference to be less than 1/e of the initial difference.
+
+    Parameters
+    ----------
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    temp : float, optional
+        Prescribed fixed temperature
+    salt : float, optional
+        Prescribed fixed salinity
+    pfunc : function, optional
+        Function to calculate pertubation and baseline
+
+    Returns
+    -------
+    time : int
+        Time in days to return to 1/e difference between pertubation
+        and baseline timeseries.
+    """
+    mds,mdp = pfunc(wind=wind, mldp=mldp, temp=temp, salt=salt)
+    pvec = (mdp.DIC-mds.DIC)[mpos+1:]/(mdp.DIC-mds.DIC)[mpos+1] 
+    pvec[pvec<1/np.exp(1)] = 0 
+    return np.nonzero(pvec==0)[0].min()
+
+def perturbation_recovery_matrix(parname="DIC"):
+    """Compute array with recovery timescales for different winds and mld's
+    
+    Returns an array with recovery timescales based on different 
+    Mixed Layer Depths (MLDs) and windspeeds defined by the matrices
+
+        mmat,wmat = np.mgrid[10:100:5,1:20:2]
+
+    The boxmodel is run with fixed a prescribed temperature of temp 12°C and a
+    prescribed salinity of 33 PSU. The pertubation is set to occur after
+    165 day. This values has no effect on the results since T and S are fixed.
+
+    use the function 'plot_tmat_wind_mld' to visualize the results.
+
+    Parameters
+    ----------
+    parname : string, optional
+        Name of carbon species to base the recovery (default is DIC)
+
+    Returns
+    -------
+    recovery_time_matrix : ndarray
+        matrix with wind as xaxis and MLD as y-axis
+    """
+    mpos = 165 # Time of pertubation. Has no effect since T and S are fixed
+    mmat,wmat = np.mgrid[10:100:5,1:20:2]
+    tmat = mmat * np.nan
+    for m,mldp in enumerate(mmat[:,0]):
+        for n,wind in enumerate(wmat[0,:]):
+            mds,mdp = perturbation(wind=wind, mldp=mldp, mpos=mpos)
+            tvecs = getattr(mds, parname)
+            tvecp = getattr(mdp, parname)
+            pvec = (tvecp-tvecs)[mpos+1:]/(tvecp-tvecs)[mpos+1]             
+            pvec[pvec<1/np.exp(1)] = 0 
+            try:
+              tmat[m,n] = np.nonzero(pvec==0)[0].min()
+            except:
+                print(mldp,wind)
+                pass
+    return tmat
+
+def plot_tmat_wind_mld(tmat=None, parname="DIC", color="k", text=""):
+    """Plot tmat array"""
+    mmat,wmat = np.mgrid[10:100:5,1:20:2]
+    if tmat is None:
+        tmat = perturbation_recovery_matrix(parname=parname)
+    CS = pl.contour(wmat, mmat, tmat,
+                    [10,20,30,40,50,60,80,100,150,200],
+                    colors=color)
+    pl.gca().clabel(CS, inline=1, fontsize=10, fmt="%i")
+    CS.collections[0].set_label(text)
+    pl.xlim(3,18)
+    pl.xlabel("Wind speed (m/s)")
+    pl.ylabel("Mixed Layer Depth (m)")
+    pl.title("Time to recover from 51 mmol DIC perturbation (days)")
+    return CS
+
+#In [87]: plot([0.01, 0.3, 0.4, 0.5, 1, 1.5,2, 3, 4, 5], [30,28,27,26,22,18,14, 12, 10, 9], ".-")                      
diff --git a/pco2box/timeseries.py b/pco2box/timeseries.py
new file mode 100644
index 0000000..7a0e167
--- /dev/null
+++ b/pco2box/timeseries.py
@@ -0,0 +1,74 @@
+import os
+
+import numpy as np
+import pylab as pl
+ 
+from .model import run_model, perturbation
+from .figures import plot_timeseries
+
+FIGDIR = os.path.abspath("./pCO2box_figs")
+
+def plot_perturbation(wind=7, mldp=30, pref="ca"):
+    """Run the model in pertubation mode and plot the results"""
+    pl.clf()
+    fig,axes = pl.subplots(3, 1, sharex=True, num=1,figsize=(6,6))
+    model_kws = dict(pref=pref, reg="pert",
+                temp=None, salt=None, wind=wind, mldp=mldp,
+                deepmld1={"pert":0, "tpos":165},
+                deepmld2={"pert":0, "tpos":553},
+                uppwell1={"pert":0, "tpos":165, "pdays":5},
+                uppwell2={"pert":0, "tpos":553, "pdays":5})
+    md = run_model(**model_kws)
+    plot_timeseries(md, axes=axes, alpha=0.5) 
+    if pref == "ca":
+        preftxt = "CCS"
+        model_kws["uppwell1"]["pert"] = 82.5
+        model_kws["uppwell2"]["pert"] = 165
+    else:
+        preftxt = "NWA"
+        model_kws["deepmld1"]["pert"] = 17
+        model_kws["deepmld2"]["pert"] = 34
+    md = run_model(**model_kws)
+    plot_timeseries(md, axes=axes, alpha=1)
+    pl.suptitle(
+        f"Perturbations, temp and salt for {preftxt}, wind:{wind}m/s, mld:{mldp}m")
+    pl.savefig(os.path.join(FIGDIR, f"pertubation_timeseries_{pref}.pdf"))
+
+def to_csv(wind=7, mldp=30, pref="ca"):
+    """Save model results to csv
+    
+    Parameters
+    ----------
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    pref : string, optional
+        Part of the coast [ne=North Atlantic outside Gulf of Maine]
+                           ca=California Coastal Current]    
+    """
+    ptext = "CCS" if pref == "ca" else "NWA"
+    if pref == "ca":
+        preftxt = "CCS"
+        model_kws = dict(pref=pref, reg="pert",
+                     temp=None, salt=None, wind=wind, mldp=mldp,
+                     deepmld1={"pert":0, "tpos":165},
+                     deepmld2={"pert":0, "tpos":553},
+                     uppwell1={"pert":82.5, "tpos":165, "pdays":5},
+                     uppwell2={"pert":165,  "tpos":553, "pdays":5})
+    else:
+        preftxt = "NWA"
+        model_kws = dict(pref=pref, reg="pert",
+                     temp=None, salt=None, wind=wind, mldp=mldp,
+                     deepmld1={"pert":17, "tpos":165},
+                     deepmld2={"pert":34, "tpos":553},
+                     uppwell1={"pert":0,  "tpos":165, "pdays":5},
+                     uppwell2={"pert":0,  "tpos":553, "pdays":5})
+
+    md = run_model(**model_kws)
+    md.to_csv(f"boxmodel_baseline_{ptext}__wind_{wind:02}__mldp_{mldp:03}.csv")
+    for key in ["deepmld1", "deepmld1", "uppwell1", "uppwell2"]:
+        model_kws[key] = {"pert":0, "tpos":165}
+    md = run_model(**model_kws)
+    md.to_csv(f"boxmodel_pertubation_{ptext}__wind_{wind:02}__mldp_{mldp:03}.csv")
+    
diff --git a/pco2box/timeseries.py~ b/pco2box/timeseries.py~
new file mode 100644
index 0000000..d6d9aa8
--- /dev/null
+++ b/pco2box/timeseries.py~
@@ -0,0 +1,71 @@
+
+import numpy as np
+import pylab as pl
+ 
+from model import run_model, perturbation
+from figures import plot_timeseries
+
+def plot_perturbation(wind=7, mldp=30, pref="ca"):
+    """Run the model in pertubation mode and plot the results"""
+    pl.clf()
+    fig,axes = pl.subplots(3, 1, sharex=True, num=1,figsize=(6,6))
+    model_kws = dict(pref=pref, reg="pert",
+                temp=None, salt=None, wind=wind, mldp=mldp,
+                deepmld1={"pert":0, "tpos":165},
+                deepmld2={"pert":0, "tpos":553},
+                uppwell1={"pert":0, "tpos":165, "pdays":5},
+                uppwell2={"pert":0, "tpos":553, "pdays":5})
+    md = run_model(**model_kws)
+    plot_timeseries(md, axes=axes, alpha=0.5) 
+    if pref == "ca":
+        preftxt = "CCS"
+        model_kws["uppwell1"]["pert"] = 82.5
+        model_kws["uppwell2"]["pert"] = 165
+    else:
+        preftxt = "NWA"
+        model_kws["deepmld1"]["pert"] = 17
+        model_kws["deepmld2"]["pert"] = 34
+    md = run_model(**model_kws)
+    plot_timeseries(md, axes=axes, alpha=1)
+    pl.suptitle(
+        f"Perturbations, temp and salt for {preftxt}, wind:{wind}m/s, mld:{mldp}m")
+    pl.savefig(f"figs/pertubation_timeseries_{pref}.pdf")
+
+def to_csv(wind=7, mldp=30, pref="ca"):
+    """Save model results to csv
+    
+    Parameters
+    ----------
+    wind : float, optional
+        Prescribed fixed wind
+    mldp : float, optional
+        Prescribed fixed MLD
+    pref : string, optional
+        Part of the coast [ne=North Atlantic outside Gulf of Maine]
+                           ca=California Coastal Current]    
+    """
+    ptext = "CCS" if pref == "ca" else "NWA"
+    if pref == "ca":
+        preftxt = "CCS"
+        model_kws = dict(pref=pref, reg="pert",
+                     temp=None, salt=None, wind=wind, mldp=mldp,
+                     deepmld1={"pert":0, "tpos":165},
+                     deepmld2={"pert":0, "tpos":553},
+                     uppwell1={"pert":82.5, "tpos":165, "pdays":5},
+                     uppwell2={"pert":165,  "tpos":553, "pdays":5})
+    else:
+        preftxt = "NWA"
+        model_kws = dict(pref=pref, reg="pert",
+                     temp=None, salt=None, wind=wind, mldp=mldp,
+                     deepmld1={"pert":17, "tpos":165},
+                     deepmld2={"pert":34, "tpos":553},
+                     uppwell1={"pert":0,  "tpos":165, "pdays":5},
+                     uppwell2={"pert":0,  "tpos":553, "pdays":5})
+
+    md = run_model(**model_kws)
+    md.to_csv(f"boxmodel_baseline_{ptext}__wind_{wind:02}__mldp_{mldp:03}.csv")
+    for key in ["deepmld1", "deepmld1", "uppwell1", "uppwell2"]:
+        model_kws[key] = {"pert":0, "tpos":165}
+    md = run_model(**model_kws)
+    md.to_csv(f"boxmodel_pertubation_{ptext}__wind_{wind:02}__mldp_{mldp:03}.csv")
+    
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..47f3f41
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,31 @@
+[tool.poetry]
+name = "pco2box"
+version = "1.0"
+description = "Box model simulating carbon solubility and air-sea exchange in the ocean"
+authors = ["Bror Jonsson <brorfred@gmail.com>"]
+license = "MIT"
+packages = [
+    { include = "oceanbox" },
+    { include = "pco2box" },
+]
+
+
+[tool.poetry.dependencies]
+python = "^3.7"
+numpy = "^1.17"
+scipy = "^1.3"
+co2sys = "^0.1.1"
+pandas = "^0.25.3"
+matplotlib = "^3.1"
+seawater = "^3.3"
+seaborn = "^0.9.0"
+#cartopy = "^0.17.0"
+
+tables = "^3.6"
+[tool.poetry.dev-dependencies]
+ipython = "^7.9"
+pytest = "^5.2"
+
+[build-system]
+requires = ["poetry>=0.12"]
+build-backend = "poetry.masonry.api"
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/test_model.py b/tests/test_model.py
new file mode 100644
index 0000000..192b8cc
--- /dev/null
+++ b/tests/test_model.py
@@ -0,0 +1,7 @@
+
+from pco2_east_west import model
+
+def test_setup():
+
+    md = model.BoxModel(svec=730)
+    md.run()