diff --git a/calc_ref_state.py b/calc_ref_state.py
index e8777ae..4da617b 100644
--- a/calc_ref_state.py
+++ b/calc_ref_state.py
@@ -109,6 +109,10 @@ except:
     grd.p  = plev
     grd.Up = 'dbar'
     grd.Np = len(plev)
+if grd.lon[-1]+grd.lon[0] >= 1:
+   grd.symlon = 1
+else:
+   grd.symlon = 0
 lat_xy, lon_xy          = np.meshgrid(grd.lat,grd.lon)
 p_yz, lat_yz            = np.meshgrid(grd.z, grd.lat)
 lat_xyz, lon_xyz, p_xyz = np.meshgrid(grd.lat, grd.lon, grd.p)
diff --git a/mycalc_ocean.py b/mycalc_ocean.py
index 8bc5454..a2139f0 100644
--- a/mycalc_ocean.py
+++ b/mycalc_ocean.py
@@ -2,7 +2,7 @@
 # calc_rho_from_theta
 # calc_theta_from_temp
 # gsw_g
-def calc_area_lat(grd,r_earth):
+def calc_area_lat(grd,r_earth,calc_old=0):
     # def calc_area_lat(grd,r_earth) #########################################################
     # written by    : Gabriel Wolf, g.a.wolf@reading.ac.uk
     #                 adapted from get_area_latitude.m of Remi Tailleux (09/2018)
@@ -19,16 +19,45 @@ def calc_area_lat(grd,r_earth):
     # convert deg to rad for sin/cos
     deg_to_rad  = np.pi/180.0
     lat_radians = grd.lat*deg_to_rad
-    # compute length of longitudinal and latitudinal elements
-    print 'area calculation not valid for non-equidistant input lon-grid'
-    print 'area calculation not valid for non-equidistant input lat-grid'
-    dx = r_earth*deg_to_rad*(grd.Nlon/360.0)
-    dy = r_earth*deg_to_rad*abs(grd.lat[1]-grd.lat[0])
-    # calculate area
-    area_xyz = np.zeros([grd.Nlon,grd.Nlat,grd.Nz])
-    for i_lon in range(0,grd.Nlon):
-        for i_z in range(0,grd.Nz):
-            area_xyz[i_lon,:,i_z] = (dx*dy)*np.cos(lat_radians)
+    if calc_old == 1:
+        print '*** old version of area calculation ***'
+        # compute length of longitudinal and latitudinal elements
+        print '*** area calculation not valid for non-equidistant input lon-grid ***'
+        print '*** area calculation not valid for non-equidistant input lat-grid ***'
+        dx = r_earth*deg_to_rad*(grd.Nlon/360.0)
+        dy = r_earth*deg_to_rad*abs(grd.lat[1]-grd.lat[0])
+        # calculate area
+        area_xyz = np.zeros([grd.Nlon,grd.Nlat,grd.Nz])
+        for i_lon in range(0,grd.Nlon):
+            for i_z in range(0,grd.Nz):
+                area_xyz[i_lon,:,i_z] = (dx*dy)*np.cos(lat_radians)
+    else:
+        lat_d, lon_d, z_d    = np.meshgrid(grd.lat,grd.lon,grd.z)
+        if grd.symlon == 1:
+            delta_lon = 360.0
+            index_d       = np.append(np.arange(1,grd.Nlon),0)
+            lon_e         = 0.5 * (lon_d + lon_d[index_d,:,:])
+            lon_e[-1,:,:] = lon_e[-1,:,:] + 0.5 * delta_lon 
+            index_d       = np.append(grd.Nlon-1,np.arange(0,grd.Nlon-1))
+            lon_w         = 0.5 * (lon_d + lon_d[index_d,:,:])
+            lon_w[0,:,:]  = lon_w[0,:,:] - 0.5 * delta_lon
+        else:
+            lon_d   = np.append(2*grd.lon[0] - grd.lon[1],grd.lon)
+            lon_d   = np.append(lon_d,2*grd.lon[-1] - grd.lon[-2])
+            dummy1, lon_d, dummy2 = np.meshgrid(grd.lat,lon_d,grd.z)
+            index_d = np.arange(1,grd.Nlon+1)
+            lon_e   = 0.5 * (lon_d[index_d,:,:] + delta_lon + lon_d[index_d+1,:,:])
+            lon_w   = 0.5 * (lon_d[index_d,:,:] + delta_lon + lon_d[index_d-1,:,:])
+        lat_d   = np.append(max(2*grd.lat[0] - grd.lat[1],-180.0-grd.lat[0]),grd.lat)
+        lat_d   = np.append(lat_d,min(2*grd.lat[-1] - grd.lat[-2],180.0-grd.lat[-1]))
+        lat_d, dummy1, dummy2 = np.meshgrid(lat_d,grd.lon,grd.z)
+        index_d = np.arange(1,grd.Nlat+1)
+        lat_n   = 0.5 * (lat_d[:,index_d,:] + lat_d[:,index_d+1,:])
+        lat_s   = 0.5 * (lat_d[:,index_d,:] + lat_d[:,index_d-1,:])
+        r_z     = r_earth - abs(grd.z)
+        #area_xyz  = (deg_to_rad*r_z)**2*(lon_e-lon_w)*(lat_n-lat_s)*np.cos(lat_d[:,index_d,:]*deg_to_rad)
+        area_xyz  = deg_to_rad*r_z**2*(lon_e-lon_w)*(np.sin(lat_n*deg_to_rad)-np.sin(lat_s*deg_to_rad))
+        del delta_lon, lon_d, lat_d, z_d, dummy1, dummy2, index_d, r_z
     # return horizontal area for input grid
     return area_xyz