From 046e99856c4f1d4da147b00e0fdde0e28c523cc7 Mon Sep 17 00:00:00 2001 From: Gabriel Wolf Date: Fri, 21 Sep 2018 18:10:29 +0100 Subject: [PATCH] Update, including the inclusion of creating a s-th-diagram (work in progress) --- calc_ref_state.py | 52 ++++++++++++++++---- gw_ocean_refstate.py | 75 +++++++++++++++++++++++++++- mycalc_ocean.py | 113 ++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 222 insertions(+), 18 deletions(-) diff --git a/calc_ref_state.py b/calc_ref_state.py index 4da617b..8894d48 100644 --- a/calc_ref_state.py +++ b/calc_ref_state.py @@ -30,6 +30,7 @@ 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 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 @@ -44,18 +45,22 @@ from myplot_inputs import myplot_create_pn # create plotname # variables : variables necessary for calculations # MANUAL DEFINED OUTPUTS MISSING # ############################################################################################ +print '----------------------------------------------' +print 'Manual defined input data' # information about input data -print 'MANUAL DEFINED OUTPUTS MISSING' +print ' *** MANUAL DEFINED OUTPUTS MISSING ***' dim = ('lon','lat','z','time',) 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') +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 '----------------------------------------------' + # Input 1 #################################################################################### # Physical constants # input : None @@ -174,20 +179,49 @@ print '----------------------------------------------' # step e ##################################################################################### # Calculate reference state -# input : dens, zlev, area_xyz +# input : area_xyz, dens, _t, zlev, ztarget # dens : output of step c +# i_t : MISSING - DONE LOCALLY HERE # zlev : output of step a # area_xyz : output of step d -# otput : pp_rhor -# pp_rhor : interpolant of reference density +# ztarget : target depth [m], 1d - DONE LOCALLY HERE +# otput : pp_rhor, pp_rhor_botup, pp_rhor_topdn +# pp_rhor : interpolant for ref. density (hor. av. rho) +# 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) # ############################################################################################ -beg_time = time.time() +i_t = 0 +# 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 'Insert code to calculate reference state' +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') +# Deallocation +del th, s, rho print '----------------------------------------------' -end_time = time.time() -print 'Elapsed time to calculate reference state: '+ '{:.2f}'.format(end_time-beg_time)+'s' +end_time2 = time.time() +print 'Elapsed time to calculate reference state: ' +print '-> mean density (hor. averaged): '+ '{:.2f}'.format(end_time-beg_time)+'s' +print '-> sorted density in S-Theta diagram: '+ '{:.2f}'.format(end_time2-beg_time2)+'s' print '----------------------------------------------' # step f ##################################################################################### diff --git a/gw_ocean_refstate.py b/gw_ocean_refstate.py index 0041ec7..6ac71cc 100644 --- a/gw_ocean_refstate.py +++ b/gw_ocean_refstate.py @@ -1,5 +1,5 @@ def refstate_meandensity(rho_in,grd,area_xyz_in): - # def refstate_meandensity(rho,zlev,area_xyz) ############################################ + # refstate_meandensity(rho,zlev,area_xyz) ################################################ # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from get_refstate_meandensity.m of Remi Tailleux (10.04.2013) # last modified : 20.09.2018 @@ -25,3 +25,76 @@ def refstate_meandensity(rho_in,grd,area_xyz_in): pp_rhor = interpolate.PchipInterpolator(grd.z, rhor, axis=0) # return interpolant for reference state 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) ###################################### + # 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 + # 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 + # grd : Input grid [python class] + # mask : valid points [boolean], 3d + # rho : density [kg m**-3], 3d + # s : salinity [psu], 3d + # s_bin : bin range for s to construct s-th-diagram [min,max,delta] + # th : potential temperature [deg C], 3d + # th_bin : 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) + # ######################################################################################## + # 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 + 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 = + 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') + # Allocation + rhor_topdn = np.zeros(n_target) + # 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 ------------------------------ + 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 diff --git a/mycalc_ocean.py b/mycalc_ocean.py index a2139f0..412bc51 100644 --- a/mycalc_ocean.py +++ b/mycalc_ocean.py @@ -2,17 +2,19 @@ # calc_rho_from_theta # 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) ######################################################### # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from get_area_latitude.m of Remi Tailleux (09/2018) - # last modified : 19.09.2018 + # last modified : 20.09.2018 # Content/Description #################################################################### + # Usage : from mycalc_ocean import calc_area_lat # Input : grd, r_earth # grd : Input grid [python class] # r_earth : earth radius [m] # Output : area_lat - # area_lat : hor. area for 3d grid [m**2] + # area_xyz : hor. area for 3d grid [m**2] # ######################################################################################## # import modules import numpy as np @@ -20,10 +22,10 @@ def calc_area_lat(grd,r_earth,calc_old=0): deg_to_rad = np.pi/180.0 lat_radians = grd.lat*deg_to_rad if calc_old == 1: - print '*** old version of area calculation ***' + print ' *** old version of area calculation ***' # compute length of longitudinal and latitudinal elements - print '*** area calculation not valid for non-equidistant input lon-grid ***' - print '*** area calculation not valid for non-equidistant input lat-grid ***' + print ' *** area calculation not valid for non-equidistant input lon-grid ***' + print ' *** area calculation not valid for non-equidistant input lat-grid ***' dx = r_earth*deg_to_rad*(grd.Nlon/360.0) dy = r_earth*deg_to_rad*abs(grd.lat[1]-grd.lat[0]) # calculate area @@ -166,7 +168,7 @@ def calc_theta_from_temp(s,temp,p,pr): 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?' + print ' Is this above statemtent true? First estimate?' n0 = 0 n2 = 2 # Allocate values @@ -175,7 +177,7 @@ def calc_theta_from_temp(s,temp,p,pr): p1 = p pr1 = pr # First estimate for pot. temp. - print 'Where is this equation coming from?' + print ' Where is this equation coming from?' theta = t1 + (p1-pr1)*( 8.65483913395442e-6 - \ s1*1.41636299744881e-6 - \ (p1+pr1) * 7.38286467135737e-9 + \ @@ -189,7 +191,7 @@ def calc_theta_from_temp(s,temp,p,pr): # 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?' + print ' Where are the equations in this while loop coming from?' while n <= nloops: n = n+1 theta_old = theta; @@ -429,3 +431,98 @@ def gsw_g(Ns,Nt,Np,s,temp,p): (681.370187043564 - 133.5392811916956*z)*z)))) result = (g03 + g08)*1e-16 return result + +def volume_census_levitus(th, s, grd, thbin, sbin): + # volume_census_levitus(th,s,grd,thbin,sbin) ############################################# + # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk + # adapted from volume_census_levitus.m of Remi Tailleux (08.03.2013) + # last modified : 20.09.2018 + # Content/Description #################################################################### + # Usage : from mycalc_ocean import volume_census_levitus + # Input : grd, s, sbin, th, thbin + # grd : Input grid [python class] + # s : salinity [psu], 3d + # sbin : bin range for s to construct s-th-diagram [min,max,delta] + # th : potential temperature [deg C], 3d + # thbin : bin range for th to construct s-th-diagram [min,max,delta] + # Output : pdf_vol_lev, vol_ocean_lev + # pdf_vol_lev : + # vol_ocean_lev : + # ######################################################################################## + # Load modules + import numpy as np + from myconstants import r_earth + from mydata_classes import Grid # class containing grid information + from scipy.interpolate import interpn # interpolation in 3d + # error checking + if grd.symlon == 0: + raise ValueError('*** Input grid must span full longitude grid ***') + # 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 + res_xfine = 0.25 + res_yfine = 0.25 + res_zfine = 10.0 + if grd.lon[1]-grd.lon[0] > res_xfine: + lon_fine = np.arange(res_xfine/2.0,360,res_xfine) + else: + lon_fine = grd.lon + if grd.lat[1]-grd.lat[0] > res_yfine: + lat_fine = np.arange(-90.0+res_yfine/2.0,90.0,res_yfine) + else: + lat_fine = grd.lat + z_max = max(grd.z)+max(grd.z)%(10**np.floor(np.log10(max(grd.z)))) + n_z = np.ceil(z_max/res_zfine) + n_z_sub = n_z/10.0 + z_range = z_max/n_z_sub + vol_tot_ocean = 0.0 + 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 + 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 + for i_subdom in range(0,int(n_z_sub)): + print ' Subdomain ' + str(int(i_subdom+1)) + ' of ' + str(int(n_z_sub)) + # create meshgrid for interpolation + z_beg = res_zfine/2.0 + i_subdom * z_range + z_end = res_zfine/2.0 + (i_subdom+1) * z_range - 1 + z_fine = np.arange(z_beg,z_end,res_zfine) + yfine_xyz, xfine_xyz, zfine_xyz = np.meshgrid(lat_fine, lon_fine, z_fine) + # interpolate s and th onto finer grid + reshape_form = [len(lon_fine),len(lat_fine),len(z_fine)] + xyz_fine_pt = np.array([np.reshape(xfine_xyz,[np.prod(reshape_form)]),\ + np.reshape(yfine_xyz,[np.prod(reshape_form)]),\ + np.reshape(zfine_xyz,[np.prod(reshape_form)])]).T + s_fine = np.reshape(interpn((grd.lon,grd.lat,grd.z), s, xyz_fine_pt, \ + bounds_error=False, fill_value=None),reshape_form) + th_fine = np.reshape(interpn((grd.lon,grd.lat,grd.z), th, xyz_fine_pt, \ + bounds_error=False, fill_value=None),reshape_form) + del xyz_fine_pt + print ' replace this np.nan stuff' + 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 + 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): + pdf_vol_lev[i_sd,i_thd] = pdf_vol_lev[i_sd,i_thd] + \ + vol_fine_lat[i_yf] + vol_tot_ocean = vol_tot_ocean + vol_fine_lat[i_yf] + del grd_fine + # Return + return pdf_vol_lev, vol_tot_ocean +