require, "netcdf.i" func ngv (variable, nc_name, nc_meta, swapyz=, x=, y=, z=, t=, dates=, av=, addtbase=) /* DOCUMENT nc_get_var, variable or nc_get_var, variable, nc_name or nc_get_var, variable, nc_name, nc_meta Returns a structure containing the data for the specified variable along with its attributes and coordinate lists from a netCDF file. If present, nc_name specifies the name of a previously opened file object, otherwise, the default file object is "file". If present, nc_meta is the structure containing the metadata from the netCDF file. Otherwise, the global variable nc_file is used, providing the metadata from the last opened file. The netCDF file will typically have been opened by: nc_name = openb (ncfile_name); nc_meta = nc_file; 3D variables are sometimes ordered (x,z,y) on the input file. This condition is tested for by checking the name of the 3rd dimension for the string "lat" (case insensitive). In this case, the 2nd and 3rd dimensions are transposed to give (x,y,z) data. Keywords: swapyz void, transpose if 3rd dimension name contains "lat" == "no", transposing is disabled entirely. != "no", transposing is performed regardless of dimension names x = value, y = [value], z = [first, last], t= [list], dates= [list], dates takes precedence over t if both specified coordinate values of points in hyperslice to be returned. 1 value (x=120), will return the closest point. 2 values (y = [-30.5, 10.2]) will return all points between the specified values. 3 or more values return the closest point to each value. If a coordinate list is specified (3 or more values), the corresponding dimension will be ordered exactly as requested. Otherwise, the order in the input file will be preserved. av = [0,0,0,1] a vector of dimensions to average. Nonzero indicates that the average in the corresponding dimension should be returned. Should only be used for dimensions which do not require masking or weighting. e.g. a zonal average with no missing points, or most time averges (since a time independent mask can be applied afterward). addtbase != 0 means add the base time, from the time units attribute, to the output time list and reset the units attribute to "days since 00-00-00 00:00:00" */ { /* set the default file and metadata, if they were not passed */ if (is_void(nc_name)) nc_name = file; if (!is_stream(nc_name)){ print, "not a valid file object"; return []; } if (is_void(nc_meta)) nc_meta = nc_file; /* find the index of the variable name */ vars = *nc_meta.vars; ivar = where(vars.name == variable); if (numberof(ivar) == 0) { print, variable + " not found in file " + nc_meta.filename; print, "the available variables are:", vars.name; return []; } ivar = ivar(1); // find the coordinate lists of the variable // assume that only multidimensional variables will be extracted dimlist = *nc_meta.dims; dims = *vars(ivar).dimlist; ndim = numberof(dims); if (ndim < 2 || ndim > 5) { print, "The variable must have 2, 3, or 4 dimensions, not", ndim; return []; } dims += 1; // if there is a record dimension, check if this is a record variable Record = 0; unlimited = (dimlist.size==0); if (noneof(unlimited)) { unlimited = -1; } else { unlimited = where(unlimited)(0); Record = dims(1) == unlimited; } /* set the dimension names */ local dimx, dimy, dimz, dimt; dimx = dimlist(dims(0)).name; dimy = dimlist(dims(-1)).name; dimz = dimt = []; if (Record) { if (ndim == 4) dimz = dimlist(dims(-2)).name; dimt = dimlist(dims(1)).name; } else { if (ndim == 3) dimz = dimlist(dims(-2)).name; } // Check for a time variable in a file with no records and set dimt if (unlimited < 0) { i = where(vars.name == "time"); if (numberof(i) != 0) dimt = "time"; } /* check for transposition required */ if (!is_void(dimz)) { if (is_void(swapyz)) { swapyz = strmatch(dimz, "lat", 1); } else { swapyz = !(swapyz == "no"); } } else { swapyz = 0; } if (swapyz) { tmpd = dimy; dimy = dimz; dimz = tmpd; } /* set averaging flags for each dimension */ if (is_void(av)) av = [0, 0, 0, 0]; xav = av(1); yav = av(2); zav = av(3); tav = av(4); // Set the time list if (is_void(t) || is_void(dimt)) { tlist = 1; tl = 0; // set the current time if dimt defined if (!is_void(dimt)) tl = get_member(nc_name,dimt); } else { // make a list of available times if (unlimited >=0 && (dimt=="time")) tl = get_times (nc_name); else tl = collect(nc_name,dimt); // select the requested times for the available list if (numberof(t) == 1 ) { tlist = abs(tl - t)(mnx); } else if (numberof(t) == 2 ) { if (numberof(tl) == 1) { tlist = 1; } else { tmin = t(min); tmax = t(max); tlist = where ((tmin<=tl) & (tl <= tmax)); } } else { t = t(*)(-:1:numberof(tl),); tlist= abs(t-tl)(mnx,); } tl = tl(tlist); } ntl = numberof(tl); nt = ntl if (tav) nt = 1; /* set the dimension lists for the spatial dimensions */ if (is_void(x) && is_void(y) && is_void(z)) everything = 1; nx = ncsetdims(nc_name,x,dimx,xav,xl,xrng,xlist); ny = ncsetdims(nc_name,y,dimy,yav,yl,yrng,ylist); if (is_void(dimz)) { zl = [1000.]; zrng = 1; zlist = [1]; nz = 1; } else { nz = ncsetdims(nc_name,z,dimz,zav,zl,zrng,zlist); } /* check for a date variable */ date = array(double, ntl); i = where(vars.name == "date"); isdate = (numberof(i) != 0); if (!isdate) { daymonth=[31,28,31,30,31,30,31,31,30,31,30,31]; dayyear = daymonth(cum); } i = where(vars.name == "datesec"); isdtsc = (numberof(i) != 0); /* x weights, don't really know how to do this. Should use grid edges, if available. Maybe later */ xw = array(1., numberof(xl)); xw /= sum(xw); /* y weights, if this is a latitude axis, get gauss weights if available, otherwise generate a list a delta sin(lat) if not latitude, use uniform weights */ if ((numberof(yl) > 1) && (strmatch(dimy, "lat", 1))) { igw = where(vars.name == "gw"); if (numberof(igw) != 0) { yw = get_member (nc_name, "gw")(ylist); } else { yw = sin(yl(pcen)*pi/180)(dif); } } else { yw = array(1., numberof(yl)); } yw /= sum(yw); /* z weights, don't really know how to do this. Should use grid edges, if available. Maybe later */ if (numberof(zl) > 1) { zw = zl(pcen)(dif); zw /= sum(zw); } else { zw = [1.]; } /* exchange y,z lists, if required */ if (swapyz) { tmp = ylist; ylist = zlist; zlist = tmp; } /* set number of times and initialize array for data */ slice = array (float, nx, ny, nz, nt); /* time loop for extracting data */ for (n=1; n<=ntl; n++) { // position to the required record if (Record) jr, nc_name, tlist(n); /* get the date or make one up */ if (isdate) { date(n) = get_member (nc_name, "date"); if (isdtsc) date(n) += get_member (nc_name, "datesec")/86400.; } else { // this should account for the time units (days since ...) dy = tl(n)%365; month = where(dy<=dayyear)(1); daym = dy + 1 - dayyear(month); date(n) = 10000*int(tl(n)/365.) + month*100 + daym; } /* if all dimension specifications are ranges, specify endpoints, otherwise use the full index list in each coordinate */ if (everything) { data = get_member(nc_name, variable); } else if (xrng && yrng && zrng) { if (is_void(dimz)) { data = get_member (nc_name, variable) (xlist(1):xlist(0), ylist(1):ylist(0)); data = data(,,-); } else { data = get_member (nc_name, variable) (xlist(1):xlist(0), ylist(1):ylist(0), zlist(1):zlist(0)); } } else { data = get_member (nc_name, variable)(xlist, ylist, zlist); } if (swapyz) data = transpose(data,[2,3]); /* fun with averaging */ if (noneof(av)) { // no averaging at all slice(,,,n) = data; } else if (noneof(av(1:3))) { // time averaging only slice += data; } else { // space and/or time averaging if (xav) data = (data(+,,)*xw(+))(-,,); if (yav) data = (data(,+,)*yw(+))(,-,); if (zav) data = (data(,,+)*zw(+))(,,-); if (tav) { slice += data; } else { slice(,,,n) = data; } } } if (tav) slice /= ntl; /* set up the structure */ nc_variable = NC_variable (name = variable, filename=nc_meta.filename, t = &tl, date = &date, x = &xl, y = &yl, z = &zl, xw= &xw, yw= &yw, zw= &zw, xname = dimx, yname = dimy, data = &slice); /* vertical coordinate info */ if (!is_void(dimz)) { nc_variable.zname = dimz; if (dimz == "lev" || dimz == "ilev") { i = where(vars.name == "hyai"); if (numberof(i) != 0) nc_variable.hyai = &get_member(nc_name, "hyai"); i = where(vars.name == "hybi"); if (numberof(i) != 0) nc_variable.hybi = &get_member(nc_name, "hybi"); i = where(vars.name == "hyam"); if (numberof(i) != 0) nc_variable.hyam = &get_member(nc_name, "hyam"); i = where(vars.name == "hybm"); if (numberof(i) != 0) nc_variable.hybm = &get_member(nc_name, "hybm"); } i = where(vars.name == "P0"); if (numberof(i) != 0) nc_variable.P0 = &get_member(nc_name, "P0"); } else { nc_variable.zname = "none"; } /* set other variable attributes */ attrs = *vars(ivar).attrs; attr = get_attr (attrs, "long_name") if (!is_void(attr)) nc_variable.long_name = attr; attr = get_attr (attrs, "units") if (!is_void(attr)) nc_variable.units = attr; attr = get_attr (attrs, "history") if (!is_void(attr)) nc_variable.history = attr; else nc_variable.history = " "; attr = get_attr (attrs, "missing_value") if (!is_void(attr)) nc_variable.spv = &attr; /* set global attributes */ attrs = *nc_meta.attrs; attr = get_attr (attrs, "case") if (!is_void(attr)) nc_variable.case = attr; attr = get_attr (attrs, "title") if (!is_void(attr)) nc_variable.title = attr; /* set the time attributes */ if (!is_void(dimt)) { nc_variable.tname = dimt; ivar = where(vars.name == dimt)(0); attrs = *vars(ivar).attrs; attr = get_attr (attrs, "long_name"); if (!is_void(attr)) nc_variable.tlname = attr; attr = get_attr (attrs, "units"); if (!is_void(attr)) nc_variable.tunits = attr; // reset time values and units if addtbase if (!is_void(addtbase) && addtbase) { if (!is_void(attr)) { year= month= day= hour= minute= second= 0; sread, attr, format="days since %d%d%d%d:%d:%d%d", year, month, day, hour, minute, second; month = -month; day = -day; daymonth=[31,28,31,30,31,30,31,31,30,31,30,31]; dayyear = daymonth(cum); dayfrst = grow(0.,dayyear); toff = year*365 + dayfrst(month) + day-1. + (hour*3600. + minute*60. + second) / 86400.; nc_variable.tunits = "days since 00-00-00 00:00:00"; nc_variable.t = &(tl+toff); } } } /* set the x attributes */ ivar = where(vars.name == dimx)(0); attrs = *vars(ivar).attrs; attr = get_attr (attrs, "long_name"); if (!is_void(attr)) nc_variable.xlname = attr; attr = get_attr (attrs, "units"); if (!is_void(attr)) nc_variable.xunits = attr; /* set the y attributes */ ivar = where(vars.name == dimy)(0); attrs = *vars(ivar).attrs; attr = get_attr (attrs, "long_name"); if (!is_void(attr)) nc_variable.ylname = attr; attr = get_attr (attrs, "units"); if (!is_void(attr)) nc_variable.yunits = attr; /* set the z attributes */ if (!is_void(dimz)) { ivar = where(vars.name == dimz)(0); attrs = *vars(ivar).attrs; attr = get_attr (attrs, "long_name"); if (!is_void(attr)) nc_variable.zlname = attr; attr = get_attr (attrs, "units"); if (!is_void(attr)) nc_variable.zunits = attr; } return nc_variable; } struct NC_variable { string case, title, filename, name, long_name, units, history; string xname, yname, zname, tname; string xlname, ylname, zlname, tlname; string xunits, yunits, zunits, tunits; pointer x, y, z, t, date, xw, yw, zw, hyam, hyai, hybm, hybi, P0, spv, data; } func ncsetdims(nc_name, x, dimx, xav, &xl, &xrng, &xlist) { xl = get_member (nc_name, dimx); xrng = 1; xlist = indgen(numberof(xl)); if (!is_void(x)) { if (numberof(x) == 1 ) { xlist = abs(xl - x)(mnx); } else if (numberof(x) == 2 ) { if (numberof(xl) == 1) { xlist = 1; } else { xmin = x(min); xmax = x(max); xlist = where ((xmin<=xl) & (xl <= xmax)); } } else { xrng = 0; // really have a discontinous list, not a range x = x(*)(-:1:numberof(xl),); xlist= abs(x-xl)(mnx,); } xl = xl(xlist); } nx = numberof(xl); if (xav) nx = 1; return nx; } func get_attr (attrs, atrname) /* DOCUMENT get_attr, attrs, atrname Return the value of the attribute given by atrname. This function is used primarily to convert char arrays to a string. */ { nattr = where(attrs.name == atrname); if (numberof(nattr) == 0) return []; nattr = nattr(1); attr = *attrs(nattr).data; if (typeof(attr) == "char") { attr = string(attrs(nattr).data) } return attr; } func nmin(var) /* DOCUMENT val = nmin(var) Return the minimum of the data segment of the structure var, excluding special values */ { /* set contour list for difference */ if (is_void(*var.spv)) { dmin = min(*var.data); } else { list = where(*var.data != *var.spv); dmin = min((*var.data)(list)); } return dmin; } func nmax(var) /* DOCUMENT val = nmax(var) Return the maximum of the data segment of the structure var, excluding special values */ { /* set contour list for difference */ if (is_void(*var.spv)) { dmax = max(*var.data); } else { list = where(*var.data != *var.spv); dmax = max((*var.data)(list)); } return dmax; } func ndif (v1, v2, name=, units=) /* DOCUMENT v = ndif(v1, v2, name=NAME, units=UNITS) Return the difference (v1-v2) of the data in the structures v1 and v2 obsolete SEE: nop */ { v = nop(v1,v2,"-",name=name,units=units); return v; } func nadd (v1, v2, name=, units=) /* DOCUMENT v= nadd(v1, v2, name=NAME, units=UNITS) Return the sum (v1+v2) of the data in the structures v1 and v2 obsolete SEE: nop */ { v = nop(v1,v2,"+",name=name,units=units); return v; } func nmul (v1, v2, name=, units=) /* DOCUMENT v = nmul(v1, v2, name=NAME, units=UNITS) Return the product (v1*v2) of the data in the structures v1 and v2 obsolete SEE: nop */ { v = nop(v1,v2,"*",name=name,units=units); return v; } func nop (v1, v2, op, name=, long_name=, units=) /* DOCUMENT v = nop(v1, v2, op, name=NAME, units=UNITS, long_name=LNM) Return the structure (v1 "op" v2) where at least one of v1 and v2, or both, are structures of type NC_variable. The arithmetic operator, op, must be in ["+", "-", "*", "/", "^"]. NAME and UNITS are new values for the output variable */ { /* make sure that an input variable is a proper structure */ v1_struct = (nameof(structof(v1)) == "NC_variable"); v2_struct = (nameof(structof(v2)) == "NC_variable"); if (!(v1_struct+v2_struct)) { print, "neither input variable is a structure of type NC_variable"; info, v1; info, v2; return; } /* get the data arrays */ if (v1_struct) { d1 = v1.data; v1spv = *v1.spv } else { d1 = &v1; } if (v2_struct) { d2 = v2.data; v2spv = *v2.spv } else { d2 = &v2; } /* get the dimensions and expand to 4 if necessary */ dim1 = dimsof(*d1); dim2 = dimsof(*d2); while (dim1(1) < 4) {dim1 = grow(dim1, 1); dim1(1)++;} while (dim2(1) < 4) {dim2 = grow(dim2, 1); dim2(1)++;} // Check dimensions and copy the larger input structure, // if only one is a structure, then its dimensions must be the largest // Also set lists of points where the input structures have special values if (!v2_struct) { if (nallof(dim1 >= dim2)) goto nonconf; vm = v1; if (!is_void(v1spv)) list1 = where(*d1 == v1spv); } else if (!v1_struct) { if (nallof(dim1 <= dim2)) goto nonconf; vm = v2; if (!is_void(v2spv)) list2 = where(*d2 == v2spv); } else if (allof(dim1 == dim2)) { vm = v1; if (!is_void(v1spv)) list1 = where(*d1 == v1spv); if (!is_void(v2spv)) list2 = where(*d2 == v2spv); if (is_void(v1spv) && !is_void(v2spv)) vm.spv = &v2spv; } else if (allof(dim1 >= dim2)) { vm = v1; if (!is_void(v1spv)) list1 = where(*d1 == v1spv); if (!is_void(v2spv)) { list2 = where(*d2 == v2spv); if (numberof(list2) > 0) { print, "the smaller variable may not have special values"; error; } } } else if (allof(dim1 <= dim2)) { vm = v2; if (!is_void(v1spv)) { list1 = where(*d1 == v1spv); if (numberof(list1) > 0) { print, "the smaller variable may not have special values"; error; } } if (!is_void(v2spv)) list2 = where(*d2 == v2spv); } else goto nonconf; /* operate on the data */ if (op == "+") result = *d1 + *d2; else if (op == "-") result = *d1 - *d2; else if (op == "*") result = *d1 * *d2; else if (op == "/") { listz = where(*d2 == 0.); if (numberof(listz) > 0) { if (is_void(*vm.spv)) { spv = 1.e30; vm.spv = &spv; } list2 = grow(list2,listz); tmp = *d2; tmp(listz) = 1.; result = *d1 / tmp; } else { result = *d1 / *d2; } } else if (op == "^") result = *d1 ^ *d2; else {print, "unknown operator", op; error;} // set special values where either input set has them if (numberof(list1) > 0) result(list1) = *vm.spv; if (numberof(list2) > 0) result(list2) = *vm.spv; vm.data = &result; if (!is_void(name)) vm.name = name; if (is_void(long_name)) { if (v1_struct && v2_struct) vm.long_name = v1.name + " " + op + " " + v2.name; } else { vm.long_name = long_name; } if (!is_void(units)) vm.units=units if ((v1_struct && v2_struct) && (v1.case != v2.case )) vm.case = v1.case + " " + op + " " + v2.case; if (v1_struct) { hist = v1.filename + "(" + v1.name + ")"; hatt = "v1 history" + v1.history + "\n"; } else { hist = pr1(*d1); hatt = " "; } hist = hist + " " + op + " "; if (v2_struct) { hist = hist + v2.filename + "(" + v2.name + ")\n"; hatt = hatt + "v2 history" + v2.history+ "\n"; } else { hist = hist + pr1(*d2); } vm.history = hist + hatt; return vm; /* nonconforming data */ nonconf: print,"input variables not comformable",dim1,dim2; error; } func nzavdev (v1, name=, long_name=) /* DOCUMENT v = nzavdev(v1, name=NAME, long_name=LNM) Return a structure with the deviations of the v1 from the zonal (x) mean, where the mean is taken over all valid points. v1 and v are structures of type NC_variable. NAME and UNITS are new values for the output variable */ { /* make sure that an input variable is a proper structure */ v1_struct = (nameof(structof(v1)) == "NC_variable"); if (!(v1_struct)) { print, "input variable is not a structure of type NC_variable"; info, v1; return; } // get the data array and special value d1 = *v1.data; dims = dimsof(d1); nx = dims(2); if (nx == 1) { print, "input variable has a degenerate first dimension (dims)", dims; return; } spv = *v1.spv; n = dims(3)*dims(4)*dims(5); // compute deviation from zonal mean if (is_void(spv)) { zav = d1(avg,..); d1 = d1 - zav(-:1:nx,..); } else { for (i=1; i<=n; i++) { list = where(d1(,i) != spv); if (numberof(list)>0) { zav = sum(d1(list,i))/numberof(list); d1(,i) = d1(,i) - zav; } } } vm = v1; vm.data = &d1; if (!is_void(name)) vm.name = name; if (is_void(long_name)) { vm.long_name = vm.long_name + " zonal deviation"; } else { vm.long_name = long_name; } return vm; } func ncat (v1, v2) /* DOCUMENT v = ncat(v1, v2) Return the concatenation (in time) of v1 and v2, where both are structures of type NC_variable. */ { /* make sure that the input variables are proper structures */ v1_struct = (nameof(structof(v1)) == "NC_variable"); v2_struct = (nameof(structof(v2)) == "NC_variable"); if (!v1_struct || !v2_struct) { print, "the input variables are not structures of type NC_variable"; info, v1; info, v2; exit; } /* check that the spatial dimension match (do not check coordinate lists) */ d1 = dimsof(*v1.data); d2 = dimsof(*v2.data); if (nallof(d1(1:4) == d2(1:4))) { print, "the input variables do not have the same spatial dimensions"; print, d1, d2; exit; } /* copy the first input structure and grow its data and time variables */ v = v1; v.data = &grow(*v1.data, *v2.data); v.t = &grow(*v1.t, *v2.t); v.date = &grow(*v1.date, *v2.date); return v; } func ntimavg (v1, navg, name=, long_name=, units=, op=) /* DOCUMENT v = ntimavg(v1, navg, name=NAME, long_name=LNM, units=UNITS) Return the structure v1, with the record dimension replaced by the means over navg samples. The length of the input record dimension (nt) must be a multiple of navg and the output length is nt/navg. op="avg","min","max" speficies the operator to use (default "avg") NAME, LNM and UNITS are new values for the output variable */ { /* make sure that an input variable is a proper structure */ v1_struct = (nameof(structof(v1)) == "NC_variable"); if (!(v1_struct)) { print, "input variable is not a structure of type NC_variable"; info, v1; return; } if (is_void(op)) op = "avg"; /* check the time dimension and averaging interval */ dims= dimsof(*v1.data); nt = dims(0); if (is_void(navg)) navg = 12; if (nt%navg != 0) { print, "number of time samples not a multiple of averaging interal", nt, navg; error; } /* set time weights, since ngv does not deal with month lengths */ no = nt/navg; if (navg == 12) { for (iy=1;iy<=no;iy++) { twt = grow(twt,[31,28,31,30,31,30,31,31,30,31,30,31]); } } else { twt = array(1.,nt); } /* do time averages */ date1 = *v1.date; t1 = *v1.t; d1 = *v1.data; dims(0) = no; d2 = array(double,dims); t2 = array(double,no); date2 = array(double,no); tt1 = 1; tt2 = navg; for (n=1; n<=no; n++) { // d2 = grow(d2,(d1(*,tt1:tt2)*twt(tt1:tt2))(..,sum) / twt(sum:tt1:tt2)) ; // d2 = grow(d2,d1(..,avg:tt1:tt2)); if (op == "avg") { d2(..,n) = d1(..,avg:tt1:tt2); } else if (op == "min") { d2(..,n) = d1(..,min:tt1:tt2); } else if (op == "max") { d2(..,n) = d1(..,max:tt1:tt2); } else { print, "unrecognized operator"; error; } t2(n) = t1(avg:tt1:tt2); date2(n) = date1(tt2); tt1 += navg; tt2 += navg; } /* set the output variable */ v2 = v1; v2.data = &d2; v2.date = &date2; v2.t = &t2; return v2; } func nrunavg (v1, navg, name=, long_name=, units=) /* DOCUMENT v = nrunavg(v1, navg, name=NAME, long_name=LNM, units=UNITS) Return the structure v1, with the record dimension replaced by the means over navg samples. The length of the output time dimension is the same as that of the input. The first navg/2 and last navg/2 values are the same. NAME, LNM and UNITS are new values for the output variable */ { /* make sure that an input variable is a proper structure */ v1_struct = (nameof(structof(v1)) == "NC_variable"); if (!(v1_struct)) { print, "input variable is not a structure of type NC_variable"; info, v1; return; } /* check the time dimension and averaging interval */ dims= dimsof(*v1.data); nt = dims(0); if (is_void(navg)) { print, "navg must be specified"; error; } no = nt; /* do time averages */ date1 = *v1.date; t1 = *v1.t; d1 = *v1.data; dims(0) = no; d2 = array(double,dims); tt1 = 1; tt2 = navg; for (n=navg/2+1; n<=no-navg/2-1; n++) { d2(..,n) = d1(..,avg:tt1:tt2); tt1 += 1; tt2 += 1; } d2(..,1:navg/2) = d2(..,navg/2+1); d2(..,no-navg/2:no) = d2(..,no-navg/2-1); /* set the output variable */ v2 = v1; v2.data = &d2; return v2; } func trscal (var, fac=, units=) { /* DOCUMENT trscal(var) convert trace gase values to volume mixing ratio and rescale to appropriate units. If fac is nonvoid, use it directly for scaling. */ local name; local unit, mw, fac; /* set molecular weight and unit scaling factor according to name */ name = var.name; if (is_void(fac)) { if (name == "Q") { mw = 18.; fac = 1.E6; units = "ppmv"; } else if (name == "CH4") { mw = 16.; fac = 1.E6; units = "ppmv"; } else if (name == "N2O") { mw = 44.; fac = 1.E9; units = "ppbv"; } else if (name == "CFC11") { mw = 136.; fac = 1.E12; units = "pptv"; } else if (name == "CFC12") { mw = 120.; fac = 1.E12; units = "pptv"; } else if (name == "O3VMR") { fac = 1.E6; units = "ppmv"; } else { return var; } if (!is_void(mw)) fac *= 29./mw; } vs = nop(var,fac,"*",units=units); return vs; } func month_list(file, frst=, last=) /* DOCUMENT list = month_list(file, frst=first, last=last) return a 2 dimensional array of lists of time values corresponding to the 12 months, where month number is the second dimension. Beginning with time index first and end with index last. If all of the 12 lists are not the same length, abort with an error; */ { t=get_times(file); dates = long(collect(file, "date")); if (!is_void(frst) || !is_void(last)) { t = t(frst:last); dates = dates(frst:last); } days = dates%100; // check for monthly average ccm output ccmav = allof(days==1); // with a date of the 1st of next month if (ccmav) ccmav = nallof(t(2:)==t(:-1)+1) months = (dates/100)%100; if (ccmav) { months = months -1; months(where(months==0))=12; } for (m=1; m<=12; m++) { list = where(months==m); if (m==1) { mlist = list(,-); nmon = numberof(list); } else { if (numberof(list) != nmon) { print, "month_list: number of months not all the same"; print, nmon, numberof(list), "for month ", m; error; } grow,mlist,list; } } return t(mlist); } func season_list(file, &djf, &mam, &jja, &son, &t, frst=, last=) /* DOCUMENT season_list,file, DJF, MAM, JJA, SON, t, frst=first return lists of time values corresponding to the 4 seasons, DJF, MAM, JJA, SON, and all times in the file, beginning with time index first. */ { t=get_times(file); dates = long(collect(file, "date")); if (!is_void(frst) || !is_void(last)) { t = t(frst:last); dates = dates(frst:last); } days = dates%100; // check for monthly average ccm output ccmav = allof(days==1); // with a date of the 1st of next month months = (dates/100)%100; if (ccmav) { months = months -1; months(where(months==0))=12; } djflist = where(months==12 | months==1 | months==2); mamlist = where(months==3 | months==4 | months==5); jjalist = where(months==6 | months==7 | months==8); sonlist = where(months==9 | months==10 | months==11); djf = t(djflist); mam = t(mamlist); jja = t(jjalist); son = t(sonlist); if (!is_void(verbose) && (verbose == 1)) { print, "-------------dates:-------------"; print, dates; print, "djf", dates(djflist); print, "mam", dates(mamlist); print, "jja", dates(jjalist); print, "son", dates(sonlist); } return; } func season_tlist(file, months, frst=, last=) /* DOCUMENT list = season (file, [mm,mm,...], frst=first,last=last) return a list of time values corresponding to the months specified, for all times in the file, beginning with time index "first" and ending with time index "last" */ { list = month_list(file,frst=frst,last=last); list = list(,months)(*); list = list(sort(list)); return list; } func chkfiles(Flist_old,Flist_new,dirchk=) { /* DOCUMENT Flist_new = chkfiles(Flist_old, Flist_new,dirchk=) check two list of files, that either 1. none of the files in Flist_old appear in Flist_new, or 2. the lists overlap (possibly completely) return the list of new files (those elements of Flist_new which are not in Flist_old). return void if no new files. return false (0) for mismatched lists if dirchk!=0 the full directory path must match, otherwise, only the filename is checked (the default) */ nold = numberof(Flist_old); nlst = numberof(Flist_new); // extract the file names (trim directories) from the file lists local name, name_old, name_new; if (is_void(dirchk) || !dirchk) { for (n=1;n<=numberof(Flist_old);n++) { name = Flist_old(n); do { str=strtok(name,"/"); if(strlen(str(2))!=0) name=str(2); } while (strlen(str(2))!=0); name_old = grow(name_old,name); } for (n=1;n<=numberof(Flist_new);n++) { name = Flist_new(n); do { str=strtok(name,"/"); if(strlen(str(2))!=0) name=str(2); } while (strlen(str(2))!=0); name_new = grow(name_new,name); } } else { name_new = Flist_new; name_old = Flist_old; } // check for whether the lists overlap first = where(name_new(1) == name_old); /* if no overlap, return the new list */ if (numberof(first) == 0) return Flist_new; /* overlap, make sure all files match through end of old list*/ name_ovr = name_old(first(1):); novr = numberof(name_ovr); if ((novr > nlst) || (anyof(name_new(1:novr) != name_ovr))) { print, "old part of list does not match"; Flist_old; Flist_new; name_ovr; return 0; } else { /* determine whether there are any new files */ nnew = nlst - novr; if (nnew > 0) { /* return the list of new files */ return Flist_new(novr+1:); } else { /* no new files, return void */ return []; } } } func hyb_vinteg(vhyb, ps) /* DOCUMENT v = hyb_vinteg (vhyb, ps) return the vertically integrated values of variable vhyb, which must be in ccm hybrid coordinates. The surface pressure, ps, must be provided. vhyb, ps, and the result are structures of type NC_variable. */ { /* check the input data */ if (nameof(structof(vhyb)) != "NC_variable") { print, "input variable not a structure of type NC_variable"; error; } if (nameof(structof(ps)) != "NC_variable") { print, "input surface pressure not a structure of type NC_variable"; error; } /* set the output coordinate list */ zlo = [1000.]; zwo = [1.]; /* get array of pressures of hybrid surfaces for integrals midpoint variable: use interfaces interface variable: use extended midpoints */ local PS, A, B, P; PS = *ps.data; if (vhyb.zname == "lev") { A = *vhyb.hyai; B = *vhyb.hybi; } else if (vhyb.zname == "ilev") { A = grow((*vhyb.hyai)(1), *vhyb.hyam, (*vhyb.hyai)(0)); B = grow((*vhyb.hybi)(1), *vhyb.hybm, (*vhyb.hybi)(0)); } else { print, "hyb_vinteg, not a hybrid coordinate variable, zname=", vhyb.zname; error; } P0 = *vhyb.P0; P = hybridp(PS,A,B,P0); if (is_void(P)) return; DP = P(,,dif,); /* do the vertical integral */ vint = (*vhyb.data * DP)(,,sum,)(,,-,); /* create the output structure */ vprs = vhyb; vprs.long_name = vprs.long_name + " vertical integral" vprs.data = &vint; vprs.z = &zlo; vprs.zw = &zwo; vprs.zname = "pressure"; vprs.zlname= "pressure"; vprs.zunits= "mb"; vprs.hyam = &[]; vprs.hyai = &[]; vprs.hybm = &[]; vprs.hybi = &[]; vprs.spv = &spv; return vprs; } func hyb_vinterp(vhyb, ps, vlist, T=, z=, spv=) /* DOCUMENT v = hyb_vinterp (vhyb, ps, vlist,spv=SPV) or v = hyb_vinterp (vhyb, ps, vlist, T=temperature) or v = hyb_vinterp (vhyb, ps, vlist, z=VCOORD) return the vertically interpolated values of variable vhyb at the vertical locations given by vlist. vhyb must be in ccm hybrid coordinates. The surface pressure, ps must be provided, unless VCOORD is nonnil. If VCOORD is specified, it gives the coordinate values at the input coordinate surfaces and will be used directly for interpolation. Otherwise, if T is not void, it gives the temperature at the input coordinate surfaces and will be converted to potential temperature for interpolation. If neither z nor T is specified, then pressure interpolation will be performed. If the vlist has only a single element, then the output field will not have a height dimension variable and the variable name and long_name will be changed by adding the value of the vertical coordinate. e.g. Z3, vlist=500 will return name=Z3-500 long_name=500 mb ... spv is the special value to be used when an output surface is outside of the range of the input surfaces. If SPV is void, then such points will have the value at the last input surface (constant fill rather than extrapolation). The vhyb, ps, T, z, and the result are structures of type NC_variable. */ { /* check the input data */ if (nameof(structof(vhyb)) != "NC_variable") { print, "input variable not a structure of type NC_variable"; error; } if (is_void(z)) { if (nameof(structof(ps)) != "NC_variable") { print, "input surface pressure not a structure of type NC_variable"; error; } } else { if (nameof(structof(z)) != "NC_variable") { print, "input coordinate variable not a structure of type NC_variable"; error; } } if (!is_void(T)) { if (nameof(structof(T)) != "NC_variable") { print, "input temperature not a structure of type NC_variable"; error; } } /* set the output coordinate list */ if (numberof(vlist)==1) { zcoord = 0; zlist = [vlist(1)]; } else { zcoord = 1; zlist = vlist; } Zlist = zlist /* set the vertical coordinate array for intepolation */ if (!is_void(z)) { zc = *z.data; zname = z.name; zlname = z.long_name; zunits = z.units; } else { /* get array of pressures of hybrid surfaces */ local PS, A, B, P; PS = *ps.data; if (vhyb.zname == "lev") { A = *vhyb.hyam; B = *vhyb.hybm; } else if (vhyb.zname == "ilev") { A = *vhyb.hyai; B = *vhyb.hybi; } else { print, "hyb_vinterp, not a hybrid coordinate variable, zname=", vhyb.zname; error; } P0 = *vhyb.P0; zc = hybridp(PS,A,B,P0); if (is_void(zc)) return; if (!is_void(T)) { /* convert temperature to potential temperature */ zc = theta(*T.data, zc); zname = "theta"; zlname = "potential temperature"; zunits = "K"; } else { /* pressure: rescale the coordinate list to Pa */ Zlist = Zlist * 100; zname = "pressure"; zlname = "pressure"; zunits = "mb"; } } /* do the vertical interpolation */ vprsd = interp3(*vhyb.data, zc, Zlist, spv=spv); /* create the output structure */ vprs = vhyb; vprs.data = &vprsd; if (zcoord) { vprs.zname = zname; vprs.zlname= zlname; vprs.zunits= zunits; vprs.z = &zlist; vprs.zw = &zlist; } else { vprs.name = vprs.name + "-" + pr1(zlist(1)); vprs.long_name = pr1(zlist(1))+ " " + zunits + " " + vprs.long_name; vprs.zname = "none" vprs.zlname= ""; vprs.zunits= ""; vprs.z = &zlist; vprs.zw = &zlist; } vprs.hyam = &[]; vprs.hyai = &[]; vprs.hybm = &[]; vprs.hybi = &[]; vprs.spv = &spv; return vprs; } func hybridp (ps, A, B, P0) /* DOCUMENT hybridp (ps, A, B) return the pressures of hybrid surfaces given surface pressure (ps) and hybrid coordinate lists A and B. */ { local dimsp; if (is_void(P0)) P0 = 1.E5; dimsp = dimsof(ps); if (dimsp(1) < 2) { print, "hybridp: surface pressure field must have at least 2 dimensions"; print, "dimsof(ps):", dimsp; return; } nx = dimsp(2); ny = dimsp(3); if (numberof(A) != numberof(B)) { print, "hybridp: hybrid A and B lists do not have same lengths"; print, "numberof(A)",numberof(A),"numberof(B)",numberof(B); return; } P = ps * B(-:1:nx,-:1:ny,) + P0 * A(-:1:nx,-:1:ny,); return P; } func theta (T, P) /* DOCUMENT theta (T, P) return the potential temperature, given temperature T and pressure p (in Pa). */ { local dimsT, dimsP, kappa, P0; kappa = 0.287; P0 = 1.E5; dimsT = dimsof(T); dimsP = dimsof(P); if (anyof(dimsT != dimsP)) { print, "theta: dimensions of T and P do not match"; print, "theta: dimsof(T)", dimsT; print, "theta: dimsof(P)", dimsP; return; } Theta = T * (P0/P)^kappa; return Theta; } func interp3 (V, C, Co, spv=) /* DOCUMENT interp3 (V, C, Co, spv=) interpolate V in the 3rd dimension, where C is the spatially varying list of input coordinates (same dimensions as V) and Co is the one dimensional list of output coordinates. Used for vertical interpolation from a horizontally vary set of input surfaces to a uniform set of output surfaces. If spv is not void, points with output coordinates outside the range of the input coordinates will be set to spv. */ { /* interpolate in 3rd dimension */ dimsV = dimsof (V); dimsC = dimsof (C); if (dimsV(1) != 4) { print, "dimsof(V)", dimsv; return; } if (anyof(dimsV != dimsC)){ print, "dimsof(V)", dimsv; print, "dimsof(C)", dimsc; return; } dimsN = dimsV; dimsN(4) = numberof(Co); N = array(double, dimsN); /* set out of range values to spv */ if (!is_void(spv)) spvset = 1; /* loop over 1st, 2nd and 4th dimensions */ for (m=1; m<=dimsV(5); m++) { for (j=1; j<=dimsV(3); j++) { for (i=1; i<=dimsV(2); i++) { N(i,j,,m) = interp (V(i,j,,m),C(i,j,,m),Co); if (spvset) { list = where(Co > max(C(i,j,,m))); if (numberof(list) > 0) N(i,j,list,m) = spv; list = where(Co < min(C(i,j,,m))); if (numberof(list) > 0) N(i,j,list,m) = spv; } } } } return N; } func ncp(old_fname,new_fname) /* DOCUMENT ncp, old_fname, new_fname Make a copy of a netcdf file, duplicating all variables and records. This function is intended to allow the contents of the variables to be altered, for example, generating a changed initial or boundary condition file for the ccm. */ { /* open the original netcdf file */ old_file = openb(old_fname); nc_meta = nc_file; /* begin creating the file structure */ ncf = nc_create(new_fname); /* put all of the dimensions of the original file into the new file */ dims = *nc_meta.dims; ndim = numberof(dims); for (n=1; n<=ndim; n++){ nc_dimdef, ncf, dims(n).name, dims(n).size; } print, "dimensions", dims.name, dims.size; /* check whether there is an unlimited (history) dimension */ unlimited = (dims.size==0); if (noneof(unlimited)) { unlimited = -1; } else { unlimited = where(unlimited)(0); } print, "unlimited=", unlimited; /* put all of the attributes of the original file into the new file */ attrs = *nc_meta.attrs; na = numberof(attrs); history = 0; for (n=1; n<=na; n++){ name = attrs(n).name; data = *attrs(n).data; if (name == "history") data = data + "\n" + timestamp() + "> ncp, " + old_fname + ", " + new_fname; nc_attrdef, ncf, name, , data; } ncf; /* put all of the variables of the original file into the new file, also locate the record variables */ type_names= ["char", "char", "short", "long", "float", "double"]; vars = *nc_meta.vars nvar = numberof(vars); Record= array(int, nvar); for (iv=1; iv<=nvar; iv++){ vname = vars(iv).name; print, "variable", iv, vname; dimlist = *vars(iv).dimlist; if (is_void(dimlist)) { nc_vardef, ncf, vname, type_names(vars(iv).type); Record(iv)= 0; } else { newlist = (dimlist+1)(0:1:-1); dimnames = dims(newlist).name; print," dimensions", dimnames; Record(iv) = (newlist(0) == unlimited); if(Record(iv)) print," record variable"; nc_vardef, ncf, vname, type_names(vars(iv).type), dimnames=dimnames; } /* add variable attributes */ attrs = *vars(iv).attrs; mm = numberof(attrs); for (m=1; m<=mm; m++){ nc_attrdef, ncf, attrs(m).name, vname, *attrs(m).data; } } /* end the file definition */ new_file = nc_enddef(ncf); /* write all nonrecord variables to file, includes coordinates, except time (date is not a "coordinate") */ list = where(!Record); for (iv=1; iv<=numberof(list); iv++){ vname = vars(list(iv)).name; print, "copy nonRecord", vname; get_member(new_file,vname) = get_member(old_file,vname); } /* check for date on old file */ is_date = numberof(where(vars.name == "date")) > 0; is_dsec = numberof(where(vars.name == "datesec")) > 0; if(is_date) date = collect(old_file,"date"); if(is_dsec) datesec = collect(old_file,"datesec"); /* put all of the records of the original file into the new file */ t = get_times(old_file); nt = numberof(t); print, "times", nt, t; for (it=1; it<=nt; it++){ nc_addrec,new_file,t(it); new_file.time = t(it); if(is_date) new_file.date = date(it); if(is_dsec) new_file.datesec = datesec(it); print, "added record:",it,t(it),date(it),datesec(it); } /* go through the records, copying all of the record variables, calling a function to alter the values if necessary */ list = where(Record); for (it=1; it<=nt; it++){ jt,old_file,t(it); jt,new_file,t(it); print,"time",t(it),date(it),datesec(it); for (iv=1; iv<=numberof(list); iv++){ vname = vars(list(iv)).name; data = get_member(old_file,vname); /* change stuff here */ stat = nchgvars(vname,data); /* save variables */ get_member(new_file,vname) = data; } } return; } func nchgvars(vname,data) { print,vname, dimsof(data), min(data),max(data); return; } func nc_def(fname, .., old_ncf=, old_file=, swapyz=) /* DOCUMENT ncfile = nc_def(fname, var1, ..., old_ncf=OLDNCF, old_file=OLDFILE, swapyz=0/1) Define a new netcdf file to contain the variables var1, ... varN, which must be of type "NC_variable". The input variables must all have the same number of times. A netcdf file must be completely defined before any data is actually written to it. First calls nc_create to initialize the netcdf file, then adds the definitions of the input variables to it, then calls nc_enddef to create the actual file. Following the creation of the file, put data from the input variables into the file. Additional records may be added later using npv calls. If OLDNCF is defined, it is an existing netcdf file definition (metadata), all global attributes and nonrecord variables from the old file will be defined in the new file. If OLDFILE is defined, it must be the existing netcdf file handle (open file) corresponding to OLDNCF. The values of all nonrecord variables in the output file will be copied from this file. swapyz != 0, transpose the y an z dimensions (2 and 3rd) SEE ALSO: nc_create, nc_enddef, ngv, npv */ { ncf = nc_create(fname); history = timestamp() + "> nc_def"; if (is_void(swapyz)) swapyz = 0; if (is_void(verbose)) verbose = 0; /* check the number of input variable arguments */ narg = more_args(); if (narg == 0) { print, "no input arguments"; return; } /* loop over the input variables */ local VARS; for (iarg=1; iarg<=narg; iarg++) { var = next_arg(); /* make sure the variable is a proper structure */ if (nameof(structof(var)) != "NC_variable") { print, "input variable " + pr1(narg) + "not a structure of type NC_variable", nameof(structof(var)); print, var; error; } VARS = grow(VARS, var); if (iarg==1) nt = numberof(*var.t); if (nt != numberof(*var.t)) { print, "not all variables have same number of times"; return; } if (verbose) print, "variable: ", var.name /* put the dimensions of the variable into the new file, both as dimensions and as variables */ dims = *ncf.dims; if (!is_void(dims)) dnames = dims.name; if (noneof(dnames == var.xname)) { if (verbose) print, "adding ", var.xname; nc_dimdef, ncf, var.xname, numberof(*var.x); nc_vardef, ncf, var.xname, "double", dimnames=[var.xname]; nc_attrdef, ncf, "long_name", var.xname, var.xlname; nc_attrdef, ncf, "units", var.xname, var.xunits; } if (noneof(dnames == var.yname)) { if (verbose) print, "adding ", var.yname; nc_dimdef, ncf, var.yname, numberof(*var.y); nc_vardef, ncf, var.yname, "double", dimnames=[var.yname]; nc_attrdef, ncf, "long_name", var.yname, var.ylname; nc_attrdef, ncf, "units", var.yname, var.yunits; nc_vardef, ncf, "gw", "double", dimnames=[var.yname]; } if (noneof(dnames == var.zname) &&(var.zname != "none")) { if (verbose) print, "adding ", var.zname; nc_dimdef, ncf, var.zname, numberof(*var.z); nc_vardef, ncf, var.zname, "double", dimnames=[var.zname]; nc_attrdef, ncf, "long_name", var.zname, var.zlname; nc_attrdef, ncf, "units", var.zname, var.zunits; } if (noneof(dnames == var.tname)) { if (verbose) print, "adding ", var.tname; nc_dimdef, ncf, var.tname, 0; nc_vardef, ncf, var.tname, "double", dimnames=[var.tname]; nc_attrdef, ncf, "long_name", var.tname, var.tlname; nc_attrdef, ncf, "units", var.tname, var.tunits; } if (var.zname == "none") { dimnames = [var.xname, var.yname, var.tname]; } else { if (swapyz) { dimnames = [var.xname, var.zname, var.yname, var.tname]; } else { dimnames = [var.xname, var.yname, var.zname, var.tname]; } } /* add the date */ vars = *ncf.vars; if (noneof(vars.name == "date")) { nc_vardef, ncf, "date", "long", dimnames=["time"]; nc_attrdef, ncf, "long_name", "date", "current date as 6 digit integer (YYMMDD)"; nc_vardef, ncf, "datesec", "long", dimnames=["time"]; nc_attrdef, ncf, "long_name", "datesec", "seconds to complete current date"; } /* add the variable to the file definition */ nc_vardef, ncf, var.name, "float", dimnames=dimnames; /* add the variable's attributes */ nc_attrdef, ncf, "long_name", var.name, var.long_name; nc_attrdef, ncf, "units", var.name, var.units; if (!is_void(*var.spv)) nc_attrdef, ncf, "missing_value", var.name, *var.spv; nc_attrdef, ncf, "history", var.name, var.history; /* add global attributes */ if (iarg == 1) { nc_attrdef, ncf, "case", , var.case; nc_attrdef, ncf, "title", , var.title; } } if (verbose) print, "input variables added"; /* ALL INPUT VARIABLES ADDED */ /* add variables from old file if necessary */ if (!is_void(old_ncf)) { if (verbose) print,"adding variables from:",old_ncf; /* make sure all the dimensions of the old file are in the new file */ dims = *ncf.dims; odims = *old_ncf.dims; ndim = numberof(odims); for (n=1; n<=ndim; n++){ if (noneof(dims.name == odims.name(n))) nc_dimdef, ncf, odims(n).name, odims(n).size; } /* check whether there is an unlimited (history) dimension */ unlimited = (odims.size==0); if (noneof(unlimited)) { unlimited = -1; } else { unlimited = where(unlimited)(0); } // print, "unlimited=", unlimited; /* put all of the global attributes of the old file into the new file */ attrs = *ncf.attrs; oattrs = *old_ncf.attrs; na = numberof(oattrs); for (n=1; n<=na; n++){ name = oattrs(n).name; if (noneof(name == attrs.name)) { data = *oattrs(n).data; // update history attribute if (name == "history") data = data + "\n" + history; nc_attrdef, ncf, name, , data; } } /* put the nonrecord variables of the old file into the new file */ type_names= ["char", "char", "short", "long", "float", "double"]; newvars = *ncf.vars; newnames = newvars.name; vars = *old_ncf.vars; nvar = numberof(vars); Record= array(int, nvar); Skip = array(int, nvar); for (iv=1; iv<=nvar; iv++){ vname = vars(iv).name; if (verbose) print, "variable", iv, vname; dimlist = *vars(iv).dimlist; if (is_void(dimlist)) { dimnames = []; Record(iv)= 0; } else { newlist = (dimlist+1)(0:1:-1); dimnames = odims(newlist).name; if (verbose) print," dimensions", dimnames; Record(iv) = (newlist(0) == unlimited); if(Record(iv)) {if (verbose) print," record variable, skipping";continue;} } if (anyof(vname == newnames)) { if (verbose) print, " already in new file, skipping"; Skip(iv) = 1; continue; } nc_vardef, ncf, vname, type_names(vars(iv).type), dimnames=dimnames; /* add variable attributes */ attrs = *vars(iv).attrs; mm = numberof(attrs); for (m=1; m<=mm; m++){ nc_attrdef, ncf, attrs(m).name, vname, *attrs(m).data; } } } // make sure that the history attribute is defined attrs = *ncf.attrs; name = "history"; if (noneof(name == attrs.name)) nc_attrdef, ncf, name, , history; /* END OF THE FILE DEFINITION */ file = nc_enddef(ncf); set_filesize, file, 0x40000000 /* write all nonrecord variables from old file to new file, includes coordinates, except time (date is not a "coordinate") */ if (!is_void(old_file)) { list = where(!Record); for (iv=1; iv<=numberof(list); iv++){ kv = list(iv); if (!Skip(kv)) { vname = vars(kv).name; if (verbose) print, "copy nonRecord", vname; get_member(file,vname) = get_member(old_file,vname); } } } /* put the nonrecord data from the input variables into the file don't worry about making putting it in many times */ for (iarg=1; iarg<=narg; iarg++) { var = VARS(iarg); get_member(file, var.xname) = *var.x; get_member(file, var.yname) = *var.y; get_member(file, "gw") = *var.yw; if (var.zname != "none") get_member(file, var.zname) = *var.z; } /* make records and write the data */ /* put all of the records into the file, use times from last var */ t = *var.t; fdate = *var.date; date = int(fdate); datesec = int((fdate-date)*86400. + .1); nt = numberof(t); if (verbose) print, "times", nt, t; for (it=1; it<=nt; it++){ nc_addrec,file,t(it); if (verbose) print, "added record:",it,t(it),date(it),datesec(it); for (iarg=1; iarg<=narg; iarg++) { var = VARS(iarg); get_member(file, var.tname) = t(it); get_member(file, "date") = date(it); get_member(file, "datesec") = datesec(it); if (swapyz) { get_member(file, var.name) = transpose((*var.data)(..,it),[2,3]); } else { get_member(file, var.name) = (*var.data)(..,it); } } } return file; } func npv(file, .., swapyz=) /* DOCUMENT npv,file, var1, ..., swapyz=0/1 Add data to a previously defined (nc_def) netcdf file, adding new records, as necessary. The variables var1, ... varN must be of type "NC_variable" and already be defined in the netcdf file. swapyz must be set as in the nc_def call. This really shouldn't be necessary, but I don't want to pass the meta data around, or */ { if (is_void(verbose)) verbose = 0; // check the number of input variable arguments narg = more_args(); if (narg == 0) { print, "no input arguments"; return; } // loop over the input variables for (iarg=1; iarg<=narg; iarg++) { var = next_arg(); // make sure the variable is a proper structure if (nameof(structof(var)) != "NC_variable") { print, "input variable " + pr1(narg) + "not a structure of type NC_variable", nameof(structof(var)); print, var; error; } if (verbose) print, "variable: ", var.name; // set time and date variables t = *var.t; fdate = *var.date; date = int(fdate); datesec = int((fdate-date)*86400. + .1); nt = numberof(t); if (verbose) print, "times", nt, t,date; // loop over times, adding records if necessary and writing data for (it=1; it<=nt; it++){ // check times against existing times in file te = get_times(file); // does this record already exist? if (anyof(t(it)==te)) { // this record already exists, go to it jt, file, t(it); } else { // add new record to file, fill in time and date if (verbose) print, "adding record:",t(it),date(it),datesec(it); nc_addrec,file,t(it); get_member(file, var.tname) = t(it); get_member(file, "date") = date(it); get_member(file, "datesec") = datesec(it); } // add the data to the file if (swapyz) { get_member(file, var.name) = transpose((*var.data)(..,it),[2,3]); } else { get_member(file, var.name) = (*var.data)(..,it); } } } } func horavg (nc_var, xbnd=, ybnd=, mask=) /* DOCUMENT horavg, nc_var Return the horizontal average of the data in an NC_variable structure. If the 3 (z) dimension is not 1, then it is preserved, otherwise the result is a single (time) dimensioned array. KEYWORDS: mask, averaging mask of 0 and 1 matching the 1st 2 dimensions of nc_var. If set, 3 values are returned, 1) the complete average 2) the average where mask =1 3) the average where mask!=1 xbnd, a 2 element array giving the min and max x values to average ybnd, a 2 element array giving the min and max y values to average */ { /* construct averaging weights */ nx = numberof(*nc_var.x); wgts = (*nc_var.yw)(-:1:nx,) * (*nc_var.xw); if (!is_void(xbnd)) { if (numberof(xbnd) != 2) { print, "xbnd must be length 2, xbnd=", xbnd; return; } else { wgts(where(*nc_var.xmax(xbnd)),) = 0; } } if (!is_void(ybnd)) { if (numberof(ybnd) != 2) { print, "ybnd must be length 2, ybnd=", ybnd; return; } else { wgts(,where(*nc_var.ymax(ybnd))) = 0; } } if (!is_void(mask)) { wgtl = wgts* mask; xx = sum(wgtl); if (xx != 0) wgtl /= xx; wgto = wgts*(!mask); xx = sum(wgto); if (xx != 0) wgto /= xx; } wgts /= sum(wgts); var_g = (*nc_var.data * wgts)(sum,sum,,); if (!is_void(mask)) { var_l = (*nc_var.data * wgtl)(sum,sum,,); var_o = (*nc_var.data * wgto)(sum,sum,,); } if (dimsof(var_g)(2) == 1) { var_g = var_g(1,); if (!is_void(mask)) { var_l = var_l(1,); var_o = var_o(1,); } } if (!is_void(mask)) { return [var_g, var_l, var_o]; } else { return var_g; } } func horint (nc_var, xbnd=, ybnd=, mask=, radius=) /* DOCUMENT horint, nc_var Return the horizontal integral of the data in an NC_variable structure. If the 3 (z) dimension is not 1, then it is preserved, otherwise the result is a single (time) dimensioned array. KEYWORDS: mask, averaging mask of 0 and 1 matching the 1st 2 dimensions of nc_var. If set, 3 values are returned, 1) the complete average 2) the average where mask =1 3) the average where mask!=1 xbnd, a 2 element array giving the min and max x values to average ybnd, a 2 element array giving the min and max y values to average radius, the radius of the sphere (default is earth radius, 6.37122e+6); */ { /* construct integration weights for sphere */ if (is_void(radius)) radius = 6.37122e+6; fac = 4*pi*radius^2; nx = numberof(*nc_var.x); wgts = (*nc_var.yw)(-:1:nx,) * (*nc_var.xw) * fac; if (!is_void(xbnd)) { if (numberof(xbnd) != 2) { print, "xbnd must be length 2, xbnd=", xbnd; return; } else { wgts(where(*nc_var.xmax(xbnd)),) = 0; } } if (!is_void(ybnd)) { if (numberof(ybnd) != 2) { print, "ybnd must be length 2, ybnd=", ybnd; return; } else { wgts(,where(*nc_var.ymax(ybnd))) = 0; } } if (!is_void(mask)) { wgtl = wgts* mask; wgto = wgts*(!mask); } var_g = (*nc_var.data * wgts)(sum,sum,,); if (!is_void(mask)) { var_l = (*nc_var.data * wgtl)(sum,sum,,); var_o = (*nc_var.data * wgto)(sum,sum,,); } if (dimsof(var_g)(2) == 1) { var_g = var_g(1,); if (!is_void(mask)) { var_l = var_l(1,); var_o = var_o(1,); } } if (!is_void(mask)) { return [var_g, var_l, var_o]; } else { return var_g; } }