#!/usr/bin/env/python # The curv2curvipcc.py CDAT program can be used to transform variables # on a curvilinear grid to quasi-IPCC output format files: there will be # one output variable per CF compliant output file and the attributes # (coordinates', variable's and file's attributes) will be those # required by IPCC, but the variable will still be on a curvilinear grid # (whereas IPCC requires the output variables to be defined on # rectilinear grids). # # The program is currently set up to create monthly means (1 year with # 12 time steps) from variables on the ORCA 4 curvilinear grid (OPA # ocean model output), but it should be easy to have it work with other # curvilinear grids and a different number of time steps. # # Input: # # * File(s) with variables on 12 time steps and a curvilinear grid # (and vertical levels). # # Note: the script does not perform any time averaging. The 12 time # steps in the input files must already have the correct values # (e.g. the monthly means must have been computed with an external # program, or by using NCO operators). # # * File(s) with coordinates' bounds. # # Note: the script can be configured to look for the required # coordinates (and coordinates' bounds) in several files, including the # input variables files themselves). # # Output: # # * IPCC like output files (1 variable per file) # # Note: see also http://www-lsce.cea.fr/motif/contrib/utils/cdat/curv2curvipcc.shtml # # Jean-Yves Peterschmitt - LSCE - 09/2004 import sys, os from time import asctime import Numeric, MA, cdms, cdtime from cdms.coord import TransientVirtualAxis, TransientAxis2D from cdms.hgrid import TransientCurveGrid ################################################################### # Parameters ################################################################### # Experiment name #expname = '21K_OA' expname = '0K_OA' # List of the variables that have to be generated (1 variable / file) #v_out_list = ['thetao', 'so', 'kzo', # 'uox', 'voy', 'wo', # 'uoeiv', 'voeiv', 'woeiv'] v_out_list = ['thetao'] # Global attributes of the output files, required by IPCC # For more details, see 'Requirements for global attributes' in # http://www-pcmdi.llnl.gov/ipcc/IPCC_output_requirements.pdf # 'institution' global attribute globatt_inst = "LSCE (Laboratoire des Sciences du Climat et de l'Environnement, Gif-sur-Yvette, France)" # 'source' global attribute globatt_src = 'IPSLCM4_v1 (2004, ipsl_cm4_v1_9_1): atmosphere: LMDZ3.3 (IPSL-CM4_IPCC, 72x45x19); ocean: OPA8.2/ORCA4 (ipsl_cm4_v1_8, 92x76x31); seaice: LIM (ipsl_cm4_v1_8); land: ORCHIDEE/SECHIBA (orchidee_1_3)' # 'project_id' global attribute globatt_project = 'MOTIF' # 'table_id' global attribute template globatt_tblid_tpl = 'Table %s' # 'realization' global attribute globatt_realnum = 1 # 'contact' global attribute globatt_contact = 'Jean-Yves Peterschmitt - motifweb@lsce.saclay.cea.fr' # Dictionary of the attributes that depend on the experiment name attri_dic = {} expid = 'PMIP2_0K_OA' attri_dic['0K_OA'] = {'expid':expid, 'longexpid':'%s Pre-industrial Ocean-Atmosphere PMIP2 experiment' \ % (expid,), 'exp_name':'BR03', 'exp_years':'0051_100', 'comment':'Experiment length: 50 years (51 to 100)', 'references':'No references'} expid = 'PMIP2_21K_OA' attri_dic['21K_OA'] = {'expid':expid, 'longexpid':'%s Last Glacial Maximum Ocean-Atmosphere PMIP2 experiment' \ % (expid,), 'exp_name':'LGM04', 'exp_years':'101_110', 'comment':'Initial state of ocean temperature and salinity is 1600-09-01 of a LGM simulation of Hadley Centre. Experiment length: 10 years (101 to 110)', 'references':'No references'} # Set the attributes that depend on the experiment name if not attri_dic.has_key(expname): # Abort the program if the current experiment is not defined... raise 'Warning! The experiment %s is not defined' % (expname,) # 'experiment_id' global attribute expid = attri_dic[expname]['expid'] globatt_longexpid = attri_dic[expname]['longexpid'] # Local experiment name part of the 'history' global attribute exp_name = attri_dic[expname]['exp_name'] globatt_localexp = 'LSCE experiment name: %s' % (exp_name,) # Note: the exp_years variable will be used later to determine # the name of the input file exp_years = attri_dic[expname]['exp_years'] # 'comment' global attribute globatt_comment = attri_dic[expname]['comment'] # 'references' global attribute globatt_references = attri_dic[expname]['references'] # 'title' global attribute globatt_title = 'IPSL model output prepared for the MOTIF %s' % (globatt_longexpid,) # Name of the IPCC/MOTIF configuration table where the output # variables are (or should be!) described # Note: the 'table_id' global attribute of the output file will be based on this t_name = 'O_SE' # Output directory (no trailing '/' needed) f_out_dir = '.' # Output file name (including path) template # Note: the actual file name will use the name of the output variable f_out_name_tpl = '%s/%%s_%s_orcagrid_LSCE_%s.nc' % (f_out_dir, t_name, expid) # Dictionary of input variables' locations # (we use the original names in the file) v_locdic = {} #commondir = 'experiments/%s' % (exp_name,) commondir = '.' v_locdic['vomecrty'] = '%s/%s_SE_%s_grid_V.nc' % (commondir, exp_name, exp_years) v_locdic['vomeeivv'] = '%s/%s_SE_%s_grid_V.nc' % (commondir, exp_name, exp_years) v_locdic['vosaline'] = '%s/%s_SE_%s_grid_T.nc' % (commondir, exp_name, exp_years) v_locdic['votemper'] = '%s/%s_SE_%s_grid_T.nc' % (commondir, exp_name, exp_years) v_locdic['votkeavt'] = '%s/%s_SE_%s_grid_W.nc' % (commondir, exp_name, exp_years) v_locdic['vovecrtz'] = '%s/%s_SE_%s_grid_W.nc' % (commondir, exp_name, exp_years) v_locdic['voveeivw'] = '%s/%s_SE_%s_grid_W.nc' % (commondir, exp_name, exp_years) v_locdic['vozocrtx'] = '%s/%s_SE_%s_grid_U.nc' % (commondir, exp_name, exp_years) v_locdic['vozoeivu'] = '%s/%s_SE_%s_grid_U.nc' % (commondir, exp_name, exp_years) # List of the expected (i.e. 'allowed') input variable shapes # Note: the current script only allows one shape, but it should be # possible to deal with a varying number of levels (i.e. deal with the # 'no vertical levels at all' case) and maybe a varying number of time steps v_in_shapelist = [(12, 31, 76, 92)] # Dictionary of output variables features # Note units conversion and other variable alterations takes place later # (look for 'units conversion') v_featdic = {} v_featdic['thetao'] = {'in_name':'votemper', 'in_name_full':'votemper + 273.15', 'out_units':'K', 'out_stdname':'sea_water_potential_temperature', 'out_longname':'Ocean potential temperature', 'out_comment':''} v_featdic['so'] = {'in_name':'vosaline', 'in_name_full':'vosaline', 'out_units':'1e-3', 'out_stdname':'sea_water_salinity', 'out_longname':'Ocean salinity', 'out_comment':''} v_featdic['uox'] = {'in_name':'vozocrtx', 'in_name_full':'vozocrtx', 'out_units':'m s-1', 'out_stdname':'sea_water_x_velocity', 'out_longname':'Ocean velocity along X grid axis', 'out_comment':''} v_featdic['voy'] = {'in_name':'vomecrty', 'in_name_full':'vomecrty', 'out_units':'m s-1', 'out_stdname':'sea_water_y_velocity', 'out_longname':'Ocean velocity along Y grid axis', 'out_comment':''} v_featdic['wo'] = {'in_name':'vovecrtz', 'in_name_full':'vovecrtz', 'out_units':'m s-1', 'out_stdname':'upward_sea_water_velocity', 'out_longname':'Vertical current velocity', 'out_comment':''} v_featdic['uoeiv'] = {'in_name':'vozoeivu', 'in_name_full':'vozoeivu', 'out_units':'m s-1', 'out_stdname':'eddy induced velocity along x', 'out_longname':'Ocean EIV along x grid axis', 'out_comment':''} v_featdic['voeiv'] = {'in_name':'vomeeivv', 'in_name_full':'vomeeivv', 'out_units':'m s-1', 'out_stdname':'eddy induced velocity along y', 'out_longname':'Ocean EIV along x grid axis', 'out_comment':''} v_featdic['woeiv'] = {'in_name':'voveeivw', 'in_name_full':'voveeivw', 'out_units':'m s-1', 'out_stdname':'eddy induced vertical velocity', 'out_longname':'Vertical ocean EIV', 'out_comment':''} v_featdic['kzo'] = {'in_name':'votkeavt', 'in_name_full':'votkeavt', 'out_units':'m2 s-1', 'out_stdname':'vertical turbulent diffusivity coefficient', 'out_longname':'Vertical turbulent diffusivity coefficient', 'out_comment':''} # Value that will be used for the v_out.cell_methods attribute # v_out_cellmth value should be IMPROVED!! v_out_cellmth = 'time: mean' # Value of the 'missing value' in the output files v_out_missing = 1e20 # Axes dictionary # For each type of axis (and bounds), the dictionary specifies # what axis or variable should be used, and from which file # The program will use the first axes/variable in the 'names' list # it finds in the the 'files' list of files axes_dic = {} # Note: the first file in the list of files where we look for the axes # is the first input file we use first_in = v_locdic[v_featdic[v_out_list[0]]['in_name']] axes_dic['z'] = {'names':['deptht', 'depth'], 'files':[first_in, './extra_orca4grid_info.nc']} axes_dic['y'] = {'names':['nav_lat', 'lat'], 'files':[first_in, './extra_orca4grid_info.nc']} axes_dic['x'] = {'names':['nav_lon', 'lon'], 'files':[first_in, './extra_orca4grid_info.nc']} axes_dic['z_bnds'] = {'names':['bounds_depth', 'depth_bnds'], 'files':[first_in, './extra_orca4grid_info.nc']} axes_dic['y_bnds'] = {'names':['bounds_lat', 'lat_bnds'], 'files':[first_in, './extra_orca4grid_info.nc']} axes_dic['x_bnds'] = {'names':['bounds_lon', 'lon_bnds'], 'files':[first_in, './extra_orca4grid_info.nc']} # Origin of the time axis time_units = 'days since 0001-1-1' ################################################################### # Miscellaneous checks ################################################################### if not os.path.exists(f_out_dir) or not os.path.isdir(f_out_dir): raise 'Output directory %s does not exist, or is not a directory' % (f_out_dir,) ################################################################### # Creation of the axes ################################################################### # Function that tries to find a given axis data in the axes' dictionary def getAxisData(ad_name, ad_dic): axis_dic = ad_dic[ad_name] axis_files = axis_dic['files'] axis_names = axis_dic['names'] axis_found = 0 # Get the first available axis in the first available file for af_name in axis_files: if not os.path.exists(af_name): # File not found, try next file continue af = cdms.open(af_name) for an in axis_names: if af.axes.has_key(an): axis_data = af.axes[an][:] axis_found = 1 # Exit the axis name loop break if an in af.listvariables(): axis_data = af(an) axis_found = 1 # Exit the axis name loop break # Close the file af.close() # Exit the file loop, if the axis has been found if axis_found: break if axis_found: # Print info mmessage and return the axis data print ' | Axis data [ %-10s ] found as [ %12s ] in file\n |\t%s' % (ad_name, an, af_name) return axis_data else: # Exit with an error message raise ' | Axis data [%6s] not found as %s in files\n |\t%s' % (ad_name, str(axis_names), str(axis_files)) print 'Creating the axes for the output variables' sys.stdout.flush() # Create the 'time' time axis time_vals = Numeric.arange(15, 360, 30, Numeric.Float64) time = cdms.createAxis(time_vals, id='time') time.designateTime(calendar=cdtime.Calendar360) time.units = time_units time.standard_name = 'time' time.long_name = 'time' time.bounds = 'time_bnds' # Create a dimension that won't be associated with a variable # (i.e. it will only appear in the 'dimensions' section of a netcdf # file, not in the 'variables' section) bnds_axis = TransientVirtualAxis("bnds", 2) # Create the time_bnds variable (bounds of the time axis) time_bnds_vals = Numeric.zeros((12, 2), Numeric.Float64) time_bnds_vals[:, 0] = Numeric.arange(0, 360, 30, Numeric.Float64) time_bnds_vals[:, 1] = Numeric.arange(30, 361, 30, Numeric.Float64) time_bnds = cdms.createVariable(time_bnds_vals, axes=(time, bnds_axis), id='time_bnds') del time_bnds.name # Remove 'name' attribute created by cdms # Create the 'depth' axis # Note: we force the typecode of 'depth' to be Numeric.Float64 # Note: we assume that the depth is increasing, and positive down... depth_vals = getAxisData('z', axes_dic) depth = cdms.createAxis(depth_vals.astype(Numeric.Float64), id='depth') depth.designateLevel() depth.units = 'm' depth.standard_name = 'depth' depth.long_name = 'depth' depth.bounds = 'depth_bnds' depth.positive = 'down' # Get the bounds of the 'depth' axis and create a 'clean' depth_bnds variable # Note: we force the typecode of 'depth_bnds' to be Numeric.Float64 depth_bnds_vals = getAxisData('z_bnds', axes_dic) depth_bnds = cdms.createVariable(depth_bnds_vals.filled(), axes=(depth, bnds_axis), typecode=Numeric.Float64, id='depth_bnds') del depth_bnds.name # Remove 'name' attribute created by cdms # Create the dimensions for the curvilinear variables # We get the grid shape from the last 2 values of the first # tuple in the list of accepted shapes s_in = v_in_shapelist[0][-2:] y_axis = TransientVirtualAxis("y", s_in[-2]) x_axis = TransientVirtualAxis("x", s_in[-1]) nvert_axis = TransientVirtualAxis("nvert", 4) # Get the bounds of the 'lat' axis and create a 'clean' lat_bnds variable # Note: we force the typecode of 'lat_bnds' to be Numeric.Float64 lat_bnds_vals = getAxisData('y_bnds', axes_dic) lat_bnds = cdms.createVariable(lat_bnds_vals.filled(), axes=(y_axis, x_axis, nvert_axis), typecode=Numeric.Float64, id='lat_bnds') del lat_bnds.name # Remove 'name' attribute created by cdms # Create the curvilinear 'lat' latitude axis # Note: we force the typecode of 'lat' to be Numeric.Float64 lat_vals = getAxisData('y', axes_dic) lat = TransientAxis2D(lat_vals.filled(), #bounds=lat_bnds, axes=(y_axis, x_axis), typecode=Numeric.Float64, id='lat') lat.standard_name = 'latitude' lat.long_name = 'latitude' lat.units= 'degrees_north' lat.axis = 'Y' lat.bounds = 'lat_bnds' del lat.name # Remove 'name' attribute created by cdms # Create the curvilinear 'lon' longitude axis in # the same fashion as 'lat' lon_bnds_vals = getAxisData('x_bnds', axes_dic) lon_bnds = cdms.createVariable(lon_bnds_vals.filled(), axes=(y_axis, x_axis, nvert_axis), typecode=Numeric.Float64, id='lon_bnds') del lon_bnds.name # Remove 'name' attribute created by cdms lon_vals = getAxisData('x', axes_dic) lon = TransientAxis2D(lon_vals.filled(), #bounds=lon_bnds, axes=(y_axis, x_axis), typecode=Numeric.Float64, id='lon') lon.standard_name = 'longitude' lon.long_name = 'longitude' lon.units= 'degrees_east' lon.axis = 'X' lon.bounds = 'lon_bnds' del lon.name # Remove 'name' attribute created by cdms # Create the curvilinear grid associated with lat and lon grid_out = TransientCurveGrid(lat, lon, id='clean_curv_grid') ################################################################### # Creation of the variables ################################################################### # Loop on the requested output variables for v_out_name in v_out_list: print 'Processing output variable [%s]' % (v_out_name,) sys.stdout.flush() # Store the date/time when we start processing the variable # (this will be used in the 'history' attribute(s) timestamp = asctime() # Get the details' dictionary of the output variable out_dic = v_featdic[v_out_name] # Get the all the details associated with the output variable v_in_name = out_dic['in_name'] f_in_name = v_locdic[v_in_name] print ' | [%-15s] is based on [%15s] in file\n |\t%s' % (v_out_name, v_in_name, f_in_name) sys.stdout.flush() v_in_name_full = out_dic['in_name_full'] v_out_units = out_dic['out_units'] v_out_stdname = out_dic['out_stdname'] v_out_longname = out_dic['out_longname'] v_out_comment = out_dic['out_comment'] f_out_name = f_out_name_tpl % (v_out_name,) # Read the input variable f_in = cdms.open(f_in_name) v_in = f_in(v_in_name) f_in.close() # Save the units v_in_units = v_in.units # Change the input variable to a simple masked array (MV -> MA) # This is just to make sure that we won't inherit attributes/axes # of the old variable v_in = MA.masked_array(v_in) # Set the spacesaver attribute to 1, in case we have to # alter the values of the input variable later, to be sure no # type upgrading (Float32 -> Float64) will occur v_in.savespace(1) # Note, the missing value will be adjusted just before creating # the output variable, AFTER possible adjustment to the variable # (that may reset the missing value to its default value) # Check the shape of the input variable s_in = v_in.shape if s_in not in v_in_shapelist: raise ' ! Wrong shape %s for the input variable!' % (str(s_in),) # Patch the value of some input variables # (e.g. take care of the units conversion) if v_in_name == 'votemper': print ' | Converting units: C --> K' sys.stdout.flush() v_in = v_in + 273.15 # Set the value of the new (maybe) missing value # Note: there may be points of the original data that have the # value of the new missing value (which would be stupid, but # who knows) so we mask those points, before setting the # missing value... v_in = MA.masked_values(v_in, v_out_missing) v_out_missing = Numeric.array(v_out_missing, Numeric.Float32) v_in.set_fill_value(v_out_missing) # Create the output variable on the curvilinear grid # and make sure that the variable is of Float32 type! v_out = cdms.createVariable(v_in, typecode=Numeric.Float32, axes=(time, depth) + grid_out.getAxisList(), grid=grid_out, id=v_out_name) # Add IPCC requested attributes to the variable v_out.units = v_out_units v_out.standard_name = v_out_stdname # Note: cdms does not want to set '_FillValue'... v_out._FillValue = v_out_missing v_out.cell_methods = v_out_cellmth v_out.original_name = v_in_name_full v_out.history = '%s: Variable based on %s from file %s' % (timestamp, v_in_name, f_in_name) v_out.original_units = v_in_units v_out.long_name = v_out_longname v_out.comment = v_out_comment # Open the output file print ' | Writing variable to\n |\t%s... ' % (f_out_name,), sys.stdout.flush() f_out = cdms.open(f_out_name, 'w') # Set the global attributes of the file f_out.institution = globatt_inst f_out.source = globatt_src f_out.project_id = globatt_project f_out.table_id = globatt_tblid_tpl % (t_name,) f_out.realization = globatt_realnum f_out.experiment_id = globatt_longexpid f_out.contact = globatt_contact hline = '%s: File created by %s with %s' % (timestamp, os.getlogin(), sys.argv[0]) f_out.history = '%s\n%s' % (hline, globatt_localexp) f_out.comment = globatt_comment f_out.references = globatt_references f_out.title = globatt_title # Write the curvilinear axes' bounds and the variable to # the output file f_out.write(time_bnds) f_out.write(depth_bnds) f_out.write(lat_bnds) f_out.write(lon_bnds) # Note: the attributes below should not have to be set # and it is impossible to set _FillValue !! f_out.write(v_out, attributes=v_out.attributes) # Close the output file f_out.close() print '[done]' sys.stdout.flush() # End of the loop on requested output variables # The end