diff --git a/bgcval2/bgcval2_make_report.py b/bgcval2/bgcval2_make_report.py
index 812c555..c6daa22 100755
--- a/bgcval2/bgcval2_make_report.py
+++ b/bgcval2/bgcval2_make_report.py
@@ -1480,6 +1480,7 @@ def newImageLocation(fn):
     physicsKM = [
         'AMOC_26N',
         'ADRC_26N',
+        'GulfStream',
         'AtmosCO2',  
         'DrakePassage',
         'GlobalMeanTemperature',
diff --git a/bgcval2/bgcvaltools/pftnames.py b/bgcval2/bgcvaltools/pftnames.py
index bd233dd..58d8b3a 100644
--- a/bgcval2/bgcvaltools/pftnames.py
+++ b/bgcval2/bgcvaltools/pftnames.py
@@ -235,6 +235,8 @@ def makeLongNameDict():
     lnd['AMOC_26N_nomexico'] = "AMOC 26N (excluding Gulf of Mexico)"
     lnd['AMOC_32S'] = "AMOC 32S"
     lnd['ADRC_26N'] = "Atlantic Deep Return Current at 26N"
+    lnd['GulfStream'] = "Gulf Stream at 26N"
+
 
     lnd['ZonalCurrent'] = "Zonal Current"
     lnd['MeridionalCurrent'] = "Meridional Current"
diff --git a/bgcval2/functions/circulation.py b/bgcval2/functions/circulation.py
index 0e879f4..464553e 100644
--- a/bgcval2/functions/circulation.py
+++ b/bgcval2/functions/circulation.py
@@ -62,6 +62,9 @@
 
 #eORCA1_latslice26N = slice(227,228)
 eORCA1_latslice26Nnm = slice(228,229)
+eORCA1_lonslice_GS = slice(228,229)
+
+
 eORCA1_latslice32S = slice(137,138)
 
 eORCA025_latslice26Nnm = slice(794,795)
@@ -312,6 +315,72 @@ def TwentySixNorth(nc,keys,**kwargs):
     return atlmoc
 
 
+def gulfstream(nc, keys, **kwargs):
+    """
+    This function calculates the Gulf Stream at 26N.
+    This is the sum of Northbound current between 
+
+    nc: a netcdf openned as a dataset.
+    keys: a list of keys to use in this function.
+
+    """
+    areafile = get_kwarg_file(kwargs, 'areafile')
+    maskname = kwargs.get('maskname', 'tmask')
+    grid = kwargs.get('grid', 'eORCA1')
+    maxdepth = np.abs(kwargs.get('maxdepth', 2000. ))
+
+ 
+    if not loadedArea:
+        loadDataMask(areafile, maskname, grid)
+
+    if grid == 'eORCA1':
+        latslice26Nnm = eORCA1_latslice26Nnm
+
+        #data=[-80.5011659 , -79.50119298, -78.50121829, -77.50124181,
+        #           -76.50126349, -75.50128329, -74.50130118, -73.50131712,
+        #           -72.50133107, -71.50134301, -70.50135293, -69.50136079,
+        #           -68.50136658],
+        lonslice_70W = slice(207, 220) 
+
+        altmaskfile = get_kwarg_file(kwargs, 'altmaskfile', default = 'bgcval2/data/basinlandmask_eORCA1.nc')
+        if not loadedAltMask:
+             loadAtlanticMask(altmaskfile, maskname='tmaskatl', grid=grid)
+    elif grid == 'eORCA025':
+        latslice26Nnm = eORCA025_latslice26Nnm
+
+    else:
+        # grid not recognised.
+        raise ValueError('TwentySixNorth: grid not recognised: %s', grid)
+
+    if not loadedAltMask:
+        # Atlantic Mask not loaded
+        raise ValueError('TwentySixNorth: Mask not loaded: ed: %s', grid)
+        assert 0
+
+    lats = nc.variables['nav_lat'][latslice26Nnm, lonslice_70W]
+    lons = nc.variables['nav_lon'][latslice26Nnm, lonslice_70W]
+    vo = np.ma.array(nc.variables[keys[0]][0, :, latslice26Nnm, lonslice_70W]) # m/s
+    vo = np.ma.masked_where(vo.mask + (vo <= 0.), vo) 
+
+    thickness = nc.variables['thkcello'][0,:,latslice26Nnm, lonslice_70W] 
+    depth = np.abs(np.cumsum(thickness, axis=0))# depth array
+    #print(vo.shape, thickness.shape, e1v_AMOC26N.shape)
+    gs = 0.
+    for (z, la, lo), v in np.ndenumerate(vo):
+        if depth[z, la,lo] > maxdepth:
+            continue
+        if v <= 0:
+            continue
+        #print((z, la, lo),'depth:', depth[z, la,lo], (lats[la, lo],'N', lons[la, lo], 'E'),  'v:', v, 'thickness:', thickness[z, la, lo], 'width:', e1v_AMOC26N[la, lo])
+        gs += v * thickness[z, la, lo] * e1v_AMOC26N[la, lo] / 1.E06
+
+    print('Gulf Stream:', gs) # expecting a value of 32Sv ish.
+    # https://www.sciencedirect.com/science/article/pii/S0079661114001694
+
+    return gs
+
+
+
 def twentysixnorth025(nc,keys,**kwargs):
     """
     This function loads the AMOC array that is used for eORCA025
diff --git a/bgcval2/timeseries/timeseriesAnalysis.py b/bgcval2/timeseries/timeseriesAnalysis.py
index 1a5974d..2dfff87 100644
--- a/bgcval2/timeseries/timeseriesAnalysis.py
+++ b/bgcval2/timeseries/timeseriesAnalysis.py
@@ -139,10 +139,11 @@ def loadModel(self):
         # load and calculate the model info
         if len(glob.glob(self.shelvefn+'*')): # shelve files have .bak .dat .dir files now
             with shOpen(self.shelvefn) as sh:
-                print('Shelf opens fine:', self.shelvefn)
+                print('Opening Shelf file:', self.shelvefn)
                 print (sh.keys())
                 readFiles       = sh['readFiles']
                 modeldataD      = sh['modeldata']
+               
         else:
             print('Does not exist', self.shelvefn)
             readFiles = []
diff --git a/key_lists/physics.yml b/key_lists/physics.yml
index fd20f5b..80760a4 100644
--- a/key_lists/physics.yml
+++ b/key_lists/physics.yml
@@ -51,7 +51,7 @@ keys:
     #	ZonalCurrent: True             	# Zonal Veloctity
     #	MeridionalCurrent: True        	# Meridional Veloctity
     #	VerticalCurrent: True          	# Vertical Veloctity
-
+    GulfStream: True
 
     somesatr: True # salt flow # Work in progress in debug