def refstate_meandensity(rho_in,grd,area_xyz_in):
    # refstate_meandensity(rho,zlev,area_xyz) ################################################
    # written by    : Gabriel Wolf, g.a.wolf@reading.ac.uk
    #                 adapted from get_refstate_meandensity.m of Remi Tailleux (10.04.2013)
    # last modified : 20.09.2018
    # Content/Description ####################################################################
    # Usage         : from gw_ocean_refstate import refstate_meandensity
    # Input         : area_xyz, rho, zlev
    #     area_xyz  : hor. area for 3d grid  [m**2], 3d
    #     rho       : density                [kg m**-3], 3d
    #     grd       : Input grid             [python class]
    # Output        : pp_rhor               
    #     pp_rhor   : interpolant for reference state
    # ########################################################################################
    # Load modules
    import numpy as np
    from scipy import interpolate
    # Allocation
    rhor = np.zeros([grd.Nz,1])
    # compute hor. averaged density field
    rho_surf = area_xyz_in*rho_in
    mask     = np.isfinite(rho_surf)
    rhor     = sum(sum(mask*rho_surf,1),0)/sum(sum(area_xyz_in*mask,1),0)
    # define arrays for interpolation of density profiles
    pp_rhor  = interpolate.PchipInterpolator(grd.z, rhor, axis=0)
    # return interpolant for reference state
    return pp_rhor

def refstate_sorted_stheta(th,s,rho,mask,ztarget,grd,sbin=[0.0,43.0,0.002],\
    thbin=[-2.6,30.0,0.005]):
    # refstate_sorted_stheta(th,s,rho,mask,ztarget,grd) ######################################
    # written by    : Gabriel Wolf, g.a.wolf@reading.ac.uk
    #                 adapted from get_refstate_sorted.m of Remi Tailleux (10.04.2013)
    # last modified : 20.09.2018
    # Content/Description ####################################################################
    # Reference     : method based on Tailleux (2013), -> https://doi.org/10.1017/jfm.2013.509
    #               : summary of this framework in Saenz et al. (2015):
    #                 -> https://doi.org/10.1175/JPO-D-14-0105.1
    # Usage         : from gw_ocean_refstate import refstate_sorted_stheta
    # Input         : grd, mask, rho, s, sbin, thbin, th, ztarget
    #     grd       : Input grid             [python class]
    #     mask      : valid points           [boolean], 3d
    #     rho       : density                [kg m**-3], 3d
    #     s         : salinity               [psu], 3d
    #     s_bin     : bin range for s to construct s-th-diagram [min,max,delta]
    #     th        : potential temperature  [deg C], 3d
    #     th_bin    : bin range for th to construct s-th-diagram [min,max,delta]
    #     ztarget   : target height level    [m], 1d
    # Output        :pp_rhor_botup, pp_rhor_topdn
    #     pp_rhor_botup : interpolant for density reference state (calc. from bottom up)
    #     pp_rhor_topdn : interpolant for density reference state (calc. from bottom up)
    # ########################################################################################
    # Load modules
    import numpy as np
    from scipy import interpolate                  # for interpolation
    from mycalc_ocean import volume_census_levitus # to create s-th-pdf with vol info
    # Allocation
    n_target = len(ztarget)
    s[mask==False]   = np.nan
    th[mask==False]  = np.nan
    rho[mask==False] = np.nan
    # Compute the pdf of the salinity/theta data ---------------------------------------------
    print '  *** MISSING: computation of the pdf salinity/theta data ***'
    pdf_vol_lev, vol_tot_ocean = volume_census_levitus(th, s, grd, thbin, sbin)
    # Define properties of the binned temperature/salinity data ------------------------------
    #n_sbin, n_thbin = 
    print '  *** MISSING: property definition of binned temperature/salinity data ***'
    # define target pressure from target depths ----------------------------------------------
    print '  *** MISSING: define target pressure from target depths ***'
    # compute bounding salinity curves -------------------------------------------------------
    print '  *** MISSING: compute bounding salinity curves ***'
    # Seek the global min. and max. of density by bringing all parcels at the surface and at 
    #     the bottom and the oceans ----------------------------------------------------------
    print '  *** MISSING: density profiles ***'
    # Compute potential density reference to target pressure at surface ----------------------
    print '  *** MISSING: pot. density at surface ***'
    # Compute potential density referenced to last target pressure ---------------------------
    print '  *** MISSING: pot. density for last target pressure ***'
    # Compute global density minimum and maximum ---------------------------------------------
    print '  *** MISSING: min. and max. of global density ***'
    # Compute upper and lower salinity curves in theta/salt space ----------------------------
    print '  *** MISSING: upper and lower salinity curves in th-s space ***'
    # Compute the topography statistics ------------------------------------------------------
    print '  *** MISSING: topography statistics ***'
    # New attempt to compute the top-down profile --------------------------------------------
    print '  *** MISSING: rhor top-down profile ***'
    print('  Computing top-dwon profile for reference density')
    # Allocation
    rhor_topdn = np.zeros(n_target)
    # New attempt at computing bottom-up profile ---------------------------------------------
    print '  *** MISSING: rhor bottom-up profile ***'
    print('  Computing bottom-up profile for reference density')
    # Allocation
    rhor_botup = np.zeros(n_target)
    # Assign interpolants for the two possible reference states ------------------------------
    pp_rhor_botup = interpolate.PchipInterpolator(ztarget, rhor_botup, axis=0)
    pp_rhor_topdn = interpolate.PchipInterpolator(ztarget, rhor_topdn, axis=0)
    # return interpolant for the two reference states ----------------------------------------
    #return pp_rhor_botup, pp_rhor_topdn
    return pdf_vol_lev, vol_tot_ocean