# calc_area_xyz # calc_rho_from_theta # calc_sal_from_rho_bis # calc_theta_from_temp # gsw_g # volume_census_levitus # volume_etopo # volume_from_sth def calc_area_xyz(grd,r_earth,calc_old=0): # def calc_area_xyz(grd,r_earth) ######################################################### # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from get_area_latitude.m of Remi Tailleux (09/2018) # last modified : 20.09.2018 # Content/Description #################################################################### # Usage : from mycalc_ocean import calc_area_xyz # Input : grd, r_earth # grd : Input grid [python class] # r_earth : earth radius [m] # Output : area_lat # area_xyz : hor. area for 3d grid [m**2] # ######################################################################################## # import modules import numpy as np # convert deg to rad for sin/cos deg_to_rad = np.pi/180.0 lat_radians = grd.lat*deg_to_rad if calc_old == 1: print ' *** old version of area calculation ***' # compute length of longitudinal and latitudinal elements print ' *** area calculation not valid for non-equidistant input lon-grid ***' print ' *** area calculation not valid for non-equidistant input lat-grid ***' dx = r_earth*deg_to_rad*(grd.Nlon/360.0) dy = r_earth*deg_to_rad*abs(grd.lat[1]-grd.lat[0]) # calculate area area_xyz = np.zeros([grd.Nlon,grd.Nlat,grd.Nz]) for i_lon in range(0,grd.Nlon): for i_z in range(0,grd.Nz): area_xyz[i_lon,:,i_z] = (dx*dy)*np.cos(lat_radians) else: lat_d, lon_d, z_d = np.meshgrid(grd.lat,grd.lon,grd.z) if grd.symlon == 1: delta_lon = 360.0 index_d = np.append(np.arange(1,grd.Nlon),0) lon_e = 0.5 * (lon_d + lon_d[index_d,:,:]) lon_e[-1,:,:] = lon_e[-1,:,:] + 0.5 * delta_lon index_d = np.append(grd.Nlon-1,np.arange(0,grd.Nlon-1)) lon_w = 0.5 * (lon_d + lon_d[index_d,:,:]) lon_w[0,:,:] = lon_w[0,:,:] - 0.5 * delta_lon else: lon_d = np.append(2*grd.lon[0] - grd.lon[1],grd.lon) lon_d = np.append(lon_d,2*grd.lon[-1] - grd.lon[-2]) dummy1, lon_d, dummy2 = np.meshgrid(grd.lat,lon_d,grd.z) index_d = np.arange(1,grd.Nlon+1) lon_e = 0.5 * (lon_d[index_d,:,:] + delta_lon + lon_d[index_d+1,:,:]) lon_w = 0.5 * (lon_d[index_d,:,:] + delta_lon + lon_d[index_d-1,:,:]) lat_d = np.append(max(2*grd.lat[0] - grd.lat[1],-180.0-grd.lat[0]),grd.lat) lat_d = np.append(lat_d,min(2*grd.lat[-1] - grd.lat[-2],180.0-grd.lat[-1])) lat_d, dummy1, dummy2 = np.meshgrid(lat_d,grd.lon,grd.z) index_d = np.arange(1,grd.Nlat+1) lat_n = 0.5 * (lat_d[:,index_d,:] + lat_d[:,index_d+1,:]) lat_s = 0.5 * (lat_d[:,index_d,:] + lat_d[:,index_d-1,:]) r_z = r_earth - abs(grd.z) #area_xyz = (deg_to_rad*r_z)**2*(lon_e-lon_w)*(lat_n-lat_s)*np.cos(lat_d[:,index_d,:]*deg_to_rad) area_xyz = deg_to_rad*r_z**2*(lon_e-lon_w)*(np.sin(lat_n*deg_to_rad)-np.sin(lat_s*deg_to_rad)) del delta_lon, lon_d, lat_d, z_d, dummy1, dummy2, index_d, r_z # return horizontal area for input grid return area_xyz def calc_rho_from_theta(s, th, p, reference='jackettetal2004'): # calc_rho_from_theta(s,th,p,reference='jackettetal2004') ################################ # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from rho_from_theta.m of Remi Tailleux (09/2018) # last modified : 19.09.2018 # Gabriel Wolf : added polynomials for 'jackettandmcdougall1995' STILL MISSING # Content/Description #################################################################### # from mycalc_ocean import calc_rho_from_theta # Density from potential temperature, by the use of polynomials from reference # reference default is Jackett, McDougall, Feistel, Wright and Griffies 2004 # -> https://doi.org/10.1175/JTECH1946.1 ; equation (A2) and table A2 therein # usage: from mycalc_ocean import calc_rho_from_theta # calc_rho_from_theta(s,th,p,reference=) # defined by publication, e.g. author1year if there is only one author, # author1andauthor2year for two authors or author1etalyear for more than 2 authors # input: s, th, p, reference # s : salinity [psu] # th : potential temperature [deg C] # p : gauge pressure [bars] # (absolute pressure - 10.1325) # reference : used polynomials [string] # output: rho : density [kg m**-3] # check: calc_rho_from_theta(20,20,1000) = 1017.7288680196421 # ######################################################################################## # load necessary modules import numpy as np # simplify calculation by predefining some terms th2 = th*th sqrts = np.sqrt(s); pth = p*th; # ######################################################################################## # define polynomial # ######################################################################################## if reference == 'jackettetal2004': # https://doi.org/10.1175/JTECH1946.1 ; equation (A2) and table A2 therein # numerator of equ.(A2) eqA2_num = 9.9984085444849347e02 + \ th * (7.3471625860981584e00 + \ th * (-5.3211231792841769e-02 + \ th * 3.6492439109814549e-04)) + \ s * (2.5880571023991390e+00 + \ th * -6.7168282786692355e-03 + \ s * 1.9203202055760151e-03) + \ p * (1.1798263740430364e-02 + \ th2 * 9.8920219266399117e-08 + \ s * 4.6996642771754730e-06 + \ p * ( -2.5862187075154352e-08 + \ th2 * (-3.2921414007960662e-12))) # denominator of equ.(A2) eqA2_den = 1e00 + \ th * (7.2815210113327091e-03 + \ th * (-4.4787265461983921e-05 + \ th * ( 3.3851002965802430e-07 + \ th * 1.3651202389758572e-10))) + \ s * (1.7632126669040377e-03 + \ th * ( -8.8066583251206474e-06 + \ th2 * (-1.8832689434804897e-10)) + \ sqrts * (5.7463776745432097e-06 + \ th2 * 1.4716275472242334e-09)) + \ p * (6.7103246285651894e-06 + \ pth * (th2 * (-2.4461698007024582e-17) + \ p * (-9.1534417604289062e-18))); # equation (B5) for the equation of state rho = eqA2_num / eqA2_den # Deallocation del eqA2_num, eqA2_den elif reference == 'jackettandmcdougall1995': raise ValueError('No valid input for reference in the function calc_rho_from_theta.') else: raise ValueError('No valid input for reference in the function calc_rho_from_theta.') # ######################################################################################## # Deallocation of Variables # ######################################################################################## del th2, sqrts, pth # return density return rho def calc_sal_from_rho(rho,th,p,reference='jackettetal2004'): # calc_sal_from_rho(rho,th,p,reference='jackettetal2004') ################################ # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from sal_from_rho_bis.m of Remi Tailleux (DRJ on 10/12/2003) # last modified : 24.09.2018 # Gabriel Wolf : added polynomials for 'jackettandmcdougall1995' STILL MISSING print ' *** PLEASE CHECK THIS: conversion method different to R. Tailleux - why? ***' # Content/Description #################################################################### # Salinity from density and pot. temperature, by the use of polynomials from reference # reference default is Jackett, McDougall, Feistel, Wright and Griffies 2004 # -> https://doi.org/10.1175/JTECH1946.1 ; equation (A2) and table A2 therein # usage: from mycalc_ocean import calc_sal_from_rho # calc_sal_from_rho(rho,th,p,reference=) # defined by publication, e.g. author1year if there is only one author, # author1andauthor2year for two authors or author1etalyear for more than 2 authors # input: rho, th, p, reference # rho : denisty [kg m**-3] # th : potential temperature [deg C] # p : gauge pressure [bars] # (absolute pressure - 10.1325) # reference : used polynomials [string] # output: sal : salinity [psu] # check: calc_rho_from_theta(20,20,1000) = 1017.7288680196421 # ######################################################################################## # load necessary modules %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% import numpy as np # simplify calculation by predefining some terms %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% th2 = th*th pth = p*th; # define polynomial and calculate salinity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if reference == 'jackettetal2004': # https://doi.org/10.1175/JTECH1946.1 ; equation (A2) and table A2 therein # Define numerator of equ.(A2) ------------------------------------------------------- c_0 = 9.9984085444849347e02 c_th1 = 7.3471625860981584e00 c_th2 = -5.3211231792841769e-02 c_th3 = 3.6492439109814549e-04 c_s1 = 2.5880571023991390e+00 c_sth1 = -6.7168282786692355e-03 c_s2 = 1.9203202055760151e-03 c_p1 = 1.1798263740430364e-02 c_pth2 = 9.8920219266399117e-08 c_ps1 = 4.6996642771754730e-06 c_p2 = -2.5862187075154352e-08 c_p2th2 = -3.2921414007960662e-12 # Full equation (numerator): eqA2_num # eqA2_num = c_0 + th * (c_th1 + th * (c_th2 + th * c_th3)) +\ # \ s * ( c_s1 + th * c_sth1 + s * c_s2) +\ # p * ( c_p1 + th2 * c_pth2 + s * c_ps1 + p * ( c_p2 + c_p2th2)) # eqA2_num = eqA2_num_s0 + s * eqA2_num_s1 + s**2 * c_s2 eqA2_num_s0 = c_0 + th * (c_th1 + th * (c_th2 + th * c_th3)) +\ p * ( c_p1 + th2 * c_pth2 + p * ( c_p2 + c_p2th2)) eqA2_num_s1 = c_s1 + th * c_sth1 + c_ps1 * p # Define denominator of equ.(A2) ----------------------------------------------------- d_0 = 1e00 d_th1 = 7.2815210113327091e-03 d_th2 = -4.4787265461983921e-05 d_th3 = 3.3851002965802430e-07 d_th4 = 1.3651202389758572e-10 d_s1 = 1.7632126669040377e-03 d_sth1 = -8.8066583251206474e-06 d_sth3 = -1.8832689434804897e-10 d_s3d2 = 5.7463776745432097e-06 d_s3d2th2 = 1.4716275472242334e-09 d_p1 = 6.7103246285651894e-06 d_p2th3 = -2.4461698007024582e-17 d_p3th1 = -9.1534417604289062e-18 # Full equation of denominator (A2): eqA2_den # eqA2_den = d_0 + th * (d_th1 + th * (d_th2 + th * (d_th3 + th * d_th4))) +\ # s * (d_s1 + th * (d_sth1 + th2 * d_sth3) + sqrts * (d_s3d2 + th2 * d_s3d2th2)) +\ # p * (d_p1 + pth * (th2 * d_p2th3 + p * d_p3th1)) # eqA2_den = eqA2_den_s0 + s * (eqA2_den_s1 + sqrts * eqA2_den_s2) eqA2_den_s0 = d_0 + th * (d_th1 + th * (d_th2 + th * (d_th3 + th * d_th4))) +\ p * (d_p1 + pth * (th2 * d_p2th3 + p * d_p3th1)) eqA2_den_s1 = d_s1 + th * (d_sth1 + th2 * d_sth3) eqA2_den_s2 = d_s3d2 + th2 * d_s3d2th2 # Apply Newton's method to find value for salinity ----------------------------------- # d(f(s_i))/ds * (s_(i+1) - s_i) + f(s_i) = 0 # from this equation: s_(i+1) = s_i - f(s_i)/(d(f(s_i))/ds) s_old = 20.0*np.ones(len(th)) new_err = 1e-02 i_iter = 0 while i_iter < 100: i_iter = i_iter + 1 sqrts = np.sqrt(s_old) eqA2_num = eqA2_num_s0 + s_old * eqA2_num_s1 + s_old**2 * c_s2 eqA2_den = eqA2_den_s0 + s_old * (eqA2_den_s1 + sqrts * eqA2_den_s2) der_eqA2_num = eqA2_num_s1 + 2.0 * c_s2 * s_old der_eqA2_den = eqA2_den_s1 + 3.0/2.0 * sqrts * eqA2_den_s2 s_new = s_old - eqA2_num * eqA2_den / \ (der_eqA2_num * eqA2_den - eqA2_num * der_eqA2_den) i_err = np.max(abs(s_new-s_old)) s_old = s_new print ' Newton iteration = ', i_iter if i_err < new_err: print ' Newton method converged after ' + str(i_iter) + ' iterations' break if i_iter == 100: print ' *** Newton method did not converge for error limit ' + \ str(new_err) + ' ***' print ' *** Final error = ' + str(i_err) + ' ***' # Deallocation ----------------------------------------------------------------------- del sqrts, s_old, s_err, new_err, i_err, i_iter, eqA2_num, eqA2_den,\ der_eqA2_num, der_eqA2_den elif reference == 'jackettandmcdougall1995': raise ValueError('No valid input for reference in the function calc_sal_from_rho.') else: raise ValueError('No valid input for reference in the function calc_sal_from_rho.') # Deallocation and return variables ------------------------------------------------------ # Deallocation del pth, th2 # Return density return s_new def calc_theta_from_temp(s,temp,p,pr): # calc_theta_from_temp(s,temp,p,pr) ###################################################### # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from gsw_ptmp.m of Remi Tailleux # last modified : 14.09.2018 # Content/Description #################################################################### # from mycalc_ocean import calc_theta_from_temp # Potential temperature of seawater from specific entropy # Input : s, temp, p, pr # s : absolute salinity [g/Kg] # temp : temperature [deg C] # p : sea (gauge) pressure [dbar] # pr : reference (gauge) pressure [dbar] # Output : theta # theta : potential temperature [deg C] # ######################################################################################## # Load modules %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% import numpy as np # Error checking %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if s.shape != temp.shape: raise ValueError('**** Data shape of salinity and temperature do not agree ****') if s.shape != p.shape: raise ValueError('**** Data shape of salinity and pressure do not agree ****') if s.shape != pr.shape: raise ValueError('**** Data shape of salinity and reference pressure do not agree ****') # Allocation and Initialization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% s1 = s*35.0/35.16504 t1 = np.copy(temp) p1 = np.copy(p) pr1 = pr # First estimate for potential temperature %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # equ (A3) and (A4) in McDougall et al. (2003) with table A1 in Jackett et al. (2004) # McDougall et al. (2003): https://doi.org/10.1175/1520-0426(2003)20<730:AACEAF>2.0.CO;2 # Jackett et al. (2004): https://doi.org/10.1175/JTECH1946.1 # For further information see Appendix A in Jackett et al. (2004) p_polynom = ( 8.65483913395442e-6 - \ s1*1.41636299744881e-6 - \ (p1+pr1) * 7.38286467135737e-9 + \ t1 *(-8.38241357039698e-6 + \ s1 * 2.83933368585534e-8 + \ t1 * 1.77803965218656e-8 + \ (p1+pr1) * 1.71155619208233e-10)) theta = t1 + (p1-pr1)* p_polynom de_dt = 13.6 # Newton method to determine theta %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% n_iter = 3 # default i_iter = 0 # initialization for count-variable print ' *** Why using 3 iterations for Newton method in calc_theta_from_temp? ***' # NOTE: n_iter = 1 gives theta with a maximum error of ... # n_iter = 2 gives theta with a maximum error of ... n0 = 0 n2 = 2 print ' *** Where are the equations in this while loop coming from? ***' while i_iter <= n_iter: i_iter = i_iter + 1 theta_old = np.copy(theta) dentropy = gsw_g(0,1,0,s,theta,pr) - gsw_g(0,1,0,s,temp,p) # gsw_g(0,1,0,s,t,p) : entropy theta = theta - dentropy/de_dt theta = 0.5 * (theta + theta_old) dentropy_dt = -gsw_g(n0,n2,n0,s,theta,pr) theta = theta_old - dentropy/dentropy_dt # Define/return output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% return theta def gsw_g(Ns,Nt,Np,s,temp,p): # gsw_g(Ns,Nt,Np,s,temp,p) ############################################################### # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from gsw_g.m of Remi Tailleux # last modified : 14.09.2018 # Content/Description #################################################################### # seawater specific Gibbs free energy and derivatives up to order 2 # Input : Ns, Nt, Np, s, temp, p # Ns : order of s derivative # Nt : order of t derivative # Np : order of p derivative # s : absolute salinity [g/Kg] # temp : temperature [deg C] # p : sea (gauge) pressure [dbar] # Output : theta # result : # ######################################################################################## # Import modules import numpy as np # Allocation/Initialization sfac = 0.0248826675584615 # sfac = 1/(40*(35.16504/35)) x2 = sfac*s x = np.sqrt(x2) y = temp*0.025e0 z = p*1e-4 # Choose correct calculation (depending on Ns, Nt and Np) if Ns==0 and Nt==0 and Np==0: g03 = 101.342743139674 + z*(100015.695367145 + \ z*(-2544.5765420363 + z*(284.517778446287 + \ z*(-33.3146754253611 + (4.20263108803084 - 0.546428511471039*z)*z)))) + \ y*(5.90578347909402 + z*(-270.983805184062 + \ z*(776.153611613101 + z*(-196.51255088122 + (28.9796526294175 - 2.13290083518327*z)*z))) + \ y*(-12357.785933039 + z*(1455.0364540468 + \ z*(-756.558385769359 + z*(273.479662323528 + z*(-55.5604063817218 + 4.34420671917197*z)))) + \ y*(736.741204151612 + z*(-672.50778314507 + \ z*(499.360390819152 + z*(-239.545330654412 + (48.8012518593872 - 1.66307106208905*z)*z))) + \ y*(-148.185936433658 + z*(397.968445406972 + \ z*(-301.815380621876 + (152.196371733841 - 26.3748377232802*z)*z)) + \ y*(58.0259125842571 + z*(-194.618310617595 + \ z*(120.520654902025 + z*(-55.2723052340152 + 6.48190668077221*z))) + \ y*(-18.9843846514172 + y*(3.05081646487967 - 9.63108119393062*z) + \ z*(63.5113936641785 + z*(-22.2897317140459 + 8.17060541818112*z)))))))) g08 = x2*(1416.27648484197 + z*(-3310.49154044839 + \ z*(384.794152978599 + z*(-96.5324320107458 + (15.8408172766824 - 2.62480156590992*z)*z))) + \ x*(-2432.14662381794 + x*(2025.80115603697 + \ y*(543.835333000098 + y*(-68.5572509204491 + \ y*(49.3667694856254 + y*(-17.1397577419788 + 2.49697009569508*y))) - 22.6683558512829*z) + \ x*(-1091.66841042967 - 196.028306689776*y + \ x*(374.60123787784 - 48.5891069025409*x + 36.7571622995805*y) + 36.0284195611086*z) + \ z*(-54.7919133532887 + (-4.08193978912261 - 30.1755111971161*z)*z)) + \ z*(199.459603073901 + z*(-52.2940909281335 + (68.0444942726459 - 3.41251932441282*z)*z)) + \ y*(-493.407510141682 + z*(-175.292041186547 + (83.1923927801819 - 29.483064349429*z)*z) + \ y*(-43.0664675978042 + z*(383.058066002476 + z*(-54.1917262517112 + 25.6398487389914*z)) + \ y*(-10.0227370861875 - 460.319931801257*z + y*(0.875600661808945 + 234.565187611355*z))))) + \ y*(168.072408311545 + z*(729.116529735046 + \ z*(-343.956902961561 + z*(124.687671116248 + z*(-31.656964386073 + 7.04658803315449*z)))) + \ y*(880.031352997204 + y*(-225.267649263401 + \ y*(91.4260447751259 + y*(-21.6603240875311 + 2.13016970847183*y) + \ z*(-297.728741987187 + (74.726141138756 - 36.4872919001588*z)*z)) + \ z*(694.244814133268 + z*(-204.889641964903 + (113.561697840594 - 11.1282734326413*z)*z))) + \ z*(-860.764303783977 + z*(337.409530269367 + \ z*(-178.314556207638 + (44.2040358308 - 7.92001547211682*z)*z)))))) inds = x>0 if sum(inds)>0: g08[inds] = g08[inds] + x2[inds]*(5812.81456626732 + 851.226734946706*y[inds])*np.log(x[inds]) result = g03 + g08 elif Ns==1 and Nt==0 and Np==0: g08 = 8645.36753595126 + z*(-6620.98308089678 + \ z*(769.588305957198 + z*(-193.0648640214916 + (31.6816345533648 - 5.24960313181984*z)*z))) + \ x*(-7296.43987145382 + x*(8103.20462414788 + \ y*(2175.341332000392 + y*(-274.2290036817964 + \ y*(197.4670779425016 + y*(-68.5590309679152 + 9.98788038278032*y))) - 90.6734234051316*z) + \ x*(-5458.34205214835 - 980.14153344888*y + \ x*(2247.60742726704 - 340.1237483177863*x + 220.542973797483*y) + 180.142097805543*z) + \ z*(-219.1676534131548 + (-16.32775915649044 - 120.7020447884644*z)*z)) + \ z*(598.378809221703 + z*(-156.8822727844005 + (204.1334828179377 - 10.23755797323846*z)*z)) + \ y*(-1480.222530425046 + z*(-525.876123559641 + (249.57717834054571 - 88.449193048287*z)*z) + \ y*(-129.1994027934126 + z*(1149.174198007428 + z*(-162.5751787551336 + 76.9195462169742*z)) + \ y*(-30.0682112585625 - 1380.9597954037708*z + y*(2.626801985426835 + 703.695562834065*z))))) + \ y*(1187.3715515697959 + z*(1458.233059470092 + \ z*(-687.913805923122 + z*(249.375342232496 + z*(-63.313928772146 + 14.09317606630898*z)))) + \ y*(1760.062705994408 + y*(-450.535298526802 + \ y*(182.8520895502518 + y*(-43.3206481750622 + 4.26033941694366*y) + \ z*(-595.457483974374 + (149.452282277512 - 72.9745838003176*z)*z)) + \ z*(1388.489628266536 + z*(-409.779283929806 + (227.123395681188 - 22.2565468652826*z)*z))) + \ z*(-1721.528607567954 + z*(674.819060538734 + \ z*(-356.629112415276 + (88.4080716616 - 15.84003094423364*z)*z))))) inds = x>0 if sum(inds)>0: g08[inds] = g08[inds] + (11625.62913253464 + 1702.453469893412*y[inds])*np.log(x[inds]) inds = x==0 if sum(inds)>0: g08[inds] = np.nan; result = 0.5e3*sfac*g08; elif Ns==0 and Nt==1 and Np==0: g03 = 5.90578347909402 + z*(-270.983805184062 + \ z*(776.153611613101 + z*(-196.51255088122 + (28.9796526294175 - 2.13290083518327*z)*z))) + \ y*(-24715.571866078 + z*(2910.0729080936 + \ z*(-1513.116771538718 + z*(546.959324647056 + z*(-111.1208127634436 + 8.68841343834394*z)))) + \ y*(2210.2236124548363 + z*(-2017.52334943521 + \ z*(1498.081172457456 + z*(-718.6359919632359 + (146.4037555781616 - 4.9892131862671505*z)*z))) + \ y*(-592.743745734632 + z*(1591.873781627888 + \ z*(-1207.261522487504 + (608.785486935364 - 105.4993508931208*z)*z)) + \ y*(290.12956292128547 + z*(-973.091553087975 + \ z*(602.603274510125 + z*(-276.361526170076 + 32.40953340386105*z))) + \ y*(-113.90630790850321 + y*(21.35571525415769 - 67.41756835751434*z) + \ z*(381.06836198507096 + z*(-133.7383902842754 + 49.023632509086724*z))))))) g08 = x2*(168.072408311545 + z*(729.116529735046 + \ z*(-343.956902961561 + z*(124.687671116248 + z*(-31.656964386073 + 7.04658803315449*z)))) + \ x*(-493.407510141682 + x*(543.835333000098 + x*(-196.028306689776 + 36.7571622995805*x) + \ y*(-137.1145018408982 + y*(148.10030845687618 + y*(-68.5590309679152 + 12.4848504784754*y))) - \ 22.6683558512829*z) + z*(-175.292041186547 + (83.1923927801819 - 29.483064349429*z)*z) + \ y*(-86.1329351956084 + z*(766.116132004952 + z*(-108.3834525034224 + 51.2796974779828*z)) + \ y*(-30.0682112585625 - 1380.9597954037708*z + y*(3.50240264723578 + 938.26075044542*z)))) + \ y*(1760.062705994408 + y*(-675.802947790203 + \ y*(365.7041791005036 + y*(-108.30162043765552 + 12.78101825083098*y) + \ z*(-1190.914967948748 + (298.904564555024 - 145.9491676006352*z)*z)) + \ z*(2082.7344423998043 + z*(-614.668925894709 + (340.685093521782 - 33.3848202979239*z)*z))) + \ z*(-1721.528607567954 + z*(674.819060538734 + \ z*(-356.629112415276 + (88.4080716616 - 15.84003094423364*z)*z))))) inds = x>0 if sum(inds)>0: g08[inds] = g08[inds] + 851.226734946706*x2[inds]*np.log(x[inds]) result = (g03 + g08)*0.025e0 elif Ns==0 and Nt==0 and Np==1: g03 = 100015.695367145 + z*(-5089.1530840726 + \ z*(853.5533353388611 + z*(-133.2587017014444 + (21.0131554401542 - 3.278571068826234*z)*z))) + \ y*(-270.983805184062 + z*(1552.307223226202 + \ z*(-589.53765264366 + (115.91861051767 - 10.664504175916349*z)*z)) + \ y*(1455.0364540468 + z*(-1513.116771538718 + \ z*(820.438986970584 + z*(-222.2416255268872 + 21.72103359585985*z))) + \ y*(-672.50778314507 + z*(998.720781638304 + \ z*(-718.6359919632359 + (195.2050074375488 - 8.31535531044525*z)*z)) + \ y*(397.968445406972 + z*(-603.630761243752 + (456.589115201523 - 105.4993508931208*z)*z) + \ y*(-194.618310617595 + y*(63.5113936641785 - 9.63108119393062*y + \ z*(-44.5794634280918 + 24.511816254543362*z)) + \ z*(241.04130980405 + z*(-165.8169157020456 + 25.92762672308884*z))))))) g08 = x2*(-3310.49154044839 + z*(769.588305957198 + \ z*(-289.5972960322374 + (63.3632691067296 - 13.1240078295496*z)*z)) + \ x*(199.459603073901 + x*(-54.7919133532887 + 36.0284195611086*x - 22.6683558512829*y + \ (-8.16387957824522 - 90.52653359134831*z)*z) + \ z*(-104.588181856267 + (204.1334828179377 - 13.65007729765128*z)*z) + \ y*(-175.292041186547 + (166.3847855603638 - 88.449193048287*z)*z + \ y*(383.058066002476 + y*(-460.319931801257 + 234.565187611355*y) + \ z*(-108.3834525034224 + 76.9195462169742*z)))) + \ y*(729.116529735046 + z*(-687.913805923122 + \ z*(374.063013348744 + z*(-126.627857544292 + 35.23294016577245*z))) + \ y*(-860.764303783977 + y*(694.244814133268 + \ y*(-297.728741987187 + (149.452282277512 - 109.46187570047641*z)*z) + \ z*(-409.779283929806 + (340.685093521782 - 44.5130937305652*z)*z)) + \ z*(674.819060538734 + z*(-534.943668622914 + (176.8161433232 - 39.600077360584095*z)*z))))) result = (g03 + g08)*1e-8 elif Ns==0 and Nt==2 and Np==0: g03 = -24715.571866078 + z*(2910.0729080936 + z* \ (-1513.116771538718 + z*(546.959324647056 + z*(-111.1208127634436 + 8.68841343834394*z)))) + \ y*(4420.4472249096725 + z*(-4035.04669887042 + \ z*(2996.162344914912 + z*(-1437.2719839264719 + (292.8075111563232 - 9.978426372534301*z)*z))) + \ y*(-1778.231237203896 + z*(4775.621344883664 + \ z*(-3621.784567462512 + (1826.356460806092 - 316.49805267936244*z)*z)) + \ y*(1160.5182516851419 + z*(-3892.3662123519 + \ z*(2410.4130980405 + z*(-1105.446104680304 + 129.6381336154442*z))) + \ y*(-569.531539542516 + y*(128.13429152494615 - 404.50541014508605*z) + \ z*(1905.341809925355 + z*(-668.691951421377 + 245.11816254543362*z)))))) g08 = x2*(1760.062705994408 + x*(-86.1329351956084 + \ x*(-137.1145018408982 + y*(296.20061691375236 + y*(-205.67709290374563 + 49.9394019139016*y))) + \ z*(766.116132004952 + z*(-108.3834525034224 + 51.2796974779828*z)) + \ y*(-60.136422517125 - 2761.9195908075417*z + y*(10.50720794170734 + 2814.78225133626*z))) + \ y*(-1351.605895580406 + y*(1097.1125373015109 + y*(-433.20648175062206 + 63.905091254154904*y) + \ z*(-3572.7449038462437 + (896.713693665072 - 437.84750280190565*z)*z)) + \ z*(4165.4688847996085 + z*(-1229.337851789418 + (681.370187043564 - 66.7696405958478*z)*z))) + \ z*(-1721.528607567954 + z*(674.819060538734 + \ z*(-356.629112415276 + (88.4080716616 - 15.84003094423364*z)*z)))) result = (g03 + g08)*0.000625 elif Ns==1 and Nt==0 and Np==1: g08 = -6620.98308089678 + z*(1539.176611914396 + \ z*(-579.1945920644748 + (126.7265382134592 - 26.2480156590992*z)*z)) + \ x*(598.378809221703 + x*(-219.1676534131548 + 180.142097805543*x - 90.6734234051316*y + \ (-32.65551831298088 - 362.10613436539325*z)*z) + \ z*(-313.764545568801 + (612.4004484538132 - 40.95023189295384*z)*z) + \ y*(-525.876123559641 + (499.15435668109143 - 265.347579144861*z)*z + \ y*(1149.174198007428 + y*(-1380.9597954037708 + 703.695562834065*y) + \ z*(-325.1503575102672 + 230.7586386509226*z)))) + \ y*(1458.233059470092 + z*(-1375.827611846244 + \ z*(748.126026697488 + z*(-253.255715088584 + 70.4658803315449*z))) + \ y*(-1721.528607567954 + y*(1388.489628266536 + \ y*(-595.457483974374 + (298.904564555024 - 218.92375140095282*z)*z) + \ z*(-819.558567859612 + (681.370187043564 - 89.0261874611304*z)*z)) + \ z*(1349.638121077468 + z*(-1069.887337245828 + (353.6322866464 - 79.20015472116819*z)*z)))) g08 = g08 result = g08*sfac*0.5e-8 elif Ns==0 and Nt==1 and Np==1: g03 = -270.983805184062 + z*(1552.307223226202 + z*(-589.53765264366 + (115.91861051767 - 10.664504175916349*z)*z)) + \ y*(2910.0729080936 + z*(-3026.233543077436 + \ z*(1640.877973941168 + z*(-444.4832510537744 + 43.4420671917197*z))) + \ y*(-2017.52334943521 + z*(2996.162344914912 + \ z*(-2155.907975889708 + (585.6150223126464 - 24.946065931335752*z)*z)) + \ y*(1591.873781627888 + z*(-2414.523044975008 + (1826.356460806092 - 421.9974035724832*z)*z) + \ y*(-973.091553087975 + z*(1205.20654902025 + z*(-829.084578510228 + 129.6381336154442*z)) + \ y*(381.06836198507096 - 67.41756835751434*y + z*(-267.4767805685508 + 147.07089752726017*z)))))) g08 = x2*(729.116529735046 + z*(-687.913805923122 + \ z*(374.063013348744 + z*(-126.627857544292 + 35.23294016577245*z))) + \ x*(-175.292041186547 - 22.6683558512829*x + (166.3847855603638 - 88.449193048287*z)*z + \ y*(766.116132004952 + y*(-1380.9597954037708 + 938.26075044542*y) + \ z*(-216.7669050068448 + 153.8390924339484*z))) + \ y*(-1721.528607567954 + y*(2082.7344423998043 + \ y*(-1190.914967948748 + (597.809129110048 - 437.84750280190565*z)*z) + \ z*(-1229.337851789418 + (1022.055280565346 - 133.5392811916956*z)*z)) + \ z*(1349.638121077468 + z*(-1069.887337245828 + (353.6322866464 - 79.20015472116819*z)*z)))) result = (g03 + g08)*2.5e-10 elif Ns==0 and Nt==0 and Np==2: g03 = -5089.1530840726 + z*(1707.1066706777221 + \ z*(-399.7761051043332 + (84.0526217606168 - 16.39285534413117*z)*z)) + \ y*(1552.307223226202 + z*(-1179.07530528732 + (347.75583155301 - 42.658016703665396*z)*z) + \ y*(-1513.116771538718 + z*(1640.877973941168 + z*(-666.7248765806615 + 86.8841343834394*z)) + \ y*(998.720781638304 + z*(-1437.2719839264719 + (585.6150223126464 - 33.261421241781*z)*z) + \ y*(-603.630761243752 + (913.178230403046 - 316.49805267936244*z)*z + \ y*(241.04130980405 + y*(-44.5794634280918 + 49.023632509086724*z) + \ z*(-331.6338314040912 + 77.78288016926652*z)))))) g08 = x2*(769.588305957198 + z*(-579.1945920644748 + (190.08980732018878 - 52.4960313181984*z)*z) + \ x*(-104.588181856267 + x*(-8.16387957824522 - 181.05306718269662*z) + \ (408.2669656358754 - 40.95023189295384*z)*z + \ y*(166.3847855603638 - 176.898386096574*z + y*(-108.3834525034224 + 153.8390924339484*z))) + \ y*(-687.913805923122 + z*(748.126026697488 + z*(-379.883572632876 + 140.9317606630898*z)) + \ y*(674.819060538734 + z*(-1069.887337245828 + (530.4484299696 - 158.40030944233638*z)*z) + \ y*(-409.779283929806 + y*(149.452282277512 - 218.92375140095282*z) + \ (681.370187043564 - 133.5392811916956*z)*z)))) result = (g03 + g08)*1e-16 return result def volume_census_levitus(th, s, grd, thbin, sbin): # volume_census_levitus(th,s,grd,thbin,sbin) ############################################# # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from volume_census_levitus.m of Remi Tailleux (08.03.2013) # last modified : 20.09.2018 # Content/Description #################################################################### # Usage : from mycalc_ocean import volume_census_levitus # Input : grd, s, sbin, th, thbin # grd : Input grid [python class] # s : salinity [psu], 3d # sbin : bin range for s to construct s-th-diagram [min,max,delta] # th : potential temperature [deg C], 3d # thbin : bin range for th to construct s-th-diagram [min,max,delta] # Output : pdf_vol_lev, vol_ocean_lev # pdf_vol_lev : # vol_ocean_lev : # ######################################################################################## print ' Calculate ocean volume in pot.temp.-salinity-space' # Load modules %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% import numpy as np from myconstants import r_earth from mydata_classes import Grid # class containing grid information from scipy.interpolate import interpn # interpolation in 3d # error checking %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if grd.symlon == 0: raise ValueError('*** Input grid must span full longitude grid ***') # Allocation / Initialization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # Compute the number of points in each categories n_sbin = np.ceil((sbin[1]-sbin[0])/sbin[2]) n_thbin = np.ceil((thbin[1]-thbin[0])/thbin[2]) # Allocate pdf_vol_lev (pdf in s-th-space) pdf_vol_lev = np.zeros([n_sbin, n_thbin]) # Define a finer mesh associated with known volume elements that are used %%%%%%%%%%%%%%%% print ' -> Define fine mesh to compute ocean volumes for fine Th-s-space' res_xfine = 0.25 res_yfine = 0.25 res_zfine = 10.0 if grd.lon[1]-grd.lon[0] > res_xfine: lon_fine = np.arange(res_xfine/2.0,360,res_xfine) else: lon_fine = grd.lon if grd.lat[1]-grd.lat[0] > res_yfine: lat_fine = np.arange(-90.0+res_yfine/2.0,90.0,res_yfine) else: lat_fine = grd.lat z_max = max(grd.z)+max(grd.z)%(10**np.floor(np.log10(max(grd.z)))) n_z = np.ceil(z_max/res_zfine) n_z_sub = n_z/10.0 z_range = z_max/n_z_sub vol_tot_ocean = 0.0 # Compute the volume associted with each elements as a function of l %%%%%%%%%%%%%%%%%%%%% print ' *** MISSING: make volume z-dependent ***' grd_fine = Grid(0,lat_fine,0.0,[],1,len(lat_fine),[],[],[],[],[],[]) grd_fine.symlon = 1 vol_fine_lat = (r_earth*np.pi/180.0)**2 * np.cos(np.pi/180.0*lat_fine) * \ res_xfine * res_yfine * res_zfine # Calculate ocean-volume that can be associated with bins in th-s-space %%%%%%%%%%%%%%%%%% print ' -> Interpolation of Th and s to finer grid and compute associated ocean volumes' print ' -> Vertical levels: '+ str(int(n_z)) + ', splitted into ' +\ str(int(n_z_sub)) + ' subdomains' for i_subdom in range(0,int(n_z_sub)): print ' -> Subdomain ' + str(int(i_subdom+1)) + ' of ' + str(int(n_z_sub)) # create meshgrid for interpolation z_beg = res_zfine/2.0 + i_subdom * z_range z_end = res_zfine/2.0 + (i_subdom+1) * z_range - 1 z_fine = np.arange(z_beg,z_end,res_zfine) yfine_xyz, xfine_xyz, zfine_xyz = np.meshgrid(lat_fine, lon_fine, z_fine) # interpolate s and th onto finer grid reshape_form = [len(lon_fine),len(lat_fine),len(z_fine)] xyz_fine_pt = np.array([np.reshape(xfine_xyz,[np.prod(reshape_form)]),\ np.reshape(yfine_xyz,[np.prod(reshape_form)]),\ np.reshape(zfine_xyz,[np.prod(reshape_form)])]).T s_fine = np.reshape(interpn((grd.lon,grd.lat,grd.z), s, xyz_fine_pt, \ bounds_error=False, fill_value=None),reshape_form) th_fine = np.reshape(interpn((grd.lon,grd.lat,grd.z), th, xyz_fine_pt, \ bounds_error=False, fill_value=None),reshape_form) del xyz_fine_pt #print ' replace this np.nan stuff with correct fill value evaluation' #s_fine[s_fine<0 + s_fine>1e6] = np.nan #th_fine[th_fine<-1e6 + th_fine>1e6] = np.nan #s_fine[s_finesbin[1]] = np.nan #th_fine[th_finethbin[1]] = np.nan # associate s-bin and th-bin integers to s_fine and th_fine i_s = np.ceil((s_fine-sbin[0])/sbin[2]) i_th = np.ceil((th_fine-thbin[0])/thbin[2]) # write associated ocean volume into s- and th-bins print ' *** replace the loops by another method ***' for i_xf in range(0,reshape_form[0]): for i_yf in range(0,reshape_form[1]): for i_zf in range(0,reshape_form[2]): i_sd = i_s[i_xf,i_yf,i_zf] i_thd = i_th[i_xf,i_yf,i_zf] if (i_sd == i_sd) and (i_thd == i_thd): i_sd = int(np.max([np.min([i_sd,n_sbin-1]),0])) i_thd = int(np.max([np.min([i_thd,n_thbin-1]),0])) pdf_vol_lev[i_sd,i_thd] = pdf_vol_lev[i_sd,i_thd] + \ vol_fine_lat[i_yf] vol_tot_ocean = vol_tot_ocean + vol_fine_lat[i_yf] del grd_fine # Return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% return pdf_vol_lev, vol_tot_ocean def volume_etopo(fn_topo,dz=100.0): # volume_etopo(fn_topo) ################################################################## # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from volume_etopo.m of Remi Tailleux (09/2018) # last modified : 24.09.2018 # Content/Description #################################################################### # Usage : from mycalc_ocean import volume_etopo # Input : fn_topo, dz # dz : ver. step for interp. [m] # fn_topo : filename (to save Output) [string] # : Output will not be saved, if filename=[] # Output : pp_area, pp_volume # pp_area : z-dep. ocean area [interpolant] # pp_volume : z-dep. ocean volume [interpolant] # Test - compare with : https://www.ngdc.noaa.gov/mgg/global/etopo1_ocean_volumes.html # pp_area(0) = 361869395859664.3 for r_earth = 6371km # pp_volume(10000) = 1.3351019116929976e+18 for r_earth = 6371km # ######################################################################################## # Manual input print ' *** Why should I use a better topography representation? ***' # Load modules %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% import numpy as np from myconstants import r_earth # necessary physical constant from netCDF4 import Dataset # to read netcdf files from scipy import interpolate # to create the interpolants for area and volume # Read bathymetry data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% print ' Calculate z-dependent ocean area and volume using GEBCO data.' print ' GEBCO Bathymetry data downloaded from https://www.bodc.ac.uk.' print ' -> Read bathymetry data' # Read data fn_bathy = '/glusterfs/inspect/users/xg911182/data/bathymetry_data/GEBCO_2014_2D.nc' nc_fid = Dataset(fn_bathy,'r') h_topo = nc_fid.variables['elevation'][:].T lon_topo = nc_fid.variables['lon'][:] lat_topo = nc_fid.variables['lat'][:] # Define data properties (resolution) dlon = np.diff(lon_topo[0:2]) dlat = np.diff(lat_topo[0:2]) n_lontopo = len(lon_topo) n_lattopo = len(lat_topo) # Compute area as function of z %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% print ' -> Caclulate z-dependent ocean area' # calculate z-independent ocean area print ' *** horizontal ocean area z-independent: change this? ***' deg_to_rad = np.pi/180.0 area_y = (deg_to_rad * r_earth)**2 * dlon * dlat * np.cos(deg_to_rad * lat_topo) area_xy = np.tile(area_y,[n_lontopo,1]) del area_y # include z-dependence n_ztopo = int(np.ceil(-np.min(h_topo[:])/dz)) area_z = np.zeros(n_ztopo) topo_z = np.arange(0,n_ztopo)*dz for i_iter in range(0,n_ztopo): print ' ' + str(i_iter) + ' of ' + str(n_ztopo-1) area_dummy = area_xy[h_topo<=-topo_z[i_iter]] area_z[i_iter] = np.sum(area_dummy) del area_dummy # Compute Volume as function of z %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% print ' -> Calculate z-dependent ocean volume' vol_z = 0.5 * (area_z[1::] + area_z[0:-1]) * dz vol_z = np.append([0],np.cumsum(vol_z,axis=0)) # Save interpolants [to reduce calculation time] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% print ' -> Calculate and save interpolants for z-dep area and volume' # Define interpolants pp_area = interpolate.PchipInterpolator(topo_z, area_z, axis=0) pp_volume = interpolate.PchipInterpolator(topo_z, vol_z, axis=0) # Save interpolants if fn_topo: import pickle with open(fn_topo, 'wb') as output: pp_area_volume = [pp_area, pp_volume] pickle.dump(pp_area_volume, output, -1) else: print ' Computed z-dep. topography (vol. and area) will not be saved' # Calculate and return interpolants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% return pp_area, pp_volume def volume_from_sth(rho_i,vol_pdf,tt,pk,s_lower,sbin,reference='jackettetal2004'): # volume_from_sth(rho_i,vol_pdf,tt,pk,s_lower,sbin,reference='jackettetal2004') ########## # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from vol_ts_bis.m of Remi Tailleux (08.03.2013) # last modified : 25.09.2018 # Content/Description #################################################################### # This function computes the volume comprised between the two salinity curves s1 and s2 in # the S-Theta-diagram. s1 is given by the input s_lower, s2 is calculated herein by using # the function calc_sal_from_rho(rho,tt,pk), where pk is the target pressure value. # The volume for each potential temp./sal. classes is given by vol_pdf[0:ns-1,0:nth-1], # where the first index represents the salinity class and the second the pot. temp. class. # Indices defining the lower and upper salinity curves are obtained from the index # np.ceil((salinity_value-sbin[0])/sbin[2]), where sbin[0] is the overall minimum salinity # value to construct the pdf function vol_pdf. Salinity values that are NaN are either # reset to sbin[0] or sbin[1] (max value of pdf) values for the purpose of integration. # usage: from mycalc_ocean import volume_from_sth # (,reference=) # defined by publication, e.g. author1year if there is only one author, # author1andauthor2year for two authors or author1etalyear for more than 2 authors # input: rho_i, vol_pdf, tt, pk, s_lower, sbin, reference # pk : target gauge pressure [dbar], one value # (absolute pressure - 10.1325) # reference : polynomials to calc s [string] # rho_i : denisty [kg m**-3], one value # sbin : bin range for s of s-th-diagram [min,max,delta] # s_lower : lower salinity curve [psu], 1d # tt : pot. temp. array / axis of s-th-diagram [deg C] # vol_pdf : pdf of ocean-vol-fraction in s-th diagram [0<=value<=1], 2d # output: sal : salinity [psu] # check: calc_rho_from_theta(20,20,1000) = 1017.7288680196421 # ######################################################################################## # load necessary modules %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% import numpy as np from scipy.optimize import fsolve # to search for 0-points in a function # Allocation / Initialization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% vol_ts = 0.0 n_th = len(tt) # Compute salinity curve and replace NaN-values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% s_rho = calc_sal_from_rho(rho_i*np.ones(n_th),tt,pk*np.ones(n_th)) # Identify NaNs and set them either to min or max whether it is the first or end part that # is NaN if np.isnan(s_lower[0]): s_lower[np.isnan(s_lower)] = sbin[0] elif np.isnan(s_lower[-1]): s_lower[np.isnan(s_lower)] = sbin[1] if np.isnan(s_rho[0]): s_rho[np.isnan(s_rho)] = sbin[0] elif np.isnan(s_rho[-1]): s_rho[np.isnan(s_rho)] = sbin[1] s_rho[s_rho > sbin[1]] = sbin[1] s_lower[s_lower > sbin[1]] = sbin[1] # Compute volume in s-th-space between the two salinity curves %%%%%%%%%%%%%%%%%%%%%%%%%%% print ' *** can I get rid of this loop? ***' for i_iter in range(0,n_th): is_class = int(np.ceil((s_rho[i_iter]-sbin[0])/sbin[2])) is_lower = int(np.ceil((s_lower[i_iter]-sbin[0])/sbin[2])) if is_class > is_lower: vol_ts = vol_ts + np.sum(vol_pdf[is_lower+1:is_class+1,i_iter]) # Return volume fraction from s-th diagram %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% return vol_ts