diff --git a/calc_ref_state.py b/calc_ref_state.py index 162ebf7..f3d716e 100644 --- a/calc_ref_state.py +++ b/calc_ref_state.py @@ -27,6 +27,8 @@ import numpy as np # import time # to measure elapsed time # read data from get_data import get_data_woa2009 as read_data +# calculation modules +from mycalc_ocean import calc_rho_from_theta # calc of density from theta # plotting modules / functions from myplot import myplot_2dmap # plot of lat-lon fields from myplot import myplot_2d # plot of crossections fields @@ -44,13 +46,14 @@ from myplot_inputs import myplot_create_pn # create plotname # information about input data print 'MANUAL DEFINED OUTPUTS MISSING' dim = ('lon','lat','z','time',) -variables = {'s':{},'temp':{},'grd':[],'valid':[]} +variables = {'s':{},'theta':{},'grd':[],'valid':[]} # plotting input +plot_var = {'rho':[]} saveplot = input('Do you want to save the plots? Enter 1 for yes or 0 for no') dir_plot = '/glusterfs/inspect/users/xg911182/Code/Python/plots_test/' lonlat_reg = [25.0, 25.0+360.0, -90.0, 90.0] # plotting range for 2dmap plot_2d = 1 # plot 2d fields (lat-lon) -plot_cs = 0 # plot crosssections +plot_cs = 1 # plot crosssections cs_plot = {'lon':[-47.5,172.5],'lat':[0]} # crosssections # Input 1 #################################################################################### # Physical constants @@ -77,7 +80,7 @@ from myconstants import rho0, grav beg_time = time.time() print '' print '----------------------------------------------' -print ' Read data' +print 'Read data' grd, data = read_data(variables,dim) print '----------------------------------------------' end_time = time.time() @@ -88,11 +91,11 @@ print '----------------------------------------------' # Define grid for data from step a # input : lon, lat, zlev (all from step a) # output : lat_xy, lat_yz, lat_xyz, lon_xy, lon_xyz, p_xyz, plev, z_yz -# lat_* : meshgrid for latitude in x/y/z-dir [degrees] -# lon_* : meshgrid for longitude in x/y/z-dir [degrees] -# p_* : meshgrid for pressure in x/y/z-dir [dbar] +# lat_* : meshgrid for latitude in x/y/z-dim [degrees] +# lon_* : meshgrid for longitude in x/y/z-dim [degrees] +# p_xyz : meshgrid for pressure in xyz-dim [dbar] # plev : gauge pressure (p_abs-10.1325 dbar) [dbar] -# z_* : meshgrid for depth in x/y/z-dir [m] +# z_yz : meshgrid for depth in yz-dim [m] # ############################################################################################ beg_time = time.time() print '' @@ -118,14 +121,28 @@ print '----------------------------------------------' # input : salt, theta, plev, valid # salt : output of step a # theta : output of step a -# p_lev : output of step b +# plev : output of step b # output : dens # dens : density of water [kg m**-3] # ############################################################################################ beg_time = time.time() print '' print '----------------------------------------------' -print 'Insert code to calculate density' +print 'Calculate density' +# Theta mask for valid values +mask = data['theta']['valid'] +# Allocate density +data['rho'] = {} +data['rho']['val'] = data['theta']['val'] +data['rho']['fill_value'] = data['theta']['fill_value'] +data['rho']['standard_name'] = 'sea_water_density' +data['rho']['valid'] = mask +data['rho']['units'] = 'kg m**-3' +# calculate density +data['rho']['val'][mask] = calc_rho_from_theta(data['s']['val'][mask],\ + data['theta']['val'][mask],p_xyz[np.squeeze(mask)]) +# De;ete unnecessary variables +del mask print '----------------------------------------------' end_time = time.time() print 'Elapsed time to calculate density: '+ '{:.2f}'.format(end_time-beg_time)+'s' @@ -185,7 +202,7 @@ print '----------------------------------------------' # step g ##################################################################################### # plot of full 2d fields -# input : data, date_str, dir_plot, grd, i_dep, i_t, lonlat_reg +# input : data, date_str, dir_plot, grd, i_dep, i_t, lonlat_reg, plot_var # data : output of step a # date_str : MISSING - DONE LOCALLY HERE # dir_plot : output of Input 0 @@ -193,6 +210,7 @@ print '----------------------------------------------' # i_dep : MISSING - DONE LOCALLY HERE # i_t : MISSING - DONE LOCALLY HERE # lonlat_reg: output of Input 0 +# plot_var : output of Input 0 # output : None # saved/created : plots saved in dir_plot # ############################################################################################ @@ -200,7 +218,7 @@ beg_time = time.time() print '' print '----------------------------------------------' print 'Plot 2d fields (lat-lon)' -list_allkeys = data.keys() # all variables +list_allkeys = plot_var.keys() # all variables i_dep = 0 i_t = 0 date_str = '2009' @@ -223,7 +241,7 @@ print '----------------------------------------------' # step h ##################################################################################### # plot of crossection fields -# input : cs_plot, data, date_str, dir_plot, grd, i_t, plot_cs, saveplot +# input : cs_plot, data, date_str, dir_plot, grd, i_t, plot_cs, plot_var, saveplot # cs_plot : output of Input 0 # data : output of step a # date_str : MISSING - DONE LOCALLY HERE @@ -231,6 +249,7 @@ print '----------------------------------------------' # grd : output of step a # i_t : MISSING - DONE LOCALLY HERE # plot_cs : output of Input 0 +# plot_var : output of Input 0 # saveplot : output of Input 0 # output : None # saved/created : plots saved in dir_plot @@ -242,7 +261,7 @@ if plot_cs: print '' print '----------------------------------------------' print 'Plot lon-p and lat-p crossections' - list_allkeys = data.keys() # all variables + list_allkeys = plot_var.keys() # all variables for i_key in range(0,len(list_allkeys)): for i_cs in range(0,len(cs_plot['lat'])): # load fixed latitude and associated index for grd.lat diff --git a/mycalc_ocean.py b/mycalc_ocean.py new file mode 100644 index 0000000..e9dc226 --- /dev/null +++ b/mycalc_ocean.py @@ -0,0 +1,368 @@ +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 rho + + + +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 ****') + # Calculate first estimate for potential temperature + print 'Is this above statemtent true? First estimate?' + n0 = 0 + n2 = 2 + # Allocate values + s1 = s*35.0/35.16504 + t1 = temp + p1 = p + pr1 = pr + # First estimate for pot. temp. + print 'Where is this equation coming from?' + theta = t1 + (p1-pr1)*( 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)) + de_dt = 13.6 + # Iteration (while loop) to determine theta + nloops = 3 # default + # NOTE: nloops = 1 gives theta with a maximum error of ... + # nloops = 2 gives theta with a maximum error of ... + n = 0 + print 'Where are the equations in this while loop coming from?' + while n <= nloops: + n = n+1 + theta_old = 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 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