gabriel/Test
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Update of calc_ref_state.py and mycalc_ocean.py to calculate density from theta

Browse Source
This commit is contained in:
Gabriel Wolf 2018-09-20 12:11:07 +01:00
parent 0251f30fa0
commit e0246c86bd
2 changed files with 400 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
calc_ref_state.py
View File

@ -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

368
mycalc_ocean.py Normal file
View File

@ -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=<reference>)
# <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