From 164ed483158ac338123d6582ec0600875f63244a Mon Sep 17 00:00:00 2001 From: Gabriel Wolf Date: Wed, 26 Sep 2018 11:50:20 +0100 Subject: [PATCH] First complete version, BUT STILL WITH BUGS --- calc_ref_state.py | 186 ++++++++++++++------ get_data.py | 53 +++--- gw_ocean_refstate.py | 215 ++++++++++++++++++----- input_calc_ref_state.py | 14 ++ mycalc_ocean.py | 380 +++++++++++++++++++++++++++++++++++----- myconstants.py | 2 +- myplot.py | 3 + 7 files changed, 685 insertions(+), 168 deletions(-) create mode 100644 input_calc_ref_state.py diff --git a/calc_ref_state.py b/calc_ref_state.py index 8894d48..20d6aef 100644 --- a/calc_ref_state.py +++ b/calc_ref_state.py @@ -29,36 +29,39 @@ import time # to measure elapsed time 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 -from mycalc_ocean import calc_area_lat # calc horizontal area +from mycalc_ocean import calc_area_xyz # calc horizontal area from gw_ocean_refstate import refstate_meandensity, refstate_sorted_stheta # calc. ref. states # plotting modules / functions from myplot import myplot_2dmap # plot of lat-lon fields from myplot import myplot_2d # plot of crossections fields from myplot_inputs import myplot_create_pn # create plotname +import pickle # to save reference state and create anomalies +# Further modules +import os # to execute terminal commands # ############################################################################################ # Input 0 #################################################################################### # Manual defined inputs # input : None -# output : dir_plot, variables -# dir_plot : directory, where plots will be saved -# variables : variables necessary for calculations +# output : dim, dir_plot, saveplot, variables, REST MISSING +# dim : prefered dimensions for matrices [list] +# dir_plot : directory, where plots will be saved [string] +# saveplot : defines if plots are being saved [0 or 1] +# variables : variables necessary for calculations [MISSING] # MANUAL DEFINED OUTPUTS MISSING # ############################################################################################ + print '----------------------------------------------' -print 'Manual defined input data' +print 'Read manual defined input data' +import input_calc_ref_state # information about input data print ' *** MANUAL DEFINED OUTPUTS MISSING ***' dim = ('lon','lat','z','time',) variables = {'s':{},'theta':{},'grd':[],'valid':[]} -# plotting input -plot_var = {'rho':[]} -saveplot = input(' Save 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 = 1 # plot crosssections -cs_plot = {'lon':[-47.5,172.5],'lat':[0]} # crosssections +print ' Save plots? Enter 1 for YES or 0 for NO.' +saveplot = input(' -> USER INPUT: ') +print ' Define case with string (used to specify plots and data), e.g. \'_id000\'' +case_id = input(' -> USER INPUT: ') print '----------------------------------------------' # Input 1 #################################################################################### @@ -143,15 +146,15 @@ print 'Calculate density' mask = data['theta']['valid'] # Allocate density data['rho'] = {} -data['rho']['val'] = data['theta']['val'] +data['rho']['val'] = np.zeros_like(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 + data['theta']['val'][mask],p_xyz[np.squeeze(mask)],reference=ref_EOS) +# Delete unnecessary variables del mask print '----------------------------------------------' end_time = time.time() @@ -171,7 +174,7 @@ beg_time = time.time() print '' print '----------------------------------------------' print 'Calculate z-dependent ocean area' -area_xyz = calc_area_lat(grd,r_earth) +area_xyz = calc_area_xyz(grd,r_earth) print '----------------------------------------------' end_time = time.time() print 'Elapsed time to calc. ocean area: '+ '{:.2f}'.format(end_time-beg_time)+'s' @@ -190,33 +193,50 @@ print '----------------------------------------------' # pp_rhor_botup : interpolant for ref. density (sorted in th-s-space, calc. bottom-up) # pp_rhor_topdn : interpolant for ref. density (sorted in th-s-space, calc. top-down) # ############################################################################################ + +# THIS MUST BE DONE AUTOMATICALLY i_t = 0 + +beg_time = time.time() +print '' +print '----------------------------------------------' +print 'Calculate reference state for mean density' + +# Define input data +mask = data['rho']['valid'][:,:,:,i_t] + # define target depth n_target = grd.Nz+9 ztarget = np.zeros(n_target) ztarget[0:10] = np.arange(0,10) ztarget[10:] = grd.z[1:] -beg_time = time.time() -print '' -print '----------------------------------------------' -print 'Calculate reference state for mean density' -# Define input data -mask = data['rho']['valid'][:,:,:,i_t] + # calculation of reference profiles pp_rhor = refstate_meandensity(data['rho']['val'][:,:,:,i_t]*mask,grd,area_xyz*mask) end_time = time.time() beg_time2 = time.time() -print 'Calculate sorted density reference state' + # Define input data rho = data['rho']['val'] s = data['s']['val'] th = data['theta']['val'] mask = data['rho']['valid']*data['s']['valid']*data['theta']['valid'] + # calculation of reference profiles -pp_rhor_botup, pp_rhor_topdn = refstate_sorted_stheta(th,s,rho,mask,ztarget,grd) -input('Stop here') +print 'Calculate sorted density reference state' +pp_rhor_botup, pp_rhor_topdn = refstate_sorted_stheta(th, s, rho, mask, ztarget, grd,\ + ref_EOS=ref_EOS, sbin=[34.0,35.0,0.002], thbin=[-2.6,3.0,0.005],\ + test=True, case_str=case_id) + +# Save output of different reference states +print 'Save reference state' +with open('data_ref_state' + case_id + '.pickle', 'wb') as output: + pp_refstate = [pp_rhor_mean, pp_rhor_botup, pp_rhor_topdn] + pickle.dump(pp_refstate, output) + # Deallocation del th, s, rho + print '----------------------------------------------' end_time2 = time.time() print 'Elapsed time to calculate reference state: ' @@ -231,13 +251,58 @@ print '----------------------------------------------' # saved/created : plots saved in dir_plot # dir_plot : output of Input 0 # ############################################################################################ + beg_time = time.time() print '' print '----------------------------------------------' -print 'Insert code to plot reference state' +print 'Plot z-dependent reference state' + +# Create plotname +date_str = 'WOA2009' +pn_full = myplot_create_pn({'dir_plot':dir_plot,'beg_str':'Refstate',\ + 'varname':'rho','date':date_str,'end_str':case_id+'.png'}) +print 'Plot deviation of ref-state to a previous calculated ref-state?' +plot_dev = input('-> USER INPUT: ') +if plot_dev: + print 'Please enter one of the listed filenames as comparison reference state.' + os.system("ls data_ref_state*.pickle") + fn_refstate = input('-> USER INPUT: ') + fn_diff = [li for li in list(difflib.ndiff(fn_refstate,\ + 'data_ref_state.pickle')) if li[0] != ' '] + pn_anom = myplot_create_pn({'dir_plot':dir_plot,'beg_str':'Refstate_dev'+fn_diff,\ + 'varname':'rho','date':date_str,'end_str':case_id+'.png'}) + +# Define plot input +xxx = np.zeros([3,n_target]) +xxx[0,:] = pp_rhor_mean(ztarget) +xxx[1,:] = pp_rhor_botup(ztarget) +xxx[2,:] = pp_rhor_topdn(ztarget) +yyy = np.tile(-ztarget,[3,1]) +if plot_dev: + xxx_anom = np.zeros([3,n_target]) + with open(fn_refstate, 'rb') as data_ref: + dummy = pickle.load(data_ref) + xxx_anom[0,:] = xxx[0,:] - dummy[0](ztarget) + xxx_anom[1,:] = xxx[1,:] - dummy[1](ztarget) + xxx_anom[2,:] = xxx[2,:] - dummy[2](ztarget) + +# Plotting +myplot_1d(xxx,yyy,FS=22,col_c=('b','r','k'),saveplot=saveplot,\ + label_c=('mean','bottom-up','top-down'),xlabel_c='rho_ref [kg m**-3]',\ + ylabel_c='z_ref [m]',plotname=pn_full) +if plot_dev: + myplot_1d(xxx,yyy,FS=22,col_c=('b','r','k'),saveplot=saveplot,\ + label_c=('mean','bottom-up','top-down'),xlabel_c='rho_ref [kg m**-3]',\ + ylabel_c='z_ref [m]',plotname=pn_anom) + +# Deallocation +if plot_dev: + del plot_anom, fn_refstate, fn_diff, dummy, xxx_anom +del plot_dev, pn_full, xxx, yyy + print '----------------------------------------------' end_time = time.time() -print 'Elapsed time to plot and save reference state: '+ '{:.2f}'.format(end_time-beg_time)+'s' +print 'Elapsed time to plot reference state: '+ '{:.2f}'.format(end_time-beg_time)+'s' print '----------------------------------------------' # step g ##################################################################################### @@ -254,30 +319,41 @@ print '----------------------------------------------' # output : None # saved/created : plots saved in dir_plot # ############################################################################################ -beg_time = time.time() -print '' -print '----------------------------------------------' -print 'Plot 2d fields (lat-lon)' -list_allkeys = plot_var.keys() # all variables -i_dep = 0 -i_t = 0 -date_str = '2009' -for i_key in range(0,len(list_allkeys)): - # create plotname - plotname = myplot_create_pn({'dir_plot':dir_plot,'beg_str':'Map2d',\ - 'varname':list_allkeys[i_key],'zlev':[grd.z[i_dep],grd.Uz],\ - 'lonlat':lonlat_reg,'date':date_str,'end_str':'.png'}) - # define plot input - title_in = 'Depth='+'{:.1f}'.format(grd.z[i_dep])+grd.Uz - data_in = data[list_allkeys[i_key]]['val'][:,:,i_dep,i_t] - data_in[data_in==data[list_allkeys[i_key]]['fill_value']] = np.nan - myplot_2dmap(data_in,grd,lonlat_range=lonlat_reg,saveplot=1,\ - title_c=title_in,plotname=plotname,unit_data=data[list_allkeys[i_key]]['units']) - del data_in, title_in, plotname -print '----------------------------------------------' -end_time = time.time() -print 'Elapsed time to plot 2d fields (lat-lon): '+ '{:.2f}'.format(end_time-beg_time)+'s' -print '----------------------------------------------' + +if plot_2d: + beg_time = time.time() + print '\n----------------------------------------------' + print 'Plot 2d fields (lat-lon)' + + # Define Input for all following 2d-plots + list_allkeys = plot_var.keys() # all variables + i_dep = 0 + i_t = 0 + date_str = 'WOA2009' + + # Loop over all variables + for i_key in range(0,len(list_allkeys)): + # Create plotname + plotname = myplot_create_pn({'dir_plot':dir_plot,'beg_str':'Map2d',\ + 'varname':list_allkeys[i_key],'zlev':[grd.z[i_dep],grd.Uz],\ + 'lonlat':lonlat_reg,'date':date_str,'end_str':case_id+'.png'}) + + # Define plot input + title_in = 'Depth='+'{:.1f}'.format(grd.z[i_dep])+grd.Uz + data_in = data[list_allkeys[i_key]]['val'][:,:,i_dep,i_t] + data_in[data_in==data[list_allkeys[i_key]]['fill_value']] = np.nan + + # Plotting + myplot_2dmap(data_in,grd,lonlat_range=lonlat_reg,saveplot=1,\ + title_c=title_in,plotname=plotname,unit_data=data[list_allkeys[i_key]]['units']) + + # Deallocation + del data_in, title_in, plotname + + print '----------------------------------------------' + end_time = time.time() + print 'Elapsed time to plot 2d fields (lat-lon): '+ '{:.2f}'.format(end_time-beg_time)+'s' + print '----------------------------------------------' # step h ##################################################################################### # plot of crossection fields @@ -294,7 +370,7 @@ print '----------------------------------------------' # output : None # saved/created : plots saved in dir_plot # ############################################################################################ -date_str = '2009' +date_str = 'WOA2009' i_t = 0 if plot_cs: beg_time = time.time() @@ -310,7 +386,7 @@ if plot_cs: # create plotname plotname = myplot_create_pn({'dir_plot':dir_plot,'beg_str':'CS',\ 'varname':list_allkeys[i_key],'lonlat':[0,360,lat_fix,lat_fix],\ - 'date':date_str,'end_str':'.png'}) + 'date':date_str,'end_str':case_id+'.png'}) # define plot input title_in = 'Depth='+'{:.1f}'.format(grd.z[i_dep])+grd.Uz cbarlabel_in = data[list_allkeys[i_key]]['standard_name'] + \ @@ -328,7 +404,7 @@ if plot_cs: # create plotname plotname = myplot_create_pn({'dir_plot':dir_plot,'beg_str':'CS',\ 'varname':list_allkeys[i_key],'lonlat':[lon_fix,lon_fix,-90,90],\ - 'date':date_str,'end_str':'.png'}) + 'date':date_str,'end_str':case_id+'.png'}) # define plot input title_in = 'Depth='+'{:.1f}'.format(grd.z[i_dep])+grd.Uz cbarlabel_in = data[list_allkeys[i_key]]['standard_name'] + \ diff --git a/get_data.py b/get_data.py index a61b8fc..0b189f3 100644 --- a/get_data.py +++ b/get_data.py @@ -17,24 +17,28 @@ def get_data_woa2009(variables,grid_order=('lon','lat','z','time',)): # Ed. NOAA Atlas NESDIS 69, 184 pp. # -> PDF : https://www.nodc.noaa.gov/OC5/indpub.html#woa09 # ######################################################################################## - print 'Use of get_data_woa2009 in get_data.py' + print ' Load WOA2009 data' dir_data = '/glusterfs/inspect/users/xg911182/data/WOA2009/' grid_names = ('lon','lat','depth','time',) - # import modules + + # import modules %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% import numpy as np from netCDF4 import Dataset # to read netcdf files from mydata_classes import Grid # class containing grid information import myconstants as my_const # my own defined constants from mycalc_ocean import calc_theta_from_temp - # define grid_order (used for data permutation) + + # Allocation and Initialization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + # define grid_order (used for data permutation) find_index = lambda searchlist, elem: [[i for i, x in enumerate(searchlist) if x == e] for e in elem] grid_order = np.squeeze(np.array(find_index(('lon','lat','z','time',),grid_order))) + list_allkeys = variables.keys() # Allocation data = {} grd = [] - list_allkeys = variables.keys() - # Read grid information (for all variables identical) - print 'Check if all WOA data is using the same grid' + + # Read grid information (for all variables identical) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + print ' -> Read grid information' if 'grd' in list_allkeys: fn = dir_data + 'temperature_annual_1deg.nc' nc_fid = Dataset(fn,'r') @@ -49,9 +53,10 @@ def get_data_woa2009(variables,grid_order=('lon','lat','z','time',)): grd = Grid(LON,LAT,Z,TIME,len(LON),len(LAT),len(Z),len(TIME), Ulon, Ulat, Uz, Ut) # define data permutation to match input dimensions nc_fid.variables['t_an'].dimensions - # read salinity data + + # read salinity data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if 's' in list_allkeys: - print 'Load salinity data' + print ' -> Read salinity data' fn = dir_data + 'salinity_annual_1deg.nc' nc_fid = Dataset(fn,'r') # read salinity data @@ -65,9 +70,10 @@ def get_data_woa2009(variables,grid_order=('lon','lat','z','time',)): data[varname]['fill_value'] = DATA_d._FillValue data[varname]['standard_name'] = DATA_d.standard_name data[varname]['valid'] = MASK.transpose(data_perm) - # Read temperature data + + # Read temperature data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if 'temp' or 'theta' in list_allkeys: - print 'Load temperature data' + print ' -> Read temperature data' fn = dir_data + 'temperature_annual_1deg.nc' nc_fid = Dataset(fn,'r') # read temperature data @@ -83,32 +89,25 @@ def get_data_woa2009(variables,grid_order=('lon','lat','z','time',)): data[varname]['fill_value'] = DATA_d._FillValue data[varname]['standard_name'] = DATA_d.standard_name data[varname]['valid'] = MASK.transpose(data_perm) + + # compute potential temperature for valid points ------------------------------------- if 'theta' in list_allkeys: - # compute pot. temp. for valid points - print DATA_d.dimensions - print data_perm - THETA = TEMP.transpose(data_perm) + print ' -> Compute potential temperature from temperature' + THETA = np.copy(TEMP).transpose(data_perm) SA = data['s']['val'] SA = SA[data['s']['valid']] - print DATA_d.dimensions[data_perm[0]] - print DATA_d.dimensions[data_perm[1]] - print DATA_d.dimensions[data_perm[2]] - print DATA_d.dimensions[data_perm[3]] wvar = nc_fid.variables[DATA_d.dimensions[data_perm[0]]] xvar = nc_fid.variables[DATA_d.dimensions[data_perm[1]]] yvar = nc_fid.variables[DATA_d.dimensions[data_perm[2]]] zvar = nc_fid.variables[DATA_d.dimensions[data_perm[3]]] w3d, x3d, y3d, z3d = np.meshgrid(xvar,wvar,yvar,zvar) p = y3d*(my_const.rho0*my_const.grav/1e4) # in dbar - print 'Doesnt work automatically - p_xyz' del w3d, x3d, y3d, z3d, wvar, xvar, yvar, zvar - pr = p*0 - print 'SA.shape: ',SA.shape - print 'THETA.shape',THETA[data[varname]['valid']].shape - print 'p.shape, ',p[data[varname]['valid']].shape - dummy = input('Press enter to continue') + pr = p*0.0 + print ' *** why do I set p_ref=0? ***' MASK_mesh = data[varname]['valid'] - THETA[MASK_mesh] = calc_theta_from_temp(SA,THETA[MASK_mesh],p[MASK_mesh],pr[MASK_mesh]) + THETA[MASK_mesh] = calc_theta_from_temp(SA,THETA[MASK_mesh],\ + p[MASK_mesh],pr[MASK_mesh]) varname = 'theta' data[varname] = {} data[varname]['val'] = THETA @@ -116,8 +115,8 @@ def get_data_woa2009(variables,grid_order=('lon','lat','z','time',)): data[varname]['fill_value'] = DATA_d._FillValue data[varname]['standard_name'] = 'sea_water_potential_temperature' data[varname]['valid'] = MASK.transpose(data_perm) - print 'MISSING: Calculate pot. temp. from temperature data' - # return data + + # return data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if 'grd' in list_allkeys: return grd, data else: diff --git a/gw_ocean_refstate.py b/gw_ocean_refstate.py index 6ac71cc..a0179ec 100644 --- a/gw_ocean_refstate.py +++ b/gw_ocean_refstate.py @@ -15,6 +15,7 @@ def refstate_meandensity(rho_in,grd,area_xyz_in): # Load modules import numpy as np from scipy import interpolate + import pickle # To load interpolants of bathymetry data # Allocation rhor = np.zeros([grd.Nz,1]) # compute hor. averaged density field @@ -27,74 +28,202 @@ def refstate_meandensity(rho_in,grd,area_xyz_in): return pp_rhor def refstate_sorted_stheta(th,s,rho,mask,ztarget,grd,sbin=[0.0,43.0,0.002],\ - thbin=[-2.6,30.0,0.005]): - # refstate_sorted_stheta(th,s,rho,mask,ztarget,grd) ###################################### + thbin=[-2.6,30.0,0.005],ref_EOS='jackettetal2004',test=False,case_str=[]): + # refstate_sorted_stheta(th,s,rho,mask,ztarget,grd,test) ################################# # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from get_refstate_sorted.m of Remi Tailleux (10.04.2013) - # last modified : 20.09.2018 + # last modified : 24.09.2018 # Content/Description #################################################################### # Reference : method based on Tailleux (2013), -> https://doi.org/10.1017/jfm.2013.509 # : summary of this framework in Saenz et al. (2015): # -> https://doi.org/10.1175/JPO-D-14-0105.1 # Usage : from gw_ocean_refstate import refstate_sorted_stheta - # Input : grd, mask, rho, s, sbin, thbin, th, ztarget + # Input : grd, mask, ,ref_equofstate, rho, s, sbin, test, thbin, th, ztarget # grd : Input grid [python class] # mask : valid points [boolean], 3d + # ref_EOS : defines equ. of state [reference] + # : reference def. by publication, e.g. author1year for one author-public., + # : author1andauthor2year for 2-author-public. or author1etalyear for more + # : than 2 authors # rho : density [kg m**-3], 3d # s : salinity [psu], 3d - # s_bin : bin range for s to construct s-th-diagram [min,max,delta] + # sbin : bin range for s to construct s-th-diagram [min,max,delta] + # test : further test plotting and value checking [boolean], 1d # th : potential temperature [deg C], 3d - # th_bin : bin range for th to construct s-th-diagram [min,max,delta] + # thbin : bin range for th to construct s-th-diagram [min,max,delta] # ztarget : target height level [m], 1d # Output :pp_rhor_botup, pp_rhor_topdn # pp_rhor_botup : interpolant for density reference state (calc. from bottom up) # pp_rhor_topdn : interpolant for density reference state (calc. from bottom up) + # Testing : plot_sth_pdf + # plot_sth_pdf : plotting of sorted ocena-volumes in the S-TH-diagram # ######################################################################################## - # Load modules + + # Testing variables %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + plot_sth_pdf = 1*test + if plot_sth_pdf == 1: + print ' Create test plot for sorted ocean-volumes in S-Th-diagram' + from myplot import myplot_2d + + # Load modules %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% import numpy as np from scipy import interpolate # for interpolation - from mycalc_ocean import volume_census_levitus # to create s-th-pdf with vol info - # Allocation + from mycalc_ocean import calc_rho_from_theta, volume_census_levitus + # calc_rho_from_theta : to calculate density from pot. temperature + # volume_census_levitus : to create s-th-pdf with vol info + from myconstants import rho0, grav # necessary constants + + # Allocation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% n_target = len(ztarget) s[mask==False] = np.nan th[mask==False] = np.nan rho[mask==False] = np.nan - # Compute the pdf of the salinity/theta data --------------------------------------------- - print ' *** MISSING: computation of the pdf salinity/theta data ***' - pdf_vol_lev, vol_tot_ocean = volume_census_levitus(th, s, grd, thbin, sbin) - # Define properties of the binned temperature/salinity data ------------------------------ - #n_sbin, n_thbin = + + # Compute the pdf of the salinity/theta data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + try: + import pickle + with open('data_pdf_salth' + case_str + '.pickle', 'rb') as data_pdf: + dummy = pickle.load(data_pdf) + pdf_vol_lev = dummy[0] + vol_tot_ocean = dummy[1] + del dummy, data_pdf + pdf_vol_lev_norm = pdf_vol_lev/vol_tot_ocean + except: + pdf_vol_lev, vol_tot_ocean = volume_census_levitus(th, s, grd, thbin, sbin) + print ' Total volume ocean: ' + str(vol_tot_ocean) + ' m**3' + pdf_vol_lev_norm = pdf_vol_lev/vol_tot_ocean + import pickle + with open('data_pdf_salth' + case_str + '.pickle', 'wb') as output: + pdf_salth = [pdf_vol_lev, vol_tot_ocean] + pickle.dump(pdf_salth, output, -1) + + # Define properties of the binned temperature/salinity data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% print ' *** MISSING: property definition of binned temperature/salinity data ***' - # define target pressure from target depths ---------------------------------------------- - print ' *** MISSING: define target pressure from target depths ***' - # compute bounding salinity curves ------------------------------------------------------- - print ' *** MISSING: compute bounding salinity curves ***' - # Seek the global min. and max. of density by bringing all parcels at the surface and at - # the bottom and the oceans ---------------------------------------------------------- - print ' *** MISSING: density profiles ***' - # Compute potential density reference to target pressure at surface ---------------------- - print ' *** MISSING: pot. density at surface ***' - # Compute potential density referenced to last target pressure --------------------------- - print ' *** MISSING: pot. density for last target pressure ***' - # Compute global density minimum and maximum --------------------------------------------- - print ' *** MISSING: min. and max. of global density ***' - # Compute upper and lower salinity curves in theta/salt space ---------------------------- - print ' *** MISSING: upper and lower salinity curves in th-s space ***' - # Compute the topography statistics ------------------------------------------------------ - print ' *** MISSING: topography statistics ***' - # New attempt to compute the top-down profile -------------------------------------------- - print ' *** MISSING: rhor top-down profile ***' - print(' Computing top-dwon profile for reference density') + n_sbin, n_thbin = pdf_vol_lev.shape + s_axis = np.arange(sbin[0]+sbin[2]*0.5,sbin[1],sbin[2]) + th_axis = np.arange(thbin[0]+thbin[2]*0.5,thbin[1],thbin[2]) + # Testing pdf-output (pdf_vol_lev_norm) + if plot_sth_pdf: + d_step = 1 + xxx = np.arange(sbin[0],sbin[1],sbin[2]*d_step) + yyy = yyy=np.arange(thbin[0],thbin[1],thbin[2]*d_step) + ddd = pdf_vol_lev_norm[::d_step,::d_step] + ddd[ddd==0] = np.nan + ddd = np.log10(ddd) + print 'np.nanmax(ddd) = ',np.nanmax(ddd) + import datetime + plotname = 'Testplot_for_oceanvol_in_s_th_diagram_' + \ + str(datetime.date.today()) + case_str + '.png' + myplot_2d(xxx,yyy,ddd.T,xlabel_c='salinity [psu]',ylabel_c='potential temp. [deg C]',\ + plotname=plotname,saveplot=1,FS=14,d_xtick=0.2,title_c='S-TH',\ + cbarlabel_c='Volume freq. distribution [log_10]') + #import pdb + #pdb.set_trace() + #iprint 'stop here' + del xxx, yyy, ddd + + # define target pressure from target depths %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + ptarget = rho0 * grav * ztarget / 1e4 + n_target = len(ptarget) + print ' *** why using constant density for target pressure and not variable one? ***' + + # Seek the global min. and max. of density by bringing all parcels to the ocean surface + # (to find min.) and to ocean bottom (to find max.) -------------------------------------- # Allocation - rhor_topdn = np.zeros(n_target) - # New attempt at computing bottom-up profile --------------------------------------------- + rho_theta_surf = rho + rho_theta_deep = rho + p_surf = np.ones([grd.Nlon,grd.Nlat,grd.Nz])*ptarget[0] + p_deep = np.ones([grd.Nlon,grd.Nlat,grd.Nz])*ptarget[-1] + # Compute potential density reference to target pressure at surface + rho_theta_surf[mask] = calc_rho_from_theta(s[mask],th[mask],p_surf[mask],reference=ref_EOS) + # Compute potential density referenced to last target pressure + rho_theta_deep[mask] = calc_rho_from_theta(s[mask],th[mask],p_deep[mask],reference=ref_EOS) + # Compute global density minimum and maximum + rho_min = np.min(rho_theta_surf[mask]) + rho_max = np.max(rho_theta_deep[mask]) + + # Compute upper and lower salinity curves in theta/salt space %%%%%%%%%%%%%%%%%%%%%%%%%%%% + # Allocation + s_lower = np.zeros(n_thbin) + s_upper = np.zeros(n_thbin) + # Computation + s_lower = calc_sal_from_rho(rho_min,th_axis,ptarget[0],reference=ref_EOS) + s_upper = calc_sal_from_rho(rho_max,th_axis,ptarget[-1],reference=ref_EOS) + + # Compute the topography statistics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + print ' Load topographic data for z-dependent ocean area and volume' + # Load or calculate z-dep. ocean area/volume + fn_topo = 'interp_of_topo' + '.pickle' + try: + with open(fn_topo, 'rb') as data_interp: + dummy = pickle.load(data_interp) + pp_area = dummy[0] + pp_volume = dummy[1] + del dummy, data_interp + except: + print ' -> Data could not be loaded. Therefore it will be calculated.' + pp_area, pp_volume = volume_etopo(fn_topo) + del fn_topo + # define z-dep. volume and its fraction + vol_target = pp_volume(ztarget) + frac_vol_target = vol_target/vol_target[-1] + + # New attempt to compute the top-down profile %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + print ' Computing top-down profile for reference density' + # Allocation/Initialization + rhor_topdn = np.zeros(n_target) + rhor_topdn[0] = rho_min + rhor_topdn[-1] = rho_max + s_curves_td = np.zeros([n_target,n_thbin]) + s_curves_td[0,:] = s_lower + s_curves_td[-1,:] = s_upper + scurv = np.copy(s_lower) + rho_l = np.copy(rho_min) + # Use optimization method to converge to correct salinity curve + print ' *** this loop is not necessary ***' + for i_iter in range(0,n_target-1): + rho_u = calc_rho_from_theta(sbin[1],thbin[0],ptarget[i_iter+1]) + func = lambda rho_i : volume_from_sth(rho_i, pdf_vol_lev_norm,\ + th_axis, ptarget[i_iter+1], scurv, sbin) - \ + (frac_vol_target[i_iter+1]-frac_vol_target[i_iter]) + rhor_topdn[i_iter+1] = fsolve(func, np.array([rho_l,rho_u])) + s_curves_td[i_iter+1,:] = calc_sal_from_rho(rhor_topdn[i_iter+1]*np.ones(n_thbin),\ + th_axis,ptarget[i_iter+1]*np.ones(n_thbin)) + s_curv = s_curves_td[i_iter+1,:] + rho_l = rhor_topdn[i_iter+1] + print ' -> [iter(target) rhor_topdn(iter) rho_l rho_u]',\ + [i_iter,rhor_topdn[i_iter+1],rho_l, rho_u] + + # New attempt at computing bottom-up profile %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% print ' *** MISSING: rhor bottom-up profile ***' - print(' Computing bottom-up profile for reference density') - # Allocation - rhor_botup = np.zeros(n_target) - # Assign interpolants for the two possible reference states ------------------------------ + print ' Computing bottom-up profile for reference density' + # Allocation/Initialization + rhor_botup = np.zeros(n_target) + rhor_botup[0] = rho_min + rhor_botup[-1] = rho_max + s_curves_bu = np.zeros([n_target,n_thbin]) + s_curves_bu[0,:] = s_lower + s_curves_bu[-1,:] = s_upper + scurv = np.copy(s_upper) + rho_u = np.copy(rho_max) + # Use optimization method to converge to correct salinity curve + print ' *** this loop is not necessary ***' + for i_iter in range(n_target-1,1,-1): + rho_l = calc_rho_from_theta(sbin[0],thbin[1],ptarget[i_iter-1]) + func = lambda rho_i : vol_from_sth(rho_i, pdf_vol_lev_norm,\ + th_axis, ptarget[i_iter-1],scurv, sbin) - \ + (frac_vol_target[i_iter]-frac_vol_target[i_iter-1]) + rhor_botup[i_iter+1] = fsolve(func, np.array([rho_l,rho_u])) + s_curves_bu[i_iter-1,:] = calc_sal_from_rho(rhor_botup[i_iter-1]*np.ones(n_thbin),\ + th_axis,ptarget[i_iter-1]*np.ones(n_thbin)) + s_curv = s_curves_bu[i_iter-1,:] + rho_u = rhor_botup[i_iter-1] + print ' -> [iter(target) rhor_botup(iter) rho_l rho_u]',\ + [i_iter,rhor_botup[i_iter+1],rho_l, rho_u] + + # Assign interpolants for the two possible reference states %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% pp_rhor_botup = interpolate.PchipInterpolator(ztarget, rhor_botup, axis=0) pp_rhor_topdn = interpolate.PchipInterpolator(ztarget, rhor_topdn, axis=0) - # return interpolant for the two reference states ---------------------------------------- - #return pp_rhor_botup, pp_rhor_topdn - return pdf_vol_lev, vol_tot_ocean + + # return interpolant for the two reference states %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + return pp_rhor_botup, pp_rhor_topdn diff --git a/input_calc_ref_state.py b/input_calc_ref_state.py new file mode 100644 index 0000000..6e67d53 --- /dev/null +++ b/input_calc_ref_state.py @@ -0,0 +1,14 @@ +# Specifications for calculations ############################################################ +ref_EOS = 'jackettetal2004' # reference used to compute equ. of state + +# plotting input ############################################################################# + +# General info about plotting +dir_plot = '/glusterfs/inspect/users/xg911182/Code/Python/plots_test/' + +# variables used for lat-lon-plot and p-lat/p-lon crosssection plot +plot_var = {'rho':[]} +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 = 1 # plot crosssections +cs_plot = {'lon':[-47.5,172.5],'lat':[0]} # crosssections diff --git a/mycalc_ocean.py b/mycalc_ocean.py index 412bc51..e5d8150 100644 --- a/mycalc_ocean.py +++ b/mycalc_ocean.py @@ -1,15 +1,18 @@ -# calc_area_lat +# calc_area_xyz # calc_rho_from_theta +# calc_sal_from_rho_bis # calc_theta_from_temp # gsw_g # volume_census_levitus -def calc_area_lat(grd,r_earth,calc_old=0): - # def calc_area_lat(grd,r_earth) ######################################################### +# 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_lat + # Usage : from mycalc_ocean import calc_area_xyz # Input : grd, r_earth # grd : Input grid [python class] # r_earth : earth radius [m] @@ -140,7 +143,127 @@ def calc_rho_from_theta(s, th, p, reference='jackettetal2004'): # 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) ###################################################### @@ -158,49 +281,58 @@ def calc_theta_from_temp(s,temp,p,pr): # Output : theta # theta : potential temperature [deg C] # ######################################################################################## - # Load modules + + # Load modules %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% import numpy as np - # Error checking + + # 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 + + # Allocation and Initialization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% s1 = s*35.0/35.16504 - t1 = temp - p1 = p + t1 = np.copy(temp) + p1 = np.copy(p) pr1 = pr - # First estimate for pot. temp. - print ' Where is this equation coming from?' - theta = t1 + (p1-pr1)*( 8.65483913395442e-6 - \ + + # 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 - # 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; + + # 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 output + + # Define/return output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% return theta def gsw_g(Ns,Nt,Np,s,temp,p): @@ -449,20 +581,27 @@ def volume_census_levitus(th, s, grd, thbin, sbin): # pdf_vol_lev : # vol_ocean_lev : # ######################################################################################## - # Load modules + 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 + + # 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 + + # 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 @@ -479,16 +618,20 @@ def volume_census_levitus(th, s, grd, thbin, sbin): n_z_sub = n_z/10.0 z_range = z_max/n_z_sub vol_tot_ocean = 0.0 - print(' Vertical levels: '+ str(int(n_z)) + ', splitted into ' +\ - str(int(n_z_sub)) + ' subdomains') - # Compute the volume associted with each elements as a function of l + + # 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 = np.squeeze(calc_area_lat(grd_fine,r_earth))*res_zfine - # Split vertical domains into smaller subdomains + 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)) + 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 @@ -504,25 +647,178 @@ def volume_census_levitus(th, s, grd, thbin, sbin): 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' - s_fine[s_finesbin[1]] = np.nan - th_fine[th_finethbin[1]] = np.nan + #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] + vol_tot_ocean = vol_tot_ocean + vol_fine_lat[i_yf] del grd_fine - # Return + + # 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 diff --git a/myconstants.py b/myconstants.py index 1d24a24..0a75e21 100644 --- a/myconstants.py +++ b/myconstants.py @@ -1,5 +1,5 @@ rho0 = 1027.0 # reference density [kg m**-3] grav = 9.81 # gravity [m s**-2] -r_earth = 6.4e06 # earth radis [m] +r_earth = 6.378*1e6 # earth radis [m] gbuo = -grav/rho0 # diff --git a/myplot.py b/myplot.py index 2274619..92ddf4a 100644 --- a/myplot.py +++ b/myplot.py @@ -197,6 +197,9 @@ def myplot_2dmap(data_in,grd,FS=22,plotname='dummy.png',saveplot=0,lonlat_range= elif unit_data: cbar_label[-1] = '[' + unit_data + ']' # actual plotting + print cbar_min + print cbar_max + print cbar_range mymapf = plt.contourf(lonall, latall, data_in, cbar_range, cmap=plt.cm.Reds) #mymapf = plt.contourf(lonall, latall, data_in, 10, cmap=plt.cm.Reds, \ # vmin=cbar_min, vmax=cbar_max)