def refstate_meandensity(rho_in,grd,area_xyz_in):
    # def 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