From a6fbce9822a2de06b99ade525f6a52441e378ecb Mon Sep 17 00:00:00 2001
From: Gabriel Wolf <g.a.wolf@reading.ac.uk>
Date: Fri, 21 Sep 2018 11:40:48 +0100
Subject: [PATCH] Update calc_ref_state.py and add gw_ocean_refstate.py to
 calculate mean density ref. state

---
 gw_ocean_refstate.py | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)
 create mode 100644 gw_ocean_refstate.py

diff --git a/gw_ocean_refstate.py b/gw_ocean_refstate.py
new file mode 100644
index 0000000..0041ec7
--- /dev/null
+++ b/gw_ocean_refstate.py
@@ -0,0 +1,27 @@
+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