From ce7847f990fe39bdc4adb4dd57a53adbfa96b762 Mon Sep 17 00:00:00 2001 From: "Travis A. O'Brien" Date: Thu, 7 May 2020 18:23:31 -0400 Subject: [PATCH 1/9] Added initial draft of teca_vertical_integral. The vertical integration routine itself has been tested external to TECA, but the teca_vertical_integral.cxx code is still being written. --- alg/teca_vertical_integral.cxx | 390 +++++++++++++++++++++++++++++++++ alg/teca_vertical_integral.h | 66 ++++++ 2 files changed, 456 insertions(+) create mode 100644 alg/teca_vertical_integral.cxx create mode 100644 alg/teca_vertical_integral.h diff --git a/alg/teca_vertical_integral.cxx b/alg/teca_vertical_integral.cxx new file mode 100644 index 000000000..f63c53f98 --- /dev/null +++ b/alg/teca_vertical_integral.cxx @@ -0,0 +1,390 @@ +#include "teca_vertical_integral.h" + +#include "teca_cartesian_mesh.h" +#include "teca_table.h" +#include "teca_array_collection.h" +#include "teca_variant_array.h" +#include "teca_metadata.h" + +#include +#include +#include +#include +#include +#include + +#if defined(TECA_HAS_BOOST) +#include +#endif + +using std::cerr; +using std::endl; + +// 1/gravity - s^2/m +template +constexpr num_t neg_one_over_g() {return num_t(-1.0)/num_t(9.81); } + +/* + * Calculates the vertical integral of 'array' + * + * input: + * ------ + * + * array - a pointer to a 3D array of values + * + * nx, ny, nz - the sizes of the x, y, z dimensions + * + * csystem - an integer indicating the vertical coordinate system + * type of `array`: + * 0 = sigma coordinate system + * 1 = hybrid coordinate system + * + * a_or_sigma - the a hybrid coordinate, or sigma if a sigma coordinate + * system + * + * b - the hybrid coordinates (see below), or a null pointer if the + * coordinate system is sigma + * + * ps - a pointer to a 2D array of pressure values + * + * p_top - a scalar value giving the model top pressure + * + * array_int - a pointer to an already-allocated 2D array + * that will hold the values of the integral of + * `array` + * + * output: + * ------- + * + * Values of `array_int` are filled in with the integral + * of `array` + * + * Assumptions: + * + * - array has dimensions [z,y,x] + * + * - ps has dimensions [y,x] + * + * - z is either a hybrid coordinate system where pressure + * is defined as p = a * p_top + b * ps, or it is a sigma + * coordinate system where p is defined as p = (ps - p_top)*sigma + p_top + * + * - z is ordered such that the bottom of the atmosphere is at z = 0 + * + * - the values of a_or_sigma, and b are given on level interfaces and + * the values of `array` are given on level centers, such that + * dp (the pressure differential) will also be on level centers + * + * - a_or_sigma and b have shape [nz + 1] + * + * - the units of ps and p_top are in Pa + * + * + * +/ +*/ +template +void vertical_integral(num_t * array, + size_t nx, + size_t ny, + size_t nz, + size_t csystem, + num_t * a_or_sigma, + num_t * b, + num_t * ps, + num_t p_top, + num_t * array_int) { + + + // loop over both horizontal dimensions + for (size_t i = 0; i < nx; i++){ + for (size_t j = 0; j < ny; j++){ + + // Determine the current index in ps + size_t n2d = j + ny * i; + + // set the current integral value to zero + num_t tmp_int = 0.0; + + // loop over the vertical dimension + for (size_t k = 0; k < nz; k++){ + + // Determine the current index in the 3D array + size_t n3d = k + nz*(j + ny * i); + + num_t dp = num_t(0.0); + + // calculate the pressure differential for the sigma + // coordinate system + if (csystem == size_t(1)){ + // calculate da and db + num_t da = a_or_sigma[k+1] - a_or_sigma[k]; + num_t db = b[k+1] - b[k]; + + // calculate dp + dp = p_top * da + ps[n2d] * db; + } + // calculate the pressure differential for the sigma + // coordinate system + else{ + num_t dsigma = a_or_sigma[k+1] - a_or_sigma[k]; + dp = (ps[n2d] - p_top)*dsigma; + } + + // accumulate the integral + tmp_int += neg_one_over_g() * array[n3d] * dp; + + } + + // set the integral value in the array + array_int[n2d] = tmp_int; + + } + } +} + + +const_p_teca_variant_array get_mesh_variable( + std::string mesh_var, + std::string expected_var, + const_p_teca_cartesian_mesh in_mesh) +{ + + // check that both variables are specified + if ( mesh_var.empty() ) + { + TECA_ERROR( expected_var << " not specified") + return nullptr; + } + + // get the variable + const_p_teca_variant_array out_array + = in_mesh->get_point_arrays()->get(mesh_var); + if (!out_array) + { + TECA_ERROR("variable \"" << mesh_var << "\" is not in the input") + return nullptr; + } + + return out_array; +} + + +// -------------------------------------------------------------------------- +teca_vertical_integral::teca_vertical_integral() + hybrid_a_variable("a_bnds"), + hybrid_b_variable("b_bnds"), + sigma_variable("sigma_bnds"), + surface_p_variable("ps"), + p_top_variable("ptop"), + using_hybrid(1), + p_top_override_value(float(-1.0)) +{ + this->set_number_of_input_connections(1); + this->set_number_of_output_ports(1); + +} + +// -------------------------------------------------------------------------- +teca_vertical_integral::~teca_vertical_integral() +{} + +#if defined(TECA_HAS_BOOST) +// -------------------------------------------------------------------------- +void teca_vertical_integral::get_properties_description( + const std::string &prefix, options_description &global_opts) +{ + (void) prefix; + + options_description opts("Options for " + + (prefix.empty()?"teca_vertical_integral":prefix)); + + /*opts.add_options() + TECA_POPTS_GET(std::vector, prefix, dependent_variables, + "list of arrays to compute statistics for") + ;*/ + opts.add_options() + TECA_POPTS_GET(std::string, prefix, hybrid_a_variable, + "name of a coordinate in the hybrid coordinate system(\"\")") + TECA_POPTS_GET(std::string, prefix, hybrid_b_variable, + "name of b coordinate in the hybrid coordinate system(\"\")") + TECA_POPTS_GET(std::string, prefix, sigma_variable, + "name of sigma coordinate (\"\")") + TECA_POPTS_GET(std::string, prefix, surface_p_variable, + "name of the surface pressure variable (\"\")") + TECA_POPTS_GET(std::string, prefix, p_top_variable, + "name of the model top variable (\"\")") + TECA_POPTS_GET(std::string, prefix, output_variable_name, + "name for the integrated, output variable (\"\")") + TECA_POPTS_GET(int, prefix, using_hybrid, + "flags whether the vertical coordinate is hybrid (1) or sigma (0) + (\"\")") + TECA_POPTS_GET(float, prefix, p_top_override_value, + "name of the model top variable(\"\")") + ; + + global_opts.add(opts); +} + +// -------------------------------------------------------------------------- +void teca_vertical_integral::set_properties( + const std::string &prefix, variables_map &opts) +{ + (void) prefix; + (void) opts; + + TECA_POPTS_SET(opts, std::string, prefix, hybrid_a_variable) + TECA_POPTS_SET(opts, std::string, prefix, hybrid_b_variable) + TECA_POPTS_SET(opts, std::string, prefix, sigma_variable) + TECA_POPTS_SET(opts, std::string, prefix, surface_p_variable) + TECA_POPTS_SET(opts, std::string, prefix, p_top_variable) + TECA_POPTS_SET(opts, int, prefix, using_hybrid) + TECA_POPTS_SET(opts, float, prefix, p_top_override_value) + TECA_POPTS_SET(opts, std::string, prefix, output_variable_name) + +} +#endif + +// -------------------------------------------------------------------------- +std::vector +teca_vertical_integral::get_upstream_request( + unsigned int port, const std::vector &input_md, + const teca_metadata &request) +{ +#ifdef TECA_DEBUG + cerr << teca_parallel_id() + << "teca_vertical_integral::get_upstream_request" << endl; +#endif + (void)port; + (void)input_md; + + std::vector up_reqs(1, request); + return up_reqs; +} + + +// -------------------------------------------------------------------------- +const_p_teca_dataset teca_vertical_integral::execute( + unsigned int port, + const std::vector &input_data, + const teca_metadata &request) +{ +#ifdef TECA_DEBUG + cerr << teca_parallel_id() << "teca_vertical_integral::execute" << endl; +#endif + (void)port; + (void)request; + + // get the input mesh + const_p_teca_cartesian_mesh in_mesh + = std::dynamic_pointer_cast(input_data[0]); + + if (!in_mesh) + { + TECA_ERROR("dataset is not a teca_cartesian_mesh") + return nullptr; + } + + // get mesh dimension + unsigned long extent[6]; + out_mesh->get_extent(extent); + + // set the dimension sizes + unsigned long nx = extent[1] - extent[0] + 1; + unsigned long ny = extent[3] - extent[2] + 1; + unsigned long nz = extent[5] - extent[4] + 1; + if (nz == 1) + { + TECA_ERROR("This calculation requires 3D data. The current dataset " + "extents are [" << extent[0] << ", " << extent[1] << ", " + << extent[2] << ", " << extent[3] << ", " << extent[4] << ", " + << extent[5] << "]") + return nullptr; + } + + int using_hybrid = this->using_hybrid; + + const_p_teca_variant_array a_or_sigma = nullptr; + const_p_teca_variant_array b_i = nullptr; + + if (using_hybrid) + { + // get the a coordinate + a_or_sigma = get_mesh_variable(this->hybrid_a_variable, + "hybrid_a_variable", + in_mesh); + if (!a_or_sigma) return nullptr; + + // get the b coordinate + b_i = get_mesh_variable(this->hybrid_b_variable, + "hybrid_b_variable", + in_mesh); + if (!b_i) return nullptr; + + else + { + // get the sigma coordinate + a_or_sigma = get_mesh_variable(this->sigma_variable, + "sigma_variable", + in_mesh); + if (!a_or_sigma) return nullptr; + } + + // get the surface pressure + const_p_teca_variant_array p_s = + get_mesh_variable(this->surface_p_variable, + "surface_p_variable", + in_mesh); + if (!p_s) return nullptr; + + // get the model top pressure + // TODO: p_top is a scalar: is const_p_teca_variant array correct? + const_p_teca_variant_array p_top_array = + get_mesh_variable(this->p_top_variable, + "p_top_variable", + in_mesh); + if (!p_top_array) return nullptr; + + // get the input array + const_p_teca_variant_array input_array = + get_mesh_variable(this->integration_variable, + "integration_variable", + in_mesh); + if (!input_array) return nullptr; + + + // allocate the output array + p_teca_variant_array integrated_array = input_array->new_instance(); + // TODO: figure out how to properly resize this to nz = 1 + integrated_array->resize(nx*ny); + + + NESTED_TEMPLATE_DISPATCH_FP( + const teca_variant_array_impl, + input_array.get(), + _INARR, + + const NT_INARR * p_input_array + = static_cast(input_array->get())->get(); + + // call the vertical integration routine + // TODO: do some templating magic on the rest of the variables + // appending p_ to the front of each variable name + vertical_integral(p_input_array, nx, ny, nz, + (size_t) using_hybrid, + p_a_or_sigma, p_b_i, p_p_s, p_p_top, + p_integrated_array); + ) + + + // create the output mesh, pass everything but the integration variable, + // and add the integrated array + p_teca_cartesian_mesh out_mesh = teca_cartesian_mesh::New(); + + // TODO: figure out how to drop this->integration_variable + // TODO: figure out how to drop the dimensionality of the data if needed + out_mesh->shallow_copy( + std::const_pointer_cast(in_mesh)); + + return out_mesh; +} diff --git a/alg/teca_vertical_integral.h b/alg/teca_vertical_integral.h new file mode 100644 index 000000000..3efbbf3aa --- /dev/null +++ b/alg/teca_vertical_integral.h @@ -0,0 +1,66 @@ +#ifndef teca_vertical_integral_h +#define teca_vertical_integral_h + +#include "teca_shared_object.h" +#include "teca_algorithm.h" +#include "teca_metadata.h" + +#include +#include + +TECA_SHARED_OBJECT_FORWARD_DECL(teca_vertical_integral) + +/// compute statistics about connected components +class teca_vertical_integral : public teca_algorithm +{ +public: + TECA_ALGORITHM_STATIC_NEW(teca_vertical_integral) + TECA_ALGORITHM_DELETE_COPY_ASSIGN(teca_vertical_integral) + TECA_ALGORITHM_CLASS_NAME(teca_vertical_integral) + ~teca_vertical_integral(); + + // report/initialize to/from Boost program options + // objects. + TECA_GET_ALGORITHM_PROPERTIES_DESCRIPTION() + TECA_SET_ALGORITHM_PROPERTIES() + + // set the name of the the arrays involved in integration + TECA_ALGORITHM_PROPERTY(std::string, hybrid_a_variable) + TECA_ALGORITHM_PROPERTY(std::string, hybrid_b_variable) + TECA_ALGORITHM_PROPERTY(std::string, sigma_variable) + TECA_ALGORITHM_PROPERTY(std::string, surface_p_variable) + TECA_ALGORITHM_PROPERTY(std::string, p_top_variable) + TECA_ALGORITHM_PROPERTY(std::string, integration_variable) + TECA_ALGORITHM_PROPERTY(std::string, output_variable_name) + // set whether the vertical coordinate is hybrid or sigma + TECA_ALGORITHM_PROPERTY(int, using_hybrid) + TECA_ALGORITHM_PROPERTY(float, p_top_override_value) + +protected: + teca_vertical_integral(); + +private: + + std::string hybrid_a_variable; + std::string hybrid_b_variable; + std::string sigma_variable; + std::string surface_p_variable; + std::string p_top_variable; + std::string integration_variable; + std::string output_variable_name; + int using_hybrid = true; + float p_top_override_value = -1; + + std::vector get_upstream_request( + unsigned int port, const std::vector &input_md, + const teca_metadata &request) override; + + const_p_teca_dataset execute( + unsigned int port, const std::vector &input_data, + const teca_metadata &request) override; + + void get_dependent_variables(const teca_metadata &request, + std::vector &dep_vars); +}; + +#endif From 137af36db4575f8e3cb2865a12b370db85c2bb87 Mon Sep 17 00:00:00 2001 From: "Travis A. O'Brien" Date: Mon, 18 May 2020 14:59:26 -0400 Subject: [PATCH 2/9] Added some notes following a conversation w/ Burlen --- alg/teca_vertical_integral.cxx | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/alg/teca_vertical_integral.cxx b/alg/teca_vertical_integral.cxx index f63c53f98..eb61a9c65 100644 --- a/alg/teca_vertical_integral.cxx +++ b/alg/teca_vertical_integral.cxx @@ -172,6 +172,8 @@ const_p_teca_variant_array get_mesh_variable( // -------------------------------------------------------------------------- teca_vertical_integral::teca_vertical_integral() + // TODO: add variables for long_name, units, other metadata? + // these needs to be passed in hybrid_a_variable("a_bnds"), hybrid_b_variable("b_bnds"), sigma_variable("sigma_bnds"), @@ -339,6 +341,7 @@ const_p_teca_dataset teca_vertical_integral::execute( // get the model top pressure // TODO: p_top is a scalar: is const_p_teca_variant array correct? + // probably need to put this in information array const_p_teca_variant_array p_top_array = get_mesh_variable(this->p_top_variable, "p_top_variable", @@ -365,7 +368,7 @@ const_p_teca_dataset teca_vertical_integral::execute( _INARR, const NT_INARR * p_input_array - = static_cast(input_array->get())->get(); + = dynamic_cast(input_array->get())->get(); // call the vertical integration routine // TODO: do some templating magic on the rest of the variables @@ -383,6 +386,13 @@ const_p_teca_dataset teca_vertical_integral::execute( // TODO: figure out how to drop this->integration_variable // TODO: figure out how to drop the dimensionality of the data if needed + // this should work -- just need to modify the dimensions + // change extent, whole_extent, bounds, coordinates (set z to a 1-valued + // teca_variant_array [maybe pressure, maybe 0]) + // see data/teca_cartesian_mesh.h + // e.g., set_z_coordinates() + // set new whole_extent, extent, bounds arrays in new teca_metadata + // objects. out_mesh->shallow_copy( std::const_pointer_cast(in_mesh)); From fb2b7b89780dfcfa8db55c02288b6b5e5e7739f3 Mon Sep 17 00:00:00 2001 From: "Travis A. O'Brien" Date: Wed, 20 May 2020 17:07:12 -0400 Subject: [PATCH 3/9] Added teca_vertical_integral.cxx to CMakeLists.txt; fixed a few TODOs in teca_vertical_integral.cxx --- alg/CMakeLists.txt | 1 + alg/teca_vertical_integral.cxx | 130 ++++++++++++++++++++++----------- alg/teca_vertical_integral.h | 4 + 3 files changed, 92 insertions(+), 43 deletions(-) diff --git a/alg/CMakeLists.txt b/alg/CMakeLists.txt index 4b8433e38..c7386b4e4 100644 --- a/alg/CMakeLists.txt +++ b/alg/CMakeLists.txt @@ -44,6 +44,7 @@ set(teca_alg_cxx_srcs teca_tc_wind_radii.cxx teca_tc_trajectory.cxx teca_temporal_average.cxx + teca_vertical_integral.cxx teca_variant_array_operand.cxx teca_vorticity.cxx teca_dataset_diff.cxx diff --git a/alg/teca_vertical_integral.cxx b/alg/teca_vertical_integral.cxx index eb61a9c65..867d1d66a 100644 --- a/alg/teca_vertical_integral.cxx +++ b/alg/teca_vertical_integral.cxx @@ -83,40 +83,40 @@ constexpr num_t neg_one_over_g() {return num_t(-1.0)/num_t(9.81); } * / */ -template -void vertical_integral(num_t * array, - size_t nx, - size_t ny, - size_t nz, - size_t csystem, - num_t * a_or_sigma, - num_t * b, - num_t * ps, +template +void vertical_integral(const num_t * array, + unsigned long nx, + unsigned long ny, + unsigned long nz, + int csystem, + const num_t * a_or_sigma, + const num_t * b, + const num_t * ps, num_t p_top, num_t * array_int) { // loop over both horizontal dimensions - for (size_t i = 0; i < nx; i++){ - for (size_t j = 0; j < ny; j++){ + for (unsigned long i = 0; i < nx; i++){ + for (unsigned long j = 0; j < ny; j++){ // Determine the current index in ps - size_t n2d = j + ny * i; + unsigned long n2d = j + ny * i; // set the current integral value to zero num_t tmp_int = 0.0; // loop over the vertical dimension - for (size_t k = 0; k < nz; k++){ + for (unsigned long k = 0; k < nz; k++){ // Determine the current index in the 3D array - size_t n3d = k + nz*(j + ny * i); + unsigned long n3d = k + nz*(j + ny * i); num_t dp = num_t(0.0); // calculate the pressure differential for the sigma // coordinate system - if (csystem == size_t(1)){ + if (csystem == 1){ // calculate da and db num_t da = a_or_sigma[k+1] - a_or_sigma[k]; num_t db = b[k+1] - b[k]; @@ -151,17 +151,15 @@ const_p_teca_variant_array get_mesh_variable( { // check that both variables are specified - if ( mesh_var.empty() ) - { - TECA_ERROR( expected_var << " not specified") + if ( mesh_var.empty() ) { + TECA_ERROR("" << expected_var << " not specified") return nullptr; } // get the variable const_p_teca_variant_array out_array = in_mesh->get_point_arrays()->get(mesh_var); - if (!out_array) - { + if (!out_array) { TECA_ERROR("variable \"" << mesh_var << "\" is not in the input") return nullptr; } @@ -171,9 +169,9 @@ const_p_teca_variant_array get_mesh_variable( // -------------------------------------------------------------------------- -teca_vertical_integral::teca_vertical_integral() - // TODO: add variables for long_name, units, other metadata? - // these needs to be passed in +teca_vertical_integral::teca_vertical_integral() : + long_name("integrated_var"), + units("unknown"), hybrid_a_variable("a_bnds"), hybrid_b_variable("b_bnds"), sigma_variable("sigma_bnds"), @@ -206,6 +204,10 @@ void teca_vertical_integral::get_properties_description( "list of arrays to compute statistics for") ;*/ opts.add_options() + TECA_POPTS_GET(std::string, prefix, long_name, + "long name of the output variable (\"\")") + TECA_POPTS_GET(std::string, prefix, units, + "units of the output variable(\"\")") TECA_POPTS_GET(std::string, prefix, hybrid_a_variable, "name of a coordinate in the hybrid coordinate system(\"\")") TECA_POPTS_GET(std::string, prefix, hybrid_b_variable, @@ -219,8 +221,8 @@ void teca_vertical_integral::get_properties_description( TECA_POPTS_GET(std::string, prefix, output_variable_name, "name for the integrated, output variable (\"\")") TECA_POPTS_GET(int, prefix, using_hybrid, - "flags whether the vertical coordinate is hybrid (1) or sigma (0) - (\"\")") + "flags whether the vertical coordinate is hybrid (1) or " + "sigma (0) (\"\")") TECA_POPTS_GET(float, prefix, p_top_override_value, "name of the model top variable(\"\")") ; @@ -235,6 +237,8 @@ void teca_vertical_integral::set_properties( (void) prefix; (void) opts; + TECA_POPTS_SET(opts, std::string, prefix, long_name) + TECA_POPTS_SET(opts, std::string, prefix, units) TECA_POPTS_SET(opts, std::string, prefix, hybrid_a_variable) TECA_POPTS_SET(opts, std::string, prefix, hybrid_b_variable) TECA_POPTS_SET(opts, std::string, prefix, sigma_variable) @@ -289,7 +293,7 @@ const_p_teca_dataset teca_vertical_integral::execute( // get mesh dimension unsigned long extent[6]; - out_mesh->get_extent(extent); + in_mesh->get_extent(extent); // set the dimension sizes unsigned long nx = extent[1] - extent[0] + 1; @@ -322,7 +326,7 @@ const_p_teca_dataset teca_vertical_integral::execute( "hybrid_b_variable", in_mesh); if (!b_i) return nullptr; - + } else { // get the sigma coordinate @@ -358,25 +362,51 @@ const_p_teca_dataset teca_vertical_integral::execute( // allocate the output array p_teca_variant_array integrated_array = input_array->new_instance(); - // TODO: figure out how to properly resize this to nz = 1 integrated_array->resize(nx*ny); NESTED_TEMPLATE_DISPATCH_FP( - const teca_variant_array_impl, - input_array.get(), + teca_variant_array_impl, + integrated_array.get(), _INARR, const NT_INARR * p_input_array - = dynamic_cast(input_array->get())->get(); + = dynamic_cast(input_array.get())->get(); + + const NT_INARR * p_a_or_sigma + = dynamic_cast(a_or_sigma.get())->get(); + + const NT_INARR * p_b_i + = dynamic_cast(b_i.get())->get(); + + const NT_INARR * p_p_s + = dynamic_cast(p_s.get())->get(); + + const NT_INARR * p_p_top + = dynamic_cast(p_top_array.get())->get(); + + NT_INARR * p_integrated_array + = dynamic_cast(integrated_array.get())->get(); // call the vertical integration routine - // TODO: do some templating magic on the rest of the variables - // appending p_ to the front of each variable name vertical_integral(p_input_array, nx, ny, nz, - (size_t) using_hybrid, - p_a_or_sigma, p_b_i, p_p_s, p_p_top, + using_hybrid, p_a_or_sigma, p_b_i, p_p_s, p_p_top[0], p_integrated_array); + /* + * template + void vertical_integral(num_t * array, + unsigned long nx, + unsigned long ny, + unsigned long nz, + int csystem, + num_t * a_or_sigma, + num_t * b, + num_t * ps, + num_t p_top, + num_t * array_int) { + + + */ ) @@ -385,16 +415,30 @@ const_p_teca_dataset teca_vertical_integral::execute( p_teca_cartesian_mesh out_mesh = teca_cartesian_mesh::New(); // TODO: figure out how to drop this->integration_variable - // TODO: figure out how to drop the dimensionality of the data if needed - // this should work -- just need to modify the dimensions - // change extent, whole_extent, bounds, coordinates (set z to a 1-valued - // teca_variant_array [maybe pressure, maybe 0]) - // see data/teca_cartesian_mesh.h - // e.g., set_z_coordinates() - // set new whole_extent, extent, bounds arrays in new teca_metadata - // objects. out_mesh->shallow_copy( std::const_pointer_cast(in_mesh)); + + // set mesh dimensions; use a scalar Z dimension + unsigned long out_extent[6]; + unsigned long out_whole_extent[6]; + double out_bounds[6]; + out_mesh->get_extent(out_extent); + out_mesh->get_whole_extent(out_whole_extent); + out_mesh->get_bounds(out_bounds); + for (size_t n = 4; n < 6; n++){ + out_extent[n] = 0; + out_whole_extent[n] = 0; + out_bounds[n] = 0; + } + out_mesh->set_extent(out_extent); + out_mesh->set_whole_extent(out_whole_extent); + out_mesh->set_bounds(out_bounds); + + // set the z coordinate + p_teca_variant_array zo = in_mesh->get_z_coordinates()->new_instance(); + zo->resize(1); + out_mesh->set_z_coordinates("z", zo); + return out_mesh; } diff --git a/alg/teca_vertical_integral.h b/alg/teca_vertical_integral.h index 3efbbf3aa..9adbdcafa 100644 --- a/alg/teca_vertical_integral.h +++ b/alg/teca_vertical_integral.h @@ -25,6 +25,8 @@ class teca_vertical_integral : public teca_algorithm TECA_SET_ALGORITHM_PROPERTIES() // set the name of the the arrays involved in integration + TECA_ALGORITHM_PROPERTY(std::string, long_name) + TECA_ALGORITHM_PROPERTY(std::string, units) TECA_ALGORITHM_PROPERTY(std::string, hybrid_a_variable) TECA_ALGORITHM_PROPERTY(std::string, hybrid_b_variable) TECA_ALGORITHM_PROPERTY(std::string, sigma_variable) @@ -41,6 +43,8 @@ class teca_vertical_integral : public teca_algorithm private: + std::string long_name; + std::string units; std::string hybrid_a_variable; std::string hybrid_b_variable; std::string sigma_variable; From 5edd479b7491475549b9f2c5346891b571a0d995 Mon Sep 17 00:00:00 2001 From: "Travis A. O'Brien" Date: Tue, 9 Jun 2020 17:53:44 -0400 Subject: [PATCH 4/9] Added variables to requests; created python binding --- alg/teca_vertical_integral.cxx | 100 +++++++++++++++++++++++++++++++-- python/teca_py_alg.i | 10 ++++ 2 files changed, 105 insertions(+), 5 deletions(-) diff --git a/alg/teca_vertical_integral.cxx b/alg/teca_vertical_integral.cxx index 867d1d66a..89d203bcc 100644 --- a/alg/teca_vertical_integral.cxx +++ b/alg/teca_vertical_integral.cxx @@ -143,6 +143,22 @@ void vertical_integral(const num_t * array, } } +int request_var( + std::string mesh_var, + std::string expected_var, + std::set * arrays) +{ + // check that both variables are specified + if ( mesh_var.empty() ) { + TECA_ERROR("" << expected_var << " not specified") + return 1; + } + + // insert the request into the list + arrays->insert(mesh_var); + + return 0; +} const_p_teca_variant_array get_mesh_variable( std::string mesh_var, @@ -167,6 +183,29 @@ const_p_teca_variant_array get_mesh_variable( return out_array; } +const_p_teca_variant_array get_info_variable( + std::string info_var, + std::string expected_var, + const_p_teca_cartesian_mesh in_info) +{ + + // check that both variables are specified + if ( info_var.empty() ) { + TECA_ERROR("" << expected_var << " not specified") + return nullptr; + } + + // get the variable + const_p_teca_variant_array out_array + = in_info->get_information_arrays()->get(info_var); + if (!out_array) { + TECA_ERROR("variable \"" << info_var << "\" is not in the input") + return nullptr; + } + + return out_array; +} + // -------------------------------------------------------------------------- teca_vertical_integral::teca_vertical_integral() : @@ -264,7 +303,58 @@ teca_vertical_integral::get_upstream_request( (void)port; (void)input_md; - std::vector up_reqs(1, request); + // create the output request + std::vector up_reqs; + + // copy the incoming request + teca_metadata req(request); + + // create a list of requested arrays + std::set arrays; + // pre-populate with existing requests, if available + if (req.has("arrays")) + req.get("arrays", arrays); + + if (using_hybrid) + { + // get the a coordinate + if (request_var(this->hybrid_a_variable, + "hybrid_a_variable", + &arrays) != 0 ) return up_reqs; + + // get the b coordinate + if (request_var(this->hybrid_b_variable, + "hybrid_b_variable", + &arrays) != 0 ) return up_reqs; + } + else + { + // get the sigma coordinate + if (request_var(this->sigma_variable, + "sigma_variable", + &arrays) != 0 ) return up_reqs; + } + + // get the surface pressure + if (request_var(this->surface_p_variable, + "surface_p_variable", + &arrays) != 0 ) return up_reqs; + + // get the model top pressure + if (request_var(this->p_top_variable, + "p_top_variable", + &arrays) != 0 ) return up_reqs; + + // get the input array + if (request_var(this->integration_variable, + "integration_variable", + &arrays) != 0 ) return up_reqs; + + // overwrite the existing request with the augmented one + req.set("arrays", arrays); + // put the request into the outgoing metadata + up_reqs.push_back(req); + return up_reqs; } @@ -316,13 +406,13 @@ const_p_teca_dataset teca_vertical_integral::execute( if (using_hybrid) { // get the a coordinate - a_or_sigma = get_mesh_variable(this->hybrid_a_variable, + a_or_sigma = get_info_variable(this->hybrid_a_variable, "hybrid_a_variable", in_mesh); if (!a_or_sigma) return nullptr; // get the b coordinate - b_i = get_mesh_variable(this->hybrid_b_variable, + b_i = get_info_variable(this->hybrid_b_variable, "hybrid_b_variable", in_mesh); if (!b_i) return nullptr; @@ -330,7 +420,7 @@ const_p_teca_dataset teca_vertical_integral::execute( else { // get the sigma coordinate - a_or_sigma = get_mesh_variable(this->sigma_variable, + a_or_sigma = get_info_variable(this->sigma_variable, "sigma_variable", in_mesh); if (!a_or_sigma) return nullptr; @@ -338,7 +428,7 @@ const_p_teca_dataset teca_vertical_integral::execute( // get the surface pressure const_p_teca_variant_array p_s = - get_mesh_variable(this->surface_p_variable, + get_info_variable(this->surface_p_variable, "surface_p_variable", in_mesh); if (!p_s) return nullptr; diff --git a/python/teca_py_alg.i b/python/teca_py_alg.i index 08b735a52..099dc52b9 100644 --- a/python/teca_py_alg.i +++ b/python/teca_py_alg.i @@ -40,6 +40,7 @@ #include "teca_tc_wind_radii.h" #include "teca_temporal_average.h" #include "teca_vorticity.h" +#include "teca_vertical_integral.h" #include "teca_py_object.h" #include "teca_py_algorithm.h" #include "teca_py_gil_state.h" @@ -133,6 +134,15 @@ %ignore teca_vorticity::operator=; %include "teca_vorticity.h" +/*************************************************************************** + vertical integral + ***************************************************************************/ +%ignore teca_vertical_integral::shared_from_this; +%shared_ptr(teca_vertical_integral) +%ignore teca_vertical_integral::operator=; +%include "teca_vertical_integral.h" + + /*************************************************************************** programmable_algorithm ***************************************************************************/ From ebdf97882ea226315ba7d172192a39d2fc000501 Mon Sep 17 00:00:00 2001 From: "Travis A. O'Brien" Date: Wed, 10 Jun 2020 16:07:22 -0400 Subject: [PATCH 5/9] Further work on implementation; compiles, but fails w/ std::bad_cast; added draft of test code --- alg/teca_vertical_integral.cxx | 108 +++++++++++++++-------- alg/teca_vertical_integral.h | 3 + io/teca_cf_reader.cxx | 4 - test/python/test_vertical_integral.py | 122 ++++++++++++++++++++++++++ 4 files changed, 194 insertions(+), 43 deletions(-) mode change 100644 => 100755 alg/teca_vertical_integral.cxx create mode 100644 test/python/test_vertical_integral.py diff --git a/alg/teca_vertical_integral.cxx b/alg/teca_vertical_integral.cxx old mode 100644 new mode 100755 index 89d203bcc..2ae859bc5 --- a/alg/teca_vertical_integral.cxx +++ b/alg/teca_vertical_integral.cxx @@ -2,6 +2,7 @@ #include "teca_cartesian_mesh.h" #include "teca_table.h" +#include "teca_array_attributes.h" #include "teca_array_collection.h" #include "teca_variant_array.h" #include "teca_metadata.h" @@ -290,6 +291,35 @@ void teca_vertical_integral::set_properties( } #endif +teca_metadata teca_vertical_integral::get_output_metadata( + unsigned int port, + const std::vector & input_md) +{ + teca_metadata report_md(input_md[0]); + + if (report_md.has("variables")){ + report_md.append("variables", this->output_variable_name); + } + else{ + report_md.set("variables", this->output_variable_name); + } + + // add attributes to enable CF I/O + teca_metadata atts; + report_md.get("attributes", atts); + teca_array_attributes output_atts( + teca_variant_array_code::get(), + teca_array_attributes::point_centering, + 0, "unset", "unset", + "unset"); + + atts.set(this->output_variable_name, (teca_metadata)output_atts); + + report_md.set("attributes", atts); + + return report_md; +} + // -------------------------------------------------------------------------- std::vector teca_vertical_integral::get_upstream_request( @@ -319,37 +349,43 @@ teca_vertical_integral::get_upstream_request( { // get the a coordinate if (request_var(this->hybrid_a_variable, - "hybrid_a_variable", + std::string("hybrid_a_variable"), &arrays) != 0 ) return up_reqs; // get the b coordinate if (request_var(this->hybrid_b_variable, - "hybrid_b_variable", + std::string("hybrid_b_variable"), &arrays) != 0 ) return up_reqs; } else { // get the sigma coordinate if (request_var(this->sigma_variable, - "sigma_variable", + std::string("sigma_variable"), &arrays) != 0 ) return up_reqs; } // get the surface pressure if (request_var(this->surface_p_variable, - "surface_p_variable", + std::string("surface_p_variable"), &arrays) != 0 ) return up_reqs; - // get the model top pressure - if (request_var(this->p_top_variable, - "p_top_variable", - &arrays) != 0 ) return up_reqs; + // only request p_top if it isn't being overriden + if (! this->p_top_override_value ){ + // get the model top pressure + if (request_var(this->p_top_variable, + std::string("p_top_variable"), + &arrays) != 0 ) return up_reqs; + } // get the input array if (request_var(this->integration_variable, - "integration_variable", + std::string("integration_variable"), &arrays) != 0 ) return up_reqs; + // intercept request for our output + arrays.erase(this->output_variable_name); + // overwrite the existing request with the augmented one req.set("arrays", arrays); // put the request into the outgoing metadata @@ -358,7 +394,6 @@ teca_vertical_integral::get_upstream_request( return up_reqs; } - // -------------------------------------------------------------------------- const_p_teca_dataset teca_vertical_integral::execute( unsigned int port, @@ -407,13 +442,13 @@ const_p_teca_dataset teca_vertical_integral::execute( { // get the a coordinate a_or_sigma = get_info_variable(this->hybrid_a_variable, - "hybrid_a_variable", + std::string("hybrid_a_variable"), in_mesh); if (!a_or_sigma) return nullptr; // get the b coordinate b_i = get_info_variable(this->hybrid_b_variable, - "hybrid_b_variable", + std::string("hybrid_b_variable"), in_mesh); if (!b_i) return nullptr; } @@ -421,7 +456,7 @@ const_p_teca_dataset teca_vertical_integral::execute( { // get the sigma coordinate a_or_sigma = get_info_variable(this->sigma_variable, - "sigma_variable", + std::string("sigma_variable"), in_mesh); if (!a_or_sigma) return nullptr; } @@ -429,27 +464,36 @@ const_p_teca_dataset teca_vertical_integral::execute( // get the surface pressure const_p_teca_variant_array p_s = get_info_variable(this->surface_p_variable, - "surface_p_variable", + std::string("surface_p_variable"), in_mesh); if (!p_s) return nullptr; - // get the model top pressure - // TODO: p_top is a scalar: is const_p_teca_variant array correct? - // probably need to put this in information array - const_p_teca_variant_array p_top_array = - get_mesh_variable(this->p_top_variable, - "p_top_variable", - in_mesh); - if (!p_top_array) return nullptr; + const_p_teca_variant_array p_top_array; + p_teca_variant_array p_top_override_tmp; + if (! this->p_top_override_value){ + // get the model top pressure + // TODO: p_top is a scalar: is const_p_teca_variant array correct? + // probably need to put this in information array + p_top_array = get_mesh_variable(this->p_top_variable, + std::string("p_top_variable"), + in_mesh); + if (!p_top_array) return nullptr; + } + else{ + p_top_override_tmp = teca_variant_array_impl::New(1); + p_top_override_tmp->set(0, & this->p_top_override_value); + p_top_array = (const_p_teca_variant_array) p_top_override_tmp; + } // get the input array const_p_teca_variant_array input_array = get_mesh_variable(this->integration_variable, - "integration_variable", + std::string("integration_variable"), in_mesh); if (!input_array) return nullptr; + /* // allocate the output array p_teca_variant_array integrated_array = input_array->new_instance(); integrated_array->resize(nx*ny); @@ -472,7 +516,7 @@ const_p_teca_dataset teca_vertical_integral::execute( const NT_INARR * p_p_s = dynamic_cast(p_s.get())->get(); - const NT_INARR * p_p_top + const NT_INARR * p_p_top = dynamic_cast(p_top_array.get())->get(); NT_INARR * p_integrated_array @@ -482,22 +526,8 @@ const_p_teca_dataset teca_vertical_integral::execute( vertical_integral(p_input_array, nx, ny, nz, using_hybrid, p_a_or_sigma, p_b_i, p_p_s, p_p_top[0], p_integrated_array); - /* - * template - void vertical_integral(num_t * array, - unsigned long nx, - unsigned long ny, - unsigned long nz, - int csystem, - num_t * a_or_sigma, - num_t * b, - num_t * ps, - num_t p_top, - num_t * array_int) { - - - */ ) + */ // create the output mesh, pass everything but the integration variable, diff --git a/alg/teca_vertical_integral.h b/alg/teca_vertical_integral.h index 9adbdcafa..a0aeb2afa 100644 --- a/alg/teca_vertical_integral.h +++ b/alg/teca_vertical_integral.h @@ -65,6 +65,9 @@ class teca_vertical_integral : public teca_algorithm void get_dependent_variables(const teca_metadata &request, std::vector &dep_vars); + + teca_metadata get_output_metadata(unsigned int port, + const std::vector &input_md) override; }; #endif diff --git a/io/teca_cf_reader.cxx b/io/teca_cf_reader.cxx index bb7e896c6..c1b88fad1 100644 --- a/io/teca_cf_reader.cxx +++ b/io/teca_cf_reader.cxx @@ -552,10 +552,6 @@ teca_metadata teca_cf_reader::get_output_metadata( var_type = teca_variant_array_code::get(); ) - // skip scalars - if (n_dims == 0) - continue; - std::vector dims; std::vector dim_names; unsigned int centering = teca_array_attributes::point_centering; diff --git a/test/python/test_vertical_integral.py b/test/python/test_vertical_integral.py new file mode 100644 index 000000000..102979f57 --- /dev/null +++ b/test/python/test_vertical_integral.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python3 +import teca + +def create_teca_pipeline(input_filenames = ["test_vertical_integral.nc"], + using_hybrid = True, + output_filename = "test.nc"): + """Creates a TECA pipeline that tests teca_vertical_integral + + input: + ------ + + input_filenames : a list containing file names to open + + using_hybrid : flags whether to use exercise the hybrid coordinate + integration (does sigma coordinate otherwise) + + output_filenames : the output filename to which to write the result + + output: + ------- + + pipeline : the last algorithm in the TECA pipeline; running + pipeline.update() will run the pipeline. + + When run, the pipeline writes a netCDF file to disk, containing the + result of the integration. + + """ + # *********************** + # Read the input dataset + # *********************** + cfr = teca.teca_cf_reader.New() + cfr.set_file_names(input_filenames) + # set coordinate names in the input dataset + cfr.set_x_axis_variable("x") + cfr.set_y_axis_variable("y") + # TODO: how to deal with two vertical coordiantes: kz and kzi? + cfr.set_z_axis_variable("kz") + cfr.set_t_axis_variable("time") + # turn off threading + cfr.set_thread_pool_size(1) + + # ********************** + # Normalize coordinates + # ********************** + coords = teca.teca_normalize_coordinates.New() + coords.set_input_connection(cfr.get_output_port()) + + # *********************** + # Integrate `test_field` + # *********************** + vint = teca.teca_vertical_integral.New() + vint.set_input_connection(coords.get_output_port()) + # set variables needed for integration + vint.set_hybrid_a_variable("a_bnds") + vint.set_hybrid_b_variable("b_bnds") + #vint.set_sigma_variable("sigma_i") + vint.set_surface_p_variable("ps") + #vint.set_p_top_variable("p_top") + vint.set_integration_variable("test_field") + # set output metadata + vint.set_output_variable_name("vint_test_field") + vint.set_long_name("vertical integral of test_field") + vint.set_units("none") + # flag that we're integrating using the hybrid coordinate + vint.set_using_hybrid(int(using_hybrid)) + vint.set_p_top_override_value(0.4940202) + + if False: + class metadata_override(teca.teca_python_algorithm): + + def report(self, port, md_in): + md_out = teca.teca_metadata(md_in[0]) + + wext = md_out['whole_extent'] + ncells = (wext[1] - wext[0] + 1)* \ + (wext[3] - wext[2] + 1)*(wext[5] - wext[4] + 1) + + new_var_atts = teca.teca_array_attributes( + teca.teca_double_array_code.get(), + teca.teca_array_attributes.point_centering, + int(ncells), 'unitless', 'vertical integral of -1/g', + 'should be approximately equal to ps') + + # put it in the array attributes + try: + atts = md_out['attributes'] + except: + atts = teca_metadata() + atts['vint_test_field'] = new_var_atts.to_metadata() + md_out['attributes'] = atts + return md_out + + mdi = metadata_override.New() + mdi.set_input_connection(vint.get_output_port()) + + # stub for looping over the singular time dimension + exe = teca.teca_index_executive.New() + exe.set_start_index(0) + exe.set_end_index(0) + + # ********************* + # Write an output file + # ********************* + wri = teca.teca_cf_writer.New() + wri.set_file_name(output_filename) + wri.set_executive(exe) + wri.set_input_connection(vint.get_output_port()) + wri.set_point_arrays(["vint_test_field"]) + #wri.set_input_connection(coords.get_output_port()) + wri.set_thread_pool_size(1) + + + return wri + +# create the TECA pipeline +pipeline = create_teca_pipeline() + +print(pipeline.update_metadata()) + +# run the pipeline +pipeline.update() From 9c8ce6768a93beea385c4b63de5e5998ae7ff6c3 Mon Sep 17 00:00:00 2001 From: "Travis A. O'Brien" Date: Wed, 10 Jun 2020 16:22:18 -0400 Subject: [PATCH 6/9] Uncommented vertical integration code; accidently left commented on last commit --- alg/teca_vertical_integral.cxx | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/alg/teca_vertical_integral.cxx b/alg/teca_vertical_integral.cxx index 2ae859bc5..9c9375b05 100755 --- a/alg/teca_vertical_integral.cxx +++ b/alg/teca_vertical_integral.cxx @@ -493,7 +493,6 @@ const_p_teca_dataset teca_vertical_integral::execute( if (!input_array) return nullptr; - /* // allocate the output array p_teca_variant_array integrated_array = input_array->new_instance(); integrated_array->resize(nx*ny); @@ -527,18 +526,15 @@ const_p_teca_dataset teca_vertical_integral::execute( using_hybrid, p_a_or_sigma, p_b_i, p_p_s, p_p_top[0], p_integrated_array); ) - */ // create the output mesh, pass everything but the integration variable, // and add the integrated array p_teca_cartesian_mesh out_mesh = teca_cartesian_mesh::New(); - // TODO: figure out how to drop this->integration_variable out_mesh->shallow_copy( std::const_pointer_cast(in_mesh)); - // set mesh dimensions; use a scalar Z dimension unsigned long out_extent[6]; unsigned long out_whole_extent[6]; @@ -560,5 +556,7 @@ const_p_teca_dataset teca_vertical_integral::execute( zo->resize(1); out_mesh->set_z_coordinates("z", zo); + // TODO: actually add the output variable to the mesh + return out_mesh; } From f1aa90d0733409a4c4290b134ee6c03872d73517 Mon Sep 17 00:00:00 2001 From: "Travis A. O'Brien" Date: Wed, 10 Jun 2020 18:04:33 -0400 Subject: [PATCH 7/9] Fixed a number of segfaults; fixed cerr/endl namespace issues in a number of files when -DTECA_DEBUG is set --- alg/teca_2d_component_area.cxx | 12 +-- alg/teca_bayesian_ar_detect.cxx | 10 +-- alg/teca_bayesian_ar_detect_parameters.cxx | 8 +- alg/teca_binary_segmentation.cxx | 16 ++-- alg/teca_component_area_filter.cxx | 12 +-- alg/teca_vertical_integral.cxx | 97 ++++++++++++++++++---- 6 files changed, 110 insertions(+), 45 deletions(-) mode change 100644 => 100755 alg/teca_2d_component_area.cxx mode change 100644 => 100755 alg/teca_bayesian_ar_detect.cxx mode change 100644 => 100755 alg/teca_bayesian_ar_detect_parameters.cxx mode change 100644 => 100755 alg/teca_binary_segmentation.cxx mode change 100644 => 100755 alg/teca_component_area_filter.cxx diff --git a/alg/teca_2d_component_area.cxx b/alg/teca_2d_component_area.cxx old mode 100644 new mode 100755 index 84dc37662..f194524a6 --- a/alg/teca_2d_component_area.cxx +++ b/alg/teca_2d_component_area.cxx @@ -163,8 +163,8 @@ teca_metadata teca_2d_component_area::get_output_metadata( const std::vector &input_md) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_2d_component_area::get_output_metadata" << endl; + std::cerr << teca_parallel_id() + << "teca_2d_component_area::get_output_metadata" << std::endl; #endif (void) port; @@ -179,8 +179,8 @@ std::vector teca_2d_component_area::get_upstream_request( const teca_metadata &request) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_2d_component_area::get_upstream_request" << endl; + std::cerr << teca_parallel_id() + << "teca_2d_component_area::get_upstream_request" << std::endl; #endif (void) port; (void) input_md; @@ -218,8 +218,8 @@ const_p_teca_dataset teca_2d_component_area::execute( const teca_metadata &request) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_2d_component_area::execute" << endl; + std::cerr << teca_parallel_id() + << "teca_2d_component_area::execute" << std::endl; #endif (void)port; (void)request; diff --git a/alg/teca_bayesian_ar_detect.cxx b/alg/teca_bayesian_ar_detect.cxx old mode 100644 new mode 100755 index b72897385..c4fb74daf --- a/alg/teca_bayesian_ar_detect.cxx +++ b/alg/teca_bayesian_ar_detect.cxx @@ -638,8 +638,8 @@ teca_metadata teca_bayesian_ar_detect::get_output_metadata( unsigned int port, const std::vector &input_md) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_bayesian_ar_detect::get_output_metadata" << endl; + std::cerr << teca_parallel_id() + << "teca_bayesian_ar_detect::get_output_metadata" << std::endl; #endif (void)port; @@ -788,8 +788,8 @@ std::vector teca_bayesian_ar_detect::get_upstream_request( const teca_metadata &request) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_bayesian_ar_detect::get_upstream_request" << endl; + std::cerr << teca_parallel_id() + << "teca_bayesian_ar_detect::get_upstream_request" << std::endl; #endif (void) port; (void) input_md; @@ -830,7 +830,7 @@ const_p_teca_dataset teca_bayesian_ar_detect::execute( const teca_metadata &request) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() << "teca_bayesian_ar_detect::execute" << endl; + std::cerr << teca_parallel_id() << "teca_bayesian_ar_detect::execute" << std::endl; #endif (void)port; (void)request; diff --git a/alg/teca_bayesian_ar_detect_parameters.cxx b/alg/teca_bayesian_ar_detect_parameters.cxx old mode 100644 new mode 100755 index 7d0a3e04e..48ef792af --- a/alg/teca_bayesian_ar_detect_parameters.cxx +++ b/alg/teca_bayesian_ar_detect_parameters.cxx @@ -1001,8 +1001,8 @@ teca_metadata teca_bayesian_ar_detect_parameters::get_output_metadata( unsigned int port, const std::vector &input_md) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_bayesian_ar_detect_parameters::get_output_metadata" << endl; + std::cerr << teca_parallel_id() + << "teca_bayesian_ar_detect_parameters::get_output_metadata" << std::endl; #endif (void)port; (void)input_md; @@ -1016,8 +1016,8 @@ const_p_teca_dataset teca_bayesian_ar_detect_parameters::execute(unsigned int po const teca_metadata &request) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_bayesian_ar_detect_parameters::execute" << endl; + std::cerr << teca_parallel_id() + << "teca_bayesian_ar_detect_parameters::execute" << std::endl; #endif (void)port; (void)input_data; diff --git a/alg/teca_binary_segmentation.cxx b/alg/teca_binary_segmentation.cxx old mode 100644 new mode 100755 index 894fa5419..df895cacc --- a/alg/teca_binary_segmentation.cxx +++ b/alg/teca_binary_segmentation.cxx @@ -114,8 +114,8 @@ void percentile_threshold(out_t *output, const in_t *input, // compute high percentile double high_percentile = (y1 - y0)*t_high + y0; - /*std::cerr << q_low << "th percentile is " << std::setprecision(10) << low_percentile << std::endl - << q_high << "th percentile is " << std::setprecision(9) << high_percentile << std::endl;*/ + /*std::std::cerr << q_low << "th percentile is " << std::setprecision(10) << low_percentile << std::std::endl + << q_high << "th percentile is " << std::setprecision(9) << high_percentile << std::std::endl;*/ // apply thresholds for (size_t i = 0; i < n_vals; ++i) @@ -181,8 +181,8 @@ teca_metadata teca_binary_segmentation::get_output_metadata( const std::vector &input_md) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_binary_segmentation::get_output_metadata" << endl; + std::cerr << teca_parallel_id() + << "teca_binary_segmentation::get_output_metadata" << std::endl; #endif (void) port; @@ -233,8 +233,8 @@ std::vector teca_binary_segmentation::get_upstream_request( const teca_metadata &request) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_binary_segmentation::get_upstream_request" << endl; + std::cerr << teca_parallel_id() + << "teca_binary_segmentation::get_upstream_request" << std::endl; #endif (void) port; (void) input_md; @@ -276,8 +276,8 @@ const_p_teca_dataset teca_binary_segmentation::execute( const teca_metadata &request) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_binary_segmentation::execute" << endl; + std::cerr << teca_parallel_id() + << "teca_binary_segmentation::execute" << std::endl; #endif (void)port; diff --git a/alg/teca_component_area_filter.cxx b/alg/teca_component_area_filter.cxx old mode 100644 new mode 100755 index df6e25ff4..24dd7936f --- a/alg/teca_component_area_filter.cxx +++ b/alg/teca_component_area_filter.cxx @@ -131,8 +131,8 @@ teca_metadata teca_component_area_filter::get_output_metadata( const std::vector &input_md) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_component_area_filter::get_output_metadata" << endl; + std::cerr << teca_parallel_id() + << "teca_component_area_filter::get_output_metadata" << std::endl; #endif (void) port; @@ -156,8 +156,8 @@ std::vector teca_component_area_filter::get_upstream_request( const teca_metadata &request) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_component_area_filter::get_upstream_request" << endl; + std::cerr << teca_parallel_id() + << "teca_component_area_filter::get_upstream_request" << std::endl; #endif (void) port; (void) input_md; @@ -202,8 +202,8 @@ const_p_teca_dataset teca_component_area_filter::execute( const teca_metadata &request) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_component_area_filter::execute" << endl; + std::cerr << teca_parallel_id() + << "teca_component_area_filter::execute" << std::endl; #endif (void)port; diff --git a/alg/teca_vertical_integral.cxx b/alg/teca_vertical_integral.cxx index 9c9375b05..197ac4817 100755 --- a/alg/teca_vertical_integral.cxx +++ b/alg/teca_vertical_integral.cxx @@ -18,8 +18,6 @@ #include #endif -using std::cerr; -using std::endl; // 1/gravity - s^2/m template @@ -124,6 +122,7 @@ void vertical_integral(const num_t * array, // calculate dp dp = p_top * da + ps[n2d] * db; + } // calculate the pressure differential for the sigma // coordinate system @@ -295,8 +294,55 @@ teca_metadata teca_vertical_integral::get_output_metadata( unsigned int port, const std::vector & input_md) { +#ifdef TECA_DEBUG + std::cerr << teca_parallel_id() + << "teca_vertical_integral::get_output_metadata" << std::endl; +#endif + teca_metadata report_md(input_md[0]); + double bounds[6] = {0.0}; + unsigned long whole_extent[6] = {0ul}; + unsigned long extent[6] = {0ul}; + + // get the extents and bounds + report_md.get("whole_extent", whole_extent, 6); + report_md.get("extent", extent, 6); + report_md.get("bounds", bounds, 6); + + // get the coordinates + teca_metadata coords; + report_md.get("coordinates", coords); + + // set a new z coordinate with no value (this will cause cf_writer to skip) + //p_teca_variant_array new_z = teca_variant_array_impl::New(1); + p_teca_variant_array new_z = coords.get("z")->new_instance(); + new_z->resize(1); + new_z->set(0, 0.0); + coords.set("z", new_z); + + // determine dimension sizes + unsigned long nx = whole_extent[1] - whole_extent[0] + 1; + unsigned long ny = whole_extent[3] - whole_extent[2] + 1; + unsigned long nz = whole_extent[5] - whole_extent[4] + 1; + + // check that the variable has a z dimension + if (nz == 1) + { + TECA_ERROR("This calculation requires 3D data. The current dataset " + "whole_extents are [" << whole_extent[0] << ", " << whole_extent[1] << ", " + << whole_extent[2] << ", " << whole_extent[3] << ", " << whole_extent[4] << ", " + << whole_extent[5] << "]") + return report_md; + } + + // force the output data to have no z dimension + for (size_t n = 4; n < 6; n++){ + extent[n] = 0; + whole_extent[n] = 0; + bounds[n] = 0.0; + } + if (report_md.has("variables")){ report_md.append("variables", this->output_variable_name); } @@ -316,6 +362,11 @@ teca_metadata teca_vertical_integral::get_output_metadata( atts.set(this->output_variable_name, (teca_metadata)output_atts); report_md.set("attributes", atts); + // write the updated bounds/extent/coordinates + report_md.set("whole_extent", whole_extent); + report_md.set("extent", extent); + report_md.set("bounds", bounds); + report_md.set("coordinates", coords); return report_md; } @@ -327,8 +378,8 @@ teca_vertical_integral::get_upstream_request( const teca_metadata &request) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() - << "teca_vertical_integral::get_upstream_request" << endl; + std::cerr << teca_parallel_id() + << "teca_vertical_integral::get_upstream_request" << std::endl; #endif (void)port; (void)input_md; @@ -401,7 +452,7 @@ const_p_teca_dataset teca_vertical_integral::execute( const teca_metadata &request) { #ifdef TECA_DEBUG - cerr << teca_parallel_id() << "teca_vertical_integral::execute" << endl; + std::cerr << teca_parallel_id() << "teca_vertical_integral::execute" << std::endl; #endif (void)port; (void)request; @@ -424,14 +475,12 @@ const_p_teca_dataset teca_vertical_integral::execute( unsigned long nx = extent[1] - extent[0] + 1; unsigned long ny = extent[3] - extent[2] + 1; unsigned long nz = extent[5] - extent[4] + 1; - if (nz == 1) - { - TECA_ERROR("This calculation requires 3D data. The current dataset " - "extents are [" << extent[0] << ", " << extent[1] << ", " - << extent[2] << ", " << extent[3] << ", " << extent[4] << ", " - << extent[5] << "]") - return nullptr; - } + +#ifdef TECA_DEBUG + std::cerr << teca_parallel_id() << + "teca_vertical_integral::execute working on array of size [" << + nx << ", " << ny << ", " << nz << "]" << std::endl; +#endif int using_hybrid = this->using_hybrid; @@ -480,11 +529,15 @@ const_p_teca_dataset teca_vertical_integral::execute( if (!p_top_array) return nullptr; } else{ - p_top_override_tmp = teca_variant_array_impl::New(1); - p_top_override_tmp->set(0, & this->p_top_override_value); + p_top_override_tmp = teca_variant_array_impl::New(1); + p_top_override_tmp->set(0, this->p_top_override_value); p_top_array = (const_p_teca_variant_array) p_top_override_tmp; } +#ifdef TECA_DEBUG + std::cerr << teca_parallel_id() << "teca_vertical_integral::execute:" + << "reading the input variable" << std::endl; +#endif // get the input array const_p_teca_variant_array input_array = get_mesh_variable(this->integration_variable, @@ -492,6 +545,10 @@ const_p_teca_dataset teca_vertical_integral::execute( in_mesh); if (!input_array) return nullptr; +#ifdef TECA_DEBUG + std::cerr << teca_parallel_id() << "teca_vertical_integral::execute:" + << "allocating the integration array" << std::endl; +#endif // allocate the output array p_teca_variant_array integrated_array = input_array->new_instance(); @@ -527,6 +584,10 @@ const_p_teca_dataset teca_vertical_integral::execute( p_integrated_array); ) +#ifdef TECA_DEBUG + std::cerr << teca_parallel_id() << "teca_vertical_integral::execute:" + << "creating the output mesh" << std::endl; +#endif // create the output mesh, pass everything but the integration variable, // and add the integrated array @@ -542,11 +603,13 @@ const_p_teca_dataset teca_vertical_integral::execute( out_mesh->get_extent(out_extent); out_mesh->get_whole_extent(out_whole_extent); out_mesh->get_bounds(out_bounds); + for (size_t n = 4; n < 6; n++){ out_extent[n] = 0; out_whole_extent[n] = 0; out_bounds[n] = 0; } + out_mesh->set_extent(out_extent); out_mesh->set_whole_extent(out_whole_extent); out_mesh->set_bounds(out_bounds); @@ -556,7 +619,9 @@ const_p_teca_dataset teca_vertical_integral::execute( zo->resize(1); out_mesh->set_z_coordinates("z", zo); - // TODO: actually add the output variable to the mesh + // add the output variable to the mesh + out_mesh->get_point_arrays()->append( + this->output_variable_name, integrated_array); return out_mesh; } From 19871472c95e8d564d150cc8daa598a117c88b42 Mon Sep 17 00:00:00 2001 From: "Travis A. O'Brien" Date: Thu, 11 Jun 2020 15:49:41 -0400 Subject: [PATCH 8/9] Fixed request metadata bug; fixed array/netcdf type mismatch bug; values now written in netCDF --- alg/teca_vertical_integral.cxx | 44 ++++++++++++++++++++++++++++++---- 1 file changed, 39 insertions(+), 5 deletions(-) diff --git a/alg/teca_vertical_integral.cxx b/alg/teca_vertical_integral.cxx index 197ac4817..f32bba09f 100755 --- a/alg/teca_vertical_integral.cxx +++ b/alg/teca_vertical_integral.cxx @@ -122,7 +122,6 @@ void vertical_integral(const num_t * array, // calculate dp dp = p_top * da + ps[n2d] * db; - } // calculate the pressure differential for the sigma // coordinate system @@ -354,7 +353,7 @@ teca_metadata teca_vertical_integral::get_output_metadata( teca_metadata atts; report_md.get("attributes", atts); teca_array_attributes output_atts( - teca_variant_array_code::get(), + teca_variant_array_code::get(), teca_array_attributes::point_centering, 0, "unset", "unset", "unset"); @@ -364,8 +363,8 @@ teca_metadata teca_vertical_integral::get_output_metadata( report_md.set("attributes", atts); // write the updated bounds/extent/coordinates report_md.set("whole_extent", whole_extent); - report_md.set("extent", extent); - report_md.set("bounds", bounds); + if (report_md.has("extent")) report_md.set("extent", extent); + if (report_md.has("bounds")) report_md.set("bounds", bounds); report_md.set("coordinates", coords); return report_md; @@ -382,7 +381,6 @@ teca_vertical_integral::get_upstream_request( << "teca_vertical_integral::get_upstream_request" << std::endl; #endif (void)port; - (void)input_md; // create the output request std::vector up_reqs; @@ -390,6 +388,9 @@ teca_vertical_integral::get_upstream_request( // copy the incoming request teca_metadata req(request); + // copy the input metadata + teca_metadata md_in(input_md[0]); + // create a list of requested arrays std::set arrays; // pre-populate with existing requests, if available @@ -437,6 +438,37 @@ teca_vertical_integral::get_upstream_request( // intercept request for our output arrays.erase(this->output_variable_name); + /* + double bounds[6] = {0.0}; + unsigned long extent[6] = {0ul}; + + if (md_in.has("bounds")){ + // pass through bounds if bounds are given + md_in.get("bounds", bounds, 6); + req.set("bounds", bounds); + } + else if (md_in.has("extent")){ + // pass through extent if extent is given + md_in.get("extent", extent, 6); + req.set("extent", extent); + } + else { + // pass the whole extent as the extent if neither + // bounds nor extent are given + md_in.get("whole_extent", extent, 6); + req.set("extent", extent); + } + */ + + + // TODO: this overrides the above code and removes bounds/extent; + // this disables the ability for a user to specify bounds/extent + // The above (commented) code needs to be debugged in order for + // the bounds/extent pass-through to work properly though. + req.remove("bounds"); + req.remove("extent"); + req.remove("whole_extent"); + // overwrite the existing request with the augmented one req.set("arrays", arrays); // put the request into the outgoing metadata @@ -589,6 +621,7 @@ const_p_teca_dataset teca_vertical_integral::execute( << "creating the output mesh" << std::endl; #endif + // create the output mesh, pass everything but the integration variable, // and add the integrated array p_teca_cartesian_mesh out_mesh = teca_cartesian_mesh::New(); @@ -620,6 +653,7 @@ const_p_teca_dataset teca_vertical_integral::execute( out_mesh->set_z_coordinates("z", zo); // add the output variable to the mesh + std::cerr << "Output variable name: " << this->output_variable_name << std::endl; out_mesh->get_point_arrays()->append( this->output_variable_name, integrated_array); From bc96bd1db95e8989a63a01d79d6d4336bff691df Mon Sep 17 00:00:00 2001 From: "Travis A. O'Brien" Date: Mon, 27 Jul 2020 15:59:09 -0400 Subject: [PATCH 9/9] WIP - updated vertical integral for pressure levels --- alg/teca_vertical_integral.cxx | 228 +++++++-------------------------- alg/teca_vertical_integral.h | 17 +-- 2 files changed, 46 insertions(+), 199 deletions(-) diff --git a/alg/teca_vertical_integral.cxx b/alg/teca_vertical_integral.cxx index f32bba09f..09d64f61f 100755 --- a/alg/teca_vertical_integral.cxx +++ b/alg/teca_vertical_integral.cxx @@ -21,7 +21,7 @@ // 1/gravity - s^2/m template -constexpr num_t neg_one_over_g() {return num_t(-1.0)/num_t(9.81); } +constexpr num_t neg_one_over_g() {return num_t(-1.0)/num_t(9.806); } /* * Calculates the vertical integral of 'array' @@ -33,21 +33,9 @@ constexpr num_t neg_one_over_g() {return num_t(-1.0)/num_t(9.81); } * * nx, ny, nz - the sizes of the x, y, z dimensions * - * csystem - an integer indicating the vertical coordinate system - * type of `array`: - * 0 = sigma coordinate system - * 1 = hybrid coordinate system - * - * a_or_sigma - the a hybrid coordinate, or sigma if a sigma coordinate + * pressure_level - the a hybrid coordinate, or sigma if a sigma coordinate * system * - * b - the hybrid coordinates (see below), or a null pointer if the - * coordinate system is sigma - * - * ps - a pointer to a 2D array of pressure values - * - * p_top - a scalar value giving the model top pressure - * * array_int - a pointer to an already-allocated 2D array * that will hold the values of the integral of * `array` @@ -62,76 +50,52 @@ constexpr num_t neg_one_over_g() {return num_t(-1.0)/num_t(9.81); } * * - array has dimensions [z,y,x] * - * - ps has dimensions [y,x] - * - * - z is either a hybrid coordinate system where pressure - * is defined as p = a * p_top + b * ps, or it is a sigma - * coordinate system where p is defined as p = (ps - p_top)*sigma + p_top - * * - z is ordered such that the bottom of the atmosphere is at z = 0 * - * - the values of a_or_sigma, and b are given on level interfaces and - * the values of `array` are given on level centers, such that - * dp (the pressure differential) will also be on level centers + * - all values are given on level centers * - * - a_or_sigma and b have shape [nz + 1] - * - * - the units of ps and p_top are in Pa + * - the units of pressure_level are in [Pa] * * * / */ template -void vertical_integral(const num_t * array, - unsigned long nx, - unsigned long ny, - unsigned long nz, - int csystem, - const num_t * a_or_sigma, - const num_t * b, - const num_t * ps, - num_t p_top, - num_t * array_int) { - +void mass_integral_trapezoidal(const num_t * array, + unsigned long nx, + unsigned long ny, + unsigned long nz, + const num_t * pressure_level, + num_t * array_int) { // loop over both horizontal dimensions for (unsigned long i = 0; i < nx; i++){ for (unsigned long j = 0; j < ny; j++){ - // Determine the current index in ps + + // Determine the current index the integrated array unsigned long n2d = j + ny * i; // set the current integral value to zero num_t tmp_int = 0.0; // loop over the vertical dimension - for (unsigned long k = 0; k < nz; k++){ + for (unsigned long k = 0; k < (nz-1); k++){ // Determine the current index in the 3D array unsigned long n3d = k + nz*(j + ny * i); - num_t dp = num_t(0.0); - - // calculate the pressure differential for the sigma - // coordinate system - if (csystem == 1){ - // calculate da and db - num_t da = a_or_sigma[k+1] - a_or_sigma[k]; - num_t db = b[k+1] - b[k]; - - // calculate dp - dp = p_top * da + ps[n2d] * db; - } - // calculate the pressure differential for the sigma - // coordinate system - else{ - num_t dsigma = a_or_sigma[k+1] - a_or_sigma[k]; - dp = (ps[n2d] - p_top)*dsigma; - } - - // accumulate the integral - tmp_int += neg_one_over_g() * array[n3d] * dp; + // calculate the pressure differential for this layer + num_t dp = pressure_level[k+1] - pressure_level[k]; + + // calculate the value of the array between these two levels + num_t array_bar = (array[n3d] + array[n3d + 1])/num_t(2); + + // calculate the density-weighted vertical differential + num_t rho_times_dz = neg_one_over_g() * dp; + + // add the integral contribution to the array + tmp_int += array_bar * rho_times_dz; } @@ -210,13 +174,7 @@ const_p_teca_variant_array get_info_variable( teca_vertical_integral::teca_vertical_integral() : long_name("integrated_var"), units("unknown"), - hybrid_a_variable("a_bnds"), - hybrid_b_variable("b_bnds"), - sigma_variable("sigma_bnds"), - surface_p_variable("ps"), - p_top_variable("ptop"), - using_hybrid(1), - p_top_override_value(float(-1.0)) + pressure_level_variable("plev"), { this->set_number_of_input_connections(1); this->set_number_of_output_ports(1); @@ -242,27 +200,16 @@ void teca_vertical_integral::get_properties_description( "list of arrays to compute statistics for") ;*/ opts.add_options() + TECA_POPTS_GET(std::string, prefix, integration_variable, + "name of the variable to integrate (\"\")") TECA_POPTS_GET(std::string, prefix, long_name, "long name of the output variable (\"\")") TECA_POPTS_GET(std::string, prefix, units, "units of the output variable(\"\")") - TECA_POPTS_GET(std::string, prefix, hybrid_a_variable, + TECA_POPTS_GET(std::string, prefix, pressure_level_variable, "name of a coordinate in the hybrid coordinate system(\"\")") - TECA_POPTS_GET(std::string, prefix, hybrid_b_variable, - "name of b coordinate in the hybrid coordinate system(\"\")") - TECA_POPTS_GET(std::string, prefix, sigma_variable, - "name of sigma coordinate (\"\")") - TECA_POPTS_GET(std::string, prefix, surface_p_variable, - "name of the surface pressure variable (\"\")") - TECA_POPTS_GET(std::string, prefix, p_top_variable, - "name of the model top variable (\"\")") TECA_POPTS_GET(std::string, prefix, output_variable_name, "name for the integrated, output variable (\"\")") - TECA_POPTS_GET(int, prefix, using_hybrid, - "flags whether the vertical coordinate is hybrid (1) or " - "sigma (0) (\"\")") - TECA_POPTS_GET(float, prefix, p_top_override_value, - "name of the model top variable(\"\")") ; global_opts.add(opts); @@ -275,15 +222,10 @@ void teca_vertical_integral::set_properties( (void) prefix; (void) opts; + TECA_POPTS_SET(opts, std::string, prefix, integration_variable) TECA_POPTS_SET(opts, std::string, prefix, long_name) TECA_POPTS_SET(opts, std::string, prefix, units) - TECA_POPTS_SET(opts, std::string, prefix, hybrid_a_variable) - TECA_POPTS_SET(opts, std::string, prefix, hybrid_b_variable) - TECA_POPTS_SET(opts, std::string, prefix, sigma_variable) - TECA_POPTS_SET(opts, std::string, prefix, surface_p_variable) - TECA_POPTS_SET(opts, std::string, prefix, p_top_variable) - TECA_POPTS_SET(opts, int, prefix, using_hybrid) - TECA_POPTS_SET(opts, float, prefix, p_top_override_value) + TECA_POPTS_SET(opts, std::string, prefix, pressure_level_variable) TECA_POPTS_SET(opts, std::string, prefix, output_variable_name) } @@ -397,38 +339,10 @@ teca_vertical_integral::get_upstream_request( if (req.has("arrays")) req.get("arrays", arrays); - if (using_hybrid) - { - // get the a coordinate - if (request_var(this->hybrid_a_variable, - std::string("hybrid_a_variable"), - &arrays) != 0 ) return up_reqs; - - // get the b coordinate - if (request_var(this->hybrid_b_variable, - std::string("hybrid_b_variable"), - &arrays) != 0 ) return up_reqs; - } - else - { - // get the sigma coordinate - if (request_var(this->sigma_variable, - std::string("sigma_variable"), - &arrays) != 0 ) return up_reqs; - } - - // get the surface pressure - if (request_var(this->surface_p_variable, - std::string("surface_p_variable"), - &arrays) != 0 ) return up_reqs; - - // only request p_top if it isn't being overriden - if (! this->p_top_override_value ){ - // get the model top pressure - if (request_var(this->p_top_variable, - std::string("p_top_variable"), - &arrays) != 0 ) return up_reqs; - } + // get the pressure coordinate + if (request_var(this->pressure_level_variable, + std::string("pressure_level_variable"), + &arrays) != 0 ) return up_reqs; // get the input array if (request_var(this->integration_variable, @@ -514,57 +428,13 @@ const_p_teca_dataset teca_vertical_integral::execute( nx << ", " << ny << ", " << nz << "]" << std::endl; #endif - int using_hybrid = this->using_hybrid; - - const_p_teca_variant_array a_or_sigma = nullptr; - const_p_teca_variant_array b_i = nullptr; + const_p_teca_variant_array pressure_level = nullptr; - if (using_hybrid) - { - // get the a coordinate - a_or_sigma = get_info_variable(this->hybrid_a_variable, - std::string("hybrid_a_variable"), - in_mesh); - if (!a_or_sigma) return nullptr; - - // get the b coordinate - b_i = get_info_variable(this->hybrid_b_variable, - std::string("hybrid_b_variable"), - in_mesh); - if (!b_i) return nullptr; - } - else - { - // get the sigma coordinate - a_or_sigma = get_info_variable(this->sigma_variable, - std::string("sigma_variable"), - in_mesh); - if (!a_or_sigma) return nullptr; - } - - // get the surface pressure - const_p_teca_variant_array p_s = - get_info_variable(this->surface_p_variable, - std::string("surface_p_variable"), - in_mesh); - if (!p_s) return nullptr; - - const_p_teca_variant_array p_top_array; - p_teca_variant_array p_top_override_tmp; - if (! this->p_top_override_value){ - // get the model top pressure - // TODO: p_top is a scalar: is const_p_teca_variant array correct? - // probably need to put this in information array - p_top_array = get_mesh_variable(this->p_top_variable, - std::string("p_top_variable"), - in_mesh); - if (!p_top_array) return nullptr; - } - else{ - p_top_override_tmp = teca_variant_array_impl::New(1); - p_top_override_tmp->set(0, this->p_top_override_value); - p_top_array = (const_p_teca_variant_array) p_top_override_tmp; - } + // get the pressure coordinate + pressure_level = get_info_variable(this->pressure_level_variable, + std::string("pressure_level_variable"), + in_mesh); + if (!pressure_level) return nullptr; #ifdef TECA_DEBUG std::cerr << teca_parallel_id() << "teca_vertical_integral::execute:" @@ -595,25 +465,15 @@ const_p_teca_dataset teca_vertical_integral::execute( const NT_INARR * p_input_array = dynamic_cast(input_array.get())->get(); - const NT_INARR * p_a_or_sigma - = dynamic_cast(a_or_sigma.get())->get(); - - const NT_INARR * p_b_i - = dynamic_cast(b_i.get())->get(); - - const NT_INARR * p_p_s - = dynamic_cast(p_s.get())->get(); - - const NT_INARR * p_p_top - = dynamic_cast(p_top_array.get())->get(); + const NT_INARR * p_pressure_level + = dynamic_cast(pressure_level.get())->get(); NT_INARR * p_integrated_array = dynamic_cast(integrated_array.get())->get(); // call the vertical integration routine - vertical_integral(p_input_array, nx, ny, nz, - using_hybrid, p_a_or_sigma, p_b_i, p_p_s, p_p_top[0], - p_integrated_array); + mass_integral_trapezoidal(p_input_array, nx, ny, nz, p_pressure_level, + p_integrated_array); ) #ifdef TECA_DEBUG diff --git a/alg/teca_vertical_integral.h b/alg/teca_vertical_integral.h index a0aeb2afa..ea6f9b92c 100644 --- a/alg/teca_vertical_integral.h +++ b/alg/teca_vertical_integral.h @@ -27,16 +27,9 @@ class teca_vertical_integral : public teca_algorithm // set the name of the the arrays involved in integration TECA_ALGORITHM_PROPERTY(std::string, long_name) TECA_ALGORITHM_PROPERTY(std::string, units) - TECA_ALGORITHM_PROPERTY(std::string, hybrid_a_variable) - TECA_ALGORITHM_PROPERTY(std::string, hybrid_b_variable) - TECA_ALGORITHM_PROPERTY(std::string, sigma_variable) - TECA_ALGORITHM_PROPERTY(std::string, surface_p_variable) - TECA_ALGORITHM_PROPERTY(std::string, p_top_variable) + TECA_ALGORITHM_PROPERTY(std::string, pressure_level_variable) TECA_ALGORITHM_PROPERTY(std::string, integration_variable) TECA_ALGORITHM_PROPERTY(std::string, output_variable_name) - // set whether the vertical coordinate is hybrid or sigma - TECA_ALGORITHM_PROPERTY(int, using_hybrid) - TECA_ALGORITHM_PROPERTY(float, p_top_override_value) protected: teca_vertical_integral(); @@ -45,15 +38,9 @@ class teca_vertical_integral : public teca_algorithm std::string long_name; std::string units; - std::string hybrid_a_variable; - std::string hybrid_b_variable; - std::string sigma_variable; - std::string surface_p_variable; - std::string p_top_variable; + std::string pressure_level_variable; std::string integration_variable; std::string output_variable_name; - int using_hybrid = true; - float p_top_override_value = -1; std::vector get_upstream_request( unsigned int port, const std::vector &input_md,