diff --git a/calc_ref_state.py b/calc_ref_state.py
new file mode 100644
index 0000000..475064f
--- /dev/null
+++ b/calc_ref_state.py
@@ -0,0 +1,134 @@
+# ############################################################################################
+# Program name  : calc_ref_state.py
+# written by    : Gabriel Wolf, g.a.wolf@reading.ac.uk
+#                 adapted from compute_reference_state_woa2009_for_Juan.m of Remi Tailleux
+# last modified : 13.09.2018
+# ############################################################################################
+
+# Load necessary modules #####################################################################
+import numpy as np      # 
+import time             # to measure elapsed time
+# ############################################################################################
+
+# Input 0 ####################################################################################
+# Manual defined inputs
+# input         : None 
+# output        : dir_plot
+#     dir_plot  : directory, where plots will be saved 
+# MANUAL DEFINED OUTPUTS MISSING
+# ############################################################################################
+print 'MANUAL DEFINED OUTPUTS MISSING'
+
+# Input 1 ####################################################################################
+# Physical constants
+# output        : r_earth
+#     r_earth   : earth radis           [m]
+# PHYSICAL CONSTATS MISSING
+# ############################################################################################
+print 'PHYSICAL CONSTATS MISSING'
+
+# step a #####################################################################################
+# Read necessary data to calculate reference state
+# input         : None 
+# output        : salt, theta, lon, lat, zlev, valid
+#     salt      : salinity              [psu]
+#     theta     : potential temperature [degC]
+#     lon       : longitude             [degrees]
+#     lat       : latitude              [degrees]
+#     zlev      : depth of ocean        [m] ; values > 0
+#     valid     : mask for gridpoints with ocean data
+# ############################################################################################
+beg_time = time.time()
+print ''
+print ' ----------------------------------------------'
+print 'Insert code to read data'
+print ' ----------------------------------------------'
+end_time = time.time()
+print 'Elapsed time to read data:         '+ '{:.2f}'.format(end_time-beg_time)+'s'
+print ' ----------------------------------------------'
+
+# step b #####################################################################################
+# Define grid for data from step a
+# input         : lon, lat, zlev (all from step a) 
+# output        : plev
+#     plev      : MISSING DEFINITION
+# ############################################################################################
+beg_time = time.time()
+print ''
+print ' ----------------------------------------------'
+print 'Insert code to define grid of data'
+print ' ----------------------------------------------'
+end_time = time.time()
+print 'Elapsed time to define grid:       '+ '{:.2f}'.format(end_time-beg_time)+'s'
+print ' ----------------------------------------------'
+
+# step c #####################################################################################
+# Calculation of density
+# input         : salt, theta, plev, valid
+#     salt      : output of step a
+#     theta     : output of step a
+#     p_lev     : output of step b
+# output        : dens
+#     dens      : density of water      [kg m**-3]
+# ############################################################################################
+beg_time = time.time()
+print ''
+print ' ----------------------------------------------'
+print 'Insert code to calculate density'
+print ' ----------------------------------------------'
+end_time = time.time()
+print 'Elapsed time to calculate density: '+ '{:.2f}'.format(end_time-beg_time)+'s'
+print ' ----------------------------------------------'
+
+# step d #####################################################################################
+# Calculate z-dependent ocean area
+# input         : lat, r_earth,nlon, nlat, nz
+#     lat       : output of step a
+#     r_earth   : defined in Input 1
+#     nlon, nlat, nz : MISSING - PUT THIS INTO GRD?
+# output        : area_xyz
+#     area_xyz  : z-dep ocean area      [m**2]
+# ############################################################################################
+beg_time = time.time()
+print ''
+print ' ----------------------------------------------'
+print 'Insert code to calculate ocean area'
+print ' ----------------------------------------------'
+end_time = time.time()
+print 'Elapsed time to calc. ocean area:  '+ '{:.2f}'.format(end_time-beg_time)+'s'
+print ' ----------------------------------------------'
+
+# step e #####################################################################################
+# Calculate reference state
+# input         : dens, zlev, area_xyz
+#     dens      : output of step c
+#     zlev      : output of step a
+#     area_xyz  : output of step d
+# otput         : pp_rhor
+#     pp_rhor   : interpolant of reference density
+# ############################################################################################
+beg_time = time.time()
+print ''
+print ' ----------------------------------------------'
+print 'Insert code to calculate reference state'
+print ' ----------------------------------------------'
+end_time = time.time()
+print 'Elapsed time to calculate reference state: '+ '{:.2f}'.format(end_time-beg_time)+'s'
+print ' ----------------------------------------------'
+
+# step f #####################################################################################
+# plot reference state
+# input         : MISSING
+# output        : None
+# saved/created : plots saved in dir_plot
+#     dir_plot  : output of Input 0 
+# ############################################################################################
+beg_time = time.time()
+print ''
+print ' ----------------------------------------------'
+print 'Insert code to plot reference state'
+print ' ----------------------------------------------'
+end_time = time.time()
+print 'Elapsed time to plot and save reference state: '+ '{:.2f}'.format(end_time-beg_time)+'s'
+print ' ----------------------------------------------'
+