local hops; /* DOCUMENT hops.i --- Hyperslab operator suite List of routines * hadd_attr: add hyperslab attribute * happend: append hyperslabs to netCDF file * hatmepv: compute atmospheric potential vorticity * hatmepflux: compute atmospheric E-P flux * hattr: get hyperslab attribute * hbin: carry out bin-averaging on dimension * hcat: concatenate hyperslabs * hclose: close netCDF file * hcombine: combine hysperslabs from several files into a single file * hconform: check conformance of two or more hyperslabs * hcont: superimpose continental outlines * hcoord: return coordinate values * hcopy: copy hyperslab * hdata: read actual data into hyperslab * hdiff: difference two netCDF files/hyperslab arrays * hdimsof: return dimension counts for hyperslab components * hfold: fold time dimension to create index dimension * hfft: fourier transform selected dimension * hgather: combine several hyperslabs into a single array * hget: read variable form history file as hyperslab * hgrow: grow hyperslab array * hinterp: interpolate to new grid * hlegend: return descriptive legend string about hyperslab * hmask: apply region masking on hyperslab * hocnmask: superimpose land/depth mask (for ocean data) * hocndrho: compute gradients of sea water density * hocnrho: compute density of sea water * hop: carry out arithmetic/logical operations on hyperslabs * hopen: open netCDF file for reading hyperslabs * hplot: display hyperslab data graphically * hregrid: regrid data along dimension * hsave: save hyperslabs in netCDF file * hset_attr: set hyperslab attribute * hshift: shift data along dimension * hshop: carry out spherical harmonic operations * hshtran: carry out spherical harmonic transforms * hsplit: split all slices along a dimension as separate hyperslabs * hsprout: re-introduce eliminated dimensions * hsub: reduce hyperslab domain through slicing, averaging etc. * htbin: carry out bin-averaging on time-series * htwrite: write time-series to file * hunfold: unfold index dimension into time dimension * hvecplot: overlay vectors on plot * hver_wt: compute vertical weights for 3-D hyperslab * * Version: 1.5 * Modified: 26 Jan 1999, R. Saravanan * Web page: http://www.cgd.ucar.edu/gds/svn/hops * * SEE ALSO: tops, flexplot, ncops, yodel, lowops */ local lowops; /* DOCUMENT lowops --- Low-level functions for hyperslab operator suite List of routines * natmfile: open CSM format atmospheric netCDF file * nattr: return coordinate attributes of hyperslab * nattlist: return structure containing all attributes of a variable * nattmatch: check if attributes of two variables match * nbinop: carry out binary operation between two data arrays * ncat: concatenate hyperslabs along a dimension * ncheck_grid: check consistency of staggered grid * ncopyatt: copy variable attributes from one slab to another * ndataprec: return data precision, given data type * ndim_bounds: return subdomain bounds * ndimlist: read variable dimensions from netCDF file * ngetatt: get variable attributes from netCDF file * ngetcoord: get coordinate values from hyperslab * nget_handle: extract file handle from file structure * nhyperfile: open hyperslab netCDF file * ninit: initialize HOPS * ninterp: auxiliary routine for HINTERP * nlocate: return locator structure for variable in netCDF file * nmask: carry out masking operation on an array * nmiss_value: generate/convert missing value * nocnfile: open CSM format atmospheric netCDF file * nocnrho: initialize sea water density polynomial coefficients * nplotcontours: determine contour levels/labels for ocean fields * nplotrange: determine axis coordinate ranges for ocean fields * npsum: compute partial sums along dimension * nputatt: write variable attributes to netCDF file * nrange: select subrange of hyperslab dimension * nrecord: write record variables to netCDF file * nreduce: carry out rank-reduction operation on hyperslab dimension * nrotate: rotate periodic X dimension of hyperslab * nsave: write non-record variables to netCDF file * nset_attr: set coordinate attribute of hyperslab * nshiftmask: shift 2D bitmask from regular to interfacial grid * nsize_vec: return hyperslab size vector * nslabarr: access an effective hyperslab array * nsublegend: return string describing hyperslab subdomain * nunop: carry out unary operation on hyperslab * nvarattdef: define variable and its attributes in a netCDF file * * SEE ALSO: hops, ncops, yodel, flexplot */ require, "ncops.i"; require, "yodel.i"; require, "flexplot.i"; struct hfmt_struc{ // Hyperslab format descriptor structure string hopsroot, dimnames(5), dimintnames(5), dimaltnames(5); string coordnames(5), reduceops(6); long data, area_wt, z_bot; pointer dimlist, varlist, attdata, attglob, attfixed; pointer ignorelist, attdesc; pointer ocn_dimlist, ocn_varlist, ocn_attdesc; pointer atm_dimlist, atm_varlist, atm_attdesc; pointer ssh_dimlist, ssh_varlist, ssh_attdesc; long nattstd(3); double epsdate, epscoord, default_missing_value; long nchar, nsigma_coefs; pointer shtab; long ishtab; } struct hyperfile_struc{ // Yorick notation for hyperslab netCDF file descriptor string structure; // Structure name ("HYPERFILE") NC_file fmeta; // File structure descriptor long fnumber; // File number string fname; // File name string vars; // File variable list string std_dims(5); // Names of the standard dimensions string int_dims(5); // Names of the interfacial dimensions pointer other_dims; // pointer to a 1-D string array pointer other_labels; // pointer to a 1-D string array // (labels are stored concatenated together in one string, separated by ";") long phys_offset(5,2); // Physical offset for regular/int grid long reverse(5); // Coordinate reversal flag long recordvars; // No. of record variables (-1,0,1,...) long scratch; // Scratch file flag pointer x0, y0, z0; // pointers to 1-D double arrays pointer xint0, yint0, zint0; // pointers to 1-D double arrays pointer time0, date0; // pointers to 1-D double arrays pointer ilabel0; // pointer to a 1-D string array pointer iparam0; // pointer to a 1-D double array pointer template; // Hyperslab template } struct ocnfile_struc{ // Yorick notation for ocean history file descriptor string structure; // Structure name ("ATMFILE") NC_file fmeta; // File structure descriptor long fnumber; // File number string fname; // File name string vars; // File variable list string ftype; // "PROGNOSTIC" / "DIAGNOSTIC" string fconventions; // "CSM" / "HYPERSLAB" string std_dims(5); // Names of the standard dimensions string int_dims(5); // Names of the interfacial dimensions pointer other_dims; // pointer to a 1-D string array pointer other_labels; // pointer to a 1-D string array // (labels are stored concatenated together in one string, separated by ";") long phys_offset(5,2); // Physical offset for regular/int grid long reverse(5); // Coordinate reversal flag long recordvars; // No. of record variables (-1,0,1,...) long scratch; // Scratch file flag string area_wt_units; // Initial area weight units string area_wt_var, z_bot_var; // Area weight/Z_BOT variable names string area_wtint_var, z_botint_var; // Interfacial area weight/Z_BOT variable names pointer x0, y0, z0; // pointers to 1-D double arrays pointer xint0, yint0, zint0; // pointers to 1-D double arrays pointer time0, date0; // pointers to 1-D double arrays pointer ilabel0; // pointer to a 1-D string array pointer iparam0; // pointer to a 1-D double array pointer sigma0, sigmaint0; // pointers to 2-D double arrays pointer area_wt0, area_wtint0; // pointers to 2-D double arrays pointer z_bot0, z_botint0; // pointers to 2-D double arrays double f0, g0; // polar coriolis parameter, // surface gravity double rhoocn0, cpocn0, socn0; // density, specific heat, mean salinity // By convention, the first "extra" dimension corresponds to the // list of pre-defined horizontal regions, if any long nrmask; // No. of pre-defined horizontal regions pointer rmask; // Horizontal region mask // (2-D byte array) pointer template; // Hyperslab template } struct atmfile_struc{// Yorick notation for atmospheric history file descriptor string structure; // Structure name ("ATMFILE") NC_file fmeta; // File structure descriptor long fnumber; // File number string fname; // File name string vars; // File variable list string ftype; // "" string fconventions; // "CSM" / "HYPERSLAB" string std_dims(5); // Names of the standard dimensions string int_dims(5); // Names of the interfacial dimensions pointer other_dims; // pointer to a 1-D string array pointer other_labels; // pointer to a 1-D string array // (labels are stored concatenated together in one string, seperated by ";") long phys_offset(5,2); // Physical offset for regular/int grid long reverse(5); // Coordinate reversal flag long recordvars; // No. of record variables (-1,0,1,...) long scratch; // Scratch file flag string area_wt_units; // Initial area weight units string area_wt_var, z_bot_var; // Area weight/Z_BOT variable names pointer x0, y0, z0; // pointers to 1-D double arrays pointer zint0; // pointers to 1-D double arrays pointer time0, date0; // pointers to 1-D double arrays pointer ilabel0; // pointer to a 1-D string array pointer iparam0; // pointer to a 1-D double array pointer sigma0, sigmaint0; // pointers to 2-D double arrays pointer area_wt0; // pointer to 1-D double array pointer z_bot0; // pointers to 2-D double arrays double f0, g0; // polar coriolis parameter, // surface gravity double rdry0, cpdry0; // dry air constant, specific heat // By convention, the first "extra" dimension corresponds to the // list of pre-defined horizontal regions, if any long nrmask; // No. of pre-defined horizontal regions pointer rmask; // Horizontal region mask // (2-D byte array) pointer template; // Hyperslab template } struct locator_struc { // Yorick notation for slab data locator structure string structure; // structure name ("slabdata") pointer fstruc; // pointer to file structure string fname; // file name NC_file fmeta; // File structure descriptor float add_offset, scale_factor; // data offset, scale factor string area_wt_var, z_bot_var; // Area weight/Z_BOT variable names string type(3); // DATA/AREA_WT/Z_BOT type string string dimenstr; // dimension string long dim_data(5); // data dimension index number in file long dim_area_wt(5); // area_wt dimension index number long dim_z_bot(5); // z_bot index number long was_present(5); // Initial dimension presence code long slab_wrapcount; // hyperslab X wrap-around count pointer slab_fold; // hyperslab T fold descriptor pointer slab_offset, slab_count; // hyperslab offset/count } struct hyperslab { // Yorick notation for a hyperslab // Mandatory variables describing the hyperslab // (A coordinate variable may be null, if the corresponding dimension // was not present even in the original history file, i.e., prior to slicing, // rank-reduction operations.) // SDIM=5: no. of standard dimensions // NATT=56: total number of attributes // NIATT=7, NFATT=8, NSATT=41: numberof integer, float, string attributes string structure; // Structure name ("HYPERSLAB...") pointer x, y, z, time; // pointer to a 1-D double arrays pointer ilabel; // pointer to a 1-D string array pointer data; // pointer to a 5-D float/double array pointer missing_value; // pointer to a scalar of same type as // the data string name, long_name, units; pointer attlist; // Attribute list("var","nam","var:nam") pointer attcode; // Attribute codes (type(1-3),index) pointer iatt; // Integer attributes: pointer fatt; // Float attributes: pointer satt; // String attributes: string type(3); // DATA/AREA_WT/Z_BOT type string long dimension(5,3); // DATA/AREA_WT/Z_BOT dimensions long reduced(5); // DATA rank-reduction codes // Optional variables describing the hyperslab (may be null) pointer area_wt, z_bot; // pointers to a 5-D float/double arrays pointer date; // pointer to a 1-D double array pointer iparam; // pointer to a 1-D double array // Optional variables describing the full spatial domain grid pointer x0, y0, z0; // pointers to 1-D double arrays pointer xint0, yint0, zint0; // pointers to 1-D double arrays pointer ilabel0; // pointer to a 1-D string array pointer iparam0; // pointer to a 1-D double array } struct hyperslab_atm { // Yorick notation for a ATM-type hyperslab // Mandatory variables describing the hyperslab // (A coordinate variable may be null, if the corresponding dimension // was not present even in the original history file, i.e., prior to slicing, // rank-reduction operations.) // SDIM=5: no. of standard dimensions // NATT=56: total number of attributes // NIATT=7, NFATT=8, NSATT=41: numberof integer, float, string attributes string structure; // Structure name ("HYPERSLAB...") pointer x, y, z, time; // pointer to a 1-D double arrays pointer ilabel; // pointer to a 1-D string array pointer data; // pointer to a 5-D float/double array pointer missing_value; // pointer to a scalar of same type as // the data string name, long_name, units; pointer attlist; // Attribute list("var","nam","var:nam") pointer attcode; // Attribute codes (type(1-3),index) pointer iatt; // Integer attributes: pointer fatt; // Float attributes: pointer satt; // String attributes: string type(3); // DATA/AREA_WT/Z_BOT type string long dimension(5,3); // DATA/AREA_WT/Z_BOT dimensions long reduced(5); // DATA rank-reduction codes // Optional variables describing the hyperslab (may be null) pointer area_wt, z_bot; // pointers to a 5-D float/double arrays pointer date; // pointer to a 1-D double array pointer iparam; // pointer to a 1-D double array // Optional variables describing the full spatial domain grid pointer x0, y0, z0; // pointers to 1-D double arrays pointer xint0, yint0, zint0; // pointers to 1-D double arrays pointer ilabel0; // pointer to a 1-D string array pointer iparam0; // pointer to a 1-D double array // SPH extension fields double a0; // planetary radius pointer eqdx0, cosdy0; // pointer to 1-D double arrays pointer eqdxint0, cosdyint0; // pointer to 1-D double arrays // SIG extension fields pointer sigma0, sigmaint0; // pointers to 2-D double arrays // ATM extension fields pointer hgrid0; // pointer to 2-D byte array } struct hyperslab_ocn { // Yorick notation for a OCN-type hyperslab // Mandatory variables describing the hyperslab // (A coordinate variable may be null, if the corresponding dimension // was not present even in the original history file, i.e., prior to slicing, // rank-reduction operations.) // SDIM=5: no. of standard dimensions // NATT=56: total number of attributes // NIATT=7, NFATT=8, NSATT=41: numberof integer, float, string attributes string structure; // Structure name ("HYPERSLAB...") pointer x, y, z, time; // pointer to a 1-D double arrays pointer ilabel; // pointer to a 1-D string array pointer data; // pointer to a 5-D float/double array pointer missing_value; // pointer to a scalar of same type as // the data string name, long_name, units; pointer attlist; // Attribute list("var","nam","var:nam") pointer attcode; // Attribute codes (type(1-3),index) pointer iatt; // Integer attributes: pointer fatt; // Float attributes: pointer satt; // String attributes: string type(3); // DATA/AREA_WT/Z_BOT type string long dimension(5,3); // DATA/AREA_WT/Z_BOT dimensions long reduced(5); // DATA rank-reduction codes // Optional variables describing the hyperslab (may be null) pointer area_wt, z_bot; // pointers to a 5-D float/double arrays pointer date; // pointer to a 1-D double array pointer iparam; // pointer to a 1-D double array // Optional variables describing the full spatial domain grid pointer x0, y0, z0; // pointers to 1-D double arrays pointer xint0, yint0, zint0; // pointers to 1-D double arrays pointer ilabel0; // pointer to a 1-D string array pointer iparam0; // pointer to a 1-D double array // SPH extension fields double a0; // planetary radius pointer eqdx0, cosdy0; // pointer to 1-D double arrays pointer eqdxint0, cosdyint0; // pointer to 1-D double arrays // SIG extension fields pointer sigma0, sigmaint0; // pointers to 2-D double arrays // OCN extension fields pointer kmax0; // pointer to 2-D byte array } struct hyperslab_ssh { // Yorick notation for a SSH-type hyperslab // Mandatory variables describing the hyperslab // (A coordinate variable may be null, if the corresponding dimension // was not present even in the original history file, i.e., prior to slicing, // rank-reduction operations.) // SDIM=5: no. of standard dimensions // NATT=56: total number of attributes // NIATT=7, NFATT=8, NSATT=41: numberof integer, float, string attributes string structure; // Structure name ("HYPERSLAB...") pointer x, y, z, time; // pointer to a 1-D double arrays pointer ilabel; // pointer to a 1-D string array pointer data; // pointer to a 5-D float/double array pointer missing_value; // pointer to a scalar of same type as // the data string name, long_name, units; pointer attlist; // Attribute list("var","nam","var:nam") pointer attcode; // Attribute codes (type(1-3),index) pointer iatt; // Integer attributes: pointer fatt; // Float attributes: pointer satt; // String attributes: string type(3); // DATA/AREA_WT/Z_BOT type string long dimension(5,3); // DATA/AREA_WT/Z_BOT dimensions long reduced(5); // DATA rank-reduction codes // Optional variables describing the hyperslab (may be null) pointer area_wt, z_bot; // pointers to a 5-D float/double arrays pointer date; // pointer to a 1-D double array pointer iparam; // pointer to a 1-D double array // Optional variables describing the full spatial domain grid pointer x0, y0, z0; // pointers to 1-D double arrays pointer xint0, yint0, zint0; // pointers to 1-D double arrays pointer ilabel0; // pointer to a 1-D string array pointer iparam0; // pointer to a 1-D double array // SSH extension fields double a0; // planetary radius // SIG extension fields pointer sigma0, sigmaint0; // pointers to 2-D double arrays } struct attlist_struc { // Attribute list for a variable pointer inames, ivals; // Integer-valued attributes pointer fnames, fvals; // Float-valued attributes pointer snames, svals; // String-valued attributes } struct shtran_struc { // Grid info for spherical harmonic transforms long nlon, nlat, m, n, k; // Truncation parameters pointer x0, y0; // pointer to 1-D double arrays pointer eqdx0, cosdy0; // pointer to 1-D double arrays pointer area_wt0; // pointer to 2-D double array pointer hgrid0; // pointer to 2-D byte array pointer init; // pointer to Shtran_Init structure } func hadd_attr(slab,attribute,value,help=) /* DOCUMENT hadd_attr(slab,attribute,value,help=help) * Adds specified ATTRIBUTE to SLAB, sets it to VALUE, and returns new slab, * where ATTRIBUTE="varname:attribute_name", or ATTRIBUTE=":attribute_name" * for global attributes. * SEE ALSO: hset_attr, hattr, hcopy */ { func_name= "hadd_attr"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; if (param_set(help)) { write,""; write, "Function HADD_ATTR adds specified ATTRIBUTE to SLAB, && sets it"; write, "to VALUE, where ATTRIBUTE='varname:attribute_name', || "; write, " ATTRIBUTE=':attribute_name' for global attributes. "; write," E.g.,"; write," new_slab= hadd_attr(slab, 'data:new_attr', 'ATT_VALUE')"; write," adds the attribute 'data:new_attr', sets it to 'ATT_VALUE', &&"; write," returns the modified slab."; write," See also: hset_attr, hattr, hcopy"; write,""; write," Usage: hadd_attr(slab, 'varname:attribute', value)"; return timer_return(func_name, NULL); } if (typeof(slab) != "struct_instance") error, "Argument SLAB should be a structure"; if (!is_scalar(slab)) error, "Can only add attributes to scalar slabs"; // Copy attribute list/codes attlist= *(slab.attlist); attcode= *(slab.attcode); // Locate attribute attwhere= where(attlist(I0+2,) == attribute); // Split attribute into variable name and attribute name str2= strsplit(attribute, ":"); if (numberof(str2) != 2) error, "Invalid attribute specification - " + attribute; attvar= str2(I0); attname= str2(I0+1); new_slab= NULL; if ( ((attvar == "data") && \ is_where(where(*(HFMT.attdata) == attname))) || \ (attname == "dimension") || \ (attribute == ":structure") || \ is_where(attwhere) ) { // Attribute already exists; simply copy slab hcopy, slab, new_slab; if (is_where(attwhere)) { attno= attwhere(I0); if (attcode(I0+1,attno) < 0) { // "Undelete" attribute attcode(I0+1,attno)= -attcode(I0+1,attno); new_slab.attcode= ref(attcode); } } } else { // Copy slab, adding attribute attype= typeof(value); if (attype == "int") attype= "long"; if (attype == "float") attype= "double"; if ((attype != "long") && (attype != "double") && \ (attype != "string")) error, "Invalid attribute type - " + attype; hcopy, slab, new_slab, extra_atts=[attvar, attname, attype]; } // Set attribute value hset_attr, new_slab, attribute, value; return timer_return(func_name, new_slab); } func happend( &fstruc, //YORICKoutput: slab0, slab1, slab2, slab3, slab4, slab5, slab6, slab7, slab8, slab9, help=, overwrite=, partial=, nocheck=, create=) /* DOCUMENT happend, fstruc, * slab0, slab1, slab2, slab3, slab4, * slab5, slab6, slab7, slab8, slab9, * help=0/1, overwrite=0/1, partial=0/1, nocheck=0/1, * create= * * Append upto 9 hyperslabs (or hyperslab arrays) in the t-dimension * to the netCDF file described by file data structure FSTRUC, * updating FSTRUC to reflect the extended t-dimension. * The appended hyperslabs must have strong full conformance in all * dimensions, including reduced ones, and also variable, case, and domain * conformance, with the hyperslabs already present in the file * (except for the t-dimension, where only dimension-unit conformance * and monotonically increasing t-coordinate values are required). * If OVERWRITE==1, overwrite existing data. * If PARTIAL==1, allow partial records to be written. * If NOCHECK==1, no slab conformance checking is done, and area weights/Z_BOT * values are always written out (use this option with care). * If CREATE="filename" is specified and FSTRUC has a null value, * a new netCDF file is created and its file structure is returned as * FSTRUC. Subsequent calls to HAPPEND can append to this file. * SEE ALSO: hsave, hopen, hconform */ { func_name= "happend"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; if (param_set(help)) { write,""; write," Function HAPPEND appends a set of hyperslabs to a netCDF file"; write," created using HSAVE, && subsequently opened using HOPEN. E.g.,"; write," happend, fstruc, slab1, slab2"; write," appends to the file pointed to by file structure FSTRUC"; write," On returning, the FSTRUC is modified to reflect the appended data."; write," (The hyperslabs SLAB1 && SLAB2, which may themselves be arrays.)"; write," Tips:"; write," 1. overwrite=1 allows overwriting of existing records."; write," 2. partial=1 allows writing of partial records."; write," 3. nocheck=1 suppresses conformance checking."; write," 4. create='filename' creates a new netCDF if FSTRUC is null."; write,""; write," See also: hsave, hopen, hconform"; write,""; write," Usage: happend,fstruc,slab1,slab2,...,overwrite=0/1,partial=0/1,nocheck=0/1,create='filename'"; return timer_return(func_name); } if (is_null(fstruc)) { if (is_null(create)) error, "Null value for FSTRUC; specify create='filename'" // Create new save file, open it for appending, and return hsave, create, slab0, slab1, slab2, slab3, slab4, slab5, slab6, slab7, slab8, slab9, nocheck=nocheck; hopen, create, fstruc, alt=1, append=1, silent=1; return timer_return(func_name); } if (typeof(fstruc) != "struct_instance") error, "Argument FSTRUC should be a structure"; if (fstruc.structure != "HYPERFILE") error, "Cannot append to file of type " + fstruc.structure; if (fstruc.recordvars == 0) error, "Cannot append to file without unlimited (record) dimension"; // Get file handle fhandle= nget_handle( fstruc ); // List of variables in file filevars= strsplit(fstruc.vars, ","); // Prior time and date values time1= *(fstruc.time0); date1= deref(fstruc.date0); nt1= numberof(time1); // Number of slabs nslab= numberof(slab0) + numberof(slab1) + numberof(slab2) + numberof(slab3) + numberof(slab4) + numberof(slab5) + numberof(slab6) + numberof(slab7) + numberof(slab8) + numberof(slab9); if (nslab == 0) error, "No slabs specified to append"; // T/I dimension transpose flag ti_transp= 0; // Flags/dimensions for area weights, and Z_BOT values area_wt_flag= 0; area_wt_var= NULL; apresent0= NULL; along_name= NULL; aunits= NULL; aelements= NULL; area_wt1= NULL; z_bot_flag= 0; z_bot_var= NULL; zpresent0= NULL; zlong_name= NULL; zunits= NULL; zref= NULL; z_bot1= NULL; // Variable list newvars= array("",nslab); for (jslab=I0; jslab <= nslab-I1; jslab++) { // Get copy of each new slab (with data) tem_slab= nslabarr( jslab+I1, slab0, slab1, slab2, slab3, slab4, slab5, slab6, slab7, slab8, slab9 ); // Variable name varname= tem_slab.name; // Check for uniqueness of variable name if (strloc(newvars, varname) > 0) error, "Duplicate variable name in hyperslab list - " + varname; is_present= tem_slab.dimension(,HFMT.data); if (is_present(TDIM) <= 0) error, "T dimension not present in variable " + varname; if (is_present(IDIM) > 0) { // Set T/I dimension transpose flag ti_transp= 1; } // Add variable name to list newvars(jslab)= varname; if (tem_slab.type(HFMT.data) == "struct_instance") error, "Slab "+tem_slab.name+" does not contain actual data"; if (!param_set(nocheck)) { // Get corresponding old hyperslab from file (without data) old_slab= hget(varname, fstruc=fstruc, nodata=1); old_loc= *(old_slab.data); // Check conformance of old and new slabs (including reduced dimensions) isuperset= (dim_conf= (udim_conf= (unit_conf= NULL))); var_conf= (case_conf= (grid_conf= (domain_conf= NULL))); hconform, old_slab, tem_slab, isuperset, dim_conf, udim_conf, unit_conf, var_conf, case_conf, grid_conf, domain_conf, reduced=1; if (!allof(dim_conf(XDIM:ZDIM) == 2)) error, "XYZ-dimensions not conforming for variable " + varname; if (dim_conf(IDIM) != 2) error, "I-dimension not conforming for variable " + varname; if (!udim_conf(TDIM)) error, "T-dimension not unit conforming for variable " + varname; if ( (!var_conf) || (!case_conf) || (!domain_conf) ) error, "Variable/case/domain non-conformance for variable " + varname; if (is_null(old_slab.missing_value)) { if (!is_null(tem_slab.missing_value)) error, "Missing value attribute mismatch for variable " + varname; } else { if (!is_null(tem_slab.missing_value)) { if (*(old_slab.missing_value) != *(tem_slab.missing_value)) error, "Missing value attribute mismatch for variable " + varname; } } } if (jslab == I0) { // First slab time2= *(tem_slab.time); date2= deref(tem_slab.date); nt2= numberof(time2); if (is_null(date1) != is_null(date2)) error, "Date value presence mismatch"; if (time2(I0) > time1(nt1-I1)) { // Append time values overwrite_flag= 0; ioffset= nt1; time0= time1; date0= date1; grow, time0, time2; grow, date0, date2; // Create new slab file data structure with extended time coordinate fstruc2= hyperfile_struc( fmeta= fstruc.fmeta, other_dims= fstruc.other_dims, other_labels= fstruc.other_labels, x0= fstruc.x0, y0= fstruc.y0, z0= fstruc.z0, xint0= fstruc.xint0, yint0= fstruc.yint0, zint0= fstruc.zint0, time0= ref(time0), date0= ref(date0), ilabel0= fstruc.ilabel0, iparam0= fstruc.iparam0, template= fstruc.template ); fstruc2.structure= "HYPERFILE"; // Initialize fixed size fields fstruc2.fnumber= fstruc.fnumber; fstruc2.fname= fstruc.fname; fstruc2.vars= fstruc.vars; fstruc2.std_dims= fstruc.std_dims; fstruc2.int_dims= fstruc.int_dims; fstruc2.phys_offset= fstruc.phys_offset; fstruc2.recordvars= fstruc.recordvars; fstruc2.reverse= fstruc.reverse; } else { // Overwrite values if (!param_set(overwrite)) error, "Overlapping time values; specify overwrite=1 to overwrite records" overwrite_flag= 1; // Locate new time coordinate offset among old values ioffset= rangeloc( time1, time2(I0) ) - 1; if ( (ioffset+nt2) > nt1 ) error, "Time values to be overwritten are not fully overlapping"; if (!array_eq(time1(I0+ioffset:I0+ioffset+nt2-1), time2, epsilon=HFMT.epscoord)) error, "Time values to be overwritten are not contiguous"; if (!is_null(date1)) { if (!array_eq(date1(I0+ioffset:I0+ioffset+nt2-1), date2, epsilon=HFMT.epsdate/max([max(date2),1.0]) )) error, "Date values to be overwritten do not match"; } } } else { // Not the first slab; check time/date coordinate values if (!array_eq(time2,*(tem_slab.time),epsilon=HFMT.epscoord)) error, "Time values not the same for all slabs"; if (!is_null(date2)) { if (!array_eq(date2, deref(tem_slab.date), epsilon=HFMT.epsdate/max([max(date2),1.0]) )) error, "Date values not the same for all slabs"; } } if (tem_slab.type(HFMT.area_wt) != "") { // Area weights present if ((!area_wt_flag) && (!param_set(nocheck))) { // First set of area weights; check for variable name area_wt_flag= 1; area_wt_var= old_loc.area_wt_var; if (strloc(filevars, area_wt_var) == 0) area_wt_var= NULL; } if (!is_null(area_wt_var)) { // Area weights variable attributes if (is_null(apresent0)) { apresent0= tem_slab.dimension(,HFMT.area_wt); along_name= hattr( tem_slab, "area_wt:long_name" ); aunits= hattr( tem_slab, "area_wt:units" ); aelements= hattr( tem_slab, "area_wt:elements" ); area_wt1= *(tem_slab.area_wt); } else { if ( (!array_eq( apresent0, tem_slab.dimension(,HFMT.area_wt)) ) || \ (along_name != hattr(tem_slab, "area_wt:long_name")) || \ (aunits != hattr(tem_slab, "area_wt:units")) || \ (aelements != hattr(tem_slab, "area_wt:elements")) ) error, "Area weights dimensions/attributes are not the same"; if (!array_eq(area_wt1, *(tem_slab.area_wt), epsilon=HFMT.epscoord)) error, "Area weights array values do not match among all slabs"; } } } if (tem_slab.type(HFMT.z_bot) != "") { // Z_BOT values present if ((!z_bot_flag) && (!param_set(nocheck))) { // First set of Z_BOT values; check for variable name z_bot_flag= 1; z_bot_var= old_loc.z_bot_var; if (strloc(filevars, z_bot_var) == 0) z_bot_var= NULL; } if (!is_null(z_bot_var)) { // Z_BOT variable attributes if (is_null(zpresent0)) { zpresent0= tem_slab.dimension(,HFMT.z_bot); zlong_name= hattr( tem_slab, "z_bot:long_name" ); zunits= hattr( tem_slab, "z_bot:units" ); zref= hattr( tem_slab, "z_bot:ref" ); z_bot1= *(tem_slab.z_bot); } else { if ( (!array_eq( zpresent0, tem_slab.dimension(,HFMT.z_bot)) ) || \ (zlong_name != hattr(tem_slab, "z_bot:long_name")) || \ (zunits != hattr(tem_slab, "z_bot:units")) || \ (zref != hattr(tem_slab, "z_bot:ref")) ) error, "Z_BOT dimensions/attributes are not the same"; if (!array_eq(z_bot1, *(tem_slab.z_bot), epsilon=HFMT.epscoord)) error, "Z_BOT values do not match among all slabs"; } } } } // Area weights/Z_BOT array write flags awrite= 1; zwrite= 1; if (!is_null(area_wt_var)) { // Area weights variable iloc= strloc(newvars, area_wt_var); if (iloc == 0) error, "Area weights variable not found among slabs"; tem_slab= nslabarr( iloc, slab0, slab1, slab2, slab3, slab4, slab5, slab6, slab7, slab8, slab9 ); awrite= 0; aelements1= hattr(tem_slab, "area_wt:elements"); if ( (!array_eq( apresent0, tem_slab.dimension(,HFMT.data)) ) || \ (along_name != hattr(tem_slab, "area_wt:long_name")) || \ (aunits != hattr(tem_slab, "area_wt:units")) || \ (aelements != aelements1) ) error, "Area weights variable dimensions/attributes do not match"; if (!array_eq(area_wt1, *(tem_slab.data), epsilon=HFMT.epscoord)) error, "Area weights array values do not match variable values"; } if (!is_null(z_bot_var)) { // Z_BOT variable iloc= strloc(newvars, area_wt_var); if (iloc == 0) error, "Z_BOT variable not found among slabs"; tem_slab= nslabarr( iloc, slab0, slab1, slab2, slab3, slab4, slab5, slab6, slab7, slab8, slab9 ); zwrite= 0; zref1= hattr(tem_slab, "z_bot:ref"); if ( (!array_eq( zpresent0, tem_slab.dimension(,HFMT.data)) ) || \ (zlong_name != hattr(tem_slab, "data:long_name")) || \ (zunits != hattr(tem_slab, "data:units")) || \ (zref != zref1) ) error, "Z_BOT variable dimensions/attributes do not match"; if (!array_eq(z_bot1, *(tem_slab.data), epsilon=HFMT.epscoord)) error, "Z_BOT values do not match variable values"; } if (!param_set(partial)) { // Check to make sure that a full record is being written dimnames= HFMT.dimnames; for (ivar=I0; ivar <= numberof(filevars)-I1; ivar++) { if (strloc(newvars, filevars(ivar)) == 0) { vardims= nc_getdims( fstruc.fmeta, filevars(ivar) ); if (strloc(vardims, dimnames(TDIM)) > 0) error, "Set partial=1 to write partial record"; } } } // Write record variables for (k=1; k <= nt2; k++) { if (!overwrite_flag) { //YORICKbegin: // Add records for time dimension nc_addrec, fhandle, time2(k-I1) //YORICKend: // Write time values nc_putvar, fhandle, "time", [time2(k-I1)], offset=[ioffset+k-1], record=1; // Write date values if (!is_null(date2)) nc_putvar, fhandle, "date", [date2(k-I1)], offset=[ioffset+k-1], record=1; } // Write record variables nrecord, fhandle, ioffset+k, slab0, k, 1, ti_transp=ti_transp, area_wt=awrite, z_bot=zwrite; nrecord, fhandle, ioffset+k, slab1, k, 1, ti_transp=ti_transp, area_wt=awrite, z_bot=zwrite; nrecord, fhandle, ioffset+k, slab2, k, 1, ti_transp=ti_transp, area_wt=awrite, z_bot=zwrite; nrecord, fhandle, ioffset+k, slab3, k, 1, ti_transp=ti_transp, area_wt=awrite, z_bot=zwrite; nrecord, fhandle, ioffset+k, slab4, k, 1, ti_transp=ti_transp, area_wt=awrite, z_bot=zwrite; nrecord, fhandle, ioffset+k, slab5, k, 1, ti_transp=ti_transp, area_wt=awrite, z_bot=zwrite; nrecord, fhandle, ioffset+k, slab6, k, 1, ti_transp=ti_transp, area_wt=awrite, z_bot=zwrite; nrecord, fhandle, ioffset+k, slab7, k, 1, ti_transp=ti_transp, area_wt=awrite, z_bot=zwrite; nrecord, fhandle, ioffset+k, slab8, k, 1, ti_transp=ti_transp, area_wt=awrite, z_bot=zwrite; nrecord, fhandle, ioffset+k, slab9, k, 1, ti_transp=ti_transp, area_wt=awrite, z_bot=zwrite; } if (!overwrite_flag) { // Copy updated file structure fstruc= fstruc2; } return timer_return(func_name); } func hattr( slab, attribute, help=, index=) /* DOCUMENT hattr(slab, attribute, help=, index=) * Returns the specified ATTRIBUTE of SLAB, * where ATTRIBUTE="varname:attribute_name", or ATTRIBUTE=":attribute_name" * for global attributes. If the attribute is not defined, a null value is * returned. * If ATTRIBUTE="varname:" is specified, all attributes of" * of the variable are printed out, and a null string is returned. * (The "units" and "long_name" attributes for the five standard dimensions * may be accessed through the array components nattr("units",slab), * and nattr("long_name",slab). The "name", "long_name", "units", * and "missing_value" attributes for the data may accessed directly as * structure members slab.*) * If INDEX is non-null, return attributes for SLAB(INDEX) * If SLAB may be an array of hyperslabs, an array of attribute values * is returned. * SEE ALSO: hset_attr, hadd_attr, hcopy, hget */ { func_name= "hattr"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; if (param_set(help)) { write,""; write,"Function HATTR returns the specified ATTRIBUTE of SLAB,"; write," where ATTRIBUTE='varname:attribute_name', || "; write," ATTRIBUTE=':attribute_name' for global attributes. "; write," E.g.,"; write," case_name = hattr(slab, ':case_name')"; write," returns the global string attribute 'case_name'."; write," If SLAB is actually an array of hyperslabs, case_name would"; write," be an array of strings."; write," If ATTRIBUTE='varname:' is specified, all attributes of"; write," of variable 'varname' are printed out."; write," E.g.,"; write," write,hattr(slab, 'time:')"; write," prints out all attributes of variable 'time'."; write," See also: hset_attr, hadd_attr, hcopy, hget"; write,""; write," Usage: hattr(slab,'varname:attribute',index=...)"; return timer_return(func_name, ""); } if (typeof(slab) != "struct_instance") error, "Argument SLAB should be a structure"; if (is_null(index) && (!is_scalar(slab))) { // Array of hyperslabs; handle recursively att_array= NULL; for (j=I0; j <= I0+numberof(slab)-1; j++) { tem_att= NULL; tem_att= hattr( slab(j), attribute ); grow, att_array, tem_att; } return timer_return(func_name, att_array); } k1= I0; if (!is_null(index)) k1= index; // Copy attribute list/codes attlist= *(slab(k1).attlist); attcode= *(slab(k1).attcode); // Locate attribute attwhere= where(attlist(I0+2,) == attribute); if (is_where(attwhere)) { attno= attwhere(I0); // If "deleted" attribute, return null if (attcode(I0+1,attno) <= 0) return timer_return(func_name, NULL); // Attribute type/index attype= attcode(I0,attno); attinx= attcode(I0+1,attno) - I1; if (attype == 1) { return timer_return(func_name, (*(slab(k1).iatt))(attinx)); } else if (attype == 2) { return timer_return(func_name, (*(slab(k1).fatt))(attinx)); } else if (attype == 3) { return timer_return(func_name, (*(slab(k1).satt))(attinx)); } } // Special handling for attributes accessible as structure members if (attribute == ":structure") { return timer_return(func_name, slab(k1).structure); } nlen= strlen(attribute); if (nlen > 5 ) { if ( (strmid(attribute,0,5) == "data:") && \ anyof(*(HFMT.attdata) == strmid(attribute,5,nlen-5)) ) { // Data attribute accessible as structure member attname= strmid(attribute,5,nlen-5); if (attname == "missing_value") { return timer_return(func_name, deref(slab(k1).missing_value)); } else if (attname == "name") { return timer_return(func_name, slab(k1).name); } else if (attname == "long_name") { return timer_return(func_name, slab(k1).long_name); } else if (attname == "units") { return timer_return(func_name, slab(k1).units); } else { error, "Internal error 1"; } } } if (strmid(attribute,nlen-1,1) == ":") { // Recursively display all attributes of variable and return null string attvar= strmid(attribute,0,nlen-1); attwhere= where(attlist(I0,) == attvar); if (is_where(attwhere)) { attsel= attlist(I0+2,attwhere) ; for (j=I0; j <= numberof(attsel)-I1; j++) { write, attsel(j)+" = ", hattr(slab, attsel(j), index=index); } } return timer_return(func_name, ""); } // Attribute not found; return NULL value return timer_return(func_name, NULL); } func hatmepflux( t_slab, u_slab, v_slab, omega_slab, &epfy_slab, &epfz_slab, &epfdc_slab, //YORICKoutput: help=,fpolar=,pref=,gravit=,kappa=, nohistory=) /* DOCUMENT hatmepflux, t_slab, u_slab, v_slab, omega_slab, * epfy_slab, epfz_slab, epfdc_slab, * fpolar=,pref=,gravit=,kappa=, * nohistory=0/1 * Returns Eliassen-Palm flux components in pressure coordinates- * EPFY_SLAB: meridional component of E-P flux (in m2 s-2) * EPFZ_SLAB: vertical component of E-P flux (in m Pa s-2) * EPFDC_SLAB: cosine-latitude times E-P flux divergence (in m s-2) * given the temperatures (T_SLAB, in K), horizontal velocity components * (U_SLAB, V_SLAB, in m/s), and pressure tendency (OMEGA_SLAB, in Pa/s) * * FPOLAR is the polar value of the Coriolis parameter (in s-1) * (defaults to 1.4584E-04) * * PREF is the reference pressure (in Pa) used to compute potential * temperature (defaults to 1000.0e2) * * GRAVIT is the surface gravitational acceleration (in m/s2) * (defaults to 9.80616) * * KAPPA is tha dimensionless ratio R/C_p (defaults to 287.04/1004.64) * * SEE ALSO: hatmepv */ { func_name= "hatmepflux"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; //IDLbegin: //:error, "Spectral transforms module SHTRAN not yet implemented in IDL"; //IDLend: //YORICKbegin: require, "shtran.i" //YORICKend: if (param_set(help)) { write,""; write," Procedure HATMEPFLUX returns the Eliassen-Palm flux components,"; write," given T, U, V, && OMEGA."; write," E.g.,"; write," hatmepflux, t_slab, u_slab, v_slab, omega_slab,"; write," epfy_slab, epfz_slab, epfdc_slab"; write," See also: hatmepv"; write,""; write," Usage: hatmepflux, t_slab, u_slab, v_slab, omega_slab, epfy_slab, epfz_slab, epfdc_slab"; return timer_return(func_name); } if (is_null(t_slab) || is_null(u_slab) || \ is_null(v_slab) || is_null(omega_slab)) error, "Null operand(s)"; if ( (typeof(t_slab) != "struct_instance") || \ (typeof(u_slab) != "struct_instance") || \ (typeof(v_slab) != "struct_instance") || \ (typeof(omega_slab) != "struct_instance") ) error, "Operands not hyperslabs"; if ((!is_scalar(t_slab)) || \ (!is_scalar(u_slab)) || \ (!is_scalar(v_slab)) || \ (!is_scalar(omega_slab)) ) error, "Operands should be scalar slabs"; if ( (t_slab.structure != "HYPERSLAB1.0_SPH_SIG_ATM") || \ (u_slab.structure != "HYPERSLAB1.0_SPH_SIG_ATM") || \ (v_slab.structure != "HYPERSLAB1.0_SPH_SIG_ATM") || \ (omega_slab.structure != "HYPERSLAB1.0_SPH_SIG_ATM") ) error, "Incorrect slab structure for operands"; if ((t_slab.units != "K") && (t_slab.units != "Kelvin")) error, "Unrecognized units for T_SLAB - " + t_slab.units; if ((u_slab.units != "m/s") && (u_slab.units != "m s-1")) error, "Unrecognized units for U_SLAB - " + u_slab.units; if ((v_slab.units != "m/s") && (v_slab.units != "m s-1")) error, "Unrecognized units for V_SLAB - " + v_slab.units; if ((omega_slab.units != "Pa/s") && (omega_slab.units != "Pa s-1")) error, "Unrecognized units for OMEGA_SLAB - " + omega_slab.units; his_str= "<"+ t_slab.name+ ">,<" + u_slab.name + ">" + ">,<" + v_slab.name+ ">,<" + omega_slab.name +">"; f0= double(1.4584E-04); if (!is_null(fpolar)) { his_str= his_str + ",fpolar=" + strnum(fpolar); f0= fpolar; } pref1= 1000.0e2; if (!is_null(pref)) { his_str= his_str + ",pref=" + strnum(pref); pref1= pref; } gravit1= 9.80616; if (!is_null(gravit1)) { his_str= his_str + ",gravit=" + strnum(gravit); gravit1= gravit; } kappa1= 287.04/1004.64; if (!is_null(kappa1)) { his_str= his_str + ",kappa=" + strnum(kappa); kappa1= kappa; } // Compute sine/cosine of latitude SINLAT= hop("sin", hcoord(t_slab,"y"), name="SINLAT"); COSLAT= hop("cos", hcoord(t_slab,"y"), name="COSLAT"); // Compute planetary vorticity PLVOR= hop(SINLAT, "*", f0, name="PLVOR", units="s-1"); // Compute zonal-mean absolute vorticity RELVOR= hshtran(hgather(u_slab,v_slab),vec=1,curl=1,phys=1,preserve=1); ABSVORXAV= hop(PLVOR, "+", hsub(RELVOR,x="avg"), units="s-1"); RELVOR= NULL; // Get pressure values on pressure levels (Pa) PLEV= hcoord(t_slab,"z"); if (PLEV.units == "hPa") PLEV= hop(PLEV, "*", 100., units="Pa"); if (PLEV.units != "Pa") error, "Incorrect pressure units"; // Compute delta pressure DELTAP2= hop( hshift(PLEV,"z",count=1), "-", hshift(PLEV,"z",count=-1), name="DELTAP2"); // Compute potential temperature THETA= hop( t_slab, "*", hop( hop(pref1,"/",PLEV,units=""), "^", kappa1) , name="THETA"); // E-P flux components (in pressure coordinates) // // ( [...] denotes zonal-averaging, (") denotes deviation from zonal-mean) // // PSIind = cos(phi) [v"theta"] / (d[theta]/dp) // // EPFY = PSIind d[u]/dp - [v"u"] cos(phi) // // EPFZ = (f + [xi]) PSIind - [omega"u"] cos(phi) // Compute vertical derivative of zonal-mean THETA THETAXAV= hsub(THETA,x="avg"); DTHDPXAV= hop( hop( hshift(THETAXAV,"z",count=1), "-", hshift(THETAXAV,"z",count=-1)), "/", DELTAP2, name="DTHDPXAV"); THETAXAV= NULL; // Compute vertical derivative of zonal-mean U UXAV= hsub(U,x="avg"); DUDPXAV= hop( hop( hshift(UXAV,"z",count=1), "-", hshift(UXAV,"z",count=-1)), "/", DELTAP2, name="DUDPXAV"); UXAV= NULL; // Compute [v"theta"], [u"v"], [u"omega"] VP= hop(v_slab, "-", hsub(v_slab,x="avg",nohistory=1),nohistory=1); THP= hop(THETA, "-", hsub(THETA,x="avg")); VTHPXAV= hsub( hop(VP,"*",THP,name="VTHP", units="m K s-1"), x="avg", nohistory=1 ); THP= NULL; THETA= NULL; UP= hop(u_slab, "-", hsub(u_slab,x="avg")); UVPXAV= hsub( hop(UP,"*",VP, name="UVP", units="m2 s-2"), x="avg" ); VP= NULL; OMP= hop(omega_slab, "-", hsub(omega_slab,x="avg")); UOMPXAV= hsub( hop(UP,"*",OMP,name="UOMP", units="m Pa s-2"), x="avg" ); UP= NULL; OMP= NULL; // Compute induced streamfunction PSIIND= hop( hop(VTHPXAV,"*",COSLAT,nohistory=1), "/", DTHDPXAV, name="PSIIND", units="m Pa s-1", nohistory=1); // Compute E-P flux components epfy_slab= hop( hop(PSIIND,"*",DUDPXAV, units="m2 s-2", nohistory=1), "-", hop(COSLAT,"*",UVPXAV), name="EPFY", nohistory=1 ); epfz_slab= hop( hop(PSIIND,"*",ABSVORXAV, units="m Pa s-2", nohistory=1), "-", hop(COSLAT,"*",UOMPXAV), name="EPFZ", nohistory=1 ); epfy_slab.long_name= "E-P_flux_y_component"; epfz_slab.long_name= "E-P_flux_p_component"; UVPXAV= NULL; UOMPXAV= NULL; VTHPXAV= NULL; DUDPXAV= NULL; ABSVORXAV= NULL; // Broadcast cos(phi)*EPFY to have original X dimension EPFYC= hop(COSLAT,"*",epfy_slab); EPFYCB= NULL; hcopy, t_slab, EPFYCB, data=broadcast(*(EPFYC.data),hdimsof(t_slab)); EPFYCB.name= "EPFYC"; // Compute gradient of cos(phi)*EPFY EPFYCGRAD= hshtran(EPFYCB,grad=1,phys=1,preserve=1); EPFYCY= hsub(EPFYCGRAD(I0+1),x="avg",name="EPFYCY"); EPFYCY.units= "m s-2"; EPFYC= NULL; EPFYCB= NULL; EPFYCGRAD= NULL; // Compute vertical derivative of EPFZ DEPFZDP= hop( hop( hshift(epfz_slab,"z",count=1,nohistory=1), "-",hshift(epfz_slab,"z",count=-1,nohistory=1)), "/", DELTAP2, name="DEPFZDP", units="m s-2", nohistory=1); // Compute E-P flux divergence epfdc_slab= hop( hop(DEPFZDP,"*",COSLAT,nohistory=1), "+", EPFYCY, name="EPFDC", nohistory=1 ); epfdc_slab.long_name= "E-P_flux_divergence_cosine"; EPFYCY= NULL; DEPFZDP= NULL; if (!param_set(nohistory)) { // Append history info to output slabs hset_attr, epfy_slab, "data:history", hattr(epfy_slab,"data:history") + " hatmepflux, " + his_str + ", ...;" hset_attr, epfz_slab, "data:history", hattr(epfz_slab,"data:history") + " hatmepflux, " + his_str + ", ...;" hset_attr, epfdc_slab, "data:history", hattr(epfdc_slab,"data:history") + " hatmepflux, " + his_str + ", ...;" } return timer_return(func_name); } func hatmepv( t_slab, u_slab, v_slab, help=, fpolar=, pref=, gravit=, kappa=, tref=, nohistory=) /* DOCUMENT hatmepv(t_slab, u_slab, v_slab, fpolar=, * pref=, gravit=, kappa=, * tref=, nohistory=0/1) * Computes Ertel Potential Vorticity (EPV, in m2 K s-1 kg-1), * given the temperatures (T_SLAB, in K) and velocity components * (U_SLAB, V_SLAB, in m/s). * * FPOLAR is the polar value of the Coriolis parameter (in s-1) * (defaults to 1.4584E-04) * * PREF is the reference pressure (in Pa) used to compute potential * temperature (defaults to 1000.0e2) * * GRAVIT is the surface gravitational acceleration (in m/s2) * (defaults to 9.80616) * * KAPPA is tha dimensionless ratio R/C_p (defaults to 287.04/1004.64) * * If reference temperature TREF (in K) is specified, the PV is * scaled to remove the static potential temperature dependence, * and the resulting nondimensional PV is returned. For isothermal * barotropic flow, the scaled PV should show no height variation. * * SEE ALSO: hatmepflux */ { func_name= "hatmepv"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; //IDLbegin: //:error, "Spectral transforms module SHTRAN not yet implemented in IDL"; //IDLend: //YORICKbegin: require, "shtran.i" //YORICKend: if (param_set(help)) { write,""; write," Function HATMEPV computes the Ertel Potential Vorticity,"; write," given T, U, && V."; write," E.g.,"; write," pv_slab = hatmepv(t_slab, u_slab, v_slab)"; write," See also: hatmepflux"; write,""; write," Usage: hatmepv(t_slab, u_slab, v_slab, ...)"; return timer_return(func_name, NULL); } if (is_null(t_slab) || is_null(u_slab) || is_null(v_slab)) error, "Null operand(s)"; if ( (typeof(t_slab) != "struct_instance") || \ (typeof(u_slab) != "struct_instance") || \ (typeof(v_slab) != "struct_instance") ) error, "Operands not hyperslabs"; if ((!is_scalar(t_slab)) || \ (!is_scalar(u_slab)) || \ (!is_scalar(v_slab)) ) error, "Operands should be scalar slabs"; if ( (t_slab.structure != "HYPERSLAB1.0_SPH_SIG_ATM") || \ (u_slab.structure != "HYPERSLAB1.0_SPH_SIG_ATM") || \ (v_slab.structure != "HYPERSLAB1.0_SPH_SIG_ATM") ) error, "Incorrect slab structure for operands"; if ((t_slab.units != "K") && (t_slab.units != "Kelvin")) error, "Unrecognized units for T_SLAB - " + t_slab.units; if ((u_slab.units != "m/s") && (u_slab.units != "m s-1")) error, "Unrecognized units for U_SLAB - " + u_slab.units; if ((v_slab.units != "m/s") && (v_slab.units != "m s-1")) error, "Unrecognized units for V_SLAB - " + v_slab.units; his_str= "<"+ t_slab.name+ ">,<" + u_slab.name + ">" + ">,<" + v_slab.name +">"; f0= double(1.4584E-04); if (!is_null(fpolar)) { his_str= his_str + ",fpolar=" + strnum(fpolar); f0= fpolar; } pref1= 1000.0e2; if (!is_null(pref)) { his_str= his_str + ",pref=" + strnum(pref); pref1= pref; } gravit1= 9.80616; if (!is_null(gravit)) { his_str= his_str + ",gravit=" + strnum(gravit); gravit1= gravit; } kappa1= 287.04/1004.64; if (!is_null(kappa)) { his_str= his_str + ",kappa=" + strnum(kappa); kappa1= kappa; } if (param_set(tref)) { // Prepare to scale PV his_str= his_str + ",tref=" + strnum(tref); // Reference PV pvref= gravit1 * f0 * kappa1 * tref / pref1; } // Compute sine/cosine of latitude SINLAT= hop("sin", hcoord(t_slab,"y"), name="SINLAT"); COSLAT= hop("cos", hcoord(t_slab,"y"), name="COSLAT"); // Compute planetary vorticity PLVOR= hop(SINLAT, "*", f0, name="PLVOR", units="s-1"); // Get pressure values on model levels (Pa) PMOD= hcoord(t_slab,"z"); if (PMOD.units == "hPa") PMOD= hop(PMOD, "*", 100., units="Pa"); if (PMOD.units != "Pa") error, "Incorrect pressure units"; // Compute delta pressure DPMOD2= hop( hshift(PMOD,"z",count=1), "-", hshift(PMOD,"z",count=-1), name="DPMOD2"); // Compute potential temperature THETA= hop( t_slab, "*", hop( hop(pref1,"/",PMOD,units=""), "^", kappa1) , name="THETA", nohistory=1); PMOD= NULL; // Compute Ertel"s PV using potential temperature as the invariant quantity. // If the result is interpolated to theta ("isentropic") surfaces, then // the field is what is commonly known is isentropic potential vorticity (IPV) // // PV = -g Zeta d(theta)/dp, where Zeta is the absolute vorticity on theta // surfaces. // Zeta is approximated by a standard coordinate transformation from the // input surfaces to theta surfaces. // // PV = -g{(vor+f) d(theta)/dp - [d(theta)/dx dv/dp - d(theta)/dy du/dp]} // // where vor, d(theta)/dx, and d(theta)/dy are on the input surfaces and // are obtained through spectral transforms. // "vor" is the relative vorticity, a code defined derived field. // // The -g factor is included to obtain the result in commonly used units // m2 K s-1 kg-1. // Compute vertical derivative of THETA DTHDP= hop( hop( hshift(THETA,"z",count=1, nohistory=1), "-", hshift(THETA,"z",count=-1, nohistory=1)), "/", DPMOD2,name="DTHDP", nohistory=1); // Compute relative vorticity RELVOR= hshtran(hgather(u_slab,v_slab),vec=1,curl=1,phys=1,preserve=1); // Add planetary vorticity ABSVOR= hop(PLVOR, "+", RELVOR, units="s-1"); // Compute PV term1 PVTERM1= hop(DTHDP, "*", ABSVOR, name="PVTERM1", nohistory=1); DTHDP= NULL; RELVOR= NULL; ABSVOR= NULL; // Compute horizontal gradient of Theta THGRAD= hshtran(THETA,grad=1,phys=1,preserve=1); if (param_set(tref)) { // Compute scale factor for nondimensionalizing PV PVSCALE= hop( pvref, "*", hop( hop(THETA,"/",tref,units=""), "^", 1.+1./kappa1), name="PVSCALE", units="m2 K s-1 kg-1"); } THETA= NULL; // Compute vertical derivative of U, V DUDP= hop( hop( hshift(u_slab,"z",count=1), "-", hshift(u_slab,"z",count=-1)), "/", DPMOD2, name="DUDP"); DVDP= hop( hop( hshift(v_slab,"z",count=1), "-", hshift(v_slab,"z",count=-1)), "/", DPMOD2, name="DVDP"); // Compute PV term2 PVTERM2= hop( hop(THGRAD(I0+0),"*",DVDP), "-", hop(THGRAD(I0+1),"*",DUDP),name="PVTERM2",units="K/Pa s-1"); THGRAD= NULL; DUDP= NULL; DVDP= NULL; DPMOD2= NULL; // Compute PV pv_slab= hop( -gravit1, "*", hop(PVTERM1,"-",PVTERM2), name="PV", units="m2 K s-1 kg-1", nohistory=1); pv_slab.long_name= "Ertel potential vorticity"; PVTERM1= NULL; PVTERM2= NULL; if (param_set(tref)) { // Scale PV pv_slab= hop( pv_slab, "/", PVSCALE, name="ScaledPV", nohistory=1 ); PVSCALE= NULL; } if (!param_set(nohistory)) { // Append history info to slab hset_attr, pv_slab, "data:history", hattr(pv_slab,"data:history") + " hatmepv( " + his_str + ");" } // Return output slab return timer_return(func_name, pv_slab); } func hbin( slab, dim, bin_size, help=, nohistory=) /* DOCUMENT hbin, slab, dim, bin_size, nohistory=0/1 * "Bins" data along dimension DIM (="x"/"y"/"z"/"t"/"i") by averaging the * weighted data values and the corresponding coordinate values. * BIN_SIZE may be a single value, or a list of values, which are used * cyclically. Negative values in a list would correspond to skipped bins. * SEE ALSO: htbin, hsub, hcat */ { func_name= "hbin"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; if (param_set(help)) { write,""; write," Function HBIN bins data along a selected dimension."; write," E.g.,"; write," binned_slab = hbin(slab,'x',2)"; write," bins data in groups of 2 along the X dimension."; write," See also: htbin, hsub, hcat"; write,""; write," Usage: hbin(slab,'x/y/z/t/i',bin_size)"; return timer_return(func_name, NULL); } if (typeof(slab) != "struct_instance") error, "Argument SLAB should be a structure"; if (!is_scalar(slab)) { // Array of hyperslabs; handle recursively slab_array= NULL; for (j=I0; j <= I0+numberof(slab)-1; j++) { tem_slab= hbin( slab(j), dim, bin_size, nohistory=nohistory ); hgrow, slab_array, tem_slab, j, dimsof(slab); } return timer_return(func_name, slab_array); } // Determine dimension to be re-introduced mdim= strloc(HFMT.coordnames,dim,case_fold=1); if (mdim == 0) error, "Invalid dimension specification - "+dim; if (!is_number(bin_size)) error, "Invalid bin size specification"; // Dimension presence codes is_present= slab.dimension(,HFMT.data); if (is_present(mdim-I1) <= 0) error, "Dimension to be binned not present"; ncoord= (hdimsof(slab))(1+mdim-I1); // History string his_str= "<" + slab.name + ">,<" + dim + ">,[" + strcombine(strnum(bin_size),",") + "]"; // Copy slab, ensuring that it contains data slab_copy= hdata(slab); if ( (is_present(ZDIM) > 0) && \ (!is_null(slab_copy.z0)) && (!is_null(slab_copy.zint0)) && \ (hattr(slab_copy,"area_wt:elements") != "dxdydz") ) { // Ensure that Z weights are included slab_copy= hver_wt(slab_copy); } // Count number of bins istart= 1; jsize= 0; nsize= numberof(bin_size); bin_range= NULL; while (istart <= ncoord) { if (bin_size(I0+jsize) > 0) { // Bin to be averaged over iend= istart+bin_size(I0+jsize)-1; if (iend > ncoord) error, "Incomplete bin not permitted"; grow, bin_range, [istart, iend]; } // Move to next bin istart= istart + abs(bin_size(I0+jsize)); jsize= (jsize + 1) % nsize; } nbin= numberof(bin_range)/2; if (nbin == 0) error, "No bins to be averaged over"; reshape_array, bin_range, [2, 2, nbin]; out_slab= NULL; for (j=I0; j <= nbin-I1; j++) { // Extract bin from slab and average irange= bin_range(,j); if (mdim == XDIM+I1) { tem_slab= hsub( slab_copy, limx=irange, x="avg", subscript=1 ); } else if (mdim == YDIM+I1) { tem_slab= hsub( slab_copy, limy=irange, y="avg", subscript=1 ); } else if (mdim == ZDIM+I1) { tem_slab= hsub( slab_copy, limz=irange, z="avg", subscript=1 ); } else if (mdim == TDIM+I1) { tem_slab= hsub( slab_copy, limt=irange, t="avg", subscript=1 ); } else if (mdim == IDIM+I1) { tem_slab= hsub( slab_copy, limi=irange, i="avg", subscript=1 ); } // Re-introduce averaged dimension, with averaged coordinate value tem_slab= hsprout(tem_slab, dim, area_wt=1, z_bot=1, crange="avg"); // Append to slab array hgrow, out_slab, tem_slab, j, nbin, destroy=1; } // Concatenate bin-averages out_slab= hcat(out_slab); if (!param_set(nohistory)) { // Append history info to output slab hset_attr, out_slab, "data:history", hattr(slab_copy,"data:history") + " hbin(" + his_str + ");" } return timer_return(func_name, out_slab); } func hcat(slab1,slab2,help=,ilabel=,iparam=,name=) /* DOCUMENT hcat(slab1,slab2,help=,ilabel=,iparam=,name=) * Concatenates SLAB1 and SLAB2 along a non-conformant dimension, provided * all other dimensions have full strong conformance, and the slabs also * have variable and grid conformance. * (Monotonicity of the concatenated coordinates is required.) * SLAB2 may be omitted, in which case SLAB1 should be an array of hyperslabs, * all of which are concatenated together. * * For the special case where all slabs have four identical dimensions "xyzt", * but do not have the i-dimension, the i-dimension is created, and the slabs * concatenated along that dimension. In this case, an array of new ILABEL * strings may be specified. If not, default labels are generated using * the variables names, if the variable names are different, and/or * the CASE_NAME, if the slabs are not case-conformant, and/or also the * HOR_DOMAIN, VER_DOMAIN attributes, if the slabs are not domain-conformant. * (i.e., attributes values that are not the same for all slabs are appended * to the label string for the slab, and set to null). * IPARAM contains the (optional) parameter values associated with the * i-dimension. * NAME contains the new variable name for the concatenated slab. * (NAME must be specified if the variable names are not the same.) * SEE ALSO: hsub, hsprout, hsplit, hcopy, happend */ { func_name= "hcat"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; if (param_set(help)) { write,""; write," Function HCAT concatenates hyperslabs along a dimension."; write," E.g.,"; write," cat_slab = hcat(slab1,slab2)"; write," concatenates slab1 && slab2 along the non-conformant dimension."; write," cat_slab = hcat(slab_arr)"; write," concatenates the hyperslab array slab_arr along the non-conformant dimension."; write," Tips:"; write," 1. If the hyperslabs have four identical dimensions 'xyzt',"; write," a new i-dimension is automatically created. In this case"; write," ilabel=['label1',label2',...] may be specified to generate"; write," label strings for the i-dimension."; write," 2. iparam=[value1,value2,...] specifies i-parameter values."; write," See also: hsub, hsprout, hsplit, hcopy, happend"; write,""; write," Usage: hcat(slab1,slab2,ilabel=['label1','label2',...],iparam=...)"; return timer_return(func_name, NULL); } // Check conformance of hyperslabs isuperset= (dim_conf= (udim_conf= (unit_conf= NULL))); var_conf= (case_conf= (grid_conf= (domain_conf= NULL))); hconform, slab1, slab2, isuperset, dim_conf, udim_conf, unit_conf, var_conf, case_conf, grid_conf, domain_conf; // Check unit and grid conformance if ( (!unit_conf) || (!grid_conf) ) error, "Unit/grid conformance required for concatenation"; if ((!var_conf) && is_null(name)) error, "NAME parameter must be specified for differing variables"; // Locate non-conformant dimension dwhere= where(dim_conf != 2); if (is_where(dwhere)) { // At least one non-conformant dimension if (numberof(dwhere) > 1) error, "More than one non-conformant dimension; cannot concatenate" // Dimension to be concatenated mcatdim= dwhere(I0) + I1; // Concatenate dimension cat_slab= ncat(mcatdim, slab1, slab2, extend=(case_conf && domain_conf) ); } else { // No non-conformant dimension is_present= slab1(I0).dimension(,HFMT.data); if (is_present(IDIM) > 0) error, "No non-conformant dimension for concatenation"; // Create i-dimension and concatenate is_present(IDIM)= 1; if (is_null(slab2)) { // Array of slabs nslab= numberof(slab1); if (is_null(ilabel)) { // Create label strings ilabel0= array("",nslab); if (!alleq(slab1.name)) ilabel0= ilabel0 + slab1.name; for (islab=I0; islab <= nslab-I1; islab++) { if (!case_conf) ilabel0(islab)= ilabel0(islab)+hattr(slab1(islab), ":case_name"); if (!domain_conf) { hstr= hattr(slab1(islab), ":hor_subdomain"); vstr= hattr(slab1(islab), ":ver_subdomain"); if (hstr != "") ilabel0(islab)= ilabel0(islab) + ":" + hstr; if (vstr != "") ilabel0(islab)= ilabel0(islab) + ":" + vstr; } } } else { // Specified label strings if (numberof(ilabel) != nslab) error, "No. of label strings does not match number of slabs"; ilabel0= ilabel; } // Copy slabs, creating i-dimension slab1c= NULL; iparam1= NULL; for (islab=I0; islab <= nslab-I1; islab++) { tem_slab= NULL; if (!is_null(iparam)) iparam1= iparam(islab); hcopy, slab1, tem_slab, index1=islab, is_present=is_present, ilabel1=[ilabel0(islab)], iparam1=iparam1, ilabel0=ilabel0, iparam0=iparam; nset_attr, "subdomain", tem_slab, IDIM+I1, islab+I1; hgrow, slab1c, tem_slab, islab, [1, nslab], destroy=1; } slab2c= NULL; } else { // Two slab concatenation if (is_null(ilabel)) { // Create label strings ilabel0= array("",2); if (slab1.name != slab2.name) ilabel0= ilabel0 + [slab1.name, slab2.name]; if (!case_conf) { ilabel0= ilabel0 + [ hattr(slab1, ":case_name"), hattr(slab2, ":case_name") ]; } if (!domain_conf) { hstr1= hattr(slab1, ":hor_subdomain"); vstr1= hattr(slab1, ":ver_subdomain"); hstr2= hattr(slab2, ":hor_subdomain"); vstr2= hattr(slab2, ":ver_subdomain"); if (hstr1 != "") ilabel0(I0)= ilabel0(I0) + ":" + hstr1; if (vstr1 != "") ilabel0(I0)= ilabel0(I0) + ":" + vstr1; if (hstr2 != "") ilabel0(I0+1)= ilabel0(I0+1) + ":" + hstr2; if (vstr2 != "") ilabel0(I0+1)= ilabel0(I0+1) + ":" + vstr2; } } else { // Specified label strings if (numberof(ilabel) != 2) error, "No. of label strings does not match number of slabs"; ilabel0= ilabel; } // Copy slabs, creating i-dimension iparam1= NULL; slab1c= NULL; if (!is_null(iparam)) iparam1= iparam(I0); hcopy, slab1, slab1c, is_present=is_present, ilabel1=[ilabel0(I0)], iparam1=iparam1, ilabel0=ilabel0, iparam0=iparam; slab2c= NULL; if (!is_null(iparam)) iparam1= iparam(I0+1); hcopy, slab2, slab2c, is_present=is_present, ilabel1=[ilabel0(I0+1)], iparam1=iparam1, ilabel0=ilabel0, iparam0=iparam; // Initialize i-dimension nset_attr, "subdomain", slab1c, IDIM+I1, 1; nset_attr, "subdomain", slab2c, IDIM+I1, 2; } // Concatenate i-dimension mcatdim= IDIM + I1; cat_slab= ncat(mcatdim, slab1c, slab2c, extend=1 ); } if (!is_null(name)) cat_slab.name= name; if (!case_conf) { // No case conformance hset_attr, cat_slab, ":case_name", ""; hset_attr, cat_slab, ":case_title", ""; } if (!domain_conf) { // No domain conformance hset_attr, cat_slab, ":hor_subdomain", ""; hset_attr, cat_slab, ":ver_subdomain", ""; } return timer_return(func_name, cat_slab); } func hclose( &fstruc, //YORICKoutput: help=,keep=,command=,silent=) /* DOCUMENT hclose,fstruc,help=0/1,keep=0/1,command=,silent=0/1 * Close netCDF history file described structure FSTRUC and associated with * file handle FHANDLE * * Optional parameter: * fstruc -- history file data structure (set to NULL on output) * (if omitted, default history file is closed) * (KEYWORD PARAMETERS) * help -- help option * keep -- if true, do not delete scratch file on closing * command -- execute specified string as operating system command * after closing the file(s), with any % characters * substituted with the file name header (i.e., excluding any * suffix, and leading pathnames). * silent -- silent mode * SEE ALSO: hopen */ { func_name= "hclose"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; if (param_set(help)) { write,""; write," Procedure HCLOSE closes a previously opened netCDF history file."; write," E.g.,"; write," hclose "; write," closes the default history file."; write," hclose,fstruc2 "; write," closes the additional history file referred to by FSTRUC2, FHANDLE2."; write," Tips:"; write," 1. keep=1 prevents scratch files from being deleted"; write," 2. command='filter /tmp/%.his %.nc' applies operating system command,"; write," with occurrences of % substituted by file name header."; write,""; write," See also: hopen"; write,""; write," Usage: hclose[,fstruc][,keep=1][,command='..']"; return timer_return(func_name); } if (is_null(fstruc)) { // Close the default history file if (is_null(DEFAULT_FILE_STRUC)) return timer_return(func_name); fname= DEFAULT_FILE_STRUC.fname; fmeta= DEFAULT_FILE_STRUC.fmeta; fnumber= DEFAULT_FILE_STRUC.fnumber; fscratch= DEFAULT_FILE_STRUC.scratch; fhandle= nget_handle( DEFAULT_FILE_STRUC, erase=1 ); DEFAULT_FILE_STRUC= NULL; } else { // Close additional history file fname= fstruc.fname; fmeta= fstruc.fmeta; fnumber= fstruc.fnumber; fscratch= fstruc.scratch; fhandle= nget_handle( fstruc, erase=1 ); fstruc= NULL; } if (!param_set(silent)) { write, "hclose: closing file "+fname; } // Close netCDF file nc_close, fhandle; // Execute operating system command, if specified if (!is_null(command)) oscommand, command, filename=filename; if ((!param_set(keep)) && fscratch) { // Delete file(s) files= strsplit(fname,","); rm_command= "/bin/rm " + strcombine(files," "); oscommand, rm_command; } return timer_return(func_name); } func hcombine( &fstruc, //YORICKoutput: infiles, help=, varlist=, chunk_size=, create=, new_case=) /* DOCUMENT hcombine, fstruc, infiles, * help=help, varlist=varlist, chunk_size=chunk_size, * create=create, new_case=new_case * * Combine data from input files INFILES and append them out to output file * described by file data structure FSTRUC, * The appended hyperslabs must have strong full conformance in all * dimensions, including reduced ones, and also variable, case, and domain * conformance, with the hyperslabs already present in the output file. * * If VARLIST is specified, it should contain an array of strings, one for * each input file. Each string should contain a comma-separated list of * variables to be read from the file. Null string implies all variables * are to be read from the file. * * If VARLIST is specified as a single string (which may be a null string), * the same variables are read from all the input files, and concatenated * using the HCAT operator. * * If CHUNK_SIZE is specified, process data in chunks of CHUNK_SIZE * records each. This allows a smaller memory footprint for the operation. * * If CREATE="filename" is specified and FSTRUC has a null value, * a new netCDF file is created and its file structure is returned as * FSTRUC. Subsequent calls to HCOMBINE can append to this file. * * NEW_CASE is the (optional) new case name * * SEE ALSO: hsave, happend */ { func_name= "hcombine"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; if (param_set(help)) { write,""; write," Function HCOMBINE combines data from several input files &&"; write," writes/appends them to a single netCDF file. E.g.,"; write," created using HSAVE, && subsequently opened using HOPEN. E.g.,"; write," hcombine, fstruc, ['infile1', 'infile2', ...]"; write," appends all hyperslabs from the specified input files to the file"; write," pointed to by file structure FSTRUC"; write," On returning, the FSTRUC is modified to reflect the appended data."; write," (The hyperslabs SLAB1 && SLAB2, which may themselves be arrays.)"; write," Tips:"; write," 1. varlist=['file1var1,file1var2','file2var1',...] selects variables for input."; write," 2. chunk_size= specified processing in chunks of records."; write," 3. create='filename' creates a new netCDF if FSTRUC is null."; write,""; write," See also: hsave, happend"; write,""; write," Usage: hcombine,fstruc,infiles,varlist=['file1var1,file1var2','file2var1',...],chunk_size=,create='filename'"; return timer_return(func_name); } //IDL2YORICK: error, "Not yet implemented in IDL" ; if (is_null(fstruc) && is_null(create)) error, "Null value for FSTRUC; specify create='filename'" if (is_null(infiles)) error, "No input files specified"; nfiles= numberof(infiles); varlist1= array("",nfiles); if (!is_null(varlist)) { if (numberof(varlist) == 1) { varlist1(*)= varlist; } else { if (numberof(varlist) != nfiles) error, "Incorrect number of strings specified for VARLIST"; varlist1= varlist; } } // Open all input files instruc= NULL; for (j=I0; j <= nfiles-I1; j++) { hopen, infiles(j), fs, fh, time0, date0, filevars, alt=1; grow, instruc, fs; if (j == I0) { // First input file; count records nt= numberof(time0); } else { if (nt != numberof(time0)) error, "No. of records does not match for file "+infiles(j); } // Determine variables to be read if (varlist1(j) == "") varlist1(j)= strcombine(filevars,","); } // Concatenation flag cat_flag= allof(varlist1(*) == varlist1(I0)); // Chunk size nsize= nt; if (param_set(chunk_size)) nsize= chunk_size; // No. of chunks nchunk= 1 + (nt-1)/nsize//; for (ichunk=I0; ichunk <= nchunk-I1; ichunk++) { // Starting, ending time subscripts itime1= 1 + (ichunk-1)*nsize//; if (ichunk != nchunk-I1) { // Not the last chunk itime2= itime1+nsize-1//; } else { // Last chunk itime2= nt//; } write, format="Chunk %4d: [%5d, %5d]\n", ichunk, itime1, itime2; ALLVARS= NULL; if (cat_flag) { // Concatenate variables vars= strsplit(varlist1(1),","); for (ivar=I0; ivar <= numberof(vars)-I1; ivar++) { VARCAT= NULL; for (j=I0; j <= nfiles-I1; j++) { VAR= hget(vars(ivar), fstruc=instruc(j), limt=[itime1,itime2], subscript=1, strip="i"); grow, VARCAT, VAR; } VARCAT= hcat(VARCAT); if (param_set(new_case)) hset_attr, VARCAT, ":case_name", new_case//; grow, ALLVARS, VARCAT; } } else { // Combine variables for (j=I0; j <= nfiles-I1; j++) { // Read input chunks VARS= hget(strsplit(varlist1(j),","), fstruc=instruc(j), limt=[itime1,itime2], subscript=1); if (param_set(new_case)) hset_attr, VARS, ":case_name", new_case//; grow, ALLVARS, VARS; } } // Append variables to file happend, fstruc, ALLVARS, create=create; } // Close input files for (j=I0; j <= nfiles-I1; j++) { hclose, instruc(j); } return timer_return(func_name); } func hconform( slab1, slab2, &isuperset, &dim_conf, &udim_conf, &unit_conf, //YORICKoutput: &var_conf, &case_conf, &grid_conf, &domain_conf, //YORICKoutput: help=, reduced=, verbose=) /* DOCUMENT hconform, slab1, slab2, * isuperset, dim_conf, udim_conf, unit_conf, * var_conf, case_conf, grid_conf, domain_conf, * help=, reduced=, verbose= * Checks conformance characteristics of two or more hyperslabs * * Input parameters: * slab1, slab2 -- hyperslab data structures * (slab1 may also be an array of hyperslabs, * in which case slab2 should be set to NULL) * Output parameters: * isuperset -- = 0 => no slab has superset of all dimensions present * > 0 => all dimensions are at least weakly conformant, * and slab ISUPERSET contains superset of all * dimensions present * dim_conf(SDIM) -- dimension conformance value (for each standard dimension) * 0 => corresponding dimensions are non-conformable * >0 => corresponding dimensions are strongly conformable * <0 => corresponding dimensions are weakly conformable * [ 2 => present in all and identical, or missing in all slabs * ("strong full conformance") * 1 => identical where present, missing in one or more slabs * ("strong broadcast conformance") * -1 => same length, where present * ("weak broadcast conformance") * -2 => same length in all ] * ("weak full conformance") * udim_conf(SDIM)-- dimension units and other attributes conformance value * (for each standard dimension) * 1 => dimension units/long name/grid/reduction-ops are same * 0 => dimension units/long name/... are different * unit_conf -- true if all slabs have the same data units * var_conf -- true if all slabs have unit conformance, and the same * variable names, long names, and time representation * case_conf -- true if all slabs have the same case name attribute * grid_conf -- true if all slabs have the same full domain grid * domain_conf -- true if all slabs have grid conformance, and * have the same horizontal and vertical subdomain name * attributes * (KEYWORD PARAMETERS) * help -- help option * reduced -- if true, check conformance of reduced dimensions as well * verbose -- verbose option (prints out additional info) * SEE ALSO: hop, hcat, happend */ { func_name= "hconform"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; if (param_set(help)) { write,""; write," Procedure HCONFORM checks conformance of two || more hyperslabs."; write," E.g.,"; write," hconform, slab_a, slab_b, isuperset, dim_conf, udim_conf,"; write," unit_conf, var_conf, case_conf, grid_conf, domain_conf"; write," check conformance of slab_a && slab_a, returning the various"; write," conformance flags."; write," hconform, slab_arr, NULL, ..."; write," checks conformance of an array of hyperslabs slab_arr."; write," Tips:"; write," 1. verbose=1 option prints out conformance details."; write," 2. See source code && hyperslab structure definition for the"; write," meaning of all the output parameters."; write," See also: hop, hcat, happend"; write,""; write," Usage: hconform, slab_a, slab_b, isuperset, dim_conf, udim_conf,"; write," unit_conf, var_conf, case_conf, grid_conf, domain_conf,"; write," verbose=1"; return timer_return(func_name); } if ( (typeof(slab1) != "struct_instance") && \ (typeof(slab2) != "struct_instance") ) error, "Argument SLAB1, SLAB2 should be structures"; if (is_null(slab2)) { // Only one argument present; check no. of elements nslab= numberof(slab1); // Ensure that there are at least two slabs if (nslab <= 1) error, "At least two slab arguments required for conformance testing"; slab_array= 1; } else { // Two arguments present if ((!is_scalar(slab1)) || (!is_scalar(slab2))) error, "SLAB1 && SLAB2 should both be scalars"; nslab= 2; slab_array= 0; } // Conformance strings conf_str= ["WF", "WB", "NC", "SB", "SF"]; flag_str= ["F", "T"]; // Initialize dimension conformance array (set to strong full conformance) dim_conf= array(2,SDIM); // Initialize dimension unit conformance array (set to true) udim_conf= array(1,SDIM); // Initialize all dimensions present array (set to true) all_present= array(1,nslab); for (mdim=I0; mdim <= SDIM-I1; mdim++) { dim_presence= array(1,nslab); coord= NULL; if (!slab_array) { // SLAB1 is not an array if ( (slab1.dimension(mdim,HFMT.data) > 0) || \ ( (slab1.dimension(mdim,HFMT.data) < 0) && \ param_set(reduced) ) ) { // Dimension is/was present in slab coord= ngetcoord(slab1,mdim+I1); credop= slab1.reduced(mdim); cgrid= nattr("grid",slab1,mdim+I1); cunits= nattr("units",slab1,mdim+I1); clong_name= nattr("long_name",slab1,mdim+I1); } else { coord= NULL; // Dimension not present in slab; reset slab dimension presence flag dim_presence(1-I1)= 0; } if ( (slab2.dimension(mdim,HFMT.data) > 0) || \ ( (slab2.dimension(mdim,HFMT.data) < 0) && \ param_set(reduced) ) ) { // Dimension is/was present in slab coord2= ngetcoord(slab2,mdim+I1); if (!is_null(coord)) { // Check conformance of dimension grid, units, and long name udim_conf(mdim)= udim_conf(mdim) && \ (nattr("grid",slab2,mdim+I1) == cgrid) && \ (nattr("units",slab2,mdim+I1) == cunits) && \ (nattr("long_name",slab2,mdim+I1) == clong_name); if (param_set(reduced)) { // Check reduction operation on dimension udim_conf(mdim)= udim_conf(mdim) && \ (slab2.reduced(mdim) == credop); } if (numberof(coord2) == numberof(coord)) { // Same dimension length; check coordinate values if ( (!array_eq(coord2,coord,epsilon=HFMT.epscoord)) && \ (dim_conf(mdim) > 0) ) { // Weak full conformance dim_conf(mdim)= -2; } } else { // Different dimension lengths; non-conformance dim_conf(mdim)= 0; } } } else { // Dimension not present in slab; reset slab dimension presence flag dim_presence(2-I1)= 0; } } else { // SLAB1 is an array for (islab=I0; islab <= nslab-I1; islab++) { if ( (slab1(islab).dimension(mdim,HFMT.data) > 0) || \ ( (slab1(islab).dimension(mdim,HFMT.data) < 0) && \ param_set(reduced) ) ) { // Dimension is/was present in slab if (is_null(coord)) { coord= ngetcoord(slab1(islab),mdim+I1); credop= slab1(islab).reduced(mdim); cgrid= nattr("grid",slab1(islab),mdim+I1); cunits= nattr("units",slab1(islab),mdim+I1); clong_name= nattr("long_name",slab1(islab),mdim+I1); } else { coord2= ngetcoord(slab1(islab),mdim+I1); // Check conformance of dimension grid, units, and long name udim_conf(mdim)= udim_conf(mdim) && \ (nattr("grid",slab1(islab),mdim+I1) == cgrid) && \ (nattr("units",slab1(islab),mdim+I1) == cunits) && \ (nattr("long_name",slab1(islab),mdim+I1) == clong_name); if (param_set(reduced)) { // Check reduction operation on dimension udim_conf(mdim)= udim_conf(mdim) && \ (slab1(islab).reduced(mdim) == credop); } if (numberof(coord2) == numberof(coord)) { // Same dimension lengths; check coordinate values if ( (!array_eq(coord2,coord,epsilon=HFMT.epscoord)) && \ (dim_conf(mdim) > 0) ) { // Weak full conformance dim_conf(mdim)= -2; } } else { // Differing dimension lengths; non-conformance dim_conf(mdim)= 0; } } } else { // Dimension not present in slab; reset slab dimension presence flag dim_presence(islab)= 0; } } } if (!(allof(dim_presence) || noneof(dim_presence)) ) { // Some slabs with missing dimension; downgrade to broadcast conformance dim_conf(mdim)= dim_conf(mdim) / 2; } if (!udim_conf(mdim)) { // Dimension attributes not identical; downgrade to weak conformance dim_conf(mdim)= -abs(dim_conf(mdim)); } if (anyof(dim_presence(*))) { // Apply AND operation on combined slabs dimension presence flag all_present(*)= all_present(*) * dim_presence(*); } } // Superset slab number isuperset= 0; if ( (noneof(abs(dim_conf) == 0)) && anyof(all_present) ) { // All dimensions are at least weakly conformant, and // at least one slab has superset of all dimensions present; locate it isuperset= (where(all_present))(I0) + I1; } if (!slab_array) { // Units conformance unit_conf= (slab1.units == slab2.units); // Variable conformance var_conf= unit_conf && (slab1.name == slab2.name ) && \ (slab1.long_name == slab2.long_name) && \ (hattr(slab1,"data:time_rep") == hattr(slab2,"data:time_rep")); // Case conformance case_conf= (hattr( slab1,":case_name") == hattr( slab2,":case_name")); // Grid conformance grid_conf= (slab1.structure == slab2.structure) && \ (hattr(slab1,"x:period") == hattr(slab2,"x:period")) && \ (hattr(slab1,"x:rotated") == hattr(slab2,"x:rotated")) && \ array_eq(deref(slab1.x0), deref(slab2.x0), epsilon=HFMT.epscoord) && \ array_eq(deref(slab1.y0), deref(slab2.y0), epsilon=HFMT.epscoord) && \ array_eq(deref(slab1.z0), deref(slab2.z0), epsilon=HFMT.epscoord) && \ array_eq(deref(slab1.ilabel0), deref(slab2.ilabel0), epsilon=HFMT.epscoord) && \ array_eq(deref(slab1.iparam0), deref(slab2.iparam0), epsilon=HFMT.epscoord) && \ array_eq(deref(slab1.xint0), deref(slab2.xint0), epsilon=HFMT.epscoord) && \ array_eq(deref(slab1.yint0), deref(slab2.yint0), epsilon=HFMT.epscoord) && \ array_eq(deref(slab1.zint0), deref(slab2.zint0), epsilon=HFMT.epscoord); // Domain conformance domain_conf= grid_conf && \ (hattr(slab1,":hor_subdomain") == hattr(slab2,":hor_subdomain")) && \ (hattr(slab1,":ver_subdomain") == hattr(slab2,":ver_subdomain")); } else { // Units conformance unit_conf= alleq( slab1.units ); // Variable conformance var_conf= unit_conf && alleq(slab1.name ) && \ alleq(slab1.long_name) && \ alleq( hattr(slab1,"data:time_rep") ); // Case conformance case_conf= alleq( hattr( slab1,":case_name") ); // Grid conformance grid_conf= 0; if ( alleq(slab1.structure) ) { grid_conf= 1; x0= deref(slab1(I0).x0); y0= deref(slab1(I0).y0); z0= deref(slab1(I0).z0); ilabel0= deref(slab1(I0).ilabel0); iparam0= deref(slab1(I0).iparam0); xint0= deref(slab1(I0).xint0); yint0= deref(slab1(I0).yint0); zint0= deref(slab1(I0).zint0); for (islab=I0+1; islab <= nslab-I1; islab++) { grid_conf= grid_conf && \ alleq(hattr(slab1,"x:period")) && \ alleq(hattr(slab1,"x:rotated")) && \ array_eq( deref(slab1(islab).x0), x0, epsilon=HFMT.epscoord ) && \ array_eq( deref(slab1(islab).y0), y0, epsilon=HFMT.epscoord ) && \ array_eq( deref(slab1(islab).z0), z0, epsilon=HFMT.epscoord ) && \ array_eq( deref(slab1(islab).ilabel0), ilabel0, epsilon=HFMT.epscoord ) && \ array_eq( deref(slab1(islab).iparam0), iparam0, epsilon=HFMT.epscoord ) && \ array_eq( deref(slab1(islab).xint0), xint0, epsilon=HFMT.epscoord ) && \ array_eq( deref(slab1(islab).yint0), yint0, epsilon=HFMT.epscoord ) && \ array_eq( deref(slab1(islab).zint0), zint0, epsilon=HFMT.epscoord ); } } // Domain conformance domain_conf= grid_conf && alleq(hattr(slab1,":hor_subdomain")) && \ alleq(hattr(slab1,":ver_subdomain")); } if (param_set(verbose)) { // Print out conformance information write, "Superset-slab-no: ", strnum(isuperset); write, "Dimension-conformance: ", strcombine( conf_str(dim_conf(*)+2+I0), "," ); write, "Dimension-unit-conformance: ", strcombine( flag_str(udim_conf(*)+I0), "," ); write, "Unit-conformance: ", flag_str(unit_conf+I0); write, "Variable-conformance: ", flag_str(var_conf+I0); write, "Case-conformance: ", flag_str(case_conf+I0); write, "Grid-conformance: ", flag_str(grid_conf+I0); write, "Domain-conformance: ", flag_str(domain_conf+I0); } return timer_return(func_name); } func hcont(slab,help=,mval=,fill=,width=,scolor=, lcolor=,ltype=,proj=,ppars=,rotx=, terrain=,transp=) /* DOCUMENT hcont,slab,help=,mval=,fill=,width=,scolor=, * lcolor=,proj=,ppars=,rotx=,terrain=,transp= * Superimposes continental outlines on previous plot for * atmospheric/oceanic hyperslab data, and optionally solid-fills land area * * slab -- hyperslab containing atmospheric/oceanic data * (KEYWORD PARAMETERS) * help -- help option * mval -- for atmospheric data, ORO mask value to be contoured * for oceanic data, depth index values < MVAL are contoured * (default value is 1) * fill -- if set, solid-fill land area * width -- line thickness (default 0.5) * scolor -- color for solid fills (default: light gray) * lcolor -- color for continental outlines (default: black) * ltype -- linetype for continental outlines (default: solid) * proj -- projection ("","NHPOLAR","SHPOLAR","MOLLWEIDE",...) * (NOTE: projection names may be abbreviated) * ppars -- projection parameters ([start_lon, extreme_lat]/...) * rotx -- if defined, rotate X-coordinate by angle ROTX * terrain -- slab containing terrain elevation (in m) above * sea level to be used to draw continental outlines * (special case: if terrain==1, use 30-minute topographic * elevation values from TerrainBase dataset) * transp -- if set, interchanges X/Y axes * SEE ALSO: hplot */ { func_name= "hcont"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; if (param_set(help)) { write,""; write," Function HCONT superimposes continental outlines on previous plot"; write," for atmospheric/oceanic hyperslab data."; write," E.g.,"; write," hcont, slab"; write," superimposes continental outlines over previous plot."; write," Tips:"; write," 1. fill=1 produces solid-fills for continents."; write," 2. width=line_width, scolor=fill_color, lcolor=line_color, "; write," ltype=line_type control line/fill details."; write," 3. proj='projection', ppars=[] "; write," controls the projection."; write," 4. rotx=x_rotation_angle rotates the X domain."; write," 5. terrain=1 uses 30-minute topographic elevation database"; write," 6. mval=mask_value determines value to be contoured."; write," See also: hplot"; write,""; write," Usage: hcont,slab,mval=mask_value,width=line_width,scolor=fill_color,"; write," lcolor=line_color,proj=...,ppars=[...],rotx=rot_angle,terrain=1,fill=1"; return timer_return(func_name); } if (is_null(mval)) mval= 1; if (is_null(width)) width= 0.5; if (is_null(scolor)) scolor= 10; if (is_null(lcolor)) lcolor= 2; if (is_null(ltype)) ltype= 0; if (typeof(slab) != "struct_instance") error, "Argument SLAB should be a structure"; if (!is_scalar(slab)) error, "Argument SLAB should be a scalar"; is_present= slab.dimension(,HFMT.data); // Determine X/Y dimensions if ( (is_present(XDIM) <= 0) || (is_present(YDIM) <= 0) ) error, "X/Y dimensions not present in slab"; ddims= hdimsof(slab); nx= ddims(1+XDIM); ny= ddims(1+YDIM); // Return if too few grid points if ( (nx < 2) || (ny < 2) ) return timer_return(func_name); if (param_set(terrain)) { if (typeof(terrain) == "struct_instance") { elev_slab= terrain; } else { if (strlen(HFMT.hopsroot) == 0) error, "Set environment variable HOPS_ROOT to root directory of HOPS installation"; // Read topographic elevation data elev_slab= hget("TOPO30", file=HFMT.hopsroot+"/data/topo30.nc"); } // Reduce to subdomain like slab elev_slab= hsub(elev_slab, like=slab); //HARD-EXTENSIONS-BEGIN: // Determine continental mask if (slab.structure == "HYPERSLAB1.0_SPH_SIG_ATM") { // SPH_SIG_ATM extension elev_slab= hop(elev_slab, ">", 0.); } else if (slab.structure == "HYPERSLAB1.0_SPH_SIG_OCN") { // SPH_SIG_OCN extension elev_slab= hop(elev_slab, ">", 0.); } else { error, "Unrecognized slab type"; } //HARD-EXTENSIONS-END: if (param_set(fill)) { // Solid fill all land areas hplot, elev_slab, proj=proj, ppars=ppars, rotx=rotx, transp=transp, overplot=1, levs=[0.5, 1.5], nomask=1, mix=[scolor,scolor], fill=1; } else { // Plot continental outline hplot, elev_slab, proj=proj, ppars=ppars, rotx=rotx, transp=transp, overplot=1, levs=[0.5, 1.5], nomask=1, width=width, line_color=lcolor, type=ltype, c_labels=0; } } else { // Current horizontal subdomain offsets ix= nattr("subdomain",slab,XDIM+I1); iy= nattr("subdomain",slab,YDIM+I1); if ( (ix < 0) || (iy < 0)) error, "X && Y subdomains should be contiguous to plot continental outlines"; if (ix == 0) ix= 1; if (iy == 0) iy= 1; // Get regular grid coordinates for subdomain xc= ngetcoord(slab, XDIM+I1, grid=1); yc= ngetcoord(slab, YDIM+I1, grid=1); if (numberof(xc) < (nx+ix-1)) nx= numberof(xc) - (ix-1); if (numberof(yc) < (ny+iy-1)) ny= numberof(yc) - (iy-1); xc= xc(ix-I1:ix+nx-1-I1); yc= yc(iy-I1:iy+ny-1-I1); // Current rotation state x_rotated= hattr(slab, "x:rotated"); //HARD-EXTENSIONS-BEGIN: // Slab type, full domain land mask if (slab.structure == "HYPERSLAB1.0_SPH_SIG_ATM") { // SPH_SIG_ATM extension lmask= (*(slab.hgrid0) == mval); } else if (slab.structure == "HYPERSLAB1.0_SPH_SIG_OCN") { // SPH_SIG_OCN extension lmask= (*(slab.kmax0) < mval); } else { error, "Unrecognized slab type"; } //HARD-EXTENSIONS-END: if (x_rotated != 0) { // Apply rotation on full domain masks lmask= rangeop(lmask,"rot",1,count=-x_rotated); } // Reduce land mask to subdomain lmask= lmask(ix-I1:ix+nx-1-I1, iy-I1:iy+ny-1-I1); if (param_set(rotx) && (nattr("subdomain",slab,XDIM+I1) == 0)) { // Rotate X-coordinate xrotate,360.0,rotx,xc,irot; lmask= rangeop(lmask,"rot",1,count=irot); } // Compute grid-box boundaries xb= xc(I0) - 0.5*(xc(I0+1)-xc(I0)); grow, xb, 0.5*(xc(I0:nx-I1-1) + xc(I0+1:nx-I1)); grow, xb, xc(nx-I1) + 0.5*(xc(nx-I1) - xc(nx-I1-1)); yb= yc(I0) - 0.5*(yc(I0+1)-yc(I0)); grow, yb, 0.5*(yc(I0:ny-I1-1) + yc(I0+1:ny-I1)); grow, yb, yc(ny-I1) + 0.5*(yc(ny-I1) - yc(ny-I1-1)); nxb= numberof(xb); nyb= numberof(yb); // Generate grid-box boundary mesh xm= xb(,-:1:nyb)// //IDL2YORICK: xm= xb # replicate(1, nyb) ym= yb(-:1:nxb,)// //IDL2YORICK: ym= replicate(1, nxb) # yb if (!is_null(proj)) { // Projection specified pmesh= array(double, 2, numberof(xm) ); pmesh(I0,)= xm(*); pmesh(I0+1,)= ym(*); pmesh= fproject( pmesh, proj=proj, ppars=ppars); xm(*)= pmesh(I0,); ym(*)= pmesh(I0+1,); } if (param_set(transp)) { // Transpose X/Y coordinates lmask= transpose(lmask); temp= xm; xm= ym; ym= temp; temp= nx; nx= ny; ny= temp; } if (param_set(fill)) { // Solid fill all land areas for (j=I0; j <= ny-I1; j++) { for (i=I0; i <= nx-I1; i++) { if (lmask(i,j) != 0) { plfp,[char(scolor)],[ym(i,j), ym(i,j+1), ym(i+1,j+1), ym(i+1,j)], [xm(i,j), xm(i,j+1), xm(i+1,j+1), xm(i+1,j)], [4]; } } } } else { // Draw continental outlines for (j=I0; j <= ny-I1; j++) { for (i=I0; i <= nx-I1; i++) { if (lmask(i,j) != 0) { if (i < nx-1) { if (lmask(i+1,j) == 0) plg, [ym(i+1,j), ym(i+1,j+1)], [xm(i+1,j), xm(i+1,j+1)], width=width, type=ltype, color=char(lcolor), marks=0; } if (i > 0) { if (lmask(i-1,j) == 0) plg, [ym(i,j), ym(i,j+1)], [xm(i,j), xm(i,j+1)], width=width, type=ltype, color=char(lcolor), marks=0; } if (j < ny-1) { if (lmask(i,j+1) == 0) plg, [ym(i,j+1), ym(i+1,j+1)], [xm(i,j+1), xm(i+1,j+1)], width=width, type=ltype, color=char(lcolor), marks=0; } if (j > 0) { if (lmask(i,j-1) == 0) plg, [ym(i,j), ym(i+1,j)], [xm(i,j), xm(i+1,j)], width=width, type=ltype, color=char(lcolor), marks=0; } } } } } } return timer_return(func_name); } func hcoord(slab,dim,help=,all=,data_prec=) /* DOCUMENT hcoord(slab,dim,help=,all=0/1,data_prec=) * returns a version of hyperslab SLAB with the variable values replaced * by the coordinate values corresponding to dimension DIM * ("x"/"y"/"z"/"t"/"i"), and the variable name/units changed to that of * the coordinate. * (The coordinate values will be defined everywhere; i.e., there will be * no missing values.) * If DIM=="area_wt", the area weights are returned as a hyperslab. * If DIM=="z_bot", the bottom Z values are returned as a hyperslab. * If ALL==1, coordinate values are broadcast to have all data dimensions. * DATA_PREC, if set, specifies the data precision ("float"/"double") * for the coordinate values. (The default is to have the same precision * as the data.) * SEE ALSO: hinterp, hver_wt */ { func_name= "hcoord"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; if (param_set(help)) { write,""; write," Function HCOORD returns slab with data values substituted by the"; write," coordinate values for the specified dimension. E.g.,"; write," z_slab = hcoord(slab, 'z')"; write," returns a slab containing Z coordinate values."; write," Tips:"; write," 1. all=1 ensures that coordinates are broadcast to data dimensions."; write," See also: hinterp, hver_wt"; write,""; write," Usage: hcoord(slab, 'x/y/z/t/i', data_prec='float/double')"; return timer_return(func_name, NULL); } if (typeof(slab) != "struct_instance") error, "Argument SLAB should be a structure"; if (!is_scalar(slab)) { // Array of hyperslabs; handle recursively slab_array= NULL; for (i=I0; i <= I0+numberof(slab)-1; i++) { tem_slab= NULL; tem_slab= hcoord( slab(i), dim ); hgrow, slab_array, tem_slab, i, dimsof(slab), destroy=1; } return timer_return(func_name, slab_array); } if (slab.type(HFMT.data) == "") error, "Error - null data values"; // Determine coordinate precision if (param_set(data_prec)) { coord_prec= data_prec; } else { // Data type info data_type= slab.type(HFMT.data); if (data_type == "struct_instance") { loc= *(slab.data); data_type= loc.type(HFMT.data); } coord_prec= ndataprec(data_type); } if (strtolower(dim) == "area_wt") { // Return area weights slab if (is_null(slab.area_wt)) error, "Area weights not available in slab"; area_wt1= typeconv( coord_prec, *(slab.area_wt) ); apresent= slab.dimension(,HFMT.area_wt); // Create area weights slab new_slab= NULL; hcopy, slab, new_slab, is_present=apresent, data=area_wt1, missing_value="", area_wt1="", area_wt_dims=array(long,SDIM), z_bot1="", z_bot_dims=array(long,SDIM); // Change variable name/units new_slab.name= slab.name + "_area_wt"; new_slab.units= hattr(slab,"area_wt:units"); new_slab.long_name= "Area weights"; return timer_return(func_name, new_slab); } if (strtolower(dim) == "z_bot") { // Return Z_BOT values slab if (is_null(slab.z_bot)) error, "Bottom Z values not available in slab"; z_bot1= typeconv( coord_prec, *(slab.z_bot) ); zpresent= slab.dimension(,HFMT.z_bot); // Create Z_BOT values slab new_slab= NULL; hcopy, slab, new_slab, is_present=zpresent, data=z_bot1, area_wt1="", area_wt_dims=array(long,SDIM), z_bot1="", z_bot_dims=array(long,SDIM); // Change variable name/units new_slab.name= slab.name + "_z_bot"; new_slab.units= hattr(slab,"z_bot:units"); new_slab.long_name= "Bottom Z values"; return timer_return(func_name, new_slab); } // Determine selected dimension mdim= strloc(HFMT.coordnames,dim,case_fold=1); if (mdim == 0) error, "Invalid dimension specification - "+dim; // Dimension presence/reduction codes is_present= slab.dimension(,HFMT.data); is_reduced= slab.reduced(*); if (is_present(mdim-I1) <= 0) error, "Dimension needs to be present in the data to generate coordinates"; // Reduce all dimensions new_reduced= is_reduced; new_reduced(where(is_present > 0))= -strloc(HFMT.reduceops, "avg"); new_present= -abs(is_present); // Sigma coordinate flag sigma_coord= (nattr("units",slab,ZDIM+I1) == "hybrid_sigma_pressure") || \ (nattr("units",slab,ZDIM+I1) == "sigma_level"); if ((mdim != ZDIM+I1) || (!sigma_coord)) { // Not sigma Z coordinate // Coordinate variable name and units if (mdim == IDIM+I1) { new_name= "iparam"; new_long_name= "parameter"; } else { new_name= HFMT.dimnames(mdim-I1); new_long_name= nattr("long_name", slab, mdim); } new_units= nattr("units", slab, mdim); // Get coordinate values coord= ngetcoord(slab, mdim, iparam=1); ncoord= numberof(coord); if (ncoord == 0) error, "Coordinate values not available for dimension"; // Reintroduce selected dimension new_present(mdim-I1)= is_present(mdim-I1); // Copy coordinate values as data values data1= typeconv( coord_prec, coord ); // Reshape data array to have right dimensionality ddims= array(1,SDIM+1); ddims(I0)= SDIM; ddims(1+mdim-I1)= ncoord; reshape_array, data1, ddims; } else { // Sigma Z coordinate if (slab.type(HFMT.z_bot) == "") error, "Bottom Z values not available to compute vertical coordinate"; // Copy Z_BOT values z_bot1= *(slab.z_bot); z_bot_ref1= typeconv( coord_prec, hattr(slab, "z_bot:ref") ); // Coordinate variable name and units zb_long_name= hattr(slab, "z_bot:long_name"); new_name= "z"; new_units= hattr(slab, "z_bot:units"); new_long_name= ""; if (zb_long_name == "surface_pressure") { // Pressure coordinate new_name= "p"; new_long_name= "pressure"; } else if (zb_long_name == "bottom_depth") { // Depth coordinate new_long_name= "depth"; } // Copy Z coordinate z= *(slab.z); nz= numberof(z); // Locate Z offset iz= nattr("subdomain",slab,ZDIM+I1); if (iz < 0) error, "Cannot compute vertical weights for non-contiguous Z subdomain"; if (iz == 0) iz= 1; // Copy missing value missing_value= deref(slab.missing_value); // Get full domain sigma values if (is_present(ZDIM) == 1) { // Data on regular Z grid if (is_null(slab.sigma0)) error, "Full domain regular sigma values not available"; sigmafull= *(slab.sigma0); } else { // Data on interfacial Z grid if (is_null(slab.sigmaint0)) error, "Full domain interfacial sigma values not available"; sigmafull= *(slab.sigmaint0); } // Set precision of sigma values sigmafull= typeconv( coord_prec, sigmafull ); // List of current data dimensions ddims= hdimsof(slab); // Z_BOT dimensions new_ddims= hdimsof(slab, z_bot=1); // Re-introduce all dimensions present in Z_BOT zwhere= where(slab.dimension(,HFMT.z_bot) > 0); if (is_where(zwhere)) { new_present(zwhere)= is_present(zwhere); new_reduced(zwhere)= 0; } // Add Z dimension new_ddims(1+ZDIM)= ddims(1+ZDIM); new_present(ZDIM)= is_present(ZDIM); timer_call,"hcoord-crit"; //CRITICAL-SECTION-BEGIN: if (!is_null(missing_value)) { // Set missing bottom Z values to reference bottom value (temporarily) // (This is just for computational convenience to avoid overflows.) z_missing= where(z_bot1 == missing_value); if (is_where(z_missing)) z_bot1(z_missing)= z_bot_ref1; } // Set precision of Z_BOT values if (typeof(z_bot1) != coord_prec) z_bot1= typeconv( coord_prec, z_bot1 ); // Focus on Z dimension new_zfocus= dim_reshape( new_ddims, focus=ZDIM+I1 ); nxy= new_zfocus(I0+1); nti= new_zfocus(I0+3); // Reshape/create arrays reshape_array, z_bot1, [2, nxy, nti]; data1= array( typeconv(coord_prec,0), [3, nxy, nz, nti] ); // Compute vertical coordinate for (k=I0; k <= I0+nz-1; k++) { data1(,k,)= sigmafull(I0, k+iz-1)*z_bot_ref1 + sigmafull(I0+1, k+iz-1)*z_bot1(,); } // Reshape data array to final dimensions reshape_array, data1, new_ddims; timer_return,"hcoord-crit"; //CRITICAL-SECTION-END: } if (param_set(all)) { // Broadcast coordinate values to all dimensions present in original data dwhere= where(slab.dimension(,HFMT.data) > 0); new_present(dwhere)= is_present(dwhere); new_reduced(dwhere)= 0; data1= broadcast(data1, hdimsof(slab)); // Copy slab, replacing data values with coordinate values (no missing values) new_slab= NULL; hcopy, slab, new_slab, is_present=new_present, is_reduced=new_reduced, missing_value="", data=data1; } else { // Create slab containing coordinate values (no missing values, area weights) new_slab= NULL; hcopy, slab, new_slab, is_present=new_present, is_reduced=new_reduced, missing_value="", data=data1, area_wt1="", area_wt_dims=array(long,SDIM), z_bot1="", z_bot_dims=array(long,SDIM); } // Change variable name/units new_slab.name= new_name; new_slab.units= new_units; new_slab.long_name= new_long_name; return timer_return(func_name, new_slab); } func hcopy( slab1, &slab2, //YORICKoutput: help=, index1=, index2=, overwrite=, x1=, y1=, z1=, time1=, date1=, ilabel1=, iparam1=, data=, area_wt1=, z_bot1=, missing_value=, x0=, y0=, z0=, xint0=, yint0=, zint0=, ilabel0=, iparam0=, structure0=, reset_vars=, extra_atts=, area_wt_dims=, z_bot_dims=, is_present=, is_reduced=) /* DOCUMENT hcopy, slab1, slab2, help=help, index1=index1, index2=index2, * x1=x1, y1=y1, z1=z1, time1=time1, date1=date1, ilabel1=ilabel1, * iparam1=iparam1, data=data, area_wt1=area_wt1, z_bot1=z_bot1, * missing_value=missing_value, * x0=x0, y0=y0, z0=z0, xint0=xint0, yint0=yint0, zint0=zint0, * ilabel0=ilabel0, iparam0=iparam0, structure0=structure0, * reset_vars=reset_vars, extra_atts=extra_atts, * area_wt_dims=area_wt_dims, z_bot_dims=z_bot_dims, * is_present=is_present, is_reduced=is_reduced * Copy SLAB1 to SLAB2, with specified keyword parameter fields modified. * If SLAB2 is an array of hyperslabs, copy SLAB1 into element * SLAB2(INDEX2) of the array. * If SLAB2 is null, return in SLAB2 a new modified copy of SLAB1, * with default field values determined by the STRUCTURE keyword parameter. * If OVERWRITE == 1, overwrite values in SLAB2. * Either SLAB1 or SLAB2 may be arrays. If they are both arrays, they * should have the same number of elements. Otherwise, if SLAB1 is an array, * SLAB2 should be null. If SLAB2 is an array, * INDEX1/INDEX2 should be specified to choose the element of the array to * copy from or to copy to. * If SLAB2 is undefined, a new hyperslab is created. * * RESET_VARS is a list of variables whose attribute lists need to be reset * (i.e., standard attributes are set to "zero" values and the non-standard * attributes are "deleted"). * (*NOTE* Attributes of coordinate variables for undefined dimensions are * automatically reset.) * EXTRA_ATTS is a 3xn array of string triplets of the form * [ ["variable_name", "attribute", "data_type"], ...] * of extra attributes to be added to the hyperslab. * IS_PRESENT specifies what dimensions are present, and what their grid types * are. * IS_REDUCED specifies what reduction operations have been carried out on the * dimensions. * AREA_WT_DIMS specifies the current dimensions present for the area weights. * Z_BOT_DIMS specifies the current dimensions present for the Z_BOT array. * (See the hyperslab format definition document for more details.) * If keyword parameter IS_PRESENT is specified, the data array, area weights * array, and the z_bot array are reshaped to conform to the standard * coordinate dimensions. * **Important note** * Specifying NULL values to parameters is equivalent to omitting them, * which would result in their values being copied from the source SLAB1. * To prevent that and to actually assign a NULL value for a parameter, * set numeric parameters to a string value (e.g., x1="", missing_value="", * data=""), and string parameters to a numeric value (e.g., ilabel1=0). * * Output parameter: * slab2 -- new slab * SEE ALSO: hget, hsub, hcat, hsave, hattr, hdata, ncopyatt */ { func_name= "hcopy"; timer_call,func_name; // Hyperslab include file // hcom.pro: // Hyperslab include file // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Hyperslab COMMON block //IDLbegin: //:common hcom, HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST //IDL-only: //IDLend: //YORICKbegin: extern HFMT, DEFAULT_FILE_STRUC, FILE_HANDLE_LIST; //YORICKend: // Initialize hyperslab format descriptor structure, if necessary if (is_null(HFMT)) ninit; // No. of standard dimensions SDIM= 5; // Indices for the five standard dimensions XDIM=I0+0 ; YDIM=I0+1 ; ZDIM=I0+2 ; TDIM=I0+3 ; IDIM=I0+4; if (param_set(help)) { write,""; write," Function HCOPY creates/copies/modifies hyperslabs."; write," E.g., "; write," hcopy, slab1, slab2 "; write," copies slab1 to slab2, with no modifications."; write," In the above, slab1 && slab2 may be arrays with the same dimensions."; write," hcopy, slab1, slab2, index2=2"; write," copies 'scalar' slab1 to array element slab2(2)."; write," hcopy, slab1, slab2, data=new_data"; write," copies slab1 to slab2, with modified data array new_data."; write," hcopy, NULL, new_slab, structure0='HYPERSLAB1.0_SPH_SIG_OCN'"; write," creates a new ocean data hyperslab new_slab."; write," Tips:"; write," 1. structure0='HYPERSLAB1.0_...' may be used to define the hyperslab structure attribute (for new slabs)."; write," 2. overwrite=1 allows destination slab to be overwritten."; write," 3. x1=..., y1=..., z1=..., time1=..., date1=..., ilabel1=..., iparam1=..."; write," may be used to modify the coordinate arrays."; write," 4. is_present=... may be used to add/delete dimensions."; write," 5. is_reduced=... may be used to specify reduction operations."; write," 6. data=...,area_wt1=...,z_bot1=... may be used to modify the data, area-weight, bottom Z arrays"; write," 7. x0=...,y0=...,z0=...,xint0=...,yint0=...,zint0=...,ilabel0=...,iparam0=...,"; write," may be used to modify the full domain grid details"; write," 8. extra_atts=['attribute','data_type',...] adds extra attributes to the slab."; write," See also: hget, hsub, hcat, hsave, hattr, hdata, ncopyatt"; write,""; write," Usage: hcopy, slab1, slab2, overwrite=1, index1=index1, index2=index2, structure0='HYPERSLAB1.0_...', data=..., ..."; return timer_return(func_name); } if (param_set(overwrite)) { // Overwrite slab recursively if (!is_null(index2)) error, "index2=... incompatible with overwrite"; tem_slab= NULL; hcopy, slab1, tem_slab, index1=index1, x1=x1, y1=y1, z1=z1, time1=time1, date1=date1, ilabel1=ilabel1, iparam1=iparam1, data=data, area_wt1=area_wt1, z_bot1=z_bot1, missing_value=missing_value, x0=x0, y0=y0, z0=z0, xint0=xint0, yint0=yint0, zint0=zint0, ilabel0=ilabel0, iparam0=iparam0, structure0=structure0, reset_vars=reset_vars, extra_atts=extra_atts, area_wt_dims=area_wt_dims, z_bot_dims=z_bot_dims, is_present=is_present, is_reduced=is_reduced; slab2= tem_slab; return timer_return(func_name); } if (is_null(index1) && (!is_null(slab1)) && \ (!is_scalar(slab1)) ) { // Array of source hyperslabs; handle recursively if (!is_null(index2)) error, "Index incompatible with array copy"; if (is_null(index2) && (!is_null(slab2)) && \ (!is_scalar(slab2)) ) { // Array-to-array copy if (!dim_conform(dimsof(slab1), dimsof(slab2), trim=1)) error, "Non-conforming dimensions for array-to-array copy"; for (i=I0; i <= I0+numberof(slab1)-1; i++) { hcopy, slab1, slab2, index1=i, index2=i, x1=x1, y1=y1, z1=z1, time1=time1, date1=date1, ilabel1=ilabel1, iparam1=iparam1, data=data, area_wt1=area_wt1, z_bot1=z_bot1, missing_value=missing_value, x0=x0, y0=y0, z0=z0, xint0=xint0, yint0=yint0, zint0=zint0, ilabel0=ilabel0, iparam0=iparam0, structure0=structure0, reset_vars=reset_vars, extra_atts=extra_atts, area_wt_dims=area_wt_dims, z_bot_dims=z_bot_dims, is_present=is_present, is_redu