From ad47df4f60e169051ac6b1955bf0173e94c50753 Mon Sep 17 00:00:00 2001 From: Kaitlin Alexander Date: Fri, 11 Dec 2015 13:29:02 +1100 Subject: [PATCH] Added density calculation routines --- file_guide.txt | 32 ++++----- make_density_file.py | 85 +++++++++++++++++++++++ plot_ohc.py | 119 ------------------------------- plot_tke.py | 162 ------------------------------------------- plot_totalsalt.py | 115 ------------------------------ unesco.py | 74 ++++++++++++++++++++ 6 files changed, 173 insertions(+), 414 deletions(-) create mode 100644 make_density_file.py delete mode 100644 plot_ohc.py delete mode 100644 plot_tke.py delete mode 100644 plot_totalsalt.py create mode 100644 unesco.py diff --git a/file_guide.txt b/file_guide.txt index c3922d6..b19f0e5 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -241,25 +241,21 @@ plot_volume.py: Extracts the volume values written in the ocean.log files (you timestep length and the output frequency for ocean.log (INFOSTEP in Params.py if you are using ROMS-CICE-MCT). -plot_ohc.py: Calculates and plots the total heat content at each timestep of - the given ocean averages file. This file also contains useful code - for calculating volume averages. - To run: Open python or ipython, and type "run plot_ohc.py". - The script will prompt you for paths to the ocean averages - file and the grid file. +unesco.py: Calculates the UNESCO seawater equation of state (1980): given + temperature, salinity, and pressure (depth/10 is fine), returns + density. + To run: The function unesco is designed to be called by another + script. See make_density_file.py for an example. -plot_totalsalt.py: Calculates and plots the total salt content at each timestep - of the given ocean averages file. - To run: Open python or ipython, and type - "run plot_totalsalt.py". The script will prompt you - for paths to the ocean averages file and the grid - file. - -plot_tke.py: Calculates and plots the total kinetic energy at each timestep of - the given ocean averages file. - To run: Open python or ipython, and type "run plot_tke.py". The - script will prompt you for paths to the ocean averages file - and the grid file. +make_density_file.py: Given an ocean history or averages file with temperature + and salinity data, calculate density fields at each + timestep using the 1980 UENSCO seawater equation of + state. Save in a new file. + To run: Open python or ipython, and type + "run make_density_file.py". The script will + prompt you for the paths to the ROMS grid file, + the ocean history or averages file, and the + desired path to the new density file. ***NICE FIGURES*** diff --git a/make_density_file.py b/make_density_file.py new file mode 100644 index 0000000..6ee3d6a --- /dev/null +++ b/make_density_file.py @@ -0,0 +1,85 @@ +from netCDF4 import Dataset +from numpy import * +from calc_z import * +from unesco import * + +# Given an ocean history or averages file with temperature and salinity data, +# calculate density fields at each timestep using the 1980 UNESCO seawater +# equation of state. Save in a new file. +# Input: +# grid_file = path to ROMS grid file +# input_file = path to ocean history/averages file +# output_file = desired path to new density file +def make_density_file (grid_file, input_file, output_file): + + # Grid parameters + theta_s = 0.9 + theta_b = 4.0 + hc = 40 + N = 31 + + # Read grid variables + id = Dataset(grid_file, 'r') + h = id.variables['h'][:,:] + zice = id.variables['zice'][:,:] + lon = id.variables['lon_rho'][:,:] + lat = id.variables['lat_rho'][:,:] + id.close() + num_lon = size(lon, 1) + num_lat = size(lon, 0) + + # Get a 3D array of z-coordinates (metres) + z, sc_r, Cs_r = calc_z(h, zice, lon, lat, theta_s, theta_b, hc, N) + # Pressure is approximately equal to |z|/10 + press = abs(z)/10.0 + + # Set up output file + out_id = Dataset(output_file, 'w') + # Define dimensions + out_id.createDimension('xi_rho', num_lon) + out_id.createDimension('eta_rho', num_lat) + out_id.createDimension('s_rho', N) + out_id.createDimension('ocean_time', None) + # Define variables + out_id.createVariable('lon_rho', 'f8', ('eta_rho', 'xi_rho')) + out_id.variables['lon_rho'][:,:] = lon + out_id.createVariable('lat_rho', 'f8', ('eta_rho', 'xi_rho')) + out_id.variables['lat_rho'][:,:] = lat + out_id.createVariable('sc_r', 'f8', ('s_rho')) + out_id.variables['sc_r'].long_name = 'S-coordinate at rho-points' + out_id.variables['sc_r'][:] = sc_r + out_id.createVariable('ocean_time', 'f8', ('ocean_time')) + out_id.variables['ocean_time'].units = 'seconds' + out_id.createVariable('rho', 'f8', ('ocean_time', 's_rho', 'eta_rho', 'xi_rho')) + out_id.variables['rho'].long_name = 'density' + out_id.variables['rho'].units = 'kg/m^3' + + # Read time values from input file + in_id = Dataset(input_file, 'r') + time = in_id.variables['ocean_time'][:] + + # Process each timestep individually to conserve memory + for t in range(size(time)): + print 'Processing timestep '+str(t+1)+' of '+str(size(time)) + # Set a new time value in the output file + out_id.variables['ocean_time'][t] = time[t] + # Read temperature and salinity (convert to float128 to prevent + # overflow in UNESCO calculations) + temp = array(in_id.variables['temp'][t,:,:,:], dtype=float128) + salt = array(in_id.variables['salt'][t,:,:,:], dtype=float128) + # Magic happens here + rho = unesco(temp, salt, press) + # Save the results for this timestep + out_id.variables['rho'][t,:,:,:] = rho + + in_id.close() + out_id.close() + + +# Command-line interface +if __name__ == "__main__": + + grid_file = raw_input("Path to grid file: ") + input_file = raw_input("Path to ocean history/averages file: ") + output_file = raw_input("Desired path to new density file: ") + make_density_file(grid_file, input_file, output_file) diff --git a/plot_ohc.py b/plot_ohc.py deleted file mode 100644 index 133e645..0000000 --- a/plot_ohc.py +++ /dev/null @@ -1,119 +0,0 @@ -from netCDF4 import Dataset -from numpy import * -from matplotlib.pyplot import * -from calc_z import * - -# Calculate and plot the total heat content at each timestep of the given ocean -# averages file. -def plot_ohc (file_path, grid_path): - - # Grid parameters - theta_s = 0.9 - theta_b = 4.0 - hc = 40 - N = 31 - # Radius of the Earth in m - r = 6.371e6 - # Degrees to radians conversion factor - deg2rad = pi/180.0 - # Specific heat of polar seawater (J/K/kg) - cp = 3974.0 - # Celsius to Kelvin conversion constant - celsius2kelvin = 273.15 - - # Read grid variables - id = Dataset(grid_path, 'r') - h = id.variables['h'][:,:] - zice = id.variables['zice'][:,:] - lon = id.variables['lon_rho'][:,:] - lat = id.variables['lat_rho'][:,:] - mask = id.variables['mask_rho'][:,:] - id.close() - - # Mask lat and lon at land points - lon = ma.masked_where(mask==0, lon) - lat = ma.masked_where(mask==0, lat) - # Save dimensions - num_lat = size(lon, 0) - num_lon = size(lon, 1) - - # Add or subtract 360 from longitude values which wrap around - # so that longitude increases monotonically from west to east - i = tile(arange(1, num_lon+1), (num_lat, 1)) - index1 = nonzero((i > 1200)*(lon < 100)) - lon[index1] = lon[index1] + 360 - index2 = nonzero((i < 200)*(lon > 300)) - lon[index2] = lon[index2] - 360 - - # Interpolate to get longitude at the edges of each cell - w_bdry = (lon[:,0] + lon[:,-1] - 360)/2 - middle_lon = (lon[:,0:-1] + lon[:,1:])/2 - e_bdry = (lon[:,0] + 360 + lon[:,-1])/2 - lon_edges = ma.concatenate((w_bdry[:,None], middle_lon, e_bdry[:,None]), axis=1) - # Subtract to get the change in longitude over each cell - dlon = abs(lon_edges[:,1:] - lon_edges[:,0:-1]) - - # Similarly for latitude - s_bdry = lat[0,:] - middle_lat = (lat[0:-1,:] + lat[1:,:])/2 - n_bdry = lat[-1,:]*0 - 50 - lat_edges = ma.concatenate((s_bdry[None,:], middle_lat, n_bdry[None,:])) - dlat = lat_edges[1:,:] - lat_edges[0:-1,:] - - # Convert from spherical to Cartesian coordinates - # dy = r*dlat where dlat is converted to radians - dy_2d = r*dlat*pi/180.0 - # Copy into a 3D array, same at each depth level - dy = tile(dy_2d, (N,1,1)) - # dx = r*cos(lat)*dlon where lat and dlon are converted to radians - dx_2d = r*cos(pi*lat/180.0)*dlon*pi/180.0 - dx = tile(dx_2d, (N,1,1)) - - # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script - z, sc_r, Cs_r = calc_z(h, zice, lon, lat, theta_s, theta_b, hc, N) - # We have z at the midpoint of each cell, now find it on the top and - # bottom edges of each cell - z_edges = zeros((size(z,0)+1, size(z,1), size(z,2))) - z_edges[1:-1,:,:] = 0.5*(z[0:-1,:,:] + z[1:,:,:]) - # At surface, z = 0; at bottom, set z to be the same as the midpoint of - # the deepest cell - z_edges[0,:,:] = z[0,:,:] - # Now find dz - dz = z_edges[1:,:,:] - z_edges[0:-1,:,:] - dV = dx*dy*dz - - # Read time data and convert from seconds to years - id = Dataset(file_path, 'r') - time = id.variables['ocean_time'][:]/(365*24*60*60) - # Read reference density - rho0 = id.variables['rho0'][:] - - ohc = zeros(size(time)) - for t in range(size(time)): - # Read temp and rho at this timestep - temp = id.variables['temp'][t,:,:,:] + celsius2kelvin - rho = id.variables['rho'][t,:,:,:] + rho0 - # Integrate temp*rho*cp over volume to get OHC - ohc[t] = sum(temp*rho*cp*dV) - id.close() - - # Plot results - clf() - plot(time, ohc) - xlabel('Years') - ylabel('Southern Ocean Heat Content (J)') - show() - - -# Command-line interface -if __name__ == "__main__": - - file_path = raw_input("Path to ocean averages file: ") - grid_path = raw_input("Path to grid file: ") - plot_ohc(file_path, grid_path) - - - - - - diff --git a/plot_tke.py b/plot_tke.py deleted file mode 100644 index 563e272..0000000 --- a/plot_tke.py +++ /dev/null @@ -1,162 +0,0 @@ -from netCDF4 import Dataset -from numpy import * -from matplotlib.pyplot import * -from calc_z import * - -# Calculate and plot the total kinetic energy at each timestep of the given -# ocean averages file. -def plot_tke (file_path, grid_path): - - # Grid parameters - theta_s = 0.9 - theta_b = 4.0 - hc = 40 - N = 31 - # Radius of the Earth in m - r = 6.371e6 - # Degrees to radians conversion factor - deg2rad = pi/180.0 - - # Read grid variables - id = Dataset(grid_path, 'r') - h = id.variables['h'][:,:] - zice = id.variables['zice'][:,:] - # Interpolate h and zice to u and v grids - h_u = 0.5*(h[:,0:-1] + h[:,1:]) - h_v = 0.5*(h[0:-1,:] + h[1:,:]) - zice_u = 0.5*(zice[:,0:-1] + zice[:,1:]) - zice_v = 0.5*(zice[0:-1,:] + zice[1:,:]) - lon_u = id.variables['lon_u'][:,:] - lat_u = id.variables['lat_u'][:,:] - mask_u = id.variables['mask_u'][:,:] - lon_v = id.variables['lon_v'][:,:] - lat_v = id.variables['lat_v'][:,:] - mask_v = id.variables['mask_v'][:,:] - id.close() - - # Mask lat and lon at land points - lon_u = ma.masked_where(mask_u==0, lon_u) - lat_u = ma.masked_where(mask_u==0, lat_u) - lon_v = ma.masked_where(mask_v==0, lon_v) - lat_v = ma.masked_where(mask_v==0, lat_v) - # Save dimensions - num_lat_u = size(lon_u, 0) - num_lon_u = size(lon_u, 1) - num_lat_v = size(lon_v, 0) - num_lon_v = size(lon_v, 1) - - # Add or subtract 360 from longitude values which wrap around - # so that longitude increases monotonically from west to east - i_u = tile(arange(1, num_lon_u+1), (num_lat_u, 1)) - index1_u = nonzero((i_u > 1200)*(lon_u < 100)) - lon_u[index1_u] = lon_u[index1_u] + 360 - index2_u = nonzero((i_u < 200)*(lon_u > 300)) - lon_u[index2_u] = lon_u[index2_u] - 360 - # Repeat for v grid - i_v = tile(arange(1, num_lon_v+1), (num_lat_v, 1)) - index1_v = nonzero((i_v > 1200)*(lon_v < 100)) - lon_v[index1_v] = lon_v[index1_v] + 360 - index2_v = nonzero((i_v < 200)*(lon_v > 300)) - lon_v[index2_v] = lon_v[index2_v] - 360 - - # Interpolate to get longitude at the edges of each cell - w_bdry_u = (lon_u[:,0] + lon_u[:,-1] - 360)/2 - middle_lon_u = (lon_u[:,0:-1] + lon_u[:,1:])/2 - e_bdry_u = (lon_u[:,0] + 360 + lon_u[:,-1])/2 - lon_u_edges = ma.concatenate((w_bdry_u[:,None], middle_lon_u, e_bdry_u[:,None]), axis=1) - # Subtract to get the change in longitude over each cell - dlon_u = abs(lon_u_edges[:,1:] - lon_u_edges[:,0:-1]) - # Repeat for v grid - w_bdry_v = (lon_v[:,0] + lon_v[:,-1] - 360)/2 - middle_lon_v = (lon_v[:,0:-1] + lon_v[:,1:])/2 - e_bdry_v = (lon_v[:,0] + 360 + lon_v[:,-1])/2 - lon_v_edges = ma.concatenate((w_bdry_v[:,None], middle_lon_v, e_bdry_v[:,None]), axis=1) - dlon_v = abs(lon_v_edges[:,1:] - lon_v_edges[:,0:-1]) - - # Similarly for latitude - s_bdry_u = lat_u[0,:] - middle_lat_u = (lat_u[0:-1,:] + lat_u[1:,:])/2 - n_bdry_u = lat_u[-1,:]*0 - 50 - lat_u_edges = ma.concatenate((s_bdry_u[None,:], middle_lat_u, n_bdry_u[None,:])) - dlat_u = lat_u_edges[1:,:] - lat_u_edges[0:-1,:] - # Repeat for v grid - s_bdry_v = lat_v[0,:] - middle_lat_v = (lat_v[0:-1,:] + lat_v[1:,:])/2 - n_bdry_v = lat_v[-1,:]*0 - 50 - lat_v_edges = ma.concatenate((s_bdry_v[None,:], middle_lat_v, n_bdry_v[None,:])) - dlat_v = lat_v_edges[1:,:] - lat_v_edges[0:-1,:] - - # Convert from spherical to Cartesian coordinates - # dy = r*dlat where dlat is converted to radians - dy_u_2d = r*dlat_u*pi/180.0 - dy_v_2d = r*dlat_v*pi/180.0 - # Copy into a 3D array, same at each depth level - dy_u = tile(dy_u_2d, (N,1,1)) - dy_v = tile(dy_v_2d, (N,1,1)) - # dx = r*cos(lat)*dlon where lat and dlon are converted to radians - dx_u_2d = r*cos(pi*lat_u/180.0)*dlon_u*pi/180.0 - dx_v_2d = r*cos(pi*lat_v/180.0)*dlon_v*pi/180.0 - dx_u = tile(dx_u_2d, (N,1,1)) - dx_v = tile(dx_v_2d, (N,1,1)) - - # Get a 3D array of z-coordinates on u and v grids - # sc_r and Cs_r are unused in this script - z_u, sc_r, Cs_r = calc_z(h_u, zice_u, lon_u, lat_u, theta_s, theta_b, hc, N) - z_v, sc_r, Cs_r = calc_z(h_v, zice_v, lon_v, lat_v, theta_s, theta_b, hc, N) - # We have z at the midpoint of each cell, now find it on the top and - # bottom edges of each cell - z_u_edges = zeros((size(z_u,0)+1, size(z_u,1), size(z_u,2))) - z_u_edges[1:-1,:,:] = 0.5*(z_u[0:-1,:,:] + z_u[1:,:,:]) - # At surface, z = 0; at bottom, set z to be the same as the midpoint of - # the deepest cell - z_u_edges[0,:,:] = z_u[0,:,:] - # Now find dz - dz_u = z_u_edges[1:,:,:] - z_u_edges[0:-1,:,:] - # Repeat on the v grid - z_v_edges = zeros((size(z_v,0)+1, size(z_v,1), size(z_v,2))) - z_v_edges[1:-1,:,:] = 0.5*(z_v[0:-1,:,:] + z_v[1:,:,:]) - z_v_edges[0,:,:] = z_v[0,:,:] - dz_v = z_v_edges[1:,:,:] - z_v_edges[0:-1,:,:] - # Calculate dV for each grid - dV_u = dx_u*dy_u*dz_u - dV_v = dx_v*dy_v*dz_v - - # Read time data and convert from seconds to years - id = Dataset(file_path, 'r') - time = id.variables['ocean_time'][:]/(365*24*60*60) - # Read reference density - rho0 = id.variables['rho0'][:] - - avgke = zeros(size(time)) - for t in range(size(time)): - # Read u, v, and rho at this timestep - u = id.variables['u'][t,:,:,:] - v = id.variables['v'][t,:,:,:] - rho = id.variables['rho'][t,:,:,:] + rho0 - # Interpolate rho onto u and v grids - rho_u = 0.5*(rho[:,:,0:-1] + rho[:,:,1:]) - rho_v = 0.5*(rho[:,0:-1,:] + rho[:,1:,:]) - # Integrate 0.5*rho*vel^2 over volume on each grid, add for TKE - avgke[t] = sum(0.5*rho_u*u**2*dV_u) + sum(0.5*rho_v*v**2*dV_v) - id.close() - - # Plot results - clf() - plot(time, avgke) - xlabel('Years') - ylabel('Southern Ocean Total Kinetic Energy (J)') - show() - - -# Command-line interface -if __name__ == "__main__": - - file_path = raw_input("Path to ocean averages file: ") - grid_path = raw_input("Path to grid file: ") - plot_tke(file_path, grid_path) - - - - - - diff --git a/plot_totalsalt.py b/plot_totalsalt.py deleted file mode 100644 index 36fe6fc..0000000 --- a/plot_totalsalt.py +++ /dev/null @@ -1,115 +0,0 @@ -from netCDF4 import Dataset -from numpy import * -from matplotlib.pyplot import * -from calc_z import * - -# Calculate and plot the total salt content at each timestep of the given ocean -# averages file. -def plot_totalsalt (file_path, grid_path): - - # Grid parameters - theta_s = 0.9 - theta_b = 4.0 - hc = 40 - N = 31 - # Radius of the Earth in m - r = 6.371e6 - # Degrees to radians conversion factor - deg2rad = pi/180.0 - - # Read grid variables - id = Dataset(grid_path, 'r') - h = id.variables['h'][:,:] - zice = id.variables['zice'][:,:] - lon = id.variables['lon_rho'][:,:] - lat = id.variables['lat_rho'][:,:] - mask = id.variables['mask_rho'][:,:] - id.close() - - # Mask lat and lon at land points - lon = ma.masked_where(mask==0, lon) - lat = ma.masked_where(mask==0, lat) - # Save dimensions - num_lat = size(lon, 0) - num_lon = size(lon, 1) - - # Add or subtract 360 from longitude values which wrap around - # so that longitude increases monotonically from west to east - i = tile(arange(1, num_lon+1), (num_lat, 1)) - index1 = nonzero((i > 1200)*(lon < 100)) - lon[index1] = lon[index1] + 360 - index2 = nonzero((i < 200)*(lon > 300)) - lon[index2] = lon[index2] - 360 - - # Interpolate to get longitude at the edges of each cell - w_bdry = (lon[:,0] + lon[:,-1] - 360)/2 - middle_lon = (lon[:,0:-1] + lon[:,1:])/2 - e_bdry = (lon[:,0] + 360 + lon[:,-1])/2 - lon_edges = ma.concatenate((w_bdry[:,None], middle_lon, e_bdry[:,None]), axis=1) - # Subtract to get the change in longitude over each cell - dlon = abs(lon_edges[:,1:] - lon_edges[:,0:-1]) - - # Similarly for latitude - s_bdry = lat[0,:] - middle_lat = (lat[0:-1,:] + lat[1:,:])/2 - n_bdry = lat[-1,:]*0 - 50 - lat_edges = ma.concatenate((s_bdry[None,:], middle_lat, n_bdry[None,:])) - dlat = lat_edges[1:,:] - lat_edges[0:-1,:] - - # Convert from spherical to Cartesian coordinates - # dy = r*dlat where dlat is converted to radians - dy_2d = r*dlat*pi/180.0 - # Copy into a 3D array, same at each depth level - dy = tile(dy_2d, (N,1,1)) - # dx = r*cos(lat)*dlon where lat and dlon are converted to radians - dx_2d = r*cos(pi*lat/180.0)*dlon*pi/180.0 - dx = tile(dx_2d, (N,1,1)) - - # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script - z, sc_r, Cs_r = calc_z(h, zice, lon, lat, theta_s, theta_b, hc, N) - # We have z at the midpoint of each cell, now find it on the top and - # bottom edges of each cell - z_edges = zeros((size(z,0)+1, size(z,1), size(z,2))) - z_edges[1:-1,:,:] = 0.5*(z[0:-1,:,:] + z[1:,:,:]) - # At surface, z = 0; at bottom, set z to be the same as the midpoint of - # the deepest cell - z_edges[0,:,:] = z[0,:,:] - # Now find dz - dz = z_edges[1:,:,:] - z_edges[0:-1,:,:] - dV = dx*dy*dz - - # Read time data and convert from seconds to years - id = Dataset(file_path, 'r') - time = id.variables['ocean_time'][:]/(365*24*60*60) - # Read reference density - rho0 = id.variables['rho0'][:] - - totalsalt = zeros(size(time)) - for t in range(size(time)): - # Read salt and rho at this timestep - salt = id.variables['salt'][t,:,:,:] - rho = id.variables['rho'][t,:,:,:] + rho0 - # Integrate 1e-3*salt*rho over volume to get total mass of salt - totalsalt[t] = sum(1e-3*salt*rho*dV) - id.close() - - # Plot results - clf() - plot(time, totalsalt) - xlabel('Years') - ylabel('Southern Ocean Salt Content (kg)') - show() - - -# Command-line interface -if __name__ == "__main__": - - file_path = raw_input("Path to ocean averages file: ") - grid_path = raw_input("Path to grid file: ") - plot_totalsalt(file_path, grid_path) - - - - - - diff --git a/unesco.py b/unesco.py new file mode 100644 index 0000000..50da7fa --- /dev/null +++ b/unesco.py @@ -0,0 +1,74 @@ +from numpy import * + +# Compute the UNESCO equation of state for seawater density +# given temp = temperature (C), salt = salinity (psu), and press = pressure +# (bar - note that p = depth (in m) / 10 is acceptable). Fully vectorized. +# See http://unesdoc.unesco.org/images/0004/000461/046148eb.pdf; all +# intermediate variables are named the same as in this document. + +# Input: +# temp = array of any dimension, containing temperature values (C) +# salt = array of same dimension as temp, containing salinity values (psu) +# press = array of same dimension as temp, containing pressure values (bar) +# for which 0.1*depth (in m) is close enough +# Output: rho = array of same dimension as temp, containing density values +# (kg/m^3) +def unesco (temp, salt, press): + + # Set constants + b0 = 8.24493e-1 + b1 = -4.0899e-3 + b2 = 7.6438e-5 + b3 = -8.2467e-7 + b4 = 5.3875e-9 + c0 = -5.72466e-3 + c1 = 1.0227e-4 + c2 = -1.6546e-6 + d0 = 4.8314e-4 + a0 = 999.842594 + a1 = 6.793952e-2 + a2 = -9.095290e-3 + a3 = 1.001685e-4 + a4 = -1.120083e-6 + a5 = 6.536332e-9 + f0 = 54.6746 + f1 = -0.603459 + f2 = 1.09987e-2 + f3 = -6.1670e-5 + g0 = 7.944e-2 + g1 = 1.6483e-2 + g2 = -5.3009e-4 + i0 = 2.2838e-3 + i1 = -1.0981e-5 + i2 = -1.6078e-6 + j0 = 1.91075e-4 + m0 = -9.9348e-7 + m1 = 2.0816e-8 + m2 = 9.1697e-10 + e0 = 19652.21 + e1 = 148.4206 + e2 = -2.327105 + e3 = 1.360477e-2 + e4 = -5.155288e-5 + h0 = 3.239908 + h1 = 1.43713e-3 + h2 = 1.16092e-4 + h3 = -5.77905e-7 + k0 = 8.50935e-5 + k1 = -6.12293e-6 + k2 = 5.2787e-8 + + rho_0 = a0 + a1*temp + a2*temp**2 + a3*temp**3 + a4*temp**4 + a5*temp**5 + \ + (b0 + b1*temp + b2*temp**2 + b3*temp**3 + b4*temp**4)*salt + \ + (c0 + c1*temp + c2*temp**2)*salt**1.5 + d0*salt**2 + A = h0 + h1*temp + h2*temp**2 + h3*temp**3 + \ + (i0 + i1*temp + i2*temp**2)*salt + j0*salt**1.5 + B = k0 + k1*temp + k2*temp**2 + (m0 + m1*temp + m2*temp**2)*salt + K = e0 + e1*temp + e2*temp**2 + e3*temp**3 + e4*temp**4 + \ + (f0 + f1*temp + f2*temp**2 + f3*temp**3)*salt + \ + (g0 + g1*temp + g2*temp**2)*salt**1.5 + A*press + B*press**2 + rho = rho_0/(1 - press/K) + + return rho + +