diff --git a/calc_ref_state.py b/calc_ref_state.py
index f3d716e..e8777ae 100644
--- a/calc_ref_state.py
+++ b/calc_ref_state.py
@@ -29,6 +29,7 @@ import time             # to measure elapsed time
 from get_data import get_data_woa2009 as read_data
 # calculation modules
 from mycalc_ocean import calc_rho_from_theta # calc of density from theta
+from mycalc_ocean import calc_area_lat # calc horizontal area
 # plotting modules / functions
 from myplot import myplot_2dmap # plot of lat-lon fields
 from myplot import myplot_2d    # plot of crossections fields
@@ -64,7 +65,7 @@ cs_plot    = {'lon':[-47.5,172.5],'lat':[0]} # crosssections
 #     grav      : gravity               [m s**-2]
 #     gbuo      : 
 # ############################################################################################
-from myconstants import rho0, grav
+from myconstants import rho0, grav, r_earth
 
 # step a #####################################################################################
 # Read necessary data to calculate reference state
@@ -160,7 +161,8 @@ print '----------------------------------------------'
 beg_time = time.time()
 print ''
 print '----------------------------------------------'
-print 'Insert code to calculate ocean area'
+print 'Calculate z-dependent ocean area'
+area_xyz = calc_area_lat(grd,r_earth)
 print '----------------------------------------------'
 end_time = time.time()
 print 'Elapsed time to calc. ocean area:  '+ '{:.2f}'.format(end_time-beg_time)+'s'
diff --git a/mycalc_ocean.py b/mycalc_ocean.py
index e9dc226..8bc5454 100644
--- a/mycalc_ocean.py
+++ b/mycalc_ocean.py
@@ -1,4 +1,38 @@
-def calc_rho_from_theta(s, th, p, reference='jackettetal2004')
+# calc_area_lat
+# calc_rho_from_theta
+# calc_theta_from_temp
+# gsw_g
+def calc_area_lat(grd,r_earth):
+    # 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)
+    # last modified : 19.09.2018
+    # Content/Description ####################################################################
+    # Input         : grd, r_earth
+    #     grd       : Input grid             [python class]
+    #     r_earth   : earth radius           [m]
+    # Output        : area_lat               
+    #     area_lat  : hor. area for 3d grid  [m**2]
+    # ########################################################################################
+    # import modules
+    import numpy as np
+    # 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)
+    # return horizontal area for input grid
+    return area_xyz
+
+def calc_rho_from_theta(s, th, p, reference='jackettetal2004'):
     # calc_rho_from_theta(s,th,p,reference='jackettetal2004') ################################
     # written by    : Gabriel Wolf, g.a.wolf@reading.ac.uk
     #                 adapted from rho_from_theta.m of Remi Tailleux (09/2018)
@@ -72,7 +106,7 @@ def calc_rho_from_theta(s, th, p, reference='jackettetal2004')
     # Deallocation of Variables
     # ########################################################################################
     del th2, sqrts, pth
-
+    # return density
     return rho